!******************** start file depsub.for ; group depsub *************
!----------------------------------------------------------------
!  set up parameters for new orbit - called from dep or nutrav
!                                           or depshf
!
subroutine depsub(zlndep,zsinth,zcosth,zsinph,zcosph,iflag, &
     zphi0,zphigc,ze0,zegc,zbt0,zbtgc)
  !
  !  rga oct00 - removed levgb = 1,2
  !  dmc apr04 -- recommented & checked
  !
  !  pam stubberfield uda mod for use of btiltc(4,..) 14 july 1994
  !
  !  arguments -
  !   ** if iflag =0 or 1 **
  !   zlndep - trajectory length from trkrb starting pt (rsta,ysta)
  !            to depositon pt
  !   zsinth,zcosth,zsinph,zcosph:  sin and cos of thea, phia
  !            (trkrb direction angles)
  !   (this is how depsub is used when called from dep or nutrav)
  !
  !   ** if iflag =2 or 3 **
  !   the first 5 arguments have different meaning:
  !   in the code, the call to nb0loc is skipped and:
  !     zlndep --> rmjdep (major radius of deposition)
  !     zsinth --> ydep (deposition z location)
  !     zcosth --> zvtor (toroidal velocity)
  !     zsinph --> zvr   (velocity in r direction)
  !     zcosph --> zvz   (velocity in z direction)
  !   (this is how depsub is used when called from depfus)
  !
  !   iflag=0 or =2:  compute flr effect
  !   iflag=1 or =3:  ignore flr effect
  !   otherwise use trkrb system to include flr
  !
  !   outputs:
  !     zphi0:  m*vphi*r (ang. momementum) of neutral as deposited
  !     zphigc:  ang. momentum of ion (guiding ctr)
  !     ze0:  kinetic energy of neutral as deposited
  !     zegc:  kinetic energy of ion (guiding ctr)
  !     zbt0:  tor. field at neutral deposition point
  !     zbtgc:  tor. field at guiding ctr
  !
  !  common
  !  ------
  !
  !   input...
  !    vion = *neutral velocity*
  !
  !   output...
  !    (xiflr,thflr)  magnetic coordinates of ionization pt
  !    (xion,th)      magnetic coordinates of ion guiding ctr
  !    (rmjdep,ydep)  (r,y) coordinates of ionization pt
  !    (rmjion,yion)  (r,y) coordinates of guiding ctr
  !    vion = *ion guiding ctr velocity*  ****modified from input****
  !    xksid = vpll/vion
  !    (vion_dep,xksid_dep)  v,vpll/v at ionization pt
  !    ngc,nflr -- zone indices of gc and ionization pt., respectively
  !
  use nbi_com
  !
  !  b field components
  !============
  implicit none
  integer :: iflag,inflr,iermap,idummy,ingc
  !============
  real(fp) :: zphi0,zphigc,ze0,zegc,zbt0,zbtgc,zlndep
  real(fp) :: zvtor,zvr,zvz,zbmod,zvpll,zrgc,zrlov,zdt,zdr,zdz,zzgc
  real(fp) :: zpph1,fnbint,zpph2,zdelp,zephi0,zephi,zvion,zxbmbnd
  real(fp) :: zdrflr,zdzflr,ztop,z11,z12,z21,z22,zdet,zdrx,zdrth
  real(fp) :: zrlar,zsinth,zcosth,zsinph,zcosph,zxhor,zyhor
  !============
  real(fp) :: zbdir(4)  ! (3) used; cf nblocb subroutine
  !  vxb/|b| compenents
  real(fp) :: zvxbn(3)
  !
  !---------------------------------------------------------------
  !
  !	establish initial orbit parameters, guiding ctr
  !
  !  first get neutral position and velocity in (r,z) lab frame
  !
  if(iflag.le.1) then
     inflr=iflag
     call nb0loc(zlndep,zsinth,zcosth,zsinph,zcosph, &
          zxhor,zyhor,rmjdep,ydep,zvtor,zvr,zvz)
     x_hor=zxhor
     y_hor=zyhor
  else
     !  remap arguments to local variables
     inflr=iflag-2
     rmjdep=zlndep
     ydep=zsinth
     zvtor=zcosth
     zvr=zsinph
     zvz=zcosph
  end if
  !
  !  get magnetic coordinates and local field at ionization pt
  !
  call nblocb(rmjdep,ydep,xiflr,thflr,zbmod,zbdir,1,nfast_map, &
       iermap)              !zbmod in gauss
  if(iermap.gt.1) then
     write(nonlin,*) &
          ' ?? nubeam depsub: unexpected (r,z)->(xi,theta) map error!'
  end if
  !
  !  neutral momentum and energy
  vion_dep = vion
  zbt0=1.e-4_fp*zbmod*zbdir(1)
  zphi0=rmjdep*zvtor*rmassb
  ze0=v2toev*vion*vion
  !
  !  new vpll
  !    definition takes into account poloidal field
  !    note sign conventions cf. transp definition of vpll is relative
  !    to direction of toroidal current, not b
  !
  !  coordinate system in right-handed order:  (r^,phi^,z^);
  !    (r^,z^) are unit vectors in the poloidal plane which is typically
  !    represented to the right of the tokamak axis of symmetry (r=0);
  !    phi^ goes around toroidally, counter-clockwise (ccw) seen from above.
  !
  !  let sigb = nsnccw here = +1 if bphi is ccw from above, -1 otherwise
  !  let sigj = nsnccw*nsjdotb = +1 if jphi is ccw from above, else -1
  !  let sigjb = nsjdotb = sigb*sigj = +1 if jphi || bphi; -1 if anti-parallel
  !
  !  the "vpll" used in nubeam is generally nsjdotb*vpll i.e. positive for
  !  particles moving in the direction of jphi, negative for the others
  !
  !  similarly, here zvtor = nsjdotb*vphi
  !
  !  v  =  zvr*r^        + zvz*z^       + sigj*zvtor*phi^
  !  b  =  -sigj*{br}*r^ - sigj*{bz}*z^ + sigb*|bphi|*phi^
  !
  !  where {br} = |br|*(dr/dtheta)/(dl/dtheta)
  !  and   {bz} = |bz|*(dz/dtheta)/(dl/dtheta)
  !  are related to zbdir(2) and zbdir(3) returned from nblocb; -sigj gets
  !    the field to curl the right way around the enclosed plasma current.
  !
  !  (theta describes a flux surface going ccw looking at a poloidal plane
  !  drawn to the right of the tokamak axis of symmetry)
  !
  !  zvpll = nsjdotb*(v.b) = -sigb*(zvr*{br}+zvz*{bz}) + zvtor*|bphi|
  !
  zvpll=zvtor*zbdir(1)-nsnccw*(zvr*zbdir(2)+zvz*zbdir(3))
  !                 |bphi|/|b|           {br}/|b|     {bz}/|b|

  xksid_dep = zvpll/vion_dep

  !
  !  guiding ctr position (magnetic coords output)
  !  zvbxn(...) to contain:  v x b / |b|
  !
  !  v x b:

  !   |     r^         phi^           z^    |
  !   |    zvr      sigj*zvtor       zvz    |
  !   | -sigj*{br}  sigb*|bphi|  -sigj*{bz} |
  !
  !  = r^*[ -zvtor*{bz} - sigb*zvz*|bphi| ]      ! /|b| --> zvbxn(2)
  !   +phi^*(-sigj)*[ zvz*{br} - zvr*{bz} ]      ! /|b| --> zvbxn(1)
  !   +z^*[ sigb*zvr*|bphi| + zvtor*{br} ]       ! /|b| --> zvbxn(3)
  !
  xion=xiflr
  th=thflr
  !
  if(inflr.eq.1) then
     !
     zbtgc=zbt0
     zegc=ze0
     !
     ngc=min(ledge,int(lcentr+xion*nzones))
     nflr=ngc
     !
     zrgc=rmjdep
     !mg bug fixed 09/13/11
     zzgc=ydep
     !
  else
     zvxbn(1)= -nsjdotb*nsnccw*(zvr*zbdir(3)-zvz*zbdir(2))  ! sign of this
     !  term doesn't affect result, cf zdt**2 below...

     zvxbn(2)=-zvz*zbdir(1)*nsnccw-zvtor*zbdir(3)
     zvxbn(3)= zvtor*zbdir(2)+zvr*zbdir(1)*nsnccw
     !
     !  larmour displacement r,z components ! larmour radius / vperp:
     if(zbmod.gt.0.0_fp) then
        zrlov=abeam/(xzbeami*9.578e3_fp*zbmod) !zbmod in gauss
        zdt=zvxbn(1)*zrlov
        zdr=zvxbn(2)*zrlov
        zdz=zvxbn(3)*zrlov
     else
        zdt=0.0_fp
        zdr=0.0_fp
        zdz=0.0_fp
     end if
     !
     !  test code-- verify magnitude and direction of vxb, flr displacement
     !test	  zdz=zvxbn(1)*zrlov
     !test	  zvprp=sqrt(vion*vion-zvpll*zvpll)
     !test	  zvprp2=sqrt(zvxbn(1)**2+zvxbn(2)**2+zvxbn(3)**2)
     !test	  zdot1=zvxbn(1)*zbdir(1)-
     !test     >      nsnccw*(zvxbn(2)*zbdir(2)+zvxbn(3)*zbdir(3))
     !test	  zdot2=(zvxbn(1)*zvtor+zvxbn(2)*zvr+zvxbn(3)*zvz)/vion
     !test	  zrlar=zrlov*zvprp
     !test	  zrlar2=sqrt(zdr**2+zdz**2+zdz**2)
     !
     !       zrgc_old=rmjdep+zdr
     zrgc=sqrt((rmjdep+zdr)**2 + zdt**2) ! include displacement from r,z plane
     zzgc=ydep+zdz
     !
     call nblocb(zrgc,zzgc,xion,th,zbmod,zbdir,1,nfast_map,iermap)
     if(iermap.gt.1) then
        write(nonlin,*) &
             ' ?? nubeam depsub: unexpected (r,z)->(xi,theta) map error!'
     end if
     zbtgc=1.e-4_fp*zbmod*zbdir(1)
     !
     !  change in vpll:  momentum conservation
     !  evaluate change in e/c * psi(poloidal)
     !
     !  p_dep = m*zvtor*rmjdep
     !  p_gc = p_dep+mp*zb*(e/c)*(psi_gc-psi_dep)
     !  <vphi_gc> = p_gc/(m*r_gc) = vpll_gc*|bphi|/|b|
     !
     call psi_spline(xiflr,inflr,zpph1)
     call psi_spline(xion,ingc,zpph2)
     !
     zdelp=9.578e11_fp*rmassb/abeam* &
          xzbeami*(zpph2-zpph1)
     if(zbdir(1).gt.0.0_fp)then
        zvpll=(zphi0+zdelp)/(zrgc*zbdir(1)*rmassb)
     else
        zvpll=0.0_fp
     end if
     !
     !  change in kinetic energy due to electrostatic potential
     !
     zephi0=xzbeami*fnbint(xiflr,phiprg,1,mj,1,idummy,nflr)
     if(xiflr.gt.1.0_fp) nflr=lep1    ! bugfix dmc 5 oct 1994
     !
     zephi=xzbeami*fnbint(xion,phiprg,1,mj,1,idummy,ngc)
     !
     zegc=ze0+zephi0-zephi
     if(zegc.le.1.0_fp) then
        write(nonlin,8901) ze0,zegc
