!----------------------------------------------------------------------
!  orball -- do all monte carlo orbit loops (prev. part of fimain)
!
!   mg  march 2009 -- for mpi each process carries a single
!                     sublist of particles  instead of both sublist and
!                     whole list arrays.
!   dmc apr 2007 -- adding ability to repeat the advance of all orbits
!      by saving and retrieving the ptcl variables.  also: nlbout=t
!      i/o feature removed (was not working with mpi in any case).
!
!   dmc may 2001 -- commenting our rf resonant operator (defunct)
!      (the fpp code interface is being rewritten-- the code here
!       would also have to be redone.  to do this would require
!       evidence of use)...
!
!   dmc aug 1988 -- support multispecies beams
!
subroutine orball(iorbtot,ierr)
  !
  use nbi_com
  use nbi_random_pll
  use nbi_random
  use nbi_mpi_com
  use nbi_fbm_random
  use orbrzv, only : mcrf_on,rfdbg_pdiag,orbrzv_nrf,orbrzv_out
  use logmod
  use dep_mpi_split ! use this module to deposit cx tracks for statistics improvement
  ! does not interfer with depall
  use nbi_aep_mod, only : aep_init, rbq_init, id_ptcl, nbi_aep_clean
  use nbi_aep_mod, only : nbi_aep_write, nbi_rbq_write, imap_nkaep_use
  use nbi_aep_com, only : ilaep_init

  !
  implicit none

  integer, intent(out) :: iorbtot   ! total no. of orbits calculated.
  integer, intent(out) :: ierr      ! status code (0=ok)
  !
  !  local stuff...
  !
  integer :: isb,iorb,j,irecap,ind,jptr,ibflg,ixiav
  integer :: iersum,ib,ien,isplit,jsplit,indsave
  integer :: iuncon,itrap,jloc,idum
  integer :: inuse(mibs)    !original ptcl in use
  integer :: inuse_p(mibs)    !original ptcl in use on given processor
  real(fp) zwght1,zang,zxhw,zphiav,zeiav,zxksav,xjadep,zrrf
  real(fp) zadjust,zwght2,zw,zav,zvar2,znlimdt
  real(fp) zpw,zfrnm, zprft, ztprft
  integer :: i,nuseycx,icx,ierr_aep
  !-------------------------------
  !  max energy for time step adjust normalization

  real(fp), parameter :: zmaxe_dtc = 500.0e+3_fp

  !-------------------------------
  !  dmc apr 2007 -- testing rerun of orbit loop, controlled from within
  !  orball; results should be identical...
  !
  !  set the following variable to the number of times to repeat the calculation
  !  (normal value: 1, just do the calculation once).
  !
  !  calculation may be redone for other reasons (e.g. adjust rf wave field).
  !
  integer, parameter :: ipass_num = 10  ! >1 for debug only
  integer, parameter :: ipass_num2 = 22  ! >1 for debug only
  real(fp), parameter :: zpwr_errmax = 0.1_fp
  integer :: ipass,i_repeat
  logical :: irepeat_possible, irepeat, itmp_file, irepeat_lastone
  logical :: key_myidd_0
  real(fp) :: zpw2(ipass_num2),z_frnm(ipass_num2) ! save abs error and renormalization
  ! factor for rf operator, in case
  ! in all itteration did not reach
  ! 10% error in power (zpwr_errmax = 0.1_fp)
  integer, dimension(1) :: i_pass
  !-------------------------------

#ifdef __DEBUG
  real(fp) zwdep,zwdepp,zwold,zwoldp
#endif
  !
  logical ilrecap,ifirst,ilsplit,inew
  !
  !%      integer :: ibplot
  !%      logical lyes,inwonl
  !
  !  mpi related stuff
  !
  integer :: mpierr,iseed,indnew,indcx
  integer :: iscreen1,impitime
  !
  character(len=200) :: zorball_tmp_filename
  !
  real(fp) :: ztim1,ztim2,zrandom
  real(fp) :: ztime1,ztime2,ztime3,ztime4
  integer :: find_ptcl_ind
  integer, dimension(:), allocatable ::nuse_total! total # of mc ptcls currently in storage
  ! is using for a global control during deposition
  ! and orbiting
  integer :: iost
  real(fp) :: wbav_c, p_inj_rs,fbmcxsum
  integer :: ncx1
  real(fp) :: xncx0_use
  logical :: mcrf_on_save
  real(fp) :: tdepay_ind,t_depay
  integer :: iseed_aep
  !
  call infolog(' entering orball')
  if(allocated(nuse_total)) deallocate(nuse_total)
  allocate(nuse_total(mibs)); nuse_total = puseyp + 0.5_fp
  !
  !------------------------------------------------------------------
  !
  ierr=0

  itmp_file = .false.

  impitime=1
  ttime31=0.0_fp
  ttime32=0.0_fp
#ifdef __MPI
  ttime31=mpi_wtime()
