! written by nbigen.py
! *** do not edit *** (use nbigen.py instead)
module nbi_com
  use iso_c_binding, only: fp => c_double
  use physconst_mod, only: xpi=>pi, twopi, xe=>kb_eV
  use physconst_mod, only: xmu=>RMU0, degrad=>deg2rad
  use physconst_mod, only: zel, VC, zmp=>mp_g, zmtocm
  use physconst_mod, only: epslon, epsinv
  use physconst_mod, only: zero,half,one,two
  use nbi_dimensions
  implicit NONE
  SAVE
 
  real(fp), parameter :: zc=VC*zmtocm 
 
  !==============================================
  ! CONSTANTS...
 
  !==============================================
  ! INPUT items...
  ! block: sys
    ! basic system information for Monte Carlo code
 
  LOGICAL :: init_sys = .FALSE.
 
  character*32 :: runid = "NUBEAM"
        ! runid string -- alphanumeric string used as prefix for
        ! output files (e.g. state files, etc.).  Should generally
        ! be set to a runid string associated with the driver code.
  character*120 :: mission =  "standard"
        ! NUBEAM "mission" in presence of Plasma State data
        !   (each supported mission requires special code; at present
        !   only "standard" is supported).
        ! "standard" means:  assume responsibility for deposition
        !   and slowing down of neutral beam and fusion product
        !   fast ion species found in the Plasma State
  character*120 :: workpath =  " "
        ! work directory path -- for debug statefile...
        ! blank means current working directory
  integer :: nonlin = 6
        ! unit number for messages.
  integer :: nseed = 1
        ! random number seed; default not recommended!
  logical :: nlbout = .TRUE.
        ! .TRUE. (1) to turn on "lost particle" output data accumulation
  integer :: lunnbx = 0
        ! unit number for "lost particles" ascii output file
        ! (will be used if ref_namelist=T is set, see below).
        ! if lunnbx<=0 or lunnbx>99, an available I/O unit is
        ! selected automatically.
  integer :: lunres = 0
        ! auxilliary unit number for ascii file i/o -- supports
        ! debug restarts and restarts of "lunnbx" file.
        ! **** this LUN will be used **** even if .NOT.NLBOUT.
        ! if lunres<=0 or lunres>99, an available I/O unit is
        ! selected automatically.
  integer :: mrstrt = 1
        ! set to zero if there are *no* check-point / restart files
        ! (a property of the driver code).
        ! this changes the way the NB.OUT file is handled, and may
        ! have other effects in the future.
  integer :: nclass = 0
        ! set =1 to activate beam-by-beam full:half:third
        ! energy species resolved profiles of n, Pbe, Pbi,
        ! <v.B>, <(vcrit**3/v**3)v.B> (lab frame v.B) for NCLASS
        ! fast ion interface -- see nclass output block.
        ! *** warning *** dmc Dec 2003 -- definition and
        ! construction of nclass=1 outputs is not yet complete ***
  logical :: only_io = .FALSE.
        ! set .TRUE. (1) when calling NUBEAM from a serial driver that
        ! will use a separate MPI-parallel process to compute the
        ! fast ion evolution.  In this mode, only the input and
        ! output structures are allocated-- NOT the large internal
        ! computational structures.  The following routines are
        ! NOT available in this mode:  NBI_INTERP_PROFILES,
        ! NBSTART, DEPALL, ORBALL, FI_FINISH; instead the timestep
        ! advance sequence is:
        !    CALL NUBEAM_WRITE_INPUT(...)
        !     (invoke local procedure for calling parallel NUBEAM
        !     as separate process)
        !    CALL NUBEAM_READ_OUTPUT(...)
        ! these calls are more fully explained elsewhere.
  logical :: ref_namelist = .FALSE.
        ! NOTE: this is a legacy feature -- mothballed!
        ! set .TRUE. (1) to write a reference namelist
  logical :: quasi_check = .TRUE.
        ! .TRUE. to check quasineutrality of density input profiles
        ! (can set .FALSE. to disable test, but this can be dangerous!)
  logical :: nltest_output = .TRUE.
        ! Set .FALSE. (0) to suppress NBSTART output of interim state
        ! data for nubeam_test time step debug rerun capability
        ! (suppress e.g. to optimize NUBEAM mpi scaling).Default - .TRUE.
        ! debug_mpi= nltest_output for
        ! MPI job debugging
        ! debug_mpi   =.FALSE. for serial mode
        ! see below
        ! in TRANSP  analog of this switch is nlnbi_test_out
  integer :: BLK_MPI_BCAST_INT =  2000000
        ! default value for size of a block of integer variables to be used
        ! for MPI communications (bytes)
  integer :: BLK_MPI_BCAST_R8 =   30000000
        ! default value for size of a block of real variables to be used
        ! for MPI communications (bytes)
 
  ! block: times
 
  LOGICAL :: init_times = .FALSE.
 
  real(fp) :: tbm1 = 0.0_fp
        ! start time of current timestep  -- seconds
  real(fp) :: tbm2 = 0.0_fp
        ! stop time of current timestep  -- seconds
        ! the start time of timestep j+1 should match
        ! the stop time of timestep j (else a warning will be issued).
 
  ! block: grid
    ! basic grid information
 
  LOGICAL :: init_grid = .FALSE.
 
  integer :: nzones = 0
        ! number of radial zones (equally spaced in sqrt(tor.flux))
        ! this determines the grid for calculation of 1d outputs
        ! such as Pbe and Pbi; must be set to a number .le. mj-3
        ! where "mj" is a user supplied array dimension.
        ! NOTE: not included in Plasma_State_Input, because this
        ! value is used to **set** the plasma state rho_nbi grid.
  integer :: nzone_fb = 0
        ! number of radial zones in 4d fast ion distributions
        ! (r,theta,E,vpll/v).  nzones should be an integer multiple
        ! of nzone_fb; if not, nzone_fb is increased until this
        ! condition is satisfied, or if this is not possible,
        ! nzone_fb=nzones is imposed.  Note: if id_mcgrid is set,
        ! nzone_fb is reset to match the 2d GRID in xplasma.
  logical :: nlsym2b =  .FALSE.
        ! .TRUE. (1) to treat fast ion dist. fcn. as up-down symmetric.
        ! Note: if id_mcgrid is set, this is overridden by the 2d
        ! "irregular" grid in xplasma.
  integer :: nth0 =  2
        ! number of poloidal zones subtending [0:pi] in the first
        ! radial zone row adjacent to the magnetic axis.
        ! Note: if id_mcgrid is set, this is overridden by the 2d
        ! "irregular" grid in xplasma.
  real(fp) :: xbmbnd = 1.3
        ! boundary of xi-space
  integer :: id_mcgrid = 0
        ! xplasma id of 2d irregular "Monte Carlo" spatial grid.
        ! **must be supplied each time and shape must not change**
        ! if this is NOT set, then, nzone_fb, nlsym2b, and nth0
        ! inputs are used to set up the irregular grid.
  integer :: nthets = 64
        ! number of poloidal angle grid points in [J] and [1/B] 2d
        ! splines (cf nubeam/magcf.for)
  real(fp) :: ebdmax = 0.0_fp
        ! maximum energy in distribution function for beam ion species
        ! the maximum for fusion product species is set by the fusion
        ! product source energy (e.g. 3.5e6 eV for DT alphas-- adjusted
        ! upward to 5.0e6 allow for possible energy diffusion).
        ! **energy in eV**
  integer :: xplasma_in_memory = 0
        ! if >0, assume "old_xplasma" shared xplasma
        ! object has been prepositioned by caller to receive NUBEAM
        ! output; if =0, just use "old_xplasma" object; if <0, save
        ! old xplasma object for later restoration.
        !   ** most users should never touch this **
  integer :: nznbma = 50
        ! number of pitch angle zones for distribution functions
        ! this is the same for all NUBEAM fast ion species
  integer :: nznbme = 100
        ! number of energy zones for beam ion distribution functions
        ! note for fusion product species the number of energy zones
        ! is always equal to the maximum, mimxbe, regardless of the
        ! nznbme specification for beam species.
  integer :: nznbma0 = 10
        ! number of pitch angle zones for neutrals
  integer :: nznbme0 = 50
        ! number of energy zones for neutrals
  integer :: ngyro = 4
        ! NUMBER OF poloidal AIMING ANGLE ZONES, BEAM F-DIST. CALC.
        ! atan2(V_Z,V_R)-atan2(B_Z,B_R)
        ! must be > 0
        ! for ngyro=4
        ! zone 1 -- [-pi/4,pi/4]         always means perpendicular to the flux surface and outward
        ! zone 2 -- [pi/4,3/4pi]  is tangential to the flux surface with an inward FLR displacement to GC
        ! zone 4 -- [-pi/4,-3/4pi] is tangential to the flux surface with an outward FLR displacement to GC
        ! zone 3 -- [3pi/4,-3/4pi] is perpendicular and inward
 
  ! block: beams
    ! here *both* the beams and the thermal plasma species
    ! are defined together, because, for each beam specie
    ! there must be a thermal plasma specie for it to
    ! thermalize into.
    !
    ! except where noted, actual values must be given...
    ! in the course of a run, beams and species can be added
    ! but once added can never taken away.
 
  LOGICAL :: init_beams = .FALSE.
 
  integer :: ngmax = 0
        ! the number of thermal plasma species
  real(fp), dimension(:), allocatable :: backz
        ! the Z of each thermal specie (atomic no.)
  real(fp), dimension(:), allocatable :: aplasm
        ! the A of each thermal specie (AMU)
  integer :: nbeam = 0
        ! the number of distinct neutral beam sources
  real(fp), dimension(:), allocatable :: xzbeama
        ! Z of injected ion specie (atomic no.)
  real(fp), dimension(:), allocatable :: abeama
        ! A of injected ion specie (AMU)
  logical, dimension(:), allocatable :: nlco
        ! .TRUE. (1) if injected "co", i.e. in direction
        ! of the toroidal plasma current.
  integer, dimension(:), allocatable :: ntrace
        ! non-zero value means this beam (j) is a
        ! "trace element" of the beam indicated by
        ! the value of ntrace(j); it inherits the
        ! geometry of that beam and the interpretation
        ! of energy fractions is affected.
        !   (of course, usually, the default ntrace=0 is OK)
  real(fp), dimension(:), allocatable :: rtcena
        ! tangency major radius, cm, of centerline of beam
  real(fp), dimension(:), allocatable :: xlbtna
        ! distance, cm, from center of beam source grid to tangency point
  real(fp), dimension(:), allocatable :: xybsca
        ! elevation of center of beam source grid above vacuum vessel
        ! midplane "Z=0". (negative if located below midplane). cm.
        ! default value, 0.0_fp, means source lies on the midplane.
        ! in an axisymmetric system, rtcena, xlbtna, xybsca suffice
        ! to specify the source grid position.
  real(fp), dimension(:), allocatable :: xbzeta
        ! toroidal angle of beam source in an (R,zeta,Z)
        ! right handed coordinate system (degrees) (specification not
        ! required except to use "tube" diagnostic simulation software).
  integer, dimension(:), allocatable :: nbshapa
        ! ion source grid shape: 1=rectangular, 2=circular
  real(fp), dimension(:), allocatable :: bmwidra
        ! source grid half width (or radius if circular), cm.
  real(fp), dimension(:), allocatable :: bmwidza
        ! source grid half height (ignored if circular), cm.
        ! nbshapa,bmwidra,bmwidza specify the source shape and size.
  real(fp), dimension(:), allocatable :: divra
        ! horizontal divergence of beam (radians)
        !  (for more info see http://w3.pppl.gov/transp --
        !   TRANSP HELP -- go to beam controls, look for "DIVR")
  real(fp), dimension(:), allocatable :: divza
        ! vertical divergence of beam (radians)
  real(fp), dimension(:), allocatable :: foclra
        ! horizontal focal length of beam (cm)
  real(fp), dimension(:), allocatable :: foclza
        ! vertical focal length of beam (cm)
  integer, dimension(:), allocatable :: nbapsha
        ! shape of main aperture to vacuum vessel (1:rectangle, 2:circle)
        ! aperture cross section is presumed normal to beam centerline
        ! this aperture must be specified
  real(fp), dimension(:), allocatable :: rapedga
        ! aperture half-width (or radius if circular), cm
  real(fp), dimension(:), allocatable :: xzpedga
        ! aperture half-height (ignored if circular), cm
  real(fp), dimension(:), allocatable :: xlbapa
        ! distance, source grid to aperture, cm
  real(fp), dimension(:), allocatable :: xybapa
        ! beam centerline at aperture:  elevation above/below vacuum
        ! vessel midplane, cm
  real(fp), dimension(:), allocatable :: xrapoffa
        ! aperture horizontal offset relative to beam centerline, cm
        ! positive value: aperture shifted outward in major radius
        !   relative to beam centerline
  real(fp), dimension(:), allocatable :: xzapoffa
        ! aperture vertical offset relative to beam centerline, cm
        ! positive value: aperture shifted vertically upward
        !   relative to beam centerline
  integer, dimension(:), allocatable :: nbapsh2
        ! shape of optional 2nd aperture (1:rectangle, 2:circle)
  real(fp), dimension(:), allocatable :: rapedg2
        ! 2nd aperture half-width (or radius if circular), cm
        ! value of zero means: no 2nd aperture
  real(fp), dimension(:), allocatable :: xzpedg2
        ! 2nd aperture half-height (ignored if circular), cm
        ! value of zero means: no 2nd aperture
  real(fp), dimension(:), allocatable :: xlbapa2
        ! beam centerline distance, source grid to 2nd aperture, cm
        ! value of zero means: no 2nd aperture
  real(fp), dimension(:), allocatable :: xrapoff2
        ! 2nd aperture horizontal offset relative to beam centerline, cm
        ! positive value: aperture shifted outward in major radius
        !   relative to beam centerline
  real(fp), dimension(:), allocatable :: xzapoff2
        ! 2nd aperture vertical offset relative to beam centerline, cm
        ! positive value: aperture shifted vertically upward
        !   relative to beam centerline
 
  ! block: impurity
    ! Zx & Ax of impurity species
    ! can also use a single "model impurity" by
    ! giving Zx(x) and Ax(x) as average Z & A of
    ! impurities as profile input
    !
    ! in the course of a run, impurity species can be added
    ! but once added can never be taken away
 
  LOGICAL :: init_impurity = .FALSE.
 
  integer :: nrhix = 0
        ! number of impurities
  real(fp), dimension(:), allocatable :: xzimpx
        ! Z of (multiple) impurities (atomic no.)
  real(fp), dimension(:), allocatable :: aimpx
        ! A of (multiple) impurities (AMU)
 
  ! block: minority
    ! RF minority speciess
    ! in the course of a run, RF species can be added
    ! but once added can never be taken away
 
  LOGICAL :: init_minority = .FALSE.
 
  integer :: nmini = 0
        ! number of RF minority species
  real(fp), dimension(:), allocatable :: xzmini
        ! Z of RF minority specie (atomic no.)
  real(fp), dimension(:), allocatable :: amini
        ! A of RF minority specie (AMU)
  integer :: nthrf = 1
        ! = 1 if minority ions are currently thermal
 
  ! block: powers
    ! beam powers, voltages, and energy fractions
    ! energy fractions:  for H-isotope beams only
    !   ffulla(j) = fraction of beam current at full energy
    !   fhalfa(j) = fraction at half energy
    !   1-ffulla(j)-fhalfa(j) = fraction at 1/3 energy
    ! (for non-H beams we assume ffulla(j)=1, fhalfa(j)=0)
    !
    ! if f1=ffulla(j), f2=fhalfa(j), f3=1-f1-f2, and
    ! J = the beam current for beam (j), then, the
    ! relation
    !   einja(j)*J*(f1 + f2/2 +f3/3) = pinja(j)
    ! is satisfied.
    !
    ! generally the powers, voltages and (perhaps) the
    ! energy fractions should be updated every timestep.
 
  LOGICAL :: init_powers = .FALSE.
 
  real(fp), dimension(:), allocatable :: einjmin
        ! minimum energy ever injected (eV)
        ! use this to prevent problems due to Einj->0 as can
        ! occur when interpolating waveform data.
        ! TRANSP uses 1/2*<Einj> where <Einj> is the power
        ! weighted average injection energy, for each beam.
        ! note this can only be set once at the start of a
        ! run; if it is defaulted, (1/4)*max(Einj) is used
        ! as the lower limit, where max(Einj) is the max
        ! Einj for the beam, seen up until the current time.
        !   (see nbi_nbsums.for)
  real(fp), dimension(:), allocatable :: einja
        ! beam voltages (eV)
  real(fp), dimension(:), allocatable :: pinja
        ! beam powers (watts)
  real(fp), dimension(:), allocatable :: ffulla
        ! full energy current fraction (dimensionless)
        ! (a realistic value should be provided)
  real(fp), dimension(:), allocatable :: fhalfa
        ! half energy current fraction (dimensionless)
        ! (a realistic value should be provided)
  real(fp), dimension(:), allocatable :: orbrzv_toric_Prftot0
        ! total power calculated by toric in the previous step
  real(fp) :: orbrzv_toric_frnm = 1.0_fp
        ! toric RF-field strength renormalizer
 
  ! block: fusion
    ! fusion products
 
  LOGICAL :: init_fusion = .FALSE.
 
  logical :: nlfhe3 = .FALSE.
        ! .TRUE. (1) for He3 fusion product modeling
  logical :: nlfhe4 = .FALSE.
        ! .TRUE. (1) for He4 (alpha) fusion product modeling
  logical :: nlfst = .FALSE.
        ! .TRUE. (1) for DD triton fusion product modeling
  logical :: nlfsp = .FALSE.
        ! .TRUE. (1) for DD proton fusion product modeling
        ! (no recommended-- proton velocities exceed limiting
        ! assumptions in collision operator)
  logical :: nlfsp14 = .FALSE.
        ! .TRUE. (1) for DHe3 proton fusion product modeling
        ! (no recommended-- proton velocities exceed limiting
        ! assumptions in collision operator)
  logical :: nlusf3 = .FALSE.
        ! .TRUE. (1) to set He3 fusion product source magnitude from
        ! data instead of from past timestep fusion rates
  real(fp) :: xdatsf3 = 0.0_fp
        ! source magnitude (n/sec) used to renormalize He3 fusion
        ! product source profile from previous timestep
        ! only applies if nlusf3.and.nlfhe3 both .TRUE.
  logical :: nlusfa = .FALSE.
        ! .TRUE. (1) to set alpha fusion product source magnitude from
        ! data instead of from past timestep fusion rates
  real(fp) :: xdatsfa = 0.0_fp
        ! source magnitude (n/sec) used to renormalize alpha fusion
        ! product source profile from previous timestep
        ! only applies if nlusfa.and.nlfhe4 both .TRUE.
  logical :: nlusft = .FALSE.
        ! .TRUE. (1) to set triton fusion product source magnitude from
        ! data instead of from past timestep fusion rates
  real(fp) :: xdatsft = 0.0_fp
        ! source magnitude (n/sec) used to renormalize triton fusion
        ! product source profile from previous timestep
        ! only applies if nlusft.and.nlfst both .TRUE.
  logical :: nlusfp = .FALSE.
        ! .TRUE. (1) to set proton fusion product source magnitude from
        ! data instead of from past timestep fusion rates
  real(fp) :: xdatsfp = 0.0_fp
        ! source magnitude (n/sec) used to renormalize proton fusion
        ! product source profile from previous timestep
        ! only applies if nlusfp.and.nlfsp both .TRUE.
  logical :: nlusfp14 = .FALSE.
        ! .TRUE. (1) to set proton fusion product source magnitude from
        ! data instead of from past timestep fusion rates
  real(fp) :: xdatsp3 = 0.0_fp
        ! source magnitude (n/sec) used to renormalize proton fusion
        ! product source profile from previous timestep
        ! only applies if nlusfp.and.nlfsp both .TRUE.
  real(fp) :: plfhe3 = 100.0_fp
        ! (watts) threshhold power above which full Monte Carlo
        ! statistics (NPTCLF particles) are assigned to model the
        ! slowing down of the He3 fusion product ion population
  real(fp) :: plfhe4 = 10000.0_fp
        ! (watts) threshhold power above which full Monte Carlo
        ! statistics (NPTCLF particles) are assigned to model the
        ! slowing down of the He4 (alpha) fusion product ion population
  real(fp) :: plfst = 100.0_fp
        ! (watts) threshhold power above which full Monte Carlo
        ! statistics (NPTCLF particles) are assigned to model the
        ! slowing down of the triton fusion product ion population
  real(fp) :: plfsp = 100.0_fp
        ! (watts) threshhold power above which full Monte Carlo
        ! statistics (NPTCLF particles) are assigned to model the
        ! slowing down of the proton fusion product ion population
  real(fp) :: plfsp14 = 100.0_fp
        ! (watts) threshhold power above which full Monte Carlo
        ! statistics (NPTCLF particles) are assigned to model the
        ! slowing down of the proton fusion product ion population
  logical :: nlfatom = .TRUE.
        ! .FALSE. (0) to suppress atomic physics calculations on
        ! slowing down fusion product ions (recommendation:
        ! leave this switch alone).
 
  ! block: fpp
    ! Monte Carlo code support for Fokker Planck model
 
  LOGICAL :: init_fpp = .FALSE.
 
  logical :: nlbfpp = .FALSE.
        ! set .TRUE. (1) if beam slowing down is to be modeled using
        ! a fokker planck code, in which case the Monte Carlo code
        ! will be used to calculate deposition only
        ! NOTE: this feature NOT implemented DMC May 2007
 
  ! block: num
    ! numerical controls
 
  LOGICAL :: init_num = .FALSE.
 
  real(fp) :: cxsplt = 2.0_fp
        ! charge exchange model: splitting factor, tracks per event.
  real(fp) :: dxbsmoo = 0.05_fp
        ! smoothing parameter: half-width of triangular hat
        ! function convolution smoothing, in x=sqrt(phi/philim)
        ! applied to all output profiles
        ! NOTE: if dxbsmoo.gt.0 a spline method will be applied
        !   to integrated profiles; if dxbsmoo.le.0 a piecewise
        !   linear method will be applied.  Thus:
        !      dxbsmoo .le. 0 ==> PBI is a step function
        !      dxbsmoo .gt. 0 ==> PBI is piecewise quadratic, C1.
        !   which follows because the integrated profile is
        !   differenced to get the local value.
  real(fp) :: dxbsm_nc = 0.0_fp
        ! (Not currently active -- left in namelist for backwards
        ! compatibility and possible later use).
        ! original intent:
        ! dxbsm_nc smoothing parameter (like dxbsmoo) for nclass
        ! profile outputs only.  Actual smoothing for these
        ! profiles will be based on max(abs(dxbsmoo),abs(dxbsm_nc));
        ! spline or linear treatment still controlled by sign of
        ! dxbsmoo; effect is optional extra smoothing for nclass
        ! outputs. NO more active MG dec.09,2008
  integer :: nptcls = 1000
        ! number of Monte Carlo particles per (normal) beam ion specie
        ! affects statistics of slowing down calculation
  integer :: nptclf = 1000
        ! number of Monte Carlo particles per fusion product specie
        ! affects statistics of slowing down calculation
  integer :: nptclh = 500
        ! number of Monte Carlo particles per heavy beam ion specie
        ! affects statistics of slowing down / stripping calculation
  integer :: NFAST_MAP = 1
        ! in NUBEAM only: NFAST_MAP=1 - use bilinear maps to evaluate
        ! R,Z coordinates and 2x2 Jacobian,
        ! if NFAST=0 - use xplasma spline maps with tol=1.e-6 (~1.5 slower)
  real(fp) :: xpand_nptcl = 3.0_fp
        ! If xpand_nptcl.gt.0.5 it will be a limit
        ! of the total growth of the particle list
        ! after this limit (> (1.0+xpand_nptcl)*NPTCLS )
        ! revert to Russian Roulette.
        ! For xpand_nptcl.le.0.5          it will be a fraction by which
        ! nptcl{*} can expand during
        ! deposition (see below) to avoid Russian Roulette (min val 0.05_fp)
        ! (max val 0.5_fp).
        ! (max(nptcls,nptclf,nptclh,nptcl_max)+1000)*1.1_fp
        ! *(1+xpand_nptcl+frac_orbrr)
  real(fp) :: frac_depmax = 0.25_fp
        ! max fraction of nptcl{*} # of MC ptcls that can be
        ! deposited in a single step (min accepted val 0.1_fp)
  real(fp) :: frac_dep_lim = 1.5_fp
        ! nptcl{*} # of MC ptcls that can be
        ! deposited in a single step
        ! is max(1,nptcls*min(1,delt/<tau_sl>*frac_dep_lim))
  real(fp) :: frac_depmin = 0.01_fp
        ! min fraction of nptcl{*}, min # of MC ptcls to
        ! deposit in a single step (btw 0.0_fp and 0.1_fp)
  real(fp) :: frac_orbrr = 0.10_fp
        ! fraction of population needing to be "underweight"
        ! before application of Russian Roulette in orball routine
        ! accepted values betwe 0.05_fp and 0.5_fp.
  integer :: nptcl_max = 0
        ! max number of Monte Carlo particles of any type that will
        ! ever be requested; max(nptclh,nptclf,nptcls,nptcl_max) is
        ! used.
  logical :: nltrk_dep0 = .FALSE.
        ! .TRUE. (1) to save details of ndep0 beam neutral tracks
  logical :: NLDEP0_GATHER =  .FALSE.
        ! IF .TRUE. fast ion birth source data (in NUBEAM_COMP_EXEC)
        ! will be collected for all time steps or only for the last one
        ! if NLDEP0_GATHER = .false.
        ! The output files <RINID>_nubeam_birth.cdf<nn> will be kept in
        ! main directory from where program was started (in NUBEAM_COMP_EXEC).
  integer :: ndep0 = 500
        ! minimum number of deposition tracks per beam ion specie--
        ! also, the number of tracks gathered, if deposition
        ! distribution data is saved for output, i.e. if
        ! nltrk_dep0 is .TRUE.
        ! Regardless of nltrk_dep0, This control
        ! affects statistics of deposition calculations (e.g.
        ! cold electron source, CX neutral sources, etc.).
        ! In an MPI run this value may be adjusted to be an exact
        ! multiple of the number of processors used.
        ! When deposition track data is gathered, there is an
        ! option to restrict the gathering to a single specific
        ! beam injector, or even a specific energy fraction for
        ! single injector.  (See ndep_set_beam and ndep_set_ien
        ! controls).  But, if this is done, the beam current
        ! fraction (call it f) of the selected beam injector can
        ! come into plan in determining the track data gathering.
        ! If f < 10%, the tracks gathered will be
        !    ndep0_use = min(ndep0,f*ndep0_max)
        ! if f < 1%, the restriction to the selected beam subset
        ! is aborted, and data, ndep0 tracks, are gathered on
        ! all beam lines.  f is output as frac_cur_dep0.
  integer :: NDEP0_MAX =    100000
        ! approximate max # of tracks to compute, when
        ! track data is gathered for only a single beam or
        ! single beam energy fraction-- used if the selected
        ! beam current fraction f is < 10%.
  integer :: ndep_set_beam = 0
        ! beam source # to which deposition data will be gathered
        ! (must be in range 0 to mib).  Default is 0:
        ! 0 - means gather data for all available beame source
        ! if ndep_set_beam set 0, then ndep_set_ien must be set 0.
        ! If beam source# was incorrectly set, then a warning is
        ! printed, and ndep_set_beam=0 is enforced.
        ! If there is no particle current for requested beam or
        ! this current is less than 1% of total ptcl. current then
        ! collect track data for all beam sources and energy fractions
        ! (i.e. revert to ndep_set_beam=0 mode).
        ! If particle current for requested beam is less than 10%
        ! of total ptcl. current - the number of tracks collected
        ! is adjusted to min(NDEP0,f*NDEP0_MAX), where f is the beam
        ! current fraction for the selected beam, relative to the
        ! total summed over all beams.
  integer :: ndep_set_ien = 0
        ! energy fractions to be used when gathering deposition data
        ! for a given beam source (ndep_set_beam .gt. 0).
        ! 1 - full energy fraction
        ! 2 - half energy fraction
        ! 3 - one third energy
        ! 0 - for all energy fraction
        ! if energy fraction# was incorrect set data is collected
        ! for all energy fractions.
        ! For beams without energy fractions (i.e. Helium neutral
        ! beams or high energy negative ion based H/D/T beams),
        ! this control is ignored.
  integer :: nbbcal = 1000000
        ! number of interactions for [fast ion] - [fast ion]
        ! fusion rate Monte Carlo integrals-- affects statistics
        ! of this component of fusion product and neutron source
        ! calculation
  integer :: ncx0 = 100
        ! the total number of CX tracks, independent control of CX statistics
        ! if ncx0 = 0, then it will be no fbm_cx gathered and CX statistics will not be
        ! improved
  integer :: nper_cx = 4
        ! the  number of gyro phase angles between [0,2*pi] will be used for CX statistics improvement
        ! if ncx0 > 0 , nper_cx=max(1,nper_cx)
  integer :: LEV_NBIDEP =  1
        ! beam "deposision" model level:
        ! 1 -- legacy of PREACT atomic physics
        ! 2 -- ADAS atomic physics for H and He beams.
        !      He beam -- only ground state model available
        ! 3 -- PREACT ground state atomic physics
        !      ADAS interaction with impurities
        ! 4 -- ADAS ground state atomic physics
        !      PREACT interaction with impurities
  integer :: LEV_TTNUK =  0
        ! nuclear data for T(T,2n)4He  fusion reaction
        ! 0 -- legacy of PREACT nuclear physics for T(T,2n)4He
        ! 1 -- DPC fit (David.P.Coster fit added on 07/2018)
  integer :: LEVMOD_HALO =  0
        ! 0 -- driver program calls FRANTIC (default)
        ! 1 -- FRANTIC will be called in NUBEAM
        ! 2 -- 3d halo model will be called in NUBEAM
        !      for 3d halo model NSPLT_GEO and NSPLT_KIN have to be set
        !      or default values will be used
        !      as 3d halo model is a 'beam in box' then
        !      NBSBOX, NBEBOX, NXBOX, NYBOX, NLBOX, NDEPBOX etc. have to be set
  integer :: NSPLT =  100
        ! activated in case LEVMOD_HALO=2
        ! max allowed spliting along each neutral track (3D halo model),
        ! NSPLT must be >= NSPLT_GEO* NSPLT_KIN
        ! if NSPLT < 0 then it  maybe use to gather 3d density 0 generation of neutrals, only
  integer :: NSPLT_GEO =  1
        ! activated in case LEVMOD_HALO=2, and NSPLT > 0
        ! max allowed spliting along each neutral track (3D halo model),
        ! max allowed geometrical spliting along each neutral track (3D halo model)
        !
  integer :: NSPLT_KIN =  1
        ! activated in case LEVMOD_HALO=2, and NSPLT > 0
        ! max allowed spliting along each neutral track (3D halo model),
        ! max allowed kinetic spliting at each CX launch point; # of velocities
        ! that are sampled from eligible ion populations (3D halo model)
  integer :: NSPLT_GEN =  1
        ! # of halo generation to be splitted at the
        ! rate NSPLT_KIN * NSPLT_GEO, and for higher generation will be no additional
        ! splitting made.  If NSPLT_GEN = 0, all halo generation will be splitted at the
        ! rate NSPLT_KIN * NSPLT_GEO. Activated in case LEVMOD_HALO=2, and NSPLT > 0
  integer :: nmimp =  1
        ! 1 -> using multiple impurities
        ! 0 -> using single tokamakium impurity
  real(fp) :: AIMP =    16._fp
        ! ATOMIC WEIGHT OF IMPURITY, for tokamakium model only
  real(fp) :: XZIMP =    8._fp
        ! IMPURITY ATOMIC NUMBER, for tokamakium model only
  integer :: NTEST_MAX =  0
        ! set 1 to simulate bi-Maxwellian FI distribusion;
        !  -1 to skip re-normalization of electric field from TORIC
        ! in RF kick off operator in NUBEAM for bi-Maxwellian FI distribusion;
        ! -10 to skip re-normalization of electric field from TORIC
        ! in RF kick off operator in NUBEAM for realistic FI distribusion
        ! this switch may be used for standalone nubeam_comp_exec program
        ! and can be modified through nubeam_comp_exec namelist
        ! 10 same as default but extra output files will be saved to be used in TORIC:
        ! some variables on sqrt(norm poloidal flux)
  real(fp) :: BI_MAX_PL =  1._fp
        ! factor for an E_PLL component in bi-Maxwellian FI distribusion
  real(fp) :: BI_MAX_PR =  1._fp
        ! factor for an E_Perp component in bi-Maxwellian FI distribusion
  real(fp) :: wghta = 1.0_fp
        ! weight profile adjustment.  Larger values of WGHTA result
        ! in more particles with less weight per particle in the core
        ! and fewer particles with more weight per particle at the
        ! edge.  (dimensionless).
        !   wghta.ge.1.0 enforced; wghta.le.nzones enforced.
        !!! Larger WGHTA -> better statistics at core, worse at edge
        !!! Smaller WGHTA -> better statistics at edge, worse at core
  real(fp) :: goocon = 10.0_fp
        ! "goose" control -- strong affect on runtime!
        ! #bounces computed per min(tau(slow),tau(P.A.scatter),tau(cx)).
        ! magnitude of constraint is quadratic for diffusive processes.
        ! see also NGOOCON controls; if any of the NGOOCON controls have
        ! greater than their default values, an upper limit of 1.0e-4 is
        ! imposed on DT_ACC; otherwise upper limit of 1.0e-3 is allowed.
  integer :: NGOOCON_VPVBIN =  0
        ! min #bounces per FBM pitch zone width P.A. scatter
        ! setting NGOOCON_VPVBIN constrains tau(bounce)/tau(P.A.scatter) by
        ! a factor of 4/(NGOOCON_VPVBIN**2*NZNBMA**2)
        ! a limit that scales as 1/(nznbma**2) on goosing is imposed if .gt.0
        ! Caution: setting this may slow the code considerably.
  integer :: NGOOCON_EBIN =  1
        ! min #bounces per FBM energy zone width slowing down
        ! additional linear constraint on tau(bounce)/tau(slowing-down)
        ! a limit that scales as 1/nznbme on goosing is imposed
  integer :: NGOOCON_RBIN =  1
        ! min #bounces per zone width radial transport
        ! linear constraint on tau(bounce)/tau(radial-flow)
        ! quadratic constraint on tau(bounce)/tau(radial-diffusion)
        ! limits that scale as 1/(nzone_fb) and 1/(nzone_fb**2) are imposed.
  real(fp) :: cxpcon = 20.0_fp
        ! adjustment on frequency of calls to charge exchange
        ! operator during slowing down.  Higher number --> more calls
        ! default value generally OK.
  real(fp) :: fppcon = 8.0_fp
        ! adjustment on frequency of calls to collision operator
        ! per orbit half-bounce during slowing down.  Higher number
        ! --> more calls; default value generally OK.
        ! if nclass is set, the number is increased to at least
        ! 20.0_fp
  real(fp) :: fdtnxy = 0.5_fp
        ! reduction factor for orbit timesteps near the axis
        ! where a pseudo x-y space instead of a polar space is used
        ! to avoid effects of the axial singularity.  Default OK.
  integer :: nchdvp = 3
        ! orbit integrator options:
        !  =1: enforce exact conservation of angular momentum
        !  =2: enforce exact conservation of energy
        !  =3: (recommended) obtain good but not exact conservation
        !      of both quantities.
  real(fp) :: dtn = 0.0_fp
        ! ***orbit timestep adjustment parameter*** if left to zero
        ! an automatic timestep will be selected and adjusted based
        ! on orbit integration error analysis.  If set, constraints
        ! will be applied, and the time step will still be adjusted
        ! based on the error analysis.  General advice: leave at
        ! zero.  (dtn, seconds).
  real(fp) :: dt_acc = 1.0e-4
        ! orbit accuracy criterion -- orbit timesteps adjusted so
        ! that error delta(E)/Einj (and a similar measure for momentum)
        ! is less than dt_acc per half bounce; another criterion
        ! is that there be at least 10 steps per half bounce so
        ! that goosing & timestep adjustment algorithms have
        ! enough data to work.
        !    dt_acc<=0.0_fp => no adjustment (old constant dtn used).
        ! if dt_acc > 0, 1.0e-6 <= dt_acc <= 1.0e-3 is enforced.
        ! an upper limit of 1.0e-3 is enforced unless GOOCON controls
        ! are increased from their default values in which case the
        ! upper limit is 1.0e-4
  integer :: orbrzv_option = 0
        ! option to switch on & off ORBRZ(RZ-integrator) and ORBG
  integer :: orbrzv_rf_option = 0
        ! option to switch on & off Monte-Carlo RF calculation
  logical :: nlmcrf_on = .false.
  integer :: orbrzv_nr = 160
  integer :: orbrzv_nz = 160
        ! grid resolution for RZ-interpolation in ORBRZ
  real(fp) :: orbrzv_zzerr_con = 1.0e-4_fp
        ! relative error control criterion
        ! 1.0E-4 does very expeditious job
 
  ! block: atomic
    ! atomic physics controls
    ! and controls affecting the beam deposition model
 
  LOGICAL :: init_atomic = .FALSE.
 
  integer :: nsmax = 21
        ! (maximum) number of thermal neutral sources
        !   (used to be fixed at 3*ngmax due to TRANSP legacy).
  integer, dimension(:), allocatable :: khalo
        ! indicate which neutral sources are due to the
        ! beam halo:  1 means halo, 0 means edge source not halo,
        ! -1 means "auto-detect" based on assumption that neutral
        ! densities due to edge sources will be strongly peaked.
        ! towards the edge
  real(fp) :: dn0out = 1.0e+11_fp
        ! neutral gas density (/cm3) to assume beyond plasma boundary
        ! (drives "external charge exchange loss")
  integer :: nmsigx = 2
        ! control option: select cross section set for impurity
        ! stopping of fast neutrals (deposition and cx recapture)
        ! 1: Olson (very ancient), H beams only
        !    OLSON ET AL., P.R.L. 41,NUM 3,P. 163.
        ! 2: cross section set for impurity will be defined by
        ! LEV_NBIDEP and NSIGEXC settings
        ! 2 use to be: Standard set
        !    Currently:
        !    R.A. Phaneuf, R.K. Janev, M.S. Pindzola,
        !    Atomic Data for Fusion, Vol. 5,
        !    Collisions of Carbon and Oxygen Ions with
        !    Electrons, H, H2, and He,
        !    ORNL-6090/V5 (5th volume of ORNL-6086) (Feb. 1987).
  logical :: nlminsv = .TRUE.
        ! .TRUE. (1) to loop over all impurities when calculating
        ! impurity stopping cross section; .FALSE. to use an
        ! "average Z" single impurity simplifying approximation.
  logical :: nlebei = .TRUE.
        ! .TRUE. (1) to include beam energy in approximate
        ! expression for <sigma*v> for electron impact ionization
  integer :: nsigexc = 0
        ! NSIGEXC=0 means ground state model.
        ! NSIGEXC=1 means "appropriate applicable" excited states model:
        ! ADAS for LEV_NBIDEP=2 and
        ! Janev-Boley-Post for other values of LEV_NBIDEP=1
        ! NSIGEXC=2 means Janev-Boley-Post applied to any model including ADAS.
        ! Janev/Boley/Post article, Nucl. Fusion Vol.29, No. 12 (1989)
        ! NSIGEXC=3 means ADAS excited states model (adas310 library):
        ! ADAS for LEV_NBIDEP=2 and
  real(fp) :: xdepmod = 1.0_fp
        ! anomalous opacity adjustment
        ! this applies a multiplicative factor to all stopping
        ! cross sections.
  logical :: nlbbcx = .TRUE.
        ! .FALSE. (0) to disable beam-beam charge-exchange and
        ! beam-beam impact-ionization contribution to deposition.
        ! This concerns beam neutrals being stopped by collisions
        ! with slowing down beam ions.
  integer :: nbbcx_avg = 0
        ! use non average over the flux surface relative energy (ebb_rel_avg)
        ! and target temperature (rms_ebb_rel) for beam-beam reaction to calculate
        ! enhancement factor
        ! =1 then use average over the flux surface relative energy (ebb_rel_avgx)
        ! and target temperature (rms_ebb_relx)
        ! this switch is activated if nlbbcx = .TRUE.
  integer :: nbbcx_bb = 0
        ! integer switch. If nbbcx_bb=0 then contribution from beam-beam reactions to the enhancement
        ! factor will be included in case if NLBBCX=.TRUE.
        ! If nbbcx_bb.ne.0  then contribution from beam-beam reactions to the enhancement
        ! factor will be NOT included even if NLBBCX=.TRUE.
  integer :: nbfallgr = 0
        ! if set 1, then program will be stopped, if enhancement
        ! factor computed from ADAS data set became less than 1.0.
        ! This setting will be used for a debugging purposes.
        ! If set 0,then enhancement factor computed from ADAS data set
        ! will be forced to be 1.0 if it became less than 1.0. It means that
        ! ground state data will be used for this region.
        ! This switch will be used if nlbbcx=.true.
  logical :: nlbeamcx = .TRUE.
        ! set .FALSE. (0) to disable beam charge exchange losses
        ! (if .FALSE. beamcx subroutine is never called).  For
        ! technical reasons nlbeamcx=.FALSE. implies nlbbcx=.FALSE.
        ! but cx stopping of neutral beams on thermal species still
        ! occurs, and a beam halo neutral source is still computed.
        ! Caution: nb*sigma*v ionization rate profile for background
        ! neutrals will be ZERO if nlbeamcx=.FALSE. -- this is for
        ! debugging only!
  logical :: nlhvion = .TRUE.
        ! this switch only used for heavy neutral beam deposition
        ! (i.e. experiments using neutral beams to inject Argon or
        ! Neon-- very rare).  If set .TRUE. (1), ion collisions are used
        ! in calculating the stripping of the deposited heavy ion; if
        ! .FALSE. (0) only electron collisions are used.  For ordinary
        ! (Hydrogen or Helium) beams this switch has no effect.
 
  ! block: collid
    ! collision operator controls
 
  LOGICAL :: init_collid = .FALSE.
 
  real(fp) :: xcfanbi = 1.0_fp
        ! anomalous multiplier: electron drag and energy
        ! diffusion on beam ions
  real(fp) :: xdfanbi = 1.0_fp
        ! anomalous multiplier: ion drag and energy
        ! diffusion on beam ions
  real(fp) :: xefanbi = 1.0_fp
        ! anomalous multiplier: pitch angle scattering on
        ! beam ions
  real(fp) :: xcfafus = 1.0_fp
        ! anomalous multiplier: electron drag and energy
        ! diffusion on fusion product ions
  real(fp) :: xdfafus = 1.0_fp
        ! anomalous multiplier: ion drag and energy
        ! diffusion on fusion product ions
  real(fp) :: xefafus = 1.0_fp
        ! anomalous multiplier: pitch angle scattering on
        ! fusion product ions
  logical :: nlbcde = .true.
        ! .FALSE. (0) to turn off energy diffusion
  logical :: nlbcoh = .true.
        ! .FALSE. (0) to suppress acceleration of beam ions by
        ! toroidal electric field
  logical :: nlbcpa = .true.
        ! .FALSE. (0) to turn off pitch angle scattering
  logical :: nlorbo = .false.
        ! .TRUE. (1) for "orbit only" mode: suppress all collisions
        ! used only in specialized debugging situations.
 
  ! block: flr
    ! finite larmor radius corrections
 
  LOGICAL :: init_flr = .FALSE.
 
  logical :: nlbflr = .TRUE.
        ! .TRUE. (1) for standard FLR treatment, based on a mod(B)
        ! at guiding center and notion of a FLR circle about
        ! the guiding center, upon which a random gyro-phase
        ! angle is taken to define the current particle position.
  logical :: nlfbmflr = .TRUE.
        ! .TRUE. (1) to use the FLR localized fast ion
        ! distribution for beam-beam neutron calculations.
        ! Note that recent code changes (by Marina Gorelenkova)
        ! ensure that both an FLR and a GC localized fast ion
        ! distribution are computed and output on every run.
  logical :: nlbgflr = .FALSE.
        ! .TRUE. (1) for enhanced FLR treatment
        ! this deforms the FLR trajectory by looking at the
        ! variation of mod(B) in the guiding center vicinity,
        ! and is appropriate when the larmor radius is >~ L(B)
        ! e.g. in low field systems such as STs.  **but** this
        ! option is extremely expensive computationally!
  real(fp) :: gflr_min = 0.0_fp
        ! for nlbgflr=.TRUE., if the guiding-center estimate
        ! of larmor radius, in cm, is less than this, revert
        ! to the standard guiding-center based treatment (to
        ! save computational effort).
  real(fp) :: gflr_rl = 1.0e+6_fp
        ! nlbgflr=.TRUE. expert control (do not touch)
  real(fp) :: gflr_ll = 30.0_fp
        ! nlbgflr=.TRUE. expert control (do not touch)
  real(fp), dimension(:), allocatable :: gflr_xv
        ! nlbgflr=.TRUE. expert controls (do not touch)
  real(fp), dimension(:), allocatable :: gflr_op
        ! nlbgflr=.TRUE. expert controls (do not touch)
 
  ! block: saw
    ! sawtooth model controls
    ! note: on every fast ion step, sawtooth mixing may be computed, but
    ! not applied.  If in fact a sawtooth occurs, a call to "sawnbi" signals
    ! that the sawtooth transformation of the Monte Carlo population should
    ! be applied.
    !
    ! The "sawnbi" call cannot change the fast ion density and energy
    ! density profiles in the driver code; instead, these profiles are
    ! computed as a part of each regular fast ion timestep.  It is up
    ! to the driver code to save these profiles and use as needed.
    !
    ! When NUBEAM is run in client-server mode, the "sawnbi" call needs
    ! to happen on the server side and is requested via a namelist
    ! setting.
    !
    ! "sawnbi" is a legacy call.  Setting sawflag=.TRUE. has the same
    ! effect as a sawnbi call (cf sawnbi.for & sawflag2 internal flag...)
 
  LOGICAL :: init_saw = .FALSE.
 
  logical :: sawflag = .FALSE.
        ! set .TRUE. (1) if sawtooth occurred since prior step
        ! (this means: apply previously calculated sawtooth
        ! mixed new start coordinates when restarting partially
        ! slowed down orbiting ions).
  logical :: nlsawb = .TRUE.
        ! if .TRUE. (1) calculate sawtooth mixing of beam ions.
  logical :: nlsawf = .TRUE.
        ! if .TRUE. (1) calculate sawtooth mixing of fusion products.
  real(fp) :: xswfrac_allfast = 1.0_fp
        ! fraction multiplier -- fast ions participating
        ! in sawtooth mixing
  real(fp) :: xswfrac_beam = 1.0_fp
        ! fraction multiplier -- beam ions participating
        ! in sawtooth mixing: net xswfrac_allfast*xswfrac_beam
  real(fp) :: xswfrac_fusn = 1.0_fp
        ! fraction multiplier -- fusn prod ions participating
        ! in sawtooth mixing: net xswfrac_allfast*xswfrac_fusn
  integer :: nmix_kdsaw = 1
        ! KDSAW (Kadomtsev/Porcelli) mixing model option
        ! =1: standard Kadomtsev
        ! =2: Kadomtsev, ions fully mixed => flat density
        ! =3: Porcelli, two mixing regions (island annulus +
        !     interior axial region separately mixed).
        ! =4: Porcelli, ions fully mixed => flat density;
        ! -------> see KDSAW NTCC MODULE documentation
  real(fp) :: fporcelli = 1.0_fp
        ! (if nmix_kdsaw = 3 or 4) interior fractional
        ! width of q=1 island towards the magnetic axis;
        ! fporcelli=1 means island reaches all the way to
        ! axis => Kadomtsev-like mixing.
        !   0.0_fp <= fporcelli <= 1.0
        ! smaller Porcelli fractions => smaller mixing
        ! region, but the region from the axis to the q=1
        ! surface is always mixed.
        ! -------> see KDSAW NTCC MODULE documentation
 
  ! block: adif
    ! anomalous diffusion
    ! if this is used, an anomalous diffusivity profile must be
    ! provided via xplasma (see section on input data profiles).
    ! if no profile is provided, the anomalous diffusivity will
    ! be zero and this feature will be inactive.
    !  MOD DMC 2008: (a) anomalous radial velocity profile
 
  LOGICAL :: init_adif = .FALSE.
 
  integer :: nkaep =  0
        ! number of terms (NK) in expansion of MHD transport operator (Mario Podesta)
  real(fp), dimension(:), allocatable :: AEP_ASCALE
        ! scaling factor for mode amplitudes in  MHD transport operator (Mario Podesta)
  real(fp) :: AEP_USCALE =  1.0_fp
        ! scaling factor for electrical potential for MHD transport operator (Mario Podesta)
  real(fp), dimension(:), allocatable :: TAEP_FAC
        ! scaling factor to redefine time step for MHD transport operator (Mario Podesta)
  integer :: nmdifb = 0
        ! type of fast ion diffusivity
        ! =4: 2D fast ion diffusivity from ufiles
        ! =6: fast ion diffusivity from dep code
  integer :: nkdifb = 3
        ! =3: apply to all fast ions;
        ! =1: apply to beam ions only;
        ! =2: apply to fusion product ions only.
        ! NOTE: this variable always set to 3 when Plasma_State is
        ! used; in this case separate diffusivity profiles are
        ! provided for the beam ions and the fusion product ions.
  integer :: ndifbe = 0
        ! number of bins in profile of energy dependence for
        ! anomalous diffusion
  real(fp), dimension(:), allocatable :: fdifbe
        ! anomalous diffusivity multiplier, vs. energy
        ! all with reference to an anomalous diffusivity input
        ! profile Db(x) given as one of the code's profile inputs.
  real(fp), dimension(:), allocatable :: edifbe
        ! energies (eV) for anomalous diffusivity multiplier
        ! energy dependence:
        !   fdifbe(1) applies for energy range 0 to edifbe(1)
        !   fdifbe(2) applies for energy range edifbe(1) to edifbe(2)
        !     ...
        !   fdifbe(ndifbe) applies, edifbe(ndifbe-1) to edife(ndifbe).
        ! if not in any range, a multiplier of 1.0 applies.
  real(fp) :: x_qlim0 =  0.0_fp
        ! QLIM0 high diffusion region -- apply enhanced diffusivity
        ! to *all* fast species inside this radius (regardless of
        ! how nkdifb is set).  (difb profile still determines the
        ! actual diffusivity; ndifb, fdifbe, edifbe can still be
        ! used; velb used if set).
        ! (NOTE: always set to 0.0_fp if Plasma_State is used).
  integer :: ndifbe_2d =  1
        ! number of energy grid points for diffusivity vs. (x,E)
        ! 2d profiles (ndifbe_2d.gt.1 indicates presence of data)
        ! this should match the energy grid of the FD0 diffusivities
        ! nbi_getprofiles will check this.
  real(fp) :: dif2d_conv =  0.0_fp
        ! conversion factor to cm**2/sec
        ! (non-zero value indicates presence of data)
        ! (Plasma State has m**/sec, interface sets dif2d_conv = 1.0e+4_fp)
  real(fp) :: dif2d_dpa =  0.0_fp
        ! if non-zero, this is the target change per orbit
        ! bounce in pitch angle which affects goosing and range of
        ! search of diffusivities in pitch angle space to reduce
        ! goosing near high diffusivity regions; it matches
        ! delta(vpll/v) approximately if << 1.0.
        ! If ZERO it will be set to 1.0/goocon
  integer :: ngradd_opt =  0
        ! grad(Db) pinch cancellation option
        ! 0, default, means automatic: active for 1d diffusivity
        ! profile inputs, inactive for 2d inputs.  Set =1 to force
        ! grad(Db) correction; set to -1 to force its absence.
 
  ! block: ripple
    ! TF field ripple effects.  The available models are crude.
    ! A TFTR project to build a better model ran out of time and
    ! money.
 
  LOGICAL :: init_ripple = .FALSE.
 
  integer :: nrip = 0
        ! =0: no ripple loss;  =1: simple "loss time" model
        ! =2: adjustable model based on Goldston-White-Boozer and
        !     ripple field strength. CAUTION-- calibration against
        !     a more exact model is needed.  PPPL expert help should
        !     be sought before using this model
  real(fp) :: taurip = 0.0_fp
        ! must be set if nrip=1: lint* readMyid0,oss time (sec) for trapped ions
  integer :: id_ripple = 0
        ! xplasma ID for profile f(R,Z) of
        !   log((Bmax-Bmin)/(Bmax+Bmin))
        ! if not set, code looks for an xplasma function named
        ! "LOG_TF_RIPPLE".  Must be found if nrip=2.
  integer :: id_ripple2 = 0
        ! xplasma ID for profile f(R,Z) of
        !   log((Bmax-Bmin)/(Bmax+Bmin)) (2nd component ripple field)
        ! if not set, code looks for an xplasma function named
        ! "LOG_TF_RIPPLE2".  Must be found if nrip=2, AND ncoils_2 > 0.
  integer :: ncoils = 0
        ! number of TF coils.  Must be set if nrip=2.
  integer :: ncoils_2 = 0
        ! number of TF coils (2ndary ripple field).
        ! Must be set if nrip=2 AND a 2 component ripple field is to
        ! be used.
  real(fp) :: asrd = 0.0_fp
        ! nrip=2 adjustment factor for fusion product ions
  real(fp) :: bsrd = 0.0_fp
        ! nrip=2 adjustment factor for beam ions.
 
  ! block: outcon
    ! code output control options
 
  LOGICAL :: init_outcon = .FALSE.
 
  integer :: nerngfi = 1
        ! number of energy bins for trapping fraction output profiles
        ! generally only coarse binning is useful, e.g. nerngfi=5.
        ! the maximum allowed value is the (user settable) array
        ! dimension "miefi".
  real(fp), dimension(:,:), allocatable :: erngfi
        ! eV
        ! coarse energy bins for trapping fraction output profiles
        ! if nerngfi=1, there is one bin: [erngfi(1),infinity]
        ! if nergnfi>1, the bins are: [erngfi(1),erngfi(2)]
        !                             [erngfi(2),erngfi(3)]
        !                               ...
        !                             [erngfi(nerngfi),infinity]
  integer :: nsdbgb = 1
        ! index to fast ion specie selected for debug output
        ! (artifact-- not currently used).
 
  ! block: fishbone
    ! fast ion fishbone loss model (ad hoc)
 
  LOGICAL :: init_fishbone = .FALSE.
 
  logical :: nlfbon = .FALSE.
        ! .TRUE. to turn on fishbone loss model
        ! "fishbones" are an MHD phenomenon though to be capable
        ! of ejecting fast ions.
  real(fp) :: fbemin = 0.0_fp
        ! minimum energy (eV) of affected ions
  real(fp) :: fbemax = 0.0_fp
        ! maximum energy (eV) of affected ions
  real(fp) :: fvpvmn =  0.0_fp
        ! minimum midplane |vpll/v| of affected ions
  real(fp) :: fvpvmx =  0.0_fp
        ! maximum midplane |vpll/v| of affected ions
  real(fp) :: fbltim =  1.0e-7_fp
        ! average loss time (seconds) of affected ions
  real(fp) :: fshper =  5.0e-3_fp
        ! fishbone period (cycle time, seconds)
  real(fp) :: fshwid =  1.0e-3_fp
        ! fishbone duration (seconds)
  real(fp) :: tfshon =  -1000.0_fp
        ! fishbone onset time (seconds) default always on
  real(fp) :: tfshof =  +1.0e+34_fp
        ! fishbone stop time (seconds) default always on
  real(fp) :: xfishmin =  0.0_fp
        ! min (x) where fishbones have an effect
  real(fp) :: xfishmax =  1.0_fp
        ! max (x) where fishbones have an effect
        ! (x) the usual sqrt(tor.flux) normalized coordinate; average position
        ! data from last half orbit is used.
  integer :: nfbon_vpvopt = 0
        ! =0: control vpll/v of affected ions with {fvpvmn,fvpvmx}
        ! =1: only operate on trapped ions with a trapping depth
        !     exceeding the threshold fbtrap_depth (see descr.).
        !     Note (fbemin,fbemax) constraints still apply.
  real(fp) :: fbtrap_depth = 0.0_fp
        ! trapping depth threshold (nfbon_vpvopt = 1 only)
        ! an orbit is considered trapped if either:
        !    (a) max(|B|) over orbit half-bounce > max(|B1|,|B2|)
        !        where B1 and B2 are taken at the midplane crossings (all
        !        |B|s are guiding center values), or,
        !    (b) the orbit never reaches the midplane after a multiple of
        !        tau_bounce has been exceeded (i.e. ST equilibria off-midplane
        !        trapped orbit class).
        ! IF an orbit is deemed trapped, then, its trapping depth is defined as
        !        d = (Bmax(psi_avg)-max(|B|)_orb)/(Bmax(psi_avg)-Bmin(psi_avg))
        ! and it is considered eligible for fishbones if d >= fbtrap_depth.
        ! Here, psi_avg is the "average" flux surface seen by the orbit, defined
        ! as 1/2(psi_min+psi_max), where the min and max values are taken over
        ! the last half bounce;
        !
        ! note that 0 < d < 1 is expected, and is enforced with min & max
        ! functions, as the limits could be exceeded slightly due to the
        ! psi_avg calculation procedure.
  integer :: nfbon_species = 3
        ! species control
        !   =1 beam ions only
        !   =2 fusion product ions only
        !   =3 (default) all fast ions
 
  ! block: box
    ! beam-in-box neutral density calculation controls
 
  LOGICAL :: init_box = .FALSE.
 
  integer :: nbbox = 0
        ! number of "beam-in-box" neutral density profile
        ! calculations
  integer, dimension(:), allocatable :: nbsbox
        ! beam source numbers for neutral density calculation.
  integer, dimension(:), allocatable :: nsample_box
        ! # of lanched MC particles in each box
  real(fp), dimension(:), allocatable :: xninja_box
        ! NO. OF PTCLS COMING IN Ith BOX  N/SEC
  real(fp), dimension(:), allocatable :: xninjs_box
        ! n/sec NO. of beam ptcls.injected, by species for BOXes subset only
  real(fp) :: xninj_box = 0
        ! NO. OF PTCLS INJECTED N/SEC for BOXes subset only
  integer, dimension(:), allocatable :: nbebox
        ! beam energy fractions, before which to compute a 3d
        ! neutral density profile on a Cartesian grid alligned
        ! on the beam centerline.
        !   nbsbox(i)=j  nbebox(i)=k
        ! means the i'th neutral density calculation is done in
        ! front of beam source j, using energy fraction k (for H
        ! isotope beams k=1 means full energy, k=2 means half-
        ! energy, k=3 means 1/3 energy).
  integer :: nxbox = 1
        ! number of bins in horizontal direction (x) in front of beam
  integer :: nybox = 1
        ! number of bins in vertical direction (y) in front of beam
  integer :: nlbox = 100
        ! number of bins in the direction (l) of beam injection
  real(fp) :: xboxhw = 5.0_fp
        ! (cm) half-width of entire Cartesian box in x direction
        ! individual bins' x width = 2*xboxhw/nxbox
  real(fp) :: yboxhw = 5.0_fp
        ! (cm) half-height of entire Cartesian box in y direction
        ! individual bins' y width = 2*yboxhw/nybox
  real(fp) :: xlbox1 = 100.0_fp
        ! (cm) distance from beam source at which box starts
  real(fp) :: xlbox2 = 500.0_fp
        ! (cm) distance from beam source at which box ends
        ! individual bins' l length = (xlbox2-xlbox1)/nlbox
  integer :: ndepbox = 500
        ! Monte Carlo statistics control:  minimum number of neutral
        ! tracks to use per calculation.
 
  ! block: tube
    ! beam-in-tube neutral density variables
 
  LOGICAL :: init_tube = .FALSE.
 
  integer :: nbtube = 0
        ! number of beam-in-tube neutral density calculations
  integer, dimension(:), allocatable :: nbstube
        ! tube beam source numbers
  integer, dimension(:), allocatable :: nbetube
        ! beam energy fractions to be used when calculating neutral
        ! density
        !   nbstube(i)=j  nbetube(i)=k
        ! means the i'th neutral density calculation is done in
        ! front of beam source j, using energy fraction k (for H
        ! isotope beams k=1 means full energy, k=2 means half-
        ! energy, k=3 means 1/3 energy).
  real(fp), dimension(:), allocatable :: Rtube
        ! radial coordinate of start of tube (cm)
  real(fp), dimension(:), allocatable :: Ytube
        ! vertical (z-axis) coordinate of start of tube (cm)
  real(fp), dimension(:), allocatable :: XZETAtube
        ! angular coordinate of start of tube (degrees)
  real(fp), dimension(:), allocatable :: PHItube
        ! tube's sightline azimuthal angle from start of tube measured from
        ! the centerline with the positive direction in direction of the
        ! plasma current (degrees)
  real(fp), dimension(:), allocatable :: THETAtube
        ! tube's sightline elevation angle from start of tube measured from
        ! the midplane (degrees)
  integer, dimension(:), allocatable :: LMIDtube
        ! 0->tube data will be collected between XL1tube<s<XL2tube
        ! 1->tube data will be collected from XL1tube to XL2tube relative
        !    to the point on the tube which is closest to the beam's sightline
  real(fp), dimension(:), allocatable :: XL1tube
        ! tube data will be collected between XL1tube and XL2tube along sightline.
        ! this is the distance at which data will begin to be collected (cm)
  real(fp), dimension(:), allocatable :: XL2tube
        ! tube data will be collected between XL1tube and XL2tube along sightline.
        ! this is the distance at which data will stop being collected (cm)
  integer, dimension(:), allocatable :: NSEGtube
        ! number of data segments between XL1tube and XL2tube in which to compute
        ! the neutral density data.  This number must be less then mituben.
  real(fp) :: RHOtube =  1._fp
        ! radius of the tubes (cm) -- note that the number of trajectories used to
        ! compute the neutral particle density (tuben0) is inversely proportional
        ! to the tube radius.
  integer :: NDEPtube = 500
        ! estimated number of beam intersections per tube assuming no beam loss and a tube
        ! perfectly centered within the beam aperture.  The number of trajectories which
        ! will be attempted will be  NDEPtube*sqrt(area_aperture)/2*RHOtube
 
  ! block: misc
    ! miscellaneous
 
  LOGICAL :: init_misc = .FALSE.
 
  real(fp) :: plhgt = 0.0_fp
        ! (updown symmetric geometries only) height of plasma
        ! midplane above/below vacuum vessel midplane (cm);
        ! neutral beam component (e.g. ion source, aperture)
        ! elevations are given with respect to the vacuum vessel
        ! midplane.  (In updown asymmetric cases the vertical
        ! offset is presumed to be given in the MHD equilibrium).
        !   **time dependent quantity**
  logical :: nlcprb = .TRUE.
        ! set .FALSE. (0) to suppress compression operator and other
        ! effects of time varying MHD equilibrium on beam ions
        ! (recommendation: leave default, .TRUE., alone).
  logical :: nlpsirz = .FALSE.
        ! set .TRUE. (1) if an accurate Psi(R,Z) (free boundary) MHD
        ! equilibrium, extending beyond the plasma boundary, is
        ! available.
  real(fp) :: edbfac = 10.0_fp
        ! debug tripwire -- if an ion reaches edbfac x its initial
        ! deposition energy, kill the run. (dimensionless).
  integer :: nmcurb = 1
        ! set to 0 to suppress fast ion driven current.
        ! nmcurb=1 gives an aspect-ratio approximate neoclassical
        ! shielding calculation for the current; nmcurb=2 gives
        ! the simplified shielding correction (1-1/Zeff).
        ! nmcurb=3  gives neoclassical shielding calculation for the current
        ! by Lin-Liu formula Phys.Plasmas 4(11),p.4179, 1997. Other
        ! non-zero values of nmcurb give the same result as nmcurb=1
        ! (for now), but, a fully generalized geometry neoclassical
        ! shielding calculation will be implemented at some point
        ! in the future, and nmcurb=3 may be used to select this,
        ! when it exists.
        ! NMCURB  = 4, shielding calculation by Honda M., Kikuchi M., Azumi M.
        ! "Collisionality dependence on shielding factor of beam driven current"
        ! Nuclear Fusion, submitted by October 2011. Martix inversion method are used
        ! to include collisionality dependence on a shielding factor at arbitrary aspect ratio.
  logical :: nlfdep = .FALSE.
        ! set .TRUE. (1) to compute midplane index oriented deposition
        ! distribution function in the lab frame using the "fbm"
        ! grids for pitch and energy, with particles mapped to
        ! their outer midplane crossing (a few particles which
        ! never cross the midplane at R > Raxis will be missed!).
        ! setting this switch causes a specialized, rarely used
        ! output data structure to be created
        ! nlfdep ***DISABLED*** DMC May 2 2007
        ! deposition data is now output on original chordal tracks
 
  ! block: profiles
    ! xplasma ids of profile inputs to the Monte Carlo fast
    ! ion calculation: generally, plasma parameter profiles.
    ! the code also expects an xplasma equilibrium description
    ! with vacuum region and limiters defined.
    ! the flux coordinate "x" is normalized sqrt-toroidal-flux
    ! (0 at axis, 1 at plasma bdy).
    ! note: all xplasma quantities are MKS; temperatures are KeV
    ! (the beam code converts to cgs-eV when it interpolates the
    ! profiles to the Monte Carlo code grid).
    ! note: loop voltage profile is in volts.
    !       radial electrostatic potential is in volts.
 
  LOGICAL :: init_profiles = .FALSE.
 
  integer :: id_te = 0
        ! electron temperature vs "x"
  integer :: id_ti = 0
        ! ion temperature vs "x"
  integer :: id_tifusn = 0
        ! ion temperature for fusion calculations vs. "x"
        ! id_ti = id_tifusn is not unusual.  Sometimes id_ti
        ! is affected by impurities, and id_tifusn is supplied
        ! as slightly more accurate for the reactive non-impurity
        ! species.
  integer :: id_ne = 0
        ! electron density /m**3 in xplasma
        ! converted to /cm**3 internally
  integer :: ng = 0
        ! current number of gases
  integer, dimension(:), allocatable :: id_nh
        ! non-impurity thermal ion densities
  integer, dimension(:), allocatable :: id_nimp
        ! impurity thermal ion densities
  integer, dimension(:), allocatable :: id_nrf
        ! RF minority specie densities
  integer, dimension(:), allocatable :: id_eperpmin
        ! RF minority specie <Eperp>, lab frame
  integer, dimension(:), allocatable :: id_epllmin
        ! RF minority specie <Epll>, lab frame
  integer, dimension(:), allocatable :: id_nfast
        ! fast density profiles can be input e.g. if a fokker
        ! planck package (instead of the Monte Carlo package)
        ! is used for slowing down calculations.
  integer, dimension(:), allocatable :: id_eperpfast
        ! fast ion <Eperp>, lab frame can be input
  integer, dimension(:), allocatable :: id_epllfast
        ! fast ion <Epll>, lab frame can be input
  integer :: id_aimp = 0
        ! Aimp vs. x (allowed if there is a single "model impurity")
        ! AMU; non-integer values allowed as this represents an
        ! amalgum of impurities.
  integer :: id_zimp = 0
        ! Zimp vs. x (allowed if there is a single "model impurity")
        ! atomic number; non-integer values allowed as this
        ! represents an amalgum of impurities.
  integer :: id_omega = 0
        ! plasma angular velocity profile "omega" (1/sec)
        ! locally Vphi = R*omega
  integer :: id_epot = 0
        ! radial electrostatic potential profile (volts)
  integer :: id_vpoh = 0
        ! toroidal loop voltage profile
  integer :: id_curt = 0
        ! enclosed toroidal current vs. flux surface x (amps)
        !------------------
        ! one can either specify these profiles and use nkdifb to control application:
  integer :: id_difb = 0
        ! anomalous diffusivity profile (all species) (cm2/sec)
  integer :: id_velb = 0
        ! anomalous radial velocity profile (all species) (cm/sec)
        !------------------
        ! or specify anomalous transport separately for beam ions and fusion ions:
  integer :: id_difb_bi = 0
        ! anomalous diffusivity profile (beam ions) (cm2/sec)
  integer :: id_velb_bi = 0
        ! anomalous radial velocity profile (beam ions) (cm/sec)
  integer :: id_difb_fi = 0
        ! anomalous diffusivity profile (fusion ions) (cm2/sec)
  integer :: id_velb_fi = 0
        ! anomalous radial velocity profile (fusion ions) (cm/sec)
        !------------------
        ! or specify orbit-topology-dependent diffusivities vs. (x,E) separately for
        ! each fast ion species  (dmc: 2010 collab w/ B. Heidbrink & N. Gorelenkov).
  integer, dimension(:), allocatable :: id_fd0
        ! diffusivities: deeply trapped orbits
  integer, dimension(:), allocatable :: id_fdb
        ! diffusivities: barely trapped orbits
  integer, dimension(:), allocatable :: id_fdp
        ! diffusivities: barely passing "co" orbits
  integer, dimension(:), allocatable :: id_fdq
        ! diffusivities: deeply passing "co" orbits
  integer, dimension(:), allocatable :: id_fdr
        ! diffusivities: barely passing "ctr" orbits
  integer, dimension(:), allocatable :: id_fds
        ! diffusivities: deeply passing "ctr" orbits
        !------------------
  integer, dimension(:,:), allocatable :: id_n0
        ! neutral density profiles
        ! id_n0(is,ig) = neutral density of gas ig due to source is
        !   (because of charge exchange, a source of one specie
        !   generally produces non-zero neutral densities for several
        !   species).
        ! source numbering:
        !   is=1 -- wall source, gas ig=1
        !   is=2 -- halo source, gas ig=1
        !   is=3 -- optional 2nd wall source, gas ig=1
        !           if no 2nd wall source leave id_n0(is,1:inum_ng0)=0
        !   is=4 -- wall source, gas ig=2
        !   is=5 -- halo source, gas ig=2
        !   is=6 -- optional 2nd wall source, gas ig=2
        !    ...
        !   is=3*inum_ng0-2 -- wall source, gas ig=inum_ng0
        !   is=3*inum_ng0-1 -- halo source, gas ig=inum_ng0
        !   is=3*inum_ng0   -- optional 2nd wall source,
        !                      gas ig=inum_ng0
        !
        ! note: at present the neutrals profiles must be given
        ! as flux surface averages, but in the future, 2d
        ! profiles f(x,theta) will be allowed.
  integer, dimension(:,:), allocatable :: id_t0
        ! neutral "temperature" (or 2/3 <E>) profiles
        ! same indexing as id_n0
  integer, dimension(:,:), allocatable :: id_omeg0
        ! neutral populations' average angular velocity profiles
        ! same indexing as id_n0
 
  !==============================================
  ! INTERNAL items...
 