8901    format( &
             ' %depsub:  flr approximation in electrostatic field'/ &
             '  yielded negative energy (patched):'/ &
             '    ze0 = ',1pe11.4,'  zegc = ',1pe11.4,' ev.')
        zegc=1.0_fp
     end if
     !
  end if
  !
  !  ion position
  !
  yion=zzgc
  rmjion=zrgc
  !
  !  g.c. velocity, vpll/v
  !
  zvion=vion
  if(zegc.ne.ze0) vion=sqrt(zegc/v2toev)
  !
  xksid=zvpll/vion
  !
  if(abs(xksid).gt.1.0_fp) then
     if(abs(xksid).gt.1.001_fp) then
        write(nout,9901) xksid
9901    format(' %depsub - xksid = ',1pe13.6)
     end if
  end if
  !
  xksid=max(-0.999999_fp,min(0.999999_fp,xksid))
  !
  !  ion momentum
  !
  zphigc=rmassb*zbdir(1)*(xksid*vion)*zrgc
  !
  !  if this is a heavy ion, determine the gyrophase angle at
  !  deposition time -- needed by orbit
  !
  if(nlheavy(lsbeam)) then
     !  get 2d grid indexing
     zxbmbnd=xbmbnd
     xbmbnd=xminbm
     call nbigc(-1)
     xbmbnd=zxbmbnd
     zdrflr=rmjdep-zrgc
     zdzflr=ydep-zzgc
     if(.not.nlsym2b) then
        ztop = 1.0_fp
     else
        if(yion.ge.0.0_fp) then
           ztop=1.0_fp
        else
           ztop=-1.0_fp
        end if
     end if
     !
     !  solve system:  zdrflr=z11*zdrx+z12*zdrth
     !                 zdzflr=z21*zdrx+z22*zdrth
     !    zdrx and zdrth are the displacement in the grad(xi) direction
     !    and the flux-surface-tangent direction respectively; from these
     !    the cos and sin of the gyrophase angle are defined.  see usage
     !    of cospho and sinpho in subroutine nbiflr
     !
     z11=btiltc(5,ngc2)
     z12=ztop*btiltc(4,ngc2)
     z21=-ztop*btiltc(4,ngc2)
     z22=z11
     !
     zdet=z22*z11-z12*z21
     zdrx=(zdrflr*z22-zdzflr*z12)/zdet
     zdrth=(zdzflr*z11-zdrflr*z21)/zdet
     !
     zrlar=0.0_fp
     cospho=zdrx
     zrlar=zrlar+cospho**2
     sinpho=zdrth*nsnccw/btiltc(2,ngc2)
     zrlar=zrlar+sinpho**2
     !
     zrlar=sqrt(zrlar)
     cospho=cospho/zrlar
     sinpho=sinpho/zrlar
     !
  end if
  !
  return
  !
end subroutine depsub
!******************** end file depsub.for ; group depsub ***************
