!  depbox_all -- use dep to compute fast neutral density in a 3d box aligned
!     on the axis of a select neutral beam for all generation of neutrals  except 0.
!
subroutine depbox_all(ierr)
  use nbi_com
  use physconst_mod, only : evptocm_sec
  use nbi_halo_mod_def
  use nbi_random
  use nbi_random_pll
  use nbi_mpi_com
  use adas_mod,only: adas_sigv
  use logmod
  use nbi_strline,only: nbi_track_call_halo
  implicit none
  integer :: i,intrcpt
  integer :: iseed,iscreen1,idum,iseed1
  integer :: ngeo !geometrical splitting along each track
  integer :: nkin !splitting in velocity space
  integer :: nrequest
  integer :: idep0,i_gen,navail,ipos,ierr,ier
  real(fp) :: zavail,zrem,fsplt_adjust,p_dep
  real(fp) :: ftarg ! probability p_dep for a new neutral track
  integer ::  isize_trk !size allocated tracks for halo, 'itrakch'
  integer ::  i_start  ! start search
  integer :: ii1
  real(fp) :: r_end(3) ! point how far new neutral paticle could travel, in tokamak cartisian coord., [m,m,m]
  real(fp) :: r(3) !  r is a start points of a new neutral track in tokamak car.coord, [m,m,m]
  real(fp) :: v_xyz(3) ! velocity vector in tokamak cartesian coord system [cm/sec]
  real(fp) :: v_rel(3) ! v_parent-v_chid in tokamak cartesian coord system [cm/sec]
  integer :: idtr !new track id#
  integer ::  isize_new ! new size for a track arrays, as we a goinig to add new tracks
  ! sample on track that will be splitted and produces a new track/s
  ! each cpu has subset of '0' generation tracks
  integer ::  isize_need ! needed # of tracks 
  real(fp), dimension(:), allocatable :: itrmap ! mapping track id
  integer :: mpierr,i_sol,i_dep_parent,i_depnew
  integer :: lsbeam_parent ! parent's lsbeam, for 'fast' lsbeam<=nsbeam, 
                           ! for 'THERM' lsbeam>nsbeam and ig=lsbeam-nsbeam
  real(fp),  dimension(:,:), allocatable :: bdns2p_ptcl ! fast ion density ant ptcl position
  real(fp) :: zdl_vneut, z_kappa,deltp,dtl
  real(fp) :: xin,yin,zin,xout,yout,zout,z_length
  integer :: i1,iz,inpl
  integer :: ierr_dalloc,ierr_alloc
  real(fp) :: evptocm_sec2

  evptocm_sec2=evptocm_sec*evptocm_sec
  ierr_dalloc=0
  ierr_alloc=0

  ! decide on weight of a new track
  if(myidd.eq.0)then
     write(nonlin,*)' in depbox_all'
  endif
  call enterlog("depbox_all") ! all returns through exitlog()
  iseed=abs(14111000*nbrandm(idum))
  iscreen1=1400 !doesn't matter 
  i_gen=0
  intrcpt=0
  isize_trk=size(itrackh_rqinit)
  isize_need=isize_trk

#ifdef __MPI
  htime_all=0.0_fp
  htime_at=0.0_fp
  htime_alc=0.0_fp
  htime_add=0.0_fp
  htime_map=0.0_fp
  htime_gm=0.0_fp
  htime_aim=0.0_fp
  htime_tr=0.0_fp
  htime1=mpi_wtime()
  htime7=mpi_wtime()
#endif
 
  if(allocated(itrmap)) deallocate(itrmap,stat=ierr_dalloc)
  if(ierr_dalloc.ne.0) then
     call errorlog(' depbox_all: error dealloc itrmap ',ierr_dalloc)
     call bad_exit
  endif
  
  allocate(itrmap(isize_trk),stat=ierr_alloc)
  if(ierr_alloc.ne.0) then
     call errorlog(' depbox_all: error alloc itrmap isize_trk=',isize_trk)
     call bad_exit
  endif
#ifdef __MPI
  htime2=mpi_wtime()
  htime_alc=htime_alc+htime2-htime1
#endif
  itrmap=0
  idtr=-mdepbox !deposited # of '0' generation tracks
  do
     i_start=1
     idep0=0
     isize_trk=size(itrackh_rqinit)
     if(size(itrmap).lt.isize_trk)then

#ifdef __MPI
        htime1=mpi_wtime()
#endif

        if(allocated(itrmap)) deallocate(itrmap,stat=ierr_dalloc)
        if(ierr_dalloc.ne.0) then
           call errorlog(' depbox_all: error dealloc itrmap i_start',i_start)
           call bad_exit
        endif
        allocate(itrmap(isize_trk),stat=ierr_alloc)
        if(ierr_alloc.ne.0) then
           call errorlog(' depbox_all: error alloc itrmap i_start',i_start)
           call bad_exit
        endif

#ifdef __MPI
        htime2=mpi_wtime()
        htime_alc=htime_alc+htime2-htime1
#endif

     endif
     itrmap=0

     !find how many next generation will start search from a parent i_start, see comment bellow 
     do while(i_start.le.isize_trk) ! avoid dealing with not allocated arrays
        if(.not.itrackh_rqinit(i_start)) then
           intrcpt=itrackh(i_start)%izone_pl(2) !if >0 neutral passed through plasma region
           if(intrcpt.gt.0)then
              if(i_gen.eq.itrackh(i_start)%idgen)then
                 idep0=idep0+1
                 itrmap(idep0)=i_start
              endif
           endif
        else
           exit !this track and following ones are not initialized yet
        endif
        i_start=i_start+1 !i_start
     enddo

     if(idep0.eq.0)then
        if(allocated(itrmap)) deallocate(itrmap,stat=ierr_dalloc)
        if(ierr_dalloc.ne.0) then
           call errorlog(' depbox_all: error dealloc itrmap idep0=0')
           call bad_exit
        endif
        exit !all done
     endif
     call warnlog(' depbox_all: generation of neutrals i_gen=',i_gen)
     call warnlog(' depbox_all: # of tracks ',idep0)
     i_gen=i_gen+1 !start a new generation of neutrals
     
     call param_fbm_ptcl !get bdns2p_ptcl at ptcl position
 
     i_dep_parent=0
     i_depnew=0
     isize_new=itrmap(idep0) !already taken # of spots

     do ii1=1,idep0 !cycle over all 'parent' tracks
        i=itrmap(ii1)
        iseed1=iseed+i+idep0*myidd
        call nbi_random_init_pll(iseed,iseed1,iscreen1)
        zavail=itrackh(i)%weight/wmin_halo
        navail=zavail
        zrem=zavail-navail

        if(nbrandm_pll(idum).lt.zrem) navail=navail+1

        if(navail.eq.0) then
           cycle  !weight.le.wmin_halo  => ending this track -- has been ionized
        endif

        intrcpt=itrackh(i)%izone_pl(2) ! box zone idex when neutral exits plasma
        xout=itrackh(i)%ptb(1)
        yout=itrackh(i)%ptb(2)
        zout=itrackh(i)%ptb(3)
        xin=itrackh(i)%pta(1)
        yin=itrackh(i)%pta(2)
        zin=itrackh(i)%pta(3)
        z_length=sqrt((xout-xin)**2+(yout-yin)**2+(zout-zin)**2) !neutral track length 
        inpl=itrackh(i)%izone_pl(1)
        if(inpl.eq.0) then
           cycle ! this track did not enter plasma region
        endif
        
        if(itrackh(i)%pdep(intrcpt).eq.1.0_fp) then
           !this neutrals has not been ionized, no reason to make
           !neutral halo from it
           cycle
        endif
        if(itrackh(i)%xlngtk(intrcpt).gt.z_length) then
           !adjust this track exit point
           if(intrcpt.eq.1)then
              call infolog(' depbox_all: #of intrcpt =1 for track#',i) 
           else
              i1=0
              do iz=intrcpt-1,inpl,-1
                 if(itrackh(i)%xlngtk(iz).le.z_length) then
                    i1=iz
                    exit
                 endif
              enddo
              if(i1.eq.0) cycle !something wrong with this track, don't mess with it
              !otherwise reset
              if(itrackh(i)%xlngtk(iz).eq.z_length) then
                 intrcpt=i1
              else
                 intrcpt=i1+1
              endif
           endif
           itrackh(i)%izone_pl(2)=intrcpt ! box zone idex when neutral exits plasma
           itrackh(i)%xlngtk(intrcpt)=z_length
        endif
        
         