#endif
  if(only_io) then
     call errorlog( &
          'orball:orball(nubeam): improper call, only_io = .true.')
     if(myidd.eq.0) then
        write(nonlin,*) &
             ' ?orball(nubeam): improper call, only_io = .true.'
     end if
     ierr=666
     call errmsg_exit(' ?orball(nubeam): improper call, only_io=.t')
     goto 9999
  end if

  call cpu_time(ztim1)

  !  sawtooth related setups (dmc nov. 2005) (serial code)

  sawflag = sawflag.or.sawflag2  ! true to use pre-calculated operator

  !  pre_sawnbi_setup was called by nbstart to initialize calculation of an
  !  operator in the upcoming step...

  !
  !  see if this is the first step (all ptcls new)
  !
  ifirst=.true.
  do isb=1,nsbeam
     ifirst = ifirst.and.(nuse_total(isb).eq.nnew(isb))
  end do
  !
  !  clear sums:  beam ion population energy, momentum, and n
  !%      ibplot=nbplot
  !%        if(nbplot.gt.0) call nbplot_bdy(0)
  !%      call uread('orball:  plot only newly deposited orbits ? (y/n)$')
  !%      inwonl=lyes(0)
  !%      if(inwonl) nbplot=-1
  !
  norbtot=0
  !
  xnum_dt_adj0=0.0_fp
  dt_adj0_sums=0.0_fp
  !
  iersum = 0
  stime9=0.0_fp
  !
  iscreen1=1600+myidd
  mcrf_on_save=mcrf_on

  !     allocate arrays for splitting with default dimension
  call nbi_alloc_splt
  !  save array dimension for slitting, may be changed in nbi_spltexpnd.f90
  minbs_expnd=minbs

  !c init mhd activity now
  if(allocated(ilaep_init)) deallocate(ilaep_init)
  allocate(ilaep_init(nsbeam))
  ilaep_init=.false.
  if(maxval(abs(aep_ascale)).ne.0.0_fp) then 
     if(minval(aep_ascale).lt.0.0_fp)  then !RBQ
        call rbq_init(ierr_aep)
        nchdvp=1
     else
        call aep_init(ierr_aep)
     endif
     if(ierr_aep.eq.0) then
        call infolog(' aep active')
        do isb=1,nsbeam
           if(maxval(imap_nkaep_use(1:nkaep,isb)).gt.0) then
              ilaep_init(isb)=.true.
           endif
        enddo
     end if
  end if

  do isb=1,nsbeam
     !
     xncx0_orb(isb)=0._fp
     xncx0=0._fp
     inuse(isb) = nuse_total(isb)
     inuse_p(isb) = nuse(isb) !save nuse on previous step

     lsbeam=isb


     if(orbrzv_rf_option.eq.1 .and. mcrf_on_save) then
        if(orbrzv_toric_prftot(lsbeam) < 4.0e2_fp) then
           mcrf_on=.false. !do not apply rf operator if not enough rf power
        end if
        zpw2(:)=epsinv
        z_frnm(:)=1.0_fp
     end if
     ttcx=0.0_fp
     tscx=0.0_fp

     if(nlfbmfpp(lsbeam)) cycle
     !
     call inispc
     !
     eirefs(isb)=min(zmaxe_dtc,einjs(isb))
     xmrefs(isb)=raxis*abeams(isb)*1.6726e-24_fp* &
          sqrt(eirefs(isb)/v2toev)
     !
     if(.not.nlheavy(lsbeam)) call jbshld
     !
     zwght1=bntots(lsbeam)
     !
     enolds(lsbeam)=0.0_fp
     bpholds(lsbeam)=0.0_fp
     bntots(lsbeam)=0.0_fp
     wnlight(lsbeam)=0.0_fp



     !
     !  check that all orb outputs are zero now
     !
     call nbi_zcheck_orb
     call nbi_zcheck_ntr

     !  (parallelized) sawtooth operator...

     if(sawflag) then
        call sawnbi_exec
     end if

     !
     !  partial seed for ptcl rng sequences
     !
     zrandom=nbrandm(idum)
     iseed=abs(10000000*zrandom)


     !-----------------------------------------------
     !  set up for possible recalculation of orbit loop
     !  add logic here if recalculation due to rf is a possibility...

     !  see also the bottom of the loop where logic decides whether or not
     !  to actually repeat the calculation

     irepeat = .false.  ! repeat of calculation not decided yet...
     irepeat_lastone = .false.  ! repeat of calculation not decided yet...
     i_repeat=0
     itmp_file = .false.
     if(orbrzv_rf_option.eq.1 .and. mcrf_on) then
        irepeat_possible = ((ipass_num.gt.1))
     else
        irepeat_possible = ((ipass_num.gt.1).and.nltest_output)
     end if

     !----------------
     if(irepeat_possible) then
        call orball_loclist_write  ! save local ptcl list in temporary file
        itmp_file = .true.
     end if
     ipass = 0
     !-----------------------------------------------
     !  loop here if re-doing all orbits for a species
     !   (e.g. to recalculate after an rf wave field renormalization)

50   continue
     ipass = ipass + 1
     if(ipass.gt.1) then
        call nbi_zzero_orb       ! zero out orb mpi-reduced outputs...
        call nbi_zzero_ntr       ! zero out orb mpi-reduced outputs...
        call orball_loclist_read ! restore original ptcl list
        if(nlbout) then
           mlost_data=0
           call nbi_zzero_lst
        end if
        fbm(:,:,:,lsbeam)=0.0_fp
        if(nlfbmflr) then
           fbm_ptcl(:,:,:,lsbeam,:)=0.0_fp
        end if !nlfbmflr
        fbm_cx(:,:,:,lsbeam)=0.0_fp
        xncx0_orb(lsbeam)=0.0_fp
     end if

     !
     iorb=0
     !
     !      keep track of currently stored energy and ptcls
     !
     j=0
     nbsplt=0

#ifdef __DEBUG
     zwdepp = 0
     zwdep = 0
     zwoldp = 0
     zwold = 0
#endif

     !  cx recapture monte carlo ptcls counter
     irecap=0
     !  orbit splitting counters
     isplit=0
     jsplit=0
     ilsplit=.false.

     !  (parallelized) compression operator...
     if(nlcprb) then
        call bcmprsn(ierr)
        if(ierr.ne.0) then
           orb_code_error(lsbeam) = orb_code_error(lsbeam) + 1
           call errmsg_exit(' ?orball: error in bcmprsn' )
           goto 9999
        end if
     end if
     !--------------------------------------------------------------
     !  top of orbit loop
     !
60   continue

     !
     if(.not.ilsplit) then
        j=j+1
        !            call warnlog('not ilsplit')
        !  exit test
        if(j.gt.(nuse(isb)-irecap-isplit)) then
           goto 90
        end if

     end if
     !
     !
     !  next ptcl (from mainstream ptcl array)
     !    (note nbndex has been replaced by nbndexn)
     !    original: ind=nbndex(j,lsbeam)

     ind=nbndex(j,lsbeam)
     ind=find_ptcl_ind(numprocd,myidd,ind) !find proper ptcl.id#
     indsave=ind
     pinj_ay=pinjay(ind,lsbeam)
     nbsce=nbscay(ind,lsbeam)
     nbien=nbienay(ind,lsbeam)
     einj0=einjay(ind,lsbeam)
     xion=xiay(ind,lsbeam)
     th=thay(ind,lsbeam)
     vion=vay(ind,lsbeam)
     xksid=xksidy(ind,lsbeam)
     t_depay=tdepay(ind,lsbeam)
     wght=wghtay(ind,lsbeam)
     if(wght.le.0.0_fp) then
        call errorlog('orball: wght=0.0')
        call errorlog('orball: probable book keeping error !')
        call errorlog('lsbeam = ',lsbeam)
        call errorlog('ind = ',ind)
        call errorlog('nbndex(j,lsbeam) = ',nbndex(j,lsbeam))
        call errorlog('nuse_total(lsbeam) = ',nuse_total(lsbeam))
        call errorlog('nuse(lsbeam) = ',nuse(lsbeam))
        if(myidd.eq.0) then
           write(nonlin,*) ' ?? orball lsbeam = ',lsbeam,' ind = ' &
                ,ind,nbndex(j,lsbeam)
           write(nonlin,*) ' ?? orball nuse_total = ', &
                nuse_total(lsbeam),nuse(lsbeam)
           write(nonlin,*) ' ?? orball myidd = ',myidd,' j = ',j, &
                ' wght=0.0, probable book keeping error!'
           write(nonlin,*)nbndex(1:nuse(lsbeam),lsbeam)
        end if
        call errorlog(' emergency exiting orball')
        call bad_exit
     end if

     goose=goosay(ind,lsbeam)

     dt_adjust = dt_adj0s(lsbeam)
     if(dt_adjust.le.0.0_fp) dt_adjust=1.0_fp

     cxskip=cxpray(ind,lsbeam)
     fpskip=fppray(ind,lsbeam)
     xistate=xistay(ind,lsbeam)
     !
     !  launch time -- .gt.0 means a new orbit...
     !
     ! mg 11/2013
     !         bitime(1)=tdepay(ind,lsbeam)
     tdepay_ind=max(0.0_fp,tdepay(ind,lsbeam)-tbm1)
     bitime(1)=tdepay_ind
     if(bitime(1).eq.0.0_fp) then
        !  old orbit
        !
        !  as a precaution (e.g. due to changed mhd equilibrium) reduce
        !  goosing by 50% for old particles.  reduce by more if the macro
        !  beam timestep was shortened sufficiently (cf tfracp from nbstart).
        !

        goose=max(1.0001_fp,min(0.5_fp,tfracp)*goose)
        if(.not.ilsplit) then
           wght=wght*wghtakil(lsbeam) ! compensate for depall->nbkill rr
           wghtay(ind,lsbeam)=wght
           pinj_ay = pinj_ay*wghtakil(lsbeam)
           pinjay(ind,isb)=pinj_ay

        end if
        inew=.false.