#ifdef __MPI
  logical :: debug_mpi = .TRUE.
#else
  logical :: debug_mpi = .FALSE.
#endif
        ! as a default
        ! debug_mpi= nltest_output for
        ! MPI job debugging
        ! debug_mpi   =.FALSE. for serial mode
        ! see nltest_output.
        ! In nbi_state_io.90, if
        ! variable not in fld state file:
        ! default value will be set
  integer :: nbsteps = 0
        ! number of completed NB steps so far...
  integer :: ndepbox_use = 0
        ! ndepbox used (for deposition in beam-in-box statistics)
  integer :: mdepbox = 0
        ! # of deposition in 'beam-in-box' model for each CPU
        ! ndepbox_use = mdepbox * numprocd
  integer :: I_trk0 = 0
        ! neutral track ID for beam in box deposition
  real(fp) :: wmin_halo = zero
        ! min allowed track weight for 3D halo model, depends on total neutral current
        ! and MDEPBOX * NSPLT
  integer :: nsigexc_use = 0
        ! = NSIGEXC, but
        !if  NSIGEXC=1 or =3 but there are no data for enhancement factor,
        !then forced nsigexc_use = 0
  integer :: kmix_kdsaw = 0
        ! index of outermost zone involved in sawtooth mixing.
  real(fp) :: xmix_kdsaw = zero
        ! xi value of outer boundary of sawtooth mixing region.
  logical, dimension(:), allocatable :: flag_cursaw
        ! per species flag for computation of sawtooth
  real(fp), dimension(:,:), allocatable :: dxsm_kdsaw
        ! average-orbit-width-based smoothing
        ! for post-sawtooth profile
  real(fp), dimension(:,:), allocatable :: dxwt_kdsaw
        ! weighting for dxsm_kdsaw sums
  real(fp), dimension(:), allocatable :: wsum_kdspre
        ! pre-sawtooth mixing region weight sum
  real(fp), dimension(:), allocatable :: wsum_kdspost
        ! post-sawtooth mixing region weight sum
  real(fp), dimension(:,:), allocatable :: wts_ien
        ! weighting for EIN_IEN and DEI_IEN sums
        ! ebb_rel_co, ebb_rel2_co, rms_ebb_co, ebb_relx_co, rms_ebbx_co
        !
        ! are arrays separated by 2d spatial zone, injected energy
        ! fraction, beam line, beam species for co-injected beam populations
        !
        ! ebb_rel_cnt, ebb_rel2_cnt, rms_ebb_cnt, ebb_relx_cnt, rms_ebbx_cnt
        !
        ! same as above for counter-injected beam populations
        !
        ! ebb_rel_fn, ebb_rel2_fn, bnb_ien2_fn, rms_ebb_fn, ebb_relx_fn, rms_ebbx_fn
        !
        ! same as above for fusion product populations
        !
        ! bnb_ien2 -  total ptcls density by direction of injection and their species, on 2d grid zone rows
        ! bnb_ienx   is bnb_ien2 fixed on 2d grid zone rows
        ! bnb_fn - fusion product density on 2d grid
        ! bnb_fnx - fusion product density fixed  on 2d grid zone rows
        ! beam-beam interaction
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel_co
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel2_co
        !"ev**2; avg relative energy**2 for beam-beam reactions"
        ! AVG RELATIVE ENERGY**2 FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebb_co
        !"ev; RMS or relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_relx_co
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebbx_co
        !"ev; RMS of relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel_cnt
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel2_cnt
        !"ev**2; avg relative energy**2 for beam-beam reactions"
        ! AVG RELATIVE ENERGY**2 FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebb_cnt
        !"ev; RMS or relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_relx_cnt
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebbx_cnt
        !"ev; RMS or relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:,:), allocatable :: bnb_ien2
        !  beam ptcls density by direction of injection and
        !  fast ion species  on 2d spatial grid
  real(fp), dimension(:,:,:), allocatable :: bnb_ienx
        !  beam ptcls density by direction of injection and
        !  fast ion species fixed on 2d spatial grid
        ! beam ptcls - fusion products
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel_fn
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_rel2_fn
        !"ev**2; avg relative energy**2 for beam-beam reactions"
        ! AVG RELATIVE ENERGY**2 FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebb_fn
        !"ev; RMS or relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION
  real(fp), dimension(:,:,:,:), allocatable :: ebb_relx_fn
        ! "ev; avg relative energy for beam-beam reactions"
        ! AVG RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:,:,:), allocatable :: rms_ebbx_fn
        !"ev; RMS or relative energy for beam-beam reactions"
        ! RMS RELATIVE ENERGY FOR BEAM-BEAM REACTION, fixed on 2d grid zone rows
  real(fp), dimension(:,:), allocatable :: bnb_fn
        ! fusion product density on 2d spatial grid
  real(fp), dimension(:,:), allocatable :: bnb_fnx
        !  fusion product density fixed o on 2d spatial grid
  real(fp) :: tfracp = zero
        ! (TBM2-TBM1)[current step] /
        !   (TBM2-TBM1)[previous step]; use this for
        !   a priori goosing adjustment if necessary
  real(fp), dimension(:), allocatable :: orb_code_error
        ! orbit code error count sum
  real(fp), dimension(:), allocatable :: dep_code_error
        ! dep code error count sum
  real(fp), dimension(:), allocatable :: xnorbs
        ! orbit counter
  logical :: flag_dt_reset = .FALSE.
        ! flag to reset orbit timestep adjust
  real(fp), dimension(:), allocatable :: dt_adj0_sums
        ! orbit timestep adjust sums
  real(fp), dimension(:), allocatable :: xnum_dt_adj0
        ! timestep adjust sum weightings
  real(fp), dimension(:), allocatable :: xmcloss_o
        ! prompt loss MC ptcl count (orball)
  real(fp), dimension(:), allocatable :: xmcloss_d
        ! prompt loss MC ptcl count (depall)
  real(fp), dimension(:), allocatable :: weinjav
        ! wght for einjav calculation (depall)
  real(fp), dimension(:), allocatable :: ennew
        ! new energy injected this step
  real(fp), dimension(:), allocatable :: ewghtdep
        ! summed E*weight deposited by DEP calls
  real(fp), dimension(:), allocatable :: ewghtdep0
        ! as ewghtdep -- but for subset w/saved MC ions
  real(fp), dimension(:), allocatable :: wghtdep
        ! summed weight deposited by DEP calls
  real(fp), dimension(:), allocatable :: zidep0
        ! number of random generator calls to collect
        ! given number of tracks (ndep0). Also used
        ! for statictics improvment.
  real(fp), dimension(:), allocatable :: zitrk0
        ! number of random generator calls to collect
        ! given number of tracks (ndep0) only.
  real(fp), dimension(:), allocatable :: xhwdep0
        ! estimated full E. orbit half-width near axis
  real(fp), dimension(:,:), allocatable :: sfrac_survive
        ! deposition fraction sum (vs E-frac, beam#)
  real(fp), dimension(:,:), allocatable :: wfrac_survive
        ! wght for dep fraction sum (vs E-frac, beam#)
  real(fp), dimension(:,:), allocatable :: frac_survive
        ! deposition fraction sum from prev. timestep
  real(fp), dimension(:,:), allocatable :: frac_sur_inv
        ! 1/frac_survive
  real(fp), dimension(:,:), allocatable :: bldistrib
        ! beamline selection probability (from powers
        ! combined with frac_survive(...))
        ! NOTE bldistrib not yet in use DMC Dec 2005
  real(fp), dimension(:,:), allocatable :: wtqbxfs
        ! adjusted wght for tqbxfs sum
  integer :: ncx0_use = 0
        ! used for CX tracks, also see ncx0_orb in section mc_satastics
  real(fp) :: wght_cx = zero
        ! weight use for improving CX statistics
  real(fp) :: xncx0 = zero
        !"# OF CX events during orbiting"
        ! the number of CX tracks that we get from orbiting, from cxnsum
  integer :: norbtot = 0
        ! total number of orbits followed: cf orball.for
  logical :: nlvphi = .FALSE.
        ! .TRUE. if there is non-zero toroidal rotation data
  integer :: ndepmod = 0
        ! flag xdepmod.ne.1
  real(fp) :: abeam = zero
        ! *ATOMIC WEIGHT OF NEUTRAL BEAM NUCLEI
  real(fp), dimension(:), allocatable :: abeams
  real(fp), dimension(:), allocatable :: aimpj
        ! Average impurity atomic weight by zones
  real(fp), dimension(:,:,:), allocatable :: allame
        ! ELECTRON LOG LAMBDA TABLE
  real(fp), dimension(:,:,:,:), allocatable :: allamh
        ! PLASMA ION LOG LAMBDA TABLE CF FOKKER SUBROUTINE
  real(fp), dimension(:,:,:), allocatable :: allamx
        ! IMPURITY LOG LAMBDA TABLE
  real(fp), dimension(:,:), allocatable :: aratii
        ! inv. aspect ratio (midplane definition, extended grid)
  real(fp), dimension(:,:), allocatable :: bbnvol
        ! beam beam fusion
  real(fp), dimension(:,:), allocatable :: bdens2
        ! beam density calculated from fast ion distribution function at GC
  real(fp), dimension(:,:), allocatable :: bdenssp
        ! past timestep beam density
  real(fp), dimension(:,:), allocatable :: bdns2p
        ! BDENS2, PREV. STEP, FOR BB DEPO/RECAPTURE
  real(fp), dimension(:,:), allocatable :: be0cx
        ! AVG ENERGY OF CH-EX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: be0cx2
        ! 2D AVG ENERGY OF CX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: be20cx2
        ! 2D AVG ENERGY^2 OF CX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: be0x2p
        ! BE0CX2, PREVIOUS TIMESTEP
  real(fp), dimension(:), allocatable :: bi1g1x
        ! d(G)/dxi for orbg
  real(fp), dimension(:), allocatable :: bi1ovb
        ! 1/MOD(B) (** CGS **) (BSUB SPLINE FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bi1phi1x
        ! d(PHI)/dxi for orbg
  real(fp), dimension(:), allocatable :: bi1psi1x
        ! d(PSI)/dxi for orbg
  real(fp), dimension(:), allocatable :: bi2g2x
        ! d2(G)/d2xi for orbg
  real(fp), dimension(:), allocatable :: bi2phi2x
        ! d2(PHI)/d2xi for orbg
  real(fp), dimension(:), allocatable :: bi2psi2x
        ! d2(PSI)/d2xi for orbg
  real(fp), dimension(:), allocatable :: bi2x2t
        ! d2x/dt2 (x-y twiddle space near origin)
  real(fp), dimension(:), allocatable :: bi2y2t
        ! d2y/dt2 (x-y twiddle space near origin)
  real(fp), dimension(:), allocatable :: bia
        ! BEAM PARITCLE ATOMIC WEIGHT UNITS
  real(fp), dimension(:), allocatable :: bib2ts
        ! D2/D(THETA-STAR)2 (1/B) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bib2xi
        ! D2/D(XI)2 (1/B) (1/B) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bibdts
        ! D/D(THETA-STAR) (1/B) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bibdxi
        ! D/D(XI) (1/B) (1/B) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bibtsx
        ! D/D(XI) (D/D(TH.-ST.) (1/B)) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bicos
        ! COS(BITHST)
  real(fp), dimension(:), allocatable :: bidb
        ! change in B during orbit step
  real(fp), dimension(:), allocatable :: bidphi
        ! change in radial electric potential during orbit step
  real(fp), dimension(:), allocatable :: bidpph
        ! change in ion canonical momentum during orbit step
  real(fp), dimension(:), allocatable :: bidqdx
        ! DQ/DX OF MC ION
  real(fp), dimension(:), allocatable :: bidt
        ! INPUT=PREVIOUS DT, OUTPUT=CURRENT DT
  real(fp), dimension(:), allocatable :: bidtc
        ! BEAM ION COMPRESSION TIMESTEP
  real(fp), dimension(:), allocatable :: bidte
        ! change in ion total energy during orbit step
  real(fp) :: bidtnv = zero
        ! NORM VEL.*TIME STEP FOR 1 KEV PASSING ION
  real(fp), dimension(:), allocatable :: bidts
        ! change in THETA_STAR during XI-THETA_STAR orbit step
  real(fp), dimension(:), allocatable :: bidvpa
        ! change in V_parallel during orbit step (from V_par eq)
  real(fp), dimension(:), allocatable :: bidx
        ! change in X during X-Y orbit step
  real(fp), dimension(:), allocatable :: bidxdt
        ! dx/dt (x-y twiddle space near origin)
  real(fp), dimension(:), allocatable :: bidxi
        ! change in XI during XI-THETA_STAR orbit step
  real(fp), dimension(:), allocatable :: bidy
        ! change in Y during X-Y orbit step
  real(fp), dimension(:), allocatable :: bidydt
        ! dy/dt (x-y twiddle space near origin)
  real(fp), dimension(:), allocatable :: big
        ! G = R*Bt for orbg
  real(fp), dimension(:), allocatable :: bijac
        ! JACOBIAN, CM3, (BSUB SPLINE FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bijdts
        ! D/D(THETA-STAR) (JAC) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bijdxi
        ! D/D(XI) (JAC) (JAC) (FOR ORBG, ORBGXY)
  real(fp), dimension(:), allocatable :: bike
        ! KINETIC ENERGY OF EACH ION (ERGS, CONSERVED)**ORB output**
  real(fp), dimension(:), allocatable :: biksid
        ! V_PAR/V OF EACH ION
  real(fp), dimension(:), allocatable :: bim
        ! MASS OF EACH ION (GRAMS) **output from ORBit routines**
  real(fp), dimension(:), allocatable :: bimodb
        ! MOD(B) OF EACH ION AT ORBIT TIMESTEP, TESLA
  real(fp), dimension(:), allocatable :: bimu
        ! MAGNETIC MOMENT OF EACH ION
  real(fp), dimension(:), allocatable :: biphi
        ! Electric potential for orbg
  real(fp), dimension(:), allocatable :: bipphi
        ! BEAM ION CANONICAL ANGULAR MOMENTUM
  real(fp), dimension(:), allocatable :: bipsi
        ! Poloidal Flux with 2*pi normalization for orbg
  real(fp), dimension(:), allocatable :: biq
        ! q at guiding center
  real(fp), dimension(:), allocatable :: bisin
        ! SIN(BITHST)
  real(fp), dimension(:), allocatable :: bite
        ! total ion energy after orbit step
  real(fp), dimension(:), allocatable :: bith2t
        ! d2(THETA)/dt**2 for orbit step
  real(fp), dimension(:), allocatable :: bithdt
        ! d(THETA)/dt for orbit step
  real(fp), dimension(:), allocatable :: bithst
        ! THETA-STAR OF EACH ION
  real(fp), dimension(:), allocatable :: bitime
        ! TIME FOR EACH ION
  real(fp), dimension(:), allocatable :: biv
        ! ION VELOCITY (CM/SEC)
  real(fp), dimension(:), allocatable :: bivp2t
        ! d2(V_par)/dt**2 for orbit step
  real(fp), dimension(:), allocatable :: bivpar
        ! PARALLEL VELOCITY OF EACH ION (CM/SEC)
  real(fp), dimension(:), allocatable :: bivpdt
        ! d(V_par)/dt for orbit step
  real(fp), dimension(:), allocatable :: bivprp
        ! PERP VELOCITY OF EACH ION
  real(fp), dimension(:), allocatable :: bixi
        ! XI OF EACH ION = SQRT(PSI_T/PSI_T_LIM)
  real(fp), dimension(:), allocatable :: bixi2t
        ! d2(XI)/dt**2 for orbit step
  real(fp), dimension(:), allocatable :: bixidt
        ! d(XI)/dt for orbit step
  real(fp), dimension(:), allocatable :: biz
        ! BEAM PARITCLE CHARGE STATE
  real(fp), dimension(:), allocatable :: bmnmx
        ! min & max B on each flux surface  (T)
        ! bmnmx(1)=max B on surface lep1
        ! bmnmx(2*nzones+1) = min B on surface lep1
        ! bmnmx(2)=max B on surface ledge=lep1-1
        ! bmnmx(2*nzones) = min B on surface ledge=lep1-1
        !   ...
        ! bmnmx(nzones+1) = min B = max B = B on axis (singular surface)
  real(fp), dimension(:), allocatable :: rmnmx
        ! R at min & max mod(B) pts on each flux surface
  real(fp), dimension(:), allocatable :: bpbmnmx
        ! Bp/B at min & max B on each flux surface
  real(fp), dimension(:), allocatable :: btbmnmx
        ! Bt/B at min & max B on each flux surface
  real(fp), dimension(:), allocatable :: drmnmx
        ! <delta(r)> flux surface spacing at min & max B, each surf.
  real(fp), dimension(:), allocatable :: psmnmx
        ! psi on grid 1:2*nzones+1 like rmnmx, bmnmx -- although psi
        ! is a flux surface constant so this is not strictly necessary, but
        ! it is convenient in the old xja* subroutines.
  real(fp), dimension(:), allocatable :: bmaxm
        ! MAXIMUM OF MOD(B) ON EACH FLUX SURFACE.It is at zone boundaries.
  real(fp), dimension(:), allocatable :: bminm
        ! MINIMUM OF MOD(B) ON EACH FLUX SURFACE
  real(fp), dimension(:), allocatable :: rbminm
        ! R at which each BMINM occurs
  real(fp), dimension(:), allocatable :: rbmaxm
        ! R at which each BMAXM occurs
  real(fp), dimension(:), allocatable :: thbminm
        ! theta at which each BMINM occurs -- caution, branch cuts and
        ! jumps due to multiple minima on different surfaces could occur
  real(fp), dimension(:), allocatable :: thbmaxm
        ! theta at which each BMAXM occurs -- caution, branch cuts and
        ! jumps due to multiple maxima on different surfaces could occur
  real(fp), dimension(:), allocatable :: bmvol
        ! 2d grid zone volumes (cm**3)
  real(fp), dimension(:,:,:,:), allocatable :: bn002
        ! 2D BEAM 1ST GEN NEUTRAL DENSITY
  real(fp), dimension(:,:), allocatable :: bn0cx2
        ! 2D BEAM CX NEUTRAL DENSITY,Nth generation
  real(fp), dimension(:,:), allocatable :: bn0x2p
        ! 2D BEAM CX NEUTRAL DENS.,bn0cx2 on PREV. TIMESTEP
  real(fp), dimension(:), allocatable :: bntots
        ! BEAM ION POPULATION BY ION SPECIES, PREV STEP
  real(fp), dimension(:), allocatable :: boxn0
        ! neutral-density data
  real(fp) :: bphion = zero
        ! ANGULAR MOMENTUM (M*VPHI*R) OF ION (G.C.)
  real(fp), dimension(:,:), allocatable :: btiltc
        ! coarse B-field on 2d grid
  real(fp), dimension(:,:), allocatable :: btntn2
        ! beam target fusion
  real(fp) :: btori = zero
        ! TOROIDAL B FIELD AT ION G.C. ?
  real(fp), dimension(:,:), allocatable :: bv0cx2
        ! 2D AVG V(PHI) OF CX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: bv0x2p
        ! BV0CX2, PREVIOUS TIMESTEP
  real(fp), dimension(:,:), allocatable :: bvrhocx2
        ! 2D AVG V_rho OF CX BEAM NEUTRALS, FLUX SURFACE NORMAL DIRECTION
  real(fp), dimension(:,:), allocatable :: bvrhox2p
        ! BVrhoCX2, FLUX SURFACE NORMAL DIRECTION PREVIOUS TIMESTEP
  real(fp), dimension(:,:), allocatable :: bvthcx2
        ! 2D AVG V_theta OF CX BEAM NEUTRALS, FLUX SUFACE TANGENTIAL DIRECTION IN POLOIDAL PLANE
  real(fp), dimension(:,:), allocatable :: bvthx2p
        ! BVthCX2, PREVIOUS TIMESTEP
  real(fp), dimension(:,:), allocatable :: bvpllcx2
        ! 2D AVG V_||/V OF CX BEAM NEUTRALS, pitch angle
  real(fp), dimension(:,:), allocatable :: bvpllx2p
        ! BVPLLCX2, pitch angle PREVIOUS TIMESTEP
  real(fp), dimension(:,:), allocatable :: balpha_gyro
        ! 2D AVG alpha_gyro OF CX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: balpha_gyrop
        ! Balpha_gyro,  PREVIOUS TIMESTEP
  real(fp) :: bzxr = zero
        ! NORMALIZED B FIELD STRENGTH DATA
  real(fp), dimension(:), allocatable :: cloge
        ! electron coulomb log
  real(fp) :: cmskip = zero
        ! BEAM COMPRESSION CALCULATION STEP SKIP FACTOR
  real(fp) :: cospho = zero
        ! COS OF ORBIT GYRO-PHASE ANGLE
  real(fp), dimension(:), allocatable :: curt
        ! AREA INTEGRAL TOROIDAL CURRENT, AMPS, TO BDY J
  real(fp), dimension(:,:), allocatable :: cxpray
        ! CH-EX CALL FREQUENCY CONTROL ARRAY
  real(fp) :: cxskip = zero
        ! *BEAM: SKIP FACTOR ON BEAM ION CX CALCULATION
  real(fp) :: delt = zero
        ! LENGTH OF (preceding) BEAM TIMESTEP
  real(fp) :: tbm2p = zero
        ! stop time of (preceding) beam timestep
  real(fp), dimension(:,:,:), allocatable :: den0mn
        ! NEUTRAL DENSITY OF HYDROGENIC SPECIES DUE TO SOURCE
  real(fp), dimension(:), allocatable :: difb
        ! (generic) anomalous fast ion diffusion cm2/sec
  real(fp), dimension(:), allocatable :: velb
        ! (generic) anomalous fast ion radial velocity cm/sec
  real(fp), dimension(:), allocatable :: difb_bi
        ! anomalous beam ion diffusion cm2/sec
  real(fp), dimension(:), allocatable :: velb_bi
        ! anomalous beam ion radial velocity cm/sec
  real(fp), dimension(:), allocatable :: difb_fi
        ! anomalous fusion ion diffusion cm2/sec
  real(fp), dimension(:), allocatable :: velb_fi
        ! anomalous fusion ion radial velocity cm/sec
  real(fp), dimension(:,:), allocatable :: wgvdx
        ! weight profile, for computing average diffusivities and velocities
        ! see difbx(...,...); velbx(...,...)
  real(fp), dimension(:), allocatable :: edifbe_2d
        ! energy grid, all diffusivities w/(x,E) variation
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fd0
        ! FD0: diffusivity, deeply trapped ions
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fdb
        ! FDB: diffusivity, barely trapped ions
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fdp
        ! FDP: diffusivity, barely passing "co" ions
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fdq
        ! FDQ: diffusivity, deeply passing "co" ions
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fdr
        ! FDR: diffusivity, barely passing "ctr" ions
  real(fp), dimension(:,:,:), allocatable :: difbe_2d_fds
        ! FDS: diffusivity, deeply passing "ctr" ions
  real(fp) :: w_fd0 = zero
        ! FD0 weight, current half bounce
  real(fp) :: w_fdb = zero
        ! FDB weight, current half bounce
  real(fp) :: w_fdp = zero
        ! FDP weight, current half bounce
  real(fp) :: w_fdq = zero
        ! FDQ weight, current half bounce
  real(fp) :: w_fdr = zero
        ! FDR weight, current half bounce
  real(fp) :: w_fds = zero
        ! FDS weight, current half bounce
  logical :: l_fd0 = .FALSE.
        ! .TRUE. if w_fd0 .gt. ZERO
  logical :: l_fdb = .FALSE.
        ! .TRUE. if w_fdb .gt. ZERO
  logical :: l_fdp = .FALSE.
        ! .TRUE. if w_fdp .gt. ZERO
  logical :: l_fdq = .FALSE.
        ! .TRUE. if w_fdq .gt. ZERO
  logical :: l_fdr = .FALSE.
        ! .TRUE. if w_fdr .gt. ZERO
  logical :: l_fds = .FALSE.
        ! .TRUE. if w_fds .gt. ZERO
  integer :: nwx_fd = 0
        ! number of nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fd0
        ! FD0 weights, nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fdb
        ! FD0 weights, nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fdp
        ! FD0 weights, nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fdq
        ! FD0 weights, nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fdr
        ! FD0 weights, nearby diffusivity tests
  real(fp), dimension(:), allocatable :: wx_fds
        ! FD0 weights, nearby diffusivity tests
  logical :: lx_fd0 = .FALSE.
        ! .TRUE. if any of wx_fd0(...) .gt. ZERO
  logical :: lx_fdb = .FALSE.
        ! .TRUE. if any of wx_fdb(...) .gt. ZERO
  logical :: lx_fdp = .FALSE.
        ! .TRUE. if any of wx_fdp(...) .gt. ZERO
  logical :: lx_fdq = .FALSE.
        ! .TRUE. if any of wx_fdq(...) .gt. ZERO
  logical :: lx_fdr = .FALSE.
        ! .TRUE. if any of wx_fdr(...) .gt. ZERO
  logical :: lx_fds = .FALSE.
        ! .TRUE. if any of wx_fds(...) .gt. ZERO
  real(fp) :: dksgc = zero
        ! ACCUMULATED DELTA(KSID) (CONVECTIVE) IN 1/2 ORBIT
  real(fp) :: dksgd = zero
        ! ACCUMULATED DELTA(KSID) (DIFFUSIVE) IN 1/2 ORBIT
  real(fp) :: dlbxdtk = zero
        ! DISTANCE TRAVELED DOWN BOX AXIS (L) PER DISTANCE
  real(fp), dimension(:), allocatable :: dravizc
        ! FLUX SURFACE AVG ** 1/DR **
  real(fp) :: dtgcm = zero
        ! BEAM COMPRESSION TIMESTEP
  real(fp) :: dtgcx = zero
        ! GOOSED DT FOR CURRENT STEP ON BEAM CX CALCULATION
  real(fp) :: dtgcx_sum = zero
        ! GOOSED DT FOR CURRENT STEP ON BEAM CX CALCULATION
  real(fp) :: dtgfp = zero
        ! GOOSED DT FOR CURRENT STEP ON BEAM FP CALCULATION
  real(fp) :: dtnv = zero
        ! ORBIT TIMESTEP NORMALIZATION FACTOR
  real(fp) :: dvbox = zero
        ! volume of single zone in box grids
  real(fp) :: dvg = zero
        ! ACCUMULATED DELTA(V) ALONG A HALF-ORBIT
  real(fp), dimension(:), allocatable :: dvolzc
        ! volume of zone i (cm**3)
  real(fp), dimension(:), allocatable :: dvolzc_r8
        ! volume of zone i (cm**3) (real(fp))
  real(fp), dimension(:), allocatable :: dareazc
        ! cross sectional area of zone i (cm**2)
  real(fp), dimension(:), allocatable :: dareazc_r8
        ! cross sectional area of zone i (cm**2) (real(fp))
  real(fp) :: dxdlbox = zero
        ! RATE OF CHANGE OF X MOVING DOWN BOX AXIS (L)
  real(fp), dimension(:,:), allocatable :: dxi
        ! DIFFERENTIAL ELEMENTS OF XI
  real(fp) :: dxitk = zero
        ! xi-spacing of 2d grid zone rows
  real(fp), dimension(:), allocatable :: dxizb
        ! DIFFERENTIAL ELEMENTS OF XI
  real(fp), dimension(:), allocatable :: dxizc
        ! DIFFERENTIAL ELEMENTS OF XI
  real(fp) :: dydlbox = zero
        ! RATE OF CHANGE OF Y MOVING DOWN BOX AXIS (L)
  real(fp), dimension(:,:), allocatable :: eb2plp
        ! 2D BEAM PLL AVG ENERGY, PREV. TIMESTEP
  real(fp), dimension(:,:), allocatable :: eb2ppp
        ! 2D BEAM PERP AVG ENERGY, PREV. TIMESTEP
  real(fp), dimension(:,:), allocatable :: eba2pl
        ! <Epll>
  real(fp), dimension(:,:), allocatable :: eba2pp
        ! <Eperp>
  real(fp), dimension(:), allocatable :: ebdlina
  real(fp), dimension(:), allocatable :: ebdmaxa
  real(fp) :: ebdmaxm = zero
        ! Emax for current orbiting specie
  real(fp), dimension(:,:), allocatable :: efbm
        ! BEAM DIST. FCN ENERGY GRIDS, BY SPECIES.
  real(fp), dimension(:,:), allocatable :: efbmb
        ! BEAM DIST. FCN ENERGY ZONE BDYS, BY SPECIES.
  real(fp) :: einj = zero
  real(fp), dimension(:), allocatable :: einjs
        ! nominal injection energy
  real(fp), dimension(:), allocatable :: vmins
        ! reference minimal velocity = 1/100'th
        ! the nominalinjection velocity sqrt(einjs(j)/v2toev)
  real(fp), dimension(:), allocatable :: eirefs
        ! reference constrained injection energy
        !   min(200.0e3,einjs(j))
  real(fp), dimension(:), allocatable :: xmrefs
        ! reference injection momentum
        ! xmrefs(j)=Raxis*abeams(j)*1.6726e-24*sqrt(eirefs(j)/v2toev)
  real(fp) :: eldot = zero
        ! (dphilim/dt)/(2*philim)
  real(fp), dimension(:,:,:), allocatable :: en0mn
        ! NEUTRAL temperature OF HYDROGENIC SPECIES DUE TO SOURCE
  real(fp) :: ep = zero
        ! ENERGY OF MONTE-CARLO BEAM ION
  real(fp) :: epbdmx = zero
        ! MAX (V*V) IN BEAM DIST. FCN. (CM**2/SEC**2)
  real(fp) :: ephflr = zero
        ! BEAM ION ELECTROSTATIC POTENTIAL AT FLR
  real(fp) :: ephion = zero
        ! ELECTROSTATIC POTENTIAL OF ION G.C. RELATIVE TO
  real(fp), dimension(:,:), allocatable :: evol0
        ! AVG energy OF THERM. vol SCE NEUTRALS
  real(fp), dimension(:,:), allocatable :: evol02
        ! energy vol sce neutrals
  real(fp), dimension(:,:), allocatable :: ewall0
        ! AVG energy THERM. WALL SCE NEUTRALS
  real(fp) :: fac_e = zero
        ! collision coefficients energy variation factor (fokker & collid)
  real(fp), dimension(:,:,:,:), allocatable :: fbm
        ! the guiding center(GC) distribution function f(E,vpll/v,position,species)
        ! -- 2d spatial grid
        ! lab frame E, vpll/v
  real(fp), dimension(:,:,:,:,:), allocatable :: fbm_ptcl
        ! particle  distribution function f(E,vpll/v,position,species)
        ! -- 2d spatial grid
        ! lab frame E, vpll/v
  real(fp), dimension(:,:,:,:), allocatable :: fbm_cx
        ! the guiding center(GC) CX fast neutral source distribution function f(E,vpll/v,position,species)
        ! -- 2d spatial grid in plasma +SOL
        ! lab frame E, vpll/v
  real(fp), dimension(:,:), allocatable :: fbz
        ! ACTUAL/EXTERNAL TOR. B FIELD AT ZN/BDY J,IZ
  real(fp), dimension(:,:), allocatable :: fppray
        ! FP COLLISIONS FREQUENCY CONTROL ARRAY
  real(fp) :: fpskip = zero
        ! *BEAM: SKIP FACTOR ON BEAM ION FOKKER-PLANCK
  real(fp) :: dt_adjust = zero
        ! current timestep adjustment factor.
  real(fp) :: dt_adj0 = zero
        ! initial timestep adjustment factor (if species
        ! specific data unavailable).
  real(fp), dimension(:), allocatable :: dt_adj0s
        ! per species timestep adjustment factor
  real(fp) :: dt_adj0_sum = zero
  integer :: num_dt_adj0 = 0
        ! items used for update of dt_adj0 (cf orball)
  real(fp), dimension(:,:), allocatable :: goosay
        ! STORED GOOSE FACTOR FOR EACH M.C. PTCL.
  real(fp) :: goose = zero
        ! GOOSING FACTOR FOR CURRENT M.C. BEAM ION
  real(fp) :: gzfact = zero
        ! global variable to store zfact in nbstart
  real(fp), dimension(:,:), allocatable :: vtaub
        ! estimated 0-banana-width v*tau(bounce) for trapped orbits
        ! with Brefl = Bmin + 0.25*(Bmax-Bmin)  (1)
        !      Brefl = Bmin + 0.50*(Bmax-Bmin)  (2)
        !      Brefl = Bmin + 0.75*(Bmax-Bmin)  (3)
        !      Brefl = Bmin + 1.00*(Bmax-Bmin)  (4)  approximately
        ! (see bparms.for) -- xplasma routine is used for evaluation.
        ! this is to check goosing; normally goosing is updated every
        ! half bounce; if a full bounce's worth of time passes with no
        ! update (according to estimate), re-estimate goosing early.
        ! this will allow off-midplane trapped orbits, as seen on NSTX
        ! to have goosing updates.  dmc Oct 2003.
        !   estimated tau_bounce = vtaub/vion
        ! these estimates are computed on each surface (except the
        ! singular axis); each zone is assigned the maximum result
        ! from its bounding surfaces; Bmin & Bmax at outer surface of
        ! zone are used for reference (cf orbit.for)
  real(fp), dimension(:), allocatable :: vtaubmin
        ! minval(vtaub(1:4,mj)) -- saved for speed optimization
        !for the following 6 quantities suffix zc used for "zone centers"
  real(fp), dimension(:), allocatable :: gb1zc
        ! <B> FLUX SURF VOLUME AVG
  real(fp), dimension(:), allocatable :: gb2zc
        ! <B**2> FLUX SURF VOLUME AVG
  real(fp), dimension(:), allocatable :: gncfhzc
        ! <function of h>  FLUX SURF VOLUME AVG
  real(fp), dimension(:), allocatable :: gr2zc
        ! <R**2> FLUX SURF VOLUME AVG
  real(fp), dimension(:), allocatable :: grizc
        ! <1/R> FLUX SURF VOLUME AVG
  real(fp), dimension(:), allocatable :: gr2izc
        ! <1/R**2> FLUX SURF VOLUME AVG
  real(fp), dimension(:,:), allocatable :: gxtabl
  real(fp), dimension(:,:,:,:,:), allocatable :: hexbbcx0
        ! . EXCITED STATES EFFECT ON BEAM-BEAM CX
  real(fp), dimension(:,:,:), allocatable :: hexbbxrc
        ! BEAM-BEAM CX EXCITED STATES ENHANCEMENT
  real(fp), dimension(:,:,:), allocatable :: hexbbxrp
        ! BEAM-BEAM CX EXCITED STATES ENHANCEMENT
  real(fp), dimension(:,:,:), allocatable :: hvsigi
        ! N*SIGMA*V IN PLASMA ZONE J FOR HEAVY ION SPECIES
  integer :: init_xplasma = 0
        ! internal control, indicating whether xplasma object needs
        ! to be initialized from plasma-state (0 means "no").
  integer :: write_xplasma = 0
        ! internal control, indicating whether xplasma object needs
        ! to be written to file when output data processing is completed.
        ! (0 means "no").
  integer :: id_R = 0
        ! R(xi,theta) xplasma function id
  integer :: id_Z = 0
        ! Z(xi,theta) xplasma function id
  integer :: id_bphi = 0
        ! B field component functions
  integer :: id_br = 0
        ! B field component functions
  integer :: id_bz = 0
        ! B field component functions
  integer :: id_phi = 0
        ! electrostatic potential function
  integer :: id_psi = 0
        ! poloidal flux function
  integer :: ldep = 0
        ! NUMBER OF BOX OR TUBE CURRENTLY BEING PROCESSED
  integer :: lcentr = 0
        ! ZONE INDEX OF INNERMOST REAL ZONE
  integer :: lcm1 = 0
        ! LCENTR-1
  integer :: lcp1 = 0
        ! LCENTR+1 (=INDEX TO OUTER BDY OF INNERMOST ZONE)
  logical :: ldtnfr = .FALSE.
        ! timestep adjust switch
  integer :: ledge = 0
        ! ZONE INDEX OF OUTERMOST REAL ZONE
  integer :: lep1 = 0
        ! LEDGE+1 (=INDEX TO PLASMA OUTER BDY)
  integer :: lrunid = 0
        ! length of runid
  integer :: lsbeam = 0
        ! FAST ION SPECIES CURRENTLY BEING "FOLLOWED"
  integer :: minbs_expnd = 0
        ! array dimension for slitting at the begging minbs_expnd=minb
        ! if needed it will be changed in nutrav and splitting array will be
        ! expanded
  integer, dimension(:), allocatable :: lsbeama
        ! lsbeam index corresponding to each injected beam
  integer, dimension(:), allocatable :: lsbeamr
        ! lsbeam index for fast reagent, each beam-target reaction
  real(fp) :: alpha_gyro = zero
        !alpha_gyro poloidal aiming angle
        !   0  -> normal to mag. surfaces in an outward direction;
        !   pi -> normal to mag. surfaces in an intward direction:
        !   pi/2 -> tangential to mag. surfaces, "left".
        !  3pi/2-> tangential to mag. surfaces, "right"
        !                   ^ alpha_gyro
        !          |       /
        !           |     /^  positive direction
        !            |   /  ^
        !             | /    ^
        !             |/______^__ >PSI
        !             |
        !            |
        !           |
        !          |
  integer :: lgyro = 0
        ! lgyro index corresponding to gyro angle (see mgyro)
  integer :: mrfdeut = 0
  integer :: mrfhel3 = 0
  integer :: mrfhel4 = 0
  integer :: mrfhydr = 0
  integer, dimension(:), allocatable :: mrfmini
  integer :: mrftrit = 0
  integer :: mzones = 0
        ! # OF RADIAL ZONES, INCLUDING DUMMY ZONES
  integer, dimension(:), allocatable :: nadrbox
        ! box data addresses
  integer, dimension(:,:), allocatable :: nbbsp
        ! PTR TO BEAM-BEAM REAGENT DISTRIBUTIONS
  integer :: nbdeut = 0
  integer :: nbflg = 0
        ! FLAG INDICATING ION STATUS: E.G. THERMALIZED
        ! NBFLG=1 MEANS AN ION WHICH HAS SURVIVED
        ! FULL BEAM CODE TIMESTEP, AND IS STILL ORBITING
        ! NBFLG=2 ION IS THERMALIZED (collid)
        ! NBFLG=3 MEANS AN ION LOST ON A BAD ORBIT, BY C-X, BY RIPPLE,
        ! OR BY F-BONES
  integer :: nbhel3 = 0
  integer :: nbhel4 = 0
  integer :: nbhydr = 0
  logical :: nbi_init_flag = .FALSE.
        ! .TRUE. if module has been initialized
  logical :: need_read_ptcl = .FALSE.
        ! control switch to synchronize reading field state file
        ! by master processor only
        ! and ptcl state files by all processors. If .TRUE. -
        ! module has been initialized then  ptcl.state file
        ! need to be read in nubeam_step_child
  logical :: sawflag2 = .FALSE.
        ! Subroutine SAWNBI settable flag requesting sawtooth mixing
  integer :: nbicx = 0
        ! FLAG FOR CH-EX BEAM NEUT. STATUS,EG ESCAPED
  integer, dimension(:,:), allocatable :: nbndex
        ! M.C. BEAM PTCL STORAGE INDEXING ARRAY
  integer, dimension(:), allocatable :: nbouta
        ! workspace for NB.OUT routines
  integer :: nbs = 0
        ! SPECIES INDEX OF BEAM SPECIES
  integer, dimension(:), allocatable :: nbsbsa
        ! INDEX TO BEAM ION SCE, CX SPLITTEE
  integer, dimension(:), allocatable :: nbebsa
        ! INDEX TO BEAM ION energy fraction, CX SPLITTEE
  real(fp), dimension(:), allocatable :: einjbsa
        ! lab frame deposition energy, CX SPLITTEE
  real(fp), dimension(:), allocatable :: pinjbsa
        ! total power injected at MC partilce deposition without shine through power
        ! CX SPLITTEE
  integer, dimension(:,:), allocatable :: nbscay
        ! INDEX TO BEAM SOURCE FROM WHICH ION ORIGINATED
  integer, dimension(:,:), allocatable :: nbienay
        ! beam energy fraction (1=full, 2=half, 3=third) of original ion
  real(fp), dimension(:,:), allocatable :: einjay
        ! injection energy of original ion (lab frame)
  integer :: nbsce = 0
        ! SOURCE INDEX OF CURRENT ORBITTING BEAM ION
  integer :: nbien = 0
        ! ENERGY FRACTION INDEX OF CURRENT ORBITTING BEAM ION
  real(fp) :: einj0 = zero
        ! lab frame energy of injected neutral or newborn fusion product ion
  real(fp) :: pinj_ay = zero
        ! total power injected at MC partilce deposition without shine through power
  integer :: nbsplt = 0
        ! CURRENT # OF CH-EX BEAM NEUTRAL SPLITTEES
  integer :: nbtrit = 0
  integer :: ndeut = 0
  integer, dimension(:,:), allocatable :: nfbmpt
        ! NO. OF NON-ZERO PTS IN FAST ION DISTRIBUTIONS
  integer :: nfbzns = 0
        ! total no. of zones
  integer :: nfbznsi = 0
        ! total no. of zones inside plasma
  integer :: nfhel3 = 0
  integer :: nfhel4 = 0
  integer :: nflr = 0
        ! TRANSP ZONE NUMBER OF CURRENT BEAM ION (G.C.)
  integer :: nflr2 = 0
        ! BEAM (2D) SPATIAL ZONE NO. OF CURRENT ION (GC)
  integer :: nflrx = 0
        ! ZONE INDEX (INCLUDING BLOATED ZONES) (GC)
  integer :: n_gyro = 0
        ! poloidal aiming grid index OF CURRENT ION
  integer :: nfprot = 0
  integer :: nfprot14 = 0
  integer :: nftrit = 0
  integer :: ngc = 0
        ! TRANSP ZONE NUMBER OF CURRENT BEAM ION (FLR)
  integer :: ngc2 = 0
        ! BEAM (2D) SPATIAL ZONE NO. OF (FLR) ION
  integer :: ngcx = 0
        ! ZONE INDEX (INCLUDING BLOATED ZONES) (FLR)
  integer :: ngoose = 0
        ! COUNTER FOR CALCULATING PTCL GOOSE FACTOR
  integer :: nhel3 = 0
  integer :: nhel4 = 0
  integer :: nhydr = 0
  integer, dimension(:,:), allocatable :: njavay
        ! M.C. BEAM PTCL ORBIT-AVG J (ZONE INDEX) ARRAY
  integer, dimension(:,:), allocatable :: njavay_saw
        ! M.C. BEAM PTCL ORBIT-AVG J (ZONE INDEX) ARRAY
        ! post-sawtooth values
  logical, dimension(:), allocatable :: nlbba
        ! .TRUE. IF A CERTAIN BEAM-BEAM REACTION OCCURS
  logical :: nlddp = .FALSE.
        ! if DD(p) fusion is present
  logical :: nlfbm = .FALSE.
        ! Fb flag, always .TRUE.
  logical, dimension(:), allocatable :: nlfbmfpp
  logical, dimension(:), allocatable :: nlfprod
  logical :: nlfusn = .FALSE.
  logical, dimension(:), allocatable :: nlheavy
  integer :: nlith = 0
  logical :: nlmcany = .FALSE.
        ! if there are any MC species
  logical :: nlpout = .FALSE.
        ! IF PTCL IS OUTSIDE THE PLASMA
  integer, dimension(:), allocatable :: nmmbinv
        ! ARRAY OF NO. OF MOMENTS TO USE - 1/MOD(B)
  integer, dimension(:), allocatable :: nmmjac
        ! ARRAY OF NO. OF MOMENTS TO USE - JACOBIAN
  integer :: nout = 0
  integer :: nphipc = 0
        ! total number of phi(p) evaluation surfaces
  integer :: nreact = 0
        ! no. of known reactions (cf rdbbsg.for)
  integer :: nrfdeut = 0
  integer :: nrfhel3 = 0
  integer :: nrfhel4 = 0
  integer :: nrfhydr = 0
  integer, dimension(:), allocatable :: nrfmini
  integer :: nrftrit = 0
  integer :: nsbeam = 0
  integer :: nsfast = 0
  integer :: nsjdotb = 0
        ! DIRECTION OF TOROIDAL B WRT. PLASMA CURRENT
        ! +1 if Ip and Bt are parallel, -1 if anti-pll
  integer :: nsnccw = 0
        ! +1 if Bphi is CCW looking down on tokamak from above
        ! -1 otherwise
  integer :: nthzons = 0
  integer, dimension(:), allocatable :: nthzsm
        ! # theta zones total inside indicated row
  integer :: ntrit = 0
  integer, dimension(:), allocatable :: nuse
        ! # OF M.C. PTCLS CURRENTLY IN STORAGE
        ! on a given processor (MPI related)
  integer :: J_SPL = 0
        ! MC particles splitting during orbiting  in orball, on each CPU
        ! (MPI related)
  real(fp), dimension(:), allocatable :: puseyp
        ! Total # of MC ptcls currently in storage
        ! is using for a global control during deposition
        ! and orbiting
  integer, dimension(:), allocatable :: ntotptp
        ! total # of M.C. slots on a
        ! given processor
  integer, dimension(:), allocatable :: nnewp
        ! # OF M.C. NEW PTCLS
        ! on a given processor (MPI related)
  integer :: nwrite_global = 0
        ! (>0 would write global particle lists)
        ! (MPI related)
  integer :: nxskpb = 0
        ! # of xi zones per 2d grid zone row
  integer :: nxtbzn = 0
        ! set in BMZINI now, from XPLASMA info.
  integer, dimension(:), allocatable :: nznblea
  integer :: nznblo = 0
        ! no. of xi radial zones in extrapolated space
  integer, dimension(:), allocatable :: nznbmea
  integer :: nznbma_save = 0
        ! saved no. of pitch angle zones (not affected by shrink/rbloat)
  integer :: nznbmr = 0
        ! total no. of xi-zone rows
  integer :: nznbmri = 0
        ! no. of zone rows inside plasma
  integer :: nthzsm0 = 0
        ! no. of poloidal zones in innermost row
  integer :: nzp1 = 0
        ! NUMBER OF ACTIVE ZONES, PLUS ONE
  real(fp), dimension(:), allocatable :: omegag
        ! thermal plasma ROTATION toroidal ANGULAR VELOCITY (RADIANS/SEC)
  real(fp), dimension(:), allocatable :: omega_ExB
        ! ExB ANGULAR VELOCITY (RADIANS/SEC)
        !   satisfies d[electrostatic potential]/dx = -omega_ExB*x*TFLUX/(q*pi)
        !   it provides a frame transformation to recapture vpll=0 matching to
        !   trapped orbit bounce point correspondence, when a radial electro-
        !   static potential is present
  real(fp), dimension(:,:,:), allocatable :: omg0mn
        ! NEUTRAL rotation velocity OF HYDROGENIC SPECIES DUE TO SOURCE
  real(fp), dimension(:,:), allocatable :: ovol0
        ! AVG ROTATION OF THERM. VOL.SCE NEUTRALS
  real(fp), dimension(:,:), allocatable :: ovol02
        ! 2D OVOL0 (AVG ROTATION VOL SCE THERM NEUTRALS)
  real(fp), dimension(:,:), allocatable :: owall0
        ! AVG ROTATION OF THERM. WALL SCE NEUTRALS
  real(fp), dimension(:), allocatable :: phibsce
        ! toroidal angle (phi) location
  real(fp), dimension(:), allocatable :: phiprg
        ! Radial electrostatic potential (volts)
  real(fp), dimension(:), allocatable :: plflxg
        ! psi(x)
  real(fp), dimension(:), allocatable :: plflxi
        ! poloidal flux from iota(bar) integral, used in MC code
  real(fp), dimension(:,:), allocatable :: pmechay
        ! Pmech: mechanical mo.
  real(fp), dimension(:,:), allocatable :: pphiay
        ! Pphi: canonical momentum
  real(fp), dimension(:,:), allocatable :: pinjay
        ! total power injected at MC partilce deposition without shine through power
  real(fp), dimension(:,:), allocatable :: pmechay_saw
        ! Pmech: mechanical mo. -- post-sawtooth
  real(fp), dimension(:,:), allocatable :: pphiay_saw
        ! Pphi: canonical momentum -- post-sawtooth
  real(fp), dimension(:), allocatable :: qgeo
        ! VALUE OF Q AT BOUNDARY J-- PRIMARY DEFINITION
  real(fp) :: raxis = zero
        ! RADIUS OF MAGNETIC AXIS (CM)
  real(fp), dimension(:), allocatable :: rbsa
        ! CH-EX NEUT. SPLITTEES RADIAL COORDS.
  real(fp), dimension(:,:,:), allocatable :: rhbs
        ! BEAM DENSITY SPECIES IBS - CURRENT TI STEP
  real(fp), dimension(:,:), allocatable :: rhi
        ! TOTAL IMPURITY DENSITY AT BDY/ZONE J
  real(fp), dimension(:,:,:), allocatable :: rhimom
        ! MOMENTS OF THE IMPURITY DENSITY WRT
  real(fp), dimension(:,:,:), allocatable :: rhix
        ! MULTIPLE IMPURITY DENSITY AT BDY/ZONE J USED IN CODE
  real(fp), dimension(:,:,:), allocatable :: rhmin
        ! DENSITY OF ICRF MINORITY ION SPECIES (N/CM**3)
  real(fp) :: rho_axis = zero
        !
  real(fp) :: rho_bdy = zero
        !
  real(fp), dimension(:,:,:), allocatable :: rhob
        ! BACKGROUND PLASMA DENSITY AT BDY/ZONE J
  real(fp), dimension(:,:), allocatable :: rhoel
        ! ELECTRON DENSITY AT BDY/ZONE J
  real(fp) :: rlar = zero
        ! LARMOR RADIUS OF CURRENT BEAM ION
  real(fp) :: rmassb = zero
        ! MASS, GRAMS OF CURRENT FAST ION SPECIES
  real(fp) :: rmjcur = zero
        ! last inside-the-plasma radius value
  real(fp) :: rmjflr = zero
        ! BEAM ION ACTUAL MAJOR RADIUS (NOT G.C.)
  real(fp) :: rmjion = zero
        ! Rmajor OF BEAM ION G.C.
  real(fp), dimension(:,:), allocatable :: rmjionay
        ! R of GC, prev. step
  real(fp), dimension(:,:), allocatable :: rmjionay_saw
        ! R of GC, prev. step -- post-sawtooth
  real(fp), dimension(:,:), allocatable :: rmjrmp
        ! midplane major radius (extended grid)
  real(fp), dimension(:,:), allocatable :: rmnrmp
        ! midplane minor radius (extended grid)
  real(fp), dimension(:,:), allocatable :: sfusion
        ! fusion product source
  real(fp) :: sinpho = zero
        ! SIN OF ORBIT GYRO PHASE ANGLE
  real(fp), dimension(:,:,:,:), allocatable :: sn00i2
        ! SPATIALLY 2D BEAM-BEAM IONIZATION RATE TABLE
  real(fp), dimension(:,:,:,:), allocatable :: sn00x2
        ! SPATIALLY 2D BEAM-BEAM CX RATE TABLE neutral beam IB,
        ! NB*sigma*v
  real(fp), dimension(:,:,:), allocatable :: snb0i2
        ! neutral sink, fast ion ii
  real(fp), dimension(:,:,:), allocatable :: snbn02
        ! neutral sink, fast ion cx
  real(fp), dimension(:,:,:,:,:), allocatable :: svn3vii_h
        ! BEAM fast ion - BEAM hydrogen isotopes neutral
        ! IONIZATION RATE TABLE with <E_0>, <v0_II/v0>, alpha_gyro, fast ion species
  real(fp), dimension(:,:,:,:,:), allocatable :: svn3vii_he
        ! BEAM fast ion - BEAM helium isotopes neutral  (if any)
        ! IONIZATION RATE TABLE  <E_0>, <v0_II/v0>, alpha_gyro, fast ion species
  real(fp), dimension(:,:,:,:,:), allocatable :: svn3vcx_h
        ! BEAM fast ion - BEAM hydrogen isotopes neutral
        ! CX RATE TABLE with <E_0>, <v0_II/v0>, alpha_gyro, fast ion species
  real(fp), dimension(:,:,:,:,:), allocatable :: svn3vcx_he
        ! BEAM fast ion - BEAM helium isotopes neutral  (if any)
        ! CX RATE TABLE  <E_0>, <v0_II/v0>, alpha_gyro, fast ion species
  real(fp), dimension(:,:), allocatable :: splg
        ! MHD G SPLINE COEFF
  real(fp), dimension(:,:), allocatable :: splimhd
        ! MHD I(xi) SPLINE COEFF TIME DERIVATIVES
  real(fp), dimension(:,:), allocatable :: spliot
        ! MHD iota(bar) SPLINE COEFF
  real(fp), dimension(:,:), allocatable :: splphi
        ! MHD RADIAL ELECTRIC POTENTIAL
  real(fp), dimension(:,:,:), allocatable :: srbn02
        ! thermal neutral source
  real(fp), dimension(:,:,:), allocatable :: srbn02_dep
        ! thermal neutral source (deposition calculation only)
  real(fp) :: srd = zero
        ! DELTA/DELTA_S :GWB RIPPLE LOSS CRITERION
  real(fp), dimension(:), allocatable :: tbsa
        ! CH-EX NEUT. SPLITTEES TIME COORDINATE
  real(fp), dimension(:,:), allocatable :: tdepay
        ! STORED M.C. BEAM PTCL DEPOSITION TIME
  real(fp), dimension(:,:), allocatable :: te
        ! ELECTRON TEMPERATURE AT BDY/ZONE J
  real(fp) :: tflux = zero
        ! TOROIDAL FLUX ENCLOSED IN PLASMA
  real(fp), dimension(:,:), allocatable :: tfluxay
        ! actual flux surface (Wb)
  real(fp), dimension(:,:), allocatable :: tfluxay_saw
        ! actual flux surface (Wb) -- post-sawtooth
  real(fp) :: tfluxp = zero
        ! past timestep inclosed flux
  real(fp) :: th = zero
        ! M.C. ION THETA, POSITION COORDINATE
  real(fp), dimension(:,:), allocatable :: thay
        ! STORED M.C. PTCL ANGLE COORDINATES
  real(fp), dimension(:,:), allocatable :: thay_saw
        ! STORED M.C. PTCL ANGLE COORDINATES (post-sawtooth values)
  real(fp) :: thbdy0 = zero
        ! [0,pi] or [-pi,pi] depending on symmetry
  real(fp) :: thbdy1 = zero
        ! [0,pi] or [-pi,pi] depending on symmetry
  real(fp), dimension(:,:), allocatable :: thbmxm
        ! (XJADEP theta(Bmax) info)
  real(fp), dimension(:), allocatable :: thbsa
        ! CH-EX NEUT. SPLITTEES ANGLE COORD.
  real(fp) :: thcur = zero
        ! last inside-the-plasma theta value
  real(fp) :: thflr = zero
        ! FLR "THETA" LOCATION OF M.C. ION
  real(fp), dimension(:,:), allocatable :: thntns
        ! thermonuclear source
  real(fp) :: thsta = zero
        ! starting MHD angle coord for straight line
  real(fp), dimension(:), allocatable :: thzons
        ! th grid (cf bmzini)
  real(fp), dimension(:), allocatable :: xiblo
        ! bloated xi grid (cf bmzini)
  real(fp), dimension(:,:), allocatable :: ti
        ! ION TEMPERATURE AT BDY/ZONE J (eV)
  real(fp), dimension(:,:,:), allocatable :: timin
        ! TEMPERATRE OF ICRF MINORITY ION SPECIES (eV)
  real(fp), dimension(:), allocatable :: tmjsm
        ! time smoothed majority temperature.
  real(fp) :: tscx = zero
        ! ORBITTING ION CH-EX PROBABILITY INTEGRAL
  real(fp) :: ttcx = zero
        ! ORBITTING ION CH-EX PARAMETER
  real(fp) :: v2toev = zero
        ! CONVERSION FACTOR
  real(fp), dimension(:,:), allocatable :: vay
        ! STORED M.C. PTCL VELOCITIES
  real(fp), dimension(:,:), allocatable :: vay_saw
        ! STORED M.C. PTCL VELOCITIES  -- post-sawtooth
  real(fp), dimension(:), allocatable :: vbsa
        ! CH-EX NEUT SPLITTEES VELOCITIES
  real(fp), dimension(:,:), allocatable :: vbtor2
        ! <Vphi>
  real(fp), dimension(:,:), allocatable :: vbtr2p
        ! VBTOR2, PREV. STEP, FOR NUTRAV BB STOPPING
  real(fp), dimension(:,:), allocatable :: vcrit
        ! V(CRIT) BY POSITION, ION SPECIES CF FOKKER PLANCK
  real(fp), dimension(:,:), allocatable :: vcrit3
        ! V(CRIT)**3 BY POSITION, ION SPECIES CF FOKKER PLANCK
  real(fp), dimension(:,:,:,:), allocatable :: vcxbn0
        ! 2D 1ST GEN BEAM NEUTRAL V0 DIR.COSINES
  real(fp) :: vflr = zero
        ! BEAM ION SPEED AT FLR
  real(fp), dimension(:), allocatable :: viflr
        ! BEAM ION VELOC VECTOR (PHI,RHO,THETA) AT FLR
  real(fp) :: vion = zero
        ! M.C. BEAM ION VELOCITY (CM/SEC)
  real(fp) :: vioncur = zero
        ! last inside-the-plasma velocity (VION)
  real(fp) :: viplasa = zero
        ! BEAM ION VELOCITY, PLASMA FRAME, BOUNCE AVG
  real(fp), dimension(:,:), allocatable :: vmin3
        ! V(THERM)**3 BY POSITION, ION SPECIES
  real(fp), dimension(:,:), allocatable :: voln0
        ! AVG NEUTRAL DENSITY OF THERM. VOL SCE NEUTRALS
  real(fp), dimension(:,:), allocatable :: voln02
        ! 2D VOL. SCE THERM. NEUTRAL DENSITIES.
  real(fp) :: vperp = zero
        ! M.C. ION PERP VELOCITY (CM/SEC)
  real(fp) :: vplflr = zero
        ! BEAM ION PARALLEL VELOCITY AT FLR
  real(fp) :: vpll = zero
        ! M.C. ION PLL VELOCITY (CM/SEC)
  real(fp) :: vpllcur = zero
        ! last inside-the-plasma pll velocity (VPLL)
  real(fp), dimension(:), allocatable :: vpoh
        ! LOOP VOLTAGE as used for POH calculation, zone J
  real(fp) :: vppflr = zero
        ! BEAM ION PERP VELOCITY AT FLR
  real(fp), dimension(:,:,:), allocatable :: vpvn00
        ! AVERAGE VPLL/V OF 1D 1ST GEN BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: walln0
        ! AVG NEUTRAL DENSITY OF THERM. WALL SCE NEUTRALS
  real(fp), dimension(:,:), allocatable :: walln02
        ! 2D Wall SCE THERM. NEUTRAL DENSITIES.
  real(fp), dimension(:,:), allocatable :: walln0x
        ! neutral density, fixed on 2d grid zone rows
  real(fp) :: wght = zero
        ! M.C. WEIGHT OF CURRENT BEAM ION
  real(fp), dimension(:,:), allocatable :: wghtay
        ! STORED WEIGHT OF ALL M.C. PTCLS
  real(fp), dimension(:), allocatable :: wtbsa
        ! WEIGHT OF CH-EX NEUT. SPLITTEES
  real(fp), dimension(:,:), allocatable :: wtrap
        ! BTRAP normalizer
  real(fp), dimension(:,:), allocatable :: wtrap0
        ! BTRAP0 normalizer
  real(fp), dimension(:,:,:), allocatable :: wtrape
        ! BTRAPE normalizer
  real(fp) :: xbgrid = zero
        ! STARTING X DISPLACEMENT RELATIVE TO CENTERLINE OF BEAM
  real(fp), dimension(:,:), allocatable :: xbjaco
        ! 2x2 Jacobian
  real(fp) :: xminbm = zero
        ! boundary of xi-space
  real(fp), dimension(:,:), allocatable :: xcrb2c
        ! unshielded XCURB2
  real(fp), dimension(:,:), allocatable :: xcurb2
        ! driven current on 2d grid
  real(fp), dimension(:,:), allocatable :: xi
        ! FIXED MESH COORDINATE (NORMALIZED FLUX LABEL)
  real(fp), dimension(:,:), allocatable :: xi_r8
        ! FIXED MESH COORDINATE (NORMALIZED FLUX LABEL) (real(fp))
  real(fp), dimension(:,:), allocatable :: xiay
        ! STORED M.C. PTCL "XI" COORDINATES
  real(fp), dimension(:,:), allocatable :: xiay_saw
        ! STORED M.C. PTCL "XI" COORDINATES -- possible post-sawtooth
  real(fp) :: xibado = zero
        ! bad orbit error detection parameter
  real(fp) :: xicur = zero
        ! last inside-the-plasma XION value
  real(fp) :: xiflr = zero
        ! FLR "XI" LOCATION OF M.C. ION
  real(fp) :: xiliml = zero
        ! xi inside of which limiter loss detection not needed
  real(fp) :: xinflr = zero
        ! FLR INTERPOLATION FACTOR IN XI SPACE
  real(fp) :: xinflrx = zero
        ! DITTO, OUT INTO BLOATED ZONES
  real(fp) :: xingc = zero
        ! GC INTERPOLATION FACTOR IN XI SPACE
  real(fp) :: xingcx = zero
        ! DITTO, OUT INTO BLOATED ZONES
  real(fp) :: xion = zero
        ! XI
  real(fp), dimension(:), allocatable :: xisbsa
        ! CX RECAPTURE XI (ON FLR GYRO) USED IN RADIAL CURRENT ESTIMATE
  real(fp) :: xistate = zero
        ! ION XI AT START OF ORBIT TRACK (DEP OR ORB AVG)
  real(fp), dimension(:,:), allocatable :: xistay
        ! DEPOSITED XIFLR OR PREV. ORB. AVG. XION FOR MC ION USED
  real(fp), dimension(:,:), allocatable :: xistay_saw
        ! DEPOSITED XIFLR OR PREV. ORB. AVG. XION FOR MC ION USED
        ! post-sawtooth profile
  real(fp) :: xistop = zero
        ! ION XI AT END OF ORBIT TRACK (ORB AVG OR LOST)
  real(fp), dimension(:), allocatable :: xizb
        ! FIXED MESH COORDINATE (Zone Boundary)
  real(fp), dimension(:), allocatable :: xizb_r8
        ! FIXED MESH COORDINATE (Zone Boundary) (real(fp))
  real(fp), dimension(:), allocatable :: xizc
        ! FIXED MESH COORDINATE (Zone Centered)
  real(fp), dimension(:), allocatable :: xizc_r8
        ! FIXED MESH COORDINATE (Zone Centered) (real(fp))
  real(fp) :: xjav = zero
        ! AVG ZONE INDEX OF LAST HALF ORBIT
  real(fp), dimension(:), allocatable :: xksbsa
        ! VPAR/V OF CH-EX NEUT. SPLITTEES
  real(fp) :: xksflr = zero
        ! BEAM ION VPLL/V AT FLR
  real(fp) :: xksid = zero
        ! VPAR/V OF CURRENT M.C. BEAM PTCL
  real(fp), dimension(:,:), allocatable :: xksidy
        ! STORED VPAR/V OF ALL M.C. BEAM IONS
  real(fp), dimension(:,:), allocatable :: xksidy_saw
        ! STORED VPAR/V OF ALL M.C. BEAM IONS -- post-sawtooth
  real(fp), dimension(:,:), allocatable :: xmuay
        ! mu (vperp**2/B)
  real(fp), dimension(:,:), allocatable :: xmuay_saw
        ! mu (vperp**2/B) -- post-sawtooth
  real(fp) :: xninj = zero
        !TOT NO. OF PTCLS INJECTED N/SEC
  real(fp), dimension(:), allocatable :: xninja
        ! NO. OF PTCLS COMING IN IB'TH BEAMLINE N/SEC
  real(fp) :: xzbeami = zero
        ! CHARGE STATE OF CURRENTLY ORBITTING HEAVY MC ION
  real(fp), dimension(:), allocatable :: xzbeams
  real(fp), dimension(:,:), allocatable :: xzbmay
        ! SAVED CHARGE STATE OF ORBITTING HEAVY MC ION
  real(fp), dimension(:,:), allocatable :: xzeffc
        ! PLASMA COMP. ZEFF AT CURRENT TIMESTEP DIFFERS FROM
  real(fp), dimension(:), allocatable :: xzimpj
        ! Average impurity atomic number by zones
  real(fp), dimension(:,:), allocatable :: xzionay
        ! Z of GC, prev. step
  real(fp), dimension(:,:), allocatable :: xzionay_saw
        ! Z of GC, prev. step -- post-sawtooth
  real(fp) :: yaxis = zero
        ! ELEVATION OF MAGNETIC AXIS (CM)
  real(fp) :: ybgrid = zero
        ! STARTING Y DISPLACEMENT RELATIVE TO CENTERLINE OF BEAM
  real(fp) :: ycur = zero
        ! last inside-the-plasma radius value
  real(fp) :: yflr = zero
        ! HEIGHT FOR CURRENT BEAM ION (NOT G.C.)
  real(fp) :: yion = zero
        ! VERT DISPLACEMENT OF BEAM ION G.C.
  real(fp), dimension(:), allocatable :: ympx
        ! MIDPLANE Y LOCATION, BLOATED GRID
  logical, dimension(:,:), allocatable :: nlreact
        ! REACTION LIST FOR EACH FUSION PRODUCT (T/F)
  logical :: nlgoose_reset = .FALSE.
        ! set T if a goose=1 orbit bounce must be calculated
        ! for every particle at the start of each new beam time step
  real(fp) :: dtrap = zero
        ! trapping depth of ion on last orbit half-bounce
        !    if ximax = max xion of orbit on last half bounce
        !    and bimax = max B at orbit guiding center on last half bounce
        ! dtrap = -1 for passing orbits-- bimax occurs at midplane crossing
        ! for trapped orbits,
        !    dtrap = (Bmax(ximax) - bimax)/(Bmax(ximax)-Bmin(ximax))
  real(fp) :: dpass = zero
        ! passing depth dpass = 1.0 - Bmax(orb)/Brefl
        ! where Brefl = Bmax*(vion/vperp)**2 = est. mod(B) that would
        !   be needed for orbit reflection;
        !   Bmax(orb) << Brefl  => dpass =~1.0, "deeply passing" orbit
  integer :: midspl = 0
        ! size of DATSPL (t.b.d.)
  real(fp), dimension(:), allocatable :: datspl
        ! spline buffer
        !CP RA DATSPL(IDSPL)  spline buffer
  integer, dimension(:), allocatable :: ndatspl
        ! index buffer
        !CP IA NDATSPL(INSPL)
        !CP I  NTHSPL        *no. of points in theta grid -- cf dimensions
  integer :: LSPTHETA = 0
        !CP I  LSPTHETA       theta grid
        !C
        !C  following are addresses of bicubic spline functions:
        !CP I  LSPTHETS       Boozer theta-star(xi,theta)
        !CP I  LSPPHILB       Boozer nu(xi,theta) straight field line coordinate
        !C                    adjustor; related to field line length.
        !CP I  LSPB           mod(B) vs (xi,theta) -- may be unused
        !CP I  LSDELCR        Ripple Threshhold function
        !CP I  LSLOSS         Loss map
        !CP I  LSPUSE         last address in use
  integer :: LSPTHETS = 0
  integer :: LSPPHILB = 0
  integer :: LSPB = 0
  integer :: LSDELCR = 0
  integer :: LSLOSS = 0
  integer :: LSPUSE = 0
  integer :: NPHIPZC = 0
  integer :: NSDSPL = 0
        !C
        !C  for phi(p) representation in spline buffer -- singularities occur
        !C  so a high-ish resolution bilinear fit will be used, with singularity
        !C  data provided.  treatment at singularities not totally critical since
        !C  then ripple well trapping occurs and phi(p)'/phi(b)' model does not
        !C  accurately apply in any case.
        !C
        !C  dmc 28 June 1996, implementing "Ripple-Induced Energetic Particle
        !C    Loss in Tokamaks" by R.B. White et al. (priv. comm. May 1996)
        !C
        !CP I  NSDSPL         # of integers / surface
        !CP I  NPHIPZC        * # of phi(p) evaluation surfaces / TRANSP zone
        !  (the following are now dimension constants in NUBEAM module)
        !CP I  NTHPMIN        * minimum number of theta evaluations (near ctr)
        !CP I  NTHPMAX        * maximum number of theta evaluations (near edge)
  integer :: KDEBUG1 = 0
        ! debug flag -- prompt loss
  integer :: KDEBUG2 = 0
        ! debug flag -- prompt loss
  real(fp), dimension(:), allocatable :: XN_DEP
        ! ptcls... (was zndep)
  real(fp), dimension(:), allocatable :: XN_ORB
  real(fp), dimension(:), allocatable :: XP_DEP
        ! power... (was zpdep)
  real(fp), dimension(:), allocatable :: XP_ORB
  real(fp), dimension(:), allocatable :: XPH_DEP
        ! momentum... (was zphdep)
  real(fp), dimension(:), allocatable :: XPH_ORB
  real(fp) :: rsta = zero
        ! R coordinate -- start of neutral track
  real(fp) :: ysta = zero
        ! Y coordinate -- start of neutral track
  real(fp) :: esta = zero
        ! toroidal (zeta) coordinate -- start of neutral track (deg)
  real(fp) :: thea = zero
        ! theta aiming angle -- start of neutral track
  real(fp) :: phia = zero
        ! phi aiming angle -- start of neutral track
  real(fp) :: rosta = zero
        ! rho coordinate -- start of neutral track
  real(fp) :: xltkbf = zero
        ! distance from start point to plasma
  real(fp) :: xltkmd = zero
        ! length of vacuum segment (exit and re-entry) (if any)
  real(fp) :: xltkwl = zero
        ! approximate length of track to far wall
  integer :: nzcrtk = 0
        ! number of track zone intersections
  integer :: nztkmd = 0
        ! end of 1st segment, if track exits and re-enters plasma
  integer, dimension(:), allocatable :: nznttk
        ! intersected radial zone indices
  integer, dimension(:), allocatable :: nznbtk
        ! intersected 2d MC grid zone indices
  real(fp), dimension(:), allocatable :: xlngtk
        ! intersection lengths
  integer :: id_rhozb = 0
        ! XPLASMA rho id for boundary oriented profiles
  integer :: id_rhozc = 0
        ! XPLASMA rho id for zone oriented profiles
  integer :: num_rhozb = 0
        ! #pts in id_rhozb grid
  integer :: num_rhozc = 0
        ! #pts in id_rhozc grid
  real(fp), dimension(:), allocatable :: ebcons
        ! ENERGY CONTAINED IN BEAM IONS BY SPECIES
  real(fp), dimension(:), allocatable :: enolds
        ! ENERGY IN BEAM IONS BY ION SPECIES, PREV STEP
  real(fp), dimension(:,:), allocatable :: ereact
        ! SOURCE ENERGY FOR EACH FUSION PRODUCT REACTION
  real(fp), dimension(:,:,:), allocatable :: cfa
        ! ENHANCED FOKKER-PLANCK COEFFICIENT ARRAY
  real(fp), dimension(:,:,:), allocatable :: dfa
        ! ENHANCED FOKKER-PLANCK COEFFICIENT ARRAY
  real(fp), dimension(:,:,:), allocatable :: efa
        ! ENHANCED FOKKER-PLANCK COEFFICIENT ARRAY
  real(fp), dimension(:,:,:), allocatable :: exc0bb
        ! excitation effect on bb-i cross section (dep)
  real(fp), dimension(:,:,:), allocatable :: exc0bx
        ! excitation effect on bb-cx cross section (dep)
  real(fp), dimension(:,:,:), allocatable :: exc0cx
        ! excitation effect on th-cx cross section (dep)
  real(fp), dimension(:,:,:), allocatable :: exc0ie
        ! excitation effect on e-i cross section (dep)
  real(fp), dimension(:,:,:), allocatable :: exc0ii
        ! excitation effect on th.-i cross section (dep)
  real(fp), dimension(:,:,:), allocatable :: exc0iz
        ! excitation effect on impurity-i cross section (dep)
  real(fp), dimension(:,:), allocatable :: excrbb
        ! excitation effect on bb-i cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: excrbx
        ! excitation effect on bb-cx cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: excrcx
        ! excitation effect on th-cx cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: excrie
        ! excitation effect on e-i cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: excrii
        ! excitation effect on th.-i cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: excriz
        ! excitation effect on impurity-i cross section (nutrav)
  real(fp), dimension(:,:), allocatable :: wttmf
        ! weighted total magn. factor
  real(fp), dimension(:,:), allocatable :: vmin
        ! m/sec
        ! BEAM ION THERMALIZATION VELOCITIES
  real(fp), dimension(:,:), allocatable :: bv0cx
        ! AVG V(PHI) OF CH-EX BEAM NEUTRALS
  real(fp), dimension(:,:), allocatable :: fpixbts
        ! Newtons
        ! BEAM ION IR X BT POLOIDAL FORCE BY SPECIES
  real(fp), dimension(:,:), allocatable :: fpixbts_ntr
        ! Newtons
        ! BEAM ION IR X BT POLOIDAL FORCE BY SPECIES (part from NUTRAV only)
  real(fp), dimension(:,:), allocatable :: fpixbts_dep
        ! Newtons
        ! BEAM ION IR X BT POLOIDAL FORCE BY SPECIES (Deposition only -- FLR)
  real(fp), dimension(:,:), allocatable :: tirsaws
        ! amps
        ! BEAM ION RADIAL CURRENT: sawtooth, time avg
  real(fp), dimension(:,:), allocatable :: tirtots
        ! amps
        ! BEAM ION RADIAL CURRENT ACROSS SURFACES (ORB JXB)
  real(fp), dimension(:,:), allocatable :: tirtots_ntr
        ! amps
        ! BEAM ION RADIAL CURRENT ACROSS SURFACES (ORB JXB) (part from NUTRAV only)
  real(fp), dimension(:,:), allocatable :: tirtots_dep
        ! amps
        ! BEAM ION RADIAL CURRENT ACROSS SURFACES (Deposition only-- FLR)
  real(fp), dimension(:,:), allocatable :: fbdepco
        ! n/sec (integrated profile)
        ! FPDEP BY ZONE, SPECIES, SUMMED OVER ENERGY AND PITCH:
        ! DEBUG OUTPUT FOR RPLOT: CO DEPOSITION.
  real(fp), dimension(:,:), allocatable :: fbdepct
        ! n/sec (integrated profile)
        ! FPDEP BY ZONE, SPECIES, SUMMED OVER ENERGY AND PITCH:
        ! DEBUG OUTPUT FOR RPLOT: CTR DEPOSITION.
  real(fp), dimension(:), allocatable :: xwght0
        ! WGHTA-related axis to edge weight ratio
  real(fp), dimension(:), allocatable :: wghtsave
        ! WGHT multiplication factor for newly deposited ions
  real(fp), dimension(:), allocatable :: wkill
        ! WGHT taken by Russian Roulette in DEPALL->NBKILL
  real(fp), dimension(:), allocatable :: wnkill
        ! #of test ptcls taken by Russian Roulette, DEPALL->NBKILL
  real(fp), dimension(:), allocatable :: wghtakil
        ! WGHT multiplicaton factor for old orbiting ions
        ! (arrange that ptcl count is preserved by DEP RR)
  real(fp), dimension(:), allocatable :: wnlight
        ! #of underweight ptcls seen in previous step
  real(fp), dimension(:), allocatable :: wnlkill
        ! #of underweight ptcls killed (see also sforbrr(...))
  real(fp), dimension(:), allocatable :: wndoubl
        ! #of underweight ptcls with weight doubled
  real(fp), dimension(:,:), allocatable :: wbav
        ! AVG WEIGHT OF M.C. BEAM IONS AFTER DEPOSITION, used
        ! to control splittings
        ! if wghta = 1.0, then wbav will be constant =  wghtot/nnew
        ! sum of WGHT actually deposited for subsequent orbiting,
        ! divided by the number of MC ions actually deposited.
  real(fp), dimension(:,:), allocatable :: wbav_kill
        ! number of physical particles / Monte Carlo particle
        ! AVG WEIGHT OF M.C. BEAM IONS AFTER DEPOSITION.
        ! used to call Russian Roulette suicide
  real(fp) :: fact_wbav = zero
        !the renormalization factor for time dependence of power included to
        !estimate wbav
  real(fp) :: wghtot = zero
        ! the sum of WGHT actually deposited for subsequent orbiting, divided by the number
        ! of MC ions actually deposited WGHT/NNEW. Needed to estimate wbav.
  real(fp), dimension(:), allocatable :: pljbbm
        ! BEAM DRIVEN PLL CURRENT * B FIELD
  real(fp), dimension(:,:), allocatable :: pljbbms
        ! FAST ION DRIVEN <J.B>S, BY SPECIES
  real(fp), dimension(:,:), allocatable :: pljbbm_nb
        ! FAST ION DRIVEN <J.B>S, BY beam line
  real(fp), dimension(:,:), allocatable :: pljbcct
        ! UNSHIELDED <J.B>, CO sources ONLY
  real(fp), dimension(:), allocatable :: pljbfi
        ! FUSION ION DRIVEN PLL CURRENT * B FIELD
  real(fp), dimension(:,:,:), allocatable :: taubpa
        ! sec
        ! BEAM ION P.A. SCATTERING TIMES - OUTPUT PROFILES
  real(fp), dimension(:,:,:), allocatable :: taubsl
        ! sec
        ! BEAM ION SLOWING DOWN TIMES - OUTPUT PROFILES
  real(fp), dimension(:,:,:), allocatable :: taubpa_tmp
        ! sec
        ! stores taubpa temporarily before parallel loop
  real(fp), dimension(:,:,:), allocatable :: taubsl_tmp
        ! sec
        ! stores taubsl temporarily before parallel loop
        ! taubsl_tmp(mj,1,mibs) tau sl.down for fusion product
  real(fp), dimension(:), allocatable :: tausl_avg
        ! sec
        ! volume average TAU SLOWING DOWN
  real(fp), dimension(:,:), allocatable :: tuben0
        ! neutral density along the tube indexed by (isection, itube) -- (#/cm**3)
  real(fp), dimension(:,:), allocatable :: tubedot
        ! average cos(angle_beam_to_tube) along the tube indexed by (isection, itube)
  real(fp), dimension(:), allocatable :: tubes0
        ! distance from start to beginning of first section
  integer, dimension(:), allocatable :: itube_try
        ! number of beam particle trajectories attempted
  integer :: NBPLOT = 0
        ! flag to signal debug plotting when debug compile
  integer :: NBPLIST = 0
        ! LABEL FLAG FOR BEAM ION ORBIT PLOTTING
  real(fp), dimension(:), allocatable :: XPOS
        ! BEAM ORBIT PLOTTING ARRAY
  real(fp), dimension(:), allocatable :: YPOS
        ! BEAM ORBIT PLOTTING ARRAY
  integer :: ndep_beam_use = 0
        ! if .gt.0, restrict data collection to this beam#
  integer :: ndep_ien_use = 0
        ! if .gt.0, restrict energy fraction to this
  integer, dimension(:), allocatable :: mibs2use
        ! maps lsbeam/isb to mibs_use
  real(fp), dimension(:,:), allocatable :: TRACK_nbsce
        ! beam index
  real(fp), dimension(:,:), allocatable :: TRACK_nbien
        ! energy fraction index
  real(fp), dimension(:,:), allocatable :: XNTRACK
        ! # of track zones (FP for MPI-reducibility)
  real(fp), dimension(:,:), allocatable :: TRACK_RIN
        ! R of 1st entry point of track into plasma
  real(fp), dimension(:,:), allocatable :: TRACK_ZIN
        ! Z of 1st entry point of track into plasma
  real(fp), dimension(:,:), allocatable :: TRACK_PHIIN
        ! Phi of 1st entry point of track into plasma, radian
  real(fp), dimension(:,:), allocatable :: TRACK_ROUT
        ! R of last exit point of track from plasma
  real(fp), dimension(:,:), allocatable :: TRACK_ZOUT
        ! Z of last exit point of track from plasma
  real(fp), dimension(:,:), allocatable :: TRACK_PHIOUT
        ! Phi of last entry point of track into plasma, radian
  real(fp), dimension(:,:), allocatable :: TRACK_TIME
        ! time of track
  real(fp), dimension(:,:), allocatable :: TRACK_EINJ
        ! injection energy (eV) of track
  real(fp), dimension(:,:), allocatable :: TRACK_SINJ
        ! track weight (ptcls)
        ! actual deposition (ptcls) is TRACK_SINJ*[deposition probability]
        ! divide by 1/[time step] to get rate.  See TRACK_PDEP.
  real(fp), dimension(:), allocatable :: TRACK_XL
        ! relative distance to end of each track zone
  real(fp), dimension(:), allocatable :: TRACK_PDEP
        ! deposition probability to end of zone
        ! range 0:1
  real(fp), dimension(:), allocatable :: TRACK_DR_FLR
        ! delta(R) from dep point to GC point
  real(fp), dimension(:), allocatable :: TRACK_DZ_FLR
        ! delta(Z) from dep point to GC point
  real(fp), dimension(:), allocatable :: TRACK_DE_FLR
        ! delta(E) from dep point to GC point
  real(fp), dimension(:), allocatable :: TRACK_VPV_DEP
        ! vpll/v at deposition point (not GC).
  real(fp), dimension(:), allocatable :: TRACK_DVPV_FLR
        ! delta(vpll/v), dep point to GC point
  real(fp), dimension(:), allocatable :: TRACK_XJADEP
        ! 1st orbit characterization
        ! (xjadep estimate): 0=passing, 1=trapped, 2=passing & lost, 3=trapped & lost
  real(fp) :: RMJDEP = zero
        ! major radius at deposition point
  real(fp) :: YDEP = zero
        ! elevation at deposition point
  real(fp) :: XKSID_DEP = zero
        ! vpll/v at deposition point
  real(fp) :: VION_DEP = zero
        ! v at deposition point
  real(fp) :: X_HOR = zero
        !(ZXHOR,ZYHOR) position projected onto midplane, relative to start point
  real(fp) :: Y_HOR = zero
        !(RSTA,0); ZR = sqrt(ZXHOR**2+ZYHOR**2)
  character*10 :: oprefix = ' '
        ! prefix for all xplasma output profile names
  character*120 :: chbuf_debug_mpi = ' '
        ! string buffer for debug_mpi messages
  real(fp), dimension(:), allocatable :: gtime_lost
        ! time when ptcl. has been lost
  real(fp), dimension(:), allocatable :: ynbsce_lost
        ! beam# from what lost ptcl.came, 0 - means fusion products
  real(fp), dimension(:), allocatable :: ynbien_lost
        ! energy fraction # of lost ptcl.
  real(fp), dimension(:), allocatable :: wght_lost
        ! weight of lost ptcl. (#)
  real(fp), dimension(:), allocatable :: rmjion_lost
        ! R at what ptcl. has been lost
  real(fp), dimension(:), allocatable :: yyion_lost
        ! Z at what ptcl. has been lost
  real(fp), dimension(:), allocatable :: yze_lost
        ! energy of lost ptcl.
  real(fp), dimension(:), allocatable :: xksid_lost
        ! pitch angle of lost ptcl.
  real(fp), dimension(:), allocatable :: goose_lost
        ! goose factor of lost ptcl., if less than 1.00001_fp - prompt loss
  real(fp), dimension(:), allocatable :: species_lost
        ! species of lost ptcl.
  !==============================================
  ! OUTPUT items...
 
  ! ... block: fusion .........
    ! mireact -- reaction index dimension.
    ! xplasma ids:  io_bbntns(mireact), etc.
    !
    ! a legacy system for categorizing fusion reactions is in use.
    ! Reactions indexing is from 1 to NREACT, get NREACT via
    !      call catfus_num(NREACT)
    ! as of Jan 2002 NREACT=14.
    !
    ! Originally, these arrays held "beam-target" and "beam-beam"
    ! reactions; now these arrays also hold "FP-target","FP-beam",
    ! "RF-target", and "RF-beam" reactions, FP = fusion product,
    ! RF = RF tail ions, with separate indices allocated.
    !
    !      integer ireact,itype1
    !      character*3 reagent1,reagent2
    !      do ireact = 1,NREACT
    !         call catfus_mtype(ireact,itype1)
    !         call catfus_reagent(reagent1,reagent2)
    !         write(nonlin,*) itype1,' ',reagent1,reagent2
    !      enddo
    ! will return in sequence the first component type for each
    ! reaction index ireact: itype1=1 for beams, itype1=2 for FP,
    ! itype1=3 for RF; and the reagent codes "D  ","T  ", or "HE3".
    ! For the DD(n) reaction reagent2 = "D:N"; for the DD(p) reaction
    ! reagent2 = "D:P" are returned
    !
    ! The "catfus" calls are included with the beam module in the
    ! library "comput.a".  A working sample code and output follow.
    !
    !      program foo
    !      implicit NONE
    !
    !      integer nreact
    !      integer ireact,itype1
    !      character*3 reagent1,reagent2
    !
    !      call catfus_num(nreact)
    !      do ireact = 1,NREACT
    !         call catfus_mtype(ireact,itype1)
    !         call catfus_reagent(ireact,reagent1,reagent2)
    !         write(nonlin,*) '! ireact=',ireact,' itype1=',itype1,
    !     >      ' reagent1,2= ',reagent1,' ',reagent2
    !      enddo
    !
    !      stop
    !      end
    !
    ! ireact= 1  itype1= 1  reagent1,2= D   T
    ! ireact= 2  itype1= 1  reagent1,2= D   HE3
    ! ireact= 3  itype1= 1  reagent1,2= D   D:P
    ! ireact= 4  itype1= 1  reagent1,2= D   D:N
    ! ireact= 5  itype1= 1  reagent1,2= T   T
    ! ireact= 6  itype1= 1  reagent1,2= T   HE3
    ! ireact= 7  itype1= 1  reagent1,2= T   D
    ! ireact= 8  itype1= 1  reagent1,2= HE3 D
    ! ireact= 9  itype1= 2  reagent1,2= T   D
    ! ireact= 10  itype1= 2  reagent1,2= HE3 D
    ! ireact= 11  itype1= 3  reagent1,2= T   D
    ! ireact= 12  itype1= 3  reagent1,2= HE3 D
    ! ireact= 13  itype1= 3  reagent1,2= D   T
    ! ireact= 14  itype1= 3  reagent1,2= D   HE3
    !
    ! Here again, itype1 characterizes the "first" reagent:
    ! itype1=1 means "beam", itype1=2 means "fusion product",
    ! itype1=3 means "RF tail".  The second reagent is either
    ! a beam ion (output to "bb" arrays) or a thermal target
    ! ion (outupt to "bt" arrays).  Fast ion - Fast ion reactions
    ! where the 2nd reagent is a fusion product or RF tail ion
    ! are never computed in the current software setup.
    !
    ! caveats: As of Jan 2002, no cross sections were available for
    ! the T-HE3 reaction, so this index is a place-holder, and the
    ! reaction rates returned are always zero.  Also, regarding the
    ! RF reactions:  (a) "D RF minority" experiments (ireact=13,14)
    ! have probably never been done; (b) the Monte Carlo code does
    ! not have an RF operator and so cannot follow RF tail ions.
    ! Instead, the RF tail ion distribution function can be computed
    ! by another code (e.g. a Fokker-Planck package) and the
    ! distribution function passed into the MC package, and an MC
    ! integral computed for the RF-beam fast ion - fast ion fusion
    ! rate, output to the "bb" arrays.  The RF-(thermal) target
    ! fusion calculation is left to the Fokker-Planck package or
    ! whatever code is computing the RF tail distribution.  Thus
    ! for the "bt" arrays only the rows corresponding to ireact=1:10
    ! will ever be set.  This is the situation as of Jan 2002; the
    ! code may be updated in the future according to user requirements.
    !
    ! The cross sections for these reaction rate integrations come
    ! from the NTCC module "preact", and, originally, from:
    !
    ! H-S. Bosch, Review of Data and Formulas for Fusion Cross-sections,
    ! IPP I/252 (Sept 1990).
    !
    ! the output profiles in this block are flux surface averaged.
    ! 2d profiles showing poloidal variation are also available.
    !------------------------
    ! for "mireact" indexed quantities, the labeling is handled
    ! specially in nubeam/nbo_name1 ...
    !------------------------
  real(fp), dimension(:,:), allocatable :: bbntns
  INTEGER, dimension(:), allocatable :: io_bbntns
  CHARACTER*32, dimension(:), allocatable :: io_bbntns_name
        ! '#/S; fast ion -- fast ion fusion'
        ! BEAM-BEAM NEUTRON EMISSION 1D PROFILE, by reaction
        ! ALSO: INFO ON FUSION REACTIONS WHICH DO NOT MAKE NEUTRONS
  real(fp), dimension(:), allocatable :: bbntnt
  INTEGER :: io_bbntnt = 0
  CHARACTER*32 :: io_bbntnt_name = ' '
        ! '#/S; beam-beam neutrons'
        ! BEAM-BEAM TOTAL NEUTRON EMISSION 1D PROFILE
  real(fp), dimension(:,:), allocatable :: btntns
  INTEGER, dimension(:), allocatable :: io_btntns
  CHARACTER*32, dimension(:), allocatable :: io_btntns_name
        ! '#/S; fast ion -- target fusion'
        ! BEAM-TARGET NEUTRON EMISSION PROFILE, by reaction
        ! ALSO: INFO ON FUSION REACTIONS WHICH DO NOT MAKE NEUTRONS
  real(fp), dimension(:,:), allocatable :: btntrs
  INTEGER, dimension(:), allocatable :: io_btntrs
  CHARACTER*32, dimension(:), allocatable :: io_btntrs_name
        ! 'cm^3/S; fast ion -- target fusion, sigv'
        ! BEAM-TARGET NEUTRON sigv PROFILE, by reaction
        ! ALSO: INFO ON FUSION REACTIONS WHICH DO NOT MAKE NEUTRONS
  real(fp), dimension(:,:), allocatable :: btntws
  INTEGER, dimension(:), allocatable :: io_btntws
  CHARACTER*32, dimension(:), allocatable :: io_btntws_name
        ! '#/cm^3; fast ions react with thermal'
  real(fp), dimension(:), allocatable :: btntnt
  INTEGER :: io_btntnt = 0
  CHARACTER*32 :: io_btntnt_name = ' '
        ! '#/S; beam-target neutrons'
        ! BEAM-TARGET TOTAL NEUTRON EMISSION PROFILE
  real(fp), dimension(:), allocatable :: rbntnt
  INTEGER :: io_rbntnt = 0
  CHARACTER*32 :: io_rbntnt_name = ' '
        ! '#/S; beam-RF minority neutrons'
        ! BEAM-RF minority TOTAL NEUTRON EMISSION PROFILE
  real(fp) :: bbntot = zero
        ! '#/S; total beam-beam neutrons'
        ! TOTAL BEAM-BEAM NEUTRON PRODUCTION
  real(fp), dimension(:), allocatable :: bbnta
        ! '#/S; fast ion -- fast ion fusion'
        ! BEAM-BEAM NEUTRONS, BY FUSION REACTION
        ! ALSO: INFO ON FUSION REACTIONS WHICH DO NOT MAKE NEUTRONS
  real(fp) :: btneut = zero
        ! '#/S; total beam-target neutrons'
        ! BEAM TARGET NEUTRONS PER SECOND
  real(fp) :: rbneut = zero
        ! '#/S; total beam-RF minority neutrons'
        ! BEAM RF minority NEUTRONS PER SECOND
  real(fp), dimension(:), allocatable :: btnta
        ! '#/S; fast ion -- target fusion'
        ! BEAM TARGET NEUTRONS, BY SOURCE REACTION
        ! ALSO: INFO ON FUSION REACTIONS WHICH DO NOT MAKE NEUTRONS
  real(fp), dimension(:), allocatable :: sftot
        ! '#/S; total fusion product source'
        ! TOTAL FUSION PRODUCT SOURCE, BY REACTION INDEX
        ! includes beam-beam, beam-target, and thermonuclear contributions,
        ! as appropriate to reaction index.
  ! ... block: density .........
    ! mibs -- fast ion specie index
    ! xplasma ids:  io_bdenss(mibs)
    !
    ! post-sawtooth arrays: if sawtooth calculation is requested,
    ! these are the specie densities expected after a Kadomtsev
    ! sawtooth remix of affected fast ion species.
  real(fp), dimension(:,:), allocatable :: bdenss
  INTEGER, dimension(:), allocatable :: io_bdenss
  CHARACTER*32, dimension(:), allocatable :: io_bdenss_name
        ! '#; fast ion density'
        ! total number of fast ions, by specie.  Orbit averaged.
        ! Fusion product ion densities are also here.
  real(fp), dimension(:,:), allocatable :: mc_dens
  INTEGER, dimension(:), allocatable :: io_mc_dens
  CHARACTER*32, dimension(:), allocatable :: io_mc_dens_name
        ! '#/m^3; MC fast ion density (unsmoothed)'
        ! instantaneous density from particle list at end of time
        ! step; no orbit averaging; no smoothing.
  real(fp), dimension(:,:), allocatable :: bdenssw
  INTEGER, dimension(:), allocatable :: io_bdenssw
  CHARACTER*32, dimension(:), allocatable :: io_bdenssw_name
        ! '#; post-sawtooth fast ion density'
        ! same as bdenss, but if a sawtooth occur this is what
        ! the density profiles change to.
  real(fp), dimension(:,:), allocatable :: udenspl
  INTEGER, dimension(:), allocatable :: io_udenspl
  CHARACTER*32, dimension(:), allocatable :: io_udenspl_name
        ! 'J; pll energy density'
        ! pll energy density profiles (Joules, integrated)
  real(fp), dimension(:,:), allocatable :: udenspp
  INTEGER, dimension(:), allocatable :: io_udenspp
  CHARACTER*32, dimension(:), allocatable :: io_udenspp_name
        ! 'J; perp energy density'
        ! perp energy density profiles (Joules, integrated)
  real(fp), dimension(:,:), allocatable :: udensplw
  INTEGER, dimension(:), allocatable :: io_udensplw
  CHARACTER*32, dimension(:), allocatable :: io_udensplw_name
        ! 'J; post-sawtooth pll energy density'
        ! pll energy density profiles (post-sawtooth)
  real(fp), dimension(:,:), allocatable :: udensppw
  INTEGER, dimension(:), allocatable :: io_udensppw
  CHARACTER*32, dimension(:), allocatable :: io_udensppw_name
        ! 'J; post-sawtooth perp energy density'
        ! perp energy density profiles (post-sawtooth)
  ! ... block: deposition .........
    ! mibs -- fast ion specie index
    ! 1:3 -- full/half/third energy index
    ! xplasma ids:  io_bdepe(3,mibs), io_bdeps(mibs),etc.
    ! mig -- thermal gas index
    !
    ! total deposition profiles f, f(x) = n/sec inside surface "x".
    ! components broken down by stopping reaction type, generally
    ! given as a sum over all beams and energy fractions, as a sum
    ! over energy fraction for each beam specie, and as a separate
    ! profile for each beam specie and energy fraction.  For non-
    ! Hydrogenic beams, the 2nd and 3rd energy fractions are always
    ! zero.
  real(fp), dimension(:,:,:), allocatable :: bdepb
  INTEGER, dimension(:,:), allocatable :: io_bdepb
  CHARACTER*32, dimension(:,:), allocatable :: io_bdepb_name
        ! '#/S; deposition rate'
        ! deposition profile at each energy fraction due to each
        ! neutral beam injector.
  real(fp), dimension(:,:), allocatable :: bdepb_nb
  INTEGER, dimension(:), allocatable :: io_bdepb_nb
  CHARACTER*32, dimension(:), allocatable :: io_bdepb_nb_name
        ! 'beam deposition profile'
        ! deposition profile to each neutral beam injector.
  real(fp), dimension(:,:,:), allocatable :: bdepe
  INTEGER, dimension(:,:), allocatable :: io_bdepe
  CHARACTER*32, dimension(:,:), allocatable :: io_bdepe_name
        ! '#/S; deposition rate'
        ! total beam deposition, by beam energy fraction and specie.
        ! for non-Hydrogenic beams, the beams are assumed to be 100%
        ! full energy.
  real(fp), dimension(:,:), allocatable :: bdeps
  INTEGER, dimension(:), allocatable :: io_bdeps
  CHARACTER*32, dimension(:), allocatable :: io_bdeps_name
        ! '#/S; deposition rate'
        ! total beam deposition, summed over energy fractions
        ! profile is sum of results of integrations along beam neutral tracks.
  real(fp), dimension(:,:), allocatable :: sdep_oa
  INTEGER, dimension(:), allocatable :: io_sdep_oa
  CHARACTER*32, dimension(:), allocatable :: io_sdep_oa_name
        ! '#/S; deposition rate'
        ! the orbit averaged total deposition fast ion source
        ! NEW PARTICLES BEAM DEPOSITION PROFILES, BY SPECIE
        ! produce a quantity closely related to BDEPS
        ! but, vs. the orbit averaged position <X> (JIAV)
        ! instead of the instantaneous deposition location X (J) "BDEPS"
        ! FAST ION SCE --orbit averaged total BEAM DEPOSITION
  real(fp), dimension(:), allocatable :: SDEP_OAUC
        ! '#/S; deposition rate'
        ! the orbit averaged deposition fast ion source, unconfined
        ! same as sdep_oa(LEP1,:)
  real(fp), dimension(:,:), allocatable :: mc_newdens
  INTEGER, dimension(:), allocatable :: io_mc_newdens
  CHARACTER*32, dimension(:), allocatable :: io_mc_newdens_name
        ! '#/m^3; MC fast ions added (unsmoothed)'
        ! Monte Carlo fast ion source -- as actually deposited to MC particle
        ! list on most recent time step; no smoothing.
  real(fp), dimension(:), allocatable :: einjav
        ! 'KeV; nominal injection energy'
        ! for beams this is based on the full energy component only.
        ! if there is no power, the value for the species is retained
        ! from the last time that there was non-zero power.
  real(fp), dimension(:,:,:), allocatable :: sbbidpe
  INTEGER, dimension(:,:), allocatable :: io_sbbidpe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbbidpe_name
        ! '#/S; dep. by beam-beam ionization'
        ! BEAM *DEPOSITION* BY BEAM-BEAM II, BY E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbbidps
  INTEGER, dimension(:), allocatable :: io_sbbidps
  CHARACTER*32, dimension(:), allocatable :: io_sbbidps_name
        ! '#/S; dep. by beam-beam ionization'
        ! BEAM *DEPOSITION* BY BEAM-BEAM II, BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:), allocatable :: sbbxdpe
  INTEGER, dimension(:,:), allocatable :: io_sbbxdpe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbbxdpe_name
        ! '#/S; dep. by beam-beam charge exchange'
        ! BEAM *DEPO* BY BEAM-BEAM CX, BY E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbbxdps
  INTEGER, dimension(:), allocatable :: io_sbbxdps
  CHARACTER*32, dimension(:), allocatable :: io_sbbxdps_name
        ! '#/S; dep. by beam-beam charge exchange'
        ! BEAM *DEPOSITION* BY BEAM-BEAM CX, BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:,:), allocatable :: sbcxdpe
  INTEGER, dimension(:,:,:), allocatable :: io_sbcxdpe
  CHARACTER*32, dimension(:,:,:), allocatable :: io_sbcxdpe_name
        ! '#/S; dep. by charge exchange'
        ! BEAM *DEPO* BY THERM. CX, BY E.FRAC. & by species
  real(fp), dimension(:,:,:), allocatable :: sbcxdps
  INTEGER, dimension(:,:), allocatable :: io_sbcxdps
  CHARACTER*32, dimension(:,:), allocatable :: io_sbcxdps_name
        ! '#/S; dep. by charge exchange'
        ! BEAM *DEPOSITION* BY THERMAL CX, BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:), allocatable :: sbedepe
  INTEGER, dimension(:,:), allocatable :: io_sbedepe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbedepe_name
        ! '#/S; dep. by ion impact ionization'
        ! BEAM *DEPO* BY plasma ion II, by E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbedeps
  INTEGER, dimension(:), allocatable :: io_sbedeps
  CHARACTER*32, dimension(:), allocatable :: io_sbedeps_name
        ! '#/S; dep. by ion impact ionization'
        ! BEAM *DEPOSITION* BY plasma ion IONIZ., BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:), allocatable :: sbiedepe
  INTEGER, dimension(:,:), allocatable :: io_sbiedepe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbiedepe_name
        ! '#/S; dep. by e- impact ionization'
        ! BEAM *DEPO* BY electron II, by E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbiedeps
  INTEGER, dimension(:), allocatable :: io_sbiedeps
  CHARACTER*32, dimension(:), allocatable :: io_sbiedeps_name
        ! '#/S; dep. by e- impact ionization'
        ! BEAM *DEPOSITION* BY ELEC. IMPACT IONIZ., BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:), allocatable :: sbildepe
  INTEGER, dimension(:,:), allocatable :: io_sbildepe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbildepe_name
        ! '#/S; dep. by Lorentz ionization'
        ! BEAM *DEPO* BY Lorentz ioniz., by E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbildeps
  INTEGER, dimension(:), allocatable :: io_sbildeps
  CHARACTER*32, dimension(:), allocatable :: io_sbildeps_name
        ! '#/S; dep. by Lorentz ionization'
        ! BEAM *DEPOSITION* BY Lorentz IONIZ., BY SPECIES
        ! summed over energy fractions
  real(fp), dimension(:,:,:), allocatable :: sbizdepe
  INTEGER, dimension(:,:), allocatable :: io_sbizdepe
  CHARACTER*32, dimension(:,:), allocatable :: io_sbizdepe_name
        ! '#/S; dep. by impurity stopping'
        ! BEAM *DEPO* BY IMPURITY II, BY E.FRAC. & by species
  real(fp), dimension(:,:), allocatable :: sbizdeps
  INTEGER, dimension(:), allocatable :: io_sbizdeps
  CHARACTER*32, dimension(:), allocatable :: io_sbizdeps_name
        ! '#/S; dep. by impurity stopping'
        ! BEAM *DEPOSITION* BY IMPURITY IMPACT IONIZ., BY SPECIES
        ! summed over energy fractions
  real(fp) :: frac_cur_ndep0 = zero
        ! '-; track data gathered, this fraction of beam current'
        ! current fraction of beam(s) eligible for
        ! deposition track data gathering
  integer :: ndep0_use = 0
        ! '-; NDEP0 used (for deposition statistics)'
        ! also for deposition track data gathering
        ! if nltrk_dep0 is .TRUE.
  integer :: nbs_step = 0
        ! '-; # of nbsteps used for averaged beam ion birth statistics'
        ! # of steps to gather data from ~ OUTTIM-AVGTIM to OUTTIM
        ! list of ions actually deposited for subsequent
        ! following in orball(...).  This not depend on NDEP0
        ! but rather on our census control based on NPTCLS
        ! (it would vary from step to step based on beam power
        ! history and our census control logic).
  integer :: n_outtim = 0
        ! '-; # for a  current outtim time'
  ! ... block: nclass .........
    ! beam line & beam energy fraction resolved outputs needed by
    ! NCLASS evaluation of neoclassical effects of beam forces
    !
    ! profiles:
    !   bnb_ien(mj,1:3,1:mib) -- density
    !   pbe_ien(mj,1:3,1:mib) -- electron heating
    !   pbi_ien(mj,1:3,1:mib) -- ion heating including impurity
    !   upb_ien(mj,1:3,1:mib) -- (lab frame) <v.B>
    !   bdens_nb(mj, mib) -- beam density
    !   bnbs_ien(mj, 3, mibs) "#; no. of ions"
    !   pbes_ien(mj, 3, mibs) "W; electron heating"
    !       !   pbis_ien(mj, 3, mibs) "W; ion heating"
    ! scalars:
    !   ein_ien(1:3,1:mib) -- average lab frame injection energy
    !   dei_ien(1:3,1:mib) -- RMS variance from ein(...)
    !     (from time variation of injection voltage on beamlines).
    !    bpshna (1:mib,1:3) -- SHINE-THRU POWER, BEAMLINE BY BEAMLINE
    !
    ! by species & fraction instead of by beamline & energy fraction:
    !   (same set of data as above)
    !   (fusion product contributions for NCLASS captured here).
    !
    ! *** warning *** dmc Dec 2003 -- definition and
    ! construction of nclass=1 outputs is not yet complete ***
  real(fp), dimension(:,:,:), allocatable :: bnb_ien
  INTEGER, dimension(:,:), allocatable :: io_bnb_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_bnb_ien_name
        ! '#; no. of ions'
  real(fp), dimension(:,:), allocatable :: bdens_nb
  INTEGER, dimension(:), allocatable :: io_bdens_nb
  CHARACTER*32, dimension(:), allocatable :: io_bdens_nb_name
        ! 'beam density'
        ! to each neutral beam injector
  real(fp), dimension(:,:,:), allocatable :: pbe_ien
  INTEGER, dimension(:,:), allocatable :: io_pbe_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_pbe_ien_name
        ! 'W; electron heating'
  real(fp), dimension(:,:,:), allocatable :: pbi_ien
  INTEGER, dimension(:,:), allocatable :: io_pbi_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_pbi_ien_name
        ! 'W; ion heating'
  real(fp), dimension(:,:,:), allocatable :: upb_ien
  INTEGER, dimension(:,:), allocatable :: io_upb_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_upb_ien_name
        ! 'T*(m/sec); <v.B> lab frame'
  real(fp), dimension(:,:), allocatable :: bpshna
        ! 'SHINE-THRU POWER BY BEAMLINE and energy fraction'
  real(fp), dimension(:,:), allocatable :: ein_ien
        ! 'KeV; avg injection energy'
  real(fp), dimension(:,:), allocatable :: dei_ien
        ! 'KeV; avg inj. energy RMS variation'
        ! the following are computed for beam species only
        ! (fusion product numbers calculated already, elsewhere).
  real(fp), dimension(:,:,:), allocatable :: bnbs_ien
  INTEGER, dimension(:,:), allocatable :: io_bnbs_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_bnbs_ien_name
        ! '#; no. of ions'
  real(fp), dimension(:,:,:), allocatable :: pbes_ien
  INTEGER, dimension(:,:), allocatable :: io_pbes_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_pbes_ien_name
        ! 'W; electron heating'
  real(fp), dimension(:,:,:), allocatable :: pbis_ien
  INTEGER, dimension(:,:), allocatable :: io_pbis_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_pbis_ien_name
        ! 'W; ion heating'
        ! the following are computed for all species, but for
        ! fusion products the energy fraction index is always 1.
  real(fp), dimension(:,:,:), allocatable :: upbs_ien
  INTEGER, dimension(:,:), allocatable :: io_upbs_ien
  CHARACTER*32, dimension(:,:), allocatable :: io_upbs_ien_name
        ! 'T*(m/sec); <v.B> lab frame'
  ! ... block: n0_fast .........
    ! flux surface averaged fast neutral density profiles
    ! mibs -- specie index    1:3 energy fraction index
    ! mib -- beamline index
    !
    ! xplasma ids:  io_bn00ts(3,mibs), etc...
    !
    ! indexing generally similar to indexing for the deposition
    ! profiles; in addition separate profiles for each beam line
    ! are available.
  real(fp), dimension(:,:), allocatable :: bn00t
  INTEGER, dimension(:), allocatable :: io_bn00t
  CHARACTER*32, dimension(:), allocatable :: io_bn00t_name
        ! '#; beam neutral density'
        ! TOTAL 1ST GEN BEAM NEUTRAL DENSITIES
        ! full/half/third energy, summed over all species
  real(fp), dimension(:,:,:), allocatable :: bn00ts
  INTEGER, dimension(:,:), allocatable :: io_bn00ts
  CHARACTER*32, dimension(:,:), allocatable :: io_bn00ts_name
        ! '#; beam neutral density'
        ! TOTAL 1ST GEN BEAM NEUTRAL DENSITIES by species
        ! full/half/third energy
  real(fp), dimension(:,:,:), allocatable :: bn00
  INTEGER, dimension(:,:), allocatable :: io_bn00
  CHARACTER*32, dimension(:,:), allocatable :: io_bn00_name
        ! '#; beam neutral density'
        ! 1D 1ST GEN BEAM NEUTRAL DENSITIES
        ! full/half/third energy, by beamline
  real(fp), dimension(:,:), allocatable :: bn0tot
  INTEGER, dimension(:), allocatable :: io_bn0tot
  CHARACTER*32, dimension(:), allocatable :: io_bn0tot_name
        ! '#; beam neutral density'
        ! TOTAL 1ST GEN BEAM NEUTRAL DENSITIES BY SPECIES
  real(fp), dimension(:,:), allocatable :: bn0cx
  INTEGER, dimension(:), allocatable :: io_bn0cx
  CHARACTER*32, dimension(:), allocatable :: io_bn0cx_name
        ! '#; charge exchange fast neutral density'
        ! BEAM CX-EX fast NEUTRAL DENSITY
  ! ... block: trap_fraction .........
    ! mibs -- fast ion specie index
    ! fraction of "banana trapped" ions
  real(fp), dimension(:,:), allocatable :: btrap
  INTEGER, dimension(:), allocatable :: io_btrap
  CHARACTER*32, dimension(:), allocatable :: io_btrap_name
        ! '-; trapping fraction'
        ! Trapping fraction profile of slowing down population
  real(fp), dimension(:,:), allocatable :: btrap0
  INTEGER, dimension(:), allocatable :: io_btrap0
  CHARACTER*32, dimension(:), allocatable :: io_btrap0_name
        ! '-; trapping fraction'
        ! Trapping fraction profile: newly deposited full energy ions.
  real(fp), dimension(:), allocatable :: tbtrap
        ! '-; total trapping fraction'
        ! Summed trapping fraction (over whole profile)
  real(fp), dimension(:), allocatable :: tbtrap0
        ! '-; total trapping fraction'
        ! Summed trapping fraction (over whole profile) (new full-E ions)
  ! ... block: erngfi_output .........
    ! miefi -- energy range index
    ! mibs -- fast ion specie index
    ! densities and trapping fractions constrained to
    ! energy ranges erngfi(1:nerngfi) as specified by user input.
  real(fp), dimension(:,:,:), allocatable :: btrape
  INTEGER, dimension(:,:), allocatable :: io_btrape
  CHARACTER*32, dimension(:,:), allocatable :: io_btrape_name
        ! '-; trapping fraction'
        ! Trapping fraction profiles by energy range
  real(fp), dimension(:,:,:), allocatable :: densfi
  INTEGER, dimension(:,:), allocatable :: io_densfi
  CHARACTER*32, dimension(:,:), allocatable :: io_densfi_name
        ! '#; fast ion density'
        ! ENERGY-RESTRICTED FAST ION DENSITY PROFILES
  ! ... block: average_energies .........
    ! mibs -- fast ion specie index
    ! average parallel energy, average perp energy, average energy
    !
    ! post-sawtooth arrays: if sawtooth calculation is requested,
    ! these are the average energies expected after a Kadomtsev
    ! sawtooth remix of affected fast ion species.
  real(fp), dimension(:,:), allocatable :: ebapl2s
  INTEGER, dimension(:), allocatable :: io_ebapl2s
  CHARACTER*32, dimension(:), allocatable :: io_ebapl2s_name
        ! 'Kev; avg ion pll energy'
        ! AVG BEAM ION PLL ENERGY @CURRENT STEP BY ION SPECIES
        ! includes fusion product species.
  real(fp), dimension(:,:), allocatable :: ebapp2s
  INTEGER, dimension(:), allocatable :: io_ebapp2s
  CHARACTER*32, dimension(:), allocatable :: io_ebapp2s_name
        ! 'Kev; avg ion perpenergy'
        ! AVG BEAM ION PERP ENERGY @CURRENT STEP BY ION SPECIES
        ! includes fusion product species.
  real(fp), dimension(:,:), allocatable :: ebapl2sw
  INTEGER, dimension(:), allocatable :: io_ebapl2sw
  CHARACTER*32, dimension(:), allocatable :: io_ebapl2sw_name
        ! 'Kev; avg pll ion energy after sawtooth'
        ! AVG BEAM ION PLL ENERGY @CURRENT STEP BY ION SPECIES
        ! includes fusion product species.
  real(fp), dimension(:,:), allocatable :: ebapp2sw
  INTEGER, dimension(:), allocatable :: io_ebapp2sw
  CHARACTER*32, dimension(:), allocatable :: io_ebapp2sw_name
        ! 'Kev; avg perp ion energy after sawtooth'
        ! AVG BEAM ION PERP ENERGY @CURRENT STEP BY ION SPECIES
        ! includes fusion product species.
  real(fp), dimension(:,:), allocatable :: ebavps
  INTEGER, dimension(:), allocatable :: io_ebavps
  CHARACTER*32, dimension(:), allocatable :: io_ebavps_name
        ! 'Kev; avg ion energy'
        ! AVG BEAM ION ENERGY AT ZONE J, BY ION SPECIES
        ! includes fusion product species.
  ! ... block: excited_states .........
    ! mibs -- fast ion specie index
    ! 1:3 -- full/half/third energy index
    ! xplasma ids:  io_excit0(3,mibs), io_excitr(mibs).
    !
    ! excited states correction profiles (flux surface averaged).
    ! =1.0 if no excited states correction included in the
    ! calculation, as per user input.
  real(fp), dimension(:,:,:), allocatable :: excit0
  INTEGER, dimension(:,:), allocatable :: io_excit0
  CHARACTER*32, dimension(:,:), allocatable :: io_excit0_name
        ! '-; deposition adj. (state excitation)'
        ! net effect of excitation on stopping cross section
        ! fast neutrals due to beam injection
        ! (relative to ground state model)
  real(fp), dimension(:,:), allocatable :: excitr
  INTEGER, dimension(:), allocatable :: io_excitr
  CHARACTER*32, dimension(:), allocatable :: io_excitr_name
        ! '-; deposition adjustment (state excitation)'
        ! net effect of excitation on stopping cross section:
        ! fast neutrals due to fast ion charge exchange
        ! (relative to ground state model)
  ! ... block: mc_statistics .........
    ! mibs -- fast ion specie index
    !
    !   MC statistics vs. fast ion radial location (orbit averaged).
    ! these statistics are computed after deposition, before
    ! slowing down.
  real(fp), dimension(:,:), allocatable :: nmc
  INTEGER, dimension(:), allocatable :: io_nmc
  CHARACTER*32, dimension(:), allocatable :: io_nmc_name
        ! '#; number of Monte Carlo particles'
        ! integrated number of Monte Carlo model particles after deposition
        ! but before orbiting.
  real(fp), dimension(:,:), allocatable :: wnmc
  INTEGER, dimension(:), allocatable :: io_wnmc
  CHARACTER*32, dimension(:), allocatable :: io_wnmc_name
        ! '#; number of ions in MC particles'
        ! physical number of ptcls (calculated when nmc is calculated)
        ! this is different than the fast ion density in that the latter
        ! is averaged over the entire timestep; calculated after deposition
        ! but before orbiting.
  integer, dimension(:), allocatable :: nmctot
        ! '#; number of Monte Carlo particles'
        ! total number of MC ptcls by fast species; calculated after deposition
        ! but before orbiting.
  real(fp), dimension(:), allocatable :: wnmctot
        ! '#; number of ions in MC particles'
        ! total physicsl number of ptcls; calculated after deposition
        ! but before orbiting.
  integer, dimension(:), allocatable :: naeptot
        ! '#; number of Monte Carlo particles'
        ! total number of MC ptcls by fast species evolved by MHD transport;
  integer, dimension(:,:), allocatable :: nstep_aep
        ! '#; number time step for inner MonteCarlo loop in MHD transport (M.Podesta)'
        ! total number of  time steps used by MHD transport;
  integer, dimension(:), allocatable :: nmcloss
        ! '#; number of prompt loss Monte Carlo particles'
        ! total number of MC ptcls on prompt loss orbits
  integer, dimension(:), allocatable :: nnew
        ! '# OF NEW M.C. PTCLS FOR DEPOSITION'
        ! # OF NEW M.C. PTCLS FOR DEPOSITION
  real(fp), dimension(:), allocatable :: xncx0_orb
        ! '# OF CX events during orbiting'
        ! the number of CX tracks that we get from orbiting
  ! ... block: rotation .........
    ! mibs -- fast ion specie index
    ! fast ion flux surface averaged angular velocities
    ! (flux surface averaged toroidal rotation profiles)
    ! xplasma ids:  io_omegabs, io_omegabsw
  real(fp), dimension(:,:), allocatable :: omegabs
  INTEGER, dimension(:), allocatable :: io_omegabs
  CHARACTER*32, dimension(:), allocatable :: io_omegabs_name
        ! 'rad/S; avg fast ion toroidal angular velocity'
        ! APPROX TOROIDAL ANGULAR VELOCITY OF IONS, BY FAST ION SPECIES
  real(fp), dimension(:,:), allocatable :: omegabsw
  INTEGER, dimension(:), allocatable :: io_omegabsw
  CHARACTER*32, dimension(:), allocatable :: io_omegabsw_name
        ! 'rad/S; avg fast ion ang. veloc. after sawtooth'
        ! APPROX TOROIDAL ANGULAR VELOCITY OF IONS, BY FAST ION SPECIES
        ! post-sawtooth profile
        ! note: dmc 15 March 2002: not computed; for now omegabsw=omegabs
  ! ... block: compression .........
    ! mibs -- fast ion specie index
    ! energy added (or removed) from fast population by adiabatic
    ! compression (or decompression).
    ! These are informational profiles only.
  real(fp), dimension(:), allocatable :: pbcprs
  INTEGER :: io_pbcprs = 0
  CHARACTER*32 :: io_pbcprs_name = ' '
        ! 'W; compression (de)acceleration of all fast ions'
        ! HEATING OF ALL FAST IONS BY COMPRESSION
  real(fp), dimension(:,:), allocatable :: pbcprss
  INTEGER, dimension(:), allocatable :: io_pbcprss
  CHARACTER*32, dimension(:), allocatable :: io_pbcprss_name
        ! 'W; compressional power to'
        ! HEATING OF BEAM IONS BY COMPRESSION BY SPECIES
  ! ... block: powers .........
    ! mibs -- fast ion specie index
    ! heating profiles f(x) = total power, watts, inside flux surface x.
    !
    ! this includes some correction terms, and various terms separated
    ! by fast ion specie index.  However, the powers seen by the thermal
    ! plasma are:
    !   pbe -- beam electron heating (summed over beam species pbes(...))
    !   pfe -- fusion product electron heating
    !                      (summed over fusion product species pbes(...))
    !
    !   pbi -- beam ion heating (summed over beam species pbes(...))
    !   pfi -- fusion product ion heating
    !                      (summed over fusion product species pbes(...))
    !   pbth -- power in fast ion thermalization source
    !
    !   pbe_nb(mj, mib) "W; beam heating to thermal electrons"
    !                       by neutral beam injector
    !   pbi_nb(mj, mib) "W; beam heating to thermal ions"
    !                       by neutral beam injector
    !   pbth_nb(mj, mib) "W; beam thermalization power"
    !                        by neutral beam injector
    !   pbtot_nb(mj, mib) "W; total beam power"
    !                       by neutral beam injector
    ! (the correction terms are included in the appropriate places).
    !
    ! xplasma ids: io_pbe, io_pbi, etc.
  real(fp), dimension(:,:), allocatable :: pbdfbes
  INTEGER, dimension(:), allocatable :: io_pbdfbes
  CHARACTER*32, dimension(:), allocatable :: io_pbdfbes_name
        ! 'W; power to e-, anom. diffusion adjustment'
        ! POWER TO ELECTRONS, BEAM ANOM.DIFF BY SPECIES
        ! generally a small correction; included in pbe (or pfe).
  real(fp), dimension(:,:), allocatable :: pbdfbis
  INTEGER, dimension(:), allocatable :: io_pbdfbis
  CHARACTER*32, dimension(:), allocatable :: io_pbdfbis_name
        ! 'W; power to ions, anom. diffusion adjustment'
        ! POWER TO IONS- BEAM ANOM.DIFF BY SPECIES
        ! generally a small correction; included in pbi (or pfi).
  real(fp), dimension(:), allocatable :: pbe
  INTEGER :: io_pbe = 0
  CHARACTER*32 :: io_pbe_name = ' '
        ! 'W; beam electron heating'
        ! POWER TO ELECTRONS FROM BEAM HEATING -- for electron power balance
  real(fp), dimension(:,:), allocatable :: pbes
  INTEGER, dimension(:), allocatable :: io_pbes
  CHARACTER*32, dimension(:), allocatable :: io_pbes_name
        ! 'W; fast ion electron heating'
        ! POWER TO ELECTRONS FROM BEAM HEATING, BY FAST SPECIES
  real(fp), dimension(:), allocatable :: pbi
  INTEGER :: io_pbi = 0
  CHARACTER*32 :: io_pbi_name = ' '
        ! 'W; beam ion heating'
        ! POWER TO IONS FROM BEAM SLOWING DOWN -- for ion power balance
  real(fp), dimension(:,:), allocatable :: pbis
  INTEGER, dimension(:), allocatable :: io_pbis
  CHARACTER*32, dimension(:), allocatable :: io_pbis_name
        ! 'W; fast ion -> thermal ion heating'
        ! POWER TO IONS FROM BEAM HEATING, BY SPECIES
  real(fp), dimension(:,:), allocatable :: pbohs
  INTEGER, dimension(:), allocatable :: io_pbohs
  CHARACTER*32, dimension(:), allocatable :: io_pbohs_name
        ! 'W; power, OH circuit -> fast ions'
        ! POWER PROFILE - OH CIRCUIT TO BEAM IONS BY SPECIES (informational).
  real(fp), dimension(:,:), allocatable :: pbephis
  INTEGER, dimension(:), allocatable :: io_pbephis
  CHARACTER*32, dimension(:), allocatable :: io_pbephis_name
        ! 'W; electrostatic acceleration of fast ions'
        ! POWER PROFILE -- radial electrostatic field to beam ions by species
        ! (informational).
  real(fp), dimension(:,:), allocatable :: pbephis_ntr
  INTEGER, dimension(:), allocatable :: io_pbephis_ntr
  CHARACTER*32, dimension(:), allocatable :: io_pbephis_ntr_name
        ! 'W; electrostatic acceleration of fast ions'
        ! POWER PROFILE -- radial electrostatic field to beam ions by species
        ! (informational) (part from NUTRAV only).
  real(fp), dimension(:,:), allocatable :: pbephis_dep
  INTEGER, dimension(:), allocatable :: io_pbephis_dep
  CHARACTER*32, dimension(:), allocatable :: io_pbephis_dep_name
        ! 'W; electrostatic accel (dep. FLR adjust)'
        ! POWER PROFILE -- radial electrostatic field to beam ions by species
        ! (informational).  (Deposition calc. only).
  real(fp), dimension(:), allocatable :: pbth
  INTEGER :: io_pbth = 0
  CHARACTER*32 :: io_pbth_name = ' '
        ! 'W; beam thermalization power'
        ! POWER TO IONS DUE TO BEAM ION THERMALIZATIONS -- for ion power bal.
        ! note this includes fusion product thermalization as well.
  real(fp), dimension(:,:), allocatable :: pbths
  INTEGER, dimension(:), allocatable :: io_pbths
  CHARACTER*32, dimension(:), allocatable :: io_pbths_name
        ! 'W; fast ion thermalization'
        ! POWER TO IONS DUE TO BEAM ION THERMALIZATIONS BY FAST ION SPECIES
  real(fp), dimension(:,:), allocatable :: pbe_nb
  INTEGER, dimension(:), allocatable :: io_pbe_nb
  CHARACTER*32, dimension(:), allocatable :: io_pbe_nb_name
        ! 'W; beam heating to thermal electrons'
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: pbi_nb
  INTEGER, dimension(:), allocatable :: io_pbi_nb
  CHARACTER*32, dimension(:), allocatable :: io_pbi_nb_name
        ! 'W; beam heating to thermal ions'
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: pbth_nb
  INTEGER, dimension(:), allocatable :: io_pbth_nb
  CHARACTER*32, dimension(:), allocatable :: io_pbth_nb_name
        ! 'W; beam thermalization power'
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: pbtot_nb
  INTEGER, dimension(:), allocatable :: io_pbtot_nb
  CHARACTER*32, dimension(:), allocatable :: io_pbtot_nb_name
        ! 'W; total beam power'
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: pbthr2s
  INTEGER, dimension(:), allocatable :: io_pbthr2s
  CHARACTER*32, dimension(:), allocatable :: io_pbthr2s_name
        ! 'W; ion thermalization rotation adjustment'
        ! ROTATION/ASSYMETRY PERTUBATION TO PBTH BY FAST ION SPECIES
        ! a small correction term, included in pbth
  real(fp), dimension(:), allocatable :: pfe
  INTEGER :: io_pfe = 0
  CHARACTER*32 :: io_pfe_name = ' '
        ! 'W; fusion product electron heating'
        ! POWER TO ELECTRONS FROM FUSION ION HEATING -- for electron power bal.
  real(fp), dimension(:), allocatable :: pfi
  INTEGER :: io_pfi = 0
  CHARACTER*32 :: io_pfi_name = ' '
        ! 'W; fusion product thermal ion heating'
        ! POWER TO IONS FROM FUSION ION SLOWING DOWN -- for ion power bal.
  real(fp), dimension(:,:), allocatable :: orbrzv_Prf
  INTEGER, dimension(:), allocatable :: io_orbrzv_Prf
  CHARACTER*32, dimension(:), allocatable :: io_orbrzv_Prf_name
        ! 'W; RF-heating on fast ions from mcrf'
        ! power profile calculated by mcrf
  real(fp), dimension(:,:), allocatable :: orbrzv_toric_Prf
  INTEGER, dimension(:), allocatable :: io_orbrzv_toric_Prf
  CHARACTER*32, dimension(:), allocatable :: io_orbrzv_toric_Prf_name
        ! 'W; RF-heating on fast ions from toric'
        ! power profile calculated by toric
  ! ... block: currents .........
    ! mibs -- fast ion specie index
    ! driven current profiles
    !
    !   curjbeam and curjfusi (or their <J.B> form equivalents)
    !   will be needed for the parallel Ohms Law expression in the
    !   poloidal field diffusion equation; the rest are just for
    !   information.
  real(fp), dimension(:), allocatable :: curjbeam
  INTEGER :: io_curjbeam = 0
  CHARACTER*32 :: io_curjbeam_name = ' '
        ! 'A; beam ion current drive (shielded)'
  real(fp), dimension(:), allocatable :: curjfusi
  INTEGER :: io_curjfusi = 0
  CHARACTER*32 :: io_curjfusi_name = ' '
        ! 'A; fusion product ion current drive (shielded)'
  real(fp), dimension(:,:), allocatable :: curjs
  INTEGER, dimension(:), allocatable :: io_curjs
  CHARACTER*32, dimension(:), allocatable :: io_curjs_name
        ! 'A; shielded fast ion driven current'
  real(fp), dimension(:,:), allocatable :: ucurjs
  INTEGER, dimension(:), allocatable :: io_ucurjs
  CHARACTER*32, dimension(:), allocatable :: io_ucurjs_name
        ! 'A; unshielded fast ion driven current'
  real(fp), dimension(:), allocatable :: ucurjco
  INTEGER :: io_ucurjco = 0
  CHARACTER*32 :: io_ucurjco_name = ' '
        ! 'A; unshielded current driven by co-injected beams'
  real(fp), dimension(:), allocatable :: ucurjctr
  INTEGER :: io_ucurjctr = 0
  CHARACTER*32 :: io_ucurjctr_name = ' '
        ! 'A; unshielded current driven by counter-injected beams'
        ! unshielded current driven by "ctr" injected beams, amps
  real(fp), dimension(:,:), allocatable :: ucurjnb
  INTEGER, dimension(:), allocatable :: io_ucurjnb
  CHARACTER*32, dimension(:), allocatable :: io_ucurjnb_name
        ! 'A; unshielded fast ion driven current by beam'
  real(fp), dimension(:,:), allocatable :: curjnb
  INTEGER, dimension(:), allocatable :: io_curjnb
  CHARACTER*32, dimension(:), allocatable :: io_curjnb_name
        ! 'A; shielded fast ion driven current by beam'
  real(fp), dimension(:), allocatable :: xjbfac
  INTEGER :: io_xjbfac = 0
  CHARACTER*32 :: io_xjbfac_name = ' '
        ! '-; species averaged neoclassical current shielding factor'
        ! --> average fast ion CURRENT SHIELDING FACTOR (btw 0 and 1)
        ! (1-Zb/ZEFF FOR classical, + NEOCLASSICAL TERMS as requested)
        ! (classical also known as Okhawa current)
        ! cf NMCURB -- member of "miscellaneous" input control structure
        !  (averaging is weighted by unshielded current absolute value)
  real(fp), dimension(:,:), allocatable :: xjbfacs
  INTEGER, dimension(:), allocatable :: io_xjbfacs
  CHARACTER*32, dimension(:), allocatable :: io_xjbfacs_name
        ! '-; neoclassical current shielding factor, by species'
  ! ... block: neutral_sources .........
    ! mibs -- fast ion specie index
    ! mig -- thermal gas index
    ! SORBN0 is a flux surface averaged volume source that might be
    ! used by a 1d thermal neutral gas model to estimate flux surface
    ! averaged halo neutral densities.  Non-flux-surface averaged 2d
    ! information is also available, to drive a 2d neutral gas model.
  real(fp), dimension(:,:), allocatable :: sbcx0s
  INTEGER, dimension(:), allocatable :: io_sbcx0s
  CHARACTER*32, dimension(:), allocatable :: io_sbcx0s_name
        ! '#/S; CX fast neutral source'
        ! SCE FCN: FAST NEUTRALS BORN BY CX INSIDE PLASMA BY FAST ION SPECIES
        ! (information only; MC code computes recapture and loss)
  real(fp), dimension(:,:), allocatable :: sorbn0
  INTEGER, dimension(:), allocatable :: io_sorbn0
  CHARACTER*32, dimension(:), allocatable :: io_sorbn0_name
        ! '#/S; thermal neutral source'
        ! Thermal NEUTRAL SOURCE DUE TO deposition or cx recapture of fast
        ! neutrals, summed over species -- ** input to neutral gas model **
  real(fp), dimension(:,:,:), allocatable :: sorbn0s
  INTEGER, dimension(:,:), allocatable :: io_sorbn0s
  CHARACTER*32, dimension(:,:), allocatable :: io_sorbn0s_name
        ! '#/S; thermal neutral source'
        ! Thermal NEUTRAL SOURCE DUE TO deposition or cx recapture of fast
        ! neutrals, broken down by fast neutral specie
  real(fp), dimension(:,:,:), allocatable :: sorbn0s_dep
  INTEGER, dimension(:,:), allocatable :: io_sorbn0s_dep
  CHARACTER*32, dimension(:,:), allocatable :: io_sorbn0s_dep_name
        ! '#/S; thermal neutral source (dep only)'
        ! Thermal NEUTRAL SOURCE DUE TO beam deposition *only*,
        ! broken down by fast neutral specie
  real(fp), dimension(:,:,:), allocatable :: sn0s3dh
  INTEGER, dimension(:,:), allocatable :: io_sn0s3dh
  CHARACTER*32, dimension(:,:), allocatable :: io_sn0s3dh_name
        ! '#/S; 3dhalo thermal neutral source zone av'
        ! Thermal NEUTRAL SOURCE DUE TO beam deposition *only*,
        ! broken down by fast neutral specie
  ! ... block: recapture .........
    ! mibs -- fast ion specie index
    ! mig -- thermal gas index
    ! profiles showing processes for recapture of charge exchange
    ! fast neutrals.
  real(fp), dimension(:,:), allocatable :: sbcxrbs
  INTEGER, dimension(:), allocatable :: io_sbcxrbs
  CHARACTER*32, dimension(:), allocatable :: io_sbcxrbs_name
        ! '#/S; fast ion induced fast neutral recapture'
        ! FAST NEUTRAL RECAPTURE BY BEAM-BEAM CX OR IONIZIZATION
  real(fp), dimension(:,:), allocatable :: sbcxres
  INTEGER, dimension(:), allocatable :: io_sbcxres
  CHARACTER*32, dimension(:), allocatable :: io_sbcxres_name
        ! '#/S; fast neutral recapture: e- impact'
        ! FAST NEUTRAL RECAPTURE BY IONIZATION (ELECTRONS)
  real(fp), dimension(:,:), allocatable :: sbcxris
  INTEGER, dimension(:), allocatable :: io_sbcxris
  CHARACTER*32, dimension(:), allocatable :: io_sbcxris_name
        ! '#/S; fast neutral recapture: th. ion impact'
        ! FAST NEUTRAL RECAPTURE BY IONIZATION (collision w/ non-impurity
        ! thermal ions)
  real(fp), dimension(:,:), allocatable :: sbcxrls
  INTEGER, dimension(:), allocatable :: io_sbcxrls
  CHARACTER*32, dimension(:), allocatable :: io_sbcxrls_name
        ! '#/S; fast neutral recapture: Lorentz ioniz.'
        ! FAST NEUTRAL RECAPTURE BY Lorentz IONIZATION
  real(fp), dimension(:,:,:), allocatable :: sbcxrxs
  INTEGER, dimension(:,:), allocatable :: io_sbcxrxs
  CHARACTER*32, dimension(:,:), allocatable :: io_sbcxrxs_name
        ! '#/S; fast neutral recapture: th ion cx'
        ! FAST NEUTRAL RECAPTURE BY CHARGE EXCHANGE W/THERMAL ION SPECIES
  real(fp), dimension(:,:), allocatable :: sbcxrzs
  INTEGER, dimension(:), allocatable :: io_sbcxrzs
  CHARACTER*32, dimension(:), allocatable :: io_sbcxrzs_name
        ! '#/S; fast neutral recapture: impurity impact'
        ! FAST NEUTRAL RECAPTURE BY IONIZATION (IMPURITIES)
  ! ... block: neutral_sinks .........
    ! mibs -- fast ion specie index
    ! mig -- thermal gas index
    !
    ! sink rates for thermal neutral gases-- flux surface averages.
    ! These should be used by a 1d thermal neutral gas model to
    ! calculate attenuation of neutral density due to the effects
    ! of fast ions.  For many experiments, this is a dominant effect
    ! on core neutral transport!
    !
    ! the profiles sinbn0 and sinb0i are flux surface averaged and
    ! suitable for use by a "1d" or flux surface averaged thermal
    ! neutral gas model; the code also outputs this data not flux
    ! surface averaged on a 2d grid: "SNBN02" and "SNB0I2"; the 2d
    ! outputs of the fast ion code are described in detail elsewhere.
    !
    ! also provided are profiles broken down by fast ion specie, for
    ! informational purposes.
  real(fp), dimension(:,:), allocatable :: sinb0i
  INTEGER, dimension(:), allocatable :: io_sinb0i
  CHARACTER*32, dimension(:), allocatable :: io_sinb0i_name
        ! '/S; fast ion neutral ionization rate'
        ! SINK rate OF NEUTRALS DUE TO fast ion IMPACT IONIZATION (sum).
  real(fp), dimension(:,:,:), allocatable :: sinb0is
  INTEGER, dimension(:,:), allocatable :: io_sinb0is
  CHARACTER*32, dimension(:,:), allocatable :: io_sinb0is_name
        ! '/S; fast ion neutral ionization rate'
        ! SINK rate OF NEUTRALS DUE TO fast ion IMPACT IONIZATION
        ! on backgroung plasma
  real(fp), dimension(:,:), allocatable :: sinbn0
  INTEGER, dimension(:), allocatable :: io_sinbn0
  CHARACTER*32, dimension(:), allocatable :: io_sinbn0_name
        ! '/S; fast ion CX rate w/thermal neutrals'
        ! NEUTRAL SINK rate DUE TO fast ion CHARGE EXCHANGE (sum).
  real(fp), dimension(:,:,:), allocatable :: sinbn0s
  INTEGER, dimension(:,:), allocatable :: io_sinbn0s
  CHARACTER*32, dimension(:,:), allocatable :: io_sinbn0s_name
        ! '/S; fast ion CX rate with neutrals'
        ! NEUTRAL SINK rate DUE TO fast ion CHARGE EXCHANGE
        ! on backgroung plasma
  ! ... block: sources .........
    ! mibs -- fast ion specie index
    ! mig -- thermal gas index
    !
    ! electron source profile, ion source profiles f(x) = no. of
    ! ptcls per second inside surface "x".
    !
    ! transport model particle continuity equations should
    ! use SORBE and SORBH; SORBTH and the
    ! breakdown by fast species profiles are informational.
  real(fp), dimension(:), allocatable :: sorbe
  INTEGER :: io_sorbe = 0
  CHARACTER*32 :: io_sorbe_name = ' '
        ! '#/S; electron source (fast neutral ionization)'
        ! ELECTRON SOURCE DUE TO BEAM DEPOSITION & CX total
        ! electrons direct from fast neutrals only; see also SORBN0
  real(fp), dimension(:,:), allocatable :: sorbes
  INTEGER, dimension(:), allocatable :: io_sorbes
  CHARACTER*32, dimension(:), allocatable :: io_sorbes_name
        ! '#/S; electron source (fast neutral ionization)'
        ! ELECTRON SOURCE DUE TO BEAM DEPOSITION & CX, each fast specie
        ! electrons direct from fast neutrals only; see also SORBN0
  real(fp), dimension(:,:), allocatable :: sorbes_ntr
  INTEGER, dimension(:), allocatable :: io_sorbes_ntr
  CHARACTER*32, dimension(:), allocatable :: io_sorbes_ntr_name
        ! '#/S; electron source (fast neutral ionization)'
        ! ELECTRON SOURCE DUE TO BEAM DEPOSITION & CX, each fast specie
        ! electrons direct from fast neutrals only; see also SORBN0
        ! came only from NUTRAV
  real(fp), dimension(:,:), allocatable :: sorbes_dep
  INTEGER, dimension(:), allocatable :: io_sorbes_dep
  CHARACTER*32, dimension(:), allocatable :: io_sorbes_dep_name
        ! '#/S; electron source (beam deposition)'
        ! ELECTRON SOURCE DUE TO BEAM DEPOSITION only, each fast specie
        ! electrons direct from fast neutrals only; see also SORBN0
  real(fp), dimension(:,:), allocatable :: sorbh
  INTEGER, dimension(:), allocatable :: io_sorbh
  CHARACTER*32, dimension(:), allocatable :: io_sorbh_name
        ! '#/S; CX adjusted thermalization source'
        ! net ION SOURCE/SINK DUE TO BEAM thermalization, less CH-EX loss
        ! sum over all fast ion species; for recapture of thermal neutrals
        ! (SORBN0) created by charge exchange, use a thermal neutral gas model.
  real(fp), dimension(:,:,:), allocatable :: sorbhs
  INTEGER, dimension(:,:), allocatable :: io_sorbhs
  CHARACTER*32, dimension(:,:), allocatable :: io_sorbhs_name
        ! '#/S; CX adjusted thermalization source'
        ! net ion source/sink due to beam, indexed by fast ion specie
  real(fp), dimension(:,:,:), allocatable :: sorbhs_ntr
  INTEGER, dimension(:,:), allocatable :: io_sorbhs_ntr
  CHARACTER*32, dimension(:,:), allocatable :: io_sorbhs_ntr_name
        ! '#/S; CX adjusted thermalization source'
        ! net ion source/sink due to beam, indexed by fast ion specie
        ! came only from NUTRAV
  real(fp), dimension(:,:,:), allocatable :: sorbhs_dep
  INTEGER, dimension(:,:), allocatable :: io_sorbhs_dep
  CHARACTER*32, dimension(:,:), allocatable :: io_sorbhs_dep_name
        ! '#/S; CX sink due to beam deposition'
        ! net ion sink due to beam ch-ex depo., indexed by fast ion specie
  real(fp), dimension(:), allocatable :: sorbth
  INTEGER :: io_sorbth = 0
  CHARACTER*32 :: io_sorbth_name = ' '
        ! '#/S; thermalization source summed over species'
        ! ION SOURCE DUE TO BEAM THERMALIZATION, summed over all species
  real(fp), dimension(:,:), allocatable :: sorbths
  INTEGER, dimension(:), allocatable :: io_sorbths
  CHARACTER*32, dimension(:), allocatable :: io_sorbths_name
        ! '#/S; thermalization source'
        ! ION SOURCE DUE TO BEAM THERMALIZATION, indexed by fast species
  ! ... block: fokker_planck .........
    ! profiles related to the collision operator:
    !
    !   --slowing down times of newly deposited fast ions
    !   --pitch angle scattering times of newly deposited fast ions
    !
    !   different profiles are given for each category of beam
    !   ion (co-injected with Ip, ctr-injected, full/half/third
    !   energy).  For fusion product species, only taubpa0 and
    !   taubsl0 are computed.
    !
    ! also computed:
    !   fraction of ion heating going to impurities
    !      PBI*[this fraction] = power to impurity species
    !   fraction of pitch angle scattering caused by impurities
  real(fp), dimension(:,:), allocatable :: taubpa0
  INTEGER, dimension(:), allocatable :: io_taubpa0
  CHARACTER*32, dimension(:), allocatable :: io_taubpa0_name
        ! 'S; P.A. scattering time at deposition (avg)'
        ! p.a. scattering time (secs), average over deposition
  real(fp), dimension(:,:), allocatable :: taubpaco1
  INTEGER, dimension(:), allocatable :: io_taubpaco1
  CHARACTER*32, dimension(:), allocatable :: io_taubpaco1_name
        ! 'S; P.A. scattering time at dep. (full E co)'
        ! p.a. scattering -- full energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubpaco2
  INTEGER, dimension(:), allocatable :: io_taubpaco2
  CHARACTER*32, dimension(:), allocatable :: io_taubpaco2_name
        ! 'S; P.A. scattering time at dep. (half E co)'
        ! p.a. scattering -- half energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubpaco3
  INTEGER, dimension(:), allocatable :: io_taubpaco3
  CHARACTER*32, dimension(:), allocatable :: io_taubpaco3_name
        ! 'S; P.A. scattering time at dep. (1/3 E co)'
        ! p.a. scattering -- 1/3 energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubpactr1
  INTEGER, dimension(:), allocatable :: io_taubpactr1
  CHARACTER*32, dimension(:), allocatable :: io_taubpactr1_name
        ! 'S; P.A. scattering time at dep. (full E ctr)'
        ! p.a. scattering -- full energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubpactr2
  INTEGER, dimension(:), allocatable :: io_taubpactr2
  CHARACTER*32, dimension(:), allocatable :: io_taubpactr2_name
        ! 'S; P.A. scattering time at dep. (half E ctr)'
        ! p.a. scattering -- half energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubpactr3
  INTEGER, dimension(:), allocatable :: io_taubpactr3
  CHARACTER*32, dimension(:), allocatable :: io_taubpactr3_name
        ! 'S; P.A. scattering time at dep. (1/3 E ctr)'
        ! p.a. scattering -- 1/3 energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubsl0
  INTEGER, dimension(:), allocatable :: io_taubsl0
  CHARACTER*32, dimension(:), allocatable :: io_taubsl0_name
        ! 'S; slowing down time at deposition (avg)'
        ! slowing down time (secs), average over deposition
  real(fp), dimension(:,:), allocatable :: taubslco1
  INTEGER, dimension(:), allocatable :: io_taubslco1
  CHARACTER*32, dimension(:), allocatable :: io_taubslco1_name
        ! 'S; slowing down time at dep. (full E co)'
        ! slowing down -- full energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubslco2
  INTEGER, dimension(:), allocatable :: io_taubslco2
  CHARACTER*32, dimension(:), allocatable :: io_taubslco2_name
        ! 'S; slowing down time at dep. (half E co)'
        ! slowing down -- half energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubslco3
  INTEGER, dimension(:), allocatable :: io_taubslco3
  CHARACTER*32, dimension(:), allocatable :: io_taubslco3_name
        ! 'S; slowing down time at dep. (1/3 E co)'
        ! slowing down -- 1/3 energy "co" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubslctr1
  INTEGER, dimension(:), allocatable :: io_taubslctr1
  CHARACTER*32, dimension(:), allocatable :: io_taubslctr1_name
        ! 'S; slowing down time at dep. (full E ctr)'
        ! slowing down -- full energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubslctr2
  INTEGER, dimension(:), allocatable :: io_taubslctr2
  CHARACTER*32, dimension(:), allocatable :: io_taubslctr2_name
        ! 'S; slowing down time at dep. (half E ctr)'
        ! slowing down -- half energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: taubslctr3
  INTEGER, dimension(:), allocatable :: io_taubslctr3
  CHARACTER*32, dimension(:), allocatable :: io_taubslctr3_name
        ! 'S; slowing down time at dep. (1/3 E ctr)'
        ! slowing down -- 1/3 energy "ctr" beam ions, sec
  real(fp), dimension(:,:), allocatable :: paxfrc
  INTEGER, dimension(:), allocatable :: io_paxfrc
  CHARACTER*32, dimension(:), allocatable :: io_paxfrc_name
        ! '-; fraction of P.A. scattering on impurities'
        ! FRACTION OF PITCH ANGLE SCATTERING DUE TO IMPURITIES, BY FAST SPECIES
        ! covers beams and fusion product species
  real(fp), dimension(:,:), allocatable :: pbxfrc
  INTEGER, dimension(:), allocatable :: io_pbxfrc
  CHARACTER*32, dimension(:), allocatable :: io_pbxfrc_name
        ! '-; fraction of ion heating on impurities'
        ! FRACTION OF ION HEATING GOING TO IMPURITY SPECIES, BY FAST SPECIES
        ! covers beams and fusion product species
  ! ... block: radial_current .........
    ! these are boundary oriented radial current profiles
    ! these are informational outputs, and apt to be noisy.
    ! DMC May 2010: anomalous transport diffusivity and pinch
    ! velocity "as observed" output profiles added.
  real(fp), dimension(:), allocatable :: tirabd
  INTEGER :: io_tirabd = 0
  CHARACTER*32 :: io_tirabd_name = ' '
        ! 'A(radial); anomalous diffusion current'
        ! BEAM ION RADIAL CURRENT: anomalous diffusion
  real(fp), dimension(:,:), allocatable :: tirabds
  INTEGER, dimension(:), allocatable :: io_tirabds
  CHARACTER*32, dimension(:), allocatable :: io_tirabds_name
        ! 'A(radial); anomalous diffusion current'
        ! BEAM ION RADIAL CURRENT: anomalous diffusion, by species
  real(fp), dimension(:), allocatable :: tirfsh
  INTEGER :: io_tirfsh = 0
  CHARACTER*32 :: io_tirfsh_name = ' '
        ! 'A(radial); fishbone current'
        ! BEAM ION RADIAL CURRENT: fishbones
  real(fp), dimension(:,:), allocatable :: tirfshs
  INTEGER, dimension(:), allocatable :: io_tirfshs
  CHARACTER*32, dimension(:), allocatable :: io_tirfshs_name
        ! 'A(radial); fishbone current'
        ! BEAM ION RADIAL CURRENT: fishbones, by species
  real(fp), dimension(:), allocatable :: tirorb
  INTEGER :: io_tirorb = 0
  CHARACTER*32 :: io_tirorb_name = ' '
        ! 'A(radial); orbit evolution current'
        ! BEAM ION RADIAL CURRENT ACROSS SURFACES: classical orbiting,
        ! computed by differencing orbit averages
  real(fp), dimension(:,:), allocatable :: tirorbs
  INTEGER, dimension(:), allocatable :: io_tirorbs
  CHARACTER*32, dimension(:), allocatable :: io_tirorbs_name
        ! 'A(radial); orbit evolution current'
        ! BEAM ION RADIAL CURRENT ACROSS SURFACES: classical orbiting,
        ! computed by differencing orbit averages, by species
  real(fp), dimension(:), allocatable :: tirrip
  INTEGER :: io_tirrip = 0
  CHARACTER*32 :: io_tirrip_name = ' '
        ! 'A(radial); TF ripple loss current'
        ! BEAM ION RADIAL CURRENT: ripple transport
  real(fp), dimension(:,:), allocatable :: tirrips
  INTEGER, dimension(:), allocatable :: io_tirrips
  CHARACTER*32, dimension(:), allocatable :: io_tirrips_name
        ! 'A(radial); TF ripple loss current'
        ! BEAM ION RADIAL CURRENT: ripple transport, by species
  real(fp), dimension(:,:), allocatable :: difbx
  INTEGER, dimension(:), allocatable :: io_difbx
  CHARACTER*32, dimension(:), allocatable :: io_difbx_name
        ! 'cm2/sec; average anomalous diffusivity'
  real(fp), dimension(:,:), allocatable :: velbx
  INTEGER, dimension(:), allocatable :: io_velbx
  CHARACTER*32, dimension(:), allocatable :: io_velbx_name
        ! 'cm/sec; average anomalous radial velocity'
  ! ... block: torques .........
    ! mibs -- fast ion specie index
    ! generally, this block contains torques applied to thermal
    ! plasma, due to various processes.  But:
    !   tqohb is a torque from the OH circuit to *fast* ions
    !   tqrpl is a torque from JxB due to ripple loss current
    ! and the coupling of these torques to the thermal plasma
    ! are considered unknown.
    !
    ! A net total torque
    !   tqbsum = collisional torque + JxB torque + thermalization torque
    !            + anom. diffusivity transfer - charge exchange loss
    !
    !          = tqbe + tqbi + tqjxb + tqbth + tqdifb - tqbcx
    !
    ! integrated profiles, f(x) = total torque inside flux surface "x".
    !
    ! for plasma momentum balance: the recommended momentum source
    ! is   tqbsum + c1*tqrpl - c2*tqohb   where c1 and c2 are
    ! adjustable parameters (e.g. TRANSP usually assumes c1=c2=0.0_fp)
    !
    ! also in the block is a dimensionless profile "tqbxf" giving
    ! the fraction of ion collisional torque going to all impurities.
  real(fp), dimension(:), allocatable :: tqbsum
  INTEGER :: io_tqbsum = 0
  CHARACTER*32 :: io_tqbsum_name = ' '
        ! 'Nt*m; summed fast ion torque'
        ! = tqbcx + tqbe + tqbi + tqbth + tqdifb + tqjxb
        ! caution: tqohb and tqrpl ** not included **
  real(fp), dimension(:), allocatable :: tqbcx
  INTEGER :: io_tqbcx = 0
  CHARACTER*32 :: io_tqbcx_name = ' '
        ! 'Nt*m; beam charge exchange loss torque'
        ! TORQUE FROM IONS DUE TO CX WITH FAST NEUTRALS **loss term**
  real(fp), dimension(:,:), allocatable :: tqbcxs
  INTEGER, dimension(:), allocatable :: io_tqbcxs
  CHARACTER*32, dimension(:), allocatable :: io_tqbcxs_name
        ! 'Nt*m; beam charge exchange loss torque'
        ! TORQUE FROM IONS DUE TO CX WITH FAST NEUTRALS **loss term**
        ! by fast specie.
  real(fp), dimension(:,:), allocatable :: tqbcxs_dep
  INTEGER, dimension(:), allocatable :: io_tqbcxs_dep
  CHARACTER*32, dimension(:), allocatable :: io_tqbcxs_dep_name
        ! 'Nt*m; beam CX loss torque during deposition'
        ! TORQUE FROM IONS DUE TO CX WITH FAST NEUTRALS **loss term**
        ! by fast specie. (deposition only)
  real(fp), dimension(:), allocatable :: tqbe
  INTEGER :: io_tqbe = 0
  CHARACTER*32 :: io_tqbe_name = ' '
        ! 'Nt*m; beam collisional torque to electrons'
        ! TORQUE TO ELECTRONS FROM BEAM SLOWING DOWN
  real(fp), dimension(:,:), allocatable :: tqbes
  INTEGER, dimension(:), allocatable :: io_tqbes
  CHARACTER*32, dimension(:), allocatable :: io_tqbes_name
        ! 'Nt*m; beam collisional torque to electrons'
        ! TORQUE TO ELECTRONS FROM BEAM SLOWING DOWN
        ! by fast specie.
  real(fp), dimension(:), allocatable :: tqbi
  INTEGER :: io_tqbi = 0
  CHARACTER*32 :: io_tqbi_name = ' '
        ! 'Nt*m; beam collisional torque to ions'
        ! TORQUE TO IONS FROM BEAM SLOWING DOWN
  real(fp), dimension(:,:), allocatable :: tqbis
  INTEGER, dimension(:), allocatable :: io_tqbis
  CHARACTER*32, dimension(:), allocatable :: io_tqbis_name
        ! 'Nt*m; beam collisional torque to ions'
        ! TORQUE TO IONS FROM BEAM SLOWING DOWN
        ! by fast specie.
  real(fp), dimension(:,:), allocatable :: tqbco_nb
  INTEGER, dimension(:), allocatable :: io_tqbco_nb
  CHARACTER*32, dimension(:), allocatable :: io_tqbco_nb_name
        ! 'Nt*m; beam collisional torque to ions and electrons'
        ! TORQUE TO IONS AND ELECTRONS FROM BEAM SLOWING DOWN
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: tqbth_nb
  INTEGER, dimension(:), allocatable :: io_tqbth_nb
  CHARACTER*32, dimension(:), allocatable :: io_tqbth_nb_name
        ! 'Nt*m; beam thermalization torque'
        ! by neutral beam injector
  real(fp), dimension(:,:), allocatable :: tqtot_nb
  INTEGER, dimension(:), allocatable :: io_tqtot_nb
  CHARACTER*32, dimension(:), allocatable :: io_tqtot_nb_name
        ! 'Nt*m; beam total torque'
        ! tqbco_nb+tqjxb_nb+tqbth_nb
        ! by neutral beam injector
  real(fp), dimension(:), allocatable :: tqbth
  INTEGER :: io_tqbth = 0
  CHARACTER*32 :: io_tqbth_name = ' '
        ! 'Nt*m; beam thermalization torque'
        ! TORQUE TO IONS FROM BEAM THERMALIZATION
  real(fp), dimension(:,:), allocatable :: tqbths
  INTEGER, dimension(:), allocatable :: io_tqbths
  CHARACTER*32, dimension(:), allocatable :: io_tqbths_name
        ! 'Nt*m; beam thermalization torque'
        ! TORQUE TO IONS FROM BEAM THERMALIZATION
        ! by fast specie.
  real(fp), dimension(:), allocatable :: tqbxf
  INTEGER :: io_tqbxf = 0
  CHARACTER*32 :: io_tqbxf_name = ' '
        ! '-; frac of collisional torque to impurities'
        ! FRACTION OF TQBI GOING TO IMPURITY SPECIES
  real(fp), dimension(:,:), allocatable :: tqbxfs
  INTEGER, dimension(:), allocatable :: io_tqbxfs
  CHARACTER*32, dimension(:), allocatable :: io_tqbxfs_name
        ! '-; frac of collisional torque to impurities'
        ! FRACTION OF TQBI GOING TO IMPURITY SPECIES
        ! by fast specie.
  real(fp), dimension(:), allocatable :: tqdifb
  INTEGER :: io_tqdifb = 0
  CHARACTER*32 :: io_tqdifb_name = ' '
        ! 'Nt*m; torque assoc. with anom. diffusion'
        ! MOMENTUM DISSIPATION - BEAM ANOM. DIFFUSION
  real(fp), dimension(:,:), allocatable :: tqdifbs
  INTEGER, dimension(:), allocatable :: io_tqdifbs
  CHARACTER*32, dimension(:), allocatable :: io_tqdifbs_name
        ! 'Nt*m; torque assoc. with anom. diffusion'
        ! MOMENTUM DISSIPATION - BEAM ANOM. DIFF BY SPECIES
        ! by fast specie.
  real(fp), dimension(:), allocatable :: tqjxb
  INTEGER :: io_tqjxb = 0
  CHARACTER*32 :: io_tqjxb_name = ' '
        ! 'Nt*m; JxB torque'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
  real(fp), dimension(:,:), allocatable :: tqjxbs
  INTEGER, dimension(:), allocatable :: io_tqjxbs
  CHARACTER*32, dimension(:), allocatable :: io_tqjxbs_name
        ! 'Nt*m; JxB torque'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
        ! by fast specie.
  real(fp), dimension(:,:), allocatable :: tqjxbs_ntr
  INTEGER, dimension(:), allocatable :: io_tqjxbs_ntr
  CHARACTER*32, dimension(:), allocatable :: io_tqjxbs_ntr_name
        ! 'Nt*m; JxB torque'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA (part from NUTRAV only)
        ! by fast specie.
  real(fp), dimension(:,:), allocatable :: tqjxb_nb
  INTEGER, dimension(:), allocatable :: io_tqjxb_nb
  CHARACTER*32, dimension(:), allocatable :: io_tqjxb_nb_name
        ! 'Nt*m; JxB torque'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
        ! by beam line
  real(fp), dimension(:,:), allocatable :: tqjxbnb_ntr
  INTEGER, dimension(:), allocatable :: io_tqjxbnb_ntr
  CHARACTER*32, dimension(:), allocatable :: io_tqjxbnb_ntr_name
        ! 'Nt*m; JxB torque'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
        ! by beam line (part from NUTRAV only)
  real(fp), dimension(:,:), allocatable :: tqjxbs_dep
  INTEGER, dimension(:), allocatable :: io_tqjxbs_dep
  CHARACTER*32, dimension(:), allocatable :: io_tqjxbs_dep_name
        ! 'Nt*m; JxB torque (deposition FLR step)'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
        ! by fast specie.  (deposition only)
  real(fp), dimension(:,:), allocatable :: tqjxb_nb_dep
  INTEGER, dimension(:), allocatable :: io_tqjxb_nb_dep
  CHARACTER*32, dimension(:), allocatable :: io_tqjxb_nb_dep_name
        ! 'Nt*m; JxB torque (deposition FLR step)'
        ! BEAM ION DEP/ORBIT J X B TORQUE TO PLASMA
        ! by beam line (deposition only)
  real(fp), dimension(:), allocatable :: tqohb
  INTEGER :: io_tqohb = 0
  CHARACTER*32 :: io_tqohb_name = ' '
        ! 'Nt*m; torque to beam ions from OH circuit'
        ! OHMIC TORQUE TO BEAM IONS =? -TORQUE TO PLASMA
  real(fp), dimension(:,:), allocatable :: tqohbs
  INTEGER, dimension(:), allocatable :: io_tqohbs
  CHARACTER*32, dimension(:), allocatable :: io_tqohbs_name
        ! 'Nt*m; torque to beam ions from OH circuit'
        ! OHMIC TORQUE TO BEAM IONS =? -TORQUE TO PLASMA
        ! by fast specie.
  real(fp), dimension(:), allocatable :: tqrpl
  INTEGER :: io_tqrpl = 0
  CHARACTER*32 :: io_tqrpl_name = ' '
        ! 'Nt*m; torque due to TF ripple loss'
        ! TOTAL RIPPLE TORQUE PROFILE
        ! under assumption: radial current due to ripple
        ! creates a JxB force.  ** this could be a large number **
        ! ** not included in tqbsum **
  real(fp), dimension(:,:), allocatable :: tqrpls
  INTEGER, dimension(:), allocatable :: io_tqrpls
  CHARACTER*32, dimension(:), allocatable :: io_tqrpls_name
        ! 'Nt*m; torque due to TF ripple loss'
        ! TOTAL RIPPLE TORQUE PROFILE -- separated by fast species.
  ! ... block: frantic .........
  real(fp), dimension(:,:,:), allocatable :: nbi_den0mn
  INTEGER, dimension(:,:), allocatable :: io_nbi_den0mn
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_den0mn_name
        ! 'the neutral density #/m**3'
  real(fp), dimension(:,:,:), allocatable :: nbi_en0mn
  INTEGER, dimension(:,:), allocatable :: io_nbi_en0mn
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_en0mn_name
        ! 'KeV, the neutral temperature'
  real(fp), dimension(:,:,:), allocatable :: nbi_omg0mn
  INTEGER, dimension(:,:), allocatable :: io_nbi_omg0mn
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_omg0mn_name
        ! 'the av. neutral toroidal velocity, m/sec'
  real(fp), dimension(:,:,:), allocatable :: nbi_snimn
  INTEGER, dimension(:,:), allocatable :: io_nbi_snimn
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_snimn_name
        ! 'the ionization ion source, #/m**3/sec'
  real(fp), dimension(:,:), allocatable :: nbi_outflx
        ! 'outflux rate (#/m**2/sec)'
  real(fp), dimension(:,:), allocatable :: nbi_eescav
        ! 'average energy (KeV/ptcl)'
  real(fp), dimension(:,:), allocatable :: nbi_cexflx
        ! 'outflux rate (#/m**2/sec)'
  real(fp), dimension(:,:,:), allocatable :: nbi_snemn
  INTEGER, dimension(:,:), allocatable :: io_nbi_snemn
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_snemn_name
        ! '#/m**3/sec the electron source'
        ! to the plasma due to ionization of neutrals
  real(fp), dimension(:,:), allocatable :: nbi_scxmn
  INTEGER, dimension(:), allocatable :: io_nbi_scxmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_scxmn_name
        ! 'the CX reaction rate #/m**3/sec'
  real(fp), dimension(:,:), allocatable :: nbi_pcxmn
  INTEGER, dimension(:), allocatable :: io_nbi_pcxmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pcxmn_name
        ! 'the CX power W/m**3 lost to ions'
  real(fp), dimension(:,:), allocatable :: nbi_pnimn
  INTEGER, dimension(:), allocatable :: io_nbi_pnimn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pnimn_name
        ! 'the ionization source power W/m**3 to ions'
  real(fp), dimension(:,:), allocatable :: nbi_tqcxmn
  INTEGER, dimension(:), allocatable :: io_nbi_tqcxmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_tqcxmn_name
        ! 'the CX torque density Nt-m/m**3 lost to ions'
  real(fp), dimension(:,:), allocatable :: nbi_tqnimn
  INTEGER, dimension(:), allocatable :: io_nbi_tqnimn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_tqnimn_name
        ! 'the ioniz. source torque density  Nt-m/m**3'
        !  to ions
  real(fp), dimension(:,:), allocatable :: nbi_pfcxmn
  INTEGER, dimension(:), allocatable :: io_nbi_pfcxmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pfcxmn_name
        ! 'the CX source friction  W/m**3 to ions'
  real(fp), dimension(:,:), allocatable :: nbi_pfnimn
  INTEGER, dimension(:), allocatable :: io_nbi_pfnimn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pfnimn_name
        ! 'the ioniz. source friction  W/m**3 to ions'
  real(fp), dimension(:,:,:), allocatable :: nbi_sflx0m
  INTEGER, dimension(:,:), allocatable :: io_nbi_sflx0m
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_sflx0m_name
        ! '#/m**3/sec, the total div'
        !  of (electrons carried by) neutral atoms
        ! 1 -- influx
        ! 2 -- outflux
  real(fp), dimension(:,:,:), allocatable :: nbi_pflx0m
  INTEGER, dimension(:,:), allocatable :: io_nbi_pflx0m
  CHARACTER*32, dimension(:,:), allocatable :: io_nbi_pflx0m_name
        ! 'W/m**3, the total div'
        ! of (power carried by) neutral atoms
        ! 1 -- influx
        ! 2 -- outflux
  real(fp), dimension(:,:), allocatable :: nbi_tqflx0
  INTEGER, dimension(:), allocatable :: io_nbi_tqflx0
  CHARACTER*32, dimension(:), allocatable :: io_nbi_tqflx0_name
        ! 'Nt-m/m**3, the total div(flux)'
        !of momentum of neutral atoms"
  real(fp), dimension(:,:), allocatable :: nbi_pb0mn
  INTEGER, dimension(:), allocatable :: io_nbi_pb0mn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pb0mn_name
        ! 'W, power in halo neutral source'
        !  (loss from thermal ions)
  real(fp), dimension(:,:), allocatable :: nbi_pscmn
  INTEGER, dimension(:), allocatable :: io_nbi_pscmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_pscmn_name
        ! 'W, beam halo power'
        !  recapture - pb0_halo
        !  (net source/sink to thermal ions)
        !  (associated particle source/sink exists)
  real(fp), dimension(:,:), allocatable :: nbi_tqb0mn
  INTEGER, dimension(:), allocatable :: io_nbi_tqb0mn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_tqb0mn_name
        ! 'Nt*m, torque of halo neutral source'
        !  (loss from thermal ions)
  real(fp), dimension(:,:), allocatable :: nbi_tqscmn
  INTEGER, dimension(:), allocatable :: io_nbi_tqscmn
  CHARACTER*32, dimension(:), allocatable :: io_nbi_tqscmn_name
        ! 'Nt*m, beam halo torque'
        ! recapture - tqb0_halo
        !  (net source/sink to thermal ions)
        !  (associated particle source/sink exists)
  ! ... block: hi_z_beams .........
    ! for Z.gt.2 beams (e.g. neon beams): <Z> of partially stripped
    ! slowing down beam ions.
    ! (caution: high-Z beam model not operational-- can be revived if
    ! there is specific user interest; contact authors at PPPL).
  real(fp), dimension(:,:), allocatable :: xzbmav
  INTEGER, dimension(:), allocatable :: io_xzbmav
  CHARACTER*32, dimension(:), allocatable :: io_xzbmav_name
        ! '-; avg <Z> of beam injected impurity ions'
        ! AVERAGE CHARGE STATE OF ORBITTING HEAVY IONS
  ! ... block: cptimes .........
    ! (not to be MPI reduced-- time 1 cpu, assume load balance)
    ! (CPBMAX and CPBMIN -- Max and Min CPU times among all
    !  processes in an parallel run)
    ! CPBROOT - CPU time for the root processor
  real(fp) :: CPBPS = zero
        ! 'S; cpu time in plasma state data extraction'
  real(fp) :: CPBSTART = zero
        ! 'S; cpu time in setup (get profiles + NBSTART)'
        ! cpu time, seconds, in NBI_PROFGET + NBSTART
  real(fp) :: CPBDEP = zero
        ! 'S; cpu time in deposition (DEPALL)'
        ! cpu time, seconds, in DEPALL
  real(fp) :: CPBORB = zero
        ! 'S; cpu time in orbiting (ORBALL)'
        ! cpu time, seconds, in ORBALL
  real(fp) :: CPBFIN = zero
        ! 'S; cpu time in post-processing (NBFINISH)'
        ! cpu time, seconds, in NBFINISH
  real(fp) :: CPBMAX = zero
        ! 'S; Max total cpu time (MPI job)'
        ! Max cpu time, seconds, MPI Job
  real(fp) :: CPBMIN = zero
        ! 'S; Min total cpu time (MPI job)'
        ! Min cpu time, seconds, MPI Job
  real(fp) :: CPBROOT = zero
        ! 'S; total cpu time for root (MPI job)'
        ! cpu time for root, seconds, MPI Job
  ! ... block: power_balance .........
    ! data related to scalar power balance for beam species
    !
    ! e.g. loss terms:  (for specie "i", all in watts):
    !     bpls(i) -- bad orbit (limiter) loss
    !     bpshins(i) -- shine-through power loss (deposition)
    !     bpcxins(i) -- charge exchange loss (fast ions inside plasma)
    !     bpcxots(i) -- charge exchange loss (ions in scrapeoff region)
    !     bpfshs(i) -- power loss due to fishbones (ad hoc model)
    !     (profiles pbaepms(:,:,i) -- power due to AEP (MHD) activity)
    !     bprpls(i) -- power loss due to TF ripple (ad hoc model)
    !
    ! and power to the thermal plasma...
    !     bpths(i) -- thermalization
    !       (profiles pbe(...),pbi(...) -- beam elec. and ion heating)
    !       (profiles pfe(...),pfi(...) -- fusion product heating)
    !
    ! and power balance
    !
    !     injected power (pinjs(i))
    !    +(compressional input power) (profile pbcprss(...)
    !    +prfb(i)    RF power to fast ions (zero in Monte Carlo model)
    !    +bpohs(i)   power to fast ions from OH circuit (loop voltage)
    !    +bperrs(i)  non-physical power due to inexact orbit integration[*]
    !    -bpsts(i)   power stored (d/dt(stored energy))
    !    -(loss terms)
    !    -(rotational work done)[**]
    !    -(power to thermal plasma)
    !     --------
    !     residual (should be small).
    !
    !     pbdepmc(i) actual MC power deposited
    !
    ! total beam power deposited
    ! there is information on charge exchange powers:
    !     bpcxi0s(i) -- charge exchange fast neutral source (plasma)
    !     bpcxx0s(i) -- charge exchange source (scrapeoff region)
    !     bpcxris(i) -- recapture of neutrals launched from plasma
    !     bpcxrxs(i) -- recapture of neutrals launched in scrapeoff region
    !     bpcxins(i) -- charge exchange loss (fast ions inside plasma)
    !     bpcxots(i) -- charge exchange loss (ions in scrapeoff region)
    !     bpcxers(i) -- straight line tracker positioning error [*]
    !
    ! [*] these are measurable numerical errors, they should be quite small
    ! [**] see "torque_work" section
  real(fp), dimension(:), allocatable :: bpcxers
        ! 'W; fast ion CX tracker error'
        ! BEAM CX POWER ERROR - TRACK NOT FOUND, BY SPECIES by specie
  real(fp), dimension(:), allocatable :: bpcxi0s
        ! 'W; CX fast neutral source (plasma)'
        ! BEAM POWER CX-ING, SUBJECT TO RECAPTURE, by specie
  real(fp), dimension(:), allocatable :: bpcxins
        ! 'W; net CX loss from in plasma'
        ! BEAM POWER LOST TO CHARGE EXCHANGE-- from INSIDE, by specie
  real(fp), dimension(:), allocatable :: bpcxots
        ! 'W; net CX loss from SOL'
        ! BEAM POWER TO CX-- from OUTSIDE OF THE PLASMA BY specie
  real(fp), dimension(:), allocatable :: bpcxris
        ! 'W; fast neutral recapture (plasma sce)'
        ! BEAM CX POWER RECAPTURED (INTERIOR SOURCE PT.) by specie
  real(fp), dimension(:), allocatable :: bpcxrxs
        ! 'W; fast neutral recapture (SOL sce)'
        ! BEAM CX POWER RECAPTURED (EXTERIOR SOURCE PT.) by specie
  real(fp), dimension(:), allocatable :: bpcxx0s
        ! 'W; CX fast neutral source (SOL)'
        ! BEAM POWER CX-ING, SUBJ. TO RECAPTURE (EXTERIOR) by specie
  real(fp) :: bperr = zero
        ! 'W; orbit integration error'
        ! BEAM POWER ERROR - ORBIT EQN. CHECK
  real(fp), dimension(:), allocatable :: bperrs
        ! 'W; orbit integration error'
        ! BEAM POWER ERROR - ORBIT EQN. CHECK BY SPECIES
  real(fp), dimension(:,:,:), allocatable :: PBAEPMS
  INTEGER, dimension(:,:), allocatable :: io_PBAEPMS
  CHARACTER*32, dimension(:,:), allocatable :: io_PBAEPMS_name
        ! 'W; power profile from AEP (MHD) activity'
        ! BEAM POWER from AEP (MHD) activity BY MODE and SPECIES
  real(fp), dimension(:), allocatable :: bpfshs
        ! 'W; fishbone loss'
        ! BEAM POWER LOST TO "FISH-BONES" BY SPECIES
  real(fp), dimension(:), allocatable :: bpls
        ! 'W; orbit loss'
        ! BEAM POWER TO LIMITER BY BEAM ION SPECIES
  real(fp), dimension(:), allocatable :: bpls_dep
        ! 'W; orbit loss (deposition calculation)'
        ! BEAM POWER TO LIMITER BY BEAM ION SPECIES (deposition calc. only)
  real(fp), dimension(:), allocatable :: bpohs
        ! 'W; power from OH circuit to fast ions'
        ! TOTAL POWER FROM OH CIRCUIT TO BEAM IONS BY SPECIES
  real(fp), dimension(:), allocatable :: bprpls
        ! 'W; TF ripple loss'
        ! BEAM POWER LOSS FROM RIPPLE BY SPECIES
  real(fp), dimension(:), allocatable :: bpshins
        ! 'W; shine-thru'
        ! BEAM SHINE THRU POWER, BY BEAM SPECIES
  real(fp), dimension(:), allocatable :: bpsts
        ! 'W; d/dt(fast ion energy)'
        ! BPST, BY FAST ION SPECIES
  real(fp), dimension(:), allocatable :: bpths
        ! 'W; thermalization power'
        ! BEAM POWER THERMALIZING INTO PLASMA, BY SPECIES
  real(fp) :: pftota = zero
        ! 'W; fusion product source'
        ! GRAND TOTAL POWER IN FUSION PRODUCT SOURCE
  real(fp), dimension(:), allocatable :: bpprma
        ! 'W; estimated prompt orbit loss'
        ! EST. PROMPT ORBIT LOSS, BEAMLINE BY BEAMLINE
  real(fp), dimension(:), allocatable :: pinjs
        ! 'W; injected power'
        ! injected power (or fusion product source power), by species.
  real(fp), dimension(:), allocatable :: pinjs_res
        ! 'W; injected power'
        ! injected power (or fusion product source power) - shine through power,
        ! by species at the time step when injected power not equals zero
  real(fp), dimension(:), allocatable :: pbdepmc
        ! 'actual MC power deposited'
        ! for beam species only
        ! total beam power deposited
  real(fp), dimension(:), allocatable :: PBDEPBA
        ! 'W; beam power balance check'
        !for beam species only
        !in case when the # of deposited MC particles >= NDEP0
        !the balance will come out exact because
        !of the normalization in subroutine DEPNRM
        !(calculation is normalized to injected power, not current).
        !In case # of deposited MC particles < NDEP0
        !it won't be perfect due to statistical error.
  real(fp), dimension(:), allocatable :: bpephis
        ! 'W; power to fast ions from E_rad'
        ! watts -- power from radial electrostatic field
  real(fp), dimension(:), allocatable :: orbrzv_Prftot
        ! 'W; RF-power to fast ions from mcrf'
        ! total power calculated by mcrf
  real(fp), dimension(:), allocatable :: orbrzv_toric_Prftot
        ! 'W; RF-power to fast ions from toric'
        ! total power calculated by toric
  real(fp), dimension(:), allocatable :: orbrzv_toric_ufrnm
        ! 'RF-field renormalizer updated in this time step'
        ! toric field renormalizer updated in this time step
  ! ... block: momentum_balance .........
    ! quantities related to scalar angular momentum balance (Nt-m)
    !
    ! momentum balance of slowing down fast ion population, species "i":
    !
    !    +bphdepsgc(i)    MC momentum deposited at GC
    !    +bphdeps(i)    MC momentum deposited at point of ionization
    !    +bphdshins(i)  shine-through  (deposition)
    !    +bphinjs(i)     momentum injected
    !    +bphdpba(i)   balance check
    !    +bphohbs(i)     from OH circuit (loop voltage)
    !    -bphfshs(i)    -(fishbone loss, if ad hoc fishbone model is active)
    !    -bphrpls(i)    -(ripple loss, if ad hoc ripple loss model active)
    !    -bphorbs(i)    -(orbit loss; ion strikes wall or limiter)
    !    -bph0was(i)    -(charge exchange loss, neutrals to wall)
    !    -bphjxbs(i)    -(JxB torque due to net radial motion)
    !    -bphcols(i)    -(collisional torque to thermal plasma)
    !    -bphdfbs(i)    -(momentum effect of anomalous diffusion operator)
    !    -thermalization  -(thermalization ion source)
    !    +bpherrs(i)   +momentum due to error in orbit integrator[*]
    !    -bphsts(i)    -(d/dt(fast ion momentum))
    !     ------
    !     bphcks(k)    residual-- should be near zero but not to machine
    !                             precision.
    !
    ! [*] numerical error, should be modest, or else the fast ion orbit
    !     timestep is too large.
    !
    ! charge exchange momentum balance
    !    +bphcxls(i)    charge exchange fast neutral source momentum
    !    -bphrcas(i)    -(recapture in plasma)
    !    -bph0was(i)    -(charge exchange loss, neutrals to wall)
    !     ---------
    !     residual      should be near zero but may be noisy
    !
  real(fp), dimension(:), allocatable :: bph0was
        ! 'Nt*m; CX fast ion loss to walls'
        ! CX MOMENTUM LOST TO WALLS BY SPECIES
  real(fp), dimension(:), allocatable :: bphcks
        ! 'Nt*m;  momentum balance check'
        ! BEAM MOMENTUM BALANCE CHECK (SHOULD=0)
  real(fp), dimension(:), allocatable :: bphcols
        ! 'Nt*m; collisional torque to plasma'
        ! BEAM MOMENTUM TO PLASMA VIA COLLISIONS BY SPECIES
  real(fp), dimension(:), allocatable :: bphcons
        ! 'Nt*m*S; momentum in fast ions'
        ! BEAM MOMENTUM CONTENT, BY BEAM ION SPECIES
  real(fp), dimension(:), allocatable :: bphcxls
        ! 'Nt*m; CX of fast ions'
        ! CX LOSS OF MOMENTUM (SUBJ. TO RECAPTURE)
  real(fp), dimension(:), allocatable :: bphinjs
        ! 'Nt*m; momentum injected'
        ! MOMENTUM INJECTED BY BEAM SPECIES
  real(fp), dimension(:), allocatable :: bphdepsgc
        ! 'Nt*m; deposition at GC'
        ! MOMENTUM DEPOSITION BY BEAM SPECIES at GC
  real(fp), dimension(:), allocatable :: bphdeps
        ! 'Nt*m; deposition at  point of ionization'
        ! MOMENTUM DEPOSITION BY BEAM SPECIES at point of ionization
  real(fp), dimension(:), allocatable :: bphshins
        ! 'shine through momentum Nt*m; deposition'
        ! shine through MOMENTUM DEPOSITION BY BEAM SPECIES
  real(fp), dimension(:), allocatable :: bphdpba
        ! 'Nt*m; beam momentum balance check'
  real(fp), dimension(:), allocatable :: bphdfbs
        ! 'Nt*m; anomalous diffusion torque'
        ! BEAM ANOM DIFF. TORQUE TO PLASMA BY SPECIES
  real(fp) :: bpherr = zero
        ! 'Nt*m; orbit integration error'
        ! BEAM MOMENTUM ERROR - ORBIT EQN. CHECK
        ! summed over all fast species
  real(fp), dimension(:), allocatable :: bpherrs
        ! 'Nt*m; orbit integration error'
        ! BEAM MOMENTUM ERROR - ORBIT EQN. CHECK
  real(fp), dimension(:), allocatable :: bphfshs
        ! 'Nt*m; fishbone loss'
        ! BEAM ION MOMENTUM TO FISHBONES BY SPECIES
  real(fp), dimension(:), allocatable :: bphjxbs
        ! 'Nt*m; JxB torque to plasma'
        ! BEAM ANG.MOMENTUM TO JXB TORQUE, BY SPECIES
  real(fp), dimension(:), allocatable :: bphjxbs_ntr
        ! 'Nt*m; JxB torque to plasma'
        ! BEAM ANG.MOMENTUM TO JXB TORQUE, BY SPECIES, from NUTRAV only
  real(fp), dimension(:), allocatable :: bphjxbs_dep
        ! 'Nt*m; JxB torque: deposition FLR step'
        ! BEAM ANG.MOMENTUM TO JXB TORQUE, BY SPECIES (deposition calc. only)
  real(fp), dimension(:), allocatable :: bphohbs
        ! 'Nt*m; torque to fast ions from OH circuit'
        ! BEAM ION MOMENTUM FROM OH CIRCUIT BY SPECIES
  real(fp), dimension(:), allocatable :: bpholds
        ! 'Nt*m*S; prev. timestep momentum'
        ! BEAM ION MOMENTUM BY SPECIES, PREV TIMESTEP
  real(fp), dimension(:), allocatable :: bphorbs
        ! 'Nt*m; orbit loss'
        ! BEAM MOMENTUM TO LIMITER, BY SPECIES
  real(fp), dimension(:), allocatable :: bphorbs_dep
        ! 'Nt*m; orbit loss (deposition FLR)'
        ! BEAM MOMENTUM TO LIMITER, BY SPECIES (deposition calc. only)
  real(fp), dimension(:), allocatable :: bphrcas
        ! 'Nt*m; recapture of CX fast neutrals'
        ! CX RECAPTURE OF MOMENTUM BY SPECIES
  real(fp), dimension(:), allocatable :: bphrpls
        ! 'Nt*m; ripple loss'
        ! MOMENTUM LOSS FROM RIPPLE BY SPECIES
  real(fp), dimension(:), allocatable :: bphsts
        ! 'Nt*m; d/dt(angular momentum)'
        ! RATE OF STORAGE OF FAST ION MOMENTUM BY SPECIES
  real(fp), dimension(:), allocatable :: bphths
        ! 'Nt*m; thermalization'
        ! BEAM ION MOMENTUM THERMALIZING BY SPECIES
  ! ... block: torque_work .........
    ! these quantities are powers associated with work done by
    ! the beam (or other fast specie) on a rotating target plasma.
    !   e.g. work done, (torque) x (target angular velocity)
  real(fp), dimension(:), allocatable :: bptcols
        ! 'W; collisional torque work done'
        ! BEAM POWER TO ROTATING PLASMA ENERGY, TORQUE WORK
  real(fp), dimension(:), allocatable :: bptdfbs
        ! 'W; Anom. diffusion JxB torque work'
        ! BEAM ANOM DIFF. POWER TO ROTATION BY SPECIES
  real(fp), dimension(:), allocatable :: bpthros
        ! 'W; rotational energy via thermalization'
        ! BEAM POWER THERM INTO PLASMA ROTATION BY SPECIES
  real(fp), dimension(:), allocatable :: bpthsfs
        ! 'W; thermalization source friction'
        ! BEAM POWER THERM ROTATIONAL SCE-FRICTION CORR.
  real(fp), dimension(:), allocatable :: bpjxbs
        ! 'W; JxB torque work'
        ! BEAM POWER TO JXB TORQUE WORK, BY SPECIES
  real(fp), dimension(:), allocatable :: bpjxbs_ntr
        ! 'W; JxB torque work'
        ! BEAM POWER TO JXB TORQUE WORK, BY SPECIES, from NUTRAV only
  real(fp), dimension(:), allocatable :: bpjxbs_dep
        ! 'W; JxB torque work in deposition step'
        ! BEAM POWER TO JXB TORQUE WORK, BY SPECIES (deposition calc. only)
  real(fp), dimension(:), allocatable :: bpthafs
        ! 'W; thermalization asymmetry correction'
        ! BEAM POWER THERM ROT. SCE ASSYMETRY CORRECTION BY FAST SPECIE
  ! ... block: particle_balance .........
    !
    !  This block contains scalar information related to fast ion particle
    !  balance.  Some quantities are available summed over fast ion species.
    !
    !  Various combinations of these quantities can be used to check the
    !  particle balance of various aspects of the fast ion model.
    !
    !  For example, some combinations used by TRANSP:
    !
    !      -- for fast specie "i":
    !
    !  General particle balance:
    !
    !    +sfdepsc(i)    deposition source
    !    +sfxrcsc(i)    charge exchange recapture source
    !    -bxfshs(i)    -(fishbone loss, if ad hoc fishbone model is active)
    !    -bxrpls(i)    -(ripple loss, if ad hoc ripple loss model active)
    !    -bxflws(i)    -(orbit loss; ion strikes wall or limiter)
    !    -bxflws_prompt(i) -(prompt orbit loss; ion strikes wall or limiter)
    !    -snbcxs(i)    -(charge exchange loss inside plasma)
    !    -snbcxos(i)   -(charge exchange loss, scrapeoff region)
    !    -thermalization  -(thermalization ion source, cf profile sorbths)
    !    -sfdtbmc(i)   -(d/dt(fast ion population in plasma))
    !     ------
    !     (residual)    (should be near zero but not to machine precision).
    !
    !  For fast ion  deposition:
    !    +xninjs(i)     number of particles injected
    !    -xninjea(i)    energy conservation adjustment
    !    -sfshine(i)    shine-through (beams only)
    !    -sfdepmc(i)    number of particles deposited
    !     -----
    !     sfdepba(i)    residual-- should be very close to zero
    !
    !  this won't be perfect, because the calculation is normalized to
    !  injected power, not current.  A related quantity is:
    !
    !     sfdeprr(i)    Monte Carlo Russian Roulette adjustment
    !                   (will affect overall particle balance; deposition
    !                   Russian Roulette is applied to make room for new
    !                   MC ions e.g. when the fast ion population is rising).
    !
    !  which is related to an option to compute extra deposition tracks
    !  beyond the number needed for new ions, e.g. to improve statistics
    !  on fast neutral density profiles.
    !
    !  Also for beam deposition:
    !
    !     sfdepiz(i)    deposition:  impact ionization on thermal plasma
    !     sfdepcx(i)    deposition:  charge exchange with thermal plasma
    !     sfdbbiz(i)    deposition:  impact ionization on fast ions
    !     sfdbbcx(i)    deposition:  charge exchange with fast ions
    !     ----------
    !     sfdepsc(i)    total deposition source integrated over neutral tracks
    !
    !  These quantities are integrated over neutral tracks; it may be worth
    !  comparing sfdepsc(i) to sfdepmc(i); the latter is the n/sec actually
    !  added to the slowing down Monte Carlo model population.  sfdepmc(i)
    !  is noisier:  each track is collapsed to a start point for a single
    !  orbiting ion; and at the end, all weights are renormalized to give
    !  the correct injected power (watts) from input data, causing a small
    !  random error in the modeled number of particles deposited.
    !
    !  There is a slowing-down particle balance that includes Monte Carlo
    !  adjustments and which should yield a residual within machine precision
    !  of zero:
    !
    !    +sfdepmc(i)   deposition source
    !    +sfdeprr(i)   Monte Carlo Russian Roulette adjustment
    !    -bxfshs(i)    -(fishbone loss, if ad hoc fishbone model is active)
    !    -bxrpls(i)    -(ripple loss, if ad hoc ripple loss model active)
    !    -bxflws(i)    -(orbit loss; ion strikes wall or limiter)
    !    -bxflws_prompt(i) -(prompt orbit loss; ion strikes wall or limiter)
    !    -sfcxesc(i)   -charge exchange fast neutrals escaping
    !    -sfcxrr(i)    -charge exchange Monte Carlo Russian Roulette adjustment
    !    -thermalization  -(thermalization ion source, cf profile sorbths)
    !    -sfdtbmc(i)   -(d/dt(ion population in plasma))
    !    -sforbrr(i)   -orbit Monte Carlo Russian Roulette adjustment
    !     --------
    !     sforbal(i)   residual -- should be zero to REAL*4 machine precision
    !
    !  Particle balance information is available on fast ion charge exchange
    !  and fast neutral recapture:
    !
    !    +sfcx0mc(i)   Fast neutrals created (n/sec, fast specie i)
    !    -sfcxesc(i)   Fast neutrals escaping
    !    -sfcxrr(i)    Monte Carlo Russian Roulette adjustment
    !    -sfcxrmc(i)   Fast neutrals recaptured
    !     ---------
    !     sfcxbal(i)   residual -- should be zero to REAL*4 precision
    !
    !  Charge exchange recapture processes are similar to deposition processes:
    !
    !     sfxrciz(i)    recapture:  impact ionization on thermal plasma
    !     sfxrccx(i)    recapture:  charge exchange with thermal plasma
    !     sfrbbiz(i)    recapture:  impact ionization on fast ions
    !     sfrbbcx(i)    recapture:  charge exchange with fast ions
    !     ---------
    !     sfxrcsc(i)   total recapture rate
    !
    !  These quantities are integrated over charge exchange neutral tracks;
    !  it may be worth comparing sfxrcsc(i) to sfcxrmc(i); the latter will
    !  be noisier, as it reflects the reduction of each neutral track to
    !  a properly weighted random decision to recapture, or not.
    !
    !  There is analogous data on the originating charge exchange process,
    !  from the point of view of the slowing down fast ions:
    !
    !     snbxw0(i)     charge exchange on "wall" (edge source) thermal neutrals
    !     snbxv0(i)     charge exchange on "halo" (vol. source) thermal neutrals
    !     snbxbb0(i)    charge exchange on original beam neutrals
    !     snbxbb1(i)    charge exchange on fast neutrals (from charge exchange).
    !     ---------
    !     snbxtot(i)    total charge exchange sink rate
    !
    !  Which can be compared sfcx0mc(i)-- the number of fast neutrals actually
    !  launched by the Monte Carlo model (this will generally be noisier).
  real(fp), dimension(:), allocatable :: bxflws
        ! '#/S; fast ions lost to wall'
        ! BEAM IONS/SEC TO LIMITER, BY SPECIES
  real(fp), dimension(:), allocatable :: bxflws_dep
        ! '#/S; deposition FLR loss'
        ! BEAM IONS/SEC TO LIMITER, BY SPECIES (deposition calc. only)
  real(fp), dimension(:), allocatable :: bxflws_prompt
        ! '#/S;  prompt losses of fast ions to wall'
        ! BEAM IONS/SEC TO LIMITER, BY SPECIES
  real(fp), dimension(:), allocatable :: bxfshs
        ! '#/S; fishbone loss'
        ! BEAM PTCLS LOST TO "FISH-BONES" BY SPECIES
  real(fp), dimension(:), allocatable :: bxrpls
        ! '#/S; TF ripple loss'
        ! PARTICLE LOSS FROM RIPPLE BY SPECIES
  real(fp), dimension(:), allocatable :: sbcxxs
        ! '#/S; lost fast neutrals born in plasma'
        ! FAST NEUTS BORN INSIDE PLASMA WHICH ESCAPE, BY SPECIES
        ! note some neutrals may be born outside; this is the difference
        ! between sbcxx* and sfcxesc*
  real(fp), dimension(:), allocatable :: sfcx0mc
        ! '#/S; Monte Carlo fast neutral source'
        ! MONTE CARLO FAST NEUTRALS LAUNCHED, N/SEC
  real(fp), dimension(:,:), allocatable :: sfcxmc1d
  INTEGER, dimension(:), allocatable :: io_sfcxmc1d
  CHARACTER*32, dimension(:), allocatable :: io_sfcxmc1d_name
        ! '#/S; Monte Carlo fast neutral source'
        ! MONTE CARLO FAST NEUTRALS LAUNCHED, N/SEC
  real(fp), dimension(:), allocatable :: sfcxbal
        ! '#/S; CX particle balance'
        ! M.C. CHARGE EXCHANGE PARTICLE BALANCE
  real(fp), dimension(:), allocatable :: sfcxesc
        ! '#/S; CX loss to walls from plasma & SOL'
        ! FAST NEUTRALS ESCAPING, N/SEC
  real(fp), dimension(:), allocatable :: sfcxrmc
        ! '#/S; CX fast neutrals recaptured'
        ! M.C. CX RECAPTURE N/SEC
  real(fp), dimension(:), allocatable :: sfcxrr
        ! '#/S; CX recapture Russian Roulette corr.'
        ! CX RECAPTURE RUSSIAN ROULETTE "LOSS", N/SEC
  real(fp), dimension(:), allocatable :: sfdbbcx
        ! '#/S; deposition by beam-beam CX'
        ! BEAM DEPOSITION BY BEAM-BEAM CHARGE EXCHANGE
  real(fp), dimension(:), allocatable :: sfdbbiz
        ! '#/S; deposition by beam-beam ionization'
        ! BEAM DEPOSITION BY BEAM-BEAM IONIZATION
  real(fp), dimension(:), allocatable :: sfdepba
        ! '#/S; deposition ptcl balance check'
        ! M.C. DEPOSITION PARTICLE BALANCE
  real(fp), dimension(:), allocatable :: sfdepcx
        ! '#/S; deposition by CX w/ thermal ions'
        ! FAST ION SCE -- BEAM DEP. BY TH. CX
  real(fp), dimension(:), allocatable :: sfdepiz
        ! '#/S; deposition, ionization by therm. ions'
        ! FAST ION SCE -- BEAM DEP. BY TH. IONIZATION
  real(fp), dimension(:), allocatable :: sfdepmc
        ! '#/S; total Monte Carlo particle deposition'
        ! MONTE CARLO DEPOSITION, N/SEC
  real(fp), dimension(:), allocatable :: sfdeprr
        ! '#/S; deposition Russian Roulette adjustment'
        ! DEPOSITION RUSSIAN ROULETTE "LOSS", N/SEC
  real(fp), dimension(:), allocatable :: sfdepsc
        ! '#/S; deposition fast ion source'
        ! FAST ION SCE -- BEAM DEPOSITION
  real(fp), dimension(:), allocatable :: sfdtbmc
        ! '#/S; d/dt(fast ion population)'
        ! D/DT(FAST ION POPULATION) MONTE CARLO EVALUATION
  real(fp), dimension(:), allocatable :: sforbal
        ! '#/S; fast ion ptcl balance'
        ! FAST ION ORBIT MONTE CARLO CODE PTCL BALANCE
  real(fp), dimension(:), allocatable :: sforbrr
        ! '#/S; orbit Russian Roulette adjustment'
        ! FAST ION ORBIT CODE RUSSIAN ROULETTE "LOSS" N/SEC
  real(fp), dimension(:), allocatable :: sfrbbcx
        ! '#/S; fast neutral recapture by beam-beam CX'
        ! RECAPTURE BY BEAM BEAM CHARGE EXCHANGE
  real(fp), dimension(:), allocatable :: sfrbbiz
        ! '#/S; fast neutral recapture by beam-beam ii'
        ! RECAPTURE BY BEAM BEAM IONIZATION
  real(fp), dimension(:), allocatable :: sfshine
        ! '#/S; beam shine through loss'
        ! BEAM SHINE THRU LOSS, N/SEC
  real(fp), dimension(:), allocatable :: sfxrccx
        ! '#/S; fast neutral recap. by therm. CX'
        ! RECAPTURE BY TH. CHARGE EXCHANGE
  real(fp), dimension(:), allocatable :: sfxrciz
        ! '#/S; fast neutral recap. by therm. ii'
        ! RECAPTURE BY TH. IONIZATION
  real(fp), dimension(:), allocatable :: sfxrcsc
        ! '#/S; fast neutral recapture (total)'
        ! FAST ION SCE -- CX RECAPTURE
  real(fp), dimension(:), allocatable :: snbcxos
        ! '#/S; CX fast neutrals created in SOL'
        ! BEAM IONS/SEC CX OUTSIDE PLASMA, BY SPECIES
  real(fp), dimension(:), allocatable :: snbcxs
        ! '#/S; CX fast neutrals created in plasma'
        ! BEAM IONS/SECOND CHARGE EXCHANGING-- INSIDE PLASMA
  real(fp), dimension(:), allocatable :: snbxbb0
        ! '#/S; 2nd gen. fast neutrals: CX w/beam'
        ! FAST ION CX: OFF NEWLY INJECTED BEAM NEUTRALS
        ! NB*N0*sigma*v
  real(fp), dimension(:), allocatable :: snbxbb1
        ! '#/S; 3rd+ generation fast neutrals'
        ! FAST ION CX: OFF CX FAST NEUTRALS
  real(fp), dimension(:), allocatable :: snbxtot
        ! '#/S; total fast neutral source'
        ! FAST ION SINK / FAST NEUTRAL SOURCE: CX
  real(fp), dimension(:,:), allocatable :: snbx1d
  INTEGER, dimension(:), allocatable :: io_snbx1d
  CHARACTER*32, dimension(:), allocatable :: io_snbx1d_name
        ! '#/S; total fast neutral source'
        ! FAST ION SINK / FAST NEUTRAL SOURCE: CX
  real(fp), dimension(:), allocatable :: snbxv0
        ! '#/S; neutral source: CX w/halo neutrals'
        ! FAST ION CX: OFF VOL. SCE (HALO) NEUTRALS
  real(fp), dimension(:), allocatable :: snbxw0
        ! '#/S; neutral source: CX w/wall neutrals'
        ! FAST ION CX: OFF WALL NEUTRALS
  real(fp), dimension(:), allocatable :: xysbcxs
        ! '#/S; d/dt(fast ions in SOL)'
        ! D/DT(T BEAM IONS OUTSIDE PLASMA) by specie
  real(fp), dimension(:), allocatable :: xninjs
        ! '#/S; total fast ions injected'
        ! n/sec beam neutrals injected
  real(fp), dimension(:), allocatable :: xninjea
        ! '#/S; injected ions: energy conservation adjustment'
        ! n/sec beam neutrals injected
  real(fp) :: sftota = zero
        ! '#/S; total fusion product source'
        ! GRAND TOTAL FUSION PRODUCT SOURCE
        ! (sum of xninjs(...) for fusion product species only)
 
end module nbi_com