!find deposition probabilty at the end of the box 
        p_dep=1.0_fp 

        if(itrackh(i)%pdep(intrcpt).gt.0.0_fp)then
           
           z_kappa=itrackh(i)%dkappa(intrcpt)
           if(intrcpt.gt.1) then
              zdl_vneut=(itrackh(i)%xlngtk(intrcpt)-itrackh(i)%xlngtk(intrcpt-1))/itrackh(i)%vneut
           else
              zdl_vneut=itrackh(i)%xlngtk(intrcpt)/itrackh(i)%vneut
           endif
           dtl = p_delt(zdl_vneut,z_kappa,deltp)                
           p_dep=1._fp-itrackh(i)%pdep(intrcpt)*deltp !deposition probabilty at the last intercept
        endif
        if((nsplt_gen.gt.0).and.(i_gen.gt.nsplt_gen)) then
           ngeo=1
           nkin=1
        else
           ngeo=max(1,int((nsplt_geo+0.5_fp)*p_dep))
           nrequest=ngeo*nsplt_kin
           fsplt_adjust=min(1._fp,real(navail,fp)/real(nrequest,fp)) ! to reduce splitting as track weight --> wmin_halo
           ngeo=max(1,int(fsplt_adjust*ngeo)) ! so many splittings allowed along 'i'th track
           nkin=max(1,int(fsplt_adjust*nsplt_kin))! so many splittings in 'v' space allowed 
        endif

        ! use rng to sample distribution{p} to choose launch point{x} for a new track
        call nbi_halo_random_fetch(ngeo,ierr)
        if(ierr.ne.0) then
           call errorlog(' depbox_all: nbi_halo_random_fetch failed')
           call bad_exit
        endif

        ipos=1
        isize_new = isize_new + ngeo*nkin
        if(isize_trk.lt.isize_new) then !check if there are enough space for a new tracks
#ifdef __MPI
        htime1=mpi_wtime()
#endif
        ! need extra space
        ! make it with some extras
           isize_need = isize_new + (idep0-ii1)*ngeo*nkin
           
           call nbi_rqinit_expand(isize_need) ! more tracks are needed
           isize_trk = isize_need !updated # of tracks
#ifdef __MPI
        htime2=mpi_wtime()
        htime_alc=htime_alc+htime2-htime1
#endif
        endif
        i_sol=0 ! for debug purposes
        i_dep_parent=i_dep_parent+1 !# of parent's tracks used for a new track deposition
        do
           if(ipos.gt.ngeo) exit
           ftarg=drndm(ipos)
           call nbi_halo_new_tr(i,ier)
           if(ier.ne.0) then 
              call errorlog(' depbox_all: nbi_halo_new_tr error')
              ierr=ier
              call bad_exit
           endif
           ipos=ipos+1
        enddo
     enddo
     call infolog(' depbox_all: #of new tracks deposited=',i_depnew)
     call infolog('             from #of parent tracks=',i_dep_parent)
     ! new generation neutrals tracks have been started
     ! find sirvival probability  for these tracks
  enddo

#ifdef __MPI
  htime1=mpi_wtime()
#endif

  if(allocated(itrmap)) deallocate(itrmap,stat=ierr_dalloc)
  if(ierr_dalloc.ne.0) then
     call errorlog(' depbox_all: error dealloc itrmap at end')
     call bad_exit
  endif
  if(allocated(bdns2p_ptcl)) deallocate(bdns2p_ptcl,stat=ierr_dalloc)
  if(ierr_dalloc.ne.0) then
     call errorlog(' depbox_all: error dealloc bdns2p_ptcl')
     call bad_exit
  endif

#ifdef __MPI
  htime2=mpi_wtime()
  htime_alc=htime_alc+htime2-htime1
#endif
#ifdef __MPI
  htime8=mpi_wtime()
  htime_all=htime_all+htime8-htime7
#endif
  call warnlog(' depbox_all: time for allocation = ', htime_alc)
  call warnlog(' depbox_all: time spent in addn0 = ', htime_add)
  call warnlog(' depbox_all: time spent in mapping = ', htime_map)
  call warnlog(' depbox_all: time spent in atomic phys= ', htime_at)
  call warnlog(' depbox_all: time spent in geom= ', htime_gm)
  call warnlog(' depbox_all: time spent in nbi_aimingangle= ', htime_aim)
  call warnlog(' depbox_all: time spent in nbi_track_call_halo= ', htime_tr)
  call warnlog(' depbox_all: all time spent depbox_all= ', htime_all)

  call exitlog("depbox_all") ! all returns through exitlog()
  if(myidd.eq.0)then
     write(nonlin,*)' exit depbox_all'
  endif
  

  
contains
  
  subroutine nbi_halo_new_tr(itr,ier)
    integer, intent(in) :: itr ! 'parent'track id
    integer, intent(out) :: ier 
    integer :: i1,i2, iz,is
    real(fp) :: p1
    real(fp) :: p2
    real(fp) :: zfac,zfac0,zfac1,frac
    integer :: ii,ib,ien
    real(fp) ::  xx, yy
    real(fp) ::  rb(3),rb0(3),r0(3)! rb - end point in box, rb0 - start point in box , for a parent
                                 ! so  rb  is a start point for a new generation track
                                 ! r0 and r - start and end points in tokamak car.coord, for a parent
                                 ! r is a start point of a new neutral track in tokamak coord system
    real(fp) :: v_ptcl , zr,zphi,v_tor,v_rot
    integer :: n_j,n_iz !transp 1d zone and beam 2d zone
    real(fp) :: zcxfr,totcx,zsp,asp,zzn,zzns,zznsp
    real(fp) :: zs_pll
    integer :: isb,irange_bb !fast particle indexing and  table energy range sector
    integer :: jba_0,lgyro_neut !pitch angle  and gyro  indices
    logical :: ilinside_pl
    character*5 :: i_type !species type: 'FASTB','FASTF','THERM'
    integer :: ie,ia,igy !energy, pitch angle, aiming angle zones #, correspondingly
    integer :: ig ! thermal species index
    real(fp) :: zangl ! angle between v_perp and normal to flux surface
    real(fp) :: zvec(3) ! velocity vector of new  particle of a given species
                      ! fast species: in tokamak cyl. coordinates [m/sec]
                      ! thermal species: in tokamak cartesian coordinates [m/sec]
    real(fp) :: zvecb(3)! unit velocity vector of new  particle in box cartesian coordinates [1/sec]
    real(fp) :: rb_end(3) ! point how far paticle could travel, in box cartisian coord.
    integer :: idbbx(3) !work array here
    real(fp) :: za ! pitch angle
    real(fp) :: v_pll,v_perp,vi2
    real(fp) :: alpha_v !aiming angle
    real(fp) :: zlen,zvec2
    integer :: irng ! 1st guess to range # of energy grid in nbeloc.
    real(fp) :: ze_rel(1),zsigv_cx(1) !relative energy/amu, sig*v cx
    real(fp) :: fkin !kinetic adjustment to the weight
    real(fp) :: dens_idtr !density of species assigned for a new track
    integer :: ind,indp1 !energy grid index
    real(fp) :: zintri,zintrp
    real(fp) :: zeovab ! energy / nucleon of neutral in plasma frame, [ev/amu]
    real(fp) :: zbnsvis,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zexcfac,zbnsvcx(mig) !n*sig*v of reactions neutrals(fast,thermal) + thermal ion
    real(fp) :: zbbnsva(2,mibs) !n*sig*v of reactions neutrals(fast,thermal) + fast ion
    real(fp) :: zbmsexc1d, zsigvnt1d !n*sig*v cx 
    real(fp), dimension(ng+nsbeam)  :: fcx_sp
    integer :: inpl
    real(fp) :: zdl,xl_start
