!******************** start file nutrav.fnb ; group nubeam *************
!*************************************************************


subroutine nutrav(icxo,n_flr,zcospa,zsinpa,zwghts,ilouta,jsplit, &
     ireto,zrs,zys,zxs,zepoti,zphii,ierr)
  !
  !  subroutine beamcx calls this routine
  !
  !  dmc 13 may 1993 -- adding call to excited states model interface
  !
  !---------------------------------------------------------------------
  !
  !  calling arguments
  !    jsplit - follow jsplit neutral tracks
  !      all start from guiding ctr (rmjion,yion)
  !      with finite flr offset
  !    zcospa,zsinpa(1...jsplit) initial gyro phase for each of the jsplit tracks
  !    zwghts(1...jsplit) montecarlo weight for each of the test neutrals
  !    ilouta(1...jsplit) .true. if track starts outside of plasma
  !
  !    *removed* zeovab ev/nucleon of test neutrals
  !      computed locally
  !
  !    ireto-- flag whether ch-ex ptcl goes back to orbit immediately
  !    or later
  !
  !       =====================================================
  !       this subroutine moves neutrals created from beam ions
  !       across the plasma column
  !       =====================================================

  !
  !  revised d. mc cune 18 jan 1982
  !
  !   add ####nb.out output-- position and velocity coordinates,
  !   energy charge and weight of escaping beam m.c. neutrals.
  !    pos. coords-- radius, angle
  !    vel. coords-- dr/dl, dy/dl in poloidal plane of escape, dz/dl
  !     at right angles to pol. plane in direction of plasma current
  !     (xksid, vparallel/v)
  !   the transformation from the coordinates (dx/dl, dy/dl, dz(old)/dl)
  !   in the poloidal plane of origin of the neutral is:
  !    dr/dl= 1/r*[(rmajor+x)*dx/dl + z*dz(old)/dl]
  !    dy/dl doesnt change
  !   dz/dl=1/r*[-z*dx/dl +(rmajor+x)*dz(old)/dl]
  !    where r= distance of escaping neutral from machine centerline
  !      when it escapes,
  !    z=l*dz(old)/dl, x=l*dx/dl, l=pathlength from pt. of origin to
  !      point of escape, rmajor= plasma major radius
  !  replaced old definition using rt (tangency radius)
  !
  !      installed flr effects  -- rob goldston dec. 26, 1982
  !
  !  dmc february 1987
  !   rotation upgrade
  !   atomic collision interactions with thermal plasma in plasma frame
  !   momentum balance accounting
  !   flr redone
  !
  !      upgrade to allow for a limiter system (defined in limlos.for)
  !      more flexible than the previous square box.  to make
  !      communication with limlos simplest, have used rmjion
  !      in place of previous zrmaj
  !      rob goldston july 21, 1983
  !
  !  modif dmc 12 aug 83:  deleted re-initialization of "ttcx" (time
  !  to next cx)-- this is now done in beamcx.for
  !
  !  rev dmc dec 1984 *** cross section section replaced ***
  !
  !    use lookup tables bnsigv, bbnsig, precalculated in
  !    subroutine nbatom ****
  !
  !-------------------------------------------------------------------
  !  geometry upgrade dmc apr 1985 ***
  !
  !  coding reorganized along the lines of dep:
  !    call to subroutine cxline evaluates flr starting parameters;
  !      (similar to nbline call in dep)
  !    call to trkrb evaluates all zone intersections
  !
  !    ptcl source terms are evaluated all along the neutral trajectory;
  !      a monte carlo choice is made to either deposit the ptcl somewhere
  !      along this trajectory or allow it to leave the system
  !
  !  common output - parameters of recaptured orbit
  !    xi,theta  vion,xksid  time=bitime(1)  wght
  !  either returned directly to orbit or indirectly via splitting storage
  !  arrays thbsa, etc.
  !
  !
  !  tbt oct 16, 1991  nbatom upgrade. assume nlvphi=t.  tbt.
  !
  !  --------------------------------------------------------------------
  !
  use nbi_com
  use nbi_mpi_com
  use nbi_random_pll
  use nbatom_mod
  use nbi_strline,only : nbi_track_call
  use logmod
  !%      use sglib_mod
  !
  !============
  implicit none
  integer, intent(in) :: icxo ! if zero -- use only for cx statistics improvement
  integer :: jsplit,ireto,ingc,idum,js,is,irng,isb,i,j,iz,indi
  integer :: indp1,ig,itrmax,inoflr,jbsplm,ierr
  !============
  real(fp) :: zsinpa,zwghts,zrs,zys,zxs,zepoti,zphii,zcospa
  real(fp) :: zrmjio,zyion,zvion,zksid,zvperp,zvpll,zxion
  real(fp) :: zthion,zxista,zxo,ztho,zvo,zkso,zwghto,zbmod,zbtgci
  real(fp) :: zepot0,zphi0,zinc,zg,zbt0,zthea,zphia,zexp,znsigl
  real(fp) :: zlvac,zlp,zlav,zlava,zrloc,zyloc,zvtor,zvr,zvy, &
       zeovab
  real(fp) :: zintri,zintrp,zexcfac,zbnsvil,zbnsvis,zbbnsvi,znsgth
  real(fp) :: zbmsexc1d, zsigvnt1d
  real(fp) :: zbmsexc2d, zsigvnt2d
  real(fp) :: zsigvnt
  real(fp) :: znsgbb,zmfpi,zfth,zfbb,zdnsl,zexpol,zdexp,zexpdl,zdn0
  real(fp) :: zdv0,zde0,zdsrat,zomega,zescap,zrand,zprbls,zlndep,zrr
  real(fp) :: zwghtp,zphi00,zphigc,ze00,zegc,zbtgc
  real(fp) :: zsinth,zcosth,zsinph,zcosph,zxhor,zyhor
  integer, dimension(jsplit) :: n_flr
  !============
  dimension zwghts(jsplit),zcospa(jsplit),zsinpa(jsplit)
  !  following used to compare fast flr and trkrb flr displacements
  !  for new code debug purposes:
  dimension zxs(jsplit),zrs(jsplit),zys(jsplit) ! test
  !
  real(fp), dimension(:), allocatable :: zxsa,zrsa,zysa  ! debug
  !
  logical ilouta(jsplit)
  logical :: ilbms1d !ture - beam-beam rections are included and for adas beam
  !stopping coefficient and enhancement factor calculation
  !turn on
  real(fp) :: zth,zxi,z_ksid,alpha_v
  !
  !%      real(fp) :: zzx(8*mj),zzy(8*mj),zzx2(3),zzy2(3)
  !%      real(fp) :: zlpl,zlplo,zlplot,zlmid,zl1,zleps
  !%      real(fp) :: zrplo,zyplo,zdummy(2,2),zzth,zzxi
  !%      integer :: ilintv,ipl,ier
  !
  !
  !-----------------------/ modified dmc 12/1984 - see dep.for, nbatom.for
  !----------------------/ now using common arrays bnsigv(i,j,k), bbnsig()
  !

  real(fp) :: zbnsvcx(mig)    ! interpolated cx cross-section data.
  real(fp) :: zbnsvii           ! interpolated ii cross-section data.
  real(fp) :: zbnsvie         ! interpolated elec i cross section data
  real(fp) :: zbnsviz         ! interpolated impurity i cross sec data
  real(fp) :: zbnsvto         ! interpolated total cross-section data.
  real(fp) :: zbbfi         ! interpolated beam-beam ii/tot x-sect data
  real(fp) :: zbbnsv         ! interpolated beam-beam total  x-sect data
  !
  real(fp) :: zbbnsva(2,mibs)
  integer :: irng2(mibs)
  !
  !-------------------------------------
  !  excited states model interface working arrays
  !
  real(fp) :: zwkthcx(2,mig)  ! th. cx w and w/o excitation effect
  real(fp) :: zwkbb(2,2,mibs) ! bb cx and ii w and w/o excitation effect
  real(fp) :: zwkthiz(2,3)    ! th. ionizations w, w/o excitation effect
  real(fp) :: zbbxfac(mibs)   ! excited states effect on beam-beam cx
  real(fp) :: zwkfac(7)       !  output excited states enhancements
  !                   1 -> electron impact ionization
  !                   2 -> impurity impact ionization
  !                   3 -> th. plasma species impact ionization
  !                   4 -> th. plasma species charge exchange
  !                   5 -> beam beam impact ionization
  !                   6 -> beam beam charge exchange
  !                   7 -> net effect on stopping cross section
  !
  !-------------------------------------
  !
  !  atomic physics:
  !  see subroutines nbatom (set-up) and nbains and nbain2 (interpolation)
  !
  real(fp) :: zbdir(4)  ! guiding ctr field data
  !
  !  track segment deposition probability, integrated distance
  !
  real(fp) :: zprb(8*mj),zlng(8*mj)
  integer :: i_zone,i_dir
  integer :: ib
  real(fp) ::  zdvrho,zdvth !velocity components -- flux surface normal direction
  !and flux suface tangential direction in poloidal plane
  real(fp) :: zbdir_nrm,zangl
  real(fp) :: zvec(3) !vr,vphi,vz velocity of the particle (work array here)
  !
  !
  !-----------------------------------------------
  !

  ierr=0
  ilbms1d=(nlbbcx.and.(nsigexc_use.eq.1).and. &
       (lev_nbidep.eq.2) &
       .and.(nbbcx_bb.eq.0))
  !
  !%      zzx=0; zzy=0; zzx2=0; zzy2=0
  !%      zzx(1)=100.0_fp
  !%      zzy(1)=100.0_fp
  !%      zzx2(1)=1.0_fp
  !%      zzy2(1)=1.0_fp
  !
  tscx=0._fp
  nbicx=0
  !
  !    follow the neutrals
  !    loop through splittees prepared in subroutine cxsr
  !     1st ptcl returned to orbit, subsequent ptcls saved in
  !      splitting arrays for later reinsertion into orbit
  !
  !  initial parameters, all neutral tracks
  !
  zrmjio=rmjion
  zyion=yion
  !
  zvion=vion
  zksid=xksid
  zvperp=vperp
  zvpll=vpll
  !
  ingc=ngc
  zxion=xion
  zthion=th
  !
  !  if saving current ion save its parameters
  if(ireto.eq.1) then
     zxista=xistate
     zxo=xion
     ztho=th
     zvo=vion
     zkso=xksid
     zwghto=wght
  end if
  !
  !  get guiding ctr field data
  !
  call nblocb(rmjion,yion,xion,th,zbmod,zbdir,2,nfast_map,idum)
  zbtgci=1.e-4_fp*zbmod*zbdir(1)
  !
  !-------------------------------------------
  !  follow all neutrals
  !
  if(jsplit.le.0) goto 1001
  ! allocate arrays
  call alloc_sa_arrays
  do js=1,jsplit
     !
     !%      call line(0)
     !%      call symbl(0)
     !
     wght=zwghts(js)

     sfcx0mc(lsbeam)=sfcx0mc(lsbeam)+wght

     sfcxmc1d(ngc,lsbeam)=sfcxmc1d(ngc,lsbeam)+wght
     !        sfcxmc1d(n_flr(js),lsbeam)=sfcxmc1d(n_flr(js),lsbeam)+wght

     !
     vion=zvion
     xksid=zksid
     vperp=zvperp
     vpll=zvpll
     !
     rmjion=zrmjio
     yion=zyion
     !
     ngc=ingc
     xion=zxion
     !
     if(ilouta(js)) then
        bpcxx0s(lsbeam)=bpcxx0s(lsbeam)+vion*vion*v2toev*wght
     else
        bpcxi0s(lsbeam)=bpcxi0s(lsbeam)+vion*vion*v2toev*wght
     end if
     !
     !
     !  initial track parameters from flr gyro - pass phase angle
     !
     xion=zxion
     th=zthion
     call cxline(zcospa(js),zsinpa(js),zrs(js),zys(js),zxs(js), &
          zepoti,zphii,zbdir,zepot0,zphi0)
     !
     zrsa(js)=rsta
     zysa(js)=ysta
     zxsa(js)=rosta
     !
     !  bt at start point
     if(rosta.le.1.0_fp) then
        zinc=nzones*(rosta-xi(n_flr(js),1))
        zinc=max(0.0_fp,min(1.0_fp,zinc))
     else
        zinc=0.0_fp
     end if
     zg=fbz(n_flr(js),1)+zinc*(fbz(n_flr(js)+1,1)-fbz(n_flr(js),1))
     zbt0=bzxr*zg/rsta
     !
     !  gc -> flr; zepotj arguments reversed because they are potential
     !  rather than kinetic energies --
     call nb0jxb(wght,zepot0,zepoti,zphii,zphi0,zbtgci,zbt0, &
          ngc,n_flr(js),xion,xiflr,2,.true.)
     !
     bphcxls(lsbeam)=bphcxls(lsbeam)+wght*zphi0
     !
     !  follow the track-- reionization only in the plasma
     !
     call nbi_shrink
     !
     !  straln is gone...        call trkrl
     !
     call nbi_track_call(rsta, ysta, phia, thea, ierr)
     if(ierr.ne.0) return
     !
     !  restore expanded plasma zones
     !
     call nbi_rbloat
     !
     !  check if track misses plasma entirely
     !
     if(nzcrtk.eq.0) then
        !
        sfcxesc(lsbeam)=sfcxesc(lsbeam)+wght
        if(ilouta(js)) then
           bpcxots(lsbeam)=bpcxots(lsbeam)+vion*vion*v2toev*wght
           bph0was(lsbeam)=bph0was(lsbeam)+zphi0*wght
        else
           if(zxs(js).lt.0.950_fp) then
              call errorlog &
                   ('nutrav: interior track missed guiding center')
              call errorlog(' rmjion=',rmjion)
              call errorlog(' yion=',yion)
              call errorlog(' rlar=',rlar)
              call errorlog(' xion=',xion)
              if(myidd.eq.0) then
                 write(nout,8801) rmjion,yion,rlar,xion
                 do is=1,js
                    write(nout,8802) is,zrs(is),zys(is),zxs(is), &
                         zrsa(is),zysa(is),zxsa(is)

                 end do
              end if