#ifdef __DEBUG
        if(.not.ilsplit) then
           zwdepp = zwdep
           zwoldp = zwold
           zwold = zwold + wght
        end if
#endif
     else
        !  new orbit
        wght=wght*wghtsave(lsbeam) ! rescale weight
        wghtay(ind,lsbeam)=wght
        pinjay(ind,isb)=pinj_ay
        inew=.true.
#ifdef __DEBUG
        zwdepp = zwdep
        zwoldp = zwold
        zwdep = zwdep + wght
#endif
     end if
     !
     jloc=njavay(ind,lsbeam)
     !
     ! mhr 19-july-1994
     !
     srd=0.0_fp
     !
     indnew=lsbeam*15111100+ind  !int((xion+xksid)*2.e+8_fp)
     iscreen1=1600+myidd
     call nbi_random_init_pll(iseed,indnew,iscreen1)
     !
     if(nlheavy(lsbeam)) then
        if(bitime(1).gt.0.0_fp) then
           xzbeami=1.0_fp  ! new orbit
           zang=xzbmay(ind,lsbeam)  ! restore gyrophase
           cospho=cos(zang)
           sinpho=sin(zang)
           bitime(1)=min(bitime(1),0.999_fp*delt) ! force some steps
        else
           xzbeami=xzbmay(ind,lsbeam)
        end if
        call jbshld
     end if


     !
     !%          if(inwonl.and.(bitime(1).gt.0.0)) then
     !%             nbplot=ibplot
     !%          end if
     !
     !-------------------------------------------
     !  check for splitting based on weight/<avg weight> -- old ions only.
     !
     ilsplit=.false.
     wbav_c=wbav(jloc,lsbeam)
     fact_wbav = 1
     ! if the code went from a period of time with low power to a period of time with high power,
     ! low weight particles from the low power regime would persist;
     ! if their presence influences wbav it could cause new particles from the high power
     ! regime to be split;
     ! the same might happen to old particles if transitioning from a period of high power to a
     ! prolonged period of low power.
     !
     !
     p_inj_rs=pinjs_res(lsbeam)   ! injected beam power without
     ! shine through power at the current time step
     ! or at the time step when power wasn't equal zero
     if(.not.inew) then
        !
        !  old orbit eligible for splitting...
        !
        !     old paricles. do not let them make a lot of splittings. increase wbav.
        if(p_inj_rs.gt.0.0_fp.and. &
             (pinjay(ind,isb).gt.p_inj_rs)) then
           !     power drops
           !     old paricles from high power regime
           fact_wbav = pinjay(ind,isb)/p_inj_rs
           wbav_c=wbav_c*fact_wbav
        end if
        !
        ! dmc apr 2011 -- force wbav to be scaled to local wnmc/nmc
        !   this is for safety; have observed situations otherwise where there is
        !   so much mc w/<w> > 2 splitting that the particle lists overflow
        !   in one case this was due to an end-of-run shortened time step which
        !   led to less new weight deposited and the new weight deposited set
        !   wbav_c to a number that was too small.
        !
        if(nmc(jloc,lsbeam).gt.0.0) then
           wbav_c = max(wbav_c,wbav_kill(jloc,lsbeam))

           !     mg, 09.01.2011 i increased  wbav_c to lower # of splittings
           !     mg old line   >            (0.8_fp*wnmc(jloc,lsbeam)/nmc(jloc,lsbeam)) )
        end if
        if(wght.gt.2*wbav_c) then
           !
           !  w/<w> >= 2
           !

           if(nuse(isb).lt.ntotptp(isb)) then
              !
              !  slot available
              !
              ilsplit=.true.
              isplit=isplit+1
              jsplit=jsplit+1
              ! save for later...
              wghtay(ind,lsbeam)= &
                   wghtay(ind,lsbeam)-wbav_c ! save for later
              wght=wbav_c   ! make slot for this now...
              !
              nuse(isb)=nuse(isb)+1
              jptr=nuse(isb)
              ind=nbndex(nuse(isb),isb)
              ind=find_ptcl_ind(numprocd,myidd,ind) !find proper ptcl.id#
              !
              pinjay(ind,lsbeam)=pinj_ay !pass to "child" from "parent"
              nbscay(ind,lsbeam)=nbsce
              nbienay(ind,lsbeam)=nbien
              einjay(ind,lsbeam)=einj0
              xiay(ind,lsbeam)=xion
              thay(ind,lsbeam)=th
              vay(ind,lsbeam)=vion
              xksidy(ind,lsbeam)=xksid
              wghtay(ind,lsbeam)=wght
              tdepay(ind,lsbeam)=t_depay
           end if
        end if
     else
        !     new paricles.define the renormalization factor for wbav for time dependence of power included to
        !
        if(pinjay(ind,isb).gt.0.0_fp.and. &
             (p_inj_rs.gt.pinjay(ind,isb))) then
           !     power up
           fact_wbav = p_inj_rs/pinjay(ind,isb)
        end if
        !

     end if
     !-------------------------------------------
     !  sums formerly calculated in serial depall have been moved here inside
     !  the parallel section...
     !
     if(.not.nlfprod(lsbeam)) then
        if(nbsce.eq.0) then
           call errorlog('orball:')
           call errorlog('beam source tracking algorithm error !')
           if(myidd.eq.0) then
              write(nonlin,*) &
                   'orball: beam source tracking algorithm error!'
           end if
           orb_code_error(lsbeam) = orb_code_error(lsbeam) + 1
        else
           zw=wght*(delt-bitime(1))
           wts_ien(nbien,nbsce)=wts_ien(nbien,nbsce)+zw
           ein_ien(nbien,nbsce)=ein_ien(nbien,nbsce)+zw*einj0
           dei_ien(nbien,nbsce)=dei_ien(nbien,nbsce)+zw*einj0*einj0
        end if
     end if
     !
     !  nmc & wnmc sums moved to depall.for 09/04/09  mg
     !
     !         nmc(jloc,lsbeam)=nmc(jloc,lsbeam)+1.0_fp
     !         wnmc(jloc,lsbeam)=wnmc(jloc,lsbeam)+wght
     !
     !------------------------
     !  clear splitting flags
     !
     ilrecap=.false.
     if(.not.ilsplit) then
        jptr=j
     end if
     !
     !  orbit
     !
     if(ilaep_init(lsbeam)) then
        id_ptcl=ind         ! set nbi_aep_mod, for aep operator
        iseed_aep=iseed+2
        call nbi_random_init_aep(iseed_aep,indnew,iscreen1) !provides unique rng seeds for tracks
     end if
65   continue
     !
#ifdef __DEBUG
     if(.not.(ilsplit.or.ilrecap)) then
        call zbal_debug(iorb,zwdepp,zwoldp)
     end if