!mg    real(fp) :: zrho !start rho for a new neutral track
!mg    real(fp) :: vphi_neut ! v tor for therm neutrals
    !take a cleanup
!mg    vphi_neut=0._fp
    zbnsvis=0._fp
    zbnsvii=0._fp
    zbnsvie=0._fp
    zbnsviz=0._fp
    zbnsvto=0._fp
    zexcfac=1._fp
    zbnsvcx(:)=0._fp
    zbbnsva(:,:)=0._fp
    zbmsexc1d=0._fp
    zsigvnt1d=0._fp
    fcx_sp=0._fp
    ier=0
    i1=0
    inpl=itrackh(itr)%izone_pl(1)
    p2=1.0_fp-itrackh(itr)%pdep(inpl) !deposition probability
    ftarg=p_dep*(1.0_fp-ftarg)+p2*ftarg !normalized to the [p_dep, p2] interval
    do iz=intrcpt,inpl,-1
       p2=1.0_fp-itrackh(itr)%pdep(iz) !deposition probability
       if(ftarg.ge.p2) then
          i1=iz
          exit
       endif
    enddo
    if(i1.eq.0) then
       call errorlog(' depbox_all: i1=0, ftarg= ',ftarg)
       call errorlog(' depbox_all: i1=0, p_dep= ',p_dep)
       call errorlog(' depbox_all: i1=0, p2= ',p2)
       do iz=intrcpt,inpl,-1
          p2=1.0_fp-itrackh(itr)%pdep(iz) !deposition probability
          call errorlog(' depbox_all: iz= ',iz)
          call errorlog(' depbox_all: 1-p(iz)= ',p2)
       enddo
       ier=1
       call bad_exit
    endif
    zdl=0.0_fp
    if(i1.eq.1) then
       xl_start=0.0_fp
    else
       xl_start=itrackh(itr)%xlngtk(i1-1)
    endif
    if(ftarg.ne.p2) then
       zdl=log((1.0_fp-p2)/(1.0_fp-ftarg))*itrackh(itr)%vneut/itrackh(itr)%dkappa(i1)
    endif
    zfac=xl_start+zdl
    if(zfac.gt.itrackh(itr)%xlngtk(intrcpt)) zfac=itrackh(itr)%xlngtk(intrcpt)

!mg    p2=1.0_fp-itrackh(itr)%pdep(inpl) 
!mg    if(intrcpt.gt.inpl) then
!mg!mg1       p1=1.0_fp-itrackh(itr)%pdep(intrcpt) 
!mg       p1=p_dep 
!mg    else
!mg       p1=p2
!mg    endif
!mg    if(ftarg.le.p2) then !first zone in plasma
!mg       i1=inpl
!mg       i2=inpl
!mg       zfac=0.0_fp
!mg       if(p2.gt.0.0_fp) then
!mg          frac=ftarg/p2
!mg          if(inpl.eq.1)then
!mg             zfac=itrackh(itr)%xlngtk(inpl)*frac ! uniformly distributed along first zone
!mg          else
!mg             zfac=itrackh(itr)%xlngtk(inpl)*frac+itrackh(itr)%xlngtk(inpl-1)*(1.0_fp-frac) ! uniformly distributed along first zone
!mg          endif
!mg       endif
!mg    else if(ftarg.ge.p1) then !shine-trough particle, not supposed to be here
!mg       if(ftarg.gt.p1) then
       ! not supposed to be there, check mapping
!mg          call errorlog(' emergency exit depbox_all: shine-trough particle, not supposed to be here')
!mg          call errorlog(' depbox_all: ftarg.ge.p1, ftarg= ',ftarg)
!mg          ier=1
!mg#ifdef __MPI
!mg          call mpi_finalize(mpierr)
!mg#endif
!mg          return
!mg       endif
!mg       i1=intrcpt
!mg       i2=intrcpt
!mg       frac=ftarg
!mg       if(intrcpt.eq.1)then
!mg         zfac=itrackh(itr)%xlngtk(intrcpt)*frac ! uniformly distributed along last zone
!mg       else
!mg          zfac=itrackh(itr)%xlngtk(intrcpt)*frac+itrackh(itr)%xlngtk(intrcpt-1)*(1.0_fp-frac) ! uniformly distributed along last zone
!mg       endif
!mg    else
!mg       do iz=inpl+1,intrcpt
!mg          p1=p2
!mg          p2=1.0_fp-itrackh(itr)%pdep(iz)
!mg          if(ftarg.le.p2) then
!mg             i1=iz-1
!mg            i2=iz
!mg             frac=(ftarg-p1)/(p2-p1)
!mg             exit
!mg          endif
!mg       enddo
!mg       zfac1=itrackh(itr)%xlngtk(i2) !track length
!mg       zfac0=itrackh(itr)%xlngtk(i1)
!mg       zfac=zfac1*frac+zfac0*(1._fp-frac) !linear approximation; not the best 
!mg    endif

    ! in [m]
    if(z_length.gt.0.0_fp)then
       frac=zfac/z_length
       if(frac.le.0.0_fp.or.frac.gt.1.0_fp) then
          call errorlog(' emergency exit depbox_all: frac = 0 or >1')
          call errorlog('  ftarg =',ftarg)
          call errorlog('  p2 =',p2)
          call errorlog('  itr=',itr)
          call errorlog('  intrcpt=',intrcpt)
          call errorlog('  i1=',i1)
          call errorlog('  xlngtk(i1) = ',itrackh(itr)%xlngtk(i1))
          ier=1
          call bad_exit
       endif
    else
       ! not supposed to be there, check mapping
       call errorlog(' emergency exit depbox_all: track has zero length')
       call errorlog(' track has zero length, generation id =',i_gen)
       call errorlog(' track has zero length, itr =',itr)
       call errorlog(' track length',itrackh(itr)%xlngtk(itrackh(itr)%maxpts))
       call errorlog(' track zone in plasma at entrance',inpl)
       call errorlog(' track zones in plasma at exit',intrcpt)
       ier=1
       call bad_exit
    endif

    ! start coordinates for a new track {x_l} [m], in box cartesian
    rb(1) = xin*(1._fp-frac) + xout*frac
    rb(2) = yin*(1._fp-frac) + yout*frac
    rb(3) = zin*(1._fp-frac) + zout*frac

#ifdef __MPI
  htime5=mpi_wtime()
#endif
   call box2car(itrackh(itr)%idbox,rb*100.0_fp,r) ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm] 
#ifdef __MPI
  htime6=mpi_wtime()
  htime_map=htime_map+htime6-htime5
#endif
    r(:)=r(:)*0.01_fp !cm -> m, in tokamak cartisian coordinates

    ! get fast particle velocity from the  parent
    !find unit vector in r direction
    rb0(1) = xin
    rb0(2) = yin
    rb0(3) = zin