8801          format( &
                   ' ? nutrav - interior track missed.  guiding center:'/ &
                   '   rmjion=',1pe11.4,' yion=',1pe11.4,' rlar=',1pe11.4,' xion=', &
                   1pe11.4/ &
                   ' no.  fast flr (r,y,xi)                 trkrb flr (r,y,xi)')
8802          format(1x,i2,6(1x,1pe12.5))
              !
              ierr=1
              return
              !
           else
              !
              !     close call-- treat as lost track and procede
              !
              bpcxins(lsbeam)=bpcxins(lsbeam)+vion*vion*v2toev* &
                   wght
              bpcxers(lsbeam)=bpcxers(lsbeam)+vion*vion*v2toev* &
                   wght
              bph0was(lsbeam)=bph0was(lsbeam)+zphi0*wght
           end if
        end if
        !
        goto 800
     else
        !     track enters plasma - sc pt.
        !
        if(.not.ilouta(js)) then
           sbcx0s(nznttk(1),lsbeam)=sbcx0s(nznttk(1),lsbeam)+wght
        end if
     end if
     !
     !  trajectory angle info
     !
     !
     zthea=degrad*thea
     zphia=degrad*phia
     !
     zsinth=sin(zthea)
     zcosth=cos(zthea)
     !
     zsinph=sin(zphia)
     zcosph=cos(zphia)
     !
     !  track through plasma; compute deposition profiles by plasma zone;
     !  compure probability of recapture in each zone
     !
     irng = 1     ! 1st guess at range # of energy grid.
     do isb=1,nsbeam
        irng2(isb)=1
     end do
     !
     zexp=1.0_fp
     znsigl=0.0_fp
     !
     !  vacuum segment length
     zlvac=xltkbf
     !
     !  loop over track segments
     !
     do i=1,nzcrtk
        !  transp zone index
        j=nznttk(i)
        !  beam zone index
        iz=nznbtk(i)
        !
        !  distance tracked
        !
        if(i.eq.1) then
           zlp=0.0_fp
           zlng(1)=xlngtk(1)
        else
           zlp=zlng(i-1)
           zlng(i)=zlp+xlngtk(i)
        end if
        !
        !  distrance tracked thru plasma - to zone midpt.
        zlav=0.5_fp*(zlng(i)+zlp)
        !
        !  distance from starting pt. (incl. vacuum track sections)
        !   check for crossing 2nd vacuum section (beamlet exits and reenters)
        if(i.eq.(nztkmd+1).and.nztkmd.ne.0) then
           zlvac=zlvac+xltkmd
        end if
        zlava=zlav+zlvac
        !
        !  look up plasma interaction rates *now* if plasma is rotating
        !  (in which case neutral energy in plasma frame is a fcn of position)
        !
        !  fetch position and velocity in (r,y,phi) coordinate system
        call nb0loc(zlava,zsinth,zcosth,zsinph,zcosph, &
             zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)

        !  mg july 27 2011 -- spend some cpu time -- get actual b vector
        !
        call nblocb(zrloc,zyloc,zxi,zth,zbmod,zbdir,2,nfast_map,idum)
        !
        if(idum.gt.1) then
           call errorlog &
                ('nutrav nblocb: unexpected (r,z)->(xi,theta) map error!')
           if(myidd.eq.0) then
              write(nonlin,*) &
                   ' ?? nutrav nblocb: unexpected (r,z)->(xi,theta) map error!'
           end if
           ierr=1
           return
        end if
        !
        !        .fetch relative energy in plasma frame
        call nberel(zeovab,j,zrloc,zvtor,zvr,zvy)
        !        .fetch interpolation factors
        call nbeloc(zeovab,irng,indi,zintri,indp1,zintrp)
        !        .interpolate table of plasma interactions
        call nbains( &
             zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zexcfac, &
             zbmsexc1d, zsigvnt1d,j,indi,zintri,indp1,zintrp)

        if(ilbms1d) then    !include beam-beam reactions
           ! recalculate zexcfac
           zexcfac = 1._fp
           if(nlfprod(lsbeam)) then
              if(zsigvnt1d.gt.0.0_fp) then
                 zexcfac = zbmsexc1d/zsigvnt1d
              end if
           else
              if(nbbcx_avg.eq.0) then
                 i_zone=iz
              else
                 i_zone=j
              end if
              zbmsexc2d=bms_exc_2d(i_zone,nbien,nbsce)
              zsigvnt2d=sigvnt_2d(i_zone,nbien,nbsce)
              if(zbmsexc2d.eq.0.0_fp) then
                 if(zsigvnt1d.gt.0.0_fp) then
                    zexcfac = zbmsexc1d/zsigvnt1d
                 end if
              else
                 zsigvnt = zsigvnt1d + zsigvnt2d
                 if(zsigvnt.gt.0.0_fp) then
                    zexcfac = (zbmsexc1d+zbmsexc2d)/zsigvnt
                 end if
              end if
           end if
        end if

        !     interpolate table of beam-beam interactions
        call nbain2(zbbnsva,iz,zvtor,zvr,zvy,irng2)
        !
        !--------------------------------------
        !
        !  excited states model option -- adjust n*sigma*v data
        !    zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zbbnsva -- used in
        !      deposition calculation
        !      zwkfac -- model enhancement factors
        !      for output
        !
        zbnsvil=0.0_fp
        if(nsigexc_use.gt.0) then
           !
           call exc_sig(zexcfac,zbbxfac,zwkfac, &
                zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto, &
                zbbnsva,ierr)
           if(ierr.ne.0) return
           !
        end if
        !
        !--------------------------------------
        !
        !  thermal species impact ionization fraction
        zbnsvis=zbnsvii+zbnsvie+zbnsviz+zbnsvil
        !  beam species impact ionization fraction
        zbbnsv=0.0_fp
        zbbnsvi=0.0_fp
        do isb=1,nsbeam
           zbbnsv=zbbnsv+zbbnsva(1,isb)+zbbnsva(2,isb)
           zbbnsvi=zbbnsvi+zbbnsva(1,isb)
        end do
        if(zbbnsv.gt.0.0_fp) then
           zbbfi=zbbnsvi/zbbnsv
        else
           zbbfi=0.0_fp
        end if
        !
        !  local inverse mean free path
        znsgth=zbnsvto/vion
        znsgbb=zbbnsv/vion
        zmfpi=znsgth+znsgbb
        !  fractional accounting
        zfth=znsgth/zmfpi
        zfbb=znsgbb/zmfpi
        !  probability integral dl/mfp
        zdnsl=xlngtk(i)*zmfpi
        zexpol=zexp
        znsigl=znsigl+zdnsl
        !  dmc 11 aug 1991-- protect against exp(-large) which produces a
        !  bogus floating overflow on the decstation:
        if(znsigl.le.80.0_fp) then
           zexp=exp(-znsigl)
        else
           zexp=0.0_fp
        end if
        zdexp=zexpol-zexp
        zexpdl=zdexp/zmfpi
        !
        if(i.eq.1) then
           zprb(1)=zdexp
        else
           zprb(i)=zprb(i-1)+zdexp
        end if
        !  update neut. density and energy
        ! as an alternative may use btiltc(4,5) instead of zbdir(3,4)
        !  btiltc(4,iz) is r component of poloidal field (proportional to dr/dth)
        !  btiltc(5,iz) iz z component of poloidal field (proportional to dz/dth)

        zbdir_nrm=zbdir(4)*zbdir(4)+zbdir(3)*zbdir(3) !must be 1
        zdn0=wght*zexpdl/vion
        !see cxline.for (l81-l82) for formulas
        zdv0=wght*zexpdl*zvtor/vion
        if(zbdir_nrm.gt.0.0_fp) then
           zdvrho=(zbdir(4)*zvr-zbdir(3)*zvy)/zbdir_nrm ! v_rho of cx beam neutrals,
           ! flux surface normal direction
           !
           zdvth=(zbdir(3)*zvr+zbdir(4)*zvy)/zbdir_nrm ! v_theta of cx beam neutrals,
           ! flux suface tangential direction in poloidal plane
        else
           zdvrho=0.0_fp
           zdvth=0.0_fp
        end if
        !
        z_ksid=(zbdir(1)*zvtor-nsnccw*zbdir(2)*zdvth)/vion !pitch angle at present location
        zangl=atan2(zsinpa(js),zcospa(js))
        if(.not.nlfbmflr) then
           alpha_v=0.0_fp
        else
           call nbi_aimingangle &
                (rmjion,yion,zangl,zvpll,zvperp,zvec,alpha_v)
        end if
        zdvrho=wght*zexpdl*zdvrho/vion
        zdvth=wght*zexpdl*zdvth/vion
        zde0=wght*zexpdl*vion
        !  1d arrays
        bn0cx(j,lsbeam)=bn0cx(j,lsbeam)+zdn0
        bv0cx(j,lsbeam)=bv0cx(j,lsbeam)+zdv0
        be0cx(j,lsbeam)=be0cx(j,lsbeam)+zde0
        !  2d arrays
        bn0cx2(iz,lsbeam)=bn0cx2(iz,lsbeam)+zdn0
        bv0cx2(iz,lsbeam)=bv0cx2(iz,lsbeam)+zdv0
        !     velocity normal to flux surface
        bvrhocx2(iz,lsbeam)=bvrhocx2(iz,lsbeam)+zdvrho
        !     velocity tangent to flux surface in poloidal direction
        bvthcx2(iz,lsbeam)=bvthcx2(iz,lsbeam)+zdvth
        bvpllcx2(iz,lsbeam)=bvpllcx2(iz,lsbeam)+ &
             wght*zexpdl*z_ksid/vion
        balpha_gyro(iz,lsbeam)=balpha_gyro(iz,lsbeam)+wght*zexpdl* &
             alpha_v/vion
        be0cx2(iz,lsbeam)=be0cx2(iz,lsbeam)+zde0
        be20cx2(iz,lsbeam)=be20cx2(iz,lsbeam)+zde0*vion*vion
        !
        !  excited states stuff
        if(nsigexc_use.gt.0) then
           do isb=1,nsbeam
              hexbbxrc(iz,lsbeam,isb)=hexbbxrc(iz,lsbeam,isb) &
                   +zdn0*zbbxfac(isb)
           end do
           excitr(j,lsbeam)=excitr(j,lsbeam)+zwkfac(7)*zdn0
        end if
        !
        !  source functions
        sorbes_ntr(j,lsbeam)=sorbes_ntr(j,lsbeam)+wght*zdexp*xzbeami* &
             (zfth*zbnsvis  +zfbb*zbbfi)

        !  nutrav output recapture fcns
        if(.not.ilouta(js)) then
           sbcxris(j,lsbeam)=sbcxris(j,lsbeam)+wght*zdexp*zfth*zbnsvii
           sbcxres(j,lsbeam)=sbcxres(j,lsbeam)+wght*zdexp*zfth*zbnsvie
           sbcxrzs(j,lsbeam)=sbcxrzs(j,lsbeam)+wght*zdexp*zfth*zbnsviz
           sbcxrls(j,lsbeam)=sbcxrls(j,lsbeam)+wght*zdexp*zfth*zbnsvil
           sbcxrbs(j,lsbeam)=sbcxrbs(j,lsbeam)+wght*zdexp*zfbb !cx + ionization
        end if
        !
        sfxrcsc(lsbeam)=sfxrcsc(lsbeam)+wght*zdexp
        sfxrciz(lsbeam)=sfxrciz(lsbeam)+wght*zdexp*zfth*zbnsvis
        sfrbbiz(lsbeam)=sfrbbiz(lsbeam)+wght*zdexp*zfbb*zbbfi
        sfrbbcx(lsbeam)=sfrbbcx(lsbeam)+ &
             wght*zdexp*zfbb*(1.0_fp-zbbfi)
        zomega=0.5_fp*(omegag(j)+omegag(j+1))
        do ig=1,ng
           zdsrat=wght*zdexp*zfth*zbnsvcx(ig)
           sfxrccx(lsbeam)=sfxrccx(lsbeam)+zdsrat
           if(.not.ilouta(js)) then
              sbcxrxs(ig,j,lsbeam)=sbcxrxs(ig,j,lsbeam)+zdsrat
           end if
           sorbhs_ntr(j,ig,lsbeam)=sorbhs_ntr(j,ig,lsbeam)-zdsrat
           sorbn0s(j,ig,lsbeam)=sorbn0s(j,ig,lsbeam)+zdsrat

           !
           srbn02(iz,ig,lsbeam)=srbn02(iz,ig,lsbeam)+zdsrat
           !  simpleminded 1d definition of bcx antitorque for now...
           !  this is consistent with the simpleminded launch in the vol. sce.
           !  in the 1d neutral transport codes
           tqbcxs(j,lsbeam)=tqbcxs(j,lsbeam)+ &
                1.6726e-24_fp*aplasm(ig)*gr2zc(j)*zomega*zdsrat

        end do
        !
        !  end of tracking loop...

     end do
     !
     !  cx loss to walls
     !
     zescap=1.0_fp-zprb(nzcrtk)
     if(ilouta(js)) then
        bpcxots(lsbeam)=bpcxots(lsbeam)+vion*vion*v2toev*wght* &
             zescap
     else
        bpcxins(lsbeam)=bpcxins(lsbeam)+vion*vion*v2toev*wght* &
             zescap
        sbcxxs(lsbeam)=sbcxxs(lsbeam)+wght*zescap
     end if
     bph0was(lsbeam)=bph0was(lsbeam)+zphi0*wght*zescap
     !
     !
     !  decide if ptcl escapes
     !
     zrand=nbrandm_pll(idum)
     if(zrand.gt.zprb(nzcrtk)) then
        sfcxesc(lsbeam)=sfcxesc(lsbeam)+wght
        goto 800
     end if
     !
     !  recapture ptcl - find the zone
     !
     zprbls=0.0_fp
     do i=1,nzcrtk
        if(i.gt.1) zprbls=zprb(i-1)
        itrmax=i
        if(zrand.le.zprb(i)) exit

     end do
     !  source zone indices
     iz=nznbtk(itrmax)
     j=nznttk(itrmax)
     !  length from source to deposition point
     zlndep=zlng(itrmax)-xlngtk(itrmax)*(zprb(itrmax)-zrand)/ &
          (zprb(itrmax)-zprbls)+xltkbf
     !  include intermediate vacuum segment if necessary
     if(itrmax.gt.nztkmd.and.nztkmd.ne.0) then
        zlndep=zlndep+xltkmd
        !%      if(icxo.ne.0) then
        !%        ilintv=4
        !%        zlmid=xltkmd
        !%        zl1=zlng(nztkmd)
        !%      end if
        !%      else
        !%      if(icxo.ne.0) then
        !%        ilintv=2
        !%        zlmid=0.0
        !%        zl1=epsinv
        !%      end if
     end if
     !  plot trajectory to point of recapture
     !%      if(icxo.ne.0) then
     !%      if(nbplot.gt.0) then
     !%        zlplot=zlndep-xltkbf-zlmid
     !%        zleps=0.00001_fp*zlplot
     !%        zlplot=zlplot-ilintv*zleps
     !%        do 440 ipl=1,100
     !%          zlplo=xltkbf+zleps+(ipl-1)*zlplot/99.0_fp
     !%          if(zlplo.gt.(zl1-zleps)) zlplo=zl1+zlmid+zleps
     !%          if(nbplot.eq.1) then
     !%            call nbi_tkxith(1,zlplo,zxi,zth,ier)
     !%            if(ier.ne.0) call bad_exit
     !%            call xmomry(zxi,zth,0,zrplo,zyplo,zdummy)
     !%            zzx(ipl+1)=zrplo
     !%            zzy(ipl+1)=zyplo
     !%          else
     !c%            zcpa=tkcpa(zlplo,ier)
     !c%            if(ier.ne.0) call bad_exit
     !c%            zzx(ipl+1)=zcpa*vion
     !c%            zzy(ipl+1)=vion*sqrt(max(0.0,(1.0-zcpa*zcpa)))
     !%            zzx(ipl+1)=0.0
     !%            zzy(ipl+1)=0.0
     !%          end if
     !% 440    continue
     !%      call line(3)
     !%      call r8_cplot(zzx,zzy)
     !%      call tsend
     !%      call line(0)
     !%      end if
     !%      end if
     !
     !    russian roulette at recapture pt
     !
     if(icxo.eq.0) then
        zrr=1.0_fp             !do not apply rr
     else
        zrr=min(1.0_fp,wght/wbav(j,lsbeam)/fact_wbav)
        if(nbrandm_pll(idum).gt.zrr) then
           sfcxrr(lsbeam)=sfcxrr(lsbeam)+wght
           cycle
        end if
     end if
     !  ptcl survives
     !
     zwghtp=wght
     wght=wght/zrr

     !  recapture source and power
     sfcxrr(lsbeam)=sfcxrr(lsbeam)+(zwghtp-wght)
     sfcxrmc(lsbeam)=sfcxrmc(lsbeam)+wght
     if(ilouta(js)) then
        bpcxrxs(lsbeam)=bpcxrxs(lsbeam)+vion*vion*v2toev*wght
     else
        bpcxris(lsbeam)=bpcxris(lsbeam)+vion*vion*v2toev*wght
     end if
     bphrcas(lsbeam)=bphrcas(lsbeam)+zphi0*wght
     !
     !  get initial orbit parameters
     !
     if(nlbflr) then
        inoflr=0
     else
        inoflr=1
     end if
     !
     call depsub(zlndep,zsinth,zcosth,zsinph,zcosph,inoflr, &
          zphi00,zphigc,ze00,zegc,zbt0,zbtgc)
     !
     call nb0jxb(wght,ze00,zegc,zphi00,zphigc,zbt0,zbtgc, &
          n_flr(js),ngc,xiflr,xion,3,.true.)

     if(icxo.eq.0) cycle
     !
     !  plot guiding ctr
     !%      if(nbplot.gt.0) then
     !%      if(nbplot.eq.1) then
     !%      call xmomry(xion,th,0,zrplo,zyplo,zdummy)
     !%      zzx2(2)=zrplo
     !%      zzy2(2)=zyplo
     !%      else
     !%      zzx2(2)=xksid*vion
     !%      zzy2(2)=sqrt(max(0.0,(1.0-xksid*xksid)))
     !%      end if
     !%      call symbl(2)
     !%      call line(-1)
     !%      call r8_cplot(zzx2,zzy2)
     !%      call tsend
     !%      end if
     !
     !  return ptcl to orbit  ...  or save in splitting array
     if(ireto.eq.0) then
        ireto=1
        zxista=xiflr  ! for radial current tracking
        zxo=xion
        ztho=th
        zwghto=wght
        zkso=xksid
        zvo=vion
     else
        !  save ptcl in splitting array for later reentry into orbit
        nbsplt=nbsplt+1
        if(nbsplt.gt.minbs_expnd) then
           call nbi_spltexpnd(ierr) ! minbs_expnd will be changed in nutrav
           if(ierr.ne.0) then
              jbsplm=minbs_expnd
              call errorlog('nutrav: insufficient array space')
              call errorlog('for ptcl splitting. increase parameter minbs')
              call errorlog('in nbi_dimensions.f90 and recompile.')
              call errorlog('current value = ',jbsplm)
              if(myidd.eq.0)then
                 write(nout,9999) jbsplm
              end if
