!
!      this subroutine advances an ion in its orbit,
!      and keeps track of power deposition
!
!
subroutine orbit(iorb,ierr)
  !
  !  dmc oct 2003 -- cleanup; reorganization of goosing and timestep control
  !       making some use of f90 techniques.  old "fbtest" calculations
  !       deleted (5% speedup).
  !
  !  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)...
  !
  !  rga oct00 - removed levgb = 1,2
  !
  !  jec  12-may-99-- james.conboy@jet.uk
  !                   #ifdef __f90, use standard open statements
  !
  !  dmc 27 sep 94-- support for radial current calculation
  !                  set xistop for orbit loss events; for
  !                  thermalization and not-lost orbits, orball
  !                  will compute the orbit average position with
  !                  xjadep
  !
  !  mhr  july 94 -- ripple model coding added
  !
  !  jfs 24 mar 94-- added tests for trouble with bad orbits (stmt 800)
  !
  !  dmc mar 1994 -- goose=1.0000 is reserved for newly deposited ions
  !    use goose=1.0001 as minimal goosing for older ions.
  !
  !  tbt sep 1993 -- put in goto 850 between call collid & abdifu.
  !
  !  dmc nov 1992 -- upgrade to support heavy beam ions.  these ions
  !    are stripped down one electron at a time -- cf subroutine beamiz
  !
  !  dmc mar 2010 -- call to bbdifu added -- *ae mode radial transport operator
  !    mu & etot conserved; momentum not conserved.
  !
  !  dmc sep 1990 -- call to abdifu added -- anomalous radial diffusion
  !  of beam ions.  goosing & mc speedup numerics handled as if abdifu
  !  were part of collid (which is where the code used to be before
  !  generalization of the code for rotating plasma)
  !
  !        bitime(1)-- orbit time of current ion. (now common)
  !          was passed variable ztim
  !          at beginning of timestep, ztim=0.0; at end ztim=delt
  !          initially injected ions start with randomly spaced
  !          times in the range of 0 to dtbm
  !          (this upgrade needed for fishbone simulation
  !          rob goldston  6-10-82)
  !
  !        xjav (now common, was passed as zjav)
  !       -- avg location (index) of ion over its last complete orbit
  !          (defined as having crossed midplane twice)
  !          if ion orbit not completed this timestep, avg j over partial
  !          orbit is used
  !
  !      code starting with c% is used only for plotting purposes
  !
  !      rev. t.terpstra november, 1991
  !      added check of nlfatom.
  !
  !        rev. d. mc cune july 14 1980
  !          cleanup: remove multiple if statements
  !          bug fix-- reinstate suppressed "swsq"
  !
  !        revised goose-- midplane test: sin(theta) is compared
  !          previous sin(theta) in
  !          section; if sign has changed, midplane has been crossed
  !          (test changed because goose section is in range of
  !          fpskip collision-skip section and old test might have
  !          missed some midplane crossings --- d. mc cune oct 1980
  !
  !      revised to include finite larmor radius effects,
  !      and to extract collisions into a separate subroutine
  !      12/15/82      rob goldston
  !
  !      upgrade to allow ions to travel outside the plasma, but
  !      inside a limiter system, defined in limlos.for.
  !      at the same time we can accumulate the beam ion distribution
  !      function out to a minor radius of rminbm, rather than just
  !      rminor, as previously.
  !      rob goldston july 21, 1983
  !
  !---------------------------------------------------------------
  !  geometry generalization - dmc, b**2 ... 1985
  !
  !   orbit calculation broken up into subsidiary routines for use
  !   in varying geometries far from or close to the magnetic axis:
  !
  !   these 3 are bob's routines:
  !    orbc - circular geometry orbit solver far from axis - no radial e-field
  !    orbs - shifted circular geometry orbit solver with radial electric field
  !    orbg - generalized geometry orbit solver far from axis with everything
  !
  !    orbpre - set up common variables for bob's routines
  !    orbout - copy outputs of bob's routines to variables used
  !             in rest of beam code.
  !
  !----------
  !xx    bob's compression routines are also added:
  !
  !    dmc nov 2001:  compression routines removed: replaced by
  !      flux position / mu / pphi conservation operator in orball.
  !
  !xx    cmprsg - compression operator (arb. geometry)
  !
  !xx  ** compression routines also compute xi <--> eldot terms in
  !xx     beam ion distribution **
  !
  !xx    cmppre - load common for compression operators
  !xx    cmpout - copy output of bob's compression operators
  !
  !-----------------------------------------------------------------
  !
  use nbi_com
  use physconst_mod, only : evptocm_sec
  use map_larmor_mod, only: map_larmor_loss, &
       map_larmor,map_larmor_compare
  use xstraln_calls
  use nbi_mpi_com
  use nbi_lost_ptcl, only: gtime,yyion,yze,ynbsce,ynbien !time, z, energy lost ptcl.,
  !source and energy fraction index
  use nbi_lost_ptcl, only: store_lost_ptcl
  use nbi_random_pll
  use orbrzv, only : mcrf_on, orbrzv_pflag,&
       orbrzv_nmidc,orbrzv_psiv,orbrzv_psix,orbrzv_bnw,orbrzv_rhopm,&
       orbrzv_rmid0,orbrzv_rc,orbrzv_zz,orbrzv_trpas,orbrzv_prmax,&
       orbrzv_bnw_lrf,orbrzv_rfdv,orbrzv_zrfd_min,il_rz_updated,nlmcrfgcon
  use logmod
  use nbi_aep_com, only : ilaep_init
  !%      use sglib_mod
  !
  implicit none
  integer,intent(in) :: iorb
  integer,intent(out) :: ierr
  !
  logical imap_larmor_loss ! .true. to us map_larmor_loss instead of limlos
  !
  !        ** local indexing scheme ** for efficient access in calculating
  !        the beam distribution function.  oct 1982 -- d. mc cune
  !
  logical ilfbm
  !
  !  explicit declarations for all local variables...
  !
  integer :: iplfp,iplcx,jbe,ievent,icol,iflri,itrak,istep,inew
  integer :: iswgoose,isumiz,icheck,istpor,istpcx,istpfp,inexcx,inexfp
  integer :: igoosc,ilim,ifkill,iripl,iflip,isuma,jba,ie,isw,iinc
  integer :: ier,iwarn,idum,ingcp,ijxb,ingc
  integer :: itrap, iuncon, ixiav

  real(fp) :: zpph1,zpph2,zeph1,zeph2,fnbint

  real(fp) :: zmaxe,zgoomax,ztimfp,zympl,ztesm,zjsum,zjtim,zdtime
  real(fp) :: ztscxp,ztslmp,ztslmg,zdph,ze1,ze2,zdltcx,zfpskip,zphp,ztesp
  real(fp) :: axicp,zstpor,zgoo,zfacfp,zfaccx,ztbtcx,zgooscx,zdvgef
  real(fp) :: zfdvg,ztbtsl,zgoosl,zfdv2g,zgsbde,ztbtsq,zgoosq,zgoosd
  real(fp) :: ztbdsq,zdelxi,zgoodf,zgoosx,zgtes,zgooso,zgtest,zgoosq2
  real(fp) :: zmaxsk,zgcomin,zgoost,zthchk,zxichk,zphiarg,zhgt
  real(fp) :: ziripl,zxip,zepp,zephip,zbphip,zbtorp,zksid
  real(fp) :: zdsdfb
  real(fp) :: zep,zdtg,zdummy,zdsxvb,ztmp,zdelxi2,zfacd
  real(fp) :: ztmtra,ztmfsh,zevolt,zfrate,zflive,zde,zxion_min
  real(fp) :: zthc2,zrc2,zvpvc2,zthc1,zrc1,zvpvc1,zbpksid,zepp_jxb
  real(fp) :: zrgdv, zrgdx, zgoorf, zgooconrf, zrhopm, zrsg
  real(fp) :: xjadep, zzv0, zxa, zxksav, zeiav, zphiav, zxhw
  real(fp) :: zpass

  !  orbit error stats -- per (half) bounce

  real(fp) :: zderrmin,zderrmax,zderrsum
  real(fp) :: zdprrmin,zdprrmax,zdprrsum
  real(fp) :: zadjust
  logical :: icheck_dt
  !
  !  debugging -- extra output to track trajectory of chosen particle
  !  (set to zero normally, to suppress debug output)
  !
  integer :: isb_trace = 0 !  1  ---- the fast specie index
  integer :: iorb_trace = 0 !  14041  ---- the selected orbit (iorb value)
  integer :: itrak_trace = 0 ! 1 -- the selected track (itrak value)
  real(fp) :: zesave
  !
  !  min and max b seen by ptcl on current orbit (since last goose eval)
  !
  real(fp) :: zminbo
  real(fp) :: zmaxbo
  !
  !  min,max xion seen by ptcl on current orbit (since last goose eval)
  !
  real(fp) :: zmaxxo,zminxo
  !
  !  info for <e>, <vpll> half-bounce averages
  !
  real(fp) :: zdtsum,zvpll_av
  !
  !  info saved or used at orbit bmax
  !
  real(fp) :: zx_bmax,zmu_bmax,zbrefl_bmax,zvpll,zomeg
  real(fp) :: zvpp  ! helper var.
  !
  !  guiding center b on previous timestep
  !
  real(fp) :: zmodb_prev
  !
  !  guiding center b on previous midplane crossing, and the one before that
  !
  real(fp) :: zmodb_midp
  real(fp) :: zmodb_midpp
  real(fp) :: zmodb,zratb          ! work variables
  real(fp) :: zbmax_fs,zbmin_fs    ! min & max b on flux surface
  real(fp) :: zxi,zxfac,bitime_save
  integer :: igcmax
  !
  !  v*dt sum -- compare to vtaub
  !
  real(fp) :: zvdtsum                      ! related to trajectory length, cm
  !
  integer :: imin_mpc               ! min desired midplane crossings
  !
  integer :: itrig                     ! reason for goose update
  !
  !  itrig=1 -- normal reason (midplane crossing)
  !  itrig=2 -- reason: bounce time expired without midplane crossing
  !
  integer :: ndxfbe_r8
  real(fp) :: alpha_v ! poloidal aiming angle
  real(fp) :: rmjion_prev,yion_prev
  !
  !%      integer :: ieplot,ibdbg,ipflr,iplsty
  !%      integer :: iplct,iplnu,jt,iregion,ioops
  !%      real(fp) :: zrelmin,zrmj
  !%      real(fp) :: zderel1,z_vperp,z_vpar
  !%      character*80 zbuf(5)
  !%      real(fp) :: zzx(2),zzy(2),zzx2(3),zzy2(3)
  !
  !--------------------------local temporary fix-rtm-----------
  !
  !%      integer, parameter :: izplot=20000
  !%      real(fp) ::  zxpos(izplot)
  !%      real(fp) ::  zypos(izplot)
  !c%      integer ::          ndvopt      !0=>descriminant positive (zdvopt.fnb)
  !                               !1=>negative descriminant
  !                               !2=>rounding error d(vpar)/dt used
  !
  !%      data  zrelmin/1.0e-6/ !min relative error in energy for plotting
  !%      data  ieplot/6/      !symbol for energy non-conservation points
  !
  !%      data ibdbg/0/   !set to 1 for mccann debug symbols
  !%      data ipflr/0/   !set to 1 or more for flr symbols
  !---------------------------------
  !
  iplfp=0                           ! plot/debugging support
  iplcx=0                           ! ditto
  !
  !        set-up single pt plotting: fp, ch-ex locations
  !%      zzx=0; zzy=0; zzx2=0; zzy2=0; zxpos=0; zypos=0
  !%      zzx(1)=1.0
  !%      zzy(1)=1.0
  !%      call stplot
  !%     &  (1,iplnu,iplct,iplcx,iplfp,iplsty,ipflr,
  !%     &   bitime(1))
  !%      jt=1
  !
  !
  !**************  initialization  **********************
  ierr=0
  !
  jbe=0                             ! init e grid search position
  imap_larmor_loss = nlbflr .and. nlbgflr .and. gflr_op(4)==0._fp
  !
  zgoomax=epsinv                    ! goomax disabled dmc sept 1990
  !
  icol=0                            ! count collid calls (debugging)
  itrak=0                           ! count tracks (debugging)
  istep=0
  il_rz_updated=.false.

  !
  !  debug energy jump errors...
  !
  zesave = v2toev*vion*vion         ! energy on entry to orbit
  zmaxe = maxval(efbmb(:,lsbeam))      ! max allowed energy
  zmaxe = max(zmaxe,edbfac*einjs(lsbeam) )

  !
  !  set flag for stepping around the gyro explicitly -- heavy
  !  beam ions only, first orbit call only... beamiz will also turn this
  !  off after first few electrons are stripped from the heavy beam ion...
  !
  if(nlbflr.and.nlheavy(lsbeam).and.(bitime(1).gt.0.0_fp)) then
     iflri=1
  else
     iflri=0
  end if
  !
  if(nlheavy(lsbeam)) then
     ilfbm=.false.
  else
     ilfbm=nlfbm
  end if
  !
  ! * re-entry point for new orbit after charge-exchange:
  !
5 continue
  itrak=itrak+1
  istep=0
  !
  if(iorb_trace.eq.iorb) then
     if((lsbeam.eq.isb_trace).and.(itrak.eq.itrak_trace)) then
        call infolog('orbit: orbit_trace iorb = ',iorb)
        call infolog(' lsbeam  = ',isb_trace)
        call infolog(' itrak = ',itrak_trace)
        if(myidd.eq.0) then
           write(nonlin,*) ' orbit_trace: iorb = ',iorb, &
                ' lsbeam,itrak = ',isb_trace,itrak_trace
        end if
     end if
  end if
  !
  icheck_dt=.false.
  if(goose.lt.1.000001_fp) then
     inew=1
     if(orbrzv_option .eq. 1) then
        dt_acc=orbrzv_zzerr_con
     end if
     if(dt_acc.gt.0.0_fp) then
        dt_acc=max(1.0e-6_fp,min(1.0e-3_fp,dt_acc))
        if(orbrzv_option .eq. 1) then
           orbrzv_zzerr_con= dt_acc
        end if
        icheck_dt=.true.
     end if
  end if
  !
  xicur=epsinv
  !
  !  the field strength has not yet been evaluated ***
  bimodb(1)=0.0_fp
  zmodb_midp = 0.0_fp
  zmodb_midpp = 0.0_fp
  !
  dtrap = -1.0_fp  ! assume passing for now...
  !
  !  the orbit timestep is not yet known
  bidt(1)=0.0_fp
  !
  !  auxilliary gc parameters and zone indices; new orbit
  !  check location...
  !
  !
  if(abs(goose-1.000000_fp).lt.1.0e-6_fp) then
     pbdepmc(lsbeam)=pbdepmc(lsbeam)+wght*vion*vion*v2toev
  end if
  call nbigc(1)
  if(ngc.eq.0) then
     icheck=min(1,istep)    ! 0 on first step-- instant orbit loss
     ilim=1
     zepp_jxb=vion*vion*v2toev
     ijxb=1
     goto 800
  end if
  !
  !  re-init goose counters
  !
  igoosc=0                          ! used to check if goose gets updated.
  ngoose=0                          ! midplane crossing counter
  !
  !      for fishbone simulation need to get two midplane crossings
  !      promptly (in tokamak time) at start of beam time step
  !         there may be other reasons for doing this; nlgoose_reset
  !         set in nbstart.for
  !
  if(nlgoose_reset) goose=min(1.0001_fp,goose)
  !
  !  clear particle status flags
  !
  nbflg=0
  ilim=0
  ifkill=0
  iripl=0
  nbicx=0
  !
  istep=0                           ! orbit step counter
  ievent=0                          ! cx flag
  iswgoose=0
  !
  zvpvc2=0.0_fp
  zrc2=0.0_fp
  zthc2=0.0_fp
  !
  if(nmdifb.eq.4.or.nmdifb.eq.6) then
     call bbdifu_init(ierr)
     if(ierr.ne.0) return
  end if
  !
  !  the above are initialized at the start of each orbit track...
  !
  !  set up goose, fp and cx acceleration schemes
  !
  !  use the contained subroutine for things that are re-initialized
  !  each half-bounce (or other operative goose re-evaluation period).
  !
  call bounce_ini(1)
  !
  !  check for ptcl past end of timestep (could result from cx splitting)
  !
  if(bitime(1).ge.delt) then
     nbflg=1
     goto 890
  end if
  !
  !****************   top of main orbit loop ******************
  !


100 continue
  istep=istep+1
  !
  !      if(jt.lt.minbplt) then      !rtm may 1990 - avoid subscript overrun
  !%      if(jt.lt.izplot) then      !rtm may 1990 - temp fix - bigger arrays
  !%        jt=jt+1
  !%      if(nbplot.gt.0) then
  !%          if(nbplot.eq.1) then
  !%          zxpos(jt)=rmjion
  !%          zypos(jt)=yion
  !%          else
  !%            zxpos(jt)=vpll
  !%            zypos(jt)=vperp
  !%          end if
  !%!turn on symbols if energy is not conserved
  !%          if(bim(1).eq.0.0) then
  !%            zderel1=0.0
  !%          else
  !%            zderel1 = abs(bidte(1))/(0.5*bim(1)*biv(1)**2)  !ignore e*phi
  !%          end if
  !%          if(zderel1.gt.zrelmin) then      !plot a symbol
  !%           if(ibdbg.gt.0) then
  !%            zzx(2) = zxpos(jt)
  !%            zzy(2) = zypos(jt)
  !%            call symbl(ieplot)
  !%            call line(-1)
  !%            call r8_cplot(zzx,zzy)
  !%           end if
  !%          else
  !%          end if
  !c%          if(ndvopt.gt.0) then      !plot a symbol
  !c%           if(ibdbg.gt.0) then
  !c%            zzx(2) = zxpos(jt)
  !c%            zzy(2) = zypos(jt)
  !c%            call symbl(ieplot+ndvopt)      !ndvopt=1,2 (neg.descr; roundoff)
  !c%            call line(-1)
  !c%            call r8_cplot(zzx,zzy)
  !c%          end if
  !c%          else
  !c%          end if
  !%!
  !%        end if
  !%      end if
  !
  !
  !      ---------------------------------------------
  !      2.0      choose which orbit subroutine to call
  !      ---------------------------------------------
  !
  !
  !
  !   ******
  !   define the calling args used by orbc, orbs, and orbg
  !   ******
  !
  zmodb_prev = bimodb(1)
  !
  ingcp=ngc
  zxip=xion
  zepp=v2toev*ep
  zepp_jxb=zepp
  zephip=ephion
  zbphip=bphion
  ijxb=0                    ! set =1 when orbit jxb evaluated
  zbtorp=btori
  zksid=xksid
  !
  !       store ion parameters before ts
  !
  call orbpre
  !
  if(orbrzv_option .eq. 1) then
     isw=0
     orbrzv_pflag=0
     call orbrzg(iorb, istep, ierr)
     il_rz_updated=.false. !this flag will be set true if r,z changed bellow this line
     if(orbrzv_pflag .ne. 0) goto 800 !bad orbit
#ifndef __MPI
#ifdef __DEBUG
     if(istep == 1 .and. mod(iorb,100) == 0) then
        print *, "iorb = ", lsbeam, nsbeam, norbtot, iorb, wght
     end if
#endif
#endif
  else
     isw=0
     call orbg(ierr)        !with everything
  end if
  if(ierr.ne.0) return
  !
  !   *****
  !   store the rsults common to orbc, orbs, and orbg
  !   *****
  !
  call orbout
  ! mg 07/25/2016 - first check if lost, as spline for g, b could be weird
  ! if far beyond boundary
  !
  !
  !  check for orbit loss - beyond limits of coordinate system
  !
  if(xion.ge.xbmbnd) then
     ilim=1
     icheck=1
     goto 800
  end if
  !
  zxion_min=min(xion,zxion_min)
  !
  if(iorb_trace.eq.iorb) then
     if(itrak_trace.eq.itrak) then
        if(isb_trace.eq.lsbeam) call orbtrace
     end if
  end if
  !
  ! kjm-rfdbg
  if((v2toev*vion*vion.gt.zmaxe.and.orbrzv_rf_option.ne.1).and. &
       (v2toev*vion*vion.gt.zmaxe.and. &
       (ntest_max.eq.0.or.abs(ntest_max).eq.10))) then
     zep=v2toev*vion*vion
     call errorlog &
          ('orbit: ptcl energy (ev) detected after orbout call ', &
          zep)
     call errorlog('at start of orbit energy(ev) was: ',zesave )
     call errorlog('species #',lsbeam)
     call errorlog(' ab = ',abeam)
     call errorlog(' zb = ',xzbeami)
     call errorlog(' xion = ',xion)
     call errorlog(' vpll = ',  vpll)
     call errorlog(' vperp = ',vperp)
     call errorlog(' vion = ',vion)
     call errorlog(' rmjion = ',rmjion)
     call errorlog(' yion = ',yion)
     call errorlog(' istep = ',istep)
     call errorlog(' iorb = ',iorb)
     if(myidd.eq.0) then
        write(nonlin,9101) zep
        write(nonlin,*) 'xion,vpll,vperp,vion : ',xion,vpll,vperp, &
             vion
        write(nonlin,*) 'rmjion,yion,istep,iorb :', rmjion,yion, &
             istep,iorb
        write(nonlin,*) &
             'at start of orbit energy was ',zesave,' ev.'
        write(nonlin,*) 'species #',lsbeam,' ab,zb=',abeam,xzbeami
     end if
9101 format(/ &
          ' ?orbit -- ',1pe11.4, &
          ' ev ptcl detected after orbout call.'/)
     ! if mcrf option is on, just report the high energy
     if(orbrzv_rf_option .ne. 1) then
        if(.not.ilaep_init(lsbeam)) then
           ierr=1
           return
        end if
     end if
  end if
  !
  !  evaluate position, larmour radius of guiding center
  call nbigc(1)
  if(ngc.eq.0) then
     icheck=1
     ilim=1
     goto 800
  end if
  !
  !   jxb torque - orbit step
  !
  zep=v2toev*ep
  call nb0jxb(wght,zepp_jxb,zep,zbphip,bphion,zbtorp,btori, &
       ingcp,ngc,zxip,xion,4,.false.)
  zepp_jxb=zep
  ijxb=1
  !
  !  add to (physical) time-since-last midplane crossing
  !  update max b seen by orbit

  !dbg      write(nonlin,8888) istep, xion,th, bimodb(1)
  !dbg 8888 format(' istep,x,th = ',i6,2(1x,1pe12.5),' mod(b)=',1pe12.5)

  ztslmp=ztslmp+bidt(1)  ! ungoosed half-bounce-time
  zvdtsum = zvdtsum + vion*bidt(1)
  !
  ! frame correction to take into account rotation/electrostatic trapping
  ! effect, approximately...
  zomeg = omegag(ngc)+xingc*(omegag(ngc+1)-omegag(ngc))
  zvpll = vpll - rmjion*zomeg/btiltc(2,ngc)
  !
  if(bimodb(1).gt.zmaxbo) then
     zmaxbo = bimodb(1)
     zx_bmax= xion
     zvpp = max(1.0_fp,vperp)
     zmu_bmax = (zvpp*zvpp)/zmaxbo
     zbrefl_bmax = max(zmaxbo,(zvpp*zvpp+zvpll*zvpll)/zmu_bmax)
  end if
  zminbo = min(zminbo,bimodb(1))
  zmaxxo = max(zmaxxo,xion)
  zminxo = min(zminxo,xion)
  !
  zdtg=bidt(1)*goose
  !
  zdtsum = zdtsum + zdtg
  zvpll_av = zvpll_av + zvpll*zdtg ! (frame corrected...)
  !
  !      treat case where ztim + zdtg > delt
  !      nbflg=1 means an ion which has survived
  !      full beam code timestep, and is still orbiting
  !
  if((bitime(1)+zdtg).ge.delt) then
     nbflg=1
     zdtg=(delt-bitime(1))
     if(zdtg.le.0.0_fp) goto 890
  end if

  !

  ztslmg=ztslmg+zdtg
  !
  !  mod dmc jan 1992 -- heavy ion beams -- first few steps integrate
  !  around gyro orbit explicitly (iflri=1) so avoid calls to limlos
  !  and gphase to avoid spurious loss and randomization of gyrophase
  !  angle
  !
  !  check for flr orbit loss -- (upgrade rjg 7/21/83)
  if(iflri.eq.0) then
     if (imap_larmor_loss .and. (rlar>max(epslon,gflr_min))) then
        call map_larmor_loss(ilim)
     else
        call xlimlos(ilim)
     end if
     if(ilim.eq.1) then
        !  if the larmor radius is equal to or larger than the distance to the
        !  limiter, we treat this as a lost orbit
        icheck=0
        goto 800
     end if
  end if
  !
  !  if doing flr, evaluate random gyrophase angle
  !
  if(nlbflr) then
     if(iflri.eq.0) call gphase
  end if
  !
  !  find the gyro position in (r,y), (xi,theta (approx)), and indices
  !  (if not doing flr, copy the g.c. coords.)
  !
  call nbiflr(zdummy,0)
  !
  !  if outside beam dist fcn zone region this is deemed an orbit loss!
  !
  if(nflr2.eq.0) then
     ilim=1
     icheck=0
     goto 800
  end if
  !
  !  check if outside plasma region on gyro
  if(.not.nlheavy(lsbeam)) then
     if(nlpout.and.(.not.nlfprod(lsbeam) .or. nlfatom)) then ! tbt
        !
        !  dmc-- force goose, etc =1 for beam orbits intersecting high-
        !   neutral density region outside plasma -- if energy is low enough
        if(iswgoose.eq.0) then
           if((ep*v2toev/abeam).le.100.0e3_fp) then
              !  force full half-orbit before goose is re-evaluated
              iswgoose=1            ! do this at most once per orbit track
              goose=min(1.0001_fp,goose)
              ngoose=0
           end if
        end if
        !  force cx calculation this step
        cxskip=1.0_fp
        inexcx=istpor
     end if
  end if
  !%
  !%           if(ipflr.ne.0) then
  !%            if (mod(istep-1,abs(ipflr)) .eq. 0) then
  !%              zzx2(1)=2.0
  !%              zzy2(1)=2.0
  !%              if (nbplot .eq. 1) then
  !%                zzx2(2)=rmjion
  !%                zzy2(2)=yion
  !%                zzx2(3)=rmjflr
  !%                zzy2(3)=yflr
  !%              else
  !%                zzx2(2)=vpll
  !%                zzy2(2)=vperp
  !%                zzx2(3)=vplflr
  !%                zzy2(3)=vppflr
  !%              end if
  !%              call symbl(0)
  !%              call line(0)
  !%              call r8_cplot(zzx2,zzy2)
  !%              if (nbplot .eq. 1) then
  !%                zzx(2) = rmjflr
  !%                zzy(2) = yflr
  !%              else
  !%                zzx(2) = vplflr
  !%                zzy(2) = vppflr
  !%              end if
  !%              call symbl(1)
  !%              call line(-1)
  !%              call r8_cplot(zzx,zzy)
  !%              if ((.not. nlbgflr) .and. ipflr<0) then
  !%               call map_larmor(ioops, iregion, zzx(2),zzy(2),
  !%     >              z_vpar, z_vperp,alpha_v, .true.)
  !%               ioops = min(11,8 + ioops)  ! change symbol with error code
  !%               call symbl(ioops)
  !%               call line(-1)
  !%               if (nbplot .eq. 2) then
  !%                zzx(2) = z_vpar
  !%                zzy(2) = z_vperp
  !%               end if
  !%               call r8_cplot(zzx,zzy)
  !%              end if
  !%              if (nlbgflr) call map_larmor_compare(istep)
  !%            end if
  !%           end if
  !%
  !
  !
  !  **************  every orbit step calculations  ******************
  !
  !
  !  count orbit timesteps
  istpor=istpor+1
  !
  !  keep track of beam code time
  zdtime = zdtime + zdtg
  bitime(1)=zjtim + zdtime
  !
  !  avg-j sum
  zjsum=zjsum+zdtg*ngc
  !
  !--------------------------------------------
  !  check for fishbone
  if(.not.nlfbon.or.ngoose.lt.2) goto 120
  ztmtra=tbm1+bitime(1)
  if(ztmtra.lt.tfshon.or.ztmtra.gt.tfshof) goto 120
  ztmfsh=ztmtra-tfshon
  if(mod(ztmfsh,fshper).gt.fshwid) goto 120
  !
  !  take action if you got a fishbone
  zevolt=abeam*ep/evptocm_sec/evptocm_sec    !*5.221e-13_fp
  call fshlos(zvpvc1,zrc1,zthc1,zvpvc2,zrc2,zthc2, &
       zevolt,zfrate)
  zflive=exp(-zdtg*zfrate)
  if(nbrandm_pll(idum).lt.zflive) goto 120
  ifkill=1
  goto 810
  !
120 continue
  !
  !  check for ripple loss - srd model
  !           if first orbit, don't check ripple loss
  if (goose.lt.1.00001_fp) go to 130
  if(nrip.lt.2) goto 130
  if(nrip.eq.2) then
     call riplos2(zksid,iripl,ierr)
     if(ierr.ne.0) return
  end if
  !
  if(iripl.eq.1) goto 820
  !
130 continue

  !--------------------------------------------
  !
  !  contribute to beam distribution function  (r.j.g. april 11, 1982)
  !      (every orbit time step, for best possible statistics !)
  !      (rminbm .ne. rminor upgrade rjg 7/21/83)
  !
  zde=bidte(1)*(1.e-7_fp/xe)
  zdph=bidpph(1)
  !
  zderrmin=min(zderrmin,zde)
  zderrmax=max(zderrmax,zde)
  zderrsum=zderrsum+abs(zde)
  !
  zdprrmin=min(zdprrmin,zdph)
  zdprrmax=max(zdprrmax,zdph)
  zdprrsum=zdprrsum+abs(zdph)
  !
  bperrs(lsbeam)=bperrs(lsbeam)+wght*zde
  bpherrs(lsbeam)=bpherrs(lsbeam)+wght*zdph


  !
  jbe=ndxfbe_r8(v2toev*ep,lsbeam,jbe)
  if(jbe <= 0) then
     !        print *, ' ?ndxfbe_r8 : out of energy',
     !     >           jbe, v2toev*ep, ebdmaxa(lsbeam), zmaxe
     jbe = 1
  end if
  !
  !  pitch angle, energy indices (nflr2 = spatial index, cf nbiflr)
  !
  jba=int((1.0_fp+xksid)*0.5_fp*real(nznbma))+1
  jba=min(nznbma,max(1,jba))
  !
  if(.not.ilfbm) goto 150
  !
  !----
  !  the distribution function...
  !mg         if(ep.ge.epbdmx) goto 150
  if(jbe.gt.nznbmea(lsbeam)) goto 150
  if(nlpout) goto 150
  !
  !  dmc feb 2009 -- based on mpi performance data, direct accumulation
  !  of fbm and fbm_ptcl is tried, rather than using zfbml_* with subsequent
  !  flipr calls.  if this helps performance, the changes will be kept.
  !
  if(nlfbmflr) then
     if((nflr2.gt.0).and.(nflr2.le.mimxbzf)) then
        fbm_ptcl(jbe,jba,nflr2,lsbeam,lgyro) = &
             fbm_ptcl(jbe,jba,nflr2,lsbeam,lgyro) + zdtg*wght
     end if
  end if !nlfbmflr
  if((ngc2.gt.0).and.(ngc2.le.mimxbzf)) then
     fbm(jbe,jba,ngc2,lsbeam) = &
          fbm(jbe,jba,ngc2,lsbeam) + zdtg*wght
  end if

150 continue
  !
  ! eval sum at guiding center
  ! commented by mg jan 30, 2012
  !         fnb(jbe,ngc,lsbeam)=fnb(jbe,ngc,lsbeam) + zdtg*wght
  !         fvdotb(jbe,ngc,lsbeam)=fvdotb(jbe,ngc,lsbeam) +
  !     >        vpll*bimodb(1)*zdtg*wght
  ! end mg
  !

  !  dmc jul 1992 -- density sums in restricted energy ranges
  !
  zep=v2toev*ep
  if(orbrzv_rf_option.eq.1) then
     if(zep.gt.efbm(nznbmea(lsbeam),lsbeam)) then
        do ie=1,nerngfi
           densfi(nflr,ie,lsbeam)=densfi(nflr,ie,lsbeam) + zdtg*wght
        end do
     end if
  else
     if(nerngfi.gt.0) then
        do ie=1,nerngfi
           ze1=erngfi(ie,lsbeam)
           if(ie.lt.nerngfi) then
              ze2=erngfi(ie+1,lsbeam)
           else
              ze2=epsinv
           end if
           if((ze1.le.zep).and.(zep.lt.ze2)) then
              densfi(nflr,ie,lsbeam)=densfi(nflr,ie,lsbeam) + &
                   zdtg*wght
           end if
        end do
     end if
  end if

  !
  !==================================================================
  !  following calculations are performed periodically, not necessarily
  !  every orbit step
  !
  !
  !***************    atomic collisions  ***************************
200 continue
  isw=0
  !
  !
  !  do cx interactions 1./cxskip (fraction) of orbit timesteps
  !      beam ion -- beam + plasma neutrals interactions
  !       (charge exchange, impact ionization)
  !
  !
  !  dmc sept 1990 -- for now, atomic physics is disabled for
  !  fusion products *** !  unless explicitly requested by nlfatom - tbt.
  !
  if(nlfprod(lsbeam) .and. .not. nlfatom) goto 300   ! tbt 11/91
  !
  if(istpor.lt.inexcx) goto 300
  !
  !  calculate when to do next collision
  iinc=cxskip
  if((cxskip-iinc).gt.nbrandm_pll(idum)) iinc=iinc+1
  inexcx=inexcx+iinc
  !
  !  evaluate goosed time-step for c-x, enhanced by cxskip
  dtgcx=zdtg*cxskip
  !
  !  increment cx calculations counter and accumulator
  istpcx=istpcx+1
  zdltcx=ttcx-ztscxp
  !
  !
  !%      if(.not.((iplcx.ne.1).or.(nbplot.le.0))) then
  !%        if(nbplot.eq.1) then
  !%        zzx(2)=rmjflr
  !%        zzy(2)=yflr
  !%        else
  !%        zzx(2)=vpll
  !%        zzy(2)=vperp
  !%        end if
  !%        call symbl(2)
  !%        call line(-1)
  !%        call r8_cplot(zzx,zzy)
  !%      end if
  !
  if((.not.nlorbo).and.(iplcx.ne.-1)) then
     !
     if(nlheavy(lsbeam)) then
        call beamiz(isw,ievent,iflri)
        if(ievent.eq.2) then
           icheck=1
           ilim=1
           goto 800                ! lost orbit
        end if
     else
        if(nlbeamcx) then
           call beamcx(isw,ievent,ierr)
           if(ierr.ne.0) return
        end if
     end if
     !
  end if
  !
  !       if ievent=1 nutrav was called -- cx took place
  !
  !      isw=1 if we are now following a new ion orbit due to
  !      charge-exchange transport
  !
  !      nbicx=1 if the particle has been split down to zero and lost
  !
  !  if ion has left as a c-x neutral, then exit from orbit code
  !
  if(nbicx.eq.1) goto 880
  !
  !  if following an old orbit, without any call to nutrav,
  !      leave c-x section
  if(ievent.eq.0) goto 300
  !
  !  deal with new orbits
  !
  if(isw.eq.1) then
     !  check for time expiration - otherwise reenter orbit like a fresh
     !  ptcl
     if(nbflg.eq.1) goto 890
     if(.not.nlheavy(lsbeam)) then
        !   jphe ! 23-feb-1994
        !dmc        if  (nlmcr2) then
        !dmc           delrs=epsinv
        !dmc           rscon=0.
        !dmc           drscon=0.
        !dmc           lrsok=.false.
        !dmc           lgosrf=.false.
        !dmc           rgosrf(1)=1.
        !dmc           rgosrf(2)=epsinv
        !dmc        end if
        !
        !   jphe ! 23-feb-1994
        !
        goto 5                     ! after cx recapture (h/he fast ion)
     else
        isumiz=isumiz+1             ! ionization level promotion: heavy ion
     end if
  end if
  !
  !  no new orbit, but nutrav was called.  therefore we must re-initialize
  !      summing for charge exchange
  !
  ztscxp=-zdltcx
  !
  !
  !*****************  fokker - planck collisions  *****************
300 continue

  !
  !  test whether to do fp this step
  if(istpor.lt.inexfp) goto 400
  !
  !  do it, but first set up when to do the one after...
  zfpskip=max(1.0_fp,fpskip)
  iinc=zfpskip
  if((zfpskip-iinc).gt.nbrandm_pll(idum)) iinc=iinc+1
  inexfp=inexfp+iinc
  !
  !  set up goosed time-step enhanced by fpskip
  dtgfp=zdtg*zfpskip
  ztimfp=ztimfp+dtgfp
  !
  !    increment fokker-planck counter
  istpfp=istpfp+1
  !%      if(.not.((iplfp.ne.1).or.(nbplot.le.0))) then
  !%        if(nbplot.eq.1) then
  !%        zzx(2)=rmjflr
  !%        zzy(2)=yflr
  !%        else
  !%        zzx(2)=vpll
  !%        zzy(2)=vperp
  !%        end if
  !%        call symbl(5)
  !%        call line(-1)
  !%        call r8_cplot(zzx,zzy)
  !%      end if
  !
  !  call fp collid subroutine, and keep track of cpu time used
  !
  !xxx      if(.not.nlorbo) then
  if(iplfp.ne.-1) then
     !
     zepp=ep                        !test
     zphp=bphion                    !test
     !----------------------------------
     icol=icol+1
     !  collisions in the plasma frame:
     bitime_save=bitime(1)
     bitime(1)=bitime(1)-zdtg+dtgfp !to pass to aep_driver
     call collid
     
     if (ilaep_init(lsbeam).and.(orbrzv_option.eq.1)) then
        il_rz_updated=.true.
     end if

     bitime(1)=bitime_save !back to original
     !
     !  if ion is thermalized, exit from orbit code
     if(nbflg.eq.2) goto 850

     !
     !  anomalous diffusion; non-orbit effects that move the g.c.
     rmjion_prev=rmjion
     yion_prev=yion
     call abdifu(jbe,zdsdfb,zdsxvb,ierr)
     if(ierr.ne.0) return
     if (orbrzv_option.eq.1) then
        if((rmjion_prev.ne.rmjion).or.(yion_prev.ne.yion)) then
           il_rz_updated=.true.
        end if
     end if
     !
     zxion_min=min(xion,zxion_min)
     !
     zepp_jxb=v2toev*vion*vion
     !
     !----------------------------------
     !
  else
     istpfp = 0                     !default for collisionless orbits
  end if
  !
  !  if ion escapes (e.g. via anomalous diffusion) exit from orbit code
  if(ngc.eq.0) then
     icheck=1
     ilim=1
     goto 800
  end if
  !
  !***************    compression  ***************************
400 continue
  !
  !  compression calculation moved to orball -- taken out of ptcl orbit loop
  !  rfkick removed altogether, until real work can be done on it...
  !
  !****************  goose factor updating  *******************
  !
  !  re-evaluate the goose factor at each midplane crossing
  !
  if(.not.bounce_chk(itrig)) goto 594
  !  at midplane crossing, control comes here
  if(itrig.eq.1 .and. orbrzv_rf_option.eq.1.and.mcrf_on) then
     ! count the # of mid-plane crossing
     orbrzv_nmidc = orbrzv_nmidc + 1
     ! calculate orbit width
     zrhopm = dsqrt(dabs(orbrzv_psiv(1))/orbrzv_psix)
     if(orbrzv_nmidc > 1) then
        orbrzv_bnw = dabs(zrhopm-orbrzv_rhopm)
        zrsg = orbrzv_rmid0 - orbrzv_rc
        zrsg = zrsg*(orbrzv_zz(1) - orbrzv_rc)
        orbrzv_trpas = sign(1.0_fp, zrsg)
     else
        ! the start of cross counting
        ! initialize the width and trpas
        orbrzv_bnw = 0.0_fp
        orbrzv_trpas = 0.0_fp
     end if
     orbrzv_bnw = orbrzv_bnw*(orbrzv_prmax-orbrzv_rc)
     orbrzv_rhopm = zrhopm
     ! update mid-plane crossing point
     orbrzv_rmid0 = orbrzv_zz(1)
  else
     ! if fail to locate the mid-plane crossing
     ! estimate orbit width conservatively (large)
     orbrzv_nmidc = 0
     orbrzv_bnw = 0.1_fp*(orbrzv_prmax-orbrzv_rc)
     orbrzv_trpas = 0.0_fp
  end if
  if(nlmcrfgcon) then
     zgooconrf = max(goocon, orbrzv_bnw_lrf)
     !zgooconrf = max(goocon, orbrzv_bnw_lrf/goocon)
  else
     zgooconrf = goocon
  end if
  !
  !
  !
  !
  !  default [floating] value of steps per half-orbit
  zstpor = real(istpor,fp)
  !  plotting gub
  !
  !  midplane crossing -- check for ripple loss
  !
  if(nrip.eq.1) then
     call riplos(ztslmg,iripl,ierr)
     if(ierr.ne.0) return
     if(iripl.eq.1) goto 820
  end if
  !
  !  remember old goose factor
  zgoo=goose
  !
  !  keep track of number of crossings
  ngoose=ngoose+itrig
  !  make sure orbit does not somehow get trapped beyond the plasma
  if((ngoose.ge.2).and.(zxion_min.gt.1.0_fp)) then
     call infolog &
          ('orbit: an orbit track is stuck entirely outside the plasma')
     call infolog(' min(xion)=',zxion_min)
     if(myidd.eq.0) then
        write(nonlin,9518) zxion_min
     end if
9518 format(' %orbit -- an orbit track is stuck entirely outside'/ &
          '  the plasma:  min(xion)=',1pe11.4)
     ilim=1                         ! treat as loss orbit
     icheck=1
     goto 800
  end if
  !
  if(nlorbo) goto 550  ! do not change goose if not slowing down
  !
  !  keep track of xksid, ep, th at last two crossing for fishbones
  if(nlfbon) then
     zvpvc1=zvpvc2
     zvpvc2=xksid
     zrc1=zrc2
     zrc2=xion
     zthc1=zthc2
     zthc2=th
  end if
  if(itrig.eq.1) then
     zmodb_midpp = zmodb_midp
     zmodb_midp = max(bimodb(1),zmodb_prev)
  end if
  !
  !  if you don't have the info to change goose, don't try to
  !
  if(ngoose.lt.2) goto 550
  !
  !  update field topology & bounce average information
  !
  zvpll_av = zvpll_av/zdtsum
  !
  dpass = -1.0_fp
  zxi=zx_bmax
  if(itrig.eq.2) then
     dtrap = 1.0_fp         ! off midplane trap
  else if(zmaxbo.eq.zminbo) then
     !   b field unchanged -- must have been 1 step orbit; assume deeply trapped
     dtrap = 1.0_fp
  else
     zxi=min(1.0_fp,zxi)
     igcmax = lcentr + zxi*nzones
     igcmax = min(ledge, igcmax)
     zxfac = (zxi-xiblo(igcmax))/ &
          (xiblo(igcmax+1)-xiblo(igcmax))
     zbmax_fs=bmaxm(igcmax) + &
          zxfac*(bmaxm(igcmax+1)-bmaxm(igcmax))
     zbmin_fs=bminm(igcmax) + &
          zxfac*(bminm(igcmax+1)-bminm(igcmax))

     zmodb=max(zmodb_midp,zmodb_midpp)

     ! passing orbit assumed if max(b) of half bounce is at fs max or
     ! at the max for the half-bounce
     zratb=max(  (zmodb-zminbo)/(zmaxbo-zminbo), &
          (zmaxbo-zbmin_fs)/(zbmax_fs-zbmin_fs) )
     zpass = 1.0_fp - zmaxbo/zbrefl_bmax ! "passing" depth
     if(zratb.gt.0.98_fp) then
        ! treat as passing orbit: max orbit |b| at midplane or flux surface maximum
        dtrap = -1.0_fp
        dpass = zpass
     else
        if(zpass.gt.0.02_fp) then
           ! override: no magnetic reflection seen
           !dbg               write(nonlin,*) ' %orbit: trap override'
           dtrap = -1.0_fp
           dpass = zpass
        else
           dtrap = (zbmax_fs-zmaxbo)/(zbmax_fs-zbmin_fs)
           dtrap = max(0.0_fp,min(1.0_fp,dtrap))
        end if
     end if
  end if
  !
  if(nmdifb.eq.4 .or. nmdifb.eq.6) then
     ! *ae mode radial transport operator
     call bbdifu(zvpll_av,ierr)
     if(ierr.ne.0) return
  end if
  !
  !-------------------------------------
  !dbg      write(nonlin,*) '  iorb, ngoose, #steps = ',iorb,ngoose,istpor
  !dbg      write(nonlin,*) '  xion,eion = ',xion,v2toev*vion*vion
  !dbg      write(nonlin,*) '   de min,max,abs_sum:  ',
  !dbg     >   zderrmin,zderrmax,zderrsum
  !dbg      write(nonlin,*) '   ---> [de abs_sum]/einj = ',
  !dbg     >   zderrsum/einjs(lsbeam)
  !dbg      write(nonlin,*) '   dph min,max,abs_sum: ',
  !dbg     >   zdprrmin,zdprrmax,zdprrsum
  !dbg      write(nonlin,*) '   ---> [dph abs_sum]/ph_ref = ',
  !dbg     >   zdprrsum/xmrefs(lsbeam)
  !-------------------------------------
  !
  if(istpfp.eq.0) then
     fpskip=min(1.0_fp,fpskip)
     cxskip=1._fp
     if(istpor.lt.4) then
        goose=1.0001_fp
     end if
     goto 550
  end if
  !
  !  remember that you have had at least one pass into changing goose
  !      (this will be checked at final exit from orbit code)
  !
  igoosc=1
  !
  !  evaluate factors which take into account whether you got the
  !      expected fraction of orbit timesteps also in collid and cxsr
  !
  zfacfp=zstpor/(max(1.0_fp,fpskip)*real(istpfp))
  zfaccx=0.0_fp
  if(istpcx.gt.0) zfaccx=zstpor/(cxskip*real(istpcx))
  !
  !  evaluate "how much" cx, colliding, and compression/ free expansion
  !       are done in a half-orbit,
  !      for use in re-evaluating goose
  !
  !  charge-exchange
  zdltcx=zfaccx*(tscx-ztscxp)
  !
  !      each ion represents wght/wbav(lcentr) "potential neutrals"
  !      for charge-exchange; ch-ex sum is calculated with this en-
  !      hancement; remove for purposes of reevaluating goose factor
  !      (d. mc cune 6 apr 1981)
  !
  if(wght.eq.0.0_fp) then
     !     ptcl. with wght==0 not supposed to be orbiting, probably came from spliting
     !     in to the region with zero density, for instance when dn0out = 0.0
     call errorlog(' orbit:: wght==0')
  end if
  ztbtcx=zdltcx*wbav(lcentr,lsbeam)*fact_wbav/wght
  if(ztbtcx.eq.0.0_fp) then
     zgooscx=zgoomax
  else
     ! kjm-rfgoose-goocon zgooscx=goose/(goocon*ztbtcx)
     zgooscx=goose/(zgooconrf*ztbtcx)
  end if
  !
  !  fokker-plancking
  if(.not.nlorbo) then
     viplasa=zfacfp*viplasa/ztimfp    ! avg plasma frame ion velocity
     zdvgef=zfacfp*dvg
     zfdvg=zdvgef/viplasa
  else
     zfdvg=0.0_fp
  endif
  !      new test 17 apr 1981  d. mc cune seperating linear and
  !      quadratic parts as per rj goldston suggestion
  dksgc=abs(zfacfp*dksgc)
  ztbtsl=max(zfdvg,dksgc)
  !
  !  re-evaluate goose just taking into account drag, convective
  !      change in xksid, and charge exchange
  if(ztbtsl.eq.0.0_fp) then
     zgoosl=zgoomax
  else
     ! kjm-rfgoose-goocon zgoosl=goose/(goocon*ztbtsl)
     zgoosl=goose/(zgooconrf*ztbtsl)
  end if
  !
  !
  !  make sure that the time step does not permit delta - e
  !      per half bounce greater than 1/2 of the energy bin
  !      width if you are calculating a beam distribution fct.
  !
  !  dmc -- make this restriction more stringent if ngoocon_ebin > 1 is set
  !
  !  dmc - c% lines to check for dead zero diffusion.  plotting vsn of
  !  code now has collisionless orbit option (june 1985)
  !
  if(ilfbm) then
     zfdv2g=(2.0_fp*viplasa + zdvgef)*zdvgef
     if(zfdv2g.eq.0.0_fp) then
        zgsbde=zgoomax
     else
        if(jbe.ge.1.and.jbe.le.nznbmea(lsbeam)) then
           zgsbde=goose*(efbmb(jbe+1,lsbeam)-efbmb(jbe,lsbeam)) &
                /(ngoocon_ebin*v2toev*zfdv2g*2.0_fp)
        else
           zgsbde=zgoomax
        end if
     end if
     zgoosl=min(zgoosl,zgsbde)
  end if
  !
  !  take into account diffusive change in xksid
  !  dmc -- condition made more stringent if ngoocon_vpvbin is set...
  !
  ztbtsq=zfacfp*dksgd
  if(ztbtsq.le.epslon) then
     zgoosq=zgoomax
  else
     ! kjm-rfgoose-goocon zgoosq=goose/(goocon**2*ztbtsq)
     zgoosq=goose/(zgooconrf**2*ztbtsq)  ! tau(90 degree scatter)
     if(ngoocon_vpvbin.gt.0) then
        ! limit according to p.a. scatter time across an fbm pitch zone
        zgoosq2=goose/(ngoocon_vpvbin**2*(nznbma**2*ztbtsq/4))
        zgoosq=min(zgoosq,zgoosq2)
     end if
     if(nmdifb.eq.4) then
        zgoosq2=goose/(ztbtsq/dif2d_dpa**2)
        zgoosq=min(zgoosq,zgoosq2)
     end if
  end if
  !
  !  take into account anomalous beam ion radial diffusion
  if(nmdifb.eq.0) then
     zgoosd=zgoomax
  else
     zfacd=zfacfp
     ztbdsq=zfacd*zdsdfb    ! abdifu
     zdelxi=0.5_fp          ! max delta(xi)/goocon bounces
     zdelxi2=dxitk/2.0_fp   ! half an fbm zone width
     if(ztbdsq.le.0.0_fp) then
        zgoosd=zgoomax
     else
        !  dont want more than one half beam zone width crossing/ half orbit
        !kjm-rfgoose-goocon zgoodf=max((goocon/zdelxi)**2,(2.0_fp/dxitk)**2)
        zgoodf=max((zgooconrf/zdelxi)**2,(ngoocon_rbin/zdelxi2)**2)
        zgoosd=goose/(zgoodf*ztbdsq)
     end if
     ztmp=zfacd*zdsxvb
     if(ztmp.gt.0.0_fp) then
        zgoodf=max(abs(goocon/zdelxi),abs(ngoocon_rbin/zdelxi2))
        zgoosd=min(zgoosd,(goose/(zgoodf*zdsxvb)))
     end if
  end if
  !
  !  d/dt(enclosed flux) (eldot)
  !  (defunct)
  zgoosx=zgoomax
  !
  goose=max(1.0001_fp, &
       min(zgooscx,zgoosl,zgoosd,zgoosq,zgoosx,zgoomax))
  !
  !  goose kept at one if heavy ion is still being stripped of electrons
  !
  if(isumiz.gt.0) then
     goose=1.0001_fp
  end if
  !
  !  careful of short orbits
  !
  if(istpor.lt.4) then
     goose=1.0001_fp
  end if
  !
  !  d. mc cune 29 mar 1983-- restrict goosing s.t. at least one
  !  orbit should be completed in a single beam timestep
  !    ... use ztslmp= (ungoosed) orbit time since last midplane crossing
  !
  imin_mpc = 3  ! #midplane crossings to have 1 full orbit guaranteed
  if(nmdifb.eq.4) imin_mpc = imin_mpc + 1
  !
  zgtes=delt/(goose*imin_mpc*(ztslmp))
  zgooso=zgtes*goose
  goose=max(1.0001_fp,min(goose,(zgooso)))
  !
  !  maximum allowable skip factor-- force avg of 3 events/ half orbit
  zmaxsk=max(1.0_fp,zstpor/3.0_fp)
  !
  !  evaluate the number of orbit steps per fokker-planck collision
  !
  fpskip=min(zmaxsk,max(1.0_fp,(zstpor/fppcon)))
  !
  !  empirical test criterion...
  zgtest=3.0_fp*xzbeami
  !
  if(goose.lt.zgtest) then
     !  if goose is small and collisional processes are fast, reduce
     !  fpskip further...
     zgcomin=min(zgoosd,zgoosq,zgoosl)
     fpskip=fpskip*min(1.0_fp,zgcomin/zgtest)
  end if
  !
  !  evaluate the number of orbit steps per c-x collision
  !      minimum of 3 ch-ex calculations per orbit required
  !
  if(istpcx.ge.3) then
     if(zdltcx.eq.0.0_fp) then
        cxskip=zmaxsk
     else
        cxskip=zstpor*zgoo/(cxpcon*zdltcx*goose)
     end if
  else
     cxskip=max(1.0_fp,fpskip)
  end if
  if(cxskip.lt.1.0_fp) then
     zgoost=cxskip*goose
     goose=max(1.0001_fp,(zgoost))
     cxskip=1.0_fp
  end if
  !
  !  upper limit
  cxskip=min(zmaxsk,cxskip)
  if(nlorbo) then
     fpskip=1.0_fp
     cxskip=1.0_fp
  endif
  !
  if(itrig.eq.2) then
     !  halve goose and reset ngoose, if this is a goose update
     !  prior to midplane crossing, due to bounce time expiration...
     ngoose=0
     goose=max(1.0001_fp,goose/2)
  else
     if(icheck_dt) then
        !  new orbits only...
        zderrsum=zderrsum/eirefs(lsbeam)
        zdprrsum=zdprrsum/xmrefs(lsbeam)
        zadjust=(max(zderrsum,zdprrsum)/dt_acc)**0.3333333_fp
        !  min/max limits allow for slow lengthening but rapid shortening
        !  of timestep...
        zadjust=max(0.8_fp,min(4.0_fp,zadjust))
        xnum_dt_adj0(lsbeam)=xnum_dt_adj0(lsbeam)+1.0_fp
        dt_adj0_sums(lsbeam)=dt_adj0_sums(lsbeam)+(1.0_fp/zadjust)
     end if
  end if
  !
550 continue
  !     kjm-rfgoose-2
  if(nlmcrfgcon) then
     if(orbrzv_rfdv(0) .gt. 0.1_fp) then
        ! update zgooconrf (based on length scale)
        ! to control rf-induced change
        zgooconrf = max(zgooconrf, 20.0_fp)
        zrgdv = (zgooconrf*orbrzv_rfdv(1))**2
        zrgdx = (zgooconrf*orbrzv_rfdv(2))**2
        zgoorf = max(zrgdv, zrgdx)
        zgoorf = zgoo/(zgoorf + epslon)
        goose = max(1.0001_fp, min(zgoorf, goose))
     end if
     if(itrig.eq.1) then
        ! keep record the # of rf-kick in the previous half-bounce
        orbrzv_rfdv(-1) = orbrzv_rfdv(0)
        orbrzv_rfdv(0:2) = 0.0_fp
        orbrzv_zrfd_min = 1.0_fp
     end if
  end if
  !
  !  re-initialize bounce sums for next goose-update cycle
  !
  call bounce_ini(0)
  !
594 continue
  !
  !*************  bottom of main orbit loop **********************
  !
  !
  !%      if((jt.ne.izplot-1).and.(ztesm*ztesp.gt.0.0)) goto 595
  !%      if(nbplot.le.0) goto 595
  !%      zxpos(1)=jt-1.0
  !%      zypos(1)=jt-1.0
  !%      call line(0)
  !%      call symbl(0)
  !%      call r8_cplot(zxpos,zypos)
  !%      call tsend
  !%      call line(-1)
  !%      zxpos(2)=zxpos(jt)
  !%      zypos(2)=zypos(jt)
  !%      jt=2
  !%      iplct=iplct+1
  !%      if(iplct.ge.iplnu) then
  !%        call stplot
  !%     & (0,iplnu,iplct,iplcx,iplfp,iplsty,ipflr,
  !%     &  bitime(1))
  !%      else
  !%      end if
  !%      if(iplct.eq.0) jt=1
  !%595      continue
  !
  !  return to orbit if bitime(1) is not now up to delt
  if(nbflg.eq.1) goto 890
  goto 100
  !
  !
  !*********************  exit points  ************************
  !
  !  particle has left beam ion distribution by
  !      bad orbit loss,
  !      or fishbones (!),
  !      or thermalization,
  !      or charge-exchange
  !       or ripple
  !
800 continue
  !
  !  this section for bad orbit ("limiter") loss
  !
  !  dmc may 2005 -- edge jxb correction associated with orbit loss
  !   per comments rob akers -- mast tokamak
  !
  if(ijxb.eq.1) then
     zxip=xion
  end if
  if(zxip.lt.1.0_fp) then
     call psi_spline(zxip,ingc,zpph1)
     call psi_spline(1.000_fp,ingc,zpph2)
     zeph1=xzbeami*fnbint(zxip,phiprg,1,mj,1,idum,ingc)
     zeph2=xzbeami*phiprg(lep1)
     !  p_lim = p_gc+mp*zb*(e/c)*(psi_lim-psi_gc)
     zbphip=bphion+9.578e11_fp*rmassb/abeam*xzbeami*(zpph2-zpph1)
     zepp=zepp_jxb+zeph1-zeph2 ! potential loss = kinetic gain
     !
     call nb0jxb(wght,zepp_jxb,zepp,bphion,zbphip,btori,btori, &
          ngc,ledge,zxip,1.0_fp,5,.false.)
     !
     xion=1.0_fp
     zepp_jxb=zepp
     bphion=zbphip
  end if
  !
  bpls(lsbeam)=bpls(lsbeam)+zepp_jxb*wght
  bxflws(lsbeam)=bxflws(lsbeam)+wght
  bphorbs(lsbeam)=bphorbs(lsbeam)+bphion*wght
  if(abs(goose-1.000000_fp).lt.1.0e-6_fp) then
     bxflws_prompt(lsbeam)=bxflws_prompt(lsbeam)+wght
     xmcloss_o(lsbeam)=xmcloss_o(lsbeam)+1.0_fp ! prompt loss
  end if
  !
  xistop=1.01_fp  ! ion has left the plasma
  !
  ! --- splines ---
  ier=0
  zphiarg=0.0_fp
  call xs_fxget(0.01_fp,zphiarg,rmjion, yion, zxichk, zthchk,iwarn)
  if ((ier.eq.0).and.(icheck.eq.1).and.(zxichk .lt. xibado)) then
     call errorlog(' orbit: bad orbit, xi = ', zxichk)
     call errorlog( ' lsbeam = ',lsbeam)
     call errorlog('   iorb = ',iorb)
     call errorlog('  itrak = ',itrak)
     call errorlog('  istep = ',istep)
     if(myidd.eq.0) then
        write (nonlin, 8003 ) zxichk
        write(nonlin,*) ' lsbeam = ',lsbeam
        write(nonlin,*) '   iorb = ',iorb
        write(nonlin,*) '  itrak = ',itrak
        write(nonlin,*) '  istep = ',istep
     end if
8003 format (/' ?orbit -- bad orbit, xi =',f10.5/)
     ierr=1
     return
  end if

  if(nlbout) then
     gtime=tbm1+bitime(1)
     yze=ep*v2toev
     if(nlsym2b) then
        zhgt=plhgt
     else
        zhgt=0.0_fp
     end if
     yyion=yion+zhgt
     ziripl=real(iripl)
     ynbsce = nbsce
     ynbien = nbien
     call store_lost_ptcl
  end if
  !
  goto 880
  !
  !  this section for fishbone loss
810 continue
  bpfshs(lsbeam)=bpfshs(lsbeam)+ep*v2toev*wght
  bxfshs(lsbeam)=bxfshs(lsbeam)+wght
  bphfshs(lsbeam)=bphfshs(lsbeam)+bphion*wght
  !
  xistop=xion       ! orbit ends at fishbone
  call radcur(xion,1.01_fp,wght,tirfshs)  ! fishbone radial current
  !
  goto 880
  !
  !  this section for ripple loss
820 continue
  bprpls(lsbeam)=bprpls(lsbeam)+ep*v2toev*wght
  bxrpls(lsbeam)=bxrpls(lsbeam)+wght
  bphrpls(lsbeam)=bphrpls(lsbeam)+bphion*wght
  !
  xistop=xion       ! orbit ends at ripple event
  call radcur(xion,1.01_fp,wght,tirrips)
  !
  if(nlbout) then
     gtime=tbm1+bitime(1)
     yze=ep*v2toev
     if(nlsym2b) then
        zhgt=plhgt
     else
        zhgt=0.0_fp
     end if
     yyion=yion+zhgt
     ziripl=real(iripl)
     ynbsce = nbsce
     ynbien = nbien
     call store_lost_ptcl
  end if
  !%       if((nrip.eq.2).and.(nbplot.gt.0)) then
  !%            if(nbplot.eq.1) then
  !%              zzx(2)=rmjflr
  !%              zzy(2)=yflr
  !%            else
  !%              zzx(2)=vpll
  !%              zzy(2)=vperp
  !%            end if
  !%            call symbl(ichar('R'))
  !%            call line(-1)
  !%            call r8_cplot(zzx,zzy)
  !%      end if
  goto 880
  !
  !  thermalization (in subroutine collid) sends you here
  !      with nbflg=2, and power flow accounted for already
  !
850 continue
  !
  !%      write(zbuf,7701)
  !%7701      format('thermalized$')
  !
  nbflg=2
  goto 890
  !
  !  bad orbit loss, c-x loss, ripple loss and f-bone loss converge here
  !      nbflg=3 means an ion lost on a bad orbit, by c-x, by ripple,
  !       or by f-bones
  !
880 continue
  nbflg=3
  !
  !%      write(zbuf,7702)
  !%7702      format('ptcl escaped$')
  !%      if(nbicx.eq.1) write(zbuf,77021)
  !%77021      format('ptcl exit:cx$')
  !%      if(ifkill.eq.1) write(zbuf,77022)
  !%77022      format('fishbone loss$')
  !%      if(iripl.eq.1) write(zbuf,77031)
  !%77031      format('ripple loss$')
  !
  !  all exits from orbit pass thru here
890 continue
  !     kjm-rfdbg
  !#ifdef __DEBUG
  !      if(orbrzv_option.eq.1 .and. orbrzv_rf_option.eq.1) then
  !        call my_rfdist_diag(iorb, istep, ierr)
  !      end if
  !#end if
  !
  !  dmc 30 sept 1994 -- if never inside plasma (last closed flux
  !  surface) call it orbit loss.
  !
  if(nbflg.le.2) then
     if(xicur.eq.epsinv) then
        goto 800  ! treat as orbit loss (ion never inside plasma)
     end if
  end if
  !
  if((nbflg.eq.1).and.(ztslmg.gt.0.0_fp)) then
     if(nrip.eq.1) call riplos(ztslmg,iripl,ierr)
     if(ierr.ne.0) return
     if(iripl.eq.1) goto 820
  end if
  !%      if(nbplot.le.0) goto 7000
  !%      if(nbflg.eq.1) write(zbuf,7703)
  !%7703      format('timestep over$')
  !%      iplsty=iplsty-15
  !%      zxpos(1)=(jt-1)
  !%      zypos(1)=(jt-1)
  !%      call symbl(0)
  !%      call line(0)
  !%      call r8_cplot(zxpos,zypos)
  !%!  plot gyro of bad orbit
  !%      if((nbflg.eq.3).and.(ilim.eq.1)) then
  !%        zrmj=rmjion
  !%        call pcircl(zrmj,yion,rlar,zxpos,zypos,40)
  !%      end if
  !%      call rmrglb(zbuf(1),iplsty)
  !%      write(zbuf,7705)
  !%7705      format('orbit time:$')
  !%      call stplab(800,iplsty,zbuf(1),bitime(1))
  !%      call line(-1)
  !%      call tsend
  !%7000      continue
  !
  !  final exit
  !
  !  reset goose factor to 1.0 if failed to change this timestep
  !
  if((.not.nlorbo).and.(igoosc.eq.0)) goose=min(1.0001_fp,goose)
  !
  !  test for excessively high kinetic energy
  if((v2toev*vion*vion.gt.zmaxe.and.orbrzv_rf_option.ne.1).and. &
       (v2toev*vion*vion.gt.zmaxe.and. &
       (ntest_max.eq.0.or.abs(ntest_max).eq.10))) then
     zep=v2toev*vion*vion
     call errorlog('orbit: ptcl energy(ev) detected on exit', &
          zep)
     if(myidd.eq.0) then
        write(nonlin,9102) zep
     end if
9102 format(/ &
          ' ?orbit -- ',1pe11.4,' ev ptcl detected on exit.'/)
     ! when rf-option is on, just report the high energy
     if(orbrzv_rf_option .ne. 1) then
        if(ilim.eq.1.and.icheck.eq.1) then
           call errorlog(' orbit: ok as ptcl lost to limiter', &
                zep)
        else
           ierr=1
           return
        end if
     end if
  end if
  !
  !  mc list density, end of step
  !
  if(nbflg.eq.1) then       !nbflg=1 means an ion which has survived
     if((lcentr.le.ngc).and.(ngc.le.ledge)) then
        mc_dens(ngc,lsbeam)=mc_dens(ngc,lsbeam) + wght
     end if
  end if
  !
  return
  !
  !-----------------------------------------------------------------------
contains
  !
  subroutine orbtrace
    use iso_c_binding, only: fp => c_double
    call infolog('orbitrace: iorb = ',iorb)
    call infolog('itrak = ',itrak)
    call infolog('istep = ', istep)
    call infolog(' e(ev) =  ',  v2toev*vion*vion)
    call infolog(' vpll/v = ',  xksid)
    call infolog(' xion  = ',xion)
    call infolog(' th = ',th)
    call infolog(' e(ev) before orbit call: ',zepp)
    if(myidd.eq.0) then
       write(nonlin,6666) iorb,itrak,istep,v2toev*vion*vion, &
            xksid,xion,th,zepp
    end if
6666 format(' ------------------------- '/ &
         ' iorb# ',i7,' itrak# ',i3,' istep# ',i7/ &
         ' e(ev)= ',1pe12.5,'  vpll/v=',1pe12.5,'  xion,th = ', &
         2(1pe12.5)/ &
         ' e(ev) before orbit call: ',1pe12.5)

    return

  end subroutine orbtrace
  !---------------------------------------------------
  subroutine bounce_ini(jcall)
    use iso_c_binding, only: fp => c_double
    !
    !  initialize local sums and controls for goose re-evaluation and update.
    !  see also bounce_chk, which decides when to do the re-eval and update.
    !
    integer, intent(in) :: jcall
    !  jcall=1 -- 1st call of orbit track
    !
    !  set up goose, fp and cx acceleration schemes
    !
    if(jcall.eq.1) isumiz=0     ! ionization count (heavy ions)
    !
    !  clear timestep counters
    !
    istpor=0
    istpcx=0
    istpfp=0
    !
    !  set up expected orbit steps to cx and to fp
    inexcx=1+nbrandm_pll(idum)*cxskip
    inexfp=1+nbrandm_pll(idum)*max(1.0_fp,fpskip)
    !
    zdsdfb=0.0_fp  ! anom diffusvity sum (abdifu)
    zdsxvb=0.0_fp  ! anom velocity sum (abdifu also)
    !
    viplasa=0.0_fp  ! avg beam ion veloc., plasma frame
    ztimfp=0.0_fp
    !
    !  save current ptcl y (sign indicates above/below midplane)
    call bounce_ymp
    !
    !      goose factor set up in main code - particle retains same
    !      goose factor for at least one half orbit
    !
    !  initialize avg-j calculation
    if(jcall.eq.1) then
       xjav=real(ngc)
       zjtim=bitime(1)
       iflip=-1
       isuma=0
       zjsum=0.00_fp
       zdtime=0.0_fp
       !  compute and reinitialize avg-j calculation
    else
       iflip=-iflip
       if(iflip.eq.1) then
          xjav=zjsum/zdtime
          zjtim=bitime(1)
          isuma=isuma+1
          zjsum=0.00_fp
          zdtime=0.0_fp
       end if
    end if
    !
    !  set up time to charge exchange and clear c-x accumulators
    if(jcall.eq.1) then
       ttcx=-log(nbrandm_pll(idum))
       tscx=0.0_fp
    else
       !
       !  rev dmc june 1984 - bug fix - reduce time-to-cx (ttcx) by time-
       !  since-cx (tscx) and reset tscx and ztscxp to 0.0.  this avoids
       !  numeric problem due to large changes in tscx sum due to external
       !  cx followed by 10**-8 smaller changes when a fat banana untraps
       !  inside in very low neutral density region in hi density plasma
       !   (avoid tscx-ztscxp=0.0 due to difference small compared to values)
       !
       ttcx=ttcx-tscx
       tscx=0.0_fp
    end if
    ztscxp=0.0_fp
    !  zero sum of time spent between midplane crossings
    ztslmp=0.0_fp
    ztslmg=0.0_fp
    !  initialize counters and accumulators that run over a half-orbit
    dvg=0.0_fp
    dksgc=0.0_fp
    dksgd=0.0_fp
    !
    !  for heavy ions, the cx numeric controls are reused to control
    !  modeling of electron stripping; if isumiz.gt.0, there were such
    !  events this bounce (dmc jan 92)
    if(isumiz.gt.0) then
       cxskip=1.0_fp
       isumiz=0
    end if
    !
    !  look for orbits that are entirely beyond the plasma boundary...
    !
    if((jcall.eq.1).or.(ngoose.ge.2)) then
       zxion_min = xion
    end if
    !
    !  check energy & momentum error accumulation / half-bounce
    !
    zderrmin=epsinv
    zderrmax=-epsinv
    zderrsum=0.0_fp
    !
    zdprrmin=epsinv
    zdprrmax=-epsinv
    zdprrsum=0.0_fp
    !
    !  min, max b seen on orbit; max xion on orbit
    !
    zmaxbo = bimodb(1)
    zminbo = zmaxbo

    zx_bmax=xion
    zvpp=max(1.0_fp,vperp)
    zmu_bmax=(zvpp*zvpp)/zmaxbo
    zbrefl_bmax = max(zmaxbo,(vion*vion/zmu_bmax))
    !
    !  min, max x on orbit; <x> = (zminxo+zmaxxo)/2 will be used
    !    <x> is a proxy for <psi> in some documentation -- est. orbit
    !    averaged position
    !
    zmaxxo = xion
    zminxo = zmaxxo
    !
    !  for computation of <v_pll>, <e>
    !
    zdtsum = 0.0_fp
    zvpll_av = 0.0_fp
    !
    !  v*dt sum
    !
    zvdtsum = 0.0_fp
    !
  end subroutine bounce_ini
  !-----------------------------------------------------------------
  subroutine bounce_ymp

    zympl=ympx(ngcx)+xingcx*(ympx(ngcx+1)-ympx(ngcx))
    ztesm=yion-zympl
    if(abs(ztesm).lt.epslon) ztesm=epslon

  end subroutine bounce_ymp
  !-----------------------------------------------------------------
  logical function bounce_chk(itrig)
    integer, intent(out) :: itrig    ! reason for true value
    ! =1 always if bounce_chk is false
    ! =1 if bounce_chk is true because of midplane crossing
    ! =2 if bounce_chk is true because bounce time expired

    real(fp) :: zbmin,zbmax,zvtaub
    real(fp) :: zintrp
    integer :: iintrp

    !  this gives the number of full orbit bounce times
    !  without a midplane crossing, needed to trigger an
    !  early regoose calculation...
    real(fp), parameter :: zorbfac = 2.0_fp

    ztesp = ztesm

    call bounce_ymp

    itrig=1
    bounce_chk=.false.
    if((ztesm*ztesp).lt.0.0_fp) then
       bounce_chk=.true.
    else if(zvdtsum.gt.zorbfac*vtaubmin(ngc)) then
       zbmin=bminm(min(lep1,ngc+1))
       zbmax=bmaxm(min(lep1,ngc+1))
       zintrp=4*(zmaxbo-zbmin)/(zbmax-zbmin)
       if(zintrp.le.1.0_fp) then
          zvtaub=vtaub(1,ngc)
       else if(zintrp.ge.4.0_fp) then
          zvtaub=vtaub(4,ngc)
       else
          iintrp=zintrp
          zintrp=zintrp-iintrp
          zvtaub=vtaub(iintrp,ngc)+ &
               zintrp*(vtaub(iintrp+1,ngc)-vtaub(iintrp,ngc))
       end if
       if(zvdtsum.gt.zorbfac*zvtaub) then
          bounce_chk=.true.
          itrig=2
       end if
    end if

  end function bounce_chk
  !-----------------------------------------------------------------------
end subroutine orbit
