subroutine tkund(np,x,px,y,py,z,g,dv,pz,al,rho,kx,eps0,coup)

      implicit none
      logical*4 coup
      integer*4 np,ndiv,i,j
      real*8 x(np),px(np),y(np),py(np),z(np),dv(np),g(np),pz(np)
      real*8 x1(np),px1(np),y1(np),py1(np)
      real*8 dp(np),aln,arho(np),theta(np)
      real*8 al,rho,kx,eps0

c       write(*,*) 'rho=',rho,'kx=',kx
        ndiv = 50

        aln = al/ndiv
c        write(*,*) al, aln, ndiv

        do 120 i=1,np
          dp(i) = (g(i) + 1.d0)*(g(i) + 1.d0) -1.d0
          arho(i)    = rho*sqrt(2.d0)*(1+dp(i))
          theta(i)   = aln/arho(i)
 120    continue

c      write(*,*) 'tkund called',dp(1)
        do 100 i=1,np
         do 110 j=1, ndiv

            x1(i)  = x(i)
            px1(i) = px(i)
            y1(i)  = y(i)
            py1(i) = py(i)

            x(i) = x(i) + px(i)*aln/2.d0
            y(i) = y(i) + py(i)*aln/2.d0


            x1(i)  = x(i)
            px1(i) = px(i)
            y1(i)  = y(i)
            py1(i)= py(i)-sinh(2.d0*kx*y(i))
     $               *aln/2.d0/arho(i)/arho(i)/kx

            x1(i) = x1(i) + px1(i)*aln/2.d0
            y1(i) = y1(i) + py1(i)*aln/2.d0

            x(i)  = x1(i)
            px(i) = px1(i)
            y(i)  = y1(i)
            py(i) = py1(i)

 110      continue  
100     continue

      return
      end
c