!-----------------------------------------------------------------------
!  nb0jxb - jxb momentum/pwr exchange, neutral track <--> ion gyro
!
subroutine nb0jxb(zwght,ze0,ze1,zphi0,zphi1,zbt0,zbt1, &
     iz0,iz1,zx0,zx1,iflag,ilntr)
  !
  !  input arguments:  zwght = weight of monte carlo particle
  !    (ze0,zphi0)  energy, momentum of ptcl in state 0 (orig.)
  !    (ze1,zphi1)  energy, momentum of ptcl in state 1 (final)
  !
  !    iz0,iz1      orig., final zzone indices of ptcl
  !    zx0,zx1      orig., final xi (flux surface) of ptcl
  !    zbt0,zbt1    orig., final btoroidal of ptcl
  !
  !  output:  add to sums
  !   bpjxbs(...) - work done by beam ions jxb torquing the plasma
  !         = beam ion energy change due to radial electric field
  !   bphjxb - momentum given to plasma from fast ions
  !   tqjxb(...) - torque profile to thermal plasma
  !
  !   tirtots(...) -- beam radial current
  !   fpixbts(...) -- bt * beam radial current
  !
  !-------------
  !
  use nbi_com
  !
  !-------------
  !  jxb / rotation effect
  !
  !
  !============
  implicit none
  integer :: iz0,iz1,iflag,izmin,izmax,iz,idum
  !============
  real(fp) :: ze0,ze1,zphi0,zphi1,zbt0,zbt1,zx0,zx1,zwght,zde,zdp
  real(fp) :: zomega,zxmin,zxmax,zpfmax,zpfmin,zf,zdf,zwghti,zsign
  real(fp) :: zswght,zfac,zbt,zyp
  !============
  logical :: ib_nbsce                 ! beam index nbsce.gt.0
  logical :: ilntr ! input, if ilntr=.true. -- contiburion from nutrav
  !                          store outputs in different arrays

  !============
  !
  ib_nbsce = (nbsce.gt.0)
  !
  zde=(ze0-ze1)*zwght
  zdp=(zphi0-zphi1)*zwght
  if(ilntr) then
     bpjxbs_ntr(lsbeam)=bpjxbs_ntr(lsbeam)+zde
     bphjxbs_ntr(lsbeam)=bphjxbs_ntr(lsbeam)+zdp
  else
     bpjxbs(lsbeam)=bpjxbs(lsbeam)+zde
     bphjxbs(lsbeam)=bphjxbs(lsbeam)+zdp
  end if
  !
  !
  !  distribute the torque
  izmin=min(iz0,iz1)
  izmax=max(iz0,iz1)
  do
     if(izmin.eq.izmax) then
        if(ilntr) then
           tqjxbs_ntr(izmin,lsbeam)=tqjxbs_ntr(izmin,lsbeam)+ &
                (zphi0-zphi1)*zwght
        else
           tqjxbs(izmin,lsbeam)=tqjxbs(izmin,lsbeam)+(zphi0-zphi1)* &
                zwght
        end if
        !
        if(ib_nbsce) then
           if(ilntr) then
              tqjxbnb_ntr(izmin,nbsce)=tqjxbnb_ntr(izmin,nbsce)+ &
                   (zphi0-zphi1)*zwght
           else
              tqjxb_nb(izmin,nbsce)=tqjxb_nb(izmin,nbsce)+ &
                   (zphi0-zphi1)*zwght
           end if
        end if
        
        !
        if(ilntr) then
           pbephis_ntr(izmin,lsbeam)=pbephis_ntr(izmin,lsbeam)+ &
                (ze0-ze1)*zwght
        else
           pbephis(izmin,lsbeam)=pbephis(izmin,lsbeam)+(ze0-ze1)*zwght
        end if
        zomega=0.5_fp*(omegag(izmin)+omegag(izmin+1))
     else
        !
        !  ptcl has crossed a surface
        !
        !  jxb torque
        !
        zxmin=min(zx0,zx1)
        zxmax=max(zx0,zx1)
        if(zxmax.eq.zxmin) then
           izmax=izmin
           cycle
        end if
        !
        if(izmax.ge.lep1) then
           zpfmax=plflxi(nzones+1)*1.0e-8_fp
           izmax=ledge
           if(izmax.eq.izmin) cycle
        else
           call psi_spline(zxmax,idum,zpfmax)
        end if
        
        call psi_spline(zxmin,idum,zpfmin)
        !
        if(zpfmax.eq.zpfmin) then
           izmax=izmin
           cycle
        end if
     !
        zf=1.0_fp/(zpfmax-zpfmin)
     !
     !  plflxi is psi(x) at surfaces from integration formula used in psi_spline
     !  index shift down 1 from std lcentr:ledge (legacy artifact)
     !
        do iz=izmin,izmax
           zdf=min(zpfmax,1.0e-8_fp*plflxi(iz))- &
                max(zpfmin,1.0e-8_fp*plflxi(iz-1))
           zwghti=zwght*zdf*zf
           if(ilntr) then
              tqjxbs_ntr(iz,lsbeam)=tqjxbs_ntr(iz,lsbeam)+ &
                   (zphi0-zphi1)*zwghti
           else
              tqjxbs(iz,lsbeam)=tqjxbs(iz,lsbeam)+(zphi0-zphi1)*zwghti
           end if
           !
           if(ib_nbsce) then
              if(ilntr) then
                 tqjxbnb_ntr(iz,nbsce)=tqjxbnb_ntr(iz,nbsce)+ &
                      (zphi0-zphi1)*zwghti
              else
                 tqjxb_nb(iz,nbsce)=tqjxb_nb(iz,nbsce)+(zphi0-zphi1)* &
                      zwghti
              end if
           end if
           !
           if(ilntr) then
              pbephis_ntr(iz,lsbeam)=pbephis_ntr(iz,lsbeam)+(ze0-ze1)* &
                   zwghti
           else
              pbephis(iz,lsbeam)=pbephis(iz,lsbeam)+(ze0-ze1)*zwghti
           end if
           
           zomega=0.5_fp*(omegag(iz)+omegag(iz+1))
           
        end do
        !
        !  radial current & bt * radial current sums
        !
        if(iz1.gt.iz0) then
           zsign=1.0_fp
        else
           zsign=-1.0_fp
        end if
        zswght=zsign*zwght
        !
        if(ilntr) then
           do iz=izmin+1,izmax
              tirtots_ntr(iz,lsbeam)=tirtots_ntr(iz,lsbeam)+xzbeami* &
                   zswght
              !
              zfac=(xi(iz,1)-zx0)/(zx1-zx0)
              zbt=zbt0+zfac*(zbt1-zbt0)
              !
              fpixbts_ntr(iz,lsbeam)=fpixbts_ntr(iz,lsbeam)+xzbeami* &
                   zswght*zbt
           end do
        else
           do iz=izmin+1,izmax
              tirtots(iz,lsbeam)=tirtots(iz,lsbeam)+xzbeami*zswght
              !
              zfac=(xi(iz,1)-zx0)/(zx1-zx0)
              zbt=zbt0+zfac*(zbt1-zbt0)
              !
              fpixbts(iz,lsbeam)=fpixbts(iz,lsbeam)+xzbeami*zswght*zbt
              
           end do
        end if
        !
     end if
     exit
  enddo
  !
  zyp=yion
  !
  return
end subroutine nb0jxb
