! dmc apr 2004 -- all xja routines use bmnmx instead of bmidp -- code
!  now in terms of min and max b on flux surfaces, not midplane values
!******************** start file xjadep.for ; group depshf *************
!-------------------------------------------------------------------
! xjadep - loop up 1st orbit adjustment to deposition j
!   see depshf subroutine where table dxjdep is set up
!
real(fp) function xjadep(zxp,zthp,zrdep,zydep,zvion,zvplldep, &
     itrap,iuncon,zxksav,zeiav,zphiav,ixiav,zxhw)
  !
  !  dmc -- apr 2004 -- improvements were made that were necessary for
  !   reliable operation with sts:
  !     (a) instead of using the midplane, orbit intercepts are found along
  !         notional surfaces at bmin(xi) and bmax(xi), where bmin and bmax
  !         are the min and max mod(b) locations on the xi flux surface,
  !         respectively.  the geometrical location of the sequence of
  !         bmin and bmax values need not be a smooth surface in real space.
  !     (b) xjavdr -- drift velocity estimates at bmin/bmax orbit intercepts --
  !         no longer assumes midplane locality and allows for b .ne. bphi.
  !     (c) conserved momentum is m*vphi*r + e/c delta(psi) (m*vphi*r replaces
  !         m*vpll*r).
  !     (d) special treatment for "bean geometries" is discarded.
  !
  ! inputs
  ! =========
  ! common:  plasma coordinate system and certain parameters
  !  e.g. plasma zone index j=lcentr, ... , ledge
  !       plflxg(j)=poloidal flux enclosed within inner bdy of
  !         plasma flux zone j, = surface j
  !       phiprg(j)=electrostatic potential at surface j
  !
  !       bmnmx(ir) -- for ir.ge.nzones+1, flux surface b minima;
  !                    for ir.le.nzones+1, flux surface b maxima
  !          (analogous to the old bmidp which implicitly expected the
  !           minima and maxima to be on the midplane).
  !
  !  also abeam = atomic weight of beam ions
  !          xe = charge of an electron
  !
  !  this list may not be complete
  !
  !  xjasub_mod module *** set up by call to subroutine xjaini
  !
  ! -----
  !  input (passed)  ** in the lab frame **
  !  ===============
  !  initial beam ion orbit parameters are passed:
  !
  !     zxp, zthp  "xi", "theta" flux coordinates of orbitting ion
  !     zrdep      major radius of orbitting ion (where deposited)
  !                (depends on zxp,zthp but is not re-evaluated for speed)
  !                cm
  !     zydep      elevation of ion above/below midplane  cm
  !     zvion      velocity of ion  cm/sec
  !     zvplldep   parallel component of velocity  cm/sec
  !
  !  output (passed) ** in the lab frame **
  !  ===============
  !
  !   itrap = 1  this is estimated to be a trapped orbit
  !   itrap = 0  estimated to be a passing orbit
  !
  !  iuncon = 0  orbit is confined
  !    then  zxksav  is vpll/v defined for fpp by "mu conservation" at
  !                  the bmin intercept of the orbit average
  !                  flux surface
  !          ixiav   is the index to the transp flux zone in which in
  !                  the orbit average flux surface lies
  !          xjadep  returns the "xi" value of the orbit average surface
  !                  (0.0 is the magnetic axis, 1.0 the plasma bdy)
  !                  the usual relative toroidal flux coordinate
  !          zeiav   returns the average energy of the ion in the orbit
  !                   ev
  !          zphiav  returns the average mechanical momentum of the ion
  !                   gram*cm**2/sec
  !
  !          zxhw    orbit half width (ximax-ximin)/2
  !
  !  iuncon = 1  orbit is unconfined
  !    then  zxksav  is left as zero
  !          ixiav   is left as lep1=ledge+1 (index to outside the plasma)
  !          xjadep  returns 1.0
  !
  !          zeiav   returns the energy of the ion leaving the plasma
  !                (i.e. crossing the xi=1 surface)
  !                   ev
  !          zphiav  returns the momentum of the ion leaving the plasma
  !                   gram*cm**2/sec
  !          zxhw    returns 0.0
  !
  !---------------------------
  !
  !  this routine is called by dep in the monte carlo beam code
  !c   and by other nubeam subroutines
  !
  !---------------------------
  !
  use nbi_com
  use physconst_mod, only : mp_g
  use xstraln_calls
  !
  use xjasub_mod
  !
  !============
  implicit none
  integer :: itrap,iuncon,ixiav,idir,ix0,ix0p1,iwarn,iout,i1,i2
  integer :: ixiavp,icout,icoutp,ierr
  !============
  real(fp) :: zthp,zrdep,zydep,zvion,zvplldep,zxksav,zeiav,zphiav
  real(fp) :: zxhw,zxp,zx,zth,zjx0,zxiav
  real(fp) :: zxfac,zphiang,zbr,zbz,zbphi,zbdepi,zbtb,zpsimc,zpsimg
  real(fp) :: zvpdp2,zr1,zr2,zvpll1,zvpll2,zvd1,xjavdr,zvd2,zvdlim
  real(fp) :: zvdi1,zvdi2,zden,zfacav,zpsi,zepav,zde,zvnew,zbout
  real(fp) :: zvn2,zvpl2,zvpllm
  real(fp) :: zxmp1,zxmp2
 !
  iuncon=0     ! rga: define in case of error
  if(.not.xjasubs_init) call xjaini
  !
  zx=zxp
  zth=zthp
  !
  !  radial interpolation set-up
  !
  zjx0=lcentr+zx*nzones
  !
  !  initial orbit point must be from inside the plasma
  !
  if(zx.lt.0.0_fp) then
     write(nonlin,*) ' ?xjadep: orbit startpoint at x < 0: ',zx
     zx=0.0_fp
  end if
  !
  if(zx.gt.1.0_fp) then
     write(nonlin,*) ' ?xjadep: orbit startpoint at x > 1: ',zx
     zx=1.0_fp
  end if
  !
  !  default outputs if orbit is lost ...
  zxksav=0.0_fp
  zxiav=1.0_fp
  ixiav=lep1
  zxhw=0.0_fp
  !
  !  sign of toroidal velocity
  !
  if(zvplldep.gt.0.0_fp) then
     idir=1
  else
     idir=-1
  end if
  !
  !  interpolation gub
  !
  ix0=zjx0
  ix0=min(ledge,ix0)
  !
  ix0p1=ix0+1
  !
  zxfac=(zx-xi(ix0,1))/dxi(ix0,2)
  !
  !  kinetic, potential energy (cf radial electrostatic potential)
  !
  zekin=zvion**2
  zepot=xzbeami* &
       (phiprg(ix0)+zxfac*(phiprg(ix0p1)-phiprg(ix0)))
  !
  !----------------
  !  compute conserved quantities:
  !    canonical momentum ("webers"):
  !
  !    magnetic moment
  !    ...get (1/b)[r,z] from fast map
  !
  zphiang=0.0_fp
  call xs_brz(0.01_fp,zphiang,zrdep,zydep,1.0_fp,zbr,zbz,zbphi, &
       iwarn)
  zbdepi=1.0_fp/sqrt(zbr**2+zbz**2+zbphi**2)
  zbtb = abs(zbphi)*zbdepi
  !
  zmass=mp_g*abeam
  !
  zpsimc=-zmass/(1.e7_fp*xe)*zbtb*zvplldep*zrdep
  zpsimg=xzbeami* &
       (plflxg(ix0)+zxfac*(plflxg(ix0p1)-plflxg(ix0)))
  !
  zpsito=zpsimg+zpsimc   ! conserved canonical momentum
  !
  zvpdp2=zvion**2-zvplldep**2
  !
  zmub=zvpdp2*zbdepi     ! conserved mu
  !----------------
  !
  !  find potential orbit midplane intercepts...
  !
  call xjaset(ix0,zvion)
  !
  iout=1
  !
  !  choose initial guess of which intercepts correspond to orbit
  !  deposited...
  !
  itrap=0
  call xjags0(zvplldep,i1,i2)
  !
  if(iloc(i1).ne.-1) then
  !
     ierr=1
     if(iloc(i2).eq.-1) then
        write(nonlin,9999) ierr
        xjadep=epsinv
        return
     endif
  !
  !  if trap intercept exists,
  !  test for trapping to outside-- based on midplane intercept fields
  !
     if(iloc(2).ne.-1) then
        if(zmub*zbi(i1).gt.zei(i1)) then
           i1=2
           i2=1
           itrap=1
        end if
     end if
  else
     i1=2
     i2=1
     itrap=1
  ! try a banana
  end if
