!******************** start file dep.fnb ; group nubeam ****************
!***************************************************************
!
!  note (dmc feb 16 2009): the subroutine name dep is changed to nbidep
!  to satisfy request of scidac facets collaboration...
!
subroutine nbidep(idepo,ierr)
  !
  use nbatom_mod
  use nbi_strline, only : nbi_track_call,nbi_tkxith
  use tube_mod, only: tube_addn0,tube_intersect_nbline
  use dep_track_counter
  use logmod
  use nbdep_debug, only : icount_depall
  !
  !
  !      this subroutine deposits new fast ions
  !
  !
  !      updated to use subroutine nbline, which simulates
  !      a fairly realistic beam-line divergeance geometry
  !
  !      also updated to use trkrb, to give more efficient
  !      tracking across the plasma, and centralized
  !      sight-line geometry calculations. permits all of
  !      beam sight-line to be used for electron source
  !      rate calculation, and beam neutral density
  !
  !
  !  dmc may 2007 -- accumulate deposition data in (mpi dpt group)
  !      track* arrays, if switch is set and 1 <= kdep0 <= ndep0.
  !      control is with local logical :: variable   iltrk0   .
  !
  !  dmc nov 2005 -- parallelization-related changes: (a) removed
  !      calling argument sum-- cf ewghtdep(...) is nbi_com now;
  !      (b) calling argument idepo=1 => new mc particle wanted;
  !      expect idepo=0 otherwise.
  !
  !  pms july 14 1994 -- use of btiltc(4,..) changed for uda
  !
  !  dmc  4 may 1994 -- use excited states fit data if nsigexc.gt.0
  !
  !  dmc 13 may 1993 -- adding call to excited states model interface
  !
  !      rob goldston  6-18-82
  !
  !      upgrade to deal with beam distribution region not equal to
  !      plasma region
  !      rob goldston  july 19, 1983
  !
  !
  !  modif dmc oct 1984-- support more variability from beamline
  !  to beamline
  !
  !  modif mar 1985 -- calculation of net atomic physics (neutral stopping
  !    cross-sections) moved to subroutine nbatom
  !    ** see subroutine nbatom **
  !
  !  modif dmc jun 1987-- compute sums to avg slowing down times,
  !                       p.a. scattering times
  !
  !  mod dmc jan 1988 --- support deposition for fpp code
  !                       if nlfbmfpp(lsbeam) we are depositing for fpp, not
  !                       monte carlo code
  !
  !  mod tbt oct 1991 --- nbatom upgrade. change code to assume nlvphi=t.
  !
  !  dmc jan 1992 -- support deposition of heavy neutral beams (neon,
  !    argon,krypton,xenon) to support analysis of tftr experiments.
  !    simplified atomic physics model is employed:  heavy neutrals are
  !    considered ionized only; the possibility of charge exchange is
  !    ignored; losses of other ion species via cx with heavy neutrals
  !    is not computed.
  !
  !  mod dmc jun 1992 -- support new option to compute 3d neutral density
  !    in a box aligned along a particular beam axis.  this mode is
  !    signalled when ndepbox.lt.0 in common.  in this mode, only the
  !    neutral density is output; all other sums are disabled.
  !
  !-----------------------------------------------------------
  !
  ! were passed, now common:
  !  xjav-- estimated avg zone index <j> for deposited ptcl's 1st orbit
  !
  !  bitime(1)-- deposition time (relative to start of step)
  !
  use nbi_com
  use nbi_mpi_com
  use nbi_random_dep
  !
  use nbdep_debug  !  debug module
  !
  !%      use sglib_mod
  !--------------------------------------
  !============
  implicit none
  integer :: idbgct,ib,ien,idum,ihit,ino,iuncbd,irng,incb,i,j,iz
  integer :: iermap,ig,indi,indp1,isb,ier,inum,ii,iuncp,iuncon
  integer :: itrap,jiav,iftan,idir,izb,ioffn,iofff,indf,indn
  integer :: itrmax,inoflr,itrk_sav1,itrk_sav2,itrk0,mtrk0_prev
  !============
  real(fp) :: zthea,zphia,ze,zltanr,zcth,zltan,zexp
  real(fp) :: znsigl,zlvac,zlava,zlp,zlav,zlvacp,zlavap,zrloc,zyloc
  real(fp) :: zvtor,zvr,zvy,zxi,zth,zvpll,zeovab,zintri,zintrp
  real(fp) :: zexcfac,zbnsvil,zbnsvis,zthnsg,zbbnsv,zbbnsvi,zbbfi
  real(fp) :: zbmsexc1d,zsigvnt1d
  real(fp) :: zbmsexc2d, zsigvnt2d
  real(fp) :: zsigvnt
  real(fp) :: zbbnsg,zmfpi,zfplas,zfbb,zdnsl,zexpol,zdexp,zexpdl
  real(fp) :: zexpav,zxl,zthl,zl2,zll,zvion,xjadep,zeiav,zphiav,zxhw
  real(fp) :: zwgpro,zweinv,zfrac1,ztop,zdn0,zfrci,zdepl,zdebb,znsv
  real(fp) :: zsigvx,zereli,zsigvx6,zsigvx8,znsvii,znii,zsigvii
  real(fp) :: zsigvcx,zdsrat,zomega,zvrel,zcf,ztausl,zef,ztaupa,zwge
  real(fp) :: zrand,zprbp,zprbls,zfaci,zlndep,zxiavl,zxiavx,zxiavi
  real(fp) :: zl0,zlmid,zl1,zl3,zphi0,zphigc,ze0,zegc,zbt0,zbtgc
  real(fp) :: zsinth,zcosth,zsinph,zcosph
  real(fp) :: zxhor,zyhor,zlntot
  logical :: jdbg
  real(fp) :: fdep_k
  real(fp) :: zx(3)          ! location of aperture for tube calculation
  real(fp) :: zvec(3)        ! direction of beam sightline for tube calculation
#ifdef __DEBUG
  real(fp) :: z_phi,xx,yy,rb(3),r(3)
  integer :: id_box(3)
