!  abdifu -- beam anomalous diffusion operator
!
subroutine abdifu(jbe,zdsdfb,zdsxvb,ierr)
  !
  !  return zdsdfb -- incremented.  sum for orbit goosing control (due to db).
  !  return zdsxvb -- incremented.  sum for orbit goosing control (due to vb).
  !
  !  rga oct00 - removed levgb = 1,2
  !
  !  updated last:
  !  09/21/93 tbt corrected pbdfbe to pbdfbi
  !
  !  dmc apr10 ...
  !    nmdifb.eq.4 version -- use pitch/energy/space dependent diffusion;
  !    modify velocity components to preserve both mu & etot=epot+ekin
  !--------------------------------
  !
  use nbi_com
  use nbi_random_pll
  !
  implicit none
  !
  !  passed arguments (input):
  !
  integer, intent(in) :: jbe          ! (fbm) energy zone index
  real(fp), intent(inout) :: zdsdfb     ! goosing sum: diffusion
  real(fp), intent(inout) :: zdsxvb     ! goosing sum: radial velocity
  integer, intent(out) :: ierr
  !
  !--------------------------------
  !
  !  this code implements a model for anomalous spatial diffusion of
  !  fast ions.
  !
  !  the guiding center of the ion is moved.
  !
  !  changes in the momentum and/or energy of the ion are accumulated
  !  in torque and/or power profiles for output
  !
  !  properties of the model controlled by namelist switch nmdifb
  !
  !
  !  nmdifb.eq.4:
  !    these are used to interpolate between:
  !        difbe_2d_fd0(x,e,specie) -- deeply trapped diffusivity, dtrap=1
  !        difbe_2d_fdb(x,e,specie) -- barely trapped diffusivity, dtrap=0
  !        difbe_2d_fdp(x,e,specie) -- barely passing (co) diffusivity, dpass=0
  !        difbe_2d_fdq(x,e,specie) -- deeply passing (co) diffusivity, dpass=1
  !        difbe_2d_fdr(x,e,specie) -- barely passing (ctr) diffusivity, dpass=0
  !        difbe_2d_fds(x,e,specie) -- deeply passing (ctr) diffusivity, dpass=1
  !
  !    also zdsdfb return is increased according to maximum diffusivity in
  !    vicinity of orbit, not just at the current orbit location...
  !
  !---------------------------------
  real(fp) :: zep,zphion,zxionp,zxi,zev,zdifb,zvelb,zveld,znorm,zdum
  real(fp) :: zintrp,zdlasq,zgxi2,zgxi,zdla,zgdr,zdr,zgdy,zdy,zdz,zr2
  real(fp) :: zr,zy,ztaui,zdddx,zdphi,zdlgxi,zdt,zd2sum,zdwt
  real(fp) :: zvbr,zvby,zvnorm
  real(fp) :: zephi,zmub
  integer :: igc,igcflg,igcflg0,i,im1,inp1,igcl,ict,ictmax,irat
  integer :: ie,ix, iok,i_repeat
  real(fp) :: zface,zfacx
  real(fp) :: zdifb_max,zde,zdx
  real(fp) :: zymp0,zymp1,zdwtp,zdsxvbp, zdsdfbp


  !---------------------------------
  !
  ierr=0
  !
  if((ngc.lt.lcentr).or.(ngc.gt.ledge)) return
  !
  if(nmdifb.le.0) return
  !
  !  code was broken out of collid dmc sept 1990 -- it executes in the
  !  lab frame, not the plasma frame.
  !
  !
  !  save current parameters
  !
  zep=ep          ! vion**2
  zev=zep*v2toev  ! kinetic energy, ev
  !
  zphion=bphion   ! orig. momentum
  zephi=ephion    ! orig. electrostatic potential
  zmub=vperp*vperp/bimodb(1)  ! orig. mu *2/ms
  !
  igc=ngc
  !
  zxionp=xion
  !
  !  height above/below midplane; abdifu will prevent this from changing
  !  sign, so that midplane crossings due to diffusion will not occur...
  !  this is for benefit of orbit goosing book keeping algorithm (dmc apr 2010).
  !
  !  (code here mimics bounce_ymp contained subroutine within orbit.for)
  !
  zymp0=ympx(ngcx)+xingcx*(ympx(ngcx+1)-ympx(ngcx))
  zymp0=yion-zymp0
  if(abs(zymp0).lt.epslon) zymp0=epslon
  !
  !  dmc may 1992 -- implement energy dependence of diffusivity
  !  this is in the lab frame
  !
  if(nmdifb.le.3) then
     zface=1.0_fp
     if(ndifbe.gt.0) then
        if(zev.lt.edifbe(1)) then
           zface=fdifbe(1)
        else if(zev.ge.edifbe(ndifbe)) then
           zface=fdifbe(ndifbe)
        else
           do i=2,ndifbe
              im1=i-1
              if((edifbe(im1).le.zev).and.(zev.le.edifbe(i))) then
                 !     piecewise linear interpolation
                 zface=fdifbe(im1)+(fdifbe(i)-fdifbe(im1))* &
                      (zev-edifbe(im1))/(edifbe(i)-edifbe(im1))
                 goto 10
              end if
           end do
        end if
     end if

  elseif(nmdifb .eq. 4) then

     ! nmdifb.eq.4:

     ! the energy interpolation factor is looked up just once;
     ! it is assumed that the energy variation due to diffusion
     ! will be small

     call eintrp(zev,ie,zface)

  end if
  !
