subroutine drift(dwtp,iend,ibeg)
c----------------------------------------------------------------------
c
      include 'param_sz.h'
      include 'bfieldcom.h'
      include 'constcom.h'
      include 'flagcom.h'
      include 'misccom.h'
      include 'pcordcom.h'
      include 'pfieldcom.h'
      include 'syscom.h'
      include 'ucom.h'
c
c--------------------------------------------------------------------------
c*
      ne=rne
      bz=bgz/gamma
      dz=bz*dcon*dwtp
      if(z+dz.ge.zloc(ne).and.ne.le.nel)then
c----particle will cross into next element during this step.
        iend=1
        dummy=0.
        if((ip.eq.999).and.(nupar.eq.1))
     *      write(nimp,*) z,dummy,dummy,dummy
        dz=zloc(ne)-z
        dwtp=dz/(bz*dcon)
      endif
c-----allow particle to go backwards
      if(ne.ge.1.and.dz.lt.0.and.z+dz.le.zloc(ne-1))then
c----particle will cross into previous element during this step
        ibeg=1
        dz=zloc(ne-1)-z
        dwtp=dz/(bz*dcon)
      endif
      if(poiflag .and. z.ge.pzmin .and. z.le.pzmax .and. ne.le.nel
     * .and. z.ge.0.) go to 10
      if(ifld.eq.0 .or. z+dz.lt.zmin .or. z.gt.zmax .or. ne.gt.nel
     * .or. z.lt.0.)then
         x=x+dz*bgx/bgz
         y=y+dz*bgy/bgz
         z=z+dz
         if(iend.eq.1)z=zloc(ne)
         if(ibeg.eq.1)z=zloc(ne-1)
         go to 2
      endif
 10   continue
      if(dwtp.eq.0.) go to 2
      beta=sqrt(1.-1./gamma**2)
      bx=bgx/gamma
      by=bgy/gamma
      dx=bx*dcon*dwtp
      dy=by*dcon*dwtp
      dl=sqrt(dx*dx+dy*dy+dz*dz)
      xp=dx/dl
      yp=dy/dl
      zp=dz/dl
      r=sqrt(x**2+y**2)
      temp=brhof*sqrt(gamma**2-1.)
      if(temp.ne.0.) then
        if(ifoclal.eq.1) then    ! foclal
           cayz=bfldlal(z)/temp     
           cayzn=bfldlal(z+dz)/temp 
           cayr=0.
        else                     ! Poisson and Helmholz
          cayz=bfld(z,r,brfld)/temp
          cayzn=bfld(z+dz,r,brfldn)/temp
           if(r.gt.1.e-10) then
             cayr=-(brfld+brfldn)*dl/(r*temp)
           else
             cayr=0.
           endif
        endif
      else
        cayz=0.
        cayzn=0.
      endif
c----drift half way.
      dl=dl*.5
      if(cayz.eq.0.)then
        x=x+dl*xp
        y=y+dl*yp
      else
        cayl=cayz*dl
        s=sin(cayl)
        c=cos(cayl)
        x = x + (s*xp + (1.-c)*yp)/cayz
        y = y + (s*yp - (1.-c)*xp)/cayz
        xxp= c*xp + s*yp
        yp = c*yp - s*xp
        xp=xxp
      endif
      z=z+dl*zp
c----apply fringe transformation
      if(ifoclal.eq.1) then
          xp=xp+.5*(cayzn-cayz)*y*zp
          yp=yp-.5*(cayzn-cayz)*x*zp
          zp=zp-.5*(cayzn-cayz)*(xp*y-yp*x)
c------allow particle to go backwards
          fac=1./sqrt(xp**2+yp**2+zp**2)
          xp=xp*fac
          yp=yp*fac
          zp=zp*fac
      else
          xp=xp+.5*cayr*y
          yp=yp-.5*cayr*x
          fac=1.-(xp**2+yp**2)
          if(fac.ge.0.) then
             zp=sign(sqrt(fac),zp)
          else
c---------particle has turned around
             xp=xp-.5*cayr*y
             yp=yp+.5*cayr*x
             zp=-zp
          endif
      endif
      bz=zp*beta
c----drift remaining distance.
      if(iend.eq.1 .or. ibeg.eq.1)then
c----adjust dwtp
        dwtp=.5*dwtp+dl*zp/(bz*dcon)
      else
c----adjust dz and check if new dz makes particle cross into next element.
        dz=.5*dwtp*bz*dcon
        if(z+dz.ge.zloc(ne))then
c----particle will cross into next element with new dz therefor adjust dz
c----and dwtp.
          iend=1
          dztemp=zloc(ne)-z
          dwtp=.5*(dwtp+dztemp/dz*dwtp)
          dl=dl*dztemp/dz
          dx=dl*xp
          dy=dl*yp
          dz=dztemp
        endif
        if(z+dz.le.zloc(ne-1))then
c----particle will cross into previous element with new dz therefor adjust dz
c----and dwtp.
          ibeg=1
          dztemp=zloc(ne-1)-z
          dwtp=.5*(dwtp+dztemp/dz*dwtp)
          dl=dl*dztemp/dz
          dx=dl*xp
          dy=dl*yp
          dz=dztemp
        endif
      endif
      if(cayzn.eq.0.)then
        x=x+xp*dl
        y=y+yp*dl
      else
        cayl=cayzn*dl
        s=sin(cayl)
        c=cos(cayl)
        x = x + (s*xp + (1.-c)*yp)/cayzn
        y = y + (s*yp - (1.-c)*xp)/cayzn
        xxp= c*xp + s*yp
        yp = c*yp - s*xp
        xp = xxp
      endif
      bx=xp*beta
      by=yp*beta
      bgx=bx*gamma
      bgy=by*gamma
      bgz=bz*gamma
 1    z=z+dz
      if(iend.eq.1)z=zloc(ne)
cbm 07/03/2001      if(ibeg.eq.1)z=zloc(ne-1)
      if((ne.gt.1).and.(ibeg.eq.1))z=zloc(ne-1)
 2    continue
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*