#endif
  !============
  !
  !----------------------------
  ! passed arguments:
  !
  integer, intent(in) :: idepo   ! =1 for new ptcl; =0 otherwise.
  integer, intent(out) :: ierr   ! exit status code, 0=ok.
  !
  !  (if ierr is set, it is a serious, unexpected error).
  !
  !----------------------------
  !
  integer :: icline

  integer :: imissed                           ! count number of misses of plasma
  integer, parameter :: imissed_abort = 1000   ! maximum allowed number of misses
  !
  real(fp) :: zprb(8*mj), zlng(8*mj), zedep(8*mj)
  !
  real(fp) :: zshf0(8*mj),zxiav(8*mj),zvpov(8*mj),zstat(8*mj)
  !
  real(fp) :: zbnsvii,zbnsvie,zbnsviz
  real(fp) :: zbnsvcx(mig),zxfrach,zxfrache
  real(fp) :: zbnsvto  !  tbt 10/17/91
  !
  real(fp) :: zbmod,zbdir(4)            ! dmc -- b vector
  !
  real(fp) :: zbbnsva(2,mibs)
  !
  real(fp) :: zcv0(3),zsign(3)
  !
  !%      real(fp) :: zzx(8*mj),zzy(8*mj),zzx2(2),zzy2(2),zdummy(2,2)
  !%      real(fp) :: zlplo,zrplo,zyplo
  !%      integer :: izork,ipl
  !
  !-------------------------------------
  !  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
  !
  !-------------------------------------
  !
  logical idepbox    ! true if this deposition is for beam in a box
  logical :: ideptube   ! true if this deposition is for beam in a tube
  logical :: iltrk0     ! true to save data on beam neutral track
  logical :: nl_beam_set ! use user's setting for beam source and energy fraction
  logical :: ilbms1d !ture - beam-beam rections are included and for adas beam
  !stopping coefficient and enhancement factor calculation
  !turn on
  !
  !%      logical :: lyes      ! uread function
  !
  integer :: i_zone,i_dir
  !------------------------------------------------------------------
  !  for debugging...
  !  serial mode
  !
  !
  !  these might be reset in a debug version:
  !
  integer, parameter :: icount_debug1 = 0
  integer, parameter :: icount_debug2 = 0
  integer, parameter :: icount_debug3 = 0
  integer, parameter :: icount_debug4 = 0
  !
  logical :: ilwgpro
  integer :: idotrk         ! related to beam in the box neutral ptcl deposition
  ! if =1 then track will be used for a 3d halo
  ! if =0 then add neutral to a box for a density only
  !
  !------------------------------------------------------------------
  !
  !  ztdep replaced by kdebug1...xph_orb in nbi_com
  !
  !  common block also used in depnrm -- testing fpp first orbit loss
  !    accounting
  !   xn_dep,xn_orb-- #/sec deposited, left orbitting
  !   xp_dep,xp_orb-- power deposited, left orbitting
  !   xph_dep,xph_orb-- ang. mo./sec deposited, left orbitting
  !
  !   differences to be accounted by:  prompt orbit losses
  !   jxb torque and associated work term
  !
  data idbgct/0/
  !
  !------------------------------------------------------------------
  !
  data zsign/1._fp,1._fp,1._fp/
  !
  !------------------------------------------------------------------
  !%      zzx=0; zzy=0; zzx2=0; zzy2=0;
  !
  ierr=0
  !
  if(nlbflr) then
     inoflr=0
  else
     inoflr=1
  end if
  !
  !     settings for code debugging
  !
  jdbg = idbg_out(0)

  kdebug1=0                 ! set =1 for debug code
  kdebug2=0                 ! set =1 for debug code
  if(idepo.eq.0) kdebug1=0
  !
  !
  ilbms1d=(nlbbcx.and.(nsigexc_use.eq.1).and. &
       (lev_nbidep.eq.2).and.(nbbcx_bb.eq.0))

  idepbox  = mdepbox<0 !it supports mpi
  ideptube = ndeptube<0
  !
  ! for beam in box option lsbeam is unknown yet
  !
  if(idepbox) then
     iltrk0 = .false. !track data will not be gathered at this time
  else

     iltrk0 = nltrk_dep0.and.(mibs2use(lsbeam).gt.0)
     iltrk0 = iltrk0 .and.(.not.ideptube)
     iltrk0 = iltrk0 .and. (kdep0.ge.0) .and. &
          (kdep0.lt.mdep0)
     !
     !
     !%      if((.not.(idepbox.or.ideptube)).and.(nbplot.gt.0)) then
     !%        zzx(1)=101.0_fp
     !%        zzy(1)=101.0_fp
     !%        zzx2(1)=1.0_fp
     !%        zzy2(1)=1.0_fp
     !%        call line(0)
     !%        call symbl(0)
     !%      end if
     !
     !---------------------------------------------------------------------
     !
     !     debug code
     !
     !
     if(kdebug1.eq.1) then
        if(ewghtdep(lsbeam).eq.0.0_fp) then
           xn_dep(lsbeam)=0.0_fp
           xn_orb(lsbeam)=0.0_fp
           xp_dep(lsbeam)=0.0_fp
           xp_orb(lsbeam)=0.0_fp
           xph_dep(lsbeam)=0.0_fp
           xph_orb(lsbeam)=0.0_fp
        end if
        xwght0(lsbeam)=1.0_fp
     end if
     !
     !---------------------------------------------------------------------
     !     error check
     !
     if(nlfprod(lsbeam)) then
        call errorlog('dep called for fusion product deposition')
        call errorlog('    should call subroutine depfus')
        if(myidd.eq.0) then
           write(nonlin,*) &
                ' ?? dep called for fusion product deposition'
           write(nonlin,*) '    should call subroutine depfus'
        end if
        ierr=1
        return
     end if
  end if
  !
  !  deposition section
  !
  !     starting parameters-- choose a chord from ion source /
  !      beamline geometry
  !
  imissed = 0