10 continue
  !
  !    cf namelist control -- nb dmc mar 2008: use of nkdifb and x_qlim0
  !                           moved to subroutine parms
  !
  !                           separate profiles now provided for beam ions
  !                           and fusion product ions.
  !
  !    d(e,x) profiles provided separately for each specie
  !
  igcflg=1                  ! *hardwire this now*
  !
  !---------------------------------------------
  !  anomalous diffusion
  !    random scatter in (r,y); inverse map back to (xi,th)
  !
  zd2sum=zero
  !
  if((igcflg.eq.1)) then
     !
     ict=0
     ictmax=1
     zdwtp = 0.0_fp
     i_repeat = 0
     do
        zdsxvbp = zdsxvb
        zdsdfbp =  zdsdfb
        ict=ict+1
        if(ict.gt.ictmax) exit

        if(ngc.gt.lcentr) then
           ztaui=1.0_fp/ &
                (xbjaco(1,1)*xbjaco(2,2)-xbjaco(1,2)*xbjaco(2,1))
           zgxi2=(ztaui*xbjaco(2,2))**2+(ztaui*xbjaco(1,2))**2
        else
           ztaui=zero
           zgxi2=(dxi(lcentr,2)*dravizc(lcentr))**2
        end if
        zgxi=sqrt(zgxi2)
        !
        !  fetch beam diffusion coeff (cf s.r. fokker) @guiding ctr
        !
        if(nmdifb.le.3) then
           !
           !  1d diffusion coefficients
           !
           inp1=ngc+1
           zintrp=xingc             ! interpolation factor from nbigc
           if(nlfprod(lsbeam)) then
              zvelb=velb_fi(ngc)+zintrp*(velb_fi(inp1)-velb_fi(ngc))
              zdifb=difb_fi(ngc)+zintrp*(difb_fi(inp1)-difb_fi(ngc))
              zdddx=(difb_fi(inp1)-difb_fi(ngc))/dxi(ngc,2)
           else
              zvelb=velb_bi(ngc)+zintrp*(velb_bi(inp1)-velb_bi(ngc))
              zdifb=difb_bi(ngc)+zintrp*(difb_bi(inp1)-difb_bi(ngc))
              zdddx=(difb_bi(inp1)-difb_bi(ngc))/dxi(ngc,2)
           end if
           ! nmdifb.le.3: grad(db) used if ngradd_opt = 0 or 1
           if(ngradd_opt.eq.-1) zdddx=zero

        else if(nmdifb .eq. 4) then
           !
           !  2d diffusion coefficients
           !     (nmdifb.eq.4)
           !  d(e,x)...
           call xintrp(xion,ix,zfacx)

           zvelb=zero

           call get_fp(zdifb,zdddx)

           ! nmdifb.ge.4: grad(db) used if ngradd_opt = 1 only
           if(ngradd_opt.le.0) zdddx=zero


        else if(nmdifb .eq. 6) then
           !
           !  2d diffusion coefficients from dep code
           !
           iok = 0
           call xeintrp(xion,zev,ix,ie,zfacx,zface,iok)

           if(iok /= 0) then
              zvelb = zero
              zdifb = zero
              zdddx = zero

           else
              zvelb=zero

              call get_fp(zdifb,zdddx)

              ! nmdifb.ge.4: grad(db) used if ngradd_opt = 1 only
              if(ngradd_opt.le.0) zdddx=zero
           end if

        end if

        zdifb=max(zero,zdifb)
        !
        !   ! energy dependence --
        !
        if(nmdifb.le.3) then
           zvelb=zface*zvelb ! radial velocity
           zdifb=zface*zdifb ! db
           zdddx=zface*zdddx ! d(db)/dx
        end if
        !
        zveld=zdddx*zgxi    ! grad(db) velocity (nmdifbe.le.3 only)
        !
        !  d=0 => no motion => exit
        !
        if((zdifb.le.zero).and.(abs(zvelb).eq.zero)) then
           zdwt = wght*dtgfp
           wgvdx(ngc,lsbeam)=wgvdx(ngc,lsbeam)+zdwt
           exit
        end if
        !
        if(abs(zvelb).gt.zero) then
           zvnorm=sqrt(xbjaco(1,2)**2+xbjaco(2,2)**2)
           zvbr= xbjaco(2,2)*zvelb/zvnorm
           zvby=-xbjaco(1,2)*zvelb/zvnorm
        else
           zvbr=zero
           zvby=zero
        end if
        !
        !  grad(d) -- used to force flux of d*grad(n), not
        !    d*grad(n) + n*grad(d)
        !
        zgdr= xbjaco(2,2)*zdddx*ztaui
        zgdy=-xbjaco(1,2)*zdddx*ztaui
        !
        !  avg cross-surface displacement (as per rob's definition)
        !
        if((zdifb.gt.zero).or.(abs(zveld+zvelb).gt.zero)) then
           !
           if(ict.eq.1) then
              !  shorten timestep if it is going to be too large...
              !
              zdlasq=2.0_fp*dtgfp*zdifb
              zdla=sqrt(zdlasq)
              !
              zdlgxi=dxi(ngc,2)/zgxi   ! local surface spacing
              !
              !  nominal step -- 1/4 local zone width
              !
              irat = 4*zdla/zdlgxi + 1
              ictmax = irat*irat ! this many steps
              !
              !  also check velocity
              !
              zdla=dtgfp*abs(zveld+zvelb)
              irat = 4*zdla/zdlgxi + 1
              ictmax = max(irat,ictmax) ! no square here
              !
              zdt = dtgfp/ictmax
              !
           end if
           !
           zdlasq=2.0_fp*zdt*zdifb
           !  sum for goosing
           !  (nmdifb.le.3)
           if(nmdifb.le.3) then
              zdsdfb=zdsdfb+zdlasq*zgxi2
           end if
           !
           !  expected <dl> for diffusive scattering
           zdla=sqrt(zdlasq)
           !
        else
           !
           zdla=zero
           if(ict.eq.1) then
              zdt = dtgfp
              ictmax = 1
           end if
           !
        end if
        !

        zdwt = wght*zdt
        wgvdx(ngc,lsbeam)=wgvdx(ngc,lsbeam)+zdwt-zdwtp
        difbx(ngc,lsbeam)=difbx(ngc,lsbeam)+(zdwt-zdwtp)*zdifb
        velbx(ngc,lsbeam)=velbx(ngc,lsbeam)+(zdwt-zdwtp)* &
             (zvelb+zveld)

        zdsxvb=zdsxvb+abs(zvelb+zveld)*zdt*zgxi
        !
        !  3d scatter of guiding ctr
        !  (r,y) & z for toroidal direction
        !
        if(zdla.gt.zero) then
           zdr= nbgauss_pll(zdla,zero) ! horizontal in pol. plane
           zdy= nbgauss_pll(zdla,zero) ! veritical in pol. plane
           zdz= nbgauss_pll(zdla,zero) ! toroidal direction
        else
           zdr=zero
           zdy=zero
           zdz=zero
        end if
        !
        zdr = zdr + zdt*(zgdr+zvbr)
        zdy = zdy + zdt*(zgdy+zvby)
        !
        !  toroidal effect
        !
        if(zdz.ne.zero) then
           zr2=rmjion**2+zdz**2
           zdr=zdr+(sqrt(zr2)-rmjion)
        end if
        !
        !  move the g.c.
        !
        call nbimov(zdr,zdy,ierr)

        if(ierr.ne.0) then
           if(ierr.eq.999.and.i_repeat.eq.0) then
              ! repeat once, use rest of the time step
              ictmax=2
              ierr=0
              zdsxvb = zdsxvbp
              zdsdfb =  zdsdfbp
              zdwtp = zdwt
              zdt = dtgfp-(ict-1)*zdt
              ict=1
              i_repeat=1
              cycle
           else
              !  ion escapes  via anomalous diffusion
              ierr = 0
              ngc=0
              xion = xbmbnd
              return
           end if
        end if
        !
        zd2sum=zd2sum+zdr*zdr+zdz*zdz
        !
        if(ict.lt.ictmax) then
           if(nfast_map.eq.1) then
              call nbmomry(xion,th,1,zr,zy,xbjaco)     ! use bilinear maps
           else
              call nbmomry_xs(xion,th,1,zr,zy,xbjaco)  !use xplasma spline maps
           end if
        end if
        !
     end do
     !
     !----------------------------------------------------------
     !  epilog
     !
     !    store new position in orbit routine variables
     !
     if(zd2sum.gt.zero) then
        !
        !  dmc: undo midplane crossing if necessary
        !
        ict=1
        do
           bixi(1) = xion
           bithst(1) = th
           !
           !  since the guiding ctr has moved, re-evaluate auxilliary parameters
           !
           call nbigc(1)
           !
           !  check midplane height
           zymp1=ympx(ngcx)+xingcx*(ympx(ngcx+1)-ympx(ngcx))
           zymp1=yion-zymp1
           if(abs(zymp1).lt.epslon) zymp1=epslon
           !
           !  prevent diffusion induced midplane crossing, by nudging the gc position
           !
           if(zymp0*zymp1.lt.zero) then
              if(zymp1.gt.zero) then
                 ! move it back below...
                 zdy = -zymp1 - raxis*1.0e-5_fp*ict
              else
                 ! move it back above...
                 zdy = -zymp1 + raxis*1.0e-5_fp*ict
              end if

              zdr = zero
              call nbimov(zdr,zdy,ierr)
              if(abs(bixi(1)-xion).le.max(abs(bixi(1)), &
                   abs(xion))*1.e-8_fp.and.abs(bithst(1)-th).le. &
                   max(abs(bithst(1)),abs(th))*1.e-8_fp) then
                 !     stepsize is too small and xs_fxget called from nbimov decides that
                 !     within its tolerances the particle doesn't have to move.
                 !     now,  move the particle by raxis*ict*1.0e-5 instead of just raxis*1.0e-5.
                 !     that increases the step each time until presumably it is enough to move it.
                 ict=ict+1
              end if
              if(ierr.ne.0) return

              cycle
           else
              exit  ! on the correct side already, as usual
           end if

        end do
        !
        !    recompute the magnetics at the new position
        !
        call bintrp
        call jintrp
        if(bijac(1).le.zero) ierr=1
        !
        !    radial current sum
        !
        call radcur(zxionp,xion,wght,tirabds)
        !
        !    new velocity; kinetic energy modified s.t. ekin+epot unchanged
        !
        zev = max(one,(zev - (ephion-zephi)))
        vion = sqrt(zev/v2toev)
        !
        vperp = sqrt(bimodb(1)*zmub)
        vperp = min(vion,vperp)
        !
        if(vpll.lt.zero) then
           vpll = -sqrt(vion*vion-vperp*vperp)
        else
           vpll = sqrt(vion*vion-vperp*vperp)
        end if
        !
        !  update gc parameters again (due to velocity change)
        !
        call nbigc(0)
        !
        !---------------------------------------------------------
        !  %dmc sept 1990 -- this version assumes no change in velocity
        !  new momentum, new energy, output sums ...
        !
        !  ** the radial electrostatic field is ignored **
        !
        if(ngc.gt.0) then
           !
           igcl=(igc+ngc)/2
           !
           zdphi=(zphion-bphion)*wght
           bphdfbs(lsbeam)=bphdfbs(lsbeam)+zdphi
           tqdifbs(igcl,lsbeam)=tqdifbs(igcl,lsbeam)+zdphi
           bptdfbs(lsbeam)=bptdfbs(lsbeam)+ &
                0.5_fp*(omegag(igcl)+omegag(igcl+1))*zdphi* &
                1.e-7_fp/xe
           ! %dmc sept 1990 -- arbitrary 50-50 splitting btw ions, electrons:
           pbdfbes(igcl,lsbeam)=pbdfbes(igcl,lsbeam) + &
                0.5_fp*(zep-ep)*v2toev*wght
           !  pbdfbi(igcl)=pbdfbe(igcl)+0.5*(zep-ep)*v2toev*wght
           pbdfbis(igcl,lsbeam)=pbdfbis(igcl,lsbeam) + &
                0.5_fp*(zep-ep)*v2toev*wght ! tbt 9/93
           !
        end if
     end if
  end if
  !
  !-----------------------
  !  nmdifb.eq.4 and 6: scan nearby configuration space for diffusivities;
  !    use the maximum in computing the goosing sum
  !
  if(nmdifb.eq.6) then

     zdifb_max = zdifb
     zde = 0.5_fp*(efbmb(jbe+1,lsbeam)-efbmb(jbe,lsbeam))
     zdx = dxitk*0.5_fp

     !  energy and x variation @ past half bounce orbit topology

     ! db(e,x-dx)
     iok = 0
     call xeintrp(xion-zdx,zev,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)

     ! db(e,x+dx)
     iok = 0
     call xeintrp(xion+zdx,zev,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)


     ! db(e-de,x)
     iok = 0
     call xeintrp(xion,zev-zde,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)

     ! db(e+de,x)
     iok = 0
     call xeintrp(xion,zev+zde,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)



     ! topology variations: @ x, x-dx, x+dx
     call getdmax(zdifb)      ! @ x
     zdifb_max = max(zdifb_max,zdifb)

     iok = 0
     call xeintrp(xion-zdx,zev,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call getdmax(zdifb)      ! @ x-dx
     zdifb_max = max(zdifb_max,zdifb)

     iok = 0
     call xeintrp(xion+zdx,zev,ix,ie,zfacx,zface,iok)
     if(iok /= 0) return

     call getdmax(zdifb)      ! @ x+dx
     zdifb_max = max(zdifb_max,zdifb)

     !-----------------------------------------
     ! ok: zdsdfb sum for output:

     zdlasq = 2.0_fp*dtgfp*zdifb_max
     zdsdfb = zdsdfb + zdlasq*zgxi2

  else if(nmdifb .eq. 4) then

     zdifb_max = zdifb
     zde = 0.5_fp*(efbmb(jbe+1,lsbeam)-efbmb(jbe,lsbeam))
     zdx = dxitk*0.5_fp

     !  energy and x variation @ past half bounce orbit topology

     call eintrp(zev,ie,zface) ! (energy may have changed)

     ! db(e,x-dx)
     call xintrp(xion-zdx,ix,zfacx)
     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)

     ! db(e,x+dx)
     call xintrp(xion+zdx,ix,zfacx)
     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)

     call xintrp(xion,ix,zfacx) ! (back to x)

     ! db(e-de,x)
     call eintrp(zev-zde,ie,zface)
     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)
     ! db(e+de,x)
     call eintrp(zev+zde,ie,zface)
     call get_fp(zdifb,zdum)
     zdifb_max = max(zdifb_max,zdifb)

     call eintrp(zev,ie,zface) ! (back to e)

     ! topology variations: @ x, x-dx, x+dx

     call getdmax(zdifb)      ! @ x
     zdifb_max = max(zdifb_max,zdifb)

     call xintrp(xion-zdx,ix,zfacx)
     call getdmax(zdifb)      ! @ x-dx
     zdifb_max = max(zdifb_max,zdifb)

     call xintrp(xion+zdx,ix,zfacx)
     call getdmax(zdifb)      ! @ x+dx
     zdifb_max = max(zdifb_max,zdifb)

     !-----------------------------------------
     ! ok: zdsdfb sum for output:

     zdlasq = 2.0_fp*dtgfp*zdifb_max
     zdsdfb = zdsdfb + zdlasq*zgxi2
  end if
  return