#ifdef __MPI
  htime5=mpi_wtime()
#endif
  call box2car(itrackh(itr)%idbox,rb0*100.0_fp,r0) ! map  x,y,z coordinates of point in box coordinates to tokamak, [m] ->[cm] 
#ifdef __MPI
  htime6=mpi_wtime()
  htime_map=htime_map+htime6-htime5
#endif
    r0=r0*0.01_fp !cm -> m, in tokamak cartisian coordinates
    rb0=r-r0
    rb0=rb0/dsqrt(sum(rb0*rb0)) ! unit vector, in tokamak cartisian coordinates
    
    v_ptcl=itrackh(itr)%vneut*100._fp ! total neutral velocity [cm/sec]
    
    v_xyz = abs(v_ptcl)*rb0 !v_x, v_y, v_z in tokamak cartesian coord system [cm/sec]
    !convert tokamak cartesian coord system to tokamak cyl. coord system
#ifdef __MPI
        htime11=mpi_wtime()
#endif
    call eq_rcyl(nsnccw*nsjdotb,r(1),r(2),zr,zphi) ![m,m] -> [m,rad]
    ! find transp zone n_j, and beam zone n_iz
    call nbi_zone_ind(zr,zphi,r(3),n_j,n_iz)
#ifdef __MPI
        htime12=mpi_wtime()
        htime_gm=htime_gm+htime12-htime11
#endif
    if(n_j.eq.lep1) then
       !sol, neutral cannot cx there, as consequences
       !cannot produce new neutral tracks
       i_sol=i_sol+1
       call infolog(' nbi_halo_new_tr: sol, neutral cannot produce new tracks there')
       call infolog(' nbi_halo_new_tr: #of entries in sol i_sol=',i_sol)
       call infolog(' nbi_halo_new_tr: #of geo splittings ngeo=',ngeo)
       return
    endif

    ilinside_pl=(n_j.lt.lep1.and.n_j.ge.lcentr.and.n_iz.gt.0.and.&
         &        n_iz.le.mimxbzf)
    if(.not.ilinside_pl) then !!track not in plasma
       ! not supposed to be there, check mapping
       call errorlog(' emergency exit depbox_all: nbi_halo_new_tr error')
       call errorlog(' nbi_halo_new_tr: nbsteps=',nbsteps)
       call errorlog(' nbi_halo_new_tr: track missed transp zone',n_j)
       call errorlog(' nbi_halo_new_tr: or track missed beam zone',n_iz)
       call errorlog(' nbi_halo_new_tr: r[m] ',zr)
       call errorlog(' nbi_halo_new_tr: phi[rad] ',zphi)
       call errorlog(' nbi_halo_new_tr: z[m] ',r(3))
       call errorlog(' nbi_halo_new_tr: track id ',itr)
       call errorlog(' nbi_halo_new_tr: # of zone intercepts ftarg= ',ftarg)
       call errorlog(' nbi_halo_new_tr: # of zone intercepts intrcpt= ',intrcpt)
       call errorlog(' nbi_halo_new_tr: new track set between zone i1= ',i1)
       call errorlog(' nbi_halo_new_tr: new track set between zone i2= ',i2)
       ier=1
       call bad_exit
    endif
    
    if(itrackh(itr)%idgen.eq.0) then
       ! this is 'FASTB' type, find and set lsbeam
       ib=nbsbox(itrackh(itr)%idbox)
       ien=nbebox(itrackh(itr)%idbox)
       do isb=1,nsbeam
          if(lsbeama(ib).eq.isb)exit
       enddo
       lsbeam=isb
       ! beam fast neutrals("0" gen) are interacted with fast ions
       ! find index for a fast ion 
       do isb=1,nsbeam
          if(bdns2p_ptcl(n_iz,isb).gt.0.0_fp) then  !check if there is some 
                                                    !fast ion density in this beam zone
             zbbnsva(1,isb)=sn00i2(ib,ien,n_iz,isb)
             zbbnsva(2,isb)=sn00x2(ib,ien,n_iz,isb) !cx
          endif
       enddo
    else
       !'gen' >'0'
       if(itrackh(itr)%type.eq.'THERM') then
          !tehrmal neutral with fast ion interaction

          call get_lsbeam(itr,lsbeam)

          ig=lsbeam-nsbeam !thermal species index

          do isb=1,nsbeam
             !thermal neutrals interact with fast ion
             if(bdns2p_ptcl(n_iz,isb).gt.0.0_fp) then  !check if there is some 
                !fast ion density in this beam zone
                zbbnsva(1,isb)=snb0i2(n_iz,ig,isb) !ii
                zbbnsva(2,isb)=snbn02(n_iz,ig,isb) !cx
             endif
          enddo
          ! i have to set lsbeam even this is thermal neutral. this action is needed to
          ! find atomic phys. data for reactions between neutral and thermal plasma. such tables
          ! are computed for an eneregy in e/amu in nbatom.for (arrays bnsv...)
          ! 
       else
          ! fast neutrals(gen>"0") are interacted with fast ions
          ! find and set lsbeam

          call get_lsbeam(itr,lsbeam)

          irange_bb=2 !for beam fast ptcl.
          jba_0=itrackh(itr)%ipitch
          lgyro_neut=itrackh(itr)%igyro
          ! in nb atomic physics tabl
          ! fast neutral (gen.>0) interact with a fast ion
          ! need to know energy of neutral in lab frame 
          ! (see cxnsum.for)
#ifdef __MPI
          htime9=mpi_wtime()
#endif
          call cxnsum_intrp(v_ptcl*v_ptcl/evptocm_sec2,-1,irange_bb,ind,indp1,zintrp) ! lookup interaction energy index 
#ifdef __MPI
          htime10=mpi_wtime()
          htime_at=htime_at+htime10-htime9
#endif
          if(itrackh(itr)%zspecies.eq.1)then !h-isotopes
             do isb=1,nsbeam
                if(bdns2p_ptcl(n_iz,isb).gt.0.0_fp) then  !check if there is some 
               
                   zbbnsva(1,isb)=(1.0_fp-zintrp)*svn3vii_h(ind,jba_0,lgyro_neut,n_iz,isb)+&
                        zintrp*svn3vii_h(indp1,jba_0,lgyro_neut,n_iz,isb)
                
                   zbbnsva(2,isb)=(1.0_fp-zintrp)*svn3vcx_h(ind,jba_0,lgyro_neut,n_iz,isb)+&
                        zintrp*svn3vcx_h(indp1,jba_0,lgyro_neut,n_iz,isb)
                endif
                
             enddo
          elseif(itrackh(itr)%zspecies.eq.2)then !he-isotopes
             
             do isb=1,nsbeam
                if(bdns2p_ptcl(n_iz,isb).gt.0.0_fp) then  !check if there is some 
               
                   zbbnsva(1,isb)=(1.0_fp-zintrp)*svn3vii_he(ind,jba_0,lgyro_neut,n_iz,isb)+&
                        zintrp*svn3vii_he(indp1,jba_0,lgyro_neut,n_iz,isb)
                
                   zbbnsva(2,isb)=(1.0_fp-zintrp)*svn3vcx_he(ind,jba_0,lgyro_neut,n_iz,isb)+&
                        zintrp*svn3vcx_he(indp1,jba_0,lgyro_neut,n_iz,isb)
                endif
                
             enddo
          endif
          
       endif
    endif
    lsbeam_parent=lsbeam
    
    ! use plasma frame for an interaction with a thermal neutrals
    ! find ptcl. energy in plasma frame,
    ! zeovab is energy of neutral in plasma frame
    call nbe_rel(zeovab,v_xyz,zr*100.0_fp,zphi,n_j)! plasma rotation is taken into account

    irng=1
    