10 continue
  !
  if(idepbox) then
     ! beam in a box 3d neutral density
     call randbl_box(ib,ien,ierr) !lsbeam and ldep will be taken from this routine
     if(ierr.ne.0) return
     call inispc
     bitime(1)=0.0_fp
  else if (ideptube) then
     ! beam in a tube 3d neutral density
     bitime(1)=0.0_fp
     ib=nbstube(ldep)
     ien=nbetube(ldep)
  else
     !
     !  normal deposition call
     !
     !  choose random start time within timestep
     bitime(1)=nbrandm_dep(idum)*delt !deposition time relative to tbm1
     bitime(1)=max(bitime(1),delt*1.0e-12_fp) !bitime=0 is considered as old mc particle
     !
     !  choose random beamline, energy fraction
     call randbl(ib,ien,ierr)
     if(ierr.ne.0) return
     zidep0(lsbeam)=zidep0(lsbeam)+1 !counter of#deposited ptcls on given cpu

     if ((ndep_beam_use.gt.0).and.iltrk0) then
        ! if ndep_beam_use set =0 --
        ! assuming all beam sources are counted
        call set_dep_ib_ien(ib,ien,nl_beam_set)

        iltrk0=iltrk0.and.nl_beam_set
     end if

     if(iltrk0) then
        kdep0=kdep0+1            ! this specie
        kdep0sum = kdep0sum + 1  ! all species
        zitrk0(lsbeam)=zitrk0(lsbeam)+1 !counter of#deposited track on given cpu
     end if
     !
  end if
  !
  !     take into account ptcl weight profile
  !
  ilwgpro=.false.
  zwgpro=1.0_fp
  zweinv=zwgpro
  if(abs(xwght0(lsbeam)-1.0_fp).ge.1.0e-10_fp) then
     ilwgpro=.true.         ! ptcl weight profile in effect
  end if
  !
  nbsce=ib
  nbien=ien
  !
  !  set up beam track starting coordinates
  !
  if(nbdep_debug_active.and.nbdep_centerline) then
     icline=1                       ! centerline track
  else
     icline=0                       ! random track per divergence, etc.
  end if
  call nbline(ib,icline)
  if(rsta.lt.0.0_fp) then
     call errorlog('nubeam dep:  error in nbline.')
     if(myidd.eq.0) then
        write(nonlin,*) ' ? nubeam dep:  error in nbline.'
     end if
     ierr=1
     return
  end if
  !
  if(jdbg) then
     call infolog( ' ** icount_debug match: lsbeam,ib,ien: ')
     call infolog('lsbeam = ',lsbeam)
     call infolog('ib = ',ib)
     call infolog('ien = ',ien)
     call infolog('rsta = ',  rsta)
     call infolog('ysta = ',  ysta)
     call infolog('thea = ',  thea)
     call infolog('phia = ',  phia)
     if(myidd.eq.0) then
        write(nonlin,*)' ** icount_debug match: lsbeam,ib,ien: ', &
             lsbeam,ib,ien

        write(nonlin,'(a,4(1x,1pe17.10))') &
             '     rsta,ysta,thea,phia: ',rsta,ysta,thea,phia
     end if
  end if
  !
  !%      if((.not.(idepbox.or.ideptube)).and.(nbplot.gt.0)) then
  !%        zzx(2)=rsta
  !%        zzy(2)=ysta
  !%      end if
  !
  if (ideptube) then
     ! check if trajectory intersected the tube and store
     ! intersection results in tube_mod
     call tube_intersect_nbline(ldep,ihit)
     if (ihit==0) return ! missed tube so exit
  end if

  zthea=degrad*thea
  zphia=degrad*phia
  !
  zsinth=sin(zthea)
  zcosth=cos(zthea)
  !
  zsinph=sin(zphia)
  zcosph=cos(zphia)
  !
  ze=eiona(ien,ib)
  !
  vion=viona(ien,ib)
  einj0=ze
  !
  !     deposition
  !
  !  follow track through plasma, compute length through each zone
  !
  !  straln is gone...      call trkrl
  !
  call nbi_track_call(rsta, ysta, phia, thea, ierr)
  if(ierr.ne.0) return
  !
  if(nbdep_debug_active) call nbdep_debug_init

  idotrk=0
  if(levmod_halo.eq.2) then
     idotrk=1
  end if

  if(idepbox) then
     call addn0(0.0_fp,1.0_fp,0.0_fp,vion,idotrk,lsbeam) !initializations only
  end if
  !
  !  check if track misses plasma entirely
  !
  if(nzcrtk.eq.0) then
     ! bug fixed 11/18/2013 by mg
     if(iltrk0) then
        kdep0 = kdep0 - 1
     end if
     if(idepbox) then
        call addn0(10.0_fp*xltkwl,1.0_fp,0.0_fp,vion,idotrk,lsbeam)
        return
     else if (ideptube) then
        call tube_addn0(0._fp,10.0_fp*xltkwl/vion,1.0_fp,0.0_fp)
        return
     else
        ewghtdep(lsbeam)=ewghtdep(lsbeam)+ze
        bpshins(lsbeam)=bpshins(lsbeam)+ze
        bpshna(ib,ien)=bpshna(ib,ien)+ze
        sfshine(lsbeam)=sfshine(lsbeam)+1.0_fp
        call nb0loc(10.0_fp*xltkwl,zsinth,zcosth,zsinph,zcosph, &
             zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
        bphshins(lsbeam)=bphshins(lsbeam)+zrloc*zvtor*rmassb
     end if
     !  get new track
     if(nbdep_debug_active.and.nbdep_centerline) then
        return                      ! centerline missed plasma
     else
        imissed = imissed+1
        if (imissed>=imissed_abort) then
           print '(a,i5,a)', '?dep: track missed the plasma ' &
                ,imissed,' consecutive times'
           ierr=1 ; return
        end if
        goto 10
     end if
  else
     imissed=0   ! in case of another goto 10
  end if
  !

90 continue
  !
  if(iltrk0) then
     if(nztkmd.gt.0) then
        xntrack(kdep0,mibs2use(lsbeam)) = nzcrtk + 2
     else
        xntrack(kdep0,mibs2use(lsbeam)) = nzcrtk
     end if
     mtrk0_prev = mtrk0
     itrk0 = mtrk0_prev
     mtrk0 = mtrk0+xntrack(kdep0,mibs2use(lsbeam))+0.1_fp

     call nbi_expand_track_zone

     call nb0loc(xltkbf,zsinth,zcosth,zsinph,zcosph, &
          zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
     track_rin(kdep0,mibs2use(lsbeam))=zrloc
     track_zin(kdep0,mibs2use(lsbeam))=zyloc + plhgt
     track_phiin(kdep0,mibs2use(lsbeam))=esta + &
          nsnccw*nsjdotb*atan2(zyhor,zxhor)/degrad

     track_sinj(kdep0,mibs2use(lsbeam)) = 1.0_fp
     track_einj(kdep0,mibs2use(lsbeam)) = einj0 !eiona(ien,ib)
     track_nbsce(kdep0,mibs2use(lsbeam))=ib
     track_nbien(kdep0,mibs2use(lsbeam))=ien
     track_time(kdep0,mibs2use(lsbeam)) = tbm1 + bitime(1)
  end if

  !
  !  trajectory length to tangency radius
  !
  !  zltan= distance to tangency pt.
  !   ino=1 if trajectory starts out headed in (towards smaller rmajor)
  !      =2 if trajectory starts headed out
  !
  !  midplane projected distance
  zltanr=max(0.0_fp,rsta*zcosph)
  if(zltanr.gt.0.0_fp) then
     zcth=zcosth
     if(zcth.le.epslon) then
        zltan=0.0_fp
        ino=2
     else
        !  distance along ptcl trajectory
        zltan=zltanr/zcth
        ino=1
     end if
  else
     zltan=0.0_fp
     ino=2
  end if
  !
  !  initialize orbit shift computation-- for each new track
  if(.not.nlheavy(lsbeam)) call xjaini
  !
  !%      izork=ib
  iuncbd=0                  ! lost orbit region bdy's counted
  !
  !
  !  track through plasma, compute deposition profiles by plasma zone,
  !  compute probability of deposition in each zone.  this probability
  !  includes a monte carlo skewing controlled by namelist input
  !  wghta-- which (for wghta.gt.1) allows the deposition of more ions
  !  of lower weight near the center, and fewer ions of greater weight
  !  towards the edge of the plasma
  !
  irng = 1      ! 1st guess to range # of energy grid in nbeloc.
  zexp=1.0_fp
  znsigl=0.0_fp
  itrk_sav1=0
  itrk_sav2=0
  !  vacuum segment length
  zlvac=xltkbf
  if(idepbox) then
     call depsub(zlvac,zsinth,zcosth,zsinph, &
          zcosph,inoflr,zphi0,zphigc,ze0,zegc,zbt0,zbtgc) !have to call to get xksid
     call addn0(zlvac,1.0_fp,0.0_fp,vion,idotrk,lsbeam)
#ifdef __DEBUG
     call nb0loc(xltkbf,zsinth,zcosth,zsinph,zcosph, &
          zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
     z_phi=esta*degrad +nsnccw*nsjdotb*atan2(zyhor,zxhor) !rad
     call eq_xcyl(1,zrloc*0.01_fp,z_phi,xx,yy)
     r(1)=xx*100._fp
     r(2)=yy*100._fp
     r(3)=zyloc
     call car2box(ldep,r,rb,id_box)
     call box2car(ldep,rb,r)
#endif
  else if (ideptube) then
     call tube_addn0(0._fp, zlvac/vion, 1.0_fp, 0.0_fp)
  end if

  if(nbdep_debug_active) then
     call nbdep_coord_record(1,zlvac,zsinth,zcosth,zsinph,zcosph)
     fbeam_nbdep(1)=1
     incb=1
  end if
  !
  !  loop thru track segments
  !
  if(nzcrtk.gt.8*mj) then
     call errorlog('dep: nzcrtk > limits')
     call errorlog('nzcrtk = ',nzcrtk)
     call errorlog('limit = ',8*mj)
     call errorlog &
          ('probably due to code error in xstraln or its srs')
     if(myidd.eq.0) then
        write(nonlin,*) ' ?? dep: nzcrtk = ',nzcrtk,' > limit ',8*mj
        write(nonlin,*) &
             '    (probably due to code error in xstraln or its srs'
     end if
     ierr=1
     return
  end if
  !
  zlava=xltkbf
  zlng=0.0_fp

  do i=1, nzcrtk
     !
     !  transp zone index
     j=nznttk(i)
     !  beam zone index
     iz=nznbtk(i)

     !
     if(j.gt.ledge.or.iz.eq.0) idotrk=0
     !  distance tracked
     !
     if(iltrk0) then
        if(itrk0.eq.size(track_xl)) then
           !               call nbi_expand_track_zone
        end if
     end if
     if(i.eq.1) then
        zlp=0.0_fp
        zlng(1)=xlngtk(1)
        if(iltrk0) then
           itrk0=mtrk0_prev+1
           track_xl(itrk0)=xlngtk(1)
        end if
     else
        zlp=zlng(i-1)
        zlng(i)=zlp+xlngtk(i)
        if(iltrk0) then
           itrk0=itrk0+1
           track_xl(itrk0)= &
                track_xl(itrk0-1)+xlngtk(i)
        end if
     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
        zlvacp=zlvac
        zlvac=zlvac+xltkmd
        zlava=zlng(i-1)+zlvac
        ino=2
        if(idepbox) then
           if(zlvac.gt.zlvacp) then
              call depsub(zlvac+zlng(i-1),zsinth,zcosth,zsinph, &
                   zcosph,inoflr,zphi0,zphigc,ze0,zegc,zbt0,zbtgc) !have to call to get xksid
              call addn0(zlvac+zlng(i-1),zexp,0.0_fp,vion,idotrk, &
                   lsbeam)
           end if
        else if (ideptube) then
           call tube_addn0((zlvacp+zlng(i-1))/vion, xltkmd/vion, &
                zexp, 0.0_fp)
        end if
        if(nbdep_debug_active) then
           call nbdep_coord_record(i+incb, &
                zlp+zlvac,zsinth,zcosth,zsinph,zcosph)
           fbeam_nbdep(i+incb)=fbeam_nbdep(i+incb-1)
           incb=incb+1
        end if
     end if
     zlavap=zlava       ! remember old value
     zlava=zlav+zlvac
     !
     if(nbdep_debug_active) then
        call nbdep_coord_record(i+incb, &
             zlng(i)+zlvac,zsinth,zcosth,zsinph,zcosph)
        lz_nbdep(i)=(lb_nbdep(i+incb)+lb_nbdep(i+incb-1))/2
     end if
     !
     !  fetch position and velocity in (r,y,phi) coordinate system
     call nb0loc(zlava,zsinth,zcosth,zsinph,zcosph, &
          zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)

     !
     !  dmc dec 4 1997 -- spend some cpu time -- get actual b vector and
     !  accurate vpll
     !
     call nblocb(zrloc,zyloc,zxi,zth,zbmod,zbdir,1,nfast_map,iermap)
     if(iermap.gt.1) then
        call errorlog &
             ('nbidep nblocb: unexpected (r,z)->(xi,theta) map error!')
        if(myidd.eq.0) then
           write(nonlin,*) &
                ' ?? nbidep nblocb: unexpected (r,z)->(xi,theta) map error!'
        end if
        ierr=1
        return
     end if

     zvpll=zvtor*zbdir(1)-nsnccw*(zvr*zbdir(2)+zvy*zbdir(3))
     !
     !  inverse local mean free path
     !  look up plasma interaction rates
     !
     call nberel(zeovab,j,zrloc,zvtor,zvr,zvy)
     !
     zedep(i)=zeovab*abeam
     !
     if(nbdep_debug_active) then
        ebplas_nbdep(i)=zeovab*abeam
        ti_nbdep(i)=ti(j,2)
        te_nbdep(i)=te(j,2)
        ne_nbdep(i)=rhoel(j,2)
        nimp_nbdep(i)=rhi(j,2)
        aimp_nbdep(i)=aimpj(j)
        zimp_nbdep(i)=xzimpj(j)
        !
        nhs_nbdep(i)=0
        ahs_nbdep(i)=0
        nhes_nbdep(i)=0
        ahes_nbdep(i)=0
        nlis_nbdep(i)=0
        alis_nbdep(i)=0
        do ig=1,ng
           if(abs(backz(ig)-1.0_fp).le.0.1_fp) then
              nhs_nbdep(i)=nhs_nbdep(i)+rhob(j,ig,2)
              ahs_nbdep(i)=ahs_nbdep(i)+rhob(j,ig,2)*aplasm(ig)
           else if(abs(backz(ig)-2.0_fp).le.0.1_fp) then
              nhes_nbdep(i)=nhes_nbdep(i)+rhob(j,ig,2)
              ahes_nbdep(i)=ahes_nbdep(i)+rhob(j,ig,2)*aplasm(ig)
           else if(abs(backz(ig)-3.0_fp).le.0.1_fp) then
              nlis_nbdep(i)=nlis_nbdep(i)+rhob(j,ig,2)
              alis_nbdep(i)=alis_nbdep(i)+rhob(j,ig,2)*aplasm(ig)
           end if
        end do
        if(nhs_nbdep(i).gt.0) ahs_nbdep(i)= &
             ahs_nbdep(i)/nhs_nbdep(i)
        if(nhes_nbdep(i).gt.0) ahes_nbdep(i)= &
             ahes_nbdep(i)/nhes_nbdep(i)
        if(nlis_nbdep(i).gt.0) alis_nbdep(i)= &
             alis_nbdep(i)/nlis_nbdep(i)
        !
        nfas_nbdep(i)=0
        do ig=1,nmini
           nfas_nbdep(i)=nfas_nbdep(i)+rhmin(j,2,ig)
        end do
        do ig=1,nsbeam
           nfas_nbdep(i)=nfas_nbdep(i)+rhbs(j,2,ig)
        end do
     end if
     !
     if(nlheavy(lsbeam)) then
        !  heavy neutral beam
        zexcfac=1.0_fp
        do ig=1,ng
           zbnsvcx(ig)=0.0_fp
        end do
        zbnsvii=1.0_fp
        zbnsvie=0.0_fp
        zbnsviz=0.0_fp
        zbnsvto=hvsigi(j,1,lsbeam)
     else
        !  normal heating beam (isotope of hydrogen or helium)
        call nbeloc(zeovab,irng,indi,zintri,indp1,zintrp)
        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

     end if
     !
     !--------------------------------------
     !
     !  compute beam-beam n*sigma*v total and ionization fraction
     !
     if(nlbbcx.and.(.not.nlheavy(lsbeam))) then
        do isb=1,nsbeam
           zbbnsva(1,isb)=sn00i2(ib,ien,iz,isb)
           zbbnsva(2,isb)=sn00x2(ib,ien,iz,isb)
        end do
     else
        do isb=1,nsbeam
           zbbnsva(1,isb)=0.0_fp
           zbbnsva(2,isb)=0.0_fp
        end do
     end if
     !
     !--------------------------------------
     !
     !  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(nlheavy(lsbeam)) then
        !  no excitation effect for neon and other hi z beams
        zbbxfac=1.0_fp
        zwkfac=1.0_fp
        if(nbdep_debug_active) excit_nbdep(i)=1.0_fp

     else if(nsigexc_use.gt.0) then
        !
        call exc_sig(zexcfac,zbbxfac,zwkfac, &
             zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto, &
             zbbnsva,ierr)
        if(ierr.ne.0) return
        if(nbdep_debug_active)then
           zexcfac = zwkfac(1) !zexcfac for ionization
           excit_nbdep(i)=zexcfac
        end if
        !

     else
        zbbxfac=1.0_fp
        zwkfac=1.0_fp
        if(nbdep_debug_active) excit_nbdep(i)=1.0_fp
     end if
     !
     !--------------------------------------
     !
     !  total ionization fraction (after adjustment)
     !  and thermal plasma stopping n*sigma
     !
     zbnsvis=zbnsvii+zbnsvie+zbnsviz+zbnsvil
     zthnsg= zbnsvto / vion
     !         zthnsg= zbnsvto / viona(ien,ib)
     !
     !  dmc apr 1990
     !  namelist option to modify plasma opacity to new fast neutrals:
     !
     if(ndepmod.eq.1) then
        zthnsg=xdepmod*zthnsg
     end if
     !
     !  net beam-beam stopping n*sigma
     !
     if(nlbbcx.and.(.not.nlheavy(lsbeam))) then
        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
     else
        zbbnsv=0.0_fp
        zbbnsvi=0.0_fp
     end if
     if(zbbnsv.gt.0.0_fp) then
        zbbfi=zbbnsvi/zbbnsv
        zbbnsg=zbbnsv/vion
        !            zbbnsg=zbbnsv/viona(ien,ib)
     else
        zbbfi=0.0_fp
        zbbnsg=0.0_fp
     end if
     !
     !  inverse mean free path:
     !
     zmfpi=zthnsg+zbbnsg
     if(nbdep_debug_active) lmfp_nbdep(i)=1/zmfpi ! mean free path
     !
     !  thermal fraction:
     zfplas=zthnsg/zmfpi
     !  beam-beam fraction:
     zfbb=zbbnsg/zmfpi
     !
     zdnsl=xlngtk(i)*zmfpi
     zexpol=zexp
     !  n*sigma*dl integral
     znsigl=znsigl+zdnsl
     zexp=exp(-znsigl)
     zdexp=zexpol-zexp
     zexpdl=zdexp/zmfpi
     !
     if(nbdep_debug_active) then
        fbeam_nbdep(i+incb)=zexp
        delf_nbdep(i)=zdexp/xlngtk(i)
     end if

     !
     !  ok have survival probability -- increment neutral density box calc.
     !
     if(idepbox) then
        zexpav=0.5_fp*(zexp+zexpol)
        ! 1/(zmfpi*vion) - mean lifetime
#ifdef __DEBUG
        call nb0loc(zlvac+zlng(i),zsinth,zcosth,zsinph,zcosph, &
             zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
        z_phi=esta*degrad +nsnccw*nsjdotb*atan2(zyhor,zxhor)
        call eq_xcyl(1,zrloc*0.01_fp,z_phi,xx,yy)
        r(1)=xx*100._fp
        r(2)=yy*100._fp
        r(3)=zyloc
        call car2box(ldep,r,rb,id_box)
        call box2car(ldep,rb,r)
#endif
        call depsub(zlvac+zlng(i),zsinth,zcosth,zsinph, &
             zcosph,inoflr,zphi0,zphigc,ze0,zegc,zbt0,zbtgc) !have to call to get xksid
        call addn0(zlvac+zlng(i),zexpol, zmfpi,vion,idotrk,lsbeam)
        cycle
     else if (ideptube) then
        call tube_addn0((zlvac+zlng(i)-xlngtk(i))/vion, &
             xlngtk(i)/vion,zexpol, zmfpi*vion)
        cycle
     end if
     !
     !  establish wghta effect on monte-carlo depo. probability
     !
     !  need estimate of current xi,theta on neutral track
     !
     !  (replacing old straln tkxi and tkth calls)
     !
     call nbi_tkxith(1,zlava,zxl,zthl,ier)
     if(ier.ne.0) then
        call errorlog(' dep -- track appears outside plasma.')
        call errorlog('  error detected at l=',zlava)
        call errorlog('  l   rho   theta  ier')
        if (myidd.eq.0)  then
           write(nonlin,*) ' ?dep -- track appears outside plasma.'
           write(nonlin,*) '  error detected at l=',zlava
           write(nonlin,*) '  l   rho   theta  ier'
        end if
        zl2=zlvac+zlng(i)
        inum=50
        do ii=1,inum
           zll=((ii-1)*zl2 + (inum-ii)*zlavap)/(inum-1)
           call nbi_tkxith(0,zll,zxl,zthl,ier)
           call  errorlog(' ' ,zll)
           call  errorlog(' ' ,zxl)
           call  errorlog(' ' ,zthl)
           call  errorlog(' ' ,ier)
           if(myidd.eq.0) then
              write(nonlin,*) zll,zxl,zthl,ier
           end if
           if(ier.eq.0) exit
        end do
        ierr=ier
        if(ierr.ne.0) then
           return
        else
           zlng(i)=(zll-zlvac)*2._fp-zlp
           zlava=zll
        end if
     end if
     !
     !  orbit shift estimated avg zone index (cf xjadep.fnb)
     !
     if(i.gt.1) iuncp=iuncon
     !
     zvion=vion
     if(nlheavy(lsbeam)) then
        zxiav(i)=zxl
        jiav=j
        iuncon=0
     else
        !
        !  ** xjadep ** call -- estimate orbit shift
        !
        zxiav(i)=xjadep(zxl,zthl,zrloc,zyloc,zvion,zvpll, &
             itrap,iuncon,zvpov(i),zeiav,zphiav,jiav,zxhw)
        !
        !
     end if
     if(iuncon.eq.1) then
        jiav=lep1 !unconfined
     end if
     !
     zshf0(i)=zxiav(i)-zxl
     !
     if((i.gt.1).and.(iuncon.ne.iuncp)) iuncbd=iuncbd+1
     !
     zstat(i)=iuncon
     !
     !  sum the bounce averaged plasma frame deposition distribution:
     !
     if(nlfbmfpp(lsbeam)) then
        call fdep(ib,zdexp,zvtor,zrloc,j,zxl, &
             jiav,zxiav(i),zvpov(i),zeiav,zphiav, &
             iuncon)
     end if
     !
     if(ilwgpro) then
        zwgpro=fdep_k(zxiav(i),xwght0(lsbeam))
        !mg            zwgpro=1.0_fp-(1.0_fp-xwght0(lsbeam))*
        !mg     >           (1.0_fp-zxiav(i)**1.5_fp)
        ! mg 01/17/18
        zwgpro=1.0_fp-(1.0_fp-xwght0(lsbeam))* &
             (1.0_fp-zxiav(i)**1.5_fp)
        zweinv=1.0_fp/zwgpro
     end if
     !
     !  zprb(i) will give the unnormalized probability of deposition at
     !    or before trkrb zone i (transp zone nznttk(i))-- including
     !    monte-carlo wghta effect
     !  the last zone's zprb(nzcrtk) is then used to normalize
     !
     if(i.eq.1) then
        zprb(1)=zweinv*zdexp
     else
        zprb(i)=zprb(i-1)+zweinv*zdexp
     end if
     !
     if(iltrk0.and.(.not.nlheavy(lsbeam))) then
        if(i.eq.1) then
           track_pdep(itrk0)=zdexp
           zlndep = xltkbf + track_xl(itrk0) &
                *0.5_fp
        else
           track_pdep(itrk0)= &
                track_pdep(itrk0-1)+zdexp
           zlndep = xltkbf + (track_xl(itrk0-1) &
                + track_xl(itrk0))*0.5_fp
        end if
        !
        !  depsub call (restore vion when done!): get flr step
        !
        call depsub(zlndep,zsinth,zcosth,zsinph,zcosph,inoflr, &
             zphi0,zphigc,ze0,zegc,zbt0,zbtgc)
        track_dr_flr(itrk0)=rmjion-rmjdep
        track_dz_flr(itrk0)=yion-ydep
        track_de_flr(itrk0)=zegc-ze0
        track_vpv_dep(itrk0)=xksid_dep
        track_dvpv_flr(itrk0)=xksid - xksid_dep
        track_xjadep(itrk0)=2*iuncon + itrap

        vion = vion_dep     ! restore vion

        if(i.eq.nztkmd) then
           itrk0=itrk0+1
           itrk_sav1 = itrk0
           track_xl(itrk0)= &
                track_xl(itrk0-1)+xltkmd*0.5_fp
           itrk0=itrk0+1
           itrk_sav2 = itrk0
           track_xl(itrk0)= &
                track_xl(itrk0-1)+xltkmd*0.5_fp
           track_pdep(itrk0-1:itrk0) = &
                track_pdep(itrk0-2)
        end if
     end if
     !
     !  check for passing tangency pt
     if((ino.eq.1).and.(zlng(i).gt.(zltan-zlvac))) then
        ino=2
        iftan=1
        zfrac1=(zltan-zlvac-zlp)/(zlng(i)-zlp)
     else
        iftan=0
     end if
     !
     !  update neut. density and velocity components
     !
     !  bn00 to be replaced with bn002 - maintained temporarily for debug
     !   but the ino index is hardwired to 1
     !
     if(.not.nlsym2b) then
        ztop=1.0_fp
     else
        if(zyloc.gt.0.0_fp) then
           ztop=1.0_fp
        else
           ztop=-1.0_fp
        end if
     end if
     !  toroidal velocity cosine
     zcv0(1)=zvtor/vion
     !  velocity cosine normal to flux surface
     zcv0(2)=(zvr*btiltc(5,iz)-zvy*ztop*btiltc(4,iz))/vion
     !  velocity cosine tangent to flux surface in poloidal direction
     zcv0(3)=(zvr*ztop*btiltc(4,iz)+zvy*btiltc(5,iz))/vion
     !
     zdn0=zexpdl/vion
     bn00(ib,ien,j)=bn00(ib,ien,j)+zdn0
     vpvn00(ib,ien,j)=vpvn00(ib,ien,j)+ &
          zcv0(1)*zsign(1)*zdn0
     !
     !  excited states stuff ...
     !
     if(nsigexc_use.gt.0) then
        excit0(j,ien,lsbeam)=excit0(j,ien,lsbeam)+zdn0*zwkfac(1)
     end if
     !
     if(iftan.eq.0) then
        !  ino=1:  ingoing neutral ino=2:  outgoing neutral
        bn002(ib,ien,ino,iz)=bn002(ib,ien,ino,iz)+zdn0
        do idir=1,3
           vcxbn0(ib,idir,ino,iz)=vcxbn0(ib,idir,ino,iz)+ &
                zcv0(idir)*zsign(idir)*zexpdl

        end do
        if(nsigexc_use.gt.0) then
           do isb=1,nsbeam
              hexbbcx0(ib,ien,ino,iz,isb)= &
                   hexbbcx0(ib,ien,ino,iz,isb)+zdn0*zbbxfac(isb)
           end do
        end if
     else
        !  zone includes tangency pt
        !   part before tangency pt ...
        bn002(ib,ien,1,iz)=bn002(ib,ien,1,iz)+zfrac1*zdn0
        do idir=1,3
           vcxbn0(ib,idir,1,iz)=vcxbn0(ib,idir,1,iz)+ &
                zcv0(idir)*zsign(idir)*zexpdl*zfrac1

        end do
        if(nsigexc_use.gt.0) then
           do isb=1,nsbeam
              hexbbcx0(ib,ien,1,iz,isb)=hexbbcx0(ib,ien,1,iz,isb) &
                   +zfrac1*zdn0*zbbxfac(isb)
           end do
        end if
        !
        !   part after tangency pt ...
        zfrci=1.0_fp-zfrac1
        bn002(ib,ien,2,iz)=bn002(ib,ien,2,iz)+zfrci*zdn0
        do idir=1,3
           vcxbn0(ib,idir,2,iz)=vcxbn0(ib,idir,2,iz)+ &
                zcv0(idir)*zsign(idir)*zexpdl*zfrci

        end do
        if(nsigexc_use.gt.0) then
           do isb=1,nsbeam
              hexbbcx0(ib,ien,2,iz,isb)=hexbbcx0(ib,ien,2,iz,isb) &
                   +zfrci*zdn0*zbbxfac(isb)
           end do
        end if
     end if
     !
     !  trapping fraction data
     !
     if((ien.eq.1).and.(.not.nlheavy(lsbeam))) then
        call addtrap0(zdexp,itrap,zeiav)
     end if
     !
     !  deposition profiles vs. transp zones
     sorbes_dep(j,lsbeam)=sorbes_dep(j,lsbeam)+xzbeami*zdexp* &
          (zfplas*zbnsvis +zfbb*zbbfi)
     !
     zdepl=zdexp*zfplas
     zdebb=zdexp*zfbb
     !
     sdep_oa(jiav,lsbeam)=sdep_oa(jiav,lsbeam)+zdexp
     sfdepsc(lsbeam)=sfdepsc(lsbeam)+zdexp
     sfdepiz(lsbeam)=sfdepiz(lsbeam)+zdepl*zbnsvis
     sfdbbcx(lsbeam)=sfdbbcx(lsbeam)+zdebb*(1.0_fp-zbbfi)
     sfdbbiz(lsbeam)=sfdbbiz(lsbeam)+zdebb*zbbfi
     !
     sbedeps(j,lsbeam)=sbedeps(j,lsbeam)+zdepl*zbnsvii
     sbiedeps(j,lsbeam)=sbiedeps(j,lsbeam)+zdepl*zbnsvie
     sbizdeps(j,lsbeam)=sbizdeps(j,lsbeam)+zdepl*zbnsviz
     sbildeps(j,lsbeam)=sbildeps(j,lsbeam)+zdepl*zbnsvil
     sbbxdps(j,lsbeam)=sbbxdps(j,lsbeam)+zdebb*(1.0_fp-zbbfi)
     sbbidps(j,lsbeam)=sbbidps(j,lsbeam)+zdebb*zbbfi
     !
     sbedepe(j,ien,lsbeam)=sbedepe(j,ien,lsbeam)+zdepl*zbnsvii
     sbiedepe(j,ien,lsbeam)=sbiedepe(j,ien,lsbeam)+zdepl*zbnsvie
     sbizdepe(j,ien,lsbeam)=sbizdepe(j,ien,lsbeam)+zdepl*zbnsviz
     sbildepe(j,ien,lsbeam)=sbildepe(j,ien,lsbeam)+zdepl*zbnsvil
     sbbxdpe(j,ien,lsbeam)= &
          sbbxdpe(j,ien,lsbeam)+zdebb*(1.0_fp-zbbfi)
     sbbidpe(j,ien,lsbeam)=sbbidpe(j,ien,lsbeam)+zdebb*zbbfi
     !
     bdepb(j,ien,nbsce)=bdepb(j,ien,nbsce)+zdexp
     !
     if(nbdep_debug_active) then
        delfie_nbdep(i)=zdepl*zbnsvie/xlngtk(i)
        delfii_nbdep(i)=zdepl*zbnsvii/xlngtk(i)
        delfix_nbdep(i)=zdepl*zbnsviz/xlngtk(i)
        delflo_nbdep(i)=zdepl*zbnsvil/xlngtk(i)
        delfcxh_nbdep(i)=0
        delfcxhe_nbdep(i)=0
        zxfrach=0
        zxfrache=0
        do ig=1,ng
           if(abs(backz(ig)-1.0_fp).le.0.1_fp) then
              delfcxh_nbdep(i)=delfcxh_nbdep(i)+zdepl*zbnsvcx(ig)
              zxfrach=zxfrach+zbnsvcx(ig)
           else if(abs(backz(ig)-2.0_fp).le.0.1_fp) then
              delfcxhe_nbdep(i)=delfcxhe_nbdep(i)+zdepl*zbnsvcx(ig)
              zxfrache=zxfrache+zbnsvcx(ig)
           end if
        end do
        delfcxh_nbdep(i)=delfcxh_nbdep(i)/xlngtk(i)
        delfcxhe_nbdep(i)=delfcxhe_nbdep(i)/xlngtk(i)
        delfbbii_nbdep(i)=zdebb*zbbfi/xlngtk(i)
        delfbbcx_nbdep(i)=zdebb*(1.0_fp-zbbfi)/xlngtk(i)
        znsv=zthnsg*vion
        !            znsv=zthnsg*viona(ien,ib)
        svie_nbdep(i)=znsv*zbnsvie/ne_nbdep(i)
        svii_nbdep(i)=znsv*zbnsvii/(nhs_nbdep(i)+ &
             nhes_nbdep(i)+nlis_nbdep(i))
        svix_nbdep(i)=znsv*zbnsviz/max(1.0_fp,nimp_nbdep(i))
        svcxh_nbdep(i)=znsv*zxfrach/max(1.0_fp,nhs_nbdep(i))
        svcxhe_nbdep(i)=znsv*zxfrache/max(1.0_fp,nhes_nbdep(i))
        izb=xzbeami+0.5_fp
        !
        call sigvte_ioniz(izb,1.0e-3_fp*te_nbdep(i),1,zsigvx)
        svie2_nbdep(i)=1.0e6_fp*zsigvx*zexcfac*xdepmod
        !
        zereli=1.0e-3_fp*(zeovab+1.5_fp*ti_nbdep(i)/aimp_nbdep(i))
        if(zimp_nbdep(i).le.6.0_fp) then
           call zstop_sigv(izb,6,zereli,1,zsigvx)
           svix2_nbdep(i)=1.0e6_fp* &
                zsigvx*zimp_nbdep(i)*zexcfac*xdepmod/6
        else if(zimp_nbdep(i).ge.8.0_fp) then
           call zstop_sigv(izb,8,zereli,1,zsigvx)
           svix2_nbdep(i)=1.0e6_fp* &
                zsigvx*zimp_nbdep(i)*zexcfac*xdepmod/8
        else
           call zstop_sigv(izb,6,zereli,1,zsigvx6)
           call zstop_sigv(izb,8,zereli,1,zsigvx8)
           zsigvx=(zsigvx6*(8-zimp_nbdep(i))+ &
                zsigvx8*(zimp_nbdep(i)-6))/2
           svix2_nbdep(i)=1.0e6_fp*zsigvx*zexcfac*xdepmod
        end if
        !
        znsvii=0
        znii=0
        if(nhs_nbdep(i).gt.0) then
           call hatom_btsigv(2,1,1.0e-3_fp*zeovab, &
                1.0e-3_fp*ti_nbdep(i)/ahs_nbdep(i),1,izb,1,zsigvii, &
                idum)
           znsvii=znsvii+1.0e6_fp*zsigvii*nhs_nbdep(i)
           znii=znii+nhs_nbdep(i)
           call hatom_btsigv(1,1,1.0e-3_fp*zeovab, &
                1.0e-3_fp*ti_nbdep(i)/ahs_nbdep(i),1,izb,1,zsigvcx, &
                idum)
           svcxh2_nbdep(i)=1.0e6_fp*zsigvcx*zexcfac*xdepmod
        else
           svcxh2_nbdep(i)=0
        end if
        if(nhes_nbdep(i).gt.0) then
           call hatom_btsigv(2,1,1.0e-3_fp*zeovab, &
                1.0e-3_fp*ti_nbdep(i)/ahes_nbdep(i),1,izb,2,zsigvii, &
                idum)
           znsvii=znsvii+1.0e6_fp*zsigvii*nhes_nbdep(i)
           znii=znii+nhes_nbdep(i)
           call hatom_btsigv(1,1,1.0e-3_fp*zeovab, &
                1.0e-3_fp*ti_nbdep(i)/ahs_nbdep(i),1,izb,2,zsigvcx, &
                idum)
           svcxhe2_nbdep(i)=1.0e6_fp*zsigvcx*zexcfac*xdepmod
        else
           svcxhe2_nbdep(i)=0
        end if
        if(nlis_nbdep(i).gt.0) then
           call hatom_btsigv(2,1,1.0e-3_fp*zeovab, &
                1.0e-3_fp*ti_nbdep(i)/alis_nbdep(i),1,izb,3,zsigvii, &
                idum)
           znsvii=znsvii+1.0e6_fp*zsigvii*nlis_nbdep(i)
           znii=znii+nlis_nbdep(i)
        end if
        svii2_nbdep(i)=znsvii*zexcfac*xdepmod/znii
        !
     end if ! end debug
     !
     zomega=0.5_fp*(omegag(j)+omegag(j+1))
     do ig=1,ng
        zdsrat=zdexp*zbnsvcx(ig)*zfplas
        sfdepcx(lsbeam)=sfdepcx(lsbeam)+zdsrat
        sorbhs_dep(j,ig,lsbeam)=sorbhs_dep(j,ig,lsbeam)-zdsrat
        sorbn0s_dep(j,ig,lsbeam)=sorbn0s_dep(j,ig,lsbeam)+zdsrat
        srbn02_dep(iz,ig,lsbeam)=srbn02_dep(iz,ig,lsbeam)+zdsrat
        sbcxdpe(j,ig,ien,lsbeam)=sbcxdpe(j,ig,ien,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_dep(j,lsbeam)=tqbcxs_dep(j,lsbeam)+ &
             1.6726e-24_fp*aplasm(ig)*gr2zc(j)*zomega*zdsrat

     end do

     !---------------------------------------------------
     !   compute tau slowing down, p.a. scattering profiles
     !   ... not for heavy neutral beams ...
     !
     if(nlheavy(lsbeam)) cycle
     !
     ioffn=6
     if(nlco(ib)) then
        iofff=0
     else
        iofff=3
     end if
     indf=ien+iofff            ! tau profile array address
     indn=indf+ioffn            ! normalizer array address
     !
     if(nlvphi) then
        zvrel=sqrt(abeam*zeovab/v2toev) ! plasma frame velocity
     else
        zvrel=vion            ! = lab frame velocity
     end if
     !
     zvrel=max(vmin(j,lsbeam),zvrel) ! nothing less than therm velocity
     !
     zcf=0.5_fp*(cfa(j,lsbeam,1)+cfa(j,lsbeam,2))
     ztausl=log((zvrel**3+vcrit3(j,lsbeam))/ &
          (vmin3(j,lsbeam)+vcrit3(j,lsbeam)))/(3.0_fp*zcf)
     !
     zef=0.5_fp*(efa(j,lsbeam,1)+efa(j,lsbeam,2))
     ztaupa=zvrel**3/(2.0_fp*zef)
     !
     taubsl(j,indf,lsbeam)=taubsl(j,indf,lsbeam)+ztausl*zdexp
     taubsl(j,indn,lsbeam)=taubsl(j,indn,lsbeam)+zdexp
     !
     taubpa(j,indf,lsbeam)=taubpa(j,indf,lsbeam)+ztaupa*zdexp
     taubpa(j,indn,lsbeam)=taubpa(j,indn,lsbeam)+zdexp
     !momentum injected
     if(.not.iltrk0) then
        if(i.eq.1)then
           zlndep = xltkbf+zlng(i)*0.5_fp
        else
           zlndep = xltkbf+(zlng(i-1)+zlng(i))*0.5_fp
        end if
        call nb0loc(zlndep,zsinth,zcosth,zsinph,zcosph, &
             zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
        bphinjs(lsbeam)=bphinjs(lsbeam)+zrloc*zvtor*rmassb*zdexp
     else
        bphinjs(lsbeam)=bphinjs(lsbeam)+zphi0*zdexp
     end if

     !
  end do

  !
  !  if dep call is for neutral density calculation only, exit now; do
  !  not try to deposit an ion.
  !
  if(idepbox) then
     call addn0(10.0_fp*xltkwl,zexp,0.0_fp,vion,idotrk,lsbeam)
     return
  else if (ideptube) then
     call tube_addn0((zlvac+zlng(nzcrtk))/vion,10.0_fp*xltkwl/vion, &
          zexp, 0.0_fp)
     return
  end if
  !
  if(iltrk0) then

     call nb0loc(zlvac+zlng(nzcrtk),zsinth,zcosth,zsinph,zcosph, &
          zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
     track_rout(kdep0,mibs2use(lsbeam))=zrloc
     track_zout(kdep0,mibs2use(lsbeam))=zyloc + plhgt
     track_phiout(kdep0,mibs2use(lsbeam))=esta + &
          nsnccw*nsjdotb*atan2(zyhor,zxhor)/degrad

     track_xl(mtrk0_prev+1:itrk0) = &
          track_xl(mtrk0_prev+1:itrk0)/ &
          track_xl(itrk0)

     if(nztkmd.gt.0) then
        ii=itrk_sav1-1
        call track_copy(itrk_sav1)
        ii=itrk_sav2+1
        call track_copy(itrk_sav2)
     end if
  end if
  !
  if(iuncbd.gt.7) then
     call infolog('dep -- track has ',iuncbd)
     call infolog(' lost orbit region bdys')
     if(myidd.eq.0) then
        write(nonlin,9001) iuncbd
     end if
9001 format(' % dep -- track has ',i2,' lost orbit region bdys')
  end if
  !
  !%      if(izork.eq.0) then
  !%        call r8_grafx2(zlng,zshf0,zstat,nzcrtk,'cm','none',
  !%     >     'transp debug test plot, run '//runid(1:lrunid),
  !%     >         'xjadep shift; (dash) orbit status',
  !%     >     '0.0 is passing, 0.1 trapped, .ge. 0.9 lost')
  !%        call r8_grafx2(zlng,zvpov,zstat,nzcrtk,'cm','none',
  !%     >     'transp debug test plot, run '//runid(1:lrunid),
  !%     >         'xjadep vpll/v; (dash) orbit status',
  !%     >     '0.0 is passing, 0.1 trapped, .ge. 0.9 lost')
  !%        call uread('dep_debug:  redo track? (y/n):$')
  !%        if(lyes(0)) goto 90
  !%      end if
  !
  !-----------------------------------------------------------
  !
  !      keep track of fraction which hits far wall
  !
  zwge=zexp*ze
  bpshins(lsbeam)=bpshins(lsbeam)+zwge
  bpshna(ib,ien)=bpshna(ib,ien)+zwge
  !
  sfshine(lsbeam)=sfshine(lsbeam)+zexp
  call nb0loc(zlvac+zlng(nzcrtk),zsinth,zcosth,zsinph,zcosph, &
       zxhor,zyhor,zrloc,zyloc,zvtor,zvr,zvy)
  bphshins(lsbeam)=bphshins(lsbeam)+zrloc*zvtor*rmassb*zexp

  !
  sfrac_survive(ien,ib)=sfrac_survive(ien,ib)+(1.0_fp-zexp)
  wfrac_survive(ien,ib)=wfrac_survive(ien,ib)+1.0_fp
  !
  !  done if depositing for fpp
  !
  if(nlfbmfpp(lsbeam).or.(idepo.eq.0)) then
     ewghtdep(lsbeam)=ewghtdep(lsbeam)+ze
     return
  end if
  !
  ewghtdep(lsbeam)=ewghtdep(lsbeam)+zwge
  !
  !-----------------------------------------------------------
  !
  !   deposit particle
  !
  !      choose initial trkrb starting zone randomly,
  !      using accumulated probability array
  !
  zrand=nbrandm_dep(idum)*zprb(nzcrtk)
  !
  zprbp=0.0_fp
  !
  do i=1, nzcrtk
     if(zrand.le.zprb(i)) then
        itrmax=i
        zprbls=zprbp
        goto 220
     end if
     zprbp=zprb(i)

  end do
  !
  call errorlog('nubeam dep -- algorithm error, zrand=',zrand)
  if(myidd.eq.0)then
     write(nonlin,*) ' ?? nubeam dep -- algorithm error, zrand=', &
          zrand
  end if
  ierr=1
  return
  !
220 continue
  !
  !      interpolate within trkrb zone and add on vacuum segments
  !  length from source to deposition point
  zfaci=(zprb(itrmax)-zrand)/(zprb(itrmax)-zprbls)
  !
  !  the bigger zfaci, the further back within the zone
  zlndep=zlng(itrmax)-xlngtk(itrmax)*zfaci+xltkbf
  !
  if(jdbg) then
     call infolog('nubeam dep --zfaci = ',zfaci)
     call infolog('zlndep = ',zlndep)
     if(myidd.eq.0) then
        write(nonlin,*) '   zfaci,zlndep: ',zfaci,zlndep
     end if
  end if
  !
  !  interpolate average xi of orbit for mc variance
  !
  zxiavl=zxiav(itrmax)
  if(zfaci.ge.0.5_fp) then
     if(itrmax.eq.1) then
        zxiavx=1.0_fp
     else
        zxiavx=zxiav(itrmax-1)
     end if
     zxiavi=zxiavl+(zfaci-0.5_fp)*(zxiavx-zxiavl)
  else
     if(itrmax.eq.nzcrtk) then
        zxiavx=1.0_fp
     else
        zxiavx=zxiav(itrmax+1)
     end if
     zxiavi=zxiavl+(0.5_fp-zfaci)*(zxiavx-zxiavl)
  end if
  !
  !  include intermediate vacuum segment in length to deposition if needed
  zl0=xltkbf
  if(itrmax.gt.nztkmd.and.nztkmd.ne.0) then
     zlndep=zlndep+xltkmd
     zlmid=xltkmd
     zl1=zlng(nztkmd)+xltkbf
     zl2=zl1+zlmid
     zl3=zlng(nzcrtk)+zlmid+xltkbf
  else
     zl1=zlng(nzcrtk)+xltkbf
     zl2=0.0_fp
     zl3=0.0_fp
     zlmid=0.0_fp
  end if
  !
  !  make sure zlndep is strictly in bounds
  if(abs(zl0-zlndep).lt.1.e-2_fp) zlndep=zl0+1.e-2_fp
  if(abs(zl1-zlndep).lt.1.e-2_fp) zlndep=zl1-1.e-2_fp
  if(abs(zl2-zlndep).lt.1.e-2_fp) zlndep=zl2+1.e-2_fp
  if(abs(zl3-zlndep).lt.1.e-2_fp) zlndep=zl3-1.e-2_fp
  !
  !%      if(nbplot.gt.0) then
  !%c  plot trajectory
  !%        do 240 ipl=1,100
  !%          zlplo=xltkbf+(ipl-1)*(zlndep-xltkbf-zlmid)/99.0_fp
  !%          if(zlplo.gt.zl1) zlplo=zlplo+zlmid
  !%          if(abs(zl0-zlplo).lt.1.e-2_fp) zlplo=zl0+1.e-2_fp
  !%          if(abs(zl1-zlplo).lt.1.e-2_fp) zlplo=zl1-1.e-2_fp
  !%          if(abs(zl2-zlplo).lt.1.e-2_fp) zlplo=zl2+1.e-2_fp
  !%          if(abs(zl3-zlplo).lt.1.e-2_fp) zlplo=zl3-1.e-2_fp
  !%          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
  !% 240    continue
  !%        zzx(1)=100.0
  !%        zzy(1)=100.0
  !%        call r8_cplot(zzx,zzy)
  !%      end if
  !
  !  get initial orbit parameters
  !
  call depsub(zlndep,zsinth,zcosth,zsinph,zcosph,inoflr, &
       zphi0,zphigc,ze0,zegc,zbt0,zbtgc)
  !
  if(jdbg) then
     call infolog('nubeam dep:  xion = ',xion)
     call infolog('               th = ',th)
     call infolog('             vion = ',vion)
     call infolog('            xksid = ',xksid)
     if(myidd.eq.0) then
        write(nonlin,'(a,4(1x,1pe17.10))') &
             '     xion,th,vion,xksid: ',xion,th,vion,xksid
     end if
  end if
  !
  !  evaluate (r,y) coords; plot guiding ctr
  !
  !%      if(nbplot.gt.0) then
  !%        call xmomry(xion,th,0,rmjion,yion,zdummy)
  !
  !%        zzx2(2)=rmjion
  !%        zzy2(2)=yion
  !%        call symbl(2)
  !%        call line(-1)
  !%        call r8_cplot(zzx2,zzy2)
  !%        call tsend
  !%      end if
  !
  !      evaluate weight of particle    (note that for consistency must
  !            use particle, not g.c. minor radius)
  !
  !  use "avg j" of first orbit
  !
  if((zstat(itrmax).gt.0.5_fp).or.(ngc.ge.lep1)) then
     xjav=lep1              ! orbit probably lost
  else
     xjav=lcentr+zxiavi*nzones      ! xjadep est. avg orbit location
  end if
  !
  if(ilwgpro) then
     zwgpro=fdep_k(zxiav(i),xwght0(lsbeam))
     !mg         zwgpro=1.0_fp-(1.0_fp-xwght0(lsbeam))*
     !mg     >        (1.0_fp-zxiav(i)**1.5_fp)
     ! mg 01/17/18
     zwgpro=1.0_fp-(1.0_fp-xwght0(lsbeam))* &
          (1.0_fp-zxiav(i)**1.5_fp)
  end if
  !
  !      average wght for a particle chosen along this chord
  !      should come out 1-zexp, in order to give equal weight
  !      to all chords
  !
  wght=zwgpro*zprb(nzcrtk)
  ewghtdep(lsbeam)=ewghtdep(lsbeam)+wght*ze
  wghtdep(lsbeam)=wghtdep(lsbeam)+wght
  !
  j=nznttk(i)
  mc_newdens(j,lsbeam) = mc_newdens(j,lsbeam) + wght
  !
  !-----------------------------------------------------------------
  !  caution: nb0jxb call increments mc sums; bphdeps is a mc sum;
  !  the normalization to be applied is the one for newly deposited
  !  particles-- ndep0 loop contribution should not be counted -- cf depnrm
  !  dmc sept 2006
  !
  call nb0jxb(wght,ze0,zegc,zphi0,zphigc,zbt0,zbtgc, &
       nflr,ngc,xiflr,xion,1,.false.)
  !
  !  momentum deposition (cf depsub)
  bphdeps(lsbeam)=bphdeps(lsbeam)+wght*zphi0  !at point of ionization
  bphdepsgc(lsbeam)=bphdepsgc(lsbeam)+wght*zphigc   !at gc
  !
  !  end caution dmc sept 2006
  !-----------------------------------------------------------------
  !
  xistate=xiflr      ! for fast ion radial current tracking
  !
  !==============================
  !
  return
  !
  !-----------------------------------
contains
  subroutine track_copy(jj)
    use iso_c_binding, only: fp => c_double

    ! copy track data into vacuum region between plasma
    ! intersection segments

    integer, intent(in) :: jj

    track_dr_flr(jj)=track_dr_flr(ii)
    track_dz_flr(jj)=track_dz_flr(ii)
    track_de_flr(jj)=track_de_flr(ii)
    track_xjadep(jj)=track_xjadep(ii)
    track_vpv_dep(jj)=track_vpv_dep(ii)
    track_dvpv_flr(jj)=track_dvpv_flr(ii)

  end subroutine track_copy

  logical function idbg_out(idum)
    integer :: idum

    idbg_out = (icount_depall.eq.icount_debug1) .or. &
         (icount_depall.eq.icount_debug2) .or. &
         (icount_depall.eq.icount_debug3) .or. &
         (icount_depall.eq.icount_debug4)

  end function idbg_out

end subroutine nbidep