#endif
     iorb=iorb+1

     !%         write(nonlin,*) 'lsbeam, nuse(lsbeam), iorb: ',
     !%     >        lsbeam,nuse(lsbeam),iorb

     call orbit(iorb,ierr)
     orb_code_error(lsbeam) = orb_code_error(lsbeam) + ierr
     if(ierr.gt.0) then
        call errmsg_exit &
             (' ?orball: orbit failed')
     end if

     if(xicur.eq.epsinv) then
        call infolog('orball: orbit never inside plasma,iorb=',iorb)
        if(myidd.eq.0) then
           write(nonlin,9011) myidd,iorb
        end if
9011    format( &
             ' %orball:  in processor ',i5,': orbit # iorb=',i10, &
             ' never inside plasma.')

     end if
     !
     !%          if(inwonl) then
     !%             if(nbplot.gt.0) call nbplot_bdy(0)
     !%             nbplot=-1
     !%          end if
     !
     !  dmc 27 sep 1994 -- new radial current calculation, based on orbit
     !  average position -- not-lost orbits only.  lost orbits had xistop
     !  set in subroutine orbit
     !
     ibflg=nbflg
     if(nbflg.eq.1) then
        if(v2toev*vioncur*vioncur.lt.100.0_fp) ibflg=2
     end if
     if(ibflg.eq.1) then
        !
        !  use last inside-the-plasma orbit coordinates
        !
        xistop=xjadep(xicur,thcur,rmjcur,ycur,vioncur,vpllcur, &
             itrap,iuncon,zxksav,zeiav,zphiav,ixiav,zxhw)
        !
        !  dmc 18 sep 1995 -- trapping fraction calculation
        !
        call addtrap(wght,itrap,zeiav)
        !
        !  end dmc 18 sep
        !
     else if(ibflg.eq.2) then
        !
        xistop=xicur                ! low energy
        !
     end if
     !
     !  only for cx loss is radial current already counted; for normal orbits
     !  and other losses add to orbit current here.
     !
     if(nbicx.eq.0) then
        call radcur(xistate,xistop,wght,tirorbs)
     end if
     !
     if(ibflg.eq.1) then
        !  zjav is avg radial location nbndex calc. subr. orbit
        jloc=xjav+0.5_fp
        jloc=max(lcentr,min(lep1,jloc))
        zrrf=wght/wbav_kill(jloc,lsbeam)
        if(zrrf.le.1.0e-4_fp) then
           ! mg 07/10/2013            if(zrrf.le.0.33333333333333_fp) then
           ! mg must be consistent with nbkill_pll.for
           wnlight(lsbeam)=wnlight(lsbeam) + 1.0_fp
        end if
     end if
     !
     !  ptcl survive to next step?
     if(ibflg.eq.1) then
        !  yes:  save parameters
        !feb92            ztsum3=ztsum3+wght
        !  (for next step compression operator)
        !
        if(nlcprb) then
           call ptclsave(xion,th,vion,xksid, &
                rmjionay(ind,lsbeam),xzionay(ind,lsbeam), &
                tfluxay(ind,lsbeam), &
                xmuay(ind,lsbeam),pphiay(ind,lsbeam), &
                pmechay(ind,lsbeam))
        end if
        !
        !  traditional state arrays
        !
        njavay(ind,lsbeam)=jloc
        xiay(ind,lsbeam)=xion
        thay(ind,lsbeam)=th
        vay(ind,lsbeam)=vion
        xksidy(ind,lsbeam)=xksid
        goosay(ind,lsbeam)=goose
        cxpray(ind,lsbeam)=cxskip
        fppray(ind,lsbeam)=fpskip
        wghtay(ind,lsbeam)=wght
        xistay(ind,lsbeam)=xistop
        !     commented by mg on 11/2013
        !            tdepay(ind,lsbeam)=0.0_fp ! old orbit now
        !
        if(nlheavy(lsbeam)) xzbmay(ind,lsbeam)=xzbeami
        !
        !  compute sawtooth operator -- shadow ptcl list
        !
        call pre_sawnbi(ind)
        !

        !*********************************************************
        !
        !***********************************************************

        !
        !  sum of surviving ptcls, energy
        enolds(lsbeam)=enolds(lsbeam)+ep*v2toev*wghtay(ind,lsbeam)
        bpholds(lsbeam)=bpholds(lsbeam)+bphion*wghtay(ind,lsbeam)
        bntots(lsbeam)=bntots(lsbeam)+wghtay(ind,lsbeam)
        !
     else
        !  no:   delete; create free ptcl slot

        wghtay(ind,lsbeam)=0.0_fp

        call idelet(nbndex(1,isb),minb,nuse(lsbeam),jptr)

        !  decrement cx-recapture ion counter if needed
        if(ilrecap) then
           irecap=irecap-1
        else if(ilsplit) then
           isplit=isplit-1
        else
           j=j-1
        end if
     end if
     !
     !  process pending cx-recaptured ions now if space permits (and there
     !  are such ions)
     !
     if((nbsplt.eq.0).or.(ilsplit)) then !no splitting occurred
        goto 60
     else
        if(nuse(isb).ge.ntotptp(isb)) then !splitting occurred but no room left
           call errorlog('orball: error nbsplt gt.0 but no room left')
           if(myidd.eq.0) then
              write(nonlin,*) &
                   ' ?orball(nubeam):nbsplt>0 but nuse(isb)>=ntotptp(isb)'
           end if
           ierr=1
           orb_code_error(lsbeam) = orb_code_error(lsbeam) + ierr
           call errmsg_exit &
                ("orball:nbsplt>0 but nuse(isb)>=ntotptp(isb)")
        end if
     end if

     !
     !  insert cx-recapture ion at end of ptcls list
     !   keep track of number of cx-recapture ions

     !
     ilrecap=.true.
     irecap=irecap+1
     !
     !  insert previously split ptcl back into mainstream
     bitime(1)=tbsa(nbsplt)
     goose=1.0001_fp  ! 1.0000 reserved for newly deposited ions

     dt_adjust = dt_adj0s(lsbeam)
     if(dt_adjust.le.0.0_fp) dt_adjust=1.0_fp

     cxskip=1.0_fp
     fpskip=1.0_fp
     pinj_ay=pinjbsa(nbsplt)
     nbsce=nbsbsa(nbsplt)
     nbien=nbebsa(nbsplt)
     einj0=einjbsa(nbsplt)
     xion=rbsa(nbsplt)
     th=thbsa(nbsplt)
     vion=vbsa(nbsplt)
     xksid=xksbsa(nbsplt)
     xistate=xisbsa(nbsplt)
     srd=0.0_fp
     wght=wtbsa(nbsplt)

     nbsplt=nbsplt-1
     !  create new ptcl slot in mainstream array
     ! ---->> examine following 3 lines carefully

     nuse(isb)=nuse(isb)+1
     jptr=nuse(isb)
     ind=nbndex(nuse(isb),isb)
     ind=find_ptcl_ind(numprocd,myidd,ind) !find proper ptcl.id#

     nbscay(ind,lsbeam)=nbsce
     nbienay(ind,lsbeam)=nbien
     einjay(ind,lsbeam)=einj0
     pinjay(ind,lsbeam)=pinj_ay
     xiay(ind,lsbeam)=xion
     thay(ind,lsbeam)=th
     vay(ind,lsbeam)=vion
     xksidy(ind,lsbeam)=xksid
     wghtay(ind,lsbeam)=wght
     tdepay(ind,lsbeam)=t_depay

     !
     !  follow orbit of this ion
     !
     goto 65
     !