9999          format(' ??? error in nutrav-- insufficient array space'/ &
                   '     for ptcl splitting.  increase parameter "minbs"'/ &
                   '     in "nbi_dimensions.f90" and recompile.'/ &
                   '     current value = ',i5)
              call bad_exit(" nutrav : not enough space to split")
           end if
        end if
        rbsa(nbsplt)=xion
        xisbsa(nbsplt)=xiflr
        thbsa(nbsplt)=th
        vbsa(nbsplt)=vion
        xksbsa(nbsplt)=xksid
        wtbsa(nbsplt)=wght
        tbsa(nbsplt)=bitime(1)
        nbsbsa(nbsplt)=nbsce
        nbebsa(nbsplt)=nbien
        einjbsa(nbsplt)=einj0
        pinjbsa(nbsplt)=pinj_ay
     end if
     cycle
     !
     !  ptcl escapes
     !
800  continue
     !%      if(nbplot.gt.0) then
     !%      zzx(1)=2.0_fp
     !%      zzy(1)=2.0_fp
     !%      if(nbplot.eq.1) then
     !%      zzx(2)=rsta
     !%      zzy(2)=ysta
     !%      else
     !%      zzx(2)=vpll
     !%      zzy(2)=vperp
     !%      end if
     !  careful - track may never have entered plasma
     !%      if(nzcrtk.eq.0) then
     !%        if(nbplot.eq.1) then
     !%          zzx(3)=rsta-raxis*nsnccw*zsinpa(js)
     !%          zzy(3)=ysta+raxis*nsnccw*zcospa(js)
     !%        else
     !%          zzx(3)=0.0
     !%          zzy(3)=0.0
     !%        end if
     !%      else
     !%        zlpl=0.99999_fp*zlng(nzcrtk)
     !%        if(nbplot.eq.1) then
     !%          call nbi_tkxith(1,zlpl,zzxi,zzth,ier)
     !%          if(ier.ne.0) call bad_exit
     !%          call xmomry(zzxi,zzth,0,zzx(3),zzy(3),zdummy)
     !%        else
     !c%          zcpa=tkcpa(zlpl,ier)  ! no longer available...
     !c%          zzx(3)=zcpa*vion
     !c%          zzy(3)=vion*sqrt(max(0.0,(1.0-zcpa*zcpa)))
     !%          zzx(3)=0.0
     !%          zzy(3)=0.0
     !%        end if
     !%      end if
     !%      call line(3)
     !%      call symbl(0)
     !%      call r8_cplot(zzx,zzy)
     !%      call tsend
     !%      call line(0)
     !%      zzx(1)=100.0_fp
     !%      zzy(1)=100.0_fp
     !%      end if
  end do
1001 continue
  !
  !==============================
  !
  if(icxo.eq.0) return
  if(ireto.eq.0) then
     !  all ptcls escaped or russian roulette suicides
     nbicx=1
     return
  else
     !  return a ptcl to orbit (velocity vion unchanged)
     xistate=zxista
     xion=zxo
     th=ztho
     wght=zwghto

     xksid=zkso
     vion=zvo
  end if
  !
  return
  !-----------------------------------
contains
  subroutine alloc_sa_arrays
    if(allocated(zxsa)) then
       if(size(zxsa).lt.jsplit)then
          deallocate(zxsa,zrsa,zysa)
          allocate(zxsa(jsplit),zrsa(jsplit),zysa(jsplit))
       end if
    else
       allocate(zxsa(jsplit),zrsa(jsplit),zysa(jsplit))
    end if
    zxsa=0.0_fp
    zrsa=0.0_fp
    zysa=0.0_fp

  end subroutine alloc_sa_arrays
end subroutine nutrav
!******************** end file nutrav.for ; group nutrav ******************