#ifdef __MPI
    htime9=mpi_wtime()
#endif
    call nbeloc(zeovab,irng,ind,zintri,indp1,zintrp) ! lookup interaction energy index in nb atomic physics tabl,
                                                     ! lsbeam in nbi_com already
    call nbains(zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zexcfac,&
         &           zbmsexc1d, zsigvnt1d,n_j,ind,zintri,indp1,zintrp)
#ifdef __MPI
    htime10=mpi_wtime()
    htime_at=htime_at+htime10-htime9
#endif
    !the total cx fraction
    fcx_sp=0.0_fp
    do ii=1,nsbeam
       fcx_sp(ii)=fcx_sp(ii)+zbbnsva(2,ii) ! cx with fast ions
    enddo
    do ii=1+nsbeam,ng+nsbeam
       fcx_sp(ii)=fcx_sp(ii)+zbnsvcx(ii-nsbeam)*zbnsvto ! cx with thermal
    enddo
    
    zcxfr=(sum(zbbnsva(2,1:nsbeam))+sum(zbnsvcx(1:ng))*zbnsvto)/&
         &        (sum(zbbnsva(1:2,1:nsbeam))+zbnsvto) !fraction of cx to total
    totcx=sum(fcx_sp)
    if(totcx.gt.0.0_fp) then
       fcx_sp=fcx_sp/totcx !fraction of cx by species
    else
       ! not supposed to be there, check atomic phys. tables
       call errorlog(' emergency exit depbox_all: nbi_halo_new_tr error')
       call errorlog(' nbi_halo_new_tr: cx fraction = 0')
       ier=1
       call bad_exit
    endif

    ! assign  species for a new track
    zzn=nbrandm_pll(idum)
    zzns=0.0_fp
    do is=1,nsbeam+ng
       if(fcx_sp(is).le.0.0_fp) cycle
       zznsp=zzns
       zzns=zzns+fcx_sp(is)
       if((zznsp.le.zzn).and.(zzns.gt.zzn)) then
          if(is.le.nsbeam) then
             ! sample {v} from fast ion distribution f_i 
             if(nlfprod(is))then
                i_type='FASTF' !fusion fast ptcl.
             else
                i_type='FASTB' !beam fast ptcl
             endif
             zsp=xzbeams(is)
             asp=abeams(is)
             dens_idtr=bdns2p(n_iz,is)/bmvol(n_iz)*1.0e6_fp
          else
             i_type='THERM'
             zsp=backz(is-nsbeam)
             asp=aplasm(is-nsbeam)
             dens_idtr=rhob(n_j,is-nsbeam,2)* 1.0e6_fp
          endif
          lsbeam=is !set lsbeam in nbi_com
          exit
       endif
    enddo
    do ii=1,nkin
       i_depnew=i_depnew+1 ! # of new tracks deposited, for statistic
       idtr=idtr+1
!indxs(ipos)+ngeo*(ii-1) !id for a new track
       if(itrackh_rqinit(idtr)) then
          
#ifdef __MPI
  htime1=mpi_wtime()
#endif
          itrackh(i)%guard = 0
          call nbi_halo_init(idtr,0,itrackh_rqinit(idtr)) ! out itrackh_rqinit(idtr)=.false.
#ifdef __MPI
  htime2=mpi_wtime()
  htime_alc=htime_alc+htime2-htime1
#endif
       endif
       itrackh(idtr)%idbox = itrackh(itr)%idbox !box id
       itrackh(idtr)%idgen = i_gen !generation
       itrackh(idtr)%pta = rb !start ponits in box coord, [m]
       itrackh(idtr)%weight = zcxfr*itrackh(itr)%weight&
            &*(1._fp-itrackh(itr)%pdep(intrcpt))/ngeo/nkin ! adjust weight to cx fraction
!mg             &/ngeo/nkin ! adjust weight to cx fraction
       itrackh(idtr)%zspecies = zsp
       itrackh(idtr)%aspecies = asp
       itrackh(idtr)%type = i_type
       
       if(i_type.eq.'THERM') then
          
          ! sample {v} from ion distribution f_i (maxwellian - for 'therm')
! the velocity components are choosen to yield a maxwellian 
! v_ti=sqrt(t_i/mi), there is should not be factor of sqrt(2) there.
! but factor v2toev=asp/1.915e12_fp i used to define zs_pll has this 2. as
! e=v**2 *v2toev
! mg 09/26/2014  fix this bug:           zs_pll=sqrt(ti(n_j,2)*1.915e12_fp/asp)
          zs_pll=sqrt(ti(n_j,2)*1.915e12_fp/asp/2.0_fp)
          !zvec - velocity vector for a given species in tokamak cartisian coord, [m/sec]
          zvec(1)=nbgauss_pll(zs_pll,0.0_fp)*0.01_fp 
          zvec(2)=nbgauss_pll(zs_pll,0.0_fp)*0.01_fp
          zvec(3)=nbgauss_pll(zs_pll,0.0_fp)*0.01_fp
!mg!vphi_neut is needed for a nbi_omg0mn 
!mg          zrho=sqrt(r(1)*r(1)+r(2)*r(2))
!mg          if(zrho.gt.0.0_fp) then
!mg             vphi_neut=(-r(2)*zvec(1)+r(1)*zvec(2))/zrho
!mg          endif
! take into account rotation
          v_rot=0.5_fp*(omegag(n_j)+omegag(n_j+1))*zr !plasma rotation,[m/sec]
!mg 01/07/15!welcome to a lab frame
          zvec(1)=zvec(1)-v_rot*sin(zphi)
          zvec(2)=zvec(2)+v_rot*cos(zphi)
          zvec2 = sum(zvec(:)*zvec(:))
          if(zvec2.gt.0.0_fp) then
             zvec(:)=zvec(:)/dsqrt(zvec2) !unit vector
             zlen = 5*max(0.01_fp*raxis,dsqrt(r(1)*r(1)+r(2)*r(2))) ![m]
             r_end(:)=r(:)+zvec(:)*zlen ![m,m,m]
          else
             r_end(:)=r(:)
          endif
       else
          
          ! sample velocity vector from fbm_ptcl
          call nbi_fastv_sample(is,n_iz,ie,ia,igy)
          ! save pitch angle and gyro angle indx
          itrackh(idtr)%ipitch = ia
          itrackh(idtr)%igyro = igy
          !
          ie=min(ie,nznbmea(is))
          ia=min(ia,nznbma)
          igy=min(igy,mgyro)
          vi2=(efbmb(ie,is)+efbmb(ie+1,is))*evptocm_sec2/2.0_fp ![cm**2/sec**2]
          za = -1.0_fp+ 2.0_fp*(ia-0.5_fp)/nznbma !pitch angle
          v_pll=za*dsqrt(vi2) !v || b component of velocity *|v|, |v| = dsqrt(efbmb(ie,is)), [cm/sec]
          v_perp=dsqrt(max(0.0_fp,vi2-v_pll*v_pll))! [cm/sec]
          !v_pll [cm/sec]
          !v_perp [cm/sec]
          zangl=twopi/mgyro*(igy-1) !get angle
#ifdef __MPI
        htime13=mpi_wtime()
#endif
          call nbi_aimingangle(zr*100._fp,r(3)*100._fp,zangl,v_pll,v_perp,zvec,alpha_v)!zvec in tokamak cyl coordinate
#ifdef __MPI
        htime14=mpi_wtime()
        htime_aim=htime_aim+htime14-htime13