9999 format(' ?? nubeam xjadep error ierr=',i10)
  !
  ! try a banana
  !
  ilconfine=.false.
  !
  ierr=2
  if(iloc(i1).eq.-1.or.iloc(i2).eq.-1) then
     write(nonlin,9999) ierr
     xjadep=epsinv
     return
  endif
  !
  iuncon=1
  if((iloc(i1).eq.0).or.(iloc(i1).ge.i2nzp1).or. &
       (iloc(i2).eq.0).or.(iloc(i2).ge.i2nzp1)) then
     zphiav=1.e7_fp*xe* &
          (xzbeami*plflxg(lep1)-zpsito)   ! exit momentum
     zeiav=v2toev*zekin+ &
          zepot-xzbeami*phiprg(lep1)   ! exit energy
     xjadep=zxiav
     return
  endif

  !
  iuncon=0
  !
  !  dmc -- save results in common
  !
  ilconfine=.true.
  iltrap = itrap.eq.1
  ilinside = iout.eq.-1
  !
  impin=i1
  impout=i2
  !
  zr1=zri(i1)
  zr2=zri(i2)
  zxmp1=zxidpi(i1)
  zxmp2=zxidpi(i2)
  zvpll1=zvplli(i1)
  zvpll2=zvplli(i2)
  !
  !  net "drift" velocities
  !  actually vert. drift velocity + e/b velocity + vert. component
  !    of pll velocity
  !  at midplane intercept pts
  !
  zvd1=xjavdr(i1)
  zvd2=xjavdr(i2)
  !
  !
  !---------------------------------------------------------------------
  !  *confined orbit*
  !
  !  orbit avg parameters
  !
  !  avg "r/a" from 1/vd weighted avg of "r/a" at midplane intercepts
  !  careful, vd can change sign or tend to zero
  !
  zvd1=abs(zvd1)
  zvd2=abs(zvd2)
  zvdlim=0.01_fp*max(zvd1,zvd2)
  zvd1=max(zvd1,zvdlim,epslon)
  zvd2=max(zvd2,zvdlim,epslon)
  !
  zvdi1=1.0_fp/zvd1
  zvdi2=1.0_fp/zvd2
  zden=(zvdi1+zvdi2)
  !
  !  avg xi="r/a"
  zxiav=(zvdi1*zxmp1+zvdi2*zxmp2)/zden
  zxhw=abs(zxmp2-zxmp1)*0.5_fp
  !
  !  indices for interpolation
  ixiav=min(ledge,int(lcentr+zxiav*nzones))
  ixiavp=ixiav+1
  zfacav=(zxiav-xi(ixiav,1))/dxi(ixiav,2)
  !
  !  avg momentum
  zpsi=xzbeami* &
       (plflxg(ixiav)+(plflxg(ixiavp)-plflxg(ixiav))*zfacav)
  zphiav=1.e7_fp*xe*(zpsi-zpsito)
  !
  !  new energy and velocity (taking into acct electrostatic potential
  !  change)
  zepav=xzbeami* &
       (phiprg(ixiav)+(phiprg(ixiavp)-phiprg(ixiav))*zfacav)
  zde=zepot-zepav   ! kinetic energy gain
  zvnew=sqrt(max(100.0_fp,(zvion**2+zde*zevtv2)))
  !
  zeiav=v2toev*zvnew*zvnew
  !
  icout=icen+iout*(ixiav-lcentr)
  icoutp=icout+iout
  !
  !  compute "avg" vpll/v: defined by mu conservation at the inner/outer
  !  midplane intercept  (inner only if inside-trapped banana, cf "iout")
  !
  zbout=bmnmx(icout)+(bmnmx(icoutp)-bmnmx(icout))*zfacav
  zvn2=zvnew**2
  zvpl2=zvn2-zmub*zbout
  !
  if(zvpl2.gt.1.e-6_fp*zvn2) then
     zvpllm=idir*sqrt(zvpl2)
     zxksav=zvpllm/zvnew
  else
     zxksav=0.001_fp
     zvpllm=zxksav*zvnew
  end if
  !
  xjadep=zxiav
  !
  return
end function xjadep
!******************** end file xjadep.for ; group depshf ***************