contains

  subroutine eintrp(zengy,ieloc,zefac)
    real(fp), intent(in) :: zengy
    integer, intent(out) :: ieloc
    real(fp), intent(out) :: zefac

    ! zengy -> (ieloc) element in edifbe_2d(1:med2d)
    !   1 <= ieloc < med2d
    !      zengy = edifbe_2d(ieloc) +
    !        zefac*(edifbe_2d(ieloc+1)-edifbe_2d(ieloc))

    ! range check applied to zengy

    !-------------
    real(fp) :: ze
    !-------------

    ze = max(edifbe_2d(1),min(edifbe_2d(med2d),zengy))

    ieloc = 1 + (ze-edifbe_2d(1))*(med2d-1)/ &
         (edifbe_2d(med2d)-edifbe_2d(1))
    if(ieloc.eq.med2d) ieloc=ieloc-1

    do
       if(edifbe_2d(ieloc).gt.ze) then
          ieloc=ieloc-1
          cycle
       else if(edifbe_2d(ieloc+1).lt.ze) then
          ieloc=ieloc+1
          cycle
       else
          exit
       end if
    end do

    zefac = (ze-edifbe_2d(ieloc))/ &
         (edifbe_2d(ieloc+1)-edifbe_2d(ieloc))
    return
  end subroutine eintrp

  subroutine xeintrp(zxin,zengy,ixloc,ieloc,zxfac,zefac,iwarn)
    real(fp), intent(in) :: zxin,zengy
    integer, intent(out) :: ixloc,ieloc
    real(fp), intent(out) :: zxfac,zefac
    integer, intent(out) :: iwarn

    !-------------
    real(fp) :: ze,zte
    !-------------

    ixloc = 0
    ieloc = 0
    zxfac = 0.0_fp
    zefac = 0.0_fp

    iwarn = 0
    call xintrp(zxin,ixloc,zxfac)

    zte = te(ixloc,1) + zxfac*(te(ixloc+1,1)-te(ixloc,1))

    if(maxval(edifbe_2d) <= 1.0e-6_fp) then
       iwarn = 1
       return
    end if

    ze = max(edifbe_2d(1),min(edifbe_2d(med2d),1.0e+3_fp*zengy/zte))

    ieloc = 1 + (ze-edifbe_2d(1))*(med2d-1)/ &
         (edifbe_2d(med2d)-edifbe_2d(1))
    if(ieloc.eq.med2d) ieloc=ieloc-1

    do
       if(edifbe_2d(ieloc).gt.ze) then
          ieloc=ieloc-1
          cycle
       else if(edifbe_2d(ieloc+1).lt.ze) then
          ieloc=ieloc+1
          cycle
       else
          exit
       end if
    end do

    zefac = (ze-edifbe_2d(ieloc))/ &
         (edifbe_2d(ieloc+1)-edifbe_2d(ieloc))
    return
  end subroutine xeintrp


  subroutine xintrp(zxin,ixloc,zxfac)
    real(fp), intent(in) :: zxin
    integer, intent(out) :: ixloc
    real(fp), intent(out) :: zxfac

    ! zxin -> (ixloc) element in xizb(lcentr:lep1)
    !   lcentr <= ixloc < lep1
    !      zxin = xizb(ixloc) + zxfac*(xizb(ixloc+1)-xizb(ixloc))

    ! range check applied to zxin

    !-------------
    real(fp) :: zx
    !-------------

    zx = max(xizb(lcentr),min(xizb(lep1),zxin))

    ixloc = lcentr + &
         (zx-xizb(lcentr))*nzones/(xizb(lep1)-xizb(lcentr))
    if(ixloc.eq.lep1) ixloc=ixloc-1

    do
       if(xizb(ixloc).gt.zx) then
          ixloc=ixloc-1
          cycle
       else if(xizb(ixloc+1).lt.zx) then
          ixloc=ixloc+1
          cycle
       else
          exit
       end if
    end do

    zxfac = (zx-xizb(ixloc))/ &
         (xizb(ixloc+1)-xizb(ixloc))
    return
  end subroutine xintrp

  subroutine get_fp(zdif,zddif)
    ! diffusivity db(e,x), current orbit topology
    real(fp), intent(out) :: zdif
    real(fp), intent(out) :: zddif

    zdif = zero

    if(l_fd0) zdif= zdif+ w_fd0*getd1(difbe_2d_fd0)
    if(l_fdb) zdif= zdif+ w_fdb*getd1(difbe_2d_fdb)
    if(l_fdp) zdif= zdif+ w_fdp*getd1(difbe_2d_fdp)
    if(l_fdq) zdif= zdif+ w_fdq*getd1(difbe_2d_fdq)
    if(l_fdr) zdif= zdif+ w_fdr*getd1(difbe_2d_fdr)
    if(l_fds) zdif= zdif+ w_fds*getd1(difbe_2d_fds)

    zddif = zero

    if(l_fd0) zddif= zddif+ w_fd0*getd1d(difbe_2d_fd0)
    if(l_fdb) zddif= zddif+ w_fdb*getd1d(difbe_2d_fdb)
    if(l_fdp) zddif= zddif+ w_fdp*getd1d(difbe_2d_fdp)
    if(l_fdq) zddif= zddif+ w_fdq*getd1d(difbe_2d_fdq)
    if(l_fdr) zddif= zddif+ w_fdr*getd1d(difbe_2d_fdr)
    if(l_fds) zddif= zddif+ w_fds*getd1d(difbe_2d_fds)

  end subroutine get_fp

  subroutine getdmax(zdif)
    ! max db(e,x), nearby orbit topologies
    real(fp), intent(out) :: zdif

    !-----------------------
    ! local:

    real(fp) :: zd_fd0,zd_fdb,zd_fdp,zd_fdq,zd_fdr,zd_fds
    real(fp) :: zdif1
    integer :: iwx

    !-----------------------
    !  base diffusivities as needed...

    if(lx_fd0) then
       zd_fd0 = getd1(difbe_2d_fd0)
    else
       zd_fd0 = zero
    end if

    if(lx_fdb) then
       zd_fdb = getd1(difbe_2d_fdb)
    else
       zd_fdb = zero
    end if

    if(lx_fdp) then
       zd_fdp = getd1(difbe_2d_fdp)
    else
       zd_fdp = zero
    end if

    if(lx_fdq) then
       zd_fdq = getd1(difbe_2d_fdq)
    else
       zd_fdq = zero
    end if

    if(lx_fdr) then
       zd_fdr = getd1(difbe_2d_fdr)
    else
       zd_fdr = zero
    end if

    if(lx_fds) then
       zd_fds = getd1(difbe_2d_fds)
    else
       zd_fds = zero
    end if

    zdif = zero

    !  max of interpolations

    do iwx=1,nwx_fd
       zdif1 = wx_fd0(iwx)*zd_fd0
       zdif1 = zdif1 + wx_fdb(iwx)*zd_fdb
       zdif1 = zdif1 + wx_fdp(iwx)*zd_fdp
       zdif1 = zdif1 + wx_fdq(iwx)*zd_fdq
       zdif1 = zdif1 + wx_fdr(iwx)*zd_fdr
       zdif1 = zdif1 + wx_fds(iwx)*zd_fds

       zdif = max(zdif,zdif1)
    end do

  end subroutine getdmax

  real(fp) function getd1(zd2d)

    ! bilinear interpolation in (e,x) space
    real(fp), dimension(:,:,:) :: zd2d

    getd1 = &
         (one-zface)* &
         ((one-zfacx)*zd2d(ix,ie,lsbeam) + &
         zfacx*zd2d(ix+1,ie,lsbeam)) + &
         zface* &
         ((one-zfacx)*zd2d(ix,ie+1,lsbeam) + &
         zfacx*zd2d(ix+1,ie+1,lsbeam))

  end function getd1

  real(fp) function getd1d(zd2d)

    !  finite difference derivative in radial coord. x, interpolated
    !  in e space

    real(fp), dimension(:,:,:) :: zd2d
    real(fp), parameter :: zeps=1.0e-12_fp

    if((ix.eq.ledge).and.(zfacx.gt.one-zeps)) then
       getd1d=zero  ! model has no spatial variation beyond plasma bdy
       return
    end if

    getd1d = ( &
         (one-zface)*(zd2d(ix+1,ie,lsbeam)-zd2d(ix,ie,lsbeam)) + &
         zface*(zd2d(ix+1,ie+1,lsbeam)-zd2d(ix,ie+1,lsbeam)) )/ &
         dxizc(ix)

  end function getd1d

end subroutine abdifu