#endif
          zvec = zvec * 0.01_fp !zvec in [m/sec]
          ! stop point in tokamak cyl coordinate, delt  is length of preceding beam timestep
          ! stop point means how far ptcl. could travel 
          zvec2 = sum(zvec(:)*zvec(:))
          if(zvec2.gt.0.0_fp) then
             zvec(:)=zvec(:)/dsqrt(zvec2) !unit vector
             zlen = 5*max(0.01_fp*raxis,zr)
#ifdef __MPI
        htime11=mpi_wtime()
#endif
             call eq_xcyl(1,zr+zvec(1)*zlen,zphi+zvec(2)*twopi*zlen/zr,xx,yy)!in tokamak cyl. to cartesian
#ifdef __MPI
        htime12=mpi_wtime()
        htime_gm=htime_gm+htime12-htime11
#endif
             r_end(:)=(/xx,yy,r(3)+zvec(3)*zlen/) !in tokamak cartesian coordinate, [m]
             ! zvec cyl.coord to cartesian coord. in tokamak coord. system
             if(zlen.gt.0.0_fp)then
                zvec(:)=(r_end(:)-r(:))/zlen !unit vector
             else
                zvec(:)=0.0_fp
             endif
          else
             r_end(:)=r(:)
          endif
          
       endif
#ifdef __MPI
       htime5=mpi_wtime()
#endif
       call car2box(itrackh(idtr)%idbox,r_end*100.00_fp,rb_end)
#ifdef __MPI
       htime6=mpi_wtime()
       htime_map=htime_map+htime6-htime5
#endif
       ! rb_end is an aproximated end point of track, it must be
       ! recalculated and reset at neutral deposition with account of atomic physics.
       ! i set it now in order to be able to calculate velocity vector laterget unit vector 
       itrackh(idtr)%ptb=rb_end*.01_fp ![cm] to [m]
       itrackh(idtr)%vneut = dsqrt(zvec2) !save neutral's velocity, [m/s]
       !
       !find kinetic adjustment fkin to the weight
       ! fkin=sig(|v_0-v_i|)*|v_0-v_i|/<sig*v> for cx
       !  have to now  relative energy
!mg as fkin =1
!mg       v_rel(:)=v_xyz(:)-zvec(:)*dsqrt(zvec2)*100._fp !v_parent-v_child in tokamak cartesian coord. system , lab frame
!mg      ze_rel(1)=sum(v_rel(:)*v_rel(:))/evptocm_sec2/1.0e+3_fp   !  kev/amu

!mg       if((lev_nbidep.eq.1).or.(lev_nbidep.eq.3)) then !preact ground state atomic physics 

!mg          call hatom_sigv(1,ze_rel,1,int(itrackh(itr)%zspecies),int(itrackh(idtr)%zspecies),zsigv_cx,ierr)

!mg       elseif((lev_nbidep.eq.2).or.(lev_nbidep.eq.4)) then !adas ground state atomic physics
!mg#ifdef __MPI
!mg          htime9=mpi_wtime()
!mg#endif
!mg          call adas_sigv(1,ze_rel,1,itrackh(itr)%zspecies,itrackh(idtr)%zspecies,zsigv_cx,ierr)
!mg#ifdef __MPI
!mg          htime10=mpi_wtime()
!mg          htime_at=htime_at+htime10-htime9
!mg#endif

!mg       endif
!mg end of: as fkin =1

       fkin=1.0_fp
!mg       if(nkin.gt.1.or.i_gen.gt.1) then
!mg          if(fcx_sp(lsbeam).gt.0.0_fp) then
!mg             fkin=zsigv_cx(1)*dens_idtr/fcx_sp(lsbeam)/totcx
!mg         endif
!mg       endif
!mg changes 12/20/2013
!mg       
       itrackh(idtr)%weight=itrackh(idtr)%weight*fkin
       ! find intersection with the irregular 2d mc grid       

       call dep_neutrals(ier)
       
       if(ier.ne.0) then
          call errorlog(' emergency exit depbox_all: nbi_halo_new_tr error')
          call errorlog(' nbi_halo_new_tr: cx fraction = 0')
          ier=1
          call bad_exit
       endif
    enddo !nkin
    !  follow track through plasma, compute length through each zone
  end subroutine nbi_halo_new_tr
  
  
  subroutine nbi_fastv_sample(ns,n_iz,ie,ia,ig)
    integer, intent(in) :: ns  ! fast ion species index
    integer, intent(in) :: n_iz  ! beam spatial zone index
    real(fp) :: dsum,zzn,de,da,d_alpha
    integer, intent(out) :: ie,ia,ig
    dsum=0.0_fp
    da = 1.0_fp/nznbma
    d_alpha = twopi/mgyro
    ig=1
    if(nlfbmflr) then
       do ie=1,nznbmea(ns)
          de=efbmb(ie+1,ns)-efbmb(ie,ns)
          do ia=1,nznbma
             do ig=1,mgyro
                dsum=dsum+fbm_ptcl(ie,ia,n_iz,ns,ig)*de*da*bmvol(n_iz)*d_alpha
             enddo
          enddo
       enddo
       if(dsum.eq.0.0_fp) then
          call errorlog(' nbi_fastv_sampl error in fi_ptcl sampling: no ptcls in beam zone=',n_iz)
       endif
       zzn=nbrandm_pll(idum)*dsum
       dsum=0.0_fp       
       do ie=1,nznbmea(ns)
          de=efbmb(ie+1,ns)-efbmb(ie,ns)
          do ia=1,nznbma
             do ig=1,mgyro
                dsum=dsum+fbm_ptcl(ie,ia,n_iz,ns,ig)*de*da*bmvol(n_iz)*d_alpha
                if(zzn.le.dsum) then
                   return
                endif
             enddo
          enddo
       enddo
    else
       do ie=1,nznbmea(ns)
          de=efbmb(ie+1,ns)-efbmb(ie,ns)
          do ia=1,nznbma
             dsum=dsum+fbm(ie,ia,n_iz,ns)*de*da*bmvol(n_iz)
          enddo
       enddo
       if(dsum.eq.0.0_fp) then
          call errorlog(' nbi_fastv_sampl error in fast ion sampling: no ptcls beam zone=',n_iz)
       endif
       zzn=nbrandm_pll(idum)*dsum
       dsum=0.0_fp       
       do ie=1,nznbmea(ns)
          de=efbmb(ie+1,ns)-efbmb(ie,ns)
          do ia=1,nznbma
             dsum=dsum+fbm(ie,ia,n_iz,ns)*de*da*bmvol(n_iz)
             if(zzn.le.dsum) then
                return
             endif
          enddo
       enddo
    endif
    
  end subroutine nbi_fastv_sample

  subroutine dep_neutrals(ier)
    integer :: ier,i,j,iz,idotrk
    real(fp) :: r_u(3),r_dep(3)
    real(fp) :: zlava,zlp,zlav,zlvac,zlvacp
    real(fp) zprb(8*mj), zlng(8*mj), zedep(8*mj)
    real(fp) zr,zphi
    real(fp) zmfpi,zthnsg,zbbnsg,znsigl,zdnsl,zexpol,zexp
    integer :: ig !thermal plasma species id, if ig=0 -- then this is not a thermal species
    integer :: irng ! 1st guess to range # of energy grid in nbeloc.
    integer :: isb  !fast particle   
    integer :: irange_bb ! table energy range sector
    integer :: jba_0,lgyro_neut !pitch angle  and gyro  indices
    real(fp) :: zbbnsv !total n*sig_v from interaction with a fast ions
    integer :: ind,indp1 !energy grid index
    real(fp) :: zintri,zintrp
    real(fp) :: zeovab ! energy / nucleon of neutral in plasma frame, [ev/amu]
    real(fp) :: zbnsvis,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zexcfac,zbnsvcx(mig) !n*sig*v of reactions neutrals(fast,thermal) + thermal ion
    real(fp) :: zbbnsva(2,mibs) !n*sig*v of reactions neutrals(fast,thermal) + fast ion
    real(fp) :: zbmsexc1d, zsigvnt1d !n*sig*v cx 
    real(fp) :: r_wrk(3),r_end_wrk(3),r_wrk1(3)
    !take a cleanup

    zbnsvis=0._fp
    zbnsvii=0._fp
    zbnsvie=0._fp
    zbnsviz=0._fp
    zbnsvto=0._fp
    zexcfac=1._fp
    zbnsvcx(:)=0._fp
    zbbnsva(:,:)=0._fp
    zbmsexc1d=0._fp
    zsigvnt1d=0._fp
    r_wrk=r
    r_end_wrk=r_end
    
    ier=0
    ! r, r_end start and end points in tokamak cartesian coord. [m,m,m]
