!  bbdifu -- beam anomalous diffusion operator
!   nmdifb.eq.4 version -- use pitch/energy/space dependent diffusion;
!
!   actual transport computed in abdifu; this routine is called once per
!   half bounce and it saves orbit-topology-dependent interpolation factors
!   which introduce pitch dependence to the diffusion operator: w_f_fp, etc.
!
!   also saved are interpolation factors for "nearby" pitch values--
!   used to probe operator strengh of nearby orbits:  wx_f_fp(...), etc.;
!   nwx_fd weight factor sets saved.
!
!   see also: bbdifu_init (separate subroutine, below) which should be
!   called, to clear bbdifu outputs, when a new orbit track integration
!   is started...
!
subroutine bbdifu(zvplla,ierr)
  !
  use nbi_com
  !
  implicit none
  !
  !  passed arguments (input):
  !
  real(fp) :: zvplla  ! half bounce averaged <vpll> -- only sign is used
  !
  !  also input:
  !  in nbi_com: dtrap ~ measure of trapping depth -- related to vpll/v
  !              dpass ~ measure of passing depth  -- related to vpll/v
  !
  !  these are computed in orbit
  !
  !    they are used to set up interpolation factors, 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
  !    (see abdifu.for)
  !
  !  passed arguments (output):
  !
  integer :: ierr   ! status code, 0=0k
  !
  !--------------------------------
  !  local...
  !
  real(fp) :: zdpa  ! est. pitch angle spread per goosed orbit
  real(fp) :: ztrap,zpass  ! nearby trapped/passing factors
  !
  !
  !---------------------------------
  !
  ierr=0
  !
  if(nmdifb.eq.4 .or. nmdifb.eq.6) then
     goto 10
  else
     return
  end if
  !
  !  pitch weight interpolation factors...
  !
