! 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
!----------------------------------------------------------------------
!  xjavdr  compute midplane intercept "drift" velocity
!
real(fp) function xjavdr(ix)
  !
  !  input ix = ptr to intercept data in xjasubs.blk common
  !  output xjavdr = drift velocity in poloidal direction
  !    = (vpll * bp/b  + c e x b/b**2  + vdrift).epol-unit-vector
  !      with vdrift=(v**2-vperp**2/2)/(omega) * (grad(b)/b)
  !      where omega=eb/m
  !
  !  dmc apr 2004 -- generalized for off-midplane bmin and bmax surfaces...
  !
  !    pll drift:  |bp|/|b|*(vpll/v)  always positive
  !    exb drift:  (|bt|/|b|)*[d(epot)/dr]/|b|  sign is that of d(epot)/d(psi)
  !  grad(b)+curv(b) drift:  (|bt|/|b|)*gbb/(vperp**2-(1/2)*vpll**2)/omega
  !    where gbb = 1/|b| * d(bmax)/dr  if at bmax on surface;
  !    where gbb = 1/|b| * d(bmin)/dr  if at bmin on surface;
  !       (sign relative to epol-unit-vector generally switches from bmin to
  !        bmax side)
  !
  !---------------------------
  !
  use nbi_com
  use physconst_mod, only : aee_amp! electron_charge/proton_mass (c*kg^-1)
  !
  use xjasub_mod
  !
  !---------------------------
  !
  !============
  implicit none
  integer :: ix,imp,i1,i2,ij
  !============
  real(fp) :: zxi,zinrho,zinmp,zdelbp,zdelep,zdelbm,zdelem,zb,zvpll
  real(fp) :: zbpb,zbtb,zvdpll,zdr,zdelv,zefld,zvexb,zgradb,zv2,zvp2
  real(fp) :: zomega,zvdrft,zv
  !============
  imp=iloc(ix)
  zxi=zxidpi(ix)
  !
  i1=imp    ! guaranteed to be .lt. 2*nzones+1 -- i.e. inside plasma
  i2=imp+1  ! guaranteed to be .le. 2*nzones+1
  !
  !  inboard/outboard sign; eletrostatic field strength
  !
  if(i1.lt.icen) then
     ij=lcentr+icen-i2
     zinrho=(zxi-xi(ij,1))/dxi(ij,2)
     zinmp=1.0_fp-zinrho
     if(i1.eq.1) then
        zdelbp=bmnmx(i1)-bmnmx(i1+1)
        zdelep=phiprg(ij+1)-phiprg(ij)
     else
        zdelbp=(bmnmx(i1-1)-bmnmx(i1+1))/2
        zdelep=(phiprg(ij+2)-phiprg(ij))/2
     end if
     zdelbm=(bmnmx(i1)-bmnmx(i1+2))/2
     if(ij.eq.lcentr) then
        zdelem=0.0_fp
     else
        zdelem=(phiprg(ij+1)-phiprg(ij-1))/2
     end if
  else
     ij=lcentr+i1-icen
     zinrho=(zxi-xi(ij,1))/dxi(ij,2)
     zinmp=zinrho
     if(i2.eq.(2*nzones+1)) then
        zdelbp=bmnmx(i2)-bmnmx(i2-1)
        zdelep=phiprg(ij+1)-phiprg(ij)
     else
        zdelbp=(bmnmx(i2+1)-bmnmx(i2-1))/2
        zdelep=(phiprg(ij+2)-phiprg(ij))/2
     end if
     zdelbm=(bmnmx(i2)-bmnmx(i2-2))/2
     if(ij.eq.lcentr) then
        zdelem=0.0_fp
     else
        zdelem=(phiprg(ij+1)-phiprg(ij-1))/2
     end if
  end if
  !
  zb=zbi(ix)
  zvpll=zvplli(ix)
  !
  zbpb=(1-zinmp)*bpbmnmx(i1)+zinmp*bpbmnmx(i2)   ! bp/b
  zbtb=(1-zinmp)*btbmnmx(i1)+zinmp*btbmnmx(i2)   ! bt/b
  !
  zvdpll=zvpll*zbpb  ! pll drift
  !
  zdr=(1-zinmp)*drmnmx(i1)+zinmp*drmnmx(i2)
  zdelv=(1-zinrho)*zdelem+zinrho*zdelep
  zefld=zdelv/zdr    ! radial field
  !
  zvexb=1.0e4_fp*zefld*zbtb/zb   ! exb drift
  !
  zgradb=((1-zinrho)*zdelbm+zinrho*zdelbp)/zdr
  !
  zv2=zei(ix)
  zvp2=max(0.0_fp,(zv2-zvpll*zvpll))
  !
  zomega=xzbeami*aee_amp*1.0e-3_fp*zb/abeam !electron_charge/proton_mass (c*g^-1)
  !
  !  generalized from midplane & 1/r b field: grad(b)+curv(b) drift
  zvdrft=((zbtb*zgradb)/zb)*(zv2-0.5_fp*zvp2)/(zomega)
  !
  !  (in old code, grad(b)/b = -1/r was assumed and the sign showed up as
  !  -zvdrft in the following line)...
  !
  zv=zvdpll+zvexb+zvdrft
  !
  xjavdr=zv
  !
  return
end function xjavdr