#ifdef __MPI
        htime15=mpi_wtime()
#endif
  
    call nbi_track_call_halo(r,r_end,ier) ! find intersection with the irregular 2d mc grid found

#ifdef __MPI
        htime16=mpi_wtime()
        htime_tr=htime_tr+htime16-htime15
#endif
    if(ier.ne.0) then
       call errorlog(' dep_neutrals error: no intersection with the irregular 2d mc grid found ')
       call bad_exit
    endif
    if(nzcrtk.gt.8*mj) then
       call errorlog('dep_neutrals: nzcrtk > limits')
       call errorlog('dep_neutrals: nzcrtk = ',nzcrtk)
       call errorlog('dep_neutrals: limit = ',8*mj)
       call errorlog('probably due to code error in xstraln or its srs') 
       if(myidd.eq.0) then
          write(nonlin,*) ' ??dep_neutrals: nzcrtk = ',nzcrtk,' > limit ',8*mj
          write(nonlin,*)'    (probably due to code error in xstraln or its srs'
       endif
       ier=1
       call bad_exit
    endif

    r_u(:)=r_end(:)-r(:) 
    r_u(:)=r_u(:)/dsqrt(sum(r_u(:)*r_u(:))) !unit vector
       
    vion=itrackh(idtr)%vneut*100._fp !velocity of neutral, [cm/s]
    idotrk=0
    if(levmod_halo.eq.2) then
       idotrk=idtr+2
    endif
    ig=0 !means that this is fast particle
    if(itrackh(idtr)%type.eq.'THERM') then
       ig=lsbeam-nsbeam !back plasma id
    endif

    znsigl=0.0_fp
    zlvac =0.0_fp  
    zlvacp=0.0_fp
#ifdef __MPI
  htime3=mpi_wtime()
#endif
    call addn0(0.0_fp,1.0_fp,0.0_fp,vion,idotrk,lsbeam_parent) !initializations only
#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
!
!  check if track misses plasma entirely
!
      if(nzcrtk.eq.0) then
         itrackh(idtr)%maxpts=0
         itrackh(idtr)%idgen=-1
! 
         return
      endif
!  vacuum segment length
      zexp=1.0_fp
      zlvac=xltkbf
      if(zlvac.gt.0.0_fp) then
#ifdef __MPI
  htime3=mpi_wtime()
#endif
         call addn0(zlvac,zexp,0.0_fp,vion,idotrk,lsbeam_parent)
#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
      endif
!  track through plasma, compute deposition profiles by plasma zone,
!  compute probability of deposition in each zone.    

!  loop thru track segments      
      zlava=xltkbf
      zlng=0.0_fp
      do i =1, nzcrtk
!  transp zone index
         j=nznttk(i)
!  beam zone index
         iz=nznbtk(i)


         if(i.eq.1) then
            zlp=0.0_fp
            zlng(1)=xlngtk(1)
         else
            zlp=zlng(i-1)
            zlng(i)=zlp+xlngtk(i)
         endif
!
!  distrance tracked thru plasma - to zone midpt.
         zlav=0.5_fp*(zlng(i)+zlp)

         if(j.eq.lep1) then
#ifdef __MPI
  htime3=mpi_wtime()
#endif
            call addn0(zlvac+zlng(i),zexp,0.0_fp,vion,idotrk,lsbeam_parent)
#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
            cycle
         elseif(j.gt.lep1.or.iz.le.0) then
#ifdef __MPI
  htime3=mpi_wtime()
#endif
            call addn0(zlvac+zlng(i),zexp,0.0_fp,vion,idotrk,lsbeam_parent)
#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
            cycle
         endif
!
!  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
            if(zlvac.gt.zlvacp) then
#ifdef __MPI
  htime3=mpi_wtime()
#endif
               call addn0(zlvac+zlng(i-1),zexp,0.0_fp,vion,idotrk,lsbeam_parent)

#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
            endif
         endif
         zlava=zlav+zlvac

!mg why         r_dep(:)= -(r(:)+r_u(:)*zlava*0.01_fp) !deposition point in tokamak cartesian coord,[m,m,m]
         r_dep(:)= r(:)+r_u(:)*zlava*0.01_fp !deposition point in tokamak cartesian coord,[m,m,m]
         v_xyz = abs(vion)*r_u !v_x, v_y, v_z in tokamak cartesian coord system [cm/sec]
         !convert tokamak cartesian coord system to tokamak cyl. coord system
#ifdef __MPI
        htime11=mpi_wtime()
#endif
         call eq_rcyl(nsnccw*nsjdotb,r_dep(1),r_dep(2),zr,zphi) ![m,m] -> [m,rad]
#ifdef __MPI
        htime12=mpi_wtime()
        htime_gm=htime_gm+htime12-htime11
#endif
         ! zeovab is energy of neutral in plasma frame

          call nbe_rel(zeovab,v_xyz,zr*100._fp,zphi,j)! plasma rotation is taken into account

         irng=1
! now get atomic physics data    
         call get_lsbeam(idtr,lsbeam)
         !thermal or fast neutrals interact with thermal plasma
#ifdef __MPI
          htime9=mpi_wtime()
#endif
         call nbeloc(zeovab,irng,ind,zintri,indp1,zintrp) ! lookup interaction energy index in nb atomic physics tabl,
                                                          ! lsbeam in nbi_com already
         call nbains(zbnsvcx,zbnsvii,zbnsvie,zbnsviz,zbnsvto,zexcfac,&
         &           zbmsexc1d, zsigvnt1d,j,ind,zintri,indp1,zintrp)
#ifdef __MPI
          htime10=mpi_wtime()
          htime_at=htime_at+htime10-htime9
#endif
         if(itrackh(idtr)%type.eq.'THERM') then
            !tehrmal neutral with fast ion interaction
            ig=lsbeam-nsbeam
            do isb=1,nsbeam
               !thermal neutrals interact with fast ion
               zbbnsva(1,isb)=snb0i2(iz,ig,isb)
               zbbnsva(2,isb)=snbn02(iz,ig,isb)
            enddo
         else
          ! fast neutrals(gen>"0") are interacted with fast ions
            irange_bb=2 !for beam fast ptcl.
          ! in nb atomic physics tabl
          ! fast neutral (gen.>0) interact with a fast ion
          ! need to know energy of neutral in lab frame 
          ! (see cxnsum.for)
#ifdef __MPI
          htime9=mpi_wtime()
#endif
            call cxnsum_intrp(vion*vion/evptocm_sec2,-1,irange_bb,ind,indp1,zintrp) ! lookup interaction energy index 
#ifdef __MPI
          htime10=mpi_wtime()
          htime_at=htime_at+htime10-htime9