90   continue

     !
     !---------------------------
     !  done following orbits of current species
     !
     xnorbs(lsbeam) = iorb

     !***************************************************
     !  mpi related code
     !  (buffer reduction and particle collection)
     !***************************************************
     j_spl=jsplit
     ! mhd activity now
     if(ilaep_init(lsbeam).and.(maxval(abs(aep_ascale)).ne.0.0)) then    
#ifdef __MPI
        call nbi_aep_reduce
#endif
        if(myidd.eq.0) then
           if(minval(aep_ascale).lt.0.0) then
               call nbi_rbq_write
          else
              call nbi_aep_write
           endif
        end if
     end if

#ifdef __MPI
     stime91=mpi_wtime()
     !
     ! that combines the values stored by each process into a single value
     ! available to master processes.(see nubeam_buf_reduce for mpi_reduce)
     !
     puseyp(lsbeam)=nuse(lsbeam)
     call nubeam_buf_reduce !ntr% arrays will not be reduce now
     if(myidd.eq.0) then
        nuse_total(lsbeam)=puseyp(lsbeam) + 0.5_fp !it is a global # mc in storage
        ! on master processor(!)
        ! because mpi_reduce not mpi_allreduce
        ! (see nubeam_buf_reduce).
        ! this value for a child processor will be
        ! updated before the next entry to depall
        ! by broadcasting of the states.
     end if

     stime92=mpi_wtime()

     stime9=stime9+(stime92-stime91)
#else
     puseyp(lsbeam)=nuse(lsbeam)
     nuse_total(lsbeam)=puseyp(lsbeam) + 0.5_fp
#endif

     !---------------------------------------------------
     ! ** decide here whether to redo orbits **
     !    (insert logic here for redo with rf wavefield renormalization)
     !    note: this needs mpi modifications so that all cpus do the same
     !    thing based on (presumably) a test of mpi-reduced rf results).
     !

     if(myidd.eq.0) then
        if(ipass.gt.ipass_num) then
           irepeat = .false.
        else
           irepeat = .true.
        end if
        if(orbrzv_prftot(lsbeam).lt.0.0_fp) then
           !     do not allow to have negative rf power from mc
           irepeat = (ipass.le.2*ipass_num)
        end if
        if(.not.irepeat.and..not.irepeat_lastone)then
           !     take the best fitting one
           irepeat_lastone=.true.
        else
           irepeat_lastone=.false.
        end if

        ! do some post-processing on rf-related outputs
        if(orbrzv_rf_option.eq.1 .and. mcrf_on) then
           ! rf-power (watt)
           ! power conversion factor erg -> watt
           zpw = 1.0e-7_fp/(tbm2-tbm1)
           orbrzv_prftot(lsbeam) = orbrzv_prftot(lsbeam)*zpw
           orbrzv_prf(:,lsbeam) = orbrzv_prf(:,lsbeam)*zpw

#ifndef __MPI
           !#ifdef __DEBUG
           ! density profiles
           orbrzv_nrf(lsbeam,:,:) = orbrzv_nrf(lsbeam,:,:)/(tbm2-tbm1)
           !#end if
#endif

           ! rf-power match test
           if(orbrzv_toric_prftot(lsbeam) < 4.0e2_fp) then
              irepeat = .false.
              irepeat_lastone=.false.
           else
              zpw = orbrzv_prftot(lsbeam)-orbrzv_toric_prftot(lsbeam)
              zpw = abs(zpw/orbrzv_toric_prftot(lsbeam))
              if(zpw < zpwr_errmax) then
                 irepeat = .false.
                 irepeat_lastone=.false.
              else
                 if(orbrzv_prftot(lsbeam).gt.0.0_fp)then
                    !                   save for future use
                    zpw2(ipass)=zpw
                    z_frnm(ipass)=orbrzv_toric_ufrnm(lsbeam)
                 end if
              end if
           end if

           if(irepeat_lastone) then
              i_pass=0
              i_pass=minloc(zpw2)
              if(i_pass(1).gt.0.and.(zpw2(i_pass(1)).lt.epsinv)) &
                   then
                 orbrzv_toric_ufrnm(lsbeam)=z_frnm(i_pass(1))
                 irepeat=.true.
              else
                 orbrzv_toric_ufrnm(lsbeam)=0.0_fp
                 orbrzv_prftot(lsbeam)=0.0_fp
                 orbrzv_prf(:,lsbeam)=0.0_fp
                 orbrzv_nrf(lsbeam,:,:)=0.0_fp
                 irepeat=.false.
              end if
              if(i_pass(1).eq.ipass) irepeat=.false. !last one was the best fit, do not repeat
           end if


           ! in the case of repeat,
           if(irepeat) then
              ! adjust toric field renormalization factor and
              if(orbrzv_prftot(lsbeam).lt.0.0_fp)then
                 zfrnm = 1.5_fp
              else
                 zfrnm =sqrt(orbrzv_toric_prftot(lsbeam)/ &
                      max(1.0_fp,orbrzv_prftot(lsbeam)))
              end if
              !     scale up is tamed to avoid too sensitive change
              !     mg               if(zfrnm .gt. 1.0_fp) zfrnm = 0.5_fp*(1.0_fp + zfrnm)

              if(.not.irepeat_lastone) then
                 orbrzv_toric_ufrnm(lsbeam) = zfrnm* &
                      orbrzv_toric_ufrnm(lsbeam)
                 !else already defined
              end if

#ifndef __MPI
              !#ifdef __DEBUG
              ! --------------------------------------------------------
              ! write rf-debug information
              do i = 1, mj
                 rfdbg_pdiag(2:6,i) = rfdbg_pdiag(2:6,i)/ &
                      (rfdbg_pdiag(1,i)+epslon)
                 rfdbg_pdiag(7,i) = rfdbg_pdiag(7,i)/(tbm2 - tbm1)
                 rfdbg_pdiag(8:9,i) = rfdbg_pdiag(8:9,i)/ &
                      (rfdbg_pdiag(1,i)+epslon)
              end do

              write(unit=orbrzv_out,fmt=*) '# ipass :', ipass
              write(unit=orbrzv_out,fmt=*) '# ion :', lsbeam, &
                   orbrzv_prftot(lsbeam), orbrzv_toric_prftot(lsbeam)
              write(unit=orbrzv_out,fmt=*) &
                   '# frnm :',orbrzv_toric_ufrnm(lsbeam)
              do i = 1, mj
                 write(unit=orbrzv_out,fmt=*) i, &
                      orbrzv_prf(i,lsbeam), &
                      rfdbg_pdiag(1:7,i), &
                      orbrzv_nrf(1,1,i), &
                      orbrzv_nrf(1,2,i), &
                      orbrzv_nrf(1,3,i), &
                      rfdbg_pdiag(8,i), rfdbg_pdiag(9,i)
              end do
              orbrzv_nrf  = 0.0_fp
              rfdbg_pdiag = 0.0_fp
              write(unit=orbrzv_out,fmt=*) ' '
              ! --------------------------------------------------------
              !#end if
