!---------------------------------------------------------------------
!  depall -- beam code:  monte carlo decisions and deposition loop
!   gen for multispecies beam and broken out of nbmain, dmc 8/88
!    **nbmain is now fimain**  -- dmc 27 jan 1997
!
!   this routine now supports deposition for fppmod as well as
!   deposition for the mc code.
!
subroutine depall(ierr)
  !
  !  jec@ccfe.ac.uk
  !                 06jan2010  -- fix use of uninitialised iokill
  !
  !  mg             march 2009 -- for mpi each process carries a single
  !                 sublist of particles  instead of both sublist and
  !                 whole list arrays. master processors (myidd=0)
  !                 keeps global control for deposited ptcls.
  !
  !  dmc            29 nov 2005 -- cleanup for parallelization
  !                                ptcl loops consolidated...
  !
  !  jfs            23 aug 1994 -- added init. of time-step magn. factor
  !  mhr              july 1994 -- ripple
  !  jfs            17 feb 1994 -- added fusion product protons
  !  tbt mods          jun 1993 -- c$   ! tbt
  !  dmc modifications feb 1992 -- adding a deposition particle balance
  !   check
  !
  use nbi_com
  use nbi_mpi_com
  use nbi_random
  use nbi_random_dep
  use tube_mod, only: deptube
  use dep_mpi_split
  use dep_track_counter
  use deptrk_io_track
  use logmod
  use nbi_bscom
  use nbdep_debug, only : icount_depall

  !
  !============
  implicit none
  integer :: iabs
  logical, save :: first_entry = .true.
  !============
  integer, intent(out) :: ierr
  !
  real(fp) :: ztest
  !
  real(fp) :: zwgrat,zwnorm,zwgfac
  real(fp) :: zbnew(mibs),zfnew(mibs),zfnewa(mibs),zeinjav(mibs),zbold
  real(fp) :: zwrats,zwght3,zcorr,zwrats_tr
  real(fp) :: zrrwfac(mibs)
  integer :: isb,iptcls(mibs),iptcl0(mibs),inlight(mibs),inuse(mibs)
  integer :: iplim,inold,innew,j,inlim,idep00,jdep00
  integer :: inloc,jxx,ind,innp1,ind1,inkill(mibs),iokill
  integer :: idep0_new,inuse_prev(mibs)
  integer :: inewmax,inusmax,inusmin,idenom,ipmin
  logical :: iorbrr(mibs)
  !
  !%      integer :: ibplot,ipln1,ipln2,iplnum,iplsta
  !%      integer :: inew,inu,ict,iplcx,iplfp,isty,ipflr
  !%      real(fp) :: ztim
  !
  real(fp) :: ztim1,ztim2
  real(fp) :: frac_dep_tausl !fraction of new ptcls. to be deposited

  integer :: iscreen1
#ifdef __MPI
  integer :: iread,iread1      ! for mpi debugging
#endif
  integer :: insbeam,iseed,idum, minb_new,minb_prev,minb_tot
  integer :: reckon_minb       ! integer function to calculate minb
  integer :: find_ptcl_ind     ! integer function to find ptcl.id#
  integer, dimension(:), allocatable :: nuse_max
  integer :: nusemax = 0
  integer :: inuse_max,isb0,max_kill
  integer :: mpierr
  !
  !  for debugging...
  !  serial mode
  !
  integer,dimension(:),allocatable :: ninddumn_tr
  !
  logical, save :: init_depall = .false.
  logical :: key_myidd_0
  integer :: nusey, nusey_tr
  integer, dimension(:), allocatable ::nuse_total ! total # of mc ptcls currently in storage
  ! is using for a global control during deposition
  ! and orbiting
  integer :: iptcl,jloc
  character(len=150) :: fname_closed !netcdf file name for track data
  logical :: iltrk0     ! true to save data on beam neutral track
  logical :: il_birth ! if .true. gather fast ion birth source
  real(fp) :: tdepay_ind
  real(fp) :: dumj(2,2)
  integer :: innew_min
  ! for modelling bi-maxwellian distribusion
  real(fp) :: zs_pll,eb_av,ewghtdep_pr,wghtdep_pr,ze,zdt
  real(fp) :: vperp_1,vperp_2

  !
  if(allocated(nuse_total)) deallocate(nuse_total)
  allocate(nuse_total(mibs)); nuse_total = puseyp + 0.5_fp
  !
  !-----------------------------------------------------------------
  !
  call infolog('entering depall: start nbsteps=',nbsteps)
  if(.not.init_depall) then
     icount_depall = 0
     init_depall = .true.
  end if
  !
#ifdef __MPI
  ttime21=0.0_fp
  ttime22=0.0_fp
  ttime21=mpi_wtime()
#endif
  insbeam=nsbeam
  ierr=0
  kdep0=0
  kdep0sum=0
  idep00=0
  !

#ifdef __MPI
  ! for mpi debugging purpose
  iread=0
  iread1=0
  if ((iread1.gt.0).and.(numprocd.eq.1)) then
     call nubeam_buf_read
     call nubeam_cxbuf_read
  end if
#endif

  if(only_io) then
     call errorlog( &
          ' emergency exit depall: improper call, only_io = .true.')
     ierr=666
     return
  end if
  !
  call cpu_time(ztim1)
  !
  minb_prev=minb
  minb_new = reckon_minb(numprocd) !find # of particles (if mpi -  per each processor)
#ifndef __MPI
  !     this will be a path for nubeam_test, when debugging mpi nubeam for one ptcl.sublist
  if(debug_mpi.and.first_entry) then   ! when debug_mpi=.true.
     ! keep only ptcl.sublist
     ! set only once at first enrty.
     !
     minb_new = minb                   ! minb,nuse have been read from
     ! <runid>_ptcl_state_cpu<myid>_<numprocd>.dat
     ! for given <myid> and <numprocd>
     nuse_total = nuse

  end if