#endif
            jba_0=itrackh(idtr)%ipitch
            lgyro_neut=itrackh(idtr)%igyro
            if(itrackh(idtr)%zspecies.eq.1)then !h-isotopes
               do isb=1,nsbeam
                  
                  zbbnsva(1,isb)=(1.0_fp-zintrp)*svn3vii_h(ind,jba_0,lgyro_neut,iz,isb)+zintrp*&
                       & svn3vii_h(indp1,jba_0,lgyro_neut,iz,isb)
                  
                  zbbnsva(2,isb)=(1.0_fp-zintrp)*svn3vcx_h(ind,jba_0,lgyro_neut,iz,isb)+zintrp*&
                       & svn3vcx_h(indp1,jba_0,lgyro_neut,iz,isb)
                  
               enddo
            elseif(itrackh(idtr)%zspecies.eq.2)then !he-isotopes
               
               do isb=1,nsbeam
                  
                  zbbnsva(1,isb)=(1.0_fp-zintrp)*svn3vii_he(ind,jba_0,lgyro_neut,iz,isb)+zintrp*&
                       & svn3vii_he(indp1,jba_0,lgyro_neut,iz,isb)
                  
                  zbbnsva(2,isb)=(1.0_fp-zintrp)*svn3vcx_he(ind,jba_0,lgyro_neut,iz,isb)+zintrp*&
                       & svn3vcx_he(indp1,jba_0,lgyro_neut,iz,isb)
                  
               enddo
            endif
         endif
         !
         !  total ionization fraction (after adjustment)
         !  and thermal plasma stopping n*sigma
         !
         zbnsvis=zbnsvii+zbnsvie+zbnsviz
         zthnsg= zbnsvto / vion
         zbbnsv=0.0_fp
         do isb=1,nsbeam
            zbbnsv=zbbnsv+zbbnsva(1,isb)+zbbnsva(2,isb)
         enddo
         if(zbbnsv.gt.0.0_fp) then
            zbbnsg=zbbnsv/vion
         else
            zbbnsg=0.0_fp
         endif
         !
         !  inverse mean free path:
         !
         zmfpi=zthnsg+zbbnsg
         zdnsl=xlngtk(i)*zmfpi
         zexpol=zexp
         !  n*sigma*dl integral
         znsigl=znsigl+zdnsl
         zexp=exp(-znsigl)
         !  ok have survival probability -- increment neutral density box calc.
#ifdef __MPI
  htime3=mpi_wtime()
#endif
         call addn0(zlvac+zlng(i),zexpol,zmfpi,vion,idotrk,lsbeam_parent)
#ifdef __MPI
  htime4=mpi_wtime()
  htime_add=htime_add+htime4-htime3
#endif
!mg  if(lsbeam_parent.le.nsbeam.and.lsbeam.gt.nsbeam)then
!mg   sn0s3dh(j,ig,lsbeam_parent)=sn0s3dh(j,ig,lsbeam_parent) + (zexpol-zexp)*itrackh(idtr)%weight 
!mg  endif
          
      enddo ! end of  loop thru track segments, i

  end subroutine dep_neutrals

  subroutine nbe_rel(ze_rel,v_xyz,z_r,z_phi,nj)
    !output argument ze_rel = energy / nucleon of neutral in plasma frame (p.f), [ev/amu]
    ! in lab frame (l.f)
    !              v_x(l.f)=v_x(p.f) - v_rot*sin(zphi)
    !              v_y(l.f)=v_y(p.f) + v_rot*cos(zphi)
    ! go to (p.f)
    ! e(p.f)=e(l.f) + v_rot^2 +2*v_rot*v_x(l.f)*sin(zphi)-2*v_rot*v_y(l.f)*cos(zphi)


    real(fp),intent(out) :: ze_rel   ![ev/amu]
    real(fp),intent(in) :: z_r,z_phi  !(r,phi) in tokamak cylindrical coor. system, [cm,rad]
    real(fp),intent(in) :: v_xyz(3) !v_x, v_y, v_z in tokamak cartesian coord system [cm/sec]
    integer,intent(in) :: nj !transp zone index
    real(fp) :: v_rot

    v_rot= 0.5_fp*(omegag(nj)+omegag(nj+1))*z_r !plasma rotation,[cm/sec]
    ze_rel=(sum(v_xyz*v_xyz)+2._fp*v_rot*(v_xyz(1)*sin(z_phi)-v_xyz(2)*cos(z_phi))&
         +v_rot**2)/evptocm_sec2 !v2toev, plasma rotation 

  end subroutine nbe_rel

  subroutine get_lsbeam(id_tr,ilsbeam)
    integer,intent(in) :: id_tr !neutral track id #
    integer,intent(out) :: ilsbeam !lsbeam for this neutral
    integer :: ig,isb
    
    ilsbeam=0
    if(itrackh(id_tr)%type.eq.'THERM') then
       do ig=1,ng
          ilsbeam=nsbeam+ig
          ! find thermal species id
          if(itrackh(id_tr)%zspecies.eq.backz(ig).and.itrackh(id_tr)%aspecies.eq.aplasm(ig))exit
       enddo
    else
       do isb=1,nsbeam
          if((itrackh(id_tr)%zspecies.eq.xzbeams(isb)).and.(itrackh(id_tr)%aspecies&
               &.eq.abeams(isb)).and.(itrackh(id_tr)%type.eq.'FASTB').eqv.(.not.nlfprod(isb))) then
             exit
          endif
       enddo
       ilsbeam=isb
    endif
    
  end subroutine get_lsbeam

  subroutine param_fbm_ptcl
    implicit none
    real(fp) :: zbmnrm,zdelebd,zdnw,zbmn
    integer :: inz,isb,ine,ina
    if(allocated(bdns2p_ptcl)) deallocate(bdns2p_ptcl)
    allocate(bdns2p_ptcl(nfbznsi,nsbeam))
    bdns2p_ptcl=0.0_fp
    zbmnrm=real(nznbma)/(delt)
    if(nlfbmflr)then
       do inz = 1, nfbznsi
          zbmn=zbmnrm/bmvol(inz)
          do isb = 1, nsbeam
             do  ine = 1, nznbmea(isb)
                zdelebd=(efbmb(ine+1,isb)-efbmb(ine,isb))
                do ina = 1, nznbma
                   zdnw=sum(fbm_ptcl(ine,ina,inz,isb,1:mgyro))
                   if(zdnw.gt.0.0_fp) then
                      bdns2p_ptcl(inz,isb)=bdns2p_ptcl(inz,isb)+zdnw*zdelebd/zbmn
                   endif
                enddo
             enddo
          enddo
       enddo
    else
       bdns2p_ptcl(1:nfbznsi,1:nsbeam)=bdns2p(1:nfbznsi,1:nsbeam) ! at gc
    endif
    return
  end subroutine param_fbm_ptcl

  !
  ! ------------------ pdelt -------------------
  ! returns the average amount of time a particle
  ! exists in a time interval delt given a decay rate of kappa
  !
  real(fp) function p_delt(delt,kappa,decay)
    real(fp), intent(in)  :: delt ! advancement time
    real(fp), intent(in)  :: kappa ! decay factor
    real(fp), intent(out) :: decay ! exp(-delt*kappa), ugly
    real(fp) :: x,z
    
    x = kappa*delt
    if (abs(x)<1.d-4) then
       z = (1.0d0-0.5d0*x*(1.-x/3.0d0))
       p_delt = delt*z
       decay = 1.-x*z
    else if (x>40) then
       p_delt = 1./kappa
       decay = 0.0d0
    else
       decay = exp(-x)
       p_delt = (1.-decay)/kappa
    end if
  end function p_delt

  
end subroutine depbox_all