#endif

           end if
        else
           irepeat = .false.
        end if
        if(orbrzv_toric_prftot(lsbeam).ge.4.0e2_fp) then
           write(nonlin,*)' orball: specie #',lsbeam, &
                ', rf repeat = ',irepeat
           write(nonlin,*)'         rf power from mc    =', &
                orbrzv_prftot(lsbeam)
           write(nonlin,*)'         rf power from toric =', &
                orbrzv_toric_prftot(lsbeam)
           if(irepeat) then
              write(nonlin,*)'      electric field rnm =', &
                   orbrzv_toric_ufrnm(lsbeam)
           end if
        end if

        i_repeat=0
        if(irepeat) i_repeat=1
        if(ntest_max.eq.-10.or.ntest_max.eq.-1) i_repeat=0
        if(ntest_max.eq.-10) then
           write(nonlin,*) &
                ' do  not apply renorm toric electric field for realistic fi'
        end if
        if(ntest_max.eq.-1) then
           write(nonlin,*) &
                ' do  not renorm apply toric electric field for maxwell fi'
        end if

     end if !myidd=0

#ifdef __MPI
     call mpi_bcast(i_repeat, 1, MPI_INTEGER, 0, nbi_comm, mpierr)
     call mpi_bcast(orbrzv_toric_ufrnm(lsbeam), 1, MPI_REAL8,0, &
          nbi_comm, mpierr)
     irepeat=(i_repeat.eq.1)
#endif
     if(abs(ntest_max).eq.1.or.ntest_max.eq.-10) then
        irepeat=.false.        !do not repeat
        mcrf_on= mcrf_on_save !restore input on rf kick
        if(ntest_max.eq.-10.or.ntest_max.eq.-1) then
           orbrzv_toric_ufrnm(lsbeam)=1.0_fp
        end if
        if(ntest_max.eq.1.and.ipass.eq.1.and.mcrf_on) then
           irepeat=.true.      !now repeat with rf kick
        end if
     end if
     if(irepeat) then
        ! zero out prf variables
        orbrzv_prftot(lsbeam) = 0.0_fp
        orbrzv_prf(:,lsbeam)  = 0.0_fp
     end if

     if(irepeat.or.(ipass.gt.1)) then
        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
           write(nonlin,*) &
                ' specie  xi             th'// &
                '             v              vpll/v   "cpu0 last ion":'

           if(iorb.eq.0) then
              write(nonlin,*) ' (species ',isb,'; cpu 0: no orbits).'

           else
              write(nonlin,8001) isb,xion,th,vion,xksid

           end if
        end if
     end if
     !
     if(irepeat) then
        call infolog(' orball: orball-repeat irepeat=t, ipass=',ipass)
        if(myidd.eq.0) then
           write(nonlin,*) &
                ' %orball: orball-repeat irepeat=t, ipass=',ipass
        end if
        goto 50
     end if




     !******************************************************
     !    mpi buffer reduction and particle collection ends
     !******************************************************


     !------------------------------------------------
     !  serial code...

     iorb = xnorbs(lsbeam) + 0.5_fp
     norbtot = norbtot + iorb

#ifdef __MPI
     call infolog(' in processor ',myidd)
     call infolog(' %mc ptcl counts, specie #',lsbeam)
     call infolog(' orig=',inuse(lsbeam))
     call infolog(' spl=',jsplit)
     call infolog(' tot=',int(puseyp(lsbeam)))
#endif
     if(myidd.eq.0) then
        write(nout,8303) myidd,lsbeam,inuse(lsbeam), &
             j_spl,int(puseyp(lsbeam))

8303    format('in processor ,',i5,' %mc ptcl counts, specie #',i1, &
             ': orig=',i10,' spl=',i10,' tot=',i10)

     end if
     !
     !%        call anmode
     !
     !  ptcl balance - related outputs
     !
     zwght2=bntots(lsbeam)
     sfdtbmc(lsbeam)=(zwght2-zwght1)/delt

     !
#ifdef __DEBUG
     call infolog('orball:  dep sums: ',sfdepmc(lsbeam))
     call infolog('         zwdep/delt= ',zwdep/delt)
     call infolog('orball:  dep rr sums: ', &
          sfdeprr(lsbeam)-sforbrr(lsbeam))
     call infolog('    (zwold-zwght1)/delt',(zwold-zwght1)/delt)
     if(myidd.eq.0) then
        write(nonlin,*) 'orball: dep sums: ', &
             sfdepmc(lsbeam),zwdep/delt

        write(nonlin,*) ' orball: dep rr sums: ', &
             sfdeprr(lsbeam)-sforbrr(lsbeam),(zwold-zwght1)/delt

     end if
#endif
     !
     !
     ! ncx0_orb # of cx tracks collected during orbiting.
     ! may not be enough for desire statistics.
     ! check if it needs to be improved cx statistics
     !
     xncx0_use=real(ncx0,fp)
     ! ifdef __MPI, in nubeam_buf_reduce  xncx0_orb(lsbeam) could be redefine to 0 if
     ! not enough cx stats

     if((xncx0_orb(lsbeam).gt.0).and. &
          xncx0_use.gt.xncx0_orb(lsbeam)) then
        !
        call cx_tracks_need
        if(ncx1.gt.0) then  !ncx1 #cx tracks per cpu

           call cxnrm       ! make a renormalization of the nutrav outputs
           call depall_mpi_split(1,ncx1*numprocd,nuseycx)
           !     parrallel part
           !     allocate random sample arrays;
           !     get sorted list of real(fp) :: random numbers
           indcx=lsbeam*123100+myidd*1.0e+6_fp  
           call nbi_random_init_pll(iseed,indcx,iscreen1)
           call nbi_fbm_random_fetch(nuseycx,ierr)
           call nbi_distrsample(nuseycx,ierr)
           !
           ! nuseycx # cx tracks on each cpu to be gathered
           if(myidd.eq.0) then
              write(nonlin,*) ' orball need ',ncx1*numprocd, &
                   ' cx tracks'
           end if
           call infolog(' orball -- start cx tracks gathering')
           do icx=1,nuseycx
              !     find weight here
              call beamcx0(icx,ierr)
              if(ierr.gt.0) then
                 call errmsg_exit &
                      (' ?orball: beamcx0 failed')
              end if

           end do
           call infolog(' orball -- done cx tracks gathering')
        end if
     end if
#ifdef __MPI
     ! reduce all data assotiated with impoved cx statistics
     call nubeam_cxbuf_reduce

#endif

     call outptb(iorb,ierr)
     if(ierr.gt.0) then
        call errmsg_exit &
             (' ?orball: outptb failed')
     end if

     orb_code_error(lsbeam) = orb_code_error(lsbeam) + ierr
     !
     call orbsmoo                  ! radial smoothing
     !
     if(nbsplt.gt.0) then
        call errorlog('orball: error nbsplt gt.0 at end')
        call errorlog('of timestep, split particles lost ****')
        call errorlog('nbsplt= ',nbsplt)
        if(myidd.eq.0) then
           write(nout,7999) nbsplt
        end if