#endif
  if(minb_new.gt.minb) then          ! check if
     ! ptcl. list arrays need to be expanded
     call infolog('depall: reallocated ptcl.list')
     call nbi_realloc_minb(minb_new) ! allocate array with new
     ! minb~max(nptcls,nptclf)/numprocd
     ! see nbi_alloc_minb
     ! minb redefined there and =minb_new
     call infolog('depall: ptcl.list dim. extended from minb =', &
          minb_prev)
     call infolog('                                to minb=',minb)
     call depall_mpi_split(minb_prev*numprocd+1, & !splitting new vacancies &
          minb*numprocd,nusey)                    !between processors
     call infolog('depall: splitting done')
     do isb=1,insbeam
        nbndex(minb_prev+1:minb,isb)=ninddumn(1:nusey)  !each processors storage indexing array 
        ntotptp(isb)=ntotptp(isb)+nusey  !each processors total # of m.c ptcls.
     end do
  end if

  !
  !      set up trkrb zones:  no pitch angle zones, no zones outside the
  !       plasma - cf subroutines shrink, rbloat (dep.for)
  !
  !  "nbi_shrink" sets nznbma=1 to speed straight line tracker, but,
  !  dep subroutine fdep needs to know the unshrunk value.  so:
  !
  nznbma_save=nznbma
  !
  call nbi_shrink
  !
  !%        ibplot=nbplot
  !
  !------------------------------------------
  !  numerical controls
  !
  frac_depmax=max(0.1_fp,min(1.0001_fp,frac_depmax))
  frac_depmin=max(0.0_fp,min(0.1_fp,frac_depmin))
  !
  !------------------------------------------
  !  adjust wbav; fetch current statistics...
  !

  do isb=1,insbeam

     lsbeam=isb

     !  save previous nuse
     inuse_prev(isb) = nuse_total(isb)
     !
     !  number of underweight particles
     inlight(isb) = wnlight(isb)+0.5_fp
     !
     !  count 1/2 ptcl each for underweight particles, when
     !  deciding how many to kill to make room for new ions coming in...
     !mg         inuse(isb) = nuse_total(isb)-inlight(isb)/2

     inuse(isb) = nuse_total(isb)
     !
     if(nmctot(isb).eq.0) cycle
     !
     !  correct wbav_kill according to summation from prior orball call
     !
     zwnorm=0
     do j=lcentr,ledge
        zwnorm=zwnorm+wnmc(j,isb)/wbav_kill(j,isb)
     end do
     zwnorm=zwnorm/nmctot(isb)
     zwnorm=max(0.8_fp,min(1.25_fp,zwnorm))
     wbav_kill(1:mj,isb)=zwnorm*wbav_kill(1:mj,isb)

  end do
  !
  !------------------------------------------
  !
  !  inverse deposition fraction, by beam & energy (prev. timestep data)...
  !  (1:3,1:mib)
  !
  frac_sur_inv = 1.0_fp/frac_survive
  zidep0=0.0_fp
  zitrk0=0.0_fp
  !
  !  beam species loop
  !
  !  serial loop -- mc census control...
  !


  !     do isb=1,nsbeam
  do isb=1,insbeam
     !
     lsbeam=isb             ! these may be used interchangeably
     call inispc
     !calculate average tau sl.down
     lsbeam=isb             ! these may be used interchangeably
     if(tausl_avg(isb).le.0.0_fp) then
        !first call: estimate approximate slowing down time
        call nbi_tausl
     end if
     !
     !xx         call nbi_zcheck_dep            ! debugging...
     !
     zeinjav(isb) = einjav(isb) ! save old value
     einjav(isb) = 0        ! initialize sum for <e>
     weinjav(isb) = 0       ! weighting for einjav sum
     xmcloss_d(isb) = 0     ! prompt loss mc ptcl counts  -- depall
     xmcloss_o(isb) = 0     ! prompt loss mc ptcl counts  -- orball
     !
     !  for heavy neutral deposition:  neutrals are assumed only singly
     !  ionized when deposited
     !
     !cc   (already set in inispc)      if(nlheavy(isb)) xzbeami=1.0_fp
     !
     !  census control depends on whether this is a beam species or a
     !  fusion product species
     !
     if(nlfbmfpp(isb)) then
        iptcls(isb)=1       ! just flag:  dep for fpp
        iptcl0(isb)=1
     else if(nlfprod(isb)) then
        !
        !     dmc 5 oct 1992 -- tau(sl) and tau(pa) calculation moved up here
        !
        call fustaus
        !
        iptcl0(isb)=nptclf
        if(isb.eq.nftrit) then
           iptcls(isb)=nptclf*min(1.0_fp,(pinjs(isb)/plfst))
        else if(isb.eq.nfprot) then
           iptcls(isb)=nptclf*min(1.0_fp,(pinjs(isb)/plfsp))
        else if(isb.eq.nfprot14) then
           iptcls(isb)=nptclf*min(1.0_fp,(pinjs(isb)/plfsp14))
        else if(isb.eq.nfhel3) then
           iptcls(isb)=nptclf*min(1.0_fp,(pinjs(isb)/plfhe3))
        else if(isb.eq.nfhel4) then
           iptcls(isb)=nptclf*min(1.0_fp,(pinjs(isb)/plfhe4))
        end if
        !  ok but don't let the census fall if the power drops transiently:
        iplim=min(nptclf,inuse(isb))
        iptcls(isb)=max(iptcls(isb),iplim)
        !  if #in use is currently zero, do not allow iptcls(isb) > 0 unless
        !  a threshold of a few ptcls is reached...
        if(nuse_total(isb).eq.0) then
           ipmin = 0.5_fp + frac_depmin*nptclf
           if(iptcls(isb).lt.ipmin) iptcls(isb)=0
        end if
     else
        iptcls(isb)=nptcls
        if(nlheavy(isb).and.(nptclh.gt.0)) iptcls(isb)=nptclh
        iptcl0(isb)=iptcls(isb)
     end if
     !
     !  if beam has been turned off but there are still
     !     beam ptcls flying around, follow these ptcls, but deposit
     !    no new ptcls.  when all beam ptcls have thermalized or
     !  left plasma, nuse=0
     !
     if((iptcls(isb).eq.0).or.(xninjs(isb).le.0.0_fp)) then
        ennew(isb)=0
        nnew(isb)=0
        inkill(isb)=0
        iorbrr(isb)=.false.
        cycle
     end if
     !
     !  deposition...
     !  monte carlo decisions *************
     !
     inusmin=100
     !do not apply rr
     !mg         ztest=inlight(isb)*1.0_fp/iptcl0(isb)
     !mg         if((ztest.ge.frac_orbrr).and.(inuse(isb).ge.inusmin)) then
     !mg            call infolog('depall: species ',isb)
     !mg            call infolog('underweight particles = ',inlight(isb))
     !mg            call infolog('russian roulette applied.')
     !mg            if(myidd.eq.0) then
     !mg               write(nonlin,*) ' %dep species ',isb,' -- ',inlight(isb),
     !mg     >              ' underweight particles:  russian roulette applied.'
     !mg           end if
     !mg            iorbrr(isb)=.true.
     !mg
     !mg         else
     !mg            call infolog('depall: species ',isb)
     !mg            call infolog('underweight particles = ',inlight(isb))
     !mg            if(myidd.eq.0) then
     !mg               write(nonlin,*) ' %dep species ',isb,' -- ',inlight(isb),
     !mg     >              ' underweight particles.'
     !mg            end if
     !mg            iorbrr(isb)=.false.
     !mg         end if

     iorbrr(isb)=.false.
     inkill(isb)=0


     ennew(isb)=pinjs(isb)*6.2418e18_fp*delt
     !
     zbold=bntots(isb)  ! weight stored before depo
     zbnew(isb)=xninjs(isb)*delt
     zfnew(isb)=zbnew(isb)/(zbnew(isb)+zbold)
     !
     !  ratio, axial weight to edge weight...
     !
     if(nlfprod(isb).or.nlfbmfpp(isb)) then
        xwght0(isb)=1.0_fp
     else
        !
        !  avg old & new values according to new injected ptcl fraction
        !
        if(wghta.gt.nzones) wghta = nzones
        wghta = max(1.0_fp,wghta)
        !
        !mg            xwght0(isb)=(1.0_fp-zfnew(isb))*xwght0(isb) +
        !mg     >           zfnew(isb)*min(1.0_fp,
        !mg     >           max(1.0_fp/wghta,xhwdep0(isb)))
        ! mg 01/17/18            xwght0(isb)=1.0_fp/wghta
        xwght0(isb)=(1.0_fp-zfnew(isb))*xwght0(isb) + &
             zfnew(isb)*min(1.0_fp, &
             max(1.0_fp/wghta,xhwdep0(isb)))
        !            xwght0(isb)=1.0_fp/wghta
        !
     end if

     !
     if(nlfbmfpp(isb)) then
        nnew(isb)=0
        nuse_total(isb)=0
        inkill(isb)=0
        iorbrr(isb)=.false.
        cycle
     end if
     !
     ! ** d. mc cune sept 29 1979 **
     !   testing new control for deposition:
     !   nnew (# of new ptcls to deposit) is related to the number
     !   of new particles coming in compared to the number already
     !   there, instead of the energy of the incoming particles
     !   compared to the energy already there.
     !
     !  bntot calculated in orbit loop
     !
     inold=0.5_fp+ (1.0_fp-zfnew(isb))*iptcls(isb) ! ideal old ptcl count
     innew=0.5_fp+ zfnew(isb)*iptcls(isb) ! ideal new ptcl count
     !
     !        russian roulette-- kill old particles to make room for new
     !
     !mg!
     frac_dep_tausl=delt/tausl_avg(isb)*frac_dep_lim
     frac_dep_tausl=min(1.0_fp,frac_dep_tausl)
     inewmax = iptcls(isb)*frac_dep_tausl +0.5_fp

     !         inewmax = iptcls(isb)*frac_depmax !disable

     inewmax = max(1,inewmax)
     if(xpand_nptcl.lt.1.0_fp.and.xpand_nptcl.gt.0.5_fp) then
        xpand_nptcl = 1.0_fp ! restrict 0.5 < xpand_nptcl < 1.0
        ! or may get trouble later when
        ! inloc bigger than ntotptp
        call infolog('depall: restrict 0.5 < xpand_nptcl < 1.0', &
             xpand_nptcl)
        if(myidd.eq.0) then
           write(nonlin,*) &
                ' %depall: restrict 0.5 < xpand_nptcl < 1.0', &
                xpand_nptcl
        end if
     end if
     !
     inusmax = iptcls(isb)*(1.0_fp+xpand_nptcl)
     !
     innew_min=1
#ifdef __MPI
     innew_min=numprocd ! allow each cpu to have at least 1 ptcl
#endif
     innew=max(innew_min,min(innew,inewmax))
     !
     !  no rr if curent number is too small
     !
     inusmax = max(inusmax,inusmin+innew)
     !
     !  russian roulette -- underweight ptcls?
     !
     !mg         if(iorbrr(isb)) then
     !mg            iokill=inlight(isb)/2
     !mg         else
     !mg            iokill=0
     !mg         end if
     !
     !  allow total number to float within a range (try to avoid russian roulette)
     !
     !
     ! now find largest value of nuse to estimate upper limit for inuse. to privent
     ! the  injection of the new partilces in amount more than {ntotptp-(nuse+nnewp)}
     ! on each processor without expanding of the particles list.
#ifdef __MPI
     if(allocated(nuse_max)) deallocate(nuse_max)
     allocate(nuse_max(numprocd))
     call mpi_allgather(nuse(isb),1,MPI_INTEGER,nuse_max,1, &
          MPI_INTEGER,nbi_comm,mpierr)
     nusemax=maxval(nuse_max,dim=1)
     deallocate(nuse_max)
#else
     nusemax=nuse(isb)
#endif
     inuse_max=nusemax*numprocd
     iokill   = 0                  ! jec 06jan10 avoid undefined variable ref in next if.. end if

     if(innew+inuse_max+minbs*numprocd.gt.minb*numprocd) then !more room needed to deposit
        ! new particles
        ! mg april 24,2014 i added 'minbs' as we allow
        ! so many recaptures in 'nutrav'
        if(innew+inuse_max+minbs*numprocd.gt.inusmax) then !apply russian roulette
           call infolog(' depall: apply rr')
           if(myidd.eq.0) then
              write(nonlin,*)' %depall: apply rr'
           end if
           iorbrr(isb)=.true.
           iokill=inlight(isb)/2
           nnew(isb)=innew
           inlim=inuse(isb)/2
           max_kill=inuse_max+innew-min(minb*numprocd,inusmax)
           inkill(isb)=min(inlim,max_kill)

        else
           minb_tot=inuse_max+innew+minbs*numprocd !total minb (on all processors) needed
           minb_tot=minb_tot+max(500,int(0.05*minb_tot))*numprocd !add some extra space to avoid trouble
           !in orball with recapture and splitting

           call split_minb(minb_tot,minb_new,numprocd)   !find minb  on each  processor (minb_new).
           !minb_tot could be increased to have factor
           !equals  'numprocd'
           minb_prev=minb
           call nbi_realloc_minb(minb_new) ! reallocate array with
           ! minb_new. minb redefined there and = minb_new
           call infolog('depall: avoid rr')
           call infolog('depall: ptcl.list was extended from minb=', &
                minb_prev)
           call infolog('                            to minb=',minb)

           if(myidd.eq.0) then
              write(nonlin,*)' %depall: avoid rr'

              write(nonlin,*) &
                   ' %depall: ptcl.list was extended from minb=', &
                   minb_prev, ' to minb=',minb

           end if
           call depall_mpi_split(minb_prev*numprocd+1, &!splitting new vacancies
                minb*numprocd,nusey) !between processors
           call infolog('depall: splitting done')

           do isb0=1,insbeam !update nbndex and ntotptp for all species
              nbndex(minb_prev+1:minb,isb0)= ninddumn(1:nusey) !each processors storage indexing array &
              ntotptp(isb0)=ntotptp(isb0)+nusey !each processors total # of m.c ptcls.
           end do

           inkill(isb)=0
           inlim=iptcls(isb)*frac_depmin
           inlim=max(2*innew,inlim)
           inlim=min(inewmax,inlim)
           nnew(isb)=max(innew,min(inlim,(iptcls(isb)-inuse(isb))))


        end if
     else
        ! mg 07/10/2013
        ! even we are tryin to avoid rr i added this part
        ! as some very light mc particles give a trouble
        ! by accelerating too much in some cases
        if(inlight(isb).gt.0) then
           iorbrr(isb)=.true.
        end if
        ! some (not all) light particles will be killed in nbkill_pll.for
        ! mg 07/10/2013
        inkill(isb)=0
        inlim=iptcls(isb)*frac_depmin
        inlim=max(2*innew,inlim)
        inlim=min(inewmax,inlim)
        nnew(isb)=max(innew,min(inlim,(iptcls(isb)-inuse(isb))))
     end if
     if(inuse(isb).gt.0) then
        idenom=inuse(isb)-inkill(isb)-iokill
        if(idenom.gt.0) then
           zrrwfac(isb)=inuse(isb)*1.0_fp/idenom
           zfnewa(isb)=nnew(isb)*1.0_fp/ &
                max(1,(nnew(isb)+inuse(isb)))
        else
           zrrwfac(isb)=1.0_fp
           zfnewa(isb)=1.0_fp
        end if
     else
        zrrwfac(isb)=1.0_fp
        zfnewa(isb)=1.0_fp
     end if
  end do

  !
  !  serial loop ends here
  !
  !---------------------------------------
  !  start parallelized code..........
  !  nbkill needs to operate on a sublist; global nuse(isb) is not known
  !    until end of pll section; nuse(isb)=nptcls only approximately not
  !    exactly at end of deposition!
  !
  !  r.r. if need more room
  !
  !     do isb=1,nsbeam
  stime4=0.0_fp
  stime5=0.0_fp
  stime6=0.0_fp
  stime7=0.0_fp
  stime8=0.0_fp
  stime41=0.0_fp
  stime42=0.0_fp
  stime51=0.0_fp
  stime52=0.0_fp
  stime61=0.0_fp
  stime62=0.0_fp
  stime71=0.0_fp
  stime72=0.0_fp
  stime81=0.0_fp
  stime82=0.0_fp
  stime10=0.0_fp
  stime101=0.0_fp
  stime102=0.0_fp
  !     check if fast ion birth source need to be gathered at this time step.
  !     this is necessary at restart
  il_birth=.false.
  if((nltrk_dep0.and.mibs_use.gt.0)) then
     il_birth=.true.
     call nbi_bsource_check(il_birth) !the presence of birth source files will be checked
  end if

  beamspeceis: do isb=1,insbeam

     !
     idep00 = 0
     lsbeam=isb             ! these may be used interchangeably

     iltrk0 = nltrk_dep0.and.(mibs2use(lsbeam).gt.0)

     call nbi_zcheck_dep

     call inispc

     !%        if(nlfprod(lsbeam)) then
     !%           write(nonlin,1001) 'fusion product',
     !%     >          abeams(lsbeam),xzbeams(lsbeam)
     !%        else
     !%           write(nonlin,1001) 'beam',
     !%     >          abeams(lsbeam),xzbeams(lsbeam)
     !%        end if
     !% 1001      format(' %depall: ',a,' species:  a= ',1pe11.4,' z= ',
     !%     >          1pe11.4)
     !%        iplnum=0
     !%        iplsta=1
     !%        call ureadi('depall:  no. of dep trajectories to plot:',
     !%     >    iplnum)
     !%        if(iplnum.gt.0) call ureadi(
     !%     >    'depall:  no. of first trajectory to plot:',iplsta)
     !------------------------
     !%        ipln1=iplsta
     !%        ipln2=iplsta+iplnum-1
     !%        if(iplnum.gt.0) then
     !%          inew=-1
     !%          call stplot(inew,inu,ict,iplcx,iplfp,isty,ipflr,
     !%     >           ztim)
     !%        end if

     !

     !     nuse_total can be redefined in nbkill_pll
     call nbkill_pll(iorbrr(isb),inkill(isb))

     !
     !  deposition-- deposit actual number needed to fill arrays
     !
     if(nnew(isb).le.0) then
        ennew(isb)=0
        nnew(isb)=0
     end if

     !----------------------------------------------------------------
     !   write particle list before orbit loop (mpi related)
     !-----------------------------------------------------------------



#ifdef __MPI
     stime61=mpi_wtime()
#endif

     !     if nuse_total(p) is zero, and nnew non-zero, we need  to split
     !     new ptcl.  between all  processors


     call depall_mpi_split(1,nnew(lsbeam),nusey)

     nnewp(lsbeam)=nusey    !# of new ptcl. on each processors

     if (nusey>0) then
        if(allocated(ninddumn_tr)) deallocate(ninddumn_tr)
        allocate( ninddumn_tr(nusey))
        ninddumn_tr=ninddumn
        if(iltrk0.and.il_birth) then
           !allocate arrays for birth source
           call nbi_bsource_alloc(nnewp(lsbeam),mibs_use)
           nbs_dep(mibs2use(lsbeam))=nnewp(lsbeam)
        end if
        !mg
        !mg            cycle
     end if


     !     if nuse(isb) is likely to be nonzero for the first time,
     !     additional initialization needs to be done
     if ((inuse_prev(isb).eq.0).and.(nnew(isb).ne.0)) then
        !     split global total # of ptcls. between all
        call depall_mpi_split(1,minb*numprocd,nusey)

        call infolog('depall: depall_mpi_split initinal done')
        if(myidd.eq.0) then
           write(nonlin,*)' %depall depall_mpi_split initinal done'
           write(nonlin,*)' %depall nuse(isb)= ',nuse(isb)
        end if
        ntotptp(isb)=nusey
        nbndex(1:ntotptp(isb),isb)=ninddumn(1:nusey)
     end if

#ifdef __MPI
     stime62=mpi_wtime()
     stime6=stime6+(stime62-stime61)
#endif


     ! intialize seed for parallel random no generator

     iseed=abs(10000000*nbrandm(idum))
     iseed=iseed*lsbeam
     iscreen1=1600+myidd

     jxx=0

     do j=1,nnewp(lsbeam)
        inloc=nuse(isb)+j
        jxx=j
        if (inloc.gt.ntotptp(lsbeam)) then
           call warnlog('depall: inloc bigger than ntotptp')
           call warnlog('inloc=',inloc)
           call warnlog('ntotptp=',ntotptp(lsbeam))
           call warnlog('nuse=',nuse(lsbeam))
           call warnlog('nnewp=',nnewp(lsbeam))
           if(myidd.eq.0) then
              write(nonlin,'(a2)') '  '
              write(nonlin,*)' !depall: inloc bigger than ntotptp:'
              write(nonlin,*) ' myidd=',myidd,' inloc=',inloc, &
                   ' ntotptp=',ntotptp(lsbeam),' nuse= ', &
                   nuse(lsbeam), '   nnewp= ',nnewp(lsbeam)
              !      nuse(lsbeam) + nnewp(lsbeam) must be .le. ntotptp(lsbeam)
              write(nonlin,'(a2)') '  '
           end if
        end if
        !%          if((ipln1.le.jxx).and.(jxx.le.ipln2)) then
        !%            nbplot=1
        !%          else
        !%            nbplot=0
        !%          end if
        !

        ind=ninddumn_tr(j)  !track id#

        call nbi_random_init_dep(iseed,ind,iscreen1) !provides unique rng seeds for tracks
        if(nlfprod(isb)) then
           icount_depall = icount_depall + 1
           idep00 = idep00 + 1
           call depfus(ierr)
           if(ierr.ne.0) then
              call errorlog('depall: depall>depfus (1)')
              call errorlog('j = ',j)
              call errorlog('icount_depall = ',icount_depall)
              if(myidd.eq.0) then
                 write(nonlin,*) &
                      ' ?depall>depfus (1) j,ierr,# = ',j,ierr, &
                      icount_depall
              end if
           end if
           dep_code_error = dep_code_error + ierr
        else
           icount_depall = icount_depall + 1
           idep00 = idep00 + 1
           if(abs(ntest_max).eq.1) then
              ewghtdep_pr=ewghtdep(isb)
              wghtdep_pr=wghtdep(isb)
           end if
           call nbidep(1,ierr)
           if(ierr.ne.0) then
              call errorlog('depall: depall>dep (1)')
              call errorlog('j = ',j)
              call errorlog('icount_depall = ',icount_depall)
              if(ierr.ne.0) then
                 if(myidd.eq.0) then
                    write(nonlin,*) &
                         ' ?depall>dep (1) j,ierr,# = ',j,ierr, &
                         icount_depall
                 end if
              end if
              dep_code_error = dep_code_error + ierr
           end if
        end if
        !
        !  stash ptcl
        !

        ind=nbndex(inloc,lsbeam) !ptcl id#
        ind=find_ptcl_ind(numprocd,myidd,ind) !find proper address for a given ptcl.id#
        nbscay(ind,isb)=nbsce
        nbienay(ind,isb)=nbien
        if(nbien.eq.1) then
           weinjav(isb)=weinjav(isb)+wght
           einjav(isb)=einjav(isb)+wght*einj0
        end if
        einjay(ind,isb)=einj0
        if(xjav.ge.lep1) then
           njavay(ind,isb)=lep1
        else
           njavay(ind,isb)=int(xjav)
        end if
        xiay(ind,isb)=xion
        thay(ind,isb)=th
        if(abs(ntest_max).eq.1) then
           ! by mg 08/08/2016 i am sampling mc particles for a  deposition  using random
           ! number generator for maxwellian distribution function
           ! ~exp(-(x-za)**2/(2.*zs**2))
           ! za = 0
           ! zs =  sqrt( c * t(j)/m_d ), where
           ! t(j)=2/3 ebavps(j) (beam average energy) in each transp j-zone.
           ! to be precise i am using ebavps(j) = max(ebavps(j), 2/3ti(j)).
           ! for bi-maxwellian i am sampling mc particles with  e_pll  and e_perp
           ! using random number generator with c_pll*t(j)  with c_perp*t(j), respectively.
           ! i compared profiles (nubeam_steps.for) for an input ebavps
           ! (energy i am sampling from) and output  (after the deposition).
           ! they look very close if good mc statistics was used.
           wght=1.0_fp
           eb_av=max(ebavps(nflr,isb), &
                0.01_fp*max(ti(nflr,2),te(nflr,2)))
           zs_pll=sqrt(.6666667_fp*eb_av/v2toev/2.0_fp) ! e=v**2 *v2toev
           ! 1/2 factor added as gaussian formula has factor 2
           ! ~exp(-(x-za)**2/(2.*zs**2))
           vperp=0.0_fp
           vpll=0.0_fp
           vperp_1=nbgauss_dep(sqrt(abs(bi_max_pr))*zs_pll,0.0_fp)

           !mg               vperp_2=nbgauss_dep(sqrt(abs(bi_max_pr))*zs_pll,0.0_fp)
           !mg               vperp=sqrt(vperp_1**2+vperp_2**2)
           vpll=nbgauss_dep(sqrt(abs(bi_max_pl))*zs_pll,0.0_fp)
           vperp=vperp_1 !sign is matter for a plots
           wght=wght*abs(vperp)
           !mg there is no sawthooth, i just borrow these names for
           ! bi-maxwellian test
           zdt=tbm2-(bitime(1)+tbm1)
           ebapl2sw(nflr,isb)=ebapl2sw(nflr,isb)+vpll*vpll* &
                v2toev*wght*zdt
           ebapp2sw(nflr,isb)=ebapp2sw(nflr,isb)+vperp*vperp* &
                v2toev*wght*zdt
           bdenss(nflr,isb)=bdenss(nflr,isb)+wght*zdt
           ze=vpll**2+vperp**2
           vion=sqrt(ze)
           ewghtdep(isb)=ewghtdep_pr+ze*v2toev*wght*zdt/delt
           wghtdep(isb)=wghtdep_pr+wght*zdt/delt
           if(vion.ne.0.0_fp) then
              xksid=vpll/vion
           else
              call errorlog &
                   (' depall: for abs(ntest_max)=1, vion = 0')
              ierr=1
              return
           end if
        end if
        vay(ind,isb)=vion
        xksidy(ind,isb)=xksid
        wghtay(ind,isb)=wght

        if (wght.le.0.0_fp) then
           call errorlog('depall: zero wght: ',wght)
           call errorlog('depall: lsbeam = ',isb)
           call errorlog('depall: ind = ',ind)
           if(myidd.eq.0) then
              write(nonlin,*) ' ?? depall: zero wght: ',wght
              write(nonlin,*) ' ?? depall: lsbeam = ',isb, &
                   ' ind = ',ind
           end if
        end if
        bitime(1)=max(epslon,bitime(1))
        !     tdepay(ind,isb)=bitime(1)      !changed by mg 11/2013 on ...
        tdepay(ind,isb)=bitime(1)+tbm1
        xistay(ind,isb)=xistate ! for radial current tracking...
        if(iltrk0.and.il_birth) then
           bs_r(j,mibs2use(lsbeam))=rmjdep
           bs_z(j,mibs2use(lsbeam))=ydep
           bs_rgc(j,mibs2use(lsbeam))=rmjion
           bs_zgc(j,mibs2use(lsbeam))=yion
           bs_xksid(j,mibs2use(lsbeam))=xksidy(ind,isb)
           bs_einj(j,mibs2use(lsbeam))=einjay(ind,isb)
           bs_wght(j,mibs2use(lsbeam))=wghtay(ind,isb)
           bs_zeta(j,mibs2use(lsbeam))=esta + &
                nsnccw*nsjdotb*atan2(y_hor,x_hor)/degrad
           if(bs_zeta(j,mibs2use(lsbeam)).lt.0.0_fp) then
              bs_zeta(j,mibs2use(lsbeam))= &
                   bs_zeta(j,mibs2use(lsbeam))+360.0_fp
           end if
           bs_time(j,mibs2use(lsbeam))=tdepay(ind,isb)
           bs_ib(j,mibs2use(lsbeam))=nbscay(ind,isb)
        end if

        !
        !  initially, ptcl gets no goosing in orbit
        goosay(ind,isb)=1.0_fp
        cxpray(ind,isb)=1.0_fp
        fppray(ind,isb)=1.0_fp
        !     for heavy ion beams -- stash the gyrophase angle here:
        !    at low charge states these orbits have large slow gyros that
        !    may require explicit integration in the orbit code
        !
        if(nlheavy(isb)) then
           xzbmay(ind,isb)=atan2(sinpho,cospho)
        end if

        !
        !  the following are calculated later in orball
        !

        if(nfast_map.eq.1) then
           call nbmomry(xiay(ind,isb),thay(ind,isb),0, &
                rmjionay(ind,isb),xzionay(ind,isb),dumj)  ! use bilinear maps
        else
           call nbmomry_xs(xiay(ind,isb),thay(ind,isb),0, &
                rmjionay(ind,isb),xzionay(ind,isb),dumj)  !use xplasma spline maps
        end if
        !mg            rmjionay(ind,isb)=0._fp
        !mg            xzionay(ind,isb)=0._fp
        tfluxay(ind,isb)=0._fp
        xmuay(ind,isb)=0._fp
        pphiay(ind,isb)=0._fp
        pmechay(ind,isb)=0._fp
        !%            if(jxx.eq.ipln2) call nbplot_bdy(0)
     end do
     !%         if(jxx.lt.ipln2) call nbplot_bdy(0)
     !
     !  portion for dep calls leading to new ptcls...
     !   [below, check nnew has not changed]
     !

     nuse_total(isb)=nuse_total(isb)+nnew(isb)
     nuse(isb)=nuse(isb)+nnewp(isb)

     !
     ewghtdep0(isb)=ewghtdep(isb)

     !
     !  additional calls to dep to improve statistics on
     !  dep output:  sorbe (elec. sce)  sorbh (cx therm ion loss term)
     !   sorbn0(cx therm neut sce term) bn00 (1st gen beam neut. dens.)
     !   bpshin (beam shine-thru power)
     !
     !   dep 'global statistics control':
     !  ndep0 defined in transp common-- input variable (ndep0_use used here)
     !  skip if already deposited that many
     !
     !  not applied to fusion products
     !

     if(nlfprod(isb)) then
        if(nnew(isb).eq.0) goto 45
        jdep00=idep00
        goto 445
     end if
     if(nnew(isb).eq.0) goto 45
     if(iltrk0) then
        if(kdep0.ge.mdep0) goto 45
        jdep00=kdep0
     else
        if(idep00.ge.mdep0) goto 45
        jdep00=idep00
     end if
     !
     !   deposit the difference, but dont save ptcls as we already
     !   have enuf to follow in orbit
     !
445  innp1=1
     if(jdep00.lt.mdep0) then
        !additional track data needed for statictics improvment
#ifdef __MPI
        stime71=mpi_wtime()
#endif
        idep0_new=nnew(lsbeam)
        add_track: do while(jdep00.lt.mdep0)

           innp1=idep0_new+1
           idep0_new=idep0_new+ndep0_use

           !-----------------------------------------------
           !        particle splitting (second loop)
           !--------------------------------------------------


           !     split btw processors for deposited ptcls.

           call depall_mpi_split(innp1,idep0_new,nusey_tr)
           if (nusey_tr>0) then
              if(allocated(ninddumn_tr)) deallocate(ninddumn_tr)
              allocate( ninddumn_tr(nusey_tr))
              ninddumn_tr=ninddumn
           end if

           iseed=abs(11111000*nbrandm(idum))

           j=1
           do
              if(j.gt.nusey_tr.or.(jdep00.ge.mdep0)) exit
              jxx=j
              !
              !%            if((ipln1.le.jxx).and.(jxx.le.ipln2)) then
              !%              nbplot=1
              !%            else
              !%              nbplot=0
              !%            end if

              ind1=ninddumn_tr(j) !track id#

              call nbi_random_init_dep(iseed,ind1,iscreen1) !provides unique rng seeds for tracks
              if(nlfprod(isb)) then
                 icount_depall = icount_depall + 1
                 idep00 = idep00 + 1
                 call depfus(ierr)
                 if(ierr.ne.0) then
                    call errorlog('depall: depall>depfus (2)')
                    call errorlog('j = ',j)
                    call errorlog('icount_depall = ',icount_depall)
                    if(myidd.eq.0) then
                       write(nonlin,*) &
                            ' ?depall>depfus (2) j,ierr,# = ',j,ierr, &
                            icount_depall
                    end if
                 end if
              else
                 icount_depall = icount_depall + 1
                 idep00 = idep00+1
                 call nbidep(0,ierr)
                 if(ierr.ne.0) then
                    call errorlog('depall: depall>dep (2) j',j)
                    call errorlog('depall: icount_depall = ', &
                         icount_depall)
                    if(myidd.eq.0)then
                       write(nonlin,*) &
                            ' ?depall>dep (2) j,ierr,# = ',j,ierr, &
                            icount_depall
                    end if
                 end if
              end if

              dep_code_error = dep_code_error + ierr
              !
              !     %            if(jxx.eq.ipln2) call nbplot_bdy(0)
              j=j+1

              jdep00=idep00
              if(.not.nlfprod(isb).and.(iltrk0)) then
                 jdep00=kdep0
              end if
           end do
        end do add_track
#ifdef __MPI
        stime72=mpi_wtime()
#endif

        stime7=stime7+(stime72-stime71)

     end if

     !
45   continue
     !     %         if(jxx.lt.ipln2) call nbplot_bdy(0)

     call depdef
     !+++++++++++++++++++++++++++++++++++++++
     !
     !     jxb torque from deposition  for each beam line (_bl)

     call depdef_bl

     kdep0 = 0
     !+++++++++++++++++++++++++++++++++++++++
     !                call mpi reduce
#ifdef __MPI
     stime81=mpi_wtime()
     ! that combines the values stored by each process into a single value
     ! available to all processes.
     ! needed to control globally all #mc ptcls.
     puseyp(isb)=nuse(isb)
     call depall_buf_reduce
     nuse_total(isb)=puseyp(isb) + 0.5_fp ! global (on all processors) # mc in storage
     stime82=mpi_wtime()
     stime8=stime8+(stime82-stime81)
     call infolog('depall: deposition buf reduced')
     if ((iread.gt.0).and.(numprocd.eq.1)) then
        call depall_buf_read
     end if
#else
     !     this will be a path for nubeam_test, when debugging mpi nubeam
     !     for one ptcl.sublist.
     !     from the working directory read in a
     !     contributions from the other processors
     !     summarized over all prosessors data from the deposition.
     if(debug_mpi.and.nltest_output) then
        if(first_entry) then  !read depall buf only once at the
           !first entry
           call infolog('depall: deposition buf read')
           call nubeam_workpath_retrieve(workpath) ! retrieve workpath, given by nubeam_test
           ! to depall_buf_read

           call depall_buf_read
        end if
     end if
     puseyp(isb)=nuse(isb)
     nuse_total(isb)=puseyp(isb) + 0.5_fp
#endif
     !mg       if(nlfprod(isb)) then
     !mg          zidep0(isb)=1       !total # ptcls.deposited
     !mg       end if
     ! set aep opertaor if available


  end do beamspeceis
  !+++++++++++++++++++++++++++++++++++++++
  !                call mpi gathe dpt% data
#ifdef __MPI
  if((nltrk_dep0.and.mibs_use.gt.0)) then  ! true gather data on beam neutral track
     stime101=mpi_wtime()
     if(.not.allocated(track_xl)) then
        call warnlog( &
             ' depall: nltrk_dep0 set but track_xl not allocated.')
     else
        call infolog( &
             ' depall: before depall_dpt_buf_gather mtrk0', &
             mtrk0)
        call depall_dpt_buf_gather
        call infolog(' depall: after depall_dpt_buf_gather')
     end if
     stime102=mpi_wtime()
     stime10=stime10+(stime102-stime101)
  end if
#endif
  !mg         return
  !---------------------------------------------------------------------
  !  serial code...
  !
  nmc=0.0_fp !storage for summation of #ptcl
  wnmc=0.0_fp!storage for summation of ptcl. weight

#ifdef __MPI
  if(myidd.eq.0.and.nltrk_dep0.and.mibs_use.gt.0) then
     if(allocated(track_sinj)) deallocate(track_sinj)
     allocate(track_sinj(ndep0_use,mibs_use))
     call deptrk_io('track_sinj',1.0_fp,track_sinj,(/1,1/),ierr) !this track data
     !must be normalized
     !have to read them
     call deptrk_io_close(fname_closed)
  end if
#endif

  do isb=1,insbeam
     lsbeam=isb

     iltrk0 = nltrk_dep0.and.(mibs2use(lsbeam).gt.0)

     wghtot = 0.0_fp
     call inispc
     !
     !  compute weight scale up for surviving particles, to compensate for
     !  weight killed by russian roulette in nbkill...

     zbold=bntots(isb)-sforbrr(isb) ! adjust for light wght ptcl kill
     if(zbold.gt.0.0_fp) then
        wghtakil(isb)=zbold/(zbold-wkill(isb))
        sfdeprr(isb)=sfdeprr(isb)+ &
             (zbold-wkill(isb))*(wghtakil(isb)-1.0_fp)
        inkill(isb)=wnkill(isb)+0.5_fp
     else
        wghtakil(isb)=1.0_fp
        inkill(isb)=0
     end if
     !

     key_myidd_0=.false.
     if(myidd.eq.0) key_myidd_0=.true.
#ifndef __MPI
     ! to allow serial code when debugging mpi
     key_myidd_0=.true.
#endif
     if(key_myidd_0) then
        if(max(inuse_prev(isb),nuse_total(isb)).gt.0) then
           write(nonlin,4802) isb,inuse_prev(isb),inkill(isb), &
                nnew(isb),nuse_total(isb)
4802       format(' %depall specie #',i1,1x,' -> ', &
                i9,' -',i9,' (killed) +',i9,' (dep) =',i9,' ptcls.')
        end if
     end if
     !
     if(nnew(isb).eq.0) goto 48
     !
     !  dep ptcl normalizer
     zwrats=zbnew(isb)/max(zidep0(isb),1._fp*nnew(isb)) !zidep0(isb) total # ptcls.deposited
     !if for beam ptcls. zidep0(isb)-
     !has been summed for all cpu at mpi_alreduce
     zwrats_tr=0.0_fp
     if(.not.nlfprod(isb)) then
        if(iltrk0) then !track_data normalizer, increased weigth compare to zwrats
           !beacause only subset data are gathered for track data

           zwrats_tr=zbnew(isb)/zitrk0(isb)
        end if
     end if
     zwgrat=ennew(isb)/ewghtdep(isb)
     !  factor for renormalizing weight of new ptcls in orbit loop
     wghtsave(isb)=ennew(isb)/ewghtdep0(isb)
     zwgfac=wghtsave(isb)/zwgrat
     if(iltrk0.and.il_birth) then
        bs_wght(:,mibs2use(lsbeam))=bs_wght(:,mibs2use(lsbeam))* &
             wghtsave(isb)
     end if

     !
     goto 49
     !
48   continue
     !
     !  pass thru here if no new ptcls are deposited; depnrm still needs to
     !  be called
     !
     zwgrat=0.0_fp
     zwrats=0.0_fp
     zwrats_tr=0.0_fp
     zwgfac=1.0_fp
     wghtsave(isb)=0.0_fp
     !
     !  proceed to renormalize dep output
     !
49   continue
     call depnrm(zwgrat,zwrats,zwrats_tr,zwgfac)
     !
     if(weinjav(isb).le.0.0_fp) then
        einjav(isb)=zeinjav(isb)
     else
        einjav(isb)=(1.0-zfnew(isb))*zeinjav(isb) + &
             zfnew(isb)*einjav(isb)/weinjav(isb)
     end if
     vmins(isb)=0.01_fp*sqrt(einjav(isb)/v2toev)
     !
     !--------------------
     !
     !  evaluate profile of average particle weights
     !
     !  if no new ptcls are deposited, let the number of m.c. ptcls decay ...
     !
     !
     !  nmc & wnmc sums moved to depall.for 09/04/09 by mg
     !
     !mpi part
     do iptcl=1,nuse(lsbeam)
        ind=nbndex(iptcl,lsbeam)
        ind=find_ptcl_ind(numprocd,myidd,ind)
        jloc=njavay(ind,lsbeam)
        nmc(jloc,lsbeam)=nmc(jloc,lsbeam)+1.0_fp
        tdepay_ind=max(0.0_fp,tdepay(ind,lsbeam)-tbm1)
        if(tdepay_ind.gt.0.0_fp) then !new orbit
           if(goosay(ind,lsbeam).ge.1.00001_fp)then
              call errorlog( &
                   ' depall: dep time for an old ptcl is wrong, t_dep= ', &
                   tdepay(ind,lsbeam))
              call errorlog(' depall: ind= ',ind)
              call errorlog(' depall: tbm1=',tbm1)
              dep_code_error = dep_code_error + 1
           end if
           wnmc(jloc,lsbeam)=wnmc(jloc,lsbeam)+wghtay(ind,lsbeam) &
                *wghtsave(lsbeam) ! rescale weight
           wghtot = wghtot + wghtay(ind,lsbeam)*wghtsave(lsbeam)
           !injected power at mc particle "birth" = injected power  - shine through power
           pinjay(ind,lsbeam)=pinjs(lsbeam)- bpshins(lsbeam)
        else
           wnmc(jloc,lsbeam)=wnmc(jloc,lsbeam)+wghtay(ind,lsbeam)
        end if
     end do
     ! save power at mc deposition time
     if(pinjs(lsbeam).gt.0.0_fp) then
        pinjs_res(lsbeam)=pinjs(lsbeam)-bpshins(lsbeam)
     end if

     !mpi part, reduce ptcl # and weight
#ifdef __MPI
     call nbi_ptcl_wght_reduce
#endif
     !the sum of actually deposited for subsequent orbiting, divided by
     !the number of mc ions actually deposited
     if(nnew(lsbeam).gt.0) then
        wghtot=wghtot/nnew(lsbeam)
     end if

     key_myidd_0=.false.
     if(myidd.eq.0) key_myidd_0=.true.
#ifdef __MPI
#else
     ! to allow serial code when debugging mpi
     key_myidd_0=.true.
#endif
     if(key_myidd_0) then

        call nmcsum

        call get_wbav

     end if


#ifdef __MPI
     call nbi_ptcl_wght_bcast
#endif


     !  follow particles in orbit and delete particles which
     !  thermalize, leave plasma, etc.
     !
     !-------
     !  ptcl balance code
     zwght3=wghtdep(isb)*wghtsave(isb)
     !  deposition of ptcls defined to include prompt loss if evaluated here
     sfdepmc(isb)=zwght3/delt + bxflws(isb)/delt
     xninjea(isb)=zwrats*max(zidep0(isb),1._fp*nnew(isb))/delt - &
          zwght3/delt - sfshine(isb)
     sfdepba(isb)=xninjs(isb)-sfshine(isb)- &
          xninjea(isb)-sfdepmc(isb)

     !
     !  end of beam species loop
     !
     !xx         call nbi_zcheck_dep
     !xx         call nbi_pput_dep
     !xx         call nbi_pget_dep
     !
     ! compute approximate slowing down time to be used for next timestep
     !
     call nbi_tausl_vol_avg
     !
  end do !species
  if(il_birth)then
     call nbi_bsource_io    !birth source to be written on each cpu
     !and nbs_step will be updated there
  end if

  if(nltrk_dep0.and.mibs_use.gt.0) then

#ifdef __MPI
     if(myidd.eq.0) then
        call deptrk_io('track_sinj',track_sinj,ierr) !this track data
        !has been normalized
        !have to write them back to cdf file
        call deptrk_io_close(fname_closed)
        if(allocated(track_sinj)) deallocate(track_sinj)
     end if
#else
     call nbi_wrdeptrk      !write deptrk data to netcdf file <runid>_dpt.cdf,
     !track data arrays will be deallocted there
#endif

  end if

  !
  zidep0 = 0.0_fp
  zitrk0 = 0.0_fp
  ierr=maxval(dep_code_error)+0.5_fp
  if(ierr.gt.0) then
     call errorlog(' emergency exit depall: ierr > 0')
     return
  end if
  !
  !-----------------------------------------------------------------------
  !  dmc june 1992 -- option to use dep to compute a fast neutral density
  !  in a 3d box along a selected beam axis

  !
  kdep0 = -1 ! do not gather "track" data
  ! if data into beam-in-box wanted and it is a time for it
  if(nbbox.gt.0) then
     xksid=0.0_fp !clean up is taken
     call depbox(ierr) !also can collect neutral density in a tube/s
     if(ierr.ne.0) then
        call errorlog(' emergency exit depall: depbox error')
        return
     end if
  else if (nbtube>0)  then
     call deptube(ierr)     ! rga, neutral density in a tube
     if(ierr.ne.0) then
        call errorlog(' emergency exiting depall: deptube error')
        return
     end if
  end if
  kdep0 = 0
  !
  !-----------------------------------------------------------------------
  !
  nbs=0
  lsbeam=0
  abeam=0.0_fp
  xzbeami=0.0_fp
  !
  !  restore expanded beam dist fcn zone grid for orbitting
  !
  call nbi_rbloat
  !
  !%      nbplot=ibplot
  !
  !  smooth output radial profiles
  !
  call depsmoo
  !
  !  get integer :: counts of mc ptcls lost
  !
  nmcloss = xmcloss_d + 0.5_fp
  !
  !  zero out nb*sigma*v sums (beam-beam interactions used in dep)
  !
  call nbsv00
  !
  call cpu_time(ztim2)

#ifdef __MPI
  ttime22=mpi_wtime()
  if (myidd.eq.0) cpbdep = cpbdep + (ttime22-ttime21)
#else
  cpbdep = cpbdep + (ztim2-ztim1)
#endif
  first_entry=.false.
  call infolog('exiting depall')
  return
end subroutine depall