10 continue

  w_fd0=zero
  w_fdb=zero
  w_fdp=zero
  w_fdq=zero
  w_fdr=zero
  w_fds=zero

  if(dtrap.lt.zero) then
     ! passing orbit
     if(zvplla.lt.zero) then
        ! ctr passing
        w_fdr = one-dpass
        w_fds = dpass

     else
        ! co passing
        w_fdp = one-dpass
        w_fdq = dpass

     end if

  else
     ! trapped orbit

     w_fd0 = dtrap
     w_fdb = one-dtrap

  end if
  !
  !  interpolation factors for nearby orbits...
  !
  zdpa = dif2d_dpa
  !
  nwx_fd = 0
  !
  wx_fd0 = zero
  wx_fdb = zero
  wx_fdp = zero
  wx_fdq = zero
  wx_fdr = zero
  wx_fds = zero

  if(dtrap.lt.zero) then
     ! passing orbit
     if(zvplla.lt.zero) then
        ! ctr passing
        if(dpass.lt.0.99_fp) then
           ! nearby orbit, more deeply passing...
           zpass = min(one,dpass+zdpa)
           nwx_fd = nwx_fd + 1
           wx_fds(nwx_fd) = zpass
           wx_fdr(nwx_fd) = one-zpass

        end if

        if(dpass.ge.zdpa) then
           ! nearby orbit, less deeply passing...
           zpass = dpass-zdpa
           nwx_fd = nwx_fd + 1
           wx_fds(nwx_fd) = zpass
           wx_fdr(nwx_fd) = one-zpass

        else
           ! nearby orbit crosses trap/passing bdy...
           nwx_fd = nwx_fd + 1
           wx_fdb(nwx_fd)=one ! barely trapped
           nwx_fd = nwx_fd + 1
           wx_fdp(nwx_fd)=one ! barely passing co
           nwx_fd = nwx_fd + 1
           wx_fdr(nwx_fd)=one ! barely passing ctr

           ztrap = zdpa - dpass
           zpass = zdpa - dpass

           ! nearby trapped orbit
           nwx_fd = nwx_fd + 1
           wx_fd0(nwx_fd) = ztrap
           wx_fdb(nwx_fd) = one-ztrap

           ! nearby co passing orbit
           nwx_fd = nwx_fd + 1
           wx_fdq(nwx_fd) = zpass
           wx_fdp(nwx_fd) = one-zpass

        end if

     else
        ! co passing
        if(dpass.lt.0.99_fp) then
           ! nearby orbit, more deeply passing...
           zpass = min(one,dpass+zdpa)
           nwx_fd = nwx_fd + 1
           wx_fdq(nwx_fd) = zpass
           wx_fdp(nwx_fd) = one-zpass

        end if

        if(dpass.ge.zdpa) then
           ! nearby orbit, less deeply passing...
           zpass = dpass-zdpa
           nwx_fd = nwx_fd + 1
           wx_fdq(nwx_fd) = zpass
           wx_fdp(nwx_fd) = one-zpass

        else
           ! nearby orbit crosses trap/passing bdy...
           nwx_fd = nwx_fd + 1
           wx_fdb(nwx_fd)=one   ! barely trapped
           nwx_fd = nwx_fd + 1
           wx_fdp(nwx_fd)=one   ! barely passing co
           nwx_fd = nwx_fd + 1
           wx_fdr(nwx_fd)=one   ! barely passing ctr

           ztrap = zdpa - dpass
           zpass = zdpa - dpass

           ! nearby trapped orbit
           nwx_fd = nwx_fd + 1
           wx_fd0(nwx_fd) = ztrap
           wx_fdb(nwx_fd) = one-ztrap

           ! nearby ctr passing orbit
           nwx_fd = nwx_fd + 1
           wx_fds(nwx_fd) = zpass
           wx_fdr(nwx_fd) = one-zpass

        end if

     end if

  else
     ! trapped orbit
     if(dtrap.lt.0.99_fp) then
        ! nearby orbit, more deeply trapped...
        ztrap = min(one,dtrap+zdpa)
        nwx_fd = nwx_fd + 1
        wx_fd0(nwx_fd) = ztrap
        wx_fdb(nwx_fd) = one-ztrap

     end if

     if(dtrap.ge.zdpa) then
        ! nearby orbit, less deeply trapped...
        ztrap = dtrap-zdpa
        nwx_fd = nwx_fd + 1
        wx_fd0(nwx_fd) = ztrap
        wx_fdb(nwx_fd) = one-ztrap

     else
        ! nearby orbit crosses trap/passing bdy...
        nwx_fd = nwx_fd + 1
        wx_fdb(nwx_fd)=one  ! barely trapped
        nwx_fd = nwx_fd + 1
        wx_fdp(nwx_fd)=one  ! barely passing co
        nwx_fd = nwx_fd + 1
        wx_fdr(nwx_fd)=one  ! barely passing ctr

        zpass = zdpa - dtrap

        ! nearby ctr passing orbit
        nwx_fd = nwx_fd + 1
        wx_fds(nwx_fd) = zpass
        wx_fdr(nwx_fd) = one-zpass

        ! nearby co passing orbit
        nwx_fd = nwx_fd + 1
        wx_fdq(nwx_fd) = zpass
        wx_fdp(nwx_fd) = one-zpass

     end if

  end if

  lx_fd0 = (maxval(wx_fd0(1:nwx_fd)).gt.zero)
  lx_fdb = (maxval(wx_fdb(1:nwx_fd)).gt.zero)
  lx_fdp = (maxval(wx_fdp(1:nwx_fd)).gt.zero)
  lx_fdq = (maxval(wx_fdq(1:nwx_fd)).gt.zero)
  lx_fdr = (maxval(wx_fdr(1:nwx_fd)).gt.zero)
  lx_fds = (maxval(wx_fds(1:nwx_fd)).gt.zero)

  l_fd0 = (w_fd0.gt.zero)
  l_fdb = (w_fdb.gt.zero)
  l_fdp = (w_fdp.gt.zero)
  l_fdq = (w_fdq.gt.zero)
  l_fdr = (w_fdr.gt.zero)
  l_fds = (w_fds.gt.zero)

  return
end subroutine bbdifu

!---------------------------------------------

subroutine bbdifu_init(ierr)
  !
  !  initialize bbdifu controls.
  !   on first half bounce, no diffusion, but for all pitch angles the
  !   diffusivities are probed, for goosing control
  !
  use nbi_com
  !
  implicit none
  !
  !  passed arguments (output):
  !
  integer :: ierr   ! status code, 0=0k
  !
  !--------------------------------
  !  local...
  !
  !--------------------------------

  ierr = 0

  w_fd0=zero
  w_fdb=zero
  w_fdp=zero
  w_fdq=zero
  w_fdr=zero
  w_fds=zero

  l_fd0 = .false.
  l_fdb = .false.
  l_fdp = .false.
  l_fdq = .false.
  l_fdr = .false.
  l_fds = .false.

  wx_fd0 = zero
  wx_fdb = zero
  wx_fdp = zero
  wx_fdq = zero
  wx_fdr = zero
  wx_fds = zero

  !  probe each diffusivity; the maximum diffusivity will
  !  determine goosing after the 1st complete half bounce; after 2nd
  !  and subsequent half bounces, bbdifu results will apply

  nwx_fd = 6

  wx_fd0(1) = one  ! f_fp on 1
  wx_fdb(2) = one  ! fdb on 2
  wx_fdp(3) = one  ! ...etc...
  wx_fdq(4) = one
  wx_fdr(5) = one
  wx_fds(6) = one

  lx_fd0 = .true.
  lx_fdb = .true.
  lx_fdp = .true.
  lx_fdq = .true.
  lx_fdr = .true.
  lx_fds = .true.

end subroutine bbdifu_init