7999    format(' ??? error in orball: nbsplt=',i10,', .gt.0 at end'/ &
             '     of timestep, split particles lost ****')
        ierr=1
        orb_code_error(lsbeam) = orb_code_error(lsbeam) + ierr
        if(ierr.gt.0) then
           call errmsg_exit &
                (' ?orball: split particles lost')
        end if

     end if
     !
     !  end of beam species loop
     !
     idum = orb_code_error(lsbeam) + 0.5_fp
     iersum = iersum + idum
     if(iersum.eq.0) then

     end if

     mcrf_on=mcrf_on_save   !restore input value

  end do !end of species loop

  if(any(ilaep_init)) then       !clean up aep or RBQ
     call nbi_aep_clean
  end if


  !     deallocate arrays for splitting
  call nbi_dalloc_splt


#ifdef __MPI
  if(nlbout) then
     call infolog('orball: lost ptcl. gather')
     if(myidd.eq.0) then
        write(nonlin,*) ' %orball: lost ptcl. gather'
     end if
     call orball_lost_gather
  else
     call infolog('orball: lost ptcl.  not gathered')
  end if
#endif



  !
  nbs=0
  lsbeam=0
  abeam=0.0_fp
  xzbeami=0.0_fp
  !
  !  update dt_adj0s(...)
  !
  do isb=1,nsbeam
     if(nlfprod(isb)) then
        znlimdt=nptclf
     else if(nlheavy(isb)) then
        znlimdt=nptclh
     else
        znlimdt=nptcls
     end if

     znlimdt = max(10.0_fp,(0.01_fp*znlimdt)) - 0.001_fp
     if(xnum_dt_adj0(isb).gt.znlimdt) then
        zadjust = dt_adj0_sums(isb)/xnum_dt_adj0(isb)
        dt_adj0s(isb) = dt_adj0s(isb)*zadjust
     end if

     if(isb.eq.1) then
        dt_adj0=dt_adj0s(isb)
     else
        dt_adj0=min(dt_adj0,dt_adj0s(isb))
     end if
  end do
  !
  flag_dt_reset = ifirst    ! global dt adjust flag for next call...
  !
  !  ein_ien calculation renorm
  !
  do ib=1,nbeam
     do ien=1,3
        zw=wts_ien(ien,ib)
        if(zw.eq.0.0_fp) then
           ein_ien(ien,ib)=0.0_fp
           dei_ien(ien,ib)=0.0_fp
        else
           zav=ein_ien(ien,ib)/zw
           ein_ien(ien,ib)=zav
           zvar2=dei_ien(ien,ib)/zw - zav*zav
           dei_ien(ien,ib)=sqrt(max(0.0_fp,zvar2)) ! rms variance
        end if
     end do
  end do
  !
  !%      nbplot=ibplot
  !
#ifndef __MPI
  ! for serial jobs only:

  if(norbtot.gt.0) then
     if(myidd.eq.0) then
        write(nonlin,*) &
             ' specie  xi             th'// &
             '             v              vpll/v     "last ion":'
     end if

     do isb=1,nsbeam
        if(nuse(isb).gt.0) then
           ind=nbndex(nuse(isb),isb)
           ind=find_ptcl_ind(numprocd,myidd,ind) !find proper ptcl.id#
           xion=xiay(ind,isb)
           th=thay(ind,isb)
           vion=vay(ind,isb)
           xksid=xksidy(ind,isb)

           if(myidd.eq.0) then
              write(nonlin,8001) isb,xion,th,vion,xksid
           end if


        end if
     end do
  end if
#endif

8001 format(3x,i1,3x,4(1x,1pe14.7))

  nmcloss = nmcloss + xmcloss_o + 0.5_fp

  !      moved to depall.for 09/04/09 mg
  !      call nmcsum

  iorbtot = norbtot
  call infolog(' orball: # of orbit = ',iorbtot)


  call cpu_time(ztim2)

  ierr = iersum

  !-------------------------------------------
  !  clear sawtooth flags -- must be reset to get another sawtooth event

  sawflag = .false.
  sawflag2 = .false.

  !----------------------------------------------------------------
  !   print execution times for various components
  !   after all species loops (mpi related)
  !----------------------------------------------------------------

#ifdef __MPI
  ttime32=mpi_wtime()
  if (myidd.eq.0) cpborb = cpborb + (ttime32-ttime31)
#else
  cpborb = cpborb + (ztim2-ztim1)
#endif

#ifdef __MPI

  if (impitime.gt.0) then

     call warnlog(' ')
     call warnlog('in nbstart: myidd= ',myidd)
     call warnlog('total time in nbstart= ',ttime12-ttime11)
     call warnlog('total time for file read= ',stime2-stime1)
     call warnlog(' ')
     call warnlog('in depall: myidd= ',myidd)
     ztime1=ttime22-ttime21
     ztime2=stime4+stime5+stime6+stime7+stime8+stime10
     call warnlog('total parallel processing time= ',ztime2)
     call warnlog('total time for calculation= ',ztime1-ztime2)
     call warnlog('total time in depall/add_track= ',stime7)
     call warnlog('time for mpi dpt_gather= ',stime10)
     call warnlog('total time in depall= ',ztime1)
     call warnlog(' ')
     call warnlog('in orball: myidd= ',myidd)
     ztime3=(ttime32-ttime31)-stime9
     call warnlog('time for mpi buf reduce= ',stime9)
     call warnlog('time for calculation= ',ztime3)
     call warnlog('total time in orball= ',ttime32-ttime31)
     call warnlog('total time in orball (using cptim)= ', &
          ztim2-ztim1)
     call warnlog(' ')
     ztime4=(ttime12-ttime11)+(ttime22-ttime21)+&
          (ttime32-ttime31)
     call warnlog('in nbstart, depall, orball: myidd= ',myidd)
     call warnlog('total time = ',ztime4)
     call warnlog(' ')
  end if

#endif

  goto 9999

  ! -- allocation error --
  ierr=777
  call errorlog('orball: allocation failure,too many particles?')
  if(myidd.eq.0) then
     write(nonlin,'(a)') '?orball: allocation failure, are ' &
          //'you using too many particles?'
  end if

9999 continue

  if(itmp_file) then
     open(unit=lunres,file=zorball_tmp_filename,status='old')
     close(unit=lunres,status='delete')
  end if
  call infolog(' exiting orball')
  return
  
contains

  subroutine cx_tracks_need

 ! find out how many tracks (ncx1) need to be evaluted on each cpu to improve cx statistics
 ! ncx0 can be adjusted here to be a multiple of numprocd = ncx0_use
    integer :: ncx_rs
    real(fp) :: zrem

    if(xncx0_orb(lsbeam).lt.1._fp) then
       ncx1 = 0
       return
    end if
    ncx_rs = int(xncx0_use - xncx0_orb(lsbeam))
    
#ifdef __MPI
 !
 !--------------------------------------------------------
 !

    if(ncx_rs>numprocd) then
       ncx1 = ncx_rs/numprocd
       if(mod(ncx_rs,numprocd).ne.0) then
          ncx1 = ncx1+1       ! round up
       end if
    else
       ncx1 = 1
    end if
    ncx0_use = int(xncx0_orb(lsbeam))+ncx1*numprocd

#else
    ncx1 = ncx_rs
    ncx0_use = int(xncx0_use)
#endif

  end subroutine cx_tracks_need

  subroutine orball_loclist_write
    character(len=6) :: zcpid
    integer :: inum, iost

 ! write ptcl arrays in tmp binary file -- do not include
 ! sawtooth variables; sawtooth operator already executed
 ! if needed.

    write(zcpid,'(i5)') myidd
    zorball_tmp_filename = trim(workpath)// &
         trim(runid)//'_orball_cpu'//trim(adjustl(zcpid))//'.tmp'
    
    open(unit=lunres,file=zorball_tmp_filename,status='unknown', &
         access='sequential',form='unformatted',iostat=iost)
    if(iost > 0) then
       call errorlog(" ?orball_loclist_write: file open failed"// &
            trim(zorball_tmp_filename))
       
       call errorlog(' emergency exiting orball_loclist_write')
       call bad_exit("orball_loclist_write : file open failed")
    end if
    
    inum = ntotptp(lsbeam)
    write(lunres) inum,nuse(lsbeam)

    write(lunres) puseyp(lsbeam),nuse_total(lsbeam)
    
    write(lunres) nbndex(1:inum,lsbeam)
    
    write(lunres) nbscay(1:inum,lsbeam)
    write(lunres) nbienay(1:inum,lsbeam)
    write(lunres) njavay(1:inum,lsbeam)
    
    write(lunres) cxpray(1:inum,lsbeam)
    write(lunres) fppray(1:inum,lsbeam)
    write(lunres) goosay(1:inum,lsbeam)
    write(lunres) einjay(1:inum,lsbeam)
    write(lunres) pmechay(1:inum,lsbeam)
    write(lunres) pinjay(1:inum,lsbeam)
    write(lunres) pphiay(1:inum,lsbeam)
    write(lunres) rmjionay(1:inum,lsbeam)
    write(lunres) tdepay(1:inum,lsbeam)
    write(lunres) tfluxay(1:inum,lsbeam)
    write(lunres) thay(1:inum,lsbeam)
    write(lunres) vay(1:inum,lsbeam)
    write(lunres) wghtay(1:inum,lsbeam)
    write(lunres) xiay(1:inum,lsbeam)
    write(lunres) xistay(1:inum,lsbeam)
    write(lunres) xksidy(1:inum,lsbeam)
    write(lunres) xmuay(1:inum,lsbeam)
    write(lunres) xzbmay(1:inum,lsbeam)
    write(lunres) xzionay(1:inum,lsbeam)
    
    close(unit=lunres)
  end subroutine orball_loclist_write

  subroutine orball_loclist_read
    integer :: inum, iost
    real(fp), dimension(:), allocatable :: dummy
    
    open(unit=lunres,file=zorball_tmp_filename,status='old', &
         access='sequential',form='unformatted',iostat=iost)
    if(iost > 0) then
       print *, " ?orball_loclist_read : file open failed"
       call errorlog(' orball_loclist_read : file open failed')
       call bad_exit("orball_loclist_read : file open failed")
    end if

    read(lunres) inum,nuse(lsbeam)
    ntotptp(lsbeam) = inum
    if(ntest_max.eq.1) then !update only weight, goose
       if(allocated(dummy)) deallocate(dummy)
       allocate(dummy(inum))
       dummy=0.0_fp
    end if

    read(lunres) puseyp(lsbeam),nuse_total(lsbeam)
    read(lunres) nbndex(1:inum,lsbeam)
    read(lunres) nbscay(1:inum,lsbeam)
    read(lunres) nbienay(1:inum,lsbeam)
    read(lunres) njavay(1:inum,lsbeam)
    if(ntest_max.eq.1) then !update only weight, goose
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
    else
       read(lunres) cxpray(1:inum,lsbeam)
       read(lunres) fppray(1:inum,lsbeam)
    end if
    read(lunres) goosay(1:inum,lsbeam)
    if(ntest_max.eq.1) then !update only weight, goose
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
    else
       read(lunres) einjay(1:inum,lsbeam)
       read(lunres) pmechay(1:inum,lsbeam)
       read(lunres) pinjay(1:inum,lsbeam)
       read(lunres) pphiay(1:inum,lsbeam)
       read(lunres) rmjionay(1:inum,lsbeam)
       read(lunres) tdepay(1:inum,lsbeam)
       read(lunres) tfluxay(1:inum,lsbeam)
       read(lunres) thay(1:inum,lsbeam)
       read(lunres) vay(1:inum,lsbeam)
    end if
    read(lunres) wghtay(1:inum,lsbeam)
    if(ntest_max.eq.1) then !update only weight, goose
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
       read(lunres) dummy(1:inum)
    else
       read(lunres) xiay(1:inum,lsbeam)
       read(lunres) xistay(1:inum,lsbeam)
       read(lunres) xksidy(1:inum,lsbeam)
       read(lunres) xmuay(1:inum,lsbeam)
       read(lunres) xzbmay(1:inum,lsbeam)
       read(lunres) xzionay(1:inum,lsbeam)
    end if
    
    close(unit=lunres)
    if(ntest_max.eq.1) then !update only weight, goose
       if(allocated(dummy)) deallocate(dummy)
    end if
  end subroutine orball_loclist_read
end subroutine orball
!=====================================================================
#ifdef __DEBUG
subroutine zbal_debug(iorb,zwdep,zwold)
  use nbi_com
  use nbi_mpi_com
  use logmod
  implicit none
  
  integer, intent(in) :: iorb
  real(fp), intent(in) :: zwdep,zwold
  
  !  debug: compute interim ptcl balance
  
  real(fp) :: zbal,zmax,bsumint,zsumth
  
  integer, parameter :: ifreq = 100
  integer, parameter :: imin = 100
  
  integer :: itest
  
  !-------------------------------
  
  
  if(iorb.lt.imin) return
  itest = iorb/ifreq
  itest = itest*ifreq
  if(itest.ne.iorb) return
  
  sfdtbmc(lsbeam)=bntots(lsbeam)-zwold
  zsumth=bsumint(sorbths(1:mj,lsbeam))
  
  zbal = zwdep &
       -bxflws(lsbeam)-zsumth &
       -bxfshs(lsbeam) &
       -sfcxesc(lsbeam)-sfcxrr(lsbeam)-bxrpls(lsbeam) &
       -sfdtbmc(lsbeam)
  zbal=abs(zbal)
  zmax=max(1.0_fp,abs(zwdep))
  zmax=max(zmax,abs(bxflws(lsbeam)))
  zmax=max(zmax,abs(zsumth))
  zmax=max(zmax,abs(bxfshs(lsbeam)))
  zmax=max(zmax,abs(sfcxesc(lsbeam)))
  zmax=max(zmax,abs(sfcxrr(lsbeam)))
  zmax=max(zmax,abs(bxrpls(lsbeam)))
  zmax=max(zmax,abs(sfdtbmc(lsbeam)))
  
  call infolog('orball: iorb=',iorb)
  call infolog('orball: ratio = ',zbal/zmax)
  if(myidd.eq.0) then
     write(nonlin,*) ' iorb=',iorb,' ratio = ',zbal/zmax
  end if
  return
end subroutine zbal_debug
#endif
