
********************************************************************
**  This help file, along with the CQL3D code, is copyrighted under
**  the GNU-GPL free-software license, as described below, and at
**  the beginning of the a_cqlp.f CQL3D source.
********************************************************************
**
**
**  TO RUN CQL3D SET NAMELIST AS DIRECTED BELOW
**  AND
**  NAME THIS FILE cqlinput
**  [Last update, Dec. 1, 2020, or later.]
**
**  Several namelist groups are defined: setup0, setup, sousetup, 
**                                 rfsetup, trsetup, frsetup,....
**
********************************************************************
**
**  Character data input is generally limited to character*8,
**  unless otherwise stated.
**
**************************************************************
**
**  This cqlinput_help file documents the namelist input variables.
**
**  A reference describing the general problem solved along 
**  with some illustrative results is:
**  "The CQL3D Fokker-Planck Code", R.W. Harvey and M.G. McCoy, 
**  appearing in Proceedings of the IAEA TCM on Advances in 
**  Simulation and Modeling of Thermonuclear Plasmas, Montreal, 
**  1992, p. 489-526, IAEA, Vienna (1993); available from U.S.
**  Dept. of Commerce, National Technical Information Service,
**  Document DE93002962.
**  See also http://www.compxco.com/cql3d.html .
**  The two dimensional (2-velocity) precedent to the code
**  is CQL: G.D. Kerbel and M.G. McCoy, Phys. Fluids 28, 3629 
**  (1985).
**  See also, file: code_overviews_cql3d_genray
**
**************************************************************
**
**   NAMELIST &setup0  SECTION
**
**   The namelist group "setup" occurs in the cqlinput file twice.
**   [Now deprecated: works, but use &setup0.  See below BH070414.]
**   The first instance of the "setup" [now setup0] section of the
**   cqlinput namelist file contains the following variables (which
**   should not be repeated in the second instance of "setup"):
**
**   mnemonic, ioutput, iuser, ibox, noplots,lnwidth,nmlstout,
**   special_calls,cqlpmod,lrz,lrzdiff,lrzmax,lrindx,
**   ls,lsmax,lsdiff,lsindx,nlwritf,nlrestrt,
**
**   BH070414:  Above use of "setup" as a "first setup namelist"
**              is depracated.   PLACE above variables in 
**              separated namelist &setup0.
**              [Modification for SWIM Integrated Plasma Simulator.]
**              Nov/2017: First setup namelist can now also be
**              be referred to as "fsetup" in the cqlinput namelist
**              file.  [An accomodation for a pathscale compiler limitation,
**              and for users who modified cql3d to accept fsetup.]
**
**   Descriptions follow immediately below.
**
*************************************************************
** 
**  mnemonic is the run designator...to help keep track of runs. 
**    for example: mnemonic="wy01"
**    mnemonic is dimensioned a48, 
**                       i.e., ok for up to 48 characters.
**    Used as prefix to output PS and netCDF files.
**    default: mnemonic="mnemonic"
**
**  ioutput(1)  [After 2020-10-21]
**    ioutput(1)=0  means the printout to screen is reduced to a minimum,
**                  leaving namelist printout, warning messages,
**                  and physics-related values used or computed by code.
**    ioutput(1)=1  The above, and more printout - 
**                  values of arrays for diagnostic purpose.
**    ioutput(1)=2 (or .ge.2) The above, and even more diagnostic printout.
**                  (default value ioutput(1)=0).
**  ioutput(2)=0 ! No usage at present [2020]
**    default: ioutput(1:2)=0  NOTE: Old cqlinput files could have 
**                                   ioutput.ne.0.
**
**  iuser ="name" of user or "number" of user (no present use).
**
**  ibox="box ***" where *** is user's box number (no present use).
**
**  Parameters are set in param.h giving maximum dimensions
**    of variables used in the code.  The code may be run
**    with dimensions which are less than the input parameters,
**    as specified below by namelist variables lrz, lrzmax,
**    ls, lsmax, (and others,) as specified below.
**
**    The parameter machinea specifies the size of integer words:
**       (=1, for 8 byte integers, =2, for 4 byte integers).
**
**  lrz is the number of radial flux surfaces in the computational
**    grid on which Fokker-Planck equation is solved.
**
**  lrzdiff="enabled", compute distribution functions on subset
**    of the flux surfaces considered. 
**    This system was setup to reduce computational time, when
**    only interested in few of the FP surfaces, but want to specify
**    the radial profile of background species on a more refined lrzmax
**    grid.  (default="disabled")
**
**  lrzmax is number of flux surfaces, including any not FP'd.
**    lrzmax=1, then this gives single flux surface runs with no input
**           of profiles
**    lrzmax.ge.4 uses input profiles of density, etc.
**    (lrzmax=2,3 to be avoided, and problems with some cubic splines.)
**    lrzmax is set =lrz, if lrzdiff.ne."enabled"
**
**  lrindx(1:lrz)=indices of flux surfaces where the FP'ing is
**    performed.
**    (if lrz=lrzmax, i.e., all flux surfaces are to be FP'ed, 
**    then it is unnecessary to specify lrzmax, lrindx, and lrzdiff).
**  
**  noplots="enabled1" inhibits all plots.
**    default: noplots="disabled" 
**    [="enabled" inhibits old unimplemented graflib plots, and was
**     necessary for code to run, but now all old graflib replaced.
**     so it should have same effect as "disabled", BH and YuP, 120122]
**    noplots="enabled1"  can be useful to load the code without 
**    a pgplot library.
**    The -Wl,-noinhibit-exec gnu loader (ld) option exists in gfortran
**    which enables loading of the code in the presence of unsatisfied
**    externals.  With noplots='enabled', the references to pgplot
**    routines in the code will be avoided.  Similar options exist
**    for several other compiler groups.
**
**  lnwidth=line-width attribute for plotting; affects lines, graph
**          markers, and text.  Specified in units of 0.005 in, and
**          must be an integer.  [default=3]
**          [Smaller value may work better for vector plots.]
**          Conversion from .ps to .pdf using ps2pdf improves the plot.
**
**
**  nmlstout="enabled", prints namelists to stdout
**           "trnscrib", transcribes the cqlinput nml file to stdout
**            else,   suppress this output.
**            default:  nmlstout="trnscrib"
**
**  special_calls="enabled", then code makes system calls to find
**                name of system and pwd, for output.   This is not
**                allowed on some systems.
**               ="external", then files uname_output and pwd_output
**                are used to printout.  These many be created with
**                an external script containing
**                uname -a > uname_output
**                pwd > pwd_output
**               ="disabled", then no system calls and no reference
**                to uname_output and pwd_output files.
**                This gets around lack of some special functions,
**                e.g., second(), with some compilers.
**               default="enabled"
**
**  nlwritf="enabled"  saves the distribution function at the end of 
**               the run (in file distrfunc).   default="disabled"
**         ="ncdfdist", don't write distn function into text file
**               distrfunc, since can use mnemonic.nc save of f.
**               This file also contains much additional data, some
**               of which is needed with some restart options.
**               Text distrfunc will contain only namelist and profiles
**               at the end of the run.
**               Don't use netcdfshort.ne."disabled", if going to 
**               use nlrestrt.
**  nlrestrt="enabled" reads the initial distribution function in 
**               text file 'distrfunc'.           default="disabled"
**               The general restart idea is that only data from end
**               of the previous run which is not specified in the
**               continuing namelist file, needs to be restored.
**          ="ncdfdist" reads the distribution function from a
**               netcdf file, 'distrfunc.nc'.  This file can be
**               the mnemonic.nc file from a previous run.
**               The text file 'distrfunc' is still needed  (NOT as of 110201), 
**               but only for reading namelist from &setup and others (except %setup0),
**               ensuring compatibility of grids, etc.  (Can change
**               mnemonic in &setup0.)
**               Reading of text distrfunc for nlrestrt='ncdfdist' or
**               'ncregrid' discontinued (BH110201).
**               BH181023: If ampfmod="enabled" then also read in prior 
**               ending value of elecfld().
**               BH+YuP191029: If iprozeff='curr_fit', then also restore 
**               prior ending value of the zeff profile.
**          ='ncregrid', same as 'ncdfdist', except the momentum-per-
**               mass grid can be reset from that used in distrfunc.nc,
**               by changing enorm, jx and xfac (only enabled changes
**               presently, not iy or lrz) in cqlinput.  If enorm/vnorm
**               is increased, the distrn for distrfunc.nc is extrapolated
**               linearly in ln(f) versus v(mom-per-mass).  Portions of the 
**               distn fnctn in distrfunc.nc within the bounds of the
**               new x() grid specified by the namelist input are 
**               interpolated onto the new grid. The FSA density of the
**               new code distribution is maintained equal to the
**               restart distribution function in the .nc file.
**               enorm/vnorm can be either increased or decreased from
**               the value in the distrfunc.nc file.
**               Reading of distrfunc text file for nlrestrt='ncdfdist' or
**               'ncregrid' discontinued (BH110201).
**               BH181023: If ampfmod="enabled" then also read in prior 
**               ending value of elecfld().
**           NB: The nlrestrt option is presently (181023) setup only
**               for a single general species (ngen=1).
**
************************************************************************
**  BH070414: Variables nlwritf and nlrestrt have been changed from   **
**            logical to character*8. Check if old NL needs adjusting.**
************************************************************************
**
**  This namelist section also contains:
**  cqlpmod,ls,lsmax,lsdiff,lsindx,nlwritf,nlrestrt,
**    lrz=1,lrzdiff,lrindx,  used for CQLP run (see below).
**
**
**    
**************************************************************
**
**   NAMELIST &setup  SECTION
**
**************************************************************
**
**   MAJOR TIME STEP CONTROLS
**   (The following variables down to the "sousetup" namelist
**     description , are in "setup" namelist.)
**
*************************************************************
**
**  dtr is the time step in seconds. 
**    (n is time step counter.  Initially n=0.  It is incremented after
**     each step).
**    dtr can be nonphysically large for implicit steady state
**    calculations, with transp="disabled": 
**    dtr= 1. --> 10. is a typical choice if implct="enabled".
**    This allows code to "jump" to steady state (for linear problems).
**    transp="enabled"/soln_method='direct' cases will require dtr 
**    values less than or of order the collision time of the nonthermal
**     particles, to get accurate results, to prevent oscillations
**     between the velocity and radial space solutions.
**    default: dtr=5.
**
**  dtr1(1:10) is array of  additional time step settings 
**    (can be shorter or longer than dtr) starting
**    at time step nondtr1(1:10), if nondtr1().gt.-1. (default: dtr1()=0.0)
**
**  nondtr1(1:10) is array of steps at which the corresponding new
**    time step dtr1(1:10) starts.
**    Choose mod(nondtr1().gt.0,nrstrt)=0, else it is reset to next
**    highest value.  (default: nondtr1(1:10)=-1).
**
**  nstop is the number of time steps (or cycles) to be run.
**    default: nstop=5
**
**  nrstrt is the number of cycles to advance on a flux
**    surface before moving on to the next flux surface.
**    Operant only if lrz > 1
**    default: nrstrt=1
**    BH100517:  Defunct functionality.  Only possibility is
**               default case, nrstrt=1.
**
**  nplot(1:nplota) gives time steps in increasing order for plotting 
**    of data for individual flux surfaces.  nplot(1)=0 is OK.
**    Also used for netcdfshort='lngshrtf' for output of distributions f.
**    default: nplot(1:nplota)=-10000
**
**  nsave(1:nsavea) For netcdfshort='lngshrtf', gives time steps in 
**    increasing order for saving the ngen particle distributions
**    to the .nc output file.  nsave(1)=0 is OK.
**    default: nsave(1:nsavea)=-10000
**
**  iactst="enabled" means certain debugging diagnostics are called:
**     subroutines diagimpd, exsweepx, and exsweept are utilized.
**     This option should be used only in debugging mode. 
**     iactst="abort" means abort the code if the computed 
**     error is .gt. 1.e-8 (relative difference between integrated
**     rhs and lhs of difference equation).
**     default: "disabled"
**
**************************************************************
**
**  MESH CONTROLS
**
**************************************************************
**
**  jx is the number speed mesh points. 
**    Speed also means momentum/rest mass for relativistic calcs.
**    default: jx=300  (in aindflt.f).
**    The basic speed mesh runs from 0. to xmax=1.0.
**
**  iy is the number of theta mesh points. Must be even number.
**    default: iy=200  (in aindflt.f).
**    The y (theta) mesh is symmetric about pi/2.
**    For urfmod.ne."disabled", iy.gt.255 not allowed with 
**      the storage scheme for ilim1,ilim2,....[BH070418: relaxed now?]
**    (Have had difficulty constructing a mesh iy as small as 30, 
**     and meshy.eq."fixed_y".)
**    The positive theta dirn (i.e., v_parallel) has positive component
**       in the positive toroidal direction (counter-clockwise, looking
**       from the top).
**    In some transport settings (transp=.ne."disabled"), the number
**       of theta (y) mesh points may reduce as a function of radius.
**       The y points can also be stacked up around the trapped-passing
**       boundary.   Thus, the y array is two dimensional, y(pitch,rho).
**
**  lz is the number of poloidal (or field line mesh points). Bounce-
**    averages are computed over this mesh. 
**    default: lz=40 (If cqlpmod.eq."enabled", lz is set equal to lsmax)
**
**  vnorm is the velocity normalization constant (cm/sec).
**    For relativistic calculations, vnorm is the maximum momentum/unit
**    rest mass on the mesh.
**    default: vnorm is set through enorm (below).
**    Note that in plots, the variable vnorm is printed as "unorm",
**    to be consistent with meaning of "u" as 
**    the relativ.moment/rest.mass.
**
**  enorm is the energy corresponding to the velocity normalization
**    constant; if kenorm corresponds to a species index, that
**    species determines vnorm through enorm, superseding vnorm. (kev).
**    For relativistic electrons, E=(gamma-1)*m*c**2=(gamma-1)*511. keV.
**                                gamma**2=1+(u/c)**2
**                                u=momentum-per-electron-rest-mass.
**    default: enorm=200., kenorm=1
**
**  tandem="enabled" facilitates the use of both ions and electrons
**    as general species. The mesh extends to a speed (momentum)
**    represented by an electron energy of enorme. However, the
**    majority of the mesh points are dedicated to a region of
**    velocity space which adequately represents ions with a
**    energy less than enormi. The ion species selected for this
**    treatment is the first general ion species listed. Most
**    of the plots are fixed so that the ion contributions are visible
**    and not squeezed into the first 1% of the plot.
**    tandem="disabled" disengages this option.
**    default: tandem="disabled"
**
**  enorme,enormi are analogues to enorm for tandem="enabled"
**    enorme corresponds to the top of the electron mesh,
**    enormi corresponds to the top of the ion mesh.
**    [defaults: enorme=enorm, enormi=enorm]
**
**  relativ="enabled" means the quasi-relativistic calculation and
**    mesh are activated.  Best for fusion plasmas.
**         ="fully" give fully relativistic calc.   Has some limitations
**           in maximum order Legendre polynomial to be used in
**           expansion for FP coefficients (mx, below, .le.3, or so).
**           See CompX report CompX-2009-1_Fully-Rel.pdf.
**    default: relativ="enabled"
**
**  xfac is the velocity mesh spacing parameter.
**    xfac=1. ==> uniform mesh;
**    xfac<1. ==> geometric mesh with greater mesh resolution at x=0;
**    xfac>1. ==> greater resolution at xmax.
**    xfac<0. ==> xpctlwr,xpctmdl,xlwr,xmdl must be input.
**    xfac<0. ==> the momentum (velocity) mesh contains three regions:
**    The region [0,xlwr] will have jx*xpctlwr evenly spaced 
**    mesh points.
**    The region [xlwr,xmdl] will have xpctmdl*jx mesh points.
**    The region [xmdl,xmax] will have the balance of the jx 
**    mesh points.
**    default: xfac=1.
**    (xfac,xpctlwr,xlwr,xpctmdl,xmdl are reset appropriately, if 
**     tandem="enabled").
**
**  tfac is the theta mesh spacing parameter.
**    tfac=1 ==> uniform mesh except in the vicinity of the
**      pass/trapped boundary region. See tbnd below.
**    tfac<1. ==> geometric mesh with points packed closer at p/t bndry.
**    tfac>1. ==> geometric mesh with points packed closer 0., pi, 
**    and pi/2.
**    default: tfac=1.
**  tfac=-1.  ==>Appropriate for radial transport, transp="enabled"
**    calculations.   This will give uniform theta mesh at each 
**    radius.   There will be slight adjustments near t-p bndry
**    at each radius, but no addition of extra points bounding the
**    trapped-passing boundary, as for tfac.gt.0. case.
**    If tfac.gt.0. and transp="enabled" and soln_method.ne."direct"
**    then tfac is reset to -1. at beginning of run.
**    
**    If tfac negative and .ne.-1, then mesh adjusted according to 
**    positive tfac above, with abs(tfac) values.  No added pts near t-p.
**  There is also a tfac.lt.0./cqlpmod="enabled" option (see below).
**
**  yreset="enabled" means pack extra theta mesh points in some
**    specified (below) region.
**    default:yreset="disabled"
**    [Not working at present, and if "enabled" will be reset to
**     "disabled" in ainsetva.f, for lrz>1.] 
**
**  numby is the number of packed theta mesh points
**    numby should be significantly less than iy/2 (or resolution
**    will be to coarse in the rest of pitch angle space).
**    Recall that the y (theta) mesh is symmetric about pi/2.
**    Only used if yreset="enabled".
**    default: numby=20
**
**  ylower and yupper specify the radian spread of the packed region
**    below pi/2. Avoid allowing [ylower,yupper] to include the p/t
**    boundary region [y(itl-1),y(itl+1)].
**    This option is automatically disabled if lrz > 1.
**    Only used if yreset="enabled".
**
**  tfacz is the field line coordinate (z) mesh spacing parameter.
**    default: tfacz=1.
**    tfacz is not operational for eqsym="none", in which case the z-mesh
**    is uniform on either side of the max-|B| point on the flux surface.
**
**  tbnd(1:lrorsa) is the width (in radians) of the passing/trapped boundary 
**    layer. There is usually no reason to change the default value.
**    If necessary the program may reset this value at some
**    flux surfaces.
**    default: tbnd(1)=.002, tbnd(2:lrorsa)=0.
**
**  rfacz is the radial mesh spacing parameter. The minimum value of
**    the mesh can be pre-specified (see roveram below). Otherwise
**    rfacz works like xfac for xfac > 0.
**    default:rfacz=.7
**
**  roveram determines the value of r/radmin = rya(1) by setting
**    rya(1)=roveram.   It works in conjunction with
**    rfacz above. If roveram < 1.e-8, the minimum flux surface is
**    determined internally with no user control. This option is
**    useful for keeping the innermost flux surface away from the
**    magnetic axis.
**    default: roveram=1.e-8
**
**  rzset="enabled" allows the user to specify the entire radial
**    rya mesh as an array of values r(l)/radmin where r(l) is the
**    radial coordinate and radmin in this case takes on
**    the meaning of the radial coodinate at the LCFS.
**    ="disabled", mesh follows rfacz directive, above.
**    default: rzset="disabled"
**
**  rya(l),l=1,lrzmax is the normalized radial mesh referred to above in
**    the rzset description. The user specifies rya(1)=.02, rya(2)=.05,...
**    rya(23)=.95, for instance, if lrzmax=23
**    Thus rya(l)=r(l)/radmin, and rya is a code array.
**    (Note: rya is dimensioned rya(0:lrza+1), so need to specify
**    first namelist input value explicitly as rya(1)=....).
**
**  radcoord="sqtorflx" (default) use normalized radial coordinate (that is,
**           flux surface label) proportional to sqrt of toroidal flux.
**           [Presently, radial diffusion eqn is only set up for this value
**            of radcoord.  Additional choices below are for no radial 
**            transport.]
**          ="sqarea" square root of flux surface area,
**          ="sqvol" square root of flux surface volume,
**          ="rminmax" 0.5*(max. major radius - 
**                          min. major radius) of flux surface
**          "polflx"  normalized poloidal flux from magnetic
**                    axis to last closed flux surface (LCFS).
**          "sqpolflx"  normalized square root of poloidal flux function.
**
*********************************************************
**
**  DIFFERENCING and SOLUTION METHOD CONTROLS
**
*************************************************************
**
**  chang= "enabled" forces Chang Cooper differencing; = "disabled" forces
**    standard centered differencing. chang="noneg" employs an additional
**    scheme to minimize resulting negative distributions in strongly
**    driven RF fields. Actually, we employ a variant of strict Chang-
**    Cooper differencing.  Typical Chang-
**    Cooper looks at dx*A/B where A is the advective coefficient and
**    B the diffusive. The closer the ratio is to zero, the closer the
**    differencing is to central. As the ratio grows larger in absolute
**    value, the differencing shifts to upwind. We have
**    empirically determined that in the presence of strong RF excitation
**    that even with the use of Chang-Cooper, the distribution can
**    be driven negative. This may be due to a 2-D effect, since the
**    RF diffusion characteristics rarely coincide with the theta and v
**    mesh directions. It was discovered that by subtracting off the
**    contribution from the RF diffusion from the total diffusion
**    coefficient and then computing the ratio above, the scheme is
**    driven in the upwind direction and the problem is vastly
**    ameliorated, in most cases, the difficulties disappear. The
**    ratio is now dx*A*B/(B-B_ql)**2 for the
**    velocity differencing. In the limit B_ql==>0, we recover
**    Chang-Cooper. For the theta differencing, the approach
**    is similar, except that if chang="noneg", strict one sided
**    differencing is enforced in the regions where the distribution has been
**    previously driven negative (during earlier time-steps). All
**    of this manipulation is a result of experimentation and
**    has not been proven to be useful rigorously....as indeed
**    neither has Chang Cooper for the 2D problem.
**    default: chang="enabled"
**
**  ineg="enabled" means negative values of the distribution function 
**    are reset to zero after time advancement; ="disabled" 
**    defeats this option.  But, see update below 03-2011 comment.
**    default: ineg="enabled" [YuP: actually, it is "disabled"]
**   ="renorm" means if f developes negative values after a time
**    step, then add 1.1*abs(min(f)) to f, and then multiply by
**    (old_line_density/new_line_density), to keep total number
**    of particles constant during this renormalization.  This option
**    is an attempt to deal with negative distributions sometimes
**    occurring with QL diffusion.  It can add unacceptable amounts of
**    high energy particles to the distribution. (BobH, sept/94)
**   ="trunc_d" is another attempt to cope with QL diffusion causing
**    negative distributions:  the QL diffusion coefficients are 
**    smoothly reduced to zero over the last 10 velocity grid points. 
**    Also, negative values of f reset to zero.
**    
**    Modification made after [03-2011] to ineg="enabled" or "trunc_d": 
**    (See diagimpd.f)
**    If distr.function f(i,j) is negative or zero 
**    at some (i,j)-point in vel. space,
**    it is set to zero for ALL jj.ge.j 
**    (and fixed i-index of pitch angle),
**    for energies greater than the maximum in the source term
**    (for example from NBI).  For lower energies, neg f is set = 0.
**    Before this modification, it was set to zero 
**    for only this specific (i,j)-point where f(i,j)<0.
**    So, no more "islands" in distribution function 
**    remaining beyond such (i,j)-points, except when they are 
**    caused by sources like NBI.
**   ="enabled1" as in ineg="enabled" [after 03-2011], except
**    f is zeroed at each pitch angle i for all j above the
**    first velocity where f(i,j)/f(1,1).lt.1.e-20.  This removes
**    spikes in the tail distribution resulting from (NBI) sources
**    injected within prompt-loss regions (e.g., lossmode(1)='simplban'
**    which may create problems for nonthermal full-wave code calculation.
**    For other purposes, the user might want to see these spikes, which
**    represent injected density for one time step, before removal.
**
**  implct= "enabled" forces implicit differencing and
**    Gaussian elimination. implct  = "disabled" forces operator 
**    splitting (can be used for time dependent calculations).
**    The implct namelist variable relates to (only) the velocity-space
**    differencing of the FP eqn, and method of solution.
**    implct="enabled", uses LAPACK linear algebra subroutines which have
**    been incorporated into the source code.
**    default: implct="enabled"
**
**  manymat="enabled" allows all the factored matrices on all flux
**    surfaces to live on disk and not disappear after the backsolve.
**    This means that the next time step will NOT require a new 
**    factorization IF THE RHS HAS NOT CHANGED. This saves lots
**    of computer time and is the recommend operating mode. 
**    manymat="disabled" does not save the factored matrices on disk
**    and is useful for single flux surface calculations where the
**    core space used for the factorization is not overwritten
**    as the code moves from flux surface to flux surface.
**      default: manymat="disabled".   Only implemented for
**      soln_method="direct"
**
**  soln_method="direct"[default],  use lapack dgbtrf/dgbtrs direct
**               LU factorization and solve of the velocity-space
**               equations.
**               v-space coeff matrix A(,) is inefficiently stored
**               in LAPACK storage system (with many zeroes).
**               In radial transport case (transp.eq."enabled") will
**               also retain the splitting method, for this setting.
**             ="itsol", obtains v-space coeffs in A(,) directly.
**               Uses SPARSKIT2 iterative GMRES solver with
**               ilut drop-tolerance preconditioner for v-space solve,
**               In radial transport case (transp.eq."enabled") will
**               also retain the splitting method, for this setting.
**             ="itsol1",uses translation of inefficient lapack
**               coeff storage of v-space coeff matrix A(,)to CSR using
**               SPARSKIT2 translation routine bndcsr.  Then solves as
**               in "itsol" case.  This option is just for cross-check
**               that CSR and LAPACK storage methods are understood.
**             ="it3dv", full 3D v-space matrix (but not
**               off-diagonal terms for transport yet) with SPARSKIT2
**               soln of full matrix simultaneously for all flux surfaces.
**               Gives better converged solution then for 
**               soln_method="itsol".  [I don't think this is compatible
**               with transp="enabled" and splitting, BH080319].
**               [Now OK, 091203].
**             ="it3drv", is full 3D-implicit solve including radial
**               transport terms, using SPARSKIT2 GMRES routine with
**               preconditioning and drop-tolerance.
**               
**  droptol=drop tolerance for the itsol preconditioner,
**          1.d-3 [default].  Off-diagonal elements smaller than
**          this factor times diagonal element will be eliminated.
**  lfil= Maximum fill-in of L and U.  If .ge. number of equations to be
**        solved, it will have no effect. See following. [default:lfil=30]
**  From itsol.f preconditioner comments:
**c---- Dual drop strategy works as follows.                             *
**c                                                                      *
**c     1) Thresholding in L and U as set by droptol. Any element whose *
**c        magnitude is less than some tolerance (relative to the abs    *
**c        value of diagonal element in u) is dropped.                   *
**c                                                                      *
**c     2) Keeping only the largest lfil elements in the i-th row of L   * 
**c        and the largest lfil elements in the i-th row of U (excluding *
**c        diagonal elements).                                           *
**c                                                                      *
**c Flexibility: one  can use  droptol=0  to get  a strategy  based on   *
**c keeping  the largest  elements in  each row  of L  and U.   Taking   *
**c droptol .ne.  0 but lfil=n will give  the usual threshold strategy   *
**c (however, fill-in is then unpredictible).                            *

**
**  lbdry(k)="conserv" means conservation boundary conditions used
**    at v=0 for general species "k" (see species descriptors).
**    Zero v-flux at v=0 is enforced in this case.
**    Also, no resetting of plasma density, so particles can
**    flux off the upper grid edge, be lost due to loss terms,
**    or be gained due to sources.
**  lbdry(k)="fixed" means zero flux boundary conditions are not 
**    enforced at v=0. Instead the distribution is held constant at 
**    v=0 during "time" advancement.  Also, if elecfld is set through 
**    eoved, then f(v=0) is only set at time step n=0, and not reset.
**    Unicity at v=0 is applied.
**  lbdry(k)="scale" is like "fixed", but at each time step after solution
**    for the distribution, f, is found, the entire distribution is scaled
**    so that the resulting field line integrated number density does
**    not change in "time". Unicity at v=0 is applied.
**    As of June 30, 2017: lbdry(k)="scale" is reset in the code
**    to "consscal".  This fixes a recent (Apr/2015) bug in cql3d, and
**    appears to be a more clear way of scaling the distribution to maintain 
**    given density.
**  lbdry(k)="consscal", then combine above "conserv" and "scale".
**            For time-dependent background plasma, nbctime.gt.0,
**            relativistic Maxwellian distributed particles at the 
**            local temperature are added.
**            For nbctime=0 case, then the whole f, the distn, is scaled.
**  lbdry(k)="conscalm", Combine above "conserv" and scale density
**                       with added (relativistic) Maxwl distrn at namelist
**                       (relativistic) temperature, including the case
**                       of nbctime=0.
**  lbdry(k)="conserv" is the ONLY option for implct="disabled"
**    default: lbdry()="conserv"  
**  NOTE:  no "enabled" or "disabled" option for lbdry().
**
**  Note: It is not always a good idea to rescale the distr.func.
**        Example: When there is a large-power NBI source, 
**        comparing to the initial background 
**        set in cqlinput. The value of xlndn00 
**        (in ratio(k,lr_)=xlndn00(k,lr_)/runden, 
**        which is the rescaling factor) is based 
**        on the initial density, so it does not include particles from NBI.
**        And the value of runden (or sden) does include all sources.
**        In such a case, it may be better to use lbdry(k)="conserv"
**
**  Following lbdry0 input only applies for implct="enabled":
**  lbdry0="enabled" means use new system computing distrn function f
**    at v=0 that takes f to be a single value.  This uses 
**    subroutine impavnc0 for implicit time advancement.  (Modifications
**    originated by Olivier Sauter). (default="enabled", implct="enabled").
**    This option only applied if lbdry(i)="conserv" or "conserv_scale".
**  lbdry0=.ne."enabled" means use old system which computed non-unique
**    f(v=0) and then set f(v=0) to density conserving average.  Uses
**    subroutine impavnc. (Applies for implct="enabled").
**
**  taunew="enabled" means use new fully numerical calculation of dtau
**          and tau times (BobH, 970210).
**          Although more accurate, there remains some work to ensure
**          compatibility of the extended theta_poloidal region of
**          dtau with the sinz,cosz,tanz determined for a single
**          orbit for each velocity grid point (BobH, 010320).
**          Oct'09: This option is the only possibility with eqsym.eq."none".
**                  Have removed extended theta_poloidal region of
**                  dtau by only integrating down the center of the
**                  pitch angle bin (nii=1); the "extended" feature is
**                  confusing [to me, BH] and doesn't add much accuracy.
**         ="disabled" gives old partially analytic method described
**          in Killeen et al. book.  (default="disabled", changed 010320).
**  nstps=related to integration for tau and dtau. The number of
**        subintervals into which each dz along the B-field is
**        divided, for int[dz/|v_parallel|]. default:nstps=100
**
********************************************************
**
**  MACHINE GEOMETRY
**
**************************************************************
**
**  machine="toroidal" or machine="mirror"
**    Use of "mirror" makes available Gary Smith's B field option
**    (see below). Use of "toroidal" allows a variety of magnetic
**    field options (see below).
**    machine="mirror" is disabled for lrz .gt. 1
**    default: machine="toroidal"
**
**  psimodel determines the model used (or the method used) to determine
**    the magnetic field variation in the poloidal direction, that
**    is psi(z)=B(z)/B(0).
**  psimodel="axitorus" means standard 1/R variation for circular cross
**    sections is used. (Only slightly sensible for eqmod.ne."enabled";
**    even with eqmod.ne."enabled", I recommend psimodel="spline"BH000416)
**  psimodel="smith" means Gary Smith's model is used for mirror problems.
**  psimodel="spline" means that if through the use of some other model,
**    we have available on some arbitrary z (field line) grid an array
**    psi(z), then we utilize splines to determine psi(z) on the
**    computational z mesh. This option is generally used only when
**    the "eq"uilibrium module is enabled and an MHD "eqdsk" file is
**    available.  (Also works with eqmod.ne."enabled",machine="toroidal".)
**    default: psimodel="axitorus",
**             But psimodel reset to "spline" if eqmod.NE."disabled".
**
**  Set gsla and gslb if machine="mirror" and psimodel="smith"
**    gsla and gslb are field scale lengths such that:
**    psi(s)=[1+(s/gsla)^2
**              +gsb*(exp(-(s-gszb)^2/gslb^2)
**              +gsb*(exp(-(s+gszb)^2/gslb^2)]/N
**    Sets gsb and gszb so psi'(zmax)=0, psi(zmax)=rmirror
**    N is chosen so psi(0)=1 and psi'(0)=0 by construction
**    default: gsla=270. gslb=35.
**
**  ephicc is the throat potential jump in kV for machine="mirror"
**    This is used to determine the loss cone.
**    default:ephicc=0.
**
**  zmax(1) is the arclength along B from midplane to throat in 
**    machine="mirror"
**    YuP: Actually, zmax is found in eqfndpsi/eqorbit 
**         as zmax(lr_)=zmax_=es_(lorbit_) (last point of field line)
**  rmirror is the mirror ratio for machine="mirror"
**    rmirror=2.
**
**  radmin is the torus minor radius (cm) for machine="toroidal"
**    default: radmin=50.   Not required, for eqmod='enabled'.
**
**  rovera(1) (for rzset.ne."enabled") is the (normalized) minor radius
**     of the toroidal flux surface
**    (r/a=r/radmin) which also serves as the toroidal drift surface
**    in the zero banana width approximation for machine="toroidal"
**    For normal operation 1.e-8 < rovera < 1.
**    If rovera is positive and less that 1.e-8, it will be reset to 1.e-8.
**    Note the rovera input option is used only for cases where lrz=1, and
**    rzset.ne."enabled".
**  rovera(1) < 0 can be utilized only if eqmod="enabled"
**    In this case an input variable povdelp is utilized to determine
**    the flux surface of interest. povdelp is described in the "eq"
**    model input below.
**    default: rovera=.1
**
**  radmaj is the torus major radius (cm) for machine="toroidal",
**    for the case that eqmod="disabled" and psimodel="axitorus"  
**    default: radmaj=100.     Not required, for eqmod='enabled'.
**
**  btor is B(toroidal) (g) at the magnetic axis for machine="toroidal"
**    for the case that eqmod="disabled" and psimodel="axitorus"  
**    default: btor=1.e+4     Not required, for eqmod='enabled'.
**
**  bth is B(poloidal) (g) at the plasma edge for machine="toroidal"
**    for the case that eqmod="disabled" and psimodel="axitorus"  
**    The poloidal coordinate at which bth is defined is pi/2.
**    default: bth=1.e+3     Not required, for eqmod='enabled'.
**
***********************************************************
**
**  SPECIES DESCRIPTORS
**
**************************************************************
**
**  In the following, the index (k) refers to the species being
**    assigned the particular value.  The following conventions are
**    assumed:
**
**    If 1 .le. k .le. ngen, then k is a general (time-advanced) species.
**    k.gt.ngen means species k is background Maxwellian.
**    General species are represented by f(theta,x,k,l_), a distribution
**    function which evolves in time, whereas the background Maxwellians
**    are represented by two quantities, temp(k,) and reden(k,).
**    In k order, the general (non-Maxwellian time advanced species) 
**    come before the background species. Electrons and
**    ions can be mixed, and any species can be represented
**    simultaneously as a general and a background species.
**    If there is more than one Maxwellian ion, then they must
**    not be separated in k by the electron species.   If more than
**    one of the ion species has the same charge number bnumb(), (e.g.,
**    H+, D+, and/or T+) then they should be placed successively at the
**    beginning of the Maxwellian ion species.
**
**********************************************************************
**
**  ngen is the number of general species. It must not exceed
**    parameter value ngena or code will abort with error message.
**    default:ngen=ngena
**
**  nmax is the number of Maxwellian species. It must not exceed
**    parameter nmaxa or the code will abort with an error message.
**    default: nmax=nmaxa
**
**  kspeci(1,k) is the name of species(k) for k=1,ntotal(=nmax+ngen).
**  kspeci(2,k) is "general" or "maxwell" for k=1,ntotal, 
**    designate species as general of Maxwellian.	
**
**  fmass(k) is the mass(k) in grams.
**    default: fmass(k)=1.e-29
**    Ions are defined by having mass .gt.1.e-27.
**
**  bnumb(k) is the atomic number  (-1. for electrons)
**    default: bnumb(k)=1.
**
**  iprote, iproti and iprone, are flags with several possible values
**     (below) governing the electron temperature, ion temperature
**     and densities of the plasma profiles for the species.
**     We have the following possibilities:
**  iprote="parabola","spline", "asdex" "prbola-t", and "spline-t"
**     (similarly for iproti and iprone).
**     Each option is described below.
**     (The "-t" are time-dependent, and are descibed in a separate
**     TIME-DEPENDENT PROFILES section below.)
**     default: "parabola" for everything
**  iprozeff, iproelec, iprovphi, ipronn:  See below.
**
**  iprote or iproti or iprone = "parabola"      
**    npwr(k) and mpwr(k) (real) specify the profiles: 
**     The  form is (center-edge)*(1-(r/a)**npwr)**mpwr+edge between the
**     specified center and edge values described below (see subr profaxis).
**  npwr(0) and mpwr(0) give reden profiles whereas npwr(k) and
**     mpwr(k) give the temp profiles for species k. npwr and mpwr
**     are real variables, not integers.
**     default: npwr(0)=2.; npwr(k)=2. for all (k=1:ngen+nmax) (k=1:ntotal).
**     default: mpwr(0)=1.; mpwr(k)=3. for all k>0.
**
**  reden(k,l) is the initial density at the midplane (particles/cm**3).
**     FOR lrz.gt.1 (and iprone.eq."parabola"): 
**     The user specifies reden(k,0) (central value) and reden(k,1)
**     (edge value). The intermediate values (l=1,lrzmax) are computed 
**     by the code (see below).  
**     FOR lrz.eq.1:
**     The user specifies reden(k,1)
**     This rule (involving lrz) is used throughout input. See temp
**     and eparc and eperc for example. 
**     default: reden(k,l))=1. (for l=1,lrzmax)
**
**  temp(k,l) is the initial temperature (Kev) of species k. This
**    applies both to general and background species. See reden above
**    for specification rule.
**    But, for time-dependent profiles, temp() can evolve in time,
**    see profiles.f.
**    default: temp(k,l)=1.
**
**  profpsi="disabled" means that the density and temperature profiles
**    will obey the npwr and mpwr directives (above).
**  profpsi="enabled" means that the profiles have the following
**    dependence (used in early days of the code):
**    reden(k,l)=reden(k,0)*((psi(l)-psimag)/psimag)**.5
**    temp(k,l)=temp(k,0)*(psi(l)-psimag)/psimag
**    default: profpsi="disabled"
**
**  iprote, iproti, iprone, iprozeff, iproelec, iprovphi, ipronn = "spline"  
**      In this case the user can specify namelist variable
**      arrays which give radial profile information for cubic splines. 
**      A spline routine will interpolate the ene, te, and/or ti data 
**      onto the computational mesh.
**      For specified profiles of ene, te, and/or ti, we have:
**      (l=1:100,k=1,ntotala=ngena+nmaxa)
**      See iprozeff below.
**    ryain(l)  (r/a)  l=1,njene (the radial coordinate)
**      (For ampfar problem spline elecfld, ryain(1) assumed = 0.)
**    enein(l,k)  the specified density of species k at ryain(l) (/cm**3)
**      (if enein(1,k)=0. for species k.gt.1, then the enein(l,k-1)/bnumb(k)
**       values will be used for species k). 
**    tein(l)   the specified electron temperature (keV)
**    tiin(l)   the specified ion temperature (same for all ion species
**                       such that bnumb(k).ne.-1.) (keV)
**    enescal, tescal, tiscal (and zeffscal, elecscal, see below)  
**    are scale factors for "spline" input.
**       (defaults = 1.0).
**    20100126:  Also added that enescal,tescal,tiscal multiply
**               all input density/temperature profiles.
**               Can be useful, e.g., when units need adjusting.
**               elecscal multiplies t-dep input for efswtch.eq."method1".
** 
**    njene must be less than njenea (= 256 is typ value) 
**       default:  0 
**
**   (There are also namelist input variables njte, njti (but no
**    njzeff) similar to ONETWO input, but for now these variables
**    have no effect, since the spline values of te, ti, zeff
**    are all assumed to be input on the same ryain(1:njene) mesh.)
**
**  iprote(or ti) (or ne) = "asdex"      
**      In this case, the code has ASDEX-type exponential profiles
**      hard-wired into the code. The user need not make further
**      specifications. See subroutines tdxinitl and tdpro.
**      Uses acoefne,acoefte.  
**      default: none
**
**  iprozeff = "disabled" (default), "parabola", "prbola-t", 
**                                   "spline-t" or "spline":
**           
**           = "disabled", then this option has no effect.
**           = "parabola","prbola-t", "spline"  or "spline-t", then this 
**             option gives ion densities so as to achieve a given zeff
**             radial profile, consistent with charge neutrality and the 
**             given electron density profile.
**
**             There must be at least one Maxwellian ion species
**             for this case.
**             If there are more than one Maxwl ion species with
**             bnumb() the same as the first (e.g., H+, D+, and/or T+),
**             then these will be treated as a unit, and their densities
**             adjusted proportionately, as required.
**             It is necessary to specify these equal-bnumb() species
**             successively at the beginning of the Maxwl ion species,
**             and the densities reden(,) should contain the relative
**             ratios of each species: e.g., reden(k,)=0.99,
**                reden(k+1,)=0.01 gives 99% 1st species.
**             Only up to two different Z ions are implemented for this
**             option.
**             If there are two different bnumb() Maxwellian ion species
**             specified, then densities of these species are adjusted 
**             to obtain the zeff-profile.
**             If only one bnumb() value (for .ge.1 Maxwl ions) is 
**             specified, then bnumb,fmass are as input above into
**             namelist, and an additional ion species is automatically
**             added with  bnumb=50., fmass=50*2.*1.67e-24.
**             The densities of the two different bnumb() species will
**             be adjusted to achieve the specified profile of zeff.
**             The first general ion species (if it exists)
**             will be given density of the first of the two Maxwellian
**             species.
**             The second general ion species (if it exists) will be given
**             the density of the second Maxwellian ion species.
**             Etc.
**             (remember that an added species will have bnumb=50,
**             fmass=50*2.*1.676e-24.)
**             In particular,
**             = "parabola", then input zeff profile by
**               zeffin(0)=central value, zeffin(1)=edge value,
**               npwrzeff,mpwrzeff (real numbers) specify parabolic
**               profile as is done above for reden.
**             = "spline", then input zeff profile data in 
**               zeffin(1:njene), as done above for enein.
**               ***NB***  MUST explicitly put zeffin(1)=,in this case.
**             If using time-dependent profiles, then need
**             iprozeff="parabola", and input profiles as specified
**             below.  iprozeff="spline" is not presently implemented
**             for time-dep. case.
**
**           = "curr_fit", then adjust zeff in such a way as to maintain
**             a target current. The target current is set 
**             with totcrt(nbctime) array. The profile of zeff(lr) is set
**             to be flat during such adjustment.  The adjustment procedure 
**             is similar to adjustment of electric field, e.g.,
**             for efswtch.eq."method4". The new value of Zeff is set as
**             Zeff_new= Zeff_old*(1 -zrelax*sign_dj*abs(djrel)**zrelax_exp)
**             where the relative "deviation" of computed current currxjtot
**             from the target current totcurtt at given time step is
**             djrel= (totcurtt-currxjtot)/max(abs(totcurtt),abs(currxjtot)) 
**             The computed current currxjtot includes the current from
**             FPE-solution (from curr() array), the bootstrap current
**             and also delta_sigma*E current 
**             [current addition from (sigma_coll-neo -sigma_banana)*E_phi]
**             The relaxation procedure for Zeff can be "tuned up" 
**             by two namelist variables: 
**           zrelax=  (default is 0.5d0), and
**           zrelax_exp= (default is 1.d0). See equation for Zeff_new above.
**             Present procedure may not work in certain instances,
**             e.g., when the runaway current becomes ~99% of the total current.
**             We might consider adding a Newton iteration in future
**             (just for time steps when the relaxation breaks down).
**             In other cases, the zeff may become unrealistically high.
**             Presently it is limited by impurity ions with bnumb=100.,
**             which are added automatically to the list of Maxwellian species;
**             their density is adjusted to keep the electron density as set by
**             the namelist (e.g. using iprone='spline-t') and to obtain the 
**             calculated value of zeff (Zeff_new in the above equation).
**
**  zeffscal=scale factor for zeffin spline input  and zeffc/zeffb
**           (default=1.0).
**
**  zeffin, npwrzeff, mpwrzeff, as above.
**
**  izeff = "backgrnd" (default), then zeff(l),l=1,lrzmax profile
**                 is calculated from all species k=1,ntotal except
**                 general electrons and Maxwellian electrons.
**                 [BH070316:  Need to check ngen.gt.1 case]
**        = "ion", then zeff(l),l=1,lrzmax calculated from all species
**                 other than general species or electrons. (Use for
**                 ion simulations, with colmodl=1,3...)
**             BH recommends checking this out a bit more in ainpla.f,
**             diaggnde.f and diaggnde2.f, if you are going to use it.
**
**  fpld(1:6,1:ngena) controls "additional" loading of the of the
**      equatorial distribution functions, at time step nfpld.
**      This set of options was instituted for purposes of study of the
**      electron knock-on runaway problem. 
**      It can also be used for INITIALIZING the 
**      Fokker-Planck'd distribution (ngen=1)
**      from externally generated distributions (see, for example,
**      the IDL procedure make_input_distn.pro, BobH, 001124).
**      (The default, fpld(1,1:ngena)=0.0,
**      gives no  loading of distribution functions beyond initial
**      (relativistic) Maxwellians derived from the specified density, 
**      temp. and zeff.)  
**      The additional loading 
**      is in the form of a internally generated vel-space-constant
**      plus a shifted "Maxwellian" distribution which is localized in
**      pitch angle, or,  distributions read in from file "fpld_dsk".
**      (See subroutine finit.)
**      fpld(1,1)= -1.0,  (Previously "fpld_dsk" option on C90)
**                 then read initial distribution function
**                 in from file "fpld_dsk" according to format
**                 given in source file finit.f.
**                 (Only works for  No additional preloading).
**      fpld(1,1)= -2.0   (Previously "fpld_ds1" option on C90)
**                then add additional loading with distribution
**                given in file "fpld_dsk", and renormalized to have
**                a fraction of the total (reden-)density 
**                as specified in fpld(2,1). Windowing of distribution is
**                given by fpld(7:10,k), as described below.
**      fpld(2,1)=in the case fpld(1,1)=-2.0, gives fraction
**                of additional loading which is in the distribution
**                specified by the "fpld_dsk" file.
**      For fpld(1,1)=real number.gt.0., (.ne.-1.0 or -2.0):
**      fpld(1,k)=fraction of flux surface averaged density of species
**                k which is in the additional preloading. (.ge.0., .le.1.)
**                (default=0.)   
**      fpld(2,k)=fraction of the additional loading which is in the
**                the shifted "Maxwellian".  (.ge.0., .le.1.)
**                (default =0.)
**      fpld(3,k)=Energy of peak of shifted "Maxwellian". (keV)
**      fpld(4,k)=Energy spread of distn. (keV).
**      fpld(5,k)=Pitch angle of peak of shifted "Max". (radians)
**      fpld(6,k)=Dispersion of distn in pitch angle (radians).
**      fpld(7,k)=Minimum energy(keV) of loading function. (deflt. 0.)
**	fpld(8,k)=Maximum energy(keV) of loading function. (deflt. 1.e10)
**      fpld(9,k)=Minimum  pitch angle(radians) of loading function (def. 0.)
**      fpld(10,k)=Maximum pitch angle(radians) of loading function (def. pi)
**          N.B.  If using Zeff (as above) increases the number of
**                ion species, must give any desired additional
**                preloading through the appropriate fpld(k,...)
**
**      nfpld = time step at which to introduce the additional loading
**              into f.
**
**************************************************************
**
**  TIME-DEPENDENT PROFILES
**
***************************************************************
**
**  Time-dependent profiles are implemented either by specifying
**    (A) central and edge values for parabolic-type profiles
**    as a function of time, or as (B) a series of spline
**    profiles as described at time-dependent splines, below.
**
**  (A) central and edge values for parabolic-type profiles
**  =======================================================
**    nbctime > 0, and profile designator (iprone, etc.) is "prbola-t"
**    The exponents of the profiles are fixed in time,
**    and are those used elsewhere in the code.
**
**    BH131029:  Presently, it does not appear their is backward
**               compatability with iprone='parabola', nbctime.gt.0.
**               Need to fix this.
**    
**  Two methods are implemented for specification of the time-
**    dependent parabolic profiles, according to tmdmeth="method1"
**    or "method2".
**  tmdmeth="method1", then specify central and edge values at a
**            sequence of times, bctime(1:nbctime), as given below.
**         ="method2"  give central and edge values at only
**            two times bctime(1:2).  The edge and central
**            values are constant up to bctime(1), vary as
**            y=y2+((y1-y2)/2.)*(1.+cos((t-t1)/(t2-t1)) up to 
**            time bctime(2), and then constant.  y1 and y2
**            are given by redenc(1:2), and so on.
**            (default="method1")
**
**  nbctime=number of times at which profile data given.
**     Set nbctime=0 for no time-dependent profiles, else
**     specify time dependent profiles.  (default=0)
**
**  bctime(1:nbctime)=times at which data given (seconds).
**    (default= 0, 1., 2., ..., float(nbctimea-1))
**    If code time .lt. bctime(1), use bctime(1) profiles.
**    If code time .gt. bctime(nbctime), use bctime(nbctime) profiles.
**
**  bctshift=time shift of bctime(array), to facilitate the present
**    configuration of cql3d which starts each run at time=0, even
**    for restarts.  Can lookup previous end time(s) in the time()
**    varaible of the  previous .nc o/p file.  Future work should 
**    enable time to be continuous, with restart.  (default=0.d0)
**
**  bctimescal=scale factor for bctime, for convenience of stretching
**    background plasma time variation.  Applied after any bctshift.
**    (default: bctime_scale=1.d0)
**
**  Setting central values =0. at the first time step
**    turns off profile generation through this namelist variable.
**    (default is all central and edge values =0.)
**    BH131029:  Looking over profiles.f, this option does not appear
**               to be implemented.
**
**  redenc(1:nbctime,k), redenb(1:nbctime,k) central and edge values
**    for each species k (cm**-3).  The powers of the parabola are
**    given by npwr(0), mpwr(0) as above, and are time-independent.
**    NOTE: for time-dependent density profiles, most sensible cases will
**        use lbdry="scale". 
**    With time dependent density and lbdry="scale", density is added
**      to achieve specified density, at either the local (in time and
**      position) temperature, or at the input value of
**      temperature as specified by temp_den.
**  temp_den=0., use local in time and radius temperature for added density.
**           .ne.0., then use this as temperature of added density (keV).
**           (default= 0.0).
**  tempc(1:nbctime,k), tempb(1:nbctime,k) central and edge values
**    for each species k (keV).  npwr(1:ntotal), mpwr(1:ntotal ) 
**    as for temp( ,k).
**  zeffc(1:nbctime),zeffb(1:nbctime), with real*8 npwrzeff,mpwrzeff,
**    specify zeff profiles as described above.
**  NOTE: if using this option (i.e., zeffc(1).ne.0.), need to have 
**        iprozeff="prbola-t"  (see iprozeff above).
**        Also, can specify two maxwellian ion species,  one
**        with high enough Z, or the code will will add a second
**        ion species with Z=50.
**  elecc(1:nbctime),elecb(1:nbctime), with real*8 npwrelec,mpwrelec
**    specify electric field profiles, unless current profile
**    control is enabled (volts/cm).  See below, D.C. ELECTRIC FIELD TERM.
**  vphic,vphib,npwrvphi,mpwrvphi similarly, specify toroidal rot (cm/sec).
**
**  Parallel Current profile (rather than specified electric field) gives
**    electric field if efswtch.eq."method2" or "method4", and
**    (xjc(1).ne.0., or totcrt.ne.0.):
**  xjc(1:nbctime),xjb(1:nbctime),npwrxj(real),mpwrxj(real) specify 
**    target current density profiles.  
**    If the following totcrt=0., but xjc(1).ne.0., 
**    then these profiles are absolute (amps/cm**2).  
**    If totcrt(1).ne.0, then these give the current density profile,
**      normalized to totcrt total plasma current. 
**  totcrt(1:nbctime) .ne.0., then total plasma current is specified
**    (amps) according to a profile given by xjc,xjb,npwrxj,mpwrxj.
**  NOTE: if using this option (i.e., xjc(1).ne.0.), need to have 
**        iprocur="prbola-t" (default = 'disabled').
**        (However, iprocur is not in ZOW version of cql3d.) (Now it is? Needs testing)
**
**  (B) a series of spline profiles
**  ===============================
**  nbctime > 0, and the respective profile indicators (iprone, etc)
**    must be set to "spline-t", as appropriate.
**    Must have tmdmeth="method1".  Not set up for above "method2" case.
**    A series of radial profiles specified at the ryain(1:nprone)
**    are given corresponding to times bctime(1:nbctime).
**    Default values for the *_t inputs below are 0.0d0.
**  enein_t(1:njene,1:ntotal,1:nbctime) input as series of density
**    profiles.   If using iprozeff = "prbola-t" or "spline-t", densities
**    are set up as described with iprozeff and izeff namelist input.
**    Can scale profiles with enescal.
**    If time-dependent spline ion densities, assume 1st in species,
**    at least.
**  tein_t(1:njene,1:nbctime) input as series of Te profiles (keV).
**    Can scale profiles with tescal.
**  tiin_t(1:njene,1:nbctime) input as series of Ti profiles (keV).
**    Applies to all ion temperatures, although coding can be readily
**    adjusted for individual species if needed.
**    Can scale profiles with tiscal.
**  zeffin_t(1:njene,1:nbctime) input as series of Zeff profiles.
**     Used if iprozeff.eq."spline-t".
**  vphiplin_t  (1:njene,1:nbctime) input as series of toroidal
**     rotation profiles (cm/sec).
**  elecin_t(1:njene,1:nbctime) input as series of Tor. E-field
**     profiles unless current profile control is enabled (volts/cm).
**     See efswtch="method1".
**  ennin_t(1:njenea,1:nbctimea,1:npaproca)  input as a series of
**     neutral density profiles (cm^-3), impurities, etc.
**     Used if ipronn='spline-t'.
**     This array can be used for calculations of ion losses due to 
**     Charge-eXchange with neutrals, but this feature is not implemented
**     in the ZOW cql3d. 
**     Also can be used for NPA diagnostics (npa_diag.ne."disabled"). 
**     Can be scaled with ennscal.
**     First index is the radial grid index, second is time-points index
**     (use  nbctime=1 and bctime(1)=0.0  for a steady-state profile).
**     Third index corresponds to type of CX (up to npaproc=5 for NPA).
**     The values of ennin_t can be specified in cqlinput, 
**     or can be read from data file. To do this, specify in cqlinput:
**       read_ennprofile_times='enabled',
**       nennprofile_files= [specify number of data files to be read]
**       Also, for the array ennprofile_files(1:nennprofile_files), 
**       specify names of data files. 
**       The files should be in text form, with 1st and 2nd line
**       being headings containing description of data,
**       followed by two columns with data: 1st column is rho_sqpolflx
**       and 2nd column is profile of neutral density [cm^-3].
**       The data is read and processed by 
**       subroutine read_convert(filename,a_new): 
**       1. Read data from files containing rho-coordinate and profile 
**       of neutral density (or other profile, in general/future devel).
**       The name of the file is given by the input 'filename'.
**       The files should be in text form, with 1st and 2nd line
**       being headings containing description of data,
**       followed by two columns with data: 1st column is rho_sqpolflx
**       and 2nd column is profile of neutral density [cm^-3] 
**       (or other profile, in general/future development).
**       2. Convert the rho_sqpolflx coordinate from file into "working"
**       rho_radcoord coordinate, which type is specified by radcoord.
**       The size of rho_radcoord is same as that of rho_sqpolflx.
**       3. Interpolate rho_radcoord into rho_new of size=njene
**       and also interpolate original density profile  
**       into the new profile a_new of size=njene.
**       The array a_new(njene) is the output of subroutine. 
**       It will be put into array (for each time index itme):
**       ennin_t(1:njene,itme,1)= a_new(1:njene)  
**       The last index (=1) is for neutrals (H or D). 
**       Higher indices (2,3,4,5) will be used for setting profiles of
**       other neutrals or impurities (not ready as of May, 2013). 
**
**
**  Parallel Current profile (rather than specified electric field) gives
**    electric field if efswtch.eq."method2" or "method4", and
**    (xjin_t(1).ne.0., or totcrt.ne.0.):
**  xjin_t(1:njene,1:nbctime) specifies target current density profiles.
**    If the following totcrt(1)=0., but xjin_t(1).ne.0., 
**    then these profiles are absolute (amps/cm**2).  
**    If totcrt(1).ne.0, then these give the current density profile,
**      normalized to totcrt() total plasma current. 
**  totcrt(1:nbctime) .ne.0., then total plasma current is specified (amps)
**    according to a profile given by xjin_t.
**  NOTE: if using this option (i.e., xjin_t(1).ne.0.), need to have 
**        iprocur="spline-t"  (default='disabled').
**  
**
**  ADDITIONAL option of exponential decay of Te or Ti, 
**  iprote="prb-expt" (and similarly, iproti="prb-expt");
**  also - similarly for 
**  iprote="spl-expt" and iproti="spl-expt".
**    This option will set the change in temperature as
**       T(t)= Tend +(T(tstart)-Tend)*exp(-(t-tstart)/tau) 
**       for electrons (use iprote=) or ions (use iproti=).
**       where exp(-(t-tstart)/tau)  is applied only at t.ge.tstart.
**       See H.M.Smith and E.Verwichte, PoP vol.15, p.072502, (2008),
**       Eqn.(7).
**    This model is useful for cases of impurity deposition, 
**    see namelist variable imp_depos_method.
** 
**    The ending (lower limit) temperature is set with 
**  temp_expt_Tend=0.010d0 [default value] [keV] final(ending) Tend after cooling.
** 
**    The start of the T-drop can be surface-dependent, through 
**  temp_expt_tstart(1:lrz)= [sec] tstart in the above Eqn., for each radial bin.
**    Default value is 0.d0 for all ll=1:lrza, which means that
**    the drop starts at t=0.0 at all magnetic surfaces.
**    This setting is suitable for imp_depos_method="instant"
**    (instant deposition of impurity at all magnetic flux surfaces).
**    In case of imp_depos_method="pellet", there is no need to set these values;
**    they will be calculated during run to match the pellet position, 
**    so that temp_expt_tstart is the time when pellet reaches a given surface.
**    Specifically, for the flux surface with major radius rpcon(ll),
**    it is calculated as 
**    temp_expt_tstart(ll)= pellet_tstart + (pellet_Rstart-rpcon(ll))/pellet_V
**
**    The decay time ("tau" in the above equation) is set by namelist variable
**  temp_expt_tau1=0.1d-3 [default value] [sec] fast decay time of T(t)
**    (Suggestion: use a typical value for Thermal Quench)
**
**    Additionally, for imp_depos_method="pellet", 
**    there is an option to add a slow decay of T, from t=pellet_tstart.
*     The reason is that the temperature profile start collapsing
**    as soon as the pellet is launched into plasma, 
**    rather than at the instant when pellet reaches a given surface.
**    Specify
**  temp_expt_tau0=3.0d-3 [default value] [sec]slow decay time of T(t) 
**    The temperature will be dropping at a rate
**    expt= exp(-(t-pellet_tstart)/temp_expt_tau0)
**    at t>pellet_tstart, simultaneously at all flux surfaces.
**    Note: to skip such a "slow decay", simply set temp_expt_tau0 
**    to a very large number.
**    Thus, with the two values set (temp_expt_tau0 and temp_expt_tau1),
**    the T profile will be slowly decaying as a whole
**    at a rate of temp_expt_tau0, and additionally it will start 
**    collapsing at a rate of temp_expt_tau1, 
**    starting from the T(t,rho) value that was obtained at that surface
**    by the time when pellet reached it.
**    
**    Similarly in case of ion temperature, we set iproti='prb-expt' ,
**    and we use same values of temp_expt_Tend, temp_expt_tau0, temp_expt_tau1,
**    and temp_expt_tstart(1:lrz).
**  
**    Note that in equation
      T(t)= Tend +(T(tstart)-Tend)*exp(-(t-tstart)/tau) 
      the value of Tend is set with namelist variable
      temp_expt_Tend, and it is one and the same for all radial surfaces
      [could be changed in future development].
      As for the value of T(start) in the above equation, 
      it is stored in arrays temp_expt_T1(k,ll) for the fast decay,
      or in temp_expt_T0(k,ll) for the slow decay.
      Presently these two arrays are set in subr. tdxinitl, 
      at t=0 (ot t=timet in a restart type of run).
      For iprote='prb-expt' (or iproti='prb-expt')
      these arrays are set using parabolic profiles for temperature, 
      based on values of temp(k,lr=0) and temp(k,lr=1) in cqlinput;
      (Same settings as for iprote.eq."parabola").
      For iprote='spl-expt' (or iproti='spl-expt')
      these arrays are set using spline profiles for temperature,
      based on values of tein_t(1:njene,1:nbctime)
      or tiin_t(1:njene,1:nbctime).
      (Similar input arrays as in case of iprote='spline-t'
       or iproti='spline-t').     
      This option ('spl-expt') is work in progress.
      Presently, for the value of T(start), it uses only one 
      time slice of tein_t(1:njene,itme)
      or tiin_t(1:njene,itme) arrays. Normally, 
      it is taken at t=0 (then, itme=1).
      In case of restart type of run, it can be other time slice,
      corresponding to the physical time of restart run.
      It means that in the present setup,
      the values of tein_t (or tiin_t) at other time indexes are not used.
      This is ok if the exp()-decay starts from t=0,
      but may fail if it starts with some delay, 
      and the profiles of T are evolved by that time.
      This issue is easy to fix in case of imp_depos_method.eq.'instant',
      but not so easy in case of imp_depos_method.eq.'pellet',
      when the temperature profile has to be "stitched" from two pieces -
      one piece is from tein_t profile where pellet did not reach yet,
      and another piece - from the exp()-decaying part.
      [Work in progress, 2022]
      
**
**************************************************************
**
**  IMPURITY MODELS  [developed in 2019]
**
***************************************************************
**
**  
**  imp_depos_method  is a method of deposition of impurity:
**    "disabled", or "pellet", or "instant".
**    Default is "disabled"
**    [There is an alternative method to read data on impurity
**     from external data files, e.g., computed by NIMROD.
**     See namelist variable  read_data= .. below]
**
**    For imp_depos_method="instant", it is assumed that the material
**    is deposited instantly at all radial bins (flux surfaces).
**    For this method, user needs to specify the profile of
**  dens0_imp(0:lrz) = Density profile of deposited impurity [1/cm^3];
**    This is just the ablated material (e.g. from pellet, flake, etc),
**    before ionization process occured. 
**    Index 0 corresponds to magnetic axis.
**    The densities of ionization states are calculated by Coronal model 
**    (see below); they can change in time, depending on T_e and n_e.
**    Also, for imp_depos_method="instant", need to specify
**  tstart_imp= Instant when impurity is deposited [sec].
**
**    For imp_depos_method="pellet", the profile of deposited impurity  
**    is calculated during run, based on parameters of pellet 
**    (mass, speed, ...;  see below). So, dens0_imp() is not needed.
**    For imp_depos_method="pellet", there is no need to specify
**    tstart_imp, it is obtained from tstart_imp=pellet_tstart 
**    (launch time of the pellet).
**
**  imp_ne_method= Method of adjusting the electron or main-ion densities
**    when an impurity is deposited, and densities of ionization states 
**    are changed in time.
**    "ni_list" or "ne_list". Default is 'ne_list'.
**    These two options correspond to the two ways of adjusting
**    electron density:
**     1. 'ni_list' option. Assume that density of main ion species 
**        (not impurity ions) is taken 
**        from input namelist (could be time-dependent); 
**        then, electron density is set to 
**          sum(n(k)*Z(k))[all k=kionm of main ions] +
**        + sum(nimp(kstate)*Zimp(kstate))
**          [all ionization states of impurity].
**     2. 'ne_list' option. Assume that electron density is 
**        taken from the input namelist, as ne(rho,t); 
**        then, reduce the density of main ions, 1st ionic species.
**        With increase of density of impurity ions, the density of 
**        main ions will go down, to maintain the value of ne from list.
**        Note: This procedure may fail if the density of main ions is low,
**        while the density of impurity is high 
**        (i.e., too many electrons are coming from impurity ions).
**        In this case, the code sets the density of main ions
**        to a very small number of 1e5 cm^-3,
**           reden(kion_main,l)= max(reden(kion_main,l),1.d5) ![cm^-3]
**        and then the electron density is calculated effectively
**        from the impurity ions, so - it cannot follow a given ne(rho,t) list.
**        For details, see profiles.f.
**
**  For imp_depos_method="pellet" (pellet propagation/ablation model),
**  the following namelist values should be set.
**    
**  pellet_Rstart= [cm] Major radius where pellet is launched.
**      Default is 230. Suggestion: Set it to rpcon(lrz), or R_LCFS radius.
**  pellet_tstart= [sec] Instant when pellet is launched, in code time.
**      Default is 0.d0.  Can be .ge.0.
**  pellet_V= [cm/s] Pellet speed. Typically 10000-900000cm/s
**      Default is 30000.d0. Assumed constant all the way through plasma.
**      Assumed that pellet travels along equatorial
**      plane, going through magnetic axis.
**  pellet_M0= [gram] Initial mass of pellet(at t=pellet_tstart)
**      Default is 30.d-3. If pellet is large, 
**      it can make to the inner border of plasma.
**      
**      Related to pellet size and ablation cloud:
**      For distributing the ablated mass among several flux surfaces,
**      assume that the ablation cloud is 5--8 times
**      larger than the pellet itself.
**      Allow for assymmetry between leading (front) 
**      and trailing (back) side of the cloud.
**  pellet_rp0= [cm] Pellet radius at t=pellet_tstart (plasma edge).
**      Default is 0.5d0
**  pellet_rcloud1= [cm] Radius of ablation cloud, leading (front) side.
**      Default is 4.d0
**  pellet_rcloud2= [cm] Radius of ablation cloud, trailing (back) side. 
**      Default is 4.d0
**      Typically pellet_rp0== rp(0) = 0.2--0.5 cm.
**      Recommended: pellet_rcloud ~(5--8)*pellet_rp0
**      Pellet radius is reduced during propagation 
**      (as the mass is reduced).
**      However, in present model, rcloud is not changed, 
**      so the cloud size remains as described by pellet_rcloud1,2 above.
**     
**      Related to description of ablation rate:
**      Assume that the ablation rate of pellet is proportional to local 
**      ne^pn * Te^pt (electron T and density in some powers),
**      and proportional to remaining_mass/pellet_M0 in some power "pm".
**      So that the local ablation rate is 
**      G[gram/s]= 
**       =Cablation* ne[cm-3]^pn *Te[keV]^pt *(Mpellet(t)/Mpellet(0))^pm
**      See REFS: "2019-03-15-Friday Science Meeting-Jie Zhang.pdf"
**      Values for those powers:
**  pellet_pn= power "pn" in the above Eqn. for G.
**      Default is 1.d0/3.d0  [REFS: should be 1/3]
**  pellet_pt= power "pt" in the above Eqn. for G
**      Default is 5.d0/3.d0  [REFS: should be 11/6, or 5/3 ]
**  pellet_pm= power "pm" in the above Eqn. for G
**      Default is 4.d0/9.d0  [REFS: should be 4/9 (so that rp^(4/3))]
**      where Mpellet(t)==Mpellet_rem is the remaining mass
**      at given radial position R(t).
**      Note that  (Mpellet(t)/Mpellet(0))^pm ~~  (rp(t)/rp(0))^(3*pm)
**      For example, when pm=2/3, we get  G~~ rp^2,
**      which means - proportional to surface area of the pellet
**      (S_pellet= 4*pi*rp^2).
**      The value of Cablation=="pellet_Cablation" is either calculated
**      during kopt=0 call of subroutine set_get_pellet, 
**      or set as a namelist value, see below.
**     
**      On calculation of pellet_Cablation value.
**      It could be either set as a namelist value, 
**      or found by an iterative procedure.
**  ipellet_method=1 [Default value]  Iterative procedure,
**      to find such value of pellet_Cablation which yields the value
**      of fraction of pellet remained at magnetic axis, 
**     pellet_fract_rem=0.5d0 [Default value] Fraction of remaining mass  
**        when pellet reaches magn.axis, i.e. it is   
**        pellet_fract_rem= (pellet_M0-dMpellet_sum(t_axis))/pellet_M0
**        where dMpellet_sum(t_axis) is the total ablated mass 
**        during the flight of the pellet from plasma edge to magn.axis.
**        The value of pellet_Cablation will be found from iterations.
**        For this ipellet_method=1, also specify these two values:
**     ipellet_iter_max=50 ! Max number of iterations
**     pellet_iter_accur=1.d-2 !Relative error (accuracy) 
**        achieved in iterations, to be compared with 
**        |pellet_fract_rem-pellet_rem_iter|/pellet_fract_rem
**  ipellet_method=2  Use the estimate (see below) for 
**      pellet_Cablation  coefficient, and let the pellet propagate
**      at such value.
**      The remaining mass of pellet when it reaches the magnetic axis
*       depends on the value of pellet_Cablation 
*       and on the initial mass pellet_M0.
*       The estimate is done from
*       pellet_Cablation= 39*(0.5)^pt *(1e-14)^pn *(pellet_rp0/0.2)^(3*pm) 
*       This is to make it consistent with expression given in
*       REFS: See "2019-03-15-Friday Science Meeting-Jie Zhang.pdf"
*       Units for pellet_Cablation : (gram/s) / (cm^-3)^pn / (keV^pt) 
**  ipellet_method=3  Specify the value of pellet_Cablation in namelist.
**      This is the easiest method. Simply specify, for example:
**    pellet_Cablation=0.001d0 [default value] [(gram/s)/(cm^-3)^pn /(keV^pt)]
**----End of imp_depos_method="pellet" (pellet propagation/ablation model) ---
**
**  For calculation of ionization states of impurity, 
**  presently based on Coronal model, the following 
**  variables should be set.
**
**  adpak='enabled' [Default value] To use ADPAK tables 
**    which are present in /ADPAK_data/ folder. 
**    If 'disabled', an alternative set of subroutines will be used, ADCDO)
**  imp_type=  Impurity type.
**    Presently ADPAK data is available for  2-Be,3-C,4-N,5-Ne,6-Ar.
**    Besides, with ADCDO subroutines (adpak='disabled') user can have 
**    the above types, and also imp_type= 1(He), 7(Xe), 8(W, Tungsten), 9(B)
**    Default is imp_type=6 (Argon).
**    For ADPAK and ADCDO subroutines/tables, need values of neutral deuterium:
**  model_dens_nD0=1  [Default value] Only one model so far for neutral D0. 
**  dens_nD0_b=1.0d10 [Default value] [1/cm^3] Edge density of neutral D0
**  dens_nD0_l=10.d0  [Default value] [cm] Scale length of exp-decay, as in  
**                           dens_nD0_b*exp[(rho-1)*radmin/dens_nD0_l]
**    Also for ADPAK subroutines/tables, need this:
**  adpak_tau_r= [sec] Characteristic time tau(r) of radial decay
**    of electron temperature T_e.
**    Default value = 1.d-3.
**    Note from A. Pigarov:
**      For disruption case, I would set tau(r)=1.e-3 sec. 
**      For Smith-like run case tau(r)=TauT, 
**      where TauT is the temperature decay time.
**      For quasi-stationary plasma it is likely about 
**      plasma confinement time ~1 s.
**
**    The remaining values for Coronal models (both ADPAK and ADCO)
**    are obtained by code from run-time values, such as 
**    electron density and temperaure, and ion temperature 
**    (it is assumed that all ions, including impurity ions, 
**     have same temperature).
**
**---End of Coronal model for impurity ---
**
**   The effects of impurity on collision operator are accessed through 
**   corrections in pitch-angle scattering coefficients and
**   energy slowing-down coefficients.  Use
**   gamafac="hesslow"  [default is gamafac="disabled"]
**   (See description of gamafac further in this file).
**   For test purposes these effects can disabled
**   by setting one or both of these variables to 0,
**    imp_bounde_collscat= 0 or 1 ! Use =1 to enable effects of scattering 
**       of electrons on partially-ionized impurity ions 
**       (Hesslow corrections, see  Hesslow et al, JPP-2018,vol.84).
**       Default value is 1.
**       Note: imp_bounde_collscat=0 will nullify gscreen(kstate,j) array
**       in subr. set_gscreen_hesslow.
**    imp_bounde_collslow= 0 or 1 ! Use =1 to enable effects of slowing down 
**       of electrons on partially-ionized impurity ions.
**       Default value is 1.
**       Note: imp_bounde_collslow=0 will nullify hbethe(kstate,j) array
**       in subr. set_gscreen_hesslow, which is  Nej*[(1/k)ln(1+hj^k) - beta^2] 
**       in the JPP-2018 paper, in Eq.(2.31).
**
**   Note: The effects of impurity are also included into KO source operator,
**   through the density of bound electrons - it is added 
**   to density of free (cold) electrons, as the target for fast electrons.
**   Also, using gamafac="hesslow" will add an energy dependence into 
**   Coulomb logarithm, according to Hesslow et al, JPP-2018, Eqs.(2.10).
**   See the ending part of subr. cfpgamma.
**
**   For effects of impurity on electron or ion temperature, presently  
**   there is an option of a prescribed T(t) exponential(in time) drop, 
**   see iprote='prb-expt' and iproti='prb-expt'.  
**   Also see iprote='spl-expt' and iproti='spl-expt'.
**
**
**   Alternative method to handle impurities - read data on impurity
**   from external data files, e.g., computed by NIMROD.
**
**  read_data='nimrod', 
**    To read data from set of files computed by NIMROD.
**    Other options can be added in future.
**    Default is 'disabled'.
**
**    Related namelist variables:
    read_data_filenames_prefix= 'rycq0'  (example; default is "notset")
    read_data_filenames_suffix= '00.dat' (example; default is "notset")
      These two variables are used by the code to form
      a template for file names. 
      Example: Consider data file 'rycq019800.dat', 
      it is made of prefix 'rycq0',
      followed by 3-digit number and suffix '00.dat'
      The code will try to find all files from  rycq000000.dat 
      to rycq099900.dat (max 1000 files, for now).
      The code will print the list of all files that are found.
      Note that the numbers in sequence of data files 
      may have gaps, like in this sequence:
      rycq020000.dat, rycq020500.dat, rycq021500.dat,
      so that not all rycq0***00.dat are present.
      The code will save the names of the above files into 
      character array read_data_filenames(1:3) [in this example].
      The data files from NIMROD contain
      FSA (or outboard-midplane) values of E field, Te, Ti, density of D+,
      densities of all ionization states of impurity (e.g., Neon), etc.
      Max number of files is 1000, for now. 
      For coupling with NIMROD, each file contains data at one time slice.
      Therefore, it is recommended to match the max number of files
      with value of nbctimea [set in param.h]. If too many, the run will be
      stopped with a warning message.
      See subr. read_data_files() for details on processing data
      from NIMROD data files.  In this subroutine,
      the data from each file is copied into CQL3D time-dependent arrays 
      (after change in units, adjustments of negative densities, etc.):
      enein_t() [for main ion species, D+, 
      and for electrons - general and Maxwellian],
      tein_t(), tiin_t() [Te and Ti profiles],
      elecin_t() [E field],
      dens_imp_t(0:nstates,1:njene,1:nbctime) [densities of each 
      ionization state]
      The radial grid ryain(1:njene) is formed automatically,
      using radcoord='rminmax' option (the only option in this case).
      The time array bctime(1:nbctime) is formed automatically
      by reading the value of time slice that is present 
      in each data file (at 1st line). 
      The value of njene (number of radial points) is determined
      automatically from the data in each file.
      The value of nbctime (number of time slices) is determined
      automatically as the number of data files.
      
    temp_min_data= Lower limit, to adjust Te and Ti data [keV]
      In NIMROD data, Te and Ti may go down to ~0.3eV; not physical.
      Default: temp_min_data=5.d-3 ![keV]

    elecfld_data= How E-field is set.
                ='readdata', then get E field directly from data file.
    elecfld_data='jspitzer', then get E as j_data/resistivity
    (the latter option was the only option in original setup).
      Default: elecfld_data='jspitzer', for backward consistency.
      (This variable, elecfld_data, is added by YuP[2022-07-13] )   
     
    ===> Note: Some other namelist values for this case 
         of read_data='nimrod' are verified and reset 
         in subr.ainsetva, particularly:
               iprone=  "spline-t"
               iprote=  "spline-t"
               iproti=  "spline-t"
               iproelec="spline-t"
               efflag="parallel"
               iprozeff='disabled'
               radcoord='rminmax'
               imp_depos_method='disabled' 
               [imp_depos_method is not applicable in this case]
         Also, it is verified that only one ion species
         (Maxwellian) is set in cqlinput, and it is D+. 

**
**************************************************************
**
**  COLLISIONAL INPUT
**
***************************************************************
**
**  ncoef is the frequency of computation of collisional coefficients.
**    We use "n" as a time step (or cycle) designator.
**    If mod(n,ncoef)=1 or n=1, then compute the coefficients.
**    default: ncoef=1
**
**  mx is the highest order polynomial used in
**    Legendre expansion of distribution and potentials for
**    purposes of computing the collisional coefficients.
**    default: mx=3 (in aindflt.f)
**
**  colmodl selects the collision model:
**  ===> NOTE:  The default value is colmodl=1 <===
**              Be careful to specify the value you want in cqlinput:
**              colmodl=3 is the common setting.
**  colmodl=0 means fully nonlinear self-collisional coefficients for
**    the general species and with the background species specified
**    in the namelist input).  Applied power can give continual energy
**    growth (in absence of collisional losses with other species,
**    or other losses such as radiation or radial transport).
**    The self-collisions conserve energy perfectly, 
**    in the limit of small grid spacing.
**    Only include background species for which it is desired to
**    include general species collsions on the background Maxwellians.
**    Thus for colmodl=0, normally, do not include a background 
**    species for the (possible multiple) general species.
**  colmodl=1 means only contributions from background included
**    in the Fokker Planck collision term.
**  colmodl=2 means only contributions from the general species are
**    included (background ignored).
**  colmodl=3 is like colmodl=0 but the p0 component (0'th order
**    Legendre polynomial) from the general species is ignored.
**    Instead, a background Maxwellian of the same species is set up
**    at a fixed density and temperature which effectively
**    supplies the p0 component for the general species. The effect
**    of this "standard" treatment is that momentum and particles
**    are conserved by the like-particle collision operator, but
**    the background Maxwellian of specified density and temperature
**    act like a heat sink, and stabilizes the Fokker-Planck'd 
**    distribution against thermal runaway.
**        This option (colmodl=3) can be useful for runaway electron 
**    and resistivity calculations.  Unlike the colmodl=0 case, 
**    which keeps accurate track of momentum
**    transfer but would allow the distribution to run away in a tau_e
**    time or so, this option keeps the bulk at the background
**    electron temperature, but still calculates momentum transfer
**    accurately. 
**        Use of colmodl=1 keeps the bulk from running away but
**    the momentum computed from the solution f will be wrong. 
**    colmodl=1 is the technique employed in the past by many codes
**    to compute DC electric field driven electron runaway rates.
**        Clearly, the colmodl=3 option is a kluge, and it does not
**    guarantee a positive definite operator. Thus if the bulk stays
**    cold but the distribution tail gets excited, the p1 component
**    might dominate the p0 component, and the velocity diffusion
**    coefficient cbl could become negative. Caveat emptor!
**        Note that the p0 component referred to above applies
**    to the contribution from a Maxwellian species at a fixed
**    temperature of the same type. For electron runaway calculations
**    or other calculations with colmodl=3,
**    one would have both general species electrons and background
**    electrons. The colmodl=3 option would "see" the contribution
**    from the background electrons as the p0'th component and the
**    other components would come from the general species electrons.
**  colmodl=4.....(undocumented and unsupported option)
**    default: colmodl=1
**
**  kfield(k)="enabled" means that species "k" is included as a field
**    (background) species in the calculation of the collision integral
**    coefficients.  kfield(k)="disabled" means it is excluded.
**    Usually, the background species are in the namelist data, in order
**    to be used as field particles in the collision operator, and 
**    thus will want to be included.
**    default: kfield(k)="enabled"
**
**  cfp_integrals='enabled' [Option added in 2020-07]
**    This option specifies how certain integrals are computed
**    which are used in subr. cfpcoefn.
**    These integrals describe contribution to BA coll. coeffs
**    from collisions of general species with the background 
**    Maxwellian species; the integrals only depend on mass
**    and local temperature of Maxwellian species.
**    Instead of calculating them over and over again
**    [as it was done originally; now accessed 
**    with cfp_integrals='disabled'], 
**    calculate them as a table over temperature grid 
**    and then reuse them by matching a local T
**    with the nearest values in the T-grid
**    [accessed with cfp_integrals='enabled'].
**    These integrals are updated if the temp() of 
**    the Maxwellian species is changed in time. 
**    Maxwellian species also include impurity ions 
**    with all ionization states [see description for 
**    namelist variable imp_depos_method].
**    This option can give some speedup in case of high-Z impurity.
**    Default is 'disabled' for backward compatibility. 
**
**  gamaset .gt. 0. means that all Coulomb logarithms are set to
**    gamaset. gamaset=0. means that the log is computed internally.
**    default: gamaset=0. **However, need to use gamaset.ne.0 for ion
**                        **general species.  This needs attention.
**  !YuP[2020-06-23] New option:
**  gamaset=-1.  means: use NRL definitions for gama(kk,k). 
**
**  gamafac = "enabled", then for electrons as (only) general species, 
**    multiply the Coulomb logarithms (e-e and e-i) in the
**    collision FP coefficients by and energy dependence as follows:
**     The energy dependent Coulomb log will be
**       alog(flamcql + min(sqrt(gamma-1),1)*(flamrp-flamcql)),
**       where flamcql is the argument of cql e-e Coulomb log,
**       and flamrp is the argument of the Rosenbluth-Putvinski
**       (Nucl. Fus. 1997) high energy electron Coulomb log.
**       The sqrt(gamma-1)-factor is chosen in accord with
**       the energy factor in the the CQL Coulomb log
**       (cf., Killeen et al. book.  The CQL Coulomb log
**       reduces to the NRL, Te.gt.10eV expression).  
**       This factor becomes significant at electron energies 
**       of order of greater than me*c**2.
** 
**  !YuP[2019-07] New option added/ still in progress [remove this line later]
**  gamafac= "hesslow", then use Eqs.2.10 in  J. Plasma Phys.(2018), vol.84
**       for Coulomb logarithms (e-e and e-i) with p-dependence (energy dependence),
**       but use our own definition of gama(kk,k) [the limit at p=0]
**       So, we set 
**          gamefac(j,k)= 1.d0 + 0.2*log(1.d0 + p_dep**2.5)/gama_ee
**       or gamefac(j,k)= 1.d0 + 0.2*log(1.d0 + p_dep**5)/gama_ei
**       where gama_ee= gama(kelecg, k=kelecg or kelecm),
**       and   gama_ei= gama(kelecg, k=kionm() or kiong()).
**       Those gama() are defined elsewhere in the code; see gamaset.
**       Also, if gamafac = "hesslow", use corrections for free_electron 
**       scattering on partially ionized ions or neutrals
**       [based on Eq.2.22 in the JPP paper],
**       and corrections for slowing down of free_e on bound electrons
**       in partially ionized ions or neutrals [based on Eq.2.31 in the JPP paper].
**  !YuP[2019-07] New option added/ still in progress [remove this line later]
**
**  qsineut = "enabled" or "maxwel" forces background electron density to
**    evolve to maintain quasi-neutrality in time. The midplane line
**    density of the electrons will equal the sum of the midplane line
**    densities of all the ionic species for qsineut="enabled".
**    For qsineut="maxwel" the sum is over background ions only.
**    Note that this option will not
**    apply if general species electrons are present. Use of this option
**    precludes use of locquas="enabled" below.
**    default: qsinuet="disabled"
**
**  locquas="enabled" means that if (1) electrons are not a general species
**    and (2) if qsineut="disabled" (above) then the background electrons are
**    not treated exactly as Maxwellians for purposes of computing their
**    contribution to the collision integral. Rather, they are assumed
**    to be local Maxwellians (at z along the orbit) with a density that
**    maintains local quasineutrality at z with all the ionic
**    species. A true Maxwellian is constant in density as a
**    function of z.
**    If locquas="disabled", this option does not apply.
**    default: "disabled"
**
**  trapmod="enabled", modifies the magnetic well on each flux surface
**    by transforming B/B_min ==>    B/B_min - trapredc*(B/B_min-1.),
**    where B is the magnetic field as a function of poloidal angle,
**    B_min is the minimum value on the flux surface.
**  trapredc is the fraction by which the magnetic well is reduced.
**    Thus, trapredc=0. gives no reduction.
**    trapredc=1. gives complete reduction, i.e., constant B (no well).
**    This option is to help in a heuristic evaluation of the
**    effect of trapping.
**    Have had problems with the theta mesh for trapredc.gt.0.99
**                                       (BobH, 981018). 
**    default: trapmod="disabled"
**  trapredc= fractional reduction in magnetic well on each flux surface.
**    default: trapredc=0.
**
**  scatmod="enabled" modifies the collisional "pitch angle scattering"
**    F coefficient, according to scatfrac, below. 
**    Setting mx=0 gives zero value to the 
**    collisional C, D, and E coefficients (the Rosenbluth g and h
**    functions become independent of pitch angle).  F is the only
**    remaining coefficient of a pitch angle derivative term.
**    F is multiplied by scatfrac, thus scatfrac=0. shuts off
**    the collisional scattering.
**    (default: scatmod="disabled").
**  scatfrac= multiplier of collisional C,E,F coeff, for testing purposes.
*     Also applied to rf coefficients.
**    (default: scatfrac=1.  Also, reset to 1. if scatmod.eq."disabled:).
**
********************************************************************
**
** ENERGY AND PARTICLE (KROOK) LOSSES
**
********************************************************************
**
**  There are three classes of particle and/or energy losses, set
**    by (1)lossmode(k), (2) torloss(k), and (3) tauegy().gt.0.:
**
**  (1)Several models are available and are invoked by setting
**  lossmode(k) to the chosen "string" in the input:
**  "snk0"==>  orbits whose energy is less than that corresponding
**    to normalized velocity x=xsink (namelist input) are loss 
**    orbits and are lost from the system with characteristic time equal to
**    the bounce time.
**  esink/xsink are alternate means for setting a loss boundary.
**    If esink.ne.0., then this energy (keV) specifies the boundary.
**    Else, use xsink (normalized velocity to vnorm). Used in cases 
**    lossmode(k)=snk0 and mirrsnk1.
**  "snk0accl"==>same as "snk0" except the loss mechanism is acceler-
**    ated by factor xsinkm (input).
**  "electron"==> forces a loss cone for general species electrons only
**    if the parallel energy is greater than eparc(k,l_) or the 
**    perpendicular energy is greater than eperc(k,l_).
**  "mirrorcc"==>potential square well model for tandem mirror central
**    cell calculations; orbits passing through a potential jump
**    at the throat of magnitude ephicc (namelist, below) are lost.
**  "mirrsnk"==> means both mirrorcc and snk0 options are used.
**  "mirrsnk1"==>means a loss cone as for mirrorcc, except
**     particles with velocity less than esink/xsink are not lost.
**    default: lossmode(k)="disabled"
**  "simplban"==> then trapped general species ions (or electrons) with
**    (banana width (poloidal gyroradius) plus gyroradius) greater than
**    distance to the (outer flux surface)*simpbfac are lost in a 
**    bounce time. See subroutine losscone for a few more details.
**    xxx As of 160602, this option uses deltarho(,,) first order orbit
**    xxx width calculation to obtain the banana width, for scrapeoff.
**    xxx ADVISEABLE to use ndeltarho='freya' or 'enabled', if using NBI.
**    xxx   [This shifts NBI birth points to their BA location, from which
**    xxx   the "simplban" loss is calculated.]
**    As of 170715: "simplban" reverted back to coding (effectively) of
**    pre-160602.  Part of problem is that the deltar coding is not presently
**    set up for k.gt.1.
**  "simplbn1" and "simplbn2" are variations on "simplban". "simplbn2"
**    was "simplban"-option up to 160403. "simplbn1" is a slight refinement
**    of "simplbn2".
**    170715: "simplbn2" now gives identical results to "simplban".
**            Presently, only recommending use of "simplban".
**            Future work can use the more accurate cql3d-hybrid-FOW code 
**            FI loss calculation to benchmark these simple models.
**  simpbfac=factor times radius of LCFS, modifying loss radius at
**    which particles are scraped off, for "simplban" losses.
**    For example, 1.05 increases the radius by 5%.  (The physical 
**    distance will depend somewhat on the value of radcoord.)
**    default: simpbfct=1.0
**  "frm_file"==> then read lossfile, specified as below.
**     [future options could be "file1",  etc. for different formats].
**  default: lossmode(k)="disabled"
**
**  lossfile(k) for lossmode(k)="frm_file", gives the
**    character*256 ascii path to a file specifying
**    region of prompt particle loss (integer 0 as function of
**    u_par,u_perp).  See file format in subroutine losscone.
**    The file gives loss-regions versus upar-uperp on a radial set
**    of flux surfaces (presently to be same as cql3d rya(1:lrz) grid).
**    Max velocity (vc_cgs) for u_par (=max(u_perp)) is specified 
**    in the file and may be different than given by enorm (for cql3d).
**    If vc_cgs.lt.vnorm, then loss regions at vc_cgs are extended
**    to vnorm at constant pitch angle.
**    Except, if data is given on twice as many flux surfaces as
**    specified by lrz, then the data set is reduced by omitting
**    data for every second flux surface.
**    default: "./prompt_loss.txt"
**
**  eparc(k,) and eperc(k,): see lossmode="electron" above.
**   They are set in the same manner as reden is set, by specifying
**   eparc(0,k) and if lrz.gt.1 eparc(1,k) also.  Similarly, eperc.
**
** (2)The torloss(k) series provides a loss term -f/tloss,
**    treated time-implicitly (for numerical stability).
** nonloss, and noffloss give time steps for start and stop
**    of these loss terms.  (defaults: nonloss=0, noffloss=10000).
**    torloss(k)="velocity" provides for an additional toroidal loss 
**    term, and the toroidal loss time will be determined by an 
**    ad hoc formula determined from input values tauloss(1-3,k):
**    tloss=tauloss(1,k)*(1.+tauloss(2,k)*(v(j))**tauloss(3,k))
**    where v(j) is the j'th velocity (momentum/mass) mesh point.
**    v(j)=x(j)*vnorm, x(j) code variable.
** tauloss(1:3,k)  loss rate coefficients (default=0.)
** torloss(k) = "flutter" gives a model similar to "velocity", but
**    dividing rather than multiplying the x part, to be able to model
**    Rechester-Rosenbluth type loss (with tauloss(3,k)=1):
**    tau = tauloss(1,k)/(1 + tauloss(2,k)*(x*vnorm)**tauloss(3,k))
**    Losses only applied to transiting particles for the following 
**       "fluttern"  models  (check the code).
** torloss(k) = "flutter1", further multiply tau(flutter) by gamma**5.
** torloss(k) = "flutter2", loss time tau is modelled on
**    electron diffusion in a stochastic magnetic field of
**    strength deltabdb, defined below.
**    The loss time is 
**       tau= (distance from plasma edge)**2/(4*v_par*D_st),
**       D_st= the magnetic field diffusion coefficient,
**           =pi*R_eff*deltabdb**2
**       1/R_eff = 1/(pi*q_safety*R_mag) + 1/lambda_mfp
** torloss(k) = "flutter3", further multiply tau(flutter2) by gamma**5.
**       default: torloss="disabled"
** deltabdb=the perpendicular fluctuating B field amplitude 
**    divided by the ambient B-field.
** torloss(k)="energy" means all particles with energy less than
**    enloss(k) (Kev) will have a loss time of tauloss(1,1)
** torloss(k)="relativ" means the loss time goes like
**   tauloss(1,k)*gamma**3
** torloss(k)="shell" is like "energy" above, but with the additional
**    effect that all particles with an energy greater than
**    .9*enorm are also lost.
** torloss(k)= "disabled" disables this option.
**    default: torloss(k)="disabled"
** enloss(k)=See torloss(k)="energy" above.
**
** (3)tauegy(k=1,ngen,0:1) .gt. 0.0 gives energy loss profiles, 
**    according to bounce average of the phenomonological 
**    energy loss term:
**            df/dt = v**-2*d/dv[1./tau(v,theta,lr_)*v**3*f/2.]
**    where   tau is energy loss time specified by
**            tau=rfact*tauegy*gamma**gamegy, v.le.vte
**                rfact*tauegy*gamma**gamegy/
**                (abs(v_par,lr_)/vte)**paregy*(v_per/vte)**peregy
**                            *(v/vte)**pegy),  v.gt.vte
**            rfact is a radially dependent modifier, as below.
**    tauegy(k,0), tauegy(k,1) are expanded out across lrzmax flux surfaces
**      parabolically, according to
**         (tauegy(k,0)-tauegy(k,1))*(1-(rho/a)**negy(k))**megy(k)+tauegy(k,1)
**  negy(1:ngen), megy(1:ngen), as above.
**  gamegy(1:ngen), paregy(1:ngen), peregy(1:ngen), pegy(1:ngen), as above.
**
**
**
**  Also have possibility of synchrotron and bremsstrahlung reaction
**    force on the electrons:
**
**  syncrad="gyro" gives synchrotron radiation energy loss
**            force on electrons according to formula of 
**            Bernstein and Baxter, Phys. Fl. Vol. 24, p. 108 (1981).
**          "geom" gives synchrotron radiation energy loss
**            force on electrons according to formula of
**            derived (BobH+Paul Parks) to account for accelerations
**            due to toroidal and poloidal following of field lines.
**          "gyrogeom" gives both "gyro" and "geom" effects.
**          (default: "disabled")
**
**  bremsrad="enabled" gives bremsstrahlung radiation energy loss 
**           force on the electrons. (combined with above phenomenological
**           energy losses).  (default:"disabled")
**           Reference: Physics Vade Mecum, H.L.Anderson, Ed. Sec. 16.07.F.
**           brfac,brfac1,brfacgm3 are additional factors which control
**           onset of an ad hoc high enhanced Bremsstrahlung which is
**           used to "catch" runaway electrons before they go off the
**           the edge of the grid.
**
**  brfac,brfac1,brfacgm3: ad hoc constants to enhance bremsstrahlung 
**    at the high energies, beyond (gamma-1)>gam4, where
**    gam4=brfacgm3*16.43/Z**(1/6), and brfacgm3>1.
**    The cross-section 'sigmarad' is evaluated as in
**    above PVM reference in the first and second of three energy ranges. 
**    In the last (upper) range, beyond (gamma-1)> brfacgm3*16.43/Z**(1/6),
**    it is multiplied by a factor (1.+brfac*gam5**brfac1) where
**    gam5 is difference (gamma-1)-gam4 (see above reference, or coding in 
**    subroutine lossegy.   The purpose of this is to force containment 
**    of high energy runaway electrons at energy greater than
**    (gam4-1.)*m*c**2, where m*c**2=510.98 keV.
**    (For investigation of knockon avalanching, etc.).  
**    (defaults: brfac=0.,brfac1=0.,brfacgm3=1.0)
**
**********************************************************************
**
**  D.C. ELECTRIC FIELD TERM
**
**********************************************************************
**
**  efswtch is a switch for the method for specifying the toroidal
**    electric field:
**    ="method1" ==> direct specification of the electric field
**                   by elecfld, eoved, or time-independent 
**                   elecin (with iproelec="spline"),or time-dependent
**                   profiles (with iproelec="prbola-t" or "spline-t"),
**                   elecc, elecb, npwrelec, mpwrelec. (default)
**    ="method2" ==> Direct (Method1) used for efield up until
**                   time step noncntrl.
**                   Then relax electric field so current converges towards
**                   specified time-dependent values of flux surface
**                   averaged parallel current (xjc,xjb,
**                   npwrxj,mpwrxj,totcrt).
**    ="method3" ==> Direct (Method1) used for efield up until 
**                   time step noncntrl
**                   Then relax electric field so current tends towards
**                   current obtained with specified elec. field up
**                   to the time step noncntrl.
**    ="method4" ==> relax electric field so current converges towards
**                   specified time-dependent values of parallel current
**                    (xjc,xjb,npwrxj,mpwrxj,totcrt).  (I.E. as in "method2").
**                   Initial electric field determined through spitzer
**                   conductivity.
**                   Presently setting this up for specified radial
**                   profile of parallel current at min B point.
**                   Could add specification of toroidal or poloidal
**                   current fairly easily (BH, 020227).
**    ="method5" ==> relax electric field so current converges towards
**                   EQDSK specified parallel current (mu_0*J_parallel=
**                   mu_0*R*(B_phi/|B|)*dp/dpsi+|B|*dF/dpsi).
**                   J_parallel is given the sign of I_plasma in the EQDSK.
**                   Initial electric field determined through spitzer
**                   conductivity. Should also set efflag="parallel";
**                   otherwise could have problems if B_phi changes
**                   sign as in an RFP.
**                   (Needs eqmod="enabled". Code checked for
**                   eqsource="eqdsk" setting, but not others.)
**    ="method6" ==> At t=0, save the value of elecfld and conductivity.
**                   At t>0, use spitzer+neoclassical conductivity
**                   to evolve electric field so that 
**                   E(t) =  E(t=0)* sigma(t=0)/sigma(t), 
**                   responding to change in temperature and/or Zeff,
**                   where sigma is based on 
**                   Sauter, Angioni and Lin-Liu, Phys.Plasmas 6, 1834 (1999) : 
**                   This method is aimed to yield j(t)=const, if no j_RE.
**                   For this method6, need to have iproelec='parabola' or 'spline'
**      
**                   (Added by YuP[2020-04-13]; see efield.f)
**
**  efswtchn="neo_hh" is switch to turn on an alternative calculation
**    of plasma current density from the code, to be used with above
**    convergence to prescribed values according to method2==>method5 
**    The calc'd current density, to be compared with the target
**    current density,  will be the sum of Hinton and Haseltine 
**    neoclassical contribution  (eta_neoclassical*elecfld) and runaway 
**    current from the electron distribution.
**    This option is provided in order to more accurately calculated current
**    in higher collisionality cases (such as in the cold plasma after
**    pellet injection or after disruption) in which the thermal particles
**    may be collisional but the tail electrons remain in the low-
**    collisionality banana regime.
**    (default="disabled").
**
**  curr_edge=An linear wedge of current varying from 0. at r/a=0.8
**            to curr_edge at r/a=1.0, which may be added to the
**            eqdsk obtained current density (efswtch="method5").
**            This is parallel current density at the min B flux surface pt.
**            Units are Amps/cm**2.  Default: curr_edge=0.
**            Might be useful for transport runs. 
**    
**  efiter="enabled" gives up to nefitera (a code parameter,  
**            presently 10) iterations for the electric field used 
**            AT EACH TIME STEP,  to achieve the target current within 
**            fractional accuracy currerr (nml below).  For
**            each iteration, the electric field is reset to
**            a new value according to the efrelax1, and
**            the Fokker-Planck equation for the time step is
**            re-solved.  "enabled" doesn't make sense for 
**            efswtch="method1".   (default="enabled".)
**
**  efrelax1=relaxation for efiter procedure. (default=0.5).
**                   
**  currerr=fractional accuracy for efiter="enabled" algorithm.
**          (default=0.1)
**
**  efrelax= relaxation parameter used in efswtch.ne."method1" or
**           efiter="enabled".
**           (efld ==> efld*( 1.-efrelax*(curr-currxj)/max(|curr|,|currxj|) ), 
**                     where currxj is the parallel flux surface area average 
**                     target current density). 
**                    ([2020-04-03] For "method4", this is generalized 
**                     using the new namelist variable efrelax_exp, see below.)
**            (See efield.f file).    
**            (default=0.5)
** 
**  efrelax_exp= power in a generalized procedure for relaxation of E field,
**            added in [2020-04-03].
**            The generalized procedure is 
**            elecfld==> elecfld*(1.-efrelax*sign(djrel*currxj)*
                                  |djrel|**efrelax_exp)
**            where djrel= (curr-currxj) / max(|curr|,|currxj|)
**            Presently this generalized procedure is only available for 
**            efswtch="method4", but could be extended to other "method#".
**            Default value is efrelax_exp=1.d0, which gives
**            the original procedure. (See efield.f file).
**            Recommended alternative value: efrelax_exp=0.5
**            which gives a faster convergence to the target current,
**            although sometimes produces "jiggles" in E(t).
**            
**
**  noncntrl= time step at which current control is turned on,
**            for efswtch.eq."method2" and "method3".
**            (default: noncntrl=0).
**
**  elecfld(l_) is the (Ohmic) electric field (volts/cm) vs flux surface
**              for machine="toroidal". Electric field values
**              are evaluated on each flux surface at the major radius
**              of the magnetic axis, and field varies on a flux
**              surface like 1/R.  That is E(R,Z)=elecfld(l_)*Rmag/R.
**    See input variable reden for lrz=1 and lrz > 1 cases.
**    Powers for "parabolic" profiles given by npwrelec,mpwrelec.
**    default: elecfld(0:1)=0.
**
**  efflag  ="toroidal", elecfld in code is in toroidal direction.
**          ="parallel", elecfld in code is parallel to magnetic field.
**           For "parallel" case, code is set up only for cqlpmod='disabled'.
**           default: "toroidal"          
**
**  iproelec="spline", then set elecfld using splines of
**            elecin(1:njena) (in manner of enein,...). 
**          ="parabola", set elecfld using elecfld(0),elecfld(1),
**           npwrelec,mpwrelec.               default: iproelec="parabola".
**          ="prbola-t" or "spline-t", time-dependent profiles
**
**  elecin(1:njena)=spline input at the ryain, of electric field (volts/cm),
**                  evaluated at major radius of the magnetic axis.
**  elecin_t(1:njena,1:nbctime)=spline input at the ryain (volts/cm) at
**    times bctime(1:nbctime), for iproelec="spline-t" 
**
**  elecscal=scale factor for elecin spline input (default=1.0).
**
**  npwrelec,mpwrelec are parameters of parabolic electric field with
**    iproelec="parabola", or time-dependent elecc,elecb above ("prbola-t").
**
**  eoved is the value of E/E-Dreicer. If eoved .le. 0 then the
**    d.c. electric field is determined from the value of
**    elecfld (above). If eoved .gt. 0, then elecfld is determined
**    from eoved by calculation. Used for runaway calculations
**    in general. Works for lrz .ge. 1.   If .gt.1, should use .ge.4,
**    since density/temperature are splined.
**    default:eoved=-.01
**    If lbdry(kelec).eq."fixed" and eoved.gt.0, then elecfld
**      is set through eoved only at time step n=0.  This avoids
**      variation of elecfld due to variation of E-Dreicer as
**      the electron density changes.
**    E-Driecer = m_e nu_0 v_Te / (2*charge), v_Te=sqrt(T_e/m_e),
**     nu_0 = 4 pi n_e charge**4 ln(lambda)/(m_e**2 v_Te**3)
**     (E-Dreicer as in Wiley and Choi, PF,p. 2193(1980) and Kulsrud et al.,
**        PRL (1973))
**
**  nonel and noffel are the on and off cycles for the electric field term
**    default: nonel=0; noffel=100000
**
**  Time-dependent profiles specified as above
**
**     
**********************************************************************
**
**  SIMULTANEOUS SOLUTION WITH AMPERE-FARADAY EQNS FOR TOROIDAL E-FLD
**
**********************************************************************
**
**    
**  ampfmod="enabled", for simultaneous soln of FP eqns and the
**                     Ampere-Faraday eqns for induced toroidal electric
**                     field in the torus, by an iterative technique
**                     using h and g functions as derived from 
**                     Kupfer,Harvey etal, PoP 1997.
**                     The initial electric field in the plasma is
**                     specified through methods specified through
**                     above D.C. ELECTRIC FIELD TERM (efswtch=method1 is
**                     the only valid choice at this time).  The Amp-Faraday
**                     solution is turned on at nonampf.
**                     General species must include electrons.
**                     (Presently, only commissioned for 1st gnrl species.)
**                     Presently, only implemented for soln_method='direct'
**                     'itsol', and 'itsol1'.  Implement for it3drv later.
**                     For the Amp-Far solution, use elecfld(lrz),
**                     elecin(), or elecb(time), depending on the value
**                     of iprone='parabola','spline', or 'spline-t'.
**                     For spline of spline-t, need 
**                     ryain(1)=0., ryain(njene)=1.
**                     Setup only for lrz=lrzmax (need smoothness).
**                     default: ampfmod="disabled".
**
**  ampfadd= "neo+bscd" [default value] Include neoclassical collisional 
**                     correction for current and also electron bootstrap current 
**                     into Amp-Far equation.
**                     The neoclassical collisional correction is calculated as
**                     j_coll_neo= (sigma_coll_neo - sigma_banana)*E_phi 
**                     where sigma_coll_neo is calculated from call_resthks
**                     with starnue>0 (calculated by sub.starnue_sptz), 
**                     while sigma_banana is calculated from call_resthks
**                     with starnue=0. This correction is useful for a collisional
**                     portion of plasma where the bounce-average approximation may fail.
**                     Other values for this variable:
**         ="add_bscd" Only add bootstrap current into Amp-Far equation,
**                     and skip the j_coll_neo correction. 
**         ="neosigma" Only add j_coll_neo correction, 
**                     and skip the bootstrap current. 
**         ="disabled" Skip both of these corrections.
**                     If "add_bscd" or "neo+bscd" option is selected
**                     then user should also set bootst="enabled"
**                     and specify jhirsh=99 (recommended) or jhirsh=88
**                     for the bootstrap model.  If user forgot to enable 
**                     the bootst="enabled" option in such a case,
**                     it will be reset in subr. ainsetva; also jhirsh=0
**                     will be reset to jhirsh=99.
**                     
**  nampfmax=max number of iterations to achieve given accuracy
**    across the torus (ampferr=max((E_l_iter+1-E_l_iter)/E_l_iter)),
**    where max is found for l=1,lrz at each iteration iter.
**    (default: nampfmax=3).
**
**  ampferr= iteration error tolerace (default=1.d-3)  [Presently 140114
**           remains to be implemented.  Use nampfmax.]
**
**  nonampf=turnon time-step for Ampere-Faraday solution.  That is,
**    the next time step, n=nonampf+1, toroidal electric field will be
**    calculated according to the Ampere-Faraday equations. 
**    (default: nonampf=0)
**
**  Note: Before the Ampere-Faraday eqns are applied, i.e, at time steps
**  earlier than nonampf, the electric field is calculated from 
**  settings according to iproelec and efswtch. 
**  Presently, the electric field before nonampf can be setup
**  as a fixed parabola (iproelec.eq."parabola"),
**  fixed spline (iproelec.eq."spline"),
**  or as a time-dependent parabola (iproelec.eq."prbola-t"),
**  with efswtch=method1, tmdmeth=method1. See subr. ampfefldb.
**  Note that for the Ampere-Faraday eqns, only the edge value
**  of electric field is needed. For the time-dependent option
**  (iproelec.eq."prbola-t") it is set through the value of elecb(1:nbctime).
**  The value of electric field at rho=0.0, 
**  which is set by namelist variable elecc(1:nbctime), 
**  is not used by Ampere-Faraday eqns, but it is used at n<nonampf
**  to setup the parabolic profile of elecfld(0:lrz) field,
**  and elecc(1:nbctime) should be ne.zero.
**  
**  
***************************************************************
**
**  TOROIDAL ROTATION:
**    Presently effects only the injection energy
**    of neutral beam particles.  (Subtracted off toroidal velocity
**    of FREYA injected particles *(R_birth/R_mag-axis), at birth point).
**    One use of this option: examine fusion neutron rate dependence
**                            on toroidal rotation rate.
**    
***************************************************************
**
**  nonvphi= time step at which toroidal rotation turned on
**           (default=10000, i.e., off)
**  noffvphi=time step for turn off.
**  iprovphi=indicates means for specifying tor rotation profile.
**           The profile is input as a flux function as for density,
**           temp, zeff.
**          ="parabola", then vphi=(vphiplin(0)-vphiplin(1))*
**                       (1-(rho/a)**npwrvphi)**mpwrvphi+vphiplin(1).
**          ="spline", then vphi given by splines of vphiplin(1:njene)
**                     at knots ryain(1:njene)
**          ="prbola-t", uses iprovphic(),iprovphib()
**          ="spline-t", uses vphiplin(1:njenea,1:nbctime) at bctime()
**            (default="disabled")
**  vphiplin(0:njenea)=gives tor vel profile (cm/sec) as specified above.
**                     (Typically 10**6-10**7  (10-100km/sec) in DIIID).
**  vphiscal= scale factor for vphiplin spline input (iprovphi="spline").
**  npwrvphi,mpwrvphi as specified above.
**
**********************************************************************
**
** PLOT, netCDF (AND PRINT) CONTROLS
**
**********************************************************************
**
**  Plotted output goes to file "mnemonic".ps
**
**  nchec is the # of time steps between sampling of data
**    for time plots (for example energy(k) vs. time).
**    default: nchec=1
**
**  pltd controls the plotting of the distributions:
**    pltd="df" plot f(n+1)-f(n); f(n+1)
**     (or "df_color" to plot in color) 
**    pltd="enabled"  plot f(n+1)
**     (or "color" to plot in color)
**    pltd="disabled" do nothing.
**    default: pltd="enabled"
**
**  pltlim,pltlimm control plotting of velocity dependent quantities
**    over a sub-portion of the mesh, and the form of the
**    independent (velocity-like) variable:
**    pltlim= "disabled",  plots versus x (u/vnorm) from
**                    x=0. to 1. (default:pltlim="disabled",pltlimm=1.)
**            "x",    plot 1d and 2d plots versus x 
**                    from 0. to pltlimm.
**            "u/c",  plot 1d and 2d plots versus u/c
**                    from 0. to pltlimm.
**            "energy", plot 1d and 2d plots verus energy (kev)
**                    from 0. to pltlimm (kev).
**
**  contrmin is the (normalized) smallest contour for contour plots
**    Contours span approximately contrmin orders of magnitude.
**    default: contrmin=1.e-12
**
**  constr is the (normalized) smallest streamline for the steady
**    state phase flow plot.
**    default:1.e-3  
**
**  ncont is the number of contours for contour plots
**    default: ncont=25
**
**    Other control of contour plots [see subr.pltcont(k,pltcase,tt_,itype)]
**    (these are hardwired, for now; see specific calls of pltcont())
**    Input itype, to identify what is plotted:
**    ! itype=1 for plots of f(),
**    !      =2 for df
**    !      =3 for source
**    !      =4 for urfb
**    !      =5 for vlfb
**    !      =6 for vlhb
**    !      =7 for rdcb
**    !      =8 for loss region (from pltlosc)
**    k is species index.
**    pltcase=1: geometric progression of plot contours,
**    pltcase>1: contours equispaced for max at temp(k,lr_).
**      (pltcase=2 is used for plots of f())
**    pltcase=0: equispaced, lin.scale from dmin to dmax [2023-03-09]
**    tt_ is contour plot heading.
** 
**  pltvecc="enabled" means vector flux plot for collisions is desired,
**    ="disabled" - not required. similarly for pltvece(electric field),
**    pltvecrf(RF field), pltvecal (sum of all fields).
**    It is necessary to use pltvecal=.ne."disabled", in order to get
**    the other vector plots in this category.
**    default: all "disabled"
**
**    veclnth=normalized maximum vector length for vector flux plots.
**      1.0=mesh scale size
**      default: veclnth=1.0
**
**    ipxy=number of perp velocities, for which vectors plotted.
**       (default: min(51,iy)).  Can be set lower.
**    jpxy=number of parallel velocities for which vectors are plotted.
**       (default: min(101,jx+1)).  Can be set lower.
**
**
**  pltend ="enabled" means time plots for density
**    energy, conservation constant, current as well as plots of
**    the local current J(v) and the local RF power prf(v) are
**    desired.
**    = "notplts" means time plots are excluded.
**    = "last" is like notplots if n .lt. nstop and is like "enabled"
**      if n=nstop.
**    = "tplts" means exclude J(v) and prf(v) plots.
**    = "disabled" means plot none of the above.
**    default: "enabled"
**
**  pltrst="enabled" means plot plasma resistivity and
**     related quantities vs. time if electrons are a general species.
**    default: pltrst="enabled"
**
**  pltvflu = "enabled" means plot normalized velocity flux (integrated
**    over theta) This gives runaway rates.
**    default: pltvflu="disabled"
**
**  pltra = "enabled" means plot runaway revelant output (<f>_theta) and
**             also PRINT some data to runaway.out.
**    default: pltra="disabled"
**
**  pltpowe = "enabled" means plot power vs time for physical energy
**    sources and sinks.
**    = "last" means do the plots only if n=nstop.
**    default: pltpowe="disabled"
**
**  pltfvs = "enabled" means plot f vs. v for some select values of theta
**           (at parallel, trapped-passing, perpendicular, anti-parallel,
**            and a pitch-angle averaged f with sin(theta0) factor).
**    default: pltfvs="disabled"
**
**  pltfofv = "enabled" means calculate f integrated over pitch angle, 
**    weighted by v_par*tau_bounce/int(ds/psi), versus u. (default="disabled")
**
**  pltprpp= "enabled" means plot reduced distribution f-parallel
**    default: pltprpp="disabled"
**  xprpmax=max perp velocity used in calc of f-parallel, normalized to vnorm.
**    (Also affects Besedin "knock-on" calculation; see below).
**    (Used for tandem.ne."enabled" and xprpmax.ne.1.0, only).
**    default: 1.
**
**  pltdn = "enabled" means plot density vs. poloidal angle for
**    various energy bins below (up to parameter negyrga bins)
**    default: pltdn="disabled"
**
**  [eegy(j=1:negyrg,1,k,lr_),eegy(j=1:negyrg,2,k,lr_)] define the 
**    energy bins (Kev) over which
**    the density is integrated for pltdn="enabled" (j'th bin
**    and k'th species) on the lr_'th radial flux surface
**    The bin associated with j=negyrga (negyrga=3 a parameter set in
**    source code) is reserved for energy bin [0,enorm].
**    For ease of input:
**    if lrzmax.gt.1 and the lr_=2 values are 0., then
**    the corresonding lr_=1 values will be duplicated to the
**    lr_=2,lrzmax values.
**  negyrg=number of energy bins.  If negyrg=0, no effect.
**    If negyrg=1, only obtain result for integration from 0 to enorm.
**    Evidently, if negyrg>1, then have to define eegy(,1:2,,) on all
**    FP'd flux surfaces.
**    (default: negyrg=0).
**
**  pltsig="enabled", then when sigmamod="enabled" give plots
**    of radial distribution of fusion power, and time-dependent
**    plots of total fusion reaction rates.
**    default: pltsig="enabled"
**
**  pltstrm = "enabled" - plot streamlines of the steady state flow.
**    default: pltstrm="disabled"
**
****Following pltflux option still needs conversion to PGPLOT****
****BH030829
**  pltflux = .ne."disabled" - plot velocity and theta components of the
**     momentum-space flux, for several component terms: 
**     pltflux1(1:7) = controls which fluxes are plotted:
**     pltflux1(1)=1. ==> sum of fluxes
**              2         collisions
**              3         parallel electric field
**              4         rf
**              5         synchrotron
**              6         Bremssstrahlung (+ phenomenological energy loss)
**              7         Krook operator slowing down
**     (default: pltflux1(1:6)=1., pltflux1(7)=0.)
**
**  pltinput="enabled" means print out the input values
**  pltinput="describe" means print out entire input-deck [defunct].
**    default: pltinput="enabled"
**
**  nrskip is a nonnegative integer. Given a Fokker-Planck'd flux surface
**    label "l", if ((l-1)/nrskip)*nrskip=l-1, then all the plots
**    described above controlled by the plt... namelist directives will
**    be plotted out for that flux surface. If nrskip=0, then the user
**    has the power to designate the flux surfaces for which the plots
**    are desired. See irzplt (below).
**    default: nrskip=5
**
** irzplt(i) is an array of length lrorsa. This variable is operative
**    if nrskip=0.       If irzplt(i)=l for some i, the FP'd
**    flux surface "l" will designate a flux surface for which
**    plots are desired, and the plt... directives above will be obeyed.
**    If irzplt(i)=0 for any i, irzplt(i) will be ignored.
**
**  plt3d="enabled" allows plotting of quantities (energy, current, etc)
**    as a function of normalized radius or normalized psi (see pltvs
**    below).
**    plt3d="enabled"
**
**  pltvs="rho" (normalized radius, see above)
**  pltvs="psi"
**    default: pltvs="rho"
**
**  pltrdc=.ne."disabled", then
**    rdcb diffusion coeffs are plotted for all flux surfaces.
**          ="one", or "enabled", plots first of nrdc diff coeff sets.
**          ="all", plot all diff coeff sets, 1:nrdc.
**          ="onecolor" Same as "one" but color plots.
**          ="allcolor" Same as "all" but color plots.
**          (default="all")
**
**  plturfb.ne."disabled" give plotting of QL B coefficient when
**          urfmod.ne."disabled" or rdcmod="enabled"
**          default:  plturfb="enabled"
**          ![2018-02-07] New: plturfb='color' for color contour plots.
**
**  pltrays='enabled'  (ACTUALLY, not presently namelist, but automatic
**                     when urfmod.eq.'enabled')
**          rays will be plotted over cross-sectional view (R,Z)
**          for urfmod.ne.'disabled' case.
**          default:  pltrays='enabled'
**
**  nplt3d(1:nplota) gives time steps for plotting of quantities
**      versus the radial variable.  
**      (Also for plots of UrfB if plturfb.ne."disabled".)
**    For such time steps,
**    and plt3d="enabled", the quantity/radius plots will appear.
**    default: nplt3d(1:nplota)=-10000
**    Also gives print out (tdoutput) for these steps.
**
**  idskf.ne."disabled", then at final time step,
**    print out meshes and distributions 
**    (((f(i,j,k,l),i=1,iy),j=1,jx),l=1,lrors), 
**    for each species k, to file named by contents of idskf.
**    idskf is declared character*8, so don't used longer filename.
**    ***Deactivated graph op as redundant, BH 000506***
*                        to file "graph  (see dsk_gr.f).
**    Also, print out more complete data to file idskf 
**    (see dskout.f): namelist, plasma profiles, meshes and
**    distributions. 
**    Note:  file name idskf should be quoted: "...".
**
**  idskrf.ne."disabled", then at final time step, 
**    print out namelist, meshes, and
**    urf quasilinear diffusion coefficients to file idskrf.
**    idskrf is declared character*8, so don't use longer filename.
**    Note:  file name idskfrf should be quoted: "...".
**
**
**  netcdfnm.ne."disabled", then create a netCDF output file
**    named according to the character*256 variable mnemonic//'.nc'
**    of code data (velocity and spatial arrays, equilibrium, 
**    plasma parameters, electric field, final distribution
**    function...,) to be used with auxiliary processing such as IDL
**    or python/matplotlib.
**    if urfmod.ne."disabled", then also write a netCDF file
**    mnemonic//'_rf.nc' containing ray data and the urfb QL coeff.
**    [But, see further qualifications according to netcdfshort.]
**    Note: // is the concatenation operator for character data.
**    (default: netcdfnm="disabled").
**
**  netcdfshort="enabled", flag to make shorter mnemonic//'.nc' 
**             and mnemonic//'_rf.nc'(if urfmod.ne."disabled")
**             and mnemonic//'_rdc.n.nc'(if urfmod.ne."disabled")  
**             netcdf files, omitting standard f(i,j,k,l) and
**             urfb/rdcb QL coeff at last time step.  
**             ="longer_f", i.e., not short, gives f at every
**                          time step (and urfb at last time step).
**                         NB: Presently (181023) code not set up for
**                             restart with this setting.
**             ="lngshrtf" (added in 2018) Save f at selected time steps
**                         only; The steps are set by specifying nplot()
**                         values in cqlinput.  The values of time
**                         are recorded into 'time_select' array 
**                         in the nc file.
**                         NB: Presently (181023) code not set up for
**                             restart with this setting.
**             ="longer_b", i.e., not short, give urfb coeff at
**                          every time step (and f at last time step).
**             ="long_jp", i.e., not short, give specific rf
**                          power pwrrf(u,r) and current currv(u,r)
**                          at each step (else only at last time step).
**                          YuP[2020-04-07]: Also added (for this option)
**                          saving favr_thet0(u,r) array, which is the value
**                          of INTEGR{f*2pi*sin(theta0)dtheta0} /4pi
**                          (Distr function at midplane averaged over pitch angle).
**             ='long_urf", urfa,urfb,urfc,urfd,urfe,urff into
**                          _rf.nc file at last time step;
**                          simularly for rdcb,rdcc,rdce,rdcf,
**                          into files _rdc.n.nc, n=1:nrdc.
**                          Otherwise, like netcdfshort="disabled".
**             default="disabled"
**
**  netcdfvecc=.ne."disabled" .and. .ne."x-theta", then at time step 
**    nstop, ouput to a netcdf file (mnemonic//'_flux_1.nc') the 
**    xpar,xperp components of the vector of velocity space fluxes 
**    due to collisions, on a normalized xpar,xperp momentum mesh 
**    also specified in the netcdf file.
**    (default: netcdfvecc="disabled")
**  netcdfvecc="x-theta", gives the u,theta-components of the velocity
**    space flux due to collisions on the code x0,theta0 grid
*     specified in the corresponding mnemonic//'.nc' file.
**    
**  netcdfvece, Similarly for electric field (mnemonic//'_flux_2.nc'),
**  netcdfvecrf, Similarly for RF field  (mnemonic//'_flux_3.nc'), 
**  netcdfvecal, Similarly for (sum of fields) (mnemonic//'_flux_4.nc').
**    default: all "disabled"
**  NOTE:  All four above netcdfvecXX setting should either by
**         "disabled" or have the same setting.
**  netcdfvecs="irzplt", then ouput netcdfvec data only at the
**                       radii designated by irzplt.
**            ="all",  then output netcdfvec data at all FP'd surfaces.
**            (default: "all")
**
**    ipxy=number of perpendicular velocities for netcdfvec o/p.
**       (default: min(51,iy)).
**    jpxy=number of parallel velocities.
**       (default: min(101,jx+1)).
**    ipxy and jpxy are the same variables as for pltvecal, etc.
**
**  f4d_out="enabled", then at time step nstop output the 2D velocity
**    distribution on a R,Z spatial grid, that is, output a 4D 
**    distribution function into a netCDF file, menmonic//'f4d.nc'.
**  f4d_out="tavg", then output time averaged 4D distribution
**    function.  Must also specify namelist tavg="enabled" and
**    associated time interval to avg over.
**
**    The f4d grids are presently specified as equispaced on the
**    following intervals:
**      Major radius grid from min to max points on outer plasma contour
**      Axial Z-grid from min to max points on outer plasma contour
**      Velocity grid from 0 to vnorm, normalized to vnorm, i.e., [0,1]
**      Pitch angle grid in [0,pi]
**  Grid dimensions are given by the namelist variables, respectively,
**  nr_f4d,nz_f4d,nv_f4d,nt_f4d.
**  The distributions are the LOCAL distributions, as calculated
**    according to the namelist ndeltarho=enabled.
**
**  tavg="enabled", then create a time avg distribution and output
**    it to the mnemonic.nc file in place of the final distn f (if
**    also have f4d_out="tavg").
**    Also, calculate NPA based on the avg distn, if f4d_out="tavg".
**    [default: "disabled"].
**
**  tavg1[1:ntavga]/tavg2[1:ntavga]= time interval (secs) of averaging,
**     for tavg=enabled. ntavga is parameter (presently =16) giving
**     the number of intervals over which the time avg is formed,
**     tavg(1) to tavg(2), etc., if not both zero. [default values: 0.d0]
**       
**

**  f3d_out="enabled", then at time step nstop
**    write 3D distribution, f3d==f(pol.angle,vpar,mu) at fixed rho,
**    to data file.
**    Original purpose: To be used by COGENT, as a boundary condition 
**    at specified rho (flux surface). Work with Ron Cohen, Mikhail Dorf.
**    See subr. f3dwrite for details [YuP August 2013].
**    In cqlinput, specify these values:
**    npol_f3d is poloidal-angle grid size,
**    nvpar_f3d, nmu_f3d are (vpar,mu) grid sizes
**    (All grids are equispaced).
**    f3d_rho is the coupling rho 
**    (in the same meaning as rovera(), 
**    and according to definition of radcoord ).
**    If not specified, the default values are
**    npol_f3d=20
**    nvpar_f3d=200
**    nmu_f3d=100 ! number of grid points in mu=0.5*m*Vperp^2/B
**    f3d_rho=0.5d0 ! f3d is saved at this rho value
**    The data is saved into file cql3d_f3d.txt, 
**    including explanation of units, in the file heading.
**    It will contain:
**    f(pol.angle,vpar,mu) at fixed rho==f3d_rho, 
**    also corresponding grids, and magnetic field B(pol.angle_grid).
**    Grids and distr.function are expressed  
**    in mks units and normalized by ref.values.
**    Example for magnetic field: [Tesla]/units_B
**    Example for f: [m^-3 (m/s)^-3]/ref_distfunc
**    Also, additional info is saved into cql3d_cogent_units.txt,
**    which is useful for setting an input file for COGENT.
**    It will contain the definition of 
**    units_B, ref_v, ref_mu, ref_distfunc, etc.
**    (These normalization units are "hard-wired" in subr.f3dwrite.)
**    Also, additional data is saved into cql3d_f3d_taub.txt,
**    which contains tau_bounce [sec] over (ipol,ivpar,imu) grid,
**    same grid that is used for the f3d distribution.
**    Default value: f3d_out="disabled"






***************************************************************
**
**  Calculation and output of first order orbit shift
**  
***************************************************************
**
**  ndeltarho=.ne."disabled" means calculate first order radial orbit
**             shifts for the first general species (ions of electrons)
**             from from given local R,Z,v,theta point to the
**             bounce-averaged position of the particle.  Output
**             is given at z-points along B-field measured from
**             the outer equatorial plane as a function of z,rho,theta0
**             (where theta0 is pitch angle at the midplane), and
**             also on a equispace R,Z,theta grid, where theta is local
**             pitch angle.  
**             These quantities with supporting data are output to the 
**             mnemonic.nc netCDF file (deltarho and deltarz arrays,
**             respectively).
**             Multiplying the arrays by v=v0 gives the shift
**             in the normalized minor radius coordinate
**             (as specified by the radcoord setting).
**             The R,Z grid is chosen about 10 percent larger than
**             the LCFS.  (default="disabled")
**
**             The ndeltarho feature is not setup to run with 
**             lrzdiff.ne."disabled".
**
**             ndeltarho can take on several values limiting the
**             application of the 1st order shift in cql3d to specified
**             portions of the code (for testing purposes):
**            ='enabled', then apply to available adjusted portions
**                of cql3d:  that is, urf diffusion, urfdamp1 ray
**                damping, freya NBI source.
**            [Following for testing purposes:]
**            ='urfb0', then just urf diffusion coeff calc
**            ='urfdamp', then just urfdamp of rays
**            ='freya', then just freya NBI source
**             
**             
**
**  nr_delta=number of R-points (default=65)
**  nz_delta=number of Z-points (default=65)
**  nt_delta=number of theta-points (Should be even, default=80)
**
**  A tri-cubic interpolation subroutine (deltar) is provided in
**  in file baviorbt.f, and can be used with the deltarz data,
**  and can be used as a template for adjusting external diagnostic
**  routines.
**
**
** 
**
***************************************************************
**
**  SOFT X-RAY DIAGNOSTIC
**  
***************************************************************
**
**  softxry="enabled" means that the diagnostic will be utilized.
**            Electron-ion and electron-electron collisions are
**            included, as described in the report CompX-02-2000.
**            Ions are treated as single species of charge Zeff,
**            density n_e.
**            X-ray spectra are added to the .nc file at the first
**            and last steps.
**          "ncdf_all" same as "enabled" except spectra for all
**            time steps is stored in the netCDF file (rather than
*             just the first and last time steps).
**          "e-ion", is a non-physical diagnostic mode which
**                   includes only e-ion collisons (no e-e).
**           default: softxry="disabled"
**
**  A right-hand cartesian x,y,z-coordinate system is assumed,
**  with the z-axis along the symmetry axis of the toroidal
**  device.
**
**  Toroidal geometry with circular or non-circular flux surfaces 
**  is assumed (depending on eqmod namelist variable).
**  The detector is assumed to be located outside the plama,
**    and without loss of generality, to be in the x-z 
**    plane (y=0.).
**  Detector location is specified by either
**    (a) rd(1:nv),thetd(1:nv): minor radius (cm) 
**                   measured from the magnetic 
**                   axis, and poloidal angle (degs) measured 
**                   from the x-axis in right-hand
**                   sense about the negative y-axis
**                   (defaults: rd(*)=100., thetd(*)=0.).
**       or, if the first elements of either of the following 
**       two inputs is .ne.0., use:
**    (b) x_sxr(1:nv),z_sxr(1:nv): the x,z-location (cm) of the detector.
**                   (defaults: x_sxr(:)=0., z_sxr(:)=0.)
**
**  nv=number of viewing cords.
**          Must be .le. nva (in param.h).
**
**  Detector viewing coords are specified by two angles,
**    (in the manner of the GENRAY and TORAY ray tracing codes):
**    thet1(1:nv)=polar angle measured from the z-axis (degs).
**                (90. deg is horizontal).  default: 90.
**    thet2(1:nv)=toroidal angle of the viewing coord, measured 
**          from the x-z plane (y=0.) in right-hand-sense with 
**          respect to the z-axis direction (degs).
**          (0. is in +R major radius dirn, 
**           180. is in -R major radius direction.)  default: 180.
**    We have detector direction vector x,y,z-components:
**          detec_x=sin(thet1)*cos(thet2)
**          detec_y=sin(thet1)*sin(thet2)
**          detec_z=cos(thet1)
**          
**  nen=number of energies at which spectral values calculated,
**          for each sightline.
**          Must be .le. nena (presently set to 30, in param.h).
**  enmin,enmax=minimum,maximum energy in the 
**              calculated spectra(keV). (defaults: 5., 50.)
**  msxr is the order of the highest polynomial used in the Legendre
**     expansion for computing SXR spectra 
**     (default value: msxr=mx in aindflt.f).
**     (Prior to 3/31/2011, msxr could not be greater than mx,
**      else reset to mx. Now, this restriction removed.  BH)
**     Values of .ge.8 are recommended.
**
**  fds= specifies the step size along the viewing cord, as a 
**     fraction of the average radial width of each radial bin 
**     (0.2 is a reasonable value).
**
***************************************************************
**
**  NEUTRAL PARTICLE (npa)  DIAGNOSTIC
**  
***************************************************************
**
**  npa_diag="enabled" means that the diagnostic will be utilized.
**           See Cql3d_npa_memo_1.pdf, CQL3D_npa_memo_2.pdf. These
**           reference the SXR description, TCV_sxr_CompX-2000-2.pdf.
**           See NPA_discussion_050826.eml for informal notes.
**           NPA spectra are added to the .nc file at the first
**           and last steps.
**          "ncdf_all" same as "enabled" except spectra for all
**            time steps is stored in the netCDF file (rather than
*             just the first and last time steps).
**           default: npa_diag="disabled"
**
**  NB: NPA presently assuming up-down symmetry in tdnpa0.
**      Fairly easy to generalize, as required [BH120328].
**
**  A right-hand cartesian x,y,z-coordinate system is assumed,
**  with the z-axis along the symmetry axis of the toroidal
**  device.
**
**  Toroidal geometry with circular or non-circular flux surfaces 
**  is assumed (depending on eqmod namelist variable).
**  The detector is assumed to be located outside the plasma,
**    and without loss of generality, to be in the x-z 
**    plane (y=0.).
**  Detector location is specified by either
**    (a) rd_npa(1:nv_npa),thetd_npa(1:nv_npa):  minor radius (cms)
**                   measured from the magnetic axis, and poloidal
**                   angle(degs) measured from the x-axis in right-hand
**                   sense about the negative y-axis
**                   (defaults: rd_npa(:)=100., thetd_npa(:)=0.).
**       or, if the first elements of either of the following 
**       two inputs is .ne.0., use:
**    (b) x_npa(1:nv_npa),z_npa(1:nv_npa): the x,z-location (cms) of 
**                   the detector. (defaults: x_npa(:)=0., z_npa(:)=0.)
**
**  nv_npa=number of viewing cords.
**          Must be .le. nva (in param.h).
**
**  Detector viewing coords are specified by two angles,
**    (in the manner of the GENRAY and TORAY ray tracing codes):
**    thet1_npa(1:nv_npa)=polar angle measured from the z-axis (degs).
**                (90. deg is horizontal).  default: 90.
**    thet2_npa(1:nv_npa)=toroidal angle of the viewing coord,
**          measured from the x-z plane (y=0.) in right-hand-sense
**          with respect to the z-axis direction (degs).
**          (0. is in +R major radius dirn, 
**           180. is in -R major radius direction.)  default: 180.
**    We have detector direction vector x,y,z-components:
**          detec_x=sin(thet1_npa)*cos(thet2_npa)
**          detec_y=sin(thet1_npa)*sin(thet2_npa)
**          detec_z=cos(thet1_npa)
**          
**  nen_npa=number of energies at which spectral values calculated,
**          for each sightline. (default =1)
**          Must be .le. nena (presently set to 30, in param.h).
**
**  enmin_npa,enmax_npa=minimum,maximum energy in the 
**              calculated spectra(keV). (defaults: 5., 50.)
**
**  m_npa [DEFUNCT] is the order of the highest polynomial used 
**     in the Legendre expansion for computing NPA spectra 
**     (.le. mxa in param.h).   default: m_npa=3
**     DEFUNCT [BH, 100521]:  No expansion of ion distribution is
**     used, since the CX cross-section is taken in the limit that
**     the ion energy is much greater than neutral energy, i.e.,
**     zero neutral temperature.
**
**  atten_npa="enabled", normal calculation of attenuation of 
**            charge exchange neutrals by reionization by e or i
**            impact.  [Could add other processes.]   (default).
**            "disabled", for numerical testing purposes.
**
**  fds_npa= specifies the step size along the viewing cord, as a 
**     fraction of the average radial width of each radial bin 
**     (0.2 is a reasonable value).
**
**  nnspec=fast ion species number, corresponding to the NPA neutral 
**         density species.  The NPA signal is derived
**         by neutralization (CX or recombination) from ion general 
**         species nnspec. Species are numbered as described above
**         for kspeci and related comments.
**         default:  nnspec=1
**
**  ipronn= "disabled", zero neutral density profile ("default")
**          "exp", specify normalized neutral density profile
**                   as exp( -(1.-rya(l))*radmin/ennl),  default
**          "spline", specify neutral density profile
**                   by values ennin(1:njene,:) at values of radius
**                   ryain(l)  (r/a)  l=1,njene (the x coordinate),
**                   as given above for enein.  njenn must be equal
**                   to njene, at this time, and is not a nml variable.
**          "spline-t", time-dependent profiles (see description of ennin_t)
**                   In this case, the profiles can be calc-ed 
*                    self-consistently by transport codes 
*                    and saved into data files (see read_ennprofile_times).
**      ipronn is a single, character*8 value.
**      ipronn has been extended to apply to up to npaproca NPA-related
**      processes [BH,20100720], below.  That is, it controls all
**      npaproc input density profiles. Presently, parameter npaproca=5.
**
**  npaproc=number of NPA processes to be considered.  Default=1
**          [If only want "cxh" and "radrecom" processes, still
**           need to set this number = 5]
**
**     Below, ennl,ennb,ennin,npa_process dimensioned to 
**     npaproc.le.npaproca.
**
**  npa_process(1:npaproc)=  Processes considered giving neutral
**                           NPA particles (set each of following in
**                           given array position): [default="notset"]
**                           (1)"cxh" gives CX of FI with hydrogenic
**                               neutral specified by 1st density 
**                               profile, ennl(1)/ennb(1) or ennin(:,1).
**                           (2)"cxb4" gives CX of FI with Boron+4 with
**                               2nd density profile, ennl(2)/ennb(2)
**                               or ennin(:,2).
**                           (3)"cxhe' CX of FI with He (not implmtd)
**                           (4)"cxc" CX of FI with C (not implmtd)
**                           (5)"radrecom" radiative recombination of
**                              FI with electrons (implmtd)
**                           [defaults: npa_process(1)="cxh", 
**                                   npa_process(2:npaproca)="notset"]
**                           Following densities are used, with storage
**                           in corresponding npaproc location, e.g., 
**                           ennin(*,1) for "cxh",
**                           ennin(*,5) for "radrecom" process.
**
**  ennl(1:npaproc)=scale length (cms), for ipronn="exp". 
**                  default: ennl()=5.
**  ennb(1:npaproc)=boundary value of neutral density (/cm**3), used for
**       normalization of above ipronn specified profiles.
**       Neutral density outside plasma may be taken equal to this.
**       default: ennb()=1.e10
**
**  ennin(1:njene,1:npaproc) = neutral density (/cm**3) profiles for
**       up to npaproca NPA-related processes.  Set values corresponding
**       to set values of npa_process().  default:  ennin(:,:)=0.0
**
**  ennscal(1:npaproc)= Scale factors for the above ennl/ennb or ennin
**       density profiles.  ennscal(i) scales ennl(i)/ennb(i), 
**       or ennin(*,i) profiles.  Default values = 1.0.
**
**  
** 
**
***************************************************************
**
**  FUSION ENERGY REACTION RATES and POWER (plus hoked up 
**    ionization and charge exchange).
**
***************************************************************
**
**  sigmamod="enabled", then calculate fusion rates between
**            relevant species, as indicated by isigmas below.
**            Reactions are between general species and itself
**            (for D+D), between general species and other Maxwellian
**            species, or between two general species.
**            If there is no general species, D, t or he3, then
**            sigmamod is reset to "disabled".
**
**  isigmas(1:6):  isigmas(1)=1, then compute d+t => alpha+n
**                        (2)=1               d+he3 => alpha+p  
**                        (3)=1               d+d => t+p
**                        (4)=1               d+d => he3+n
**                        (5)=1               hydrogenic impact ionization 
**                                            rate.
**                        (6)=1               charge exchange
**
**  mmsv= maximum order in Legendre expansion of reaction rate cross
**        section (set to mx in aindflt.f).
**     (Prior to 3/31/2011, mmsv could not be greater than mx,
**      else reset to mx. Now, this restriction removed.  BH)
**      [I don't find much sensitivity to N-dot for mmsv.gt.3. BH]
**
**  isigsgv1 = 0, Compute reactions of test (general) distribution with
**                self and others.
**             1, Also do same calc. except replacing the general species
**                with a  Maxwln with same energy and density
**                as test distn, for comparison purposes.
**             (default=0).
**
**  BH120314: Following option, isigsgv2=1 no effect on code, since
**  BH120314  coding is confused, and can't imagine a sensible application
**  BH120313: of isigsgv2=1.
**  xxxisigsgv2 = 0, Background Maxwellian distribution not included
**  xxx               in calculation of reactions,
**  xxx         = 1, Include background Maxwellian distn in calc of reactions.
**  xxx           (default=0). 
**
**  isigtst  = 1, then do sigmas(5)=1 calc using constant namelist
**                value for sigma*v = sigvi, and
**                do sigmas(6)=1 calc using constant namelist
**                value for sigma*v = sigvcx.
**  sigvi,sigvcx (cgs) as above.
**
**  Species in the code are recognized by the mass(with accuracy
**     ~1%) and charge:
**     H+ mass: 1.6726e-24 gr, bnumb=1.
**     D+ mass: 3.3436e-24 gr, bnumb=1.
**     T+ mass: 5.0074e-24 gr, bnumb=1.
**     He3++ mass: 5.0064e-24 gr, bnumb=2.
**     He4++ mass: 6.6442e-24 gr, bnumb=2.
**
**
***************************************************************
**
**  BOOTSTRAP CURRENT CALCULATION FOR NONTHERMAL DISTRIBUTIONS	
**
***************************************************************
**
**
**  bootcalc="disabled" (default) 
**          ="method1" or "method2", then compute jump conditions
**           at the trapped-passing boundary, simulating finite
**           banana width effects.  This gives a numerical 
**           calculation of bootstrap current. (Doesn't work if
**           implct.ne."enabled".or.lrz.eq.1)
**           "method1" uses radial expansion of local distribution
**             function (f=f_o-u_par/omega_c_pol*d(f)/dr)).
**           "method2" connects co-(counter-)passing region to trapped
**             particle distr. displaced half banana width
**             inwards (outwards).  (Presently assumes positive plasma
**             current.  9/93.  Needs more checking out, BH 990928)
**           See R.W.Harvey and G.Taylor, Electron Bernstein wave-bootstrap
**           current synergy in the National Spherical Torus Experiment", 
**           Physics of Plasmasa 12, 052509 (2005)
**           Only use lrzdiff.eq."disabled" with this option.
**           [YuP: Just to be clear, "method1" or "method2" produce 
**            a change in distribution function, so that the bootstrap
**            current is 'naturally' added into computation of curr() array
**            through this modification of the distribution function;
**            See functions bsl() and bsu() for details.
**            This is different from bootst="enabled" option, see below,
**            which does not produce a change in distribution function,
**            but only evaluates the bs current for general and Maxwellian
**            species based on analytical formulas; 
**            see subr.tdboothi for details, also arrays bscurm(lrz,ksp,ktype).
**            These bscurm arrays are saved into mnemonic.nc,
**            and also plotted in *.ps file (see subr.tdpltmne). Also, 
**            these arrays are used (optionally) for correction
**            of Ampere-Faraday solution, see ampfadd="add_bscd" or "neo+bscd".]
**  bootupdt="disabled" (default)
**          ="enabled" update f_o in bootstrap calculations beyond
**             initial set-up of the distribution at nonboot.
**  bootsign=+1.0 (default)
**          =a multiplier in bootstrap current calculation of the
**           jump in the distribution at trapped-passing boundary.
**           Positive bootsign gives contibution in the positive plasma
**           current direction for the usual negative radial 
**           plasma derivatives. (The direction of the plasma
**           current is derived from the eqdsk plasma current.
**           Positive is counter-clockwise viewed from above.)
**           (As of 990901, BH).
**  nonboot= time step n at which bootstrap calculation begins
**             (default=1, not 0, since method is explicit in time).
**
***************************************************************
**
**  CURRENT DRIVE DIAGNOSTICS
**
***************************************************************
**
**  bootst="enabled" means compute the bootstrap contribution to
**    the total driven current according to Hinton and Haseltine
**    formula, and add it in; It will be used
**    to update the "eqdsk" file if eqdskalt="enabled"
**    default: bootst="disabled"
**
**  jhirsh=0, use Hinton and Haseltine bootstrap model,
**            for total electron+ion bootstrap current,
**            good for all collisionalities, high aspect ratio,
**            simple electron-ion plasma.
**        =88,use Hirshman banana model, good for all aspect
**            ratios, low collisionality only. The simple electron-
**            ion plasma is extended to multiple ion species
**            using Zeff (of unknown validity). 
**            (S.P. Hirshman, Phys. Fl. 31, p. 3150-2 (1988).)
**        (default=0, to maintain backwards compatibility).
**
**        =99,use Sauter, Angioni, and Lin-Liu, PoP, 2834 (1999).
**            Corrected according to Erratum in Phys.Plasmas 2002,v.9,p.5140;
**            good now for any collisionality (coded by YuP 2019-12)
**            aspect ratios, general eqdsk data. For details,
**            see comments in subr. tdboothi.
**            This is the best jhirsh option.
**
**        jhirsh=88/99 give calc of bscurm for electrons and/or
**               ions,  based upon thermal or nonthermal 
**               energies.      
**
**  kpress(k)="enabled" means that species "k" is represented
**    in computing the total plasma pressure. This quantity is
**    used as output to the "eqdsk" file if eqdskalt="enabled".
**    default: kpress(k)="enabled" for k=1,..., ntotal
**
**************************************************************
**
**  INTERACTION WITH OTHER CODES
**
**************************************************************
**
**    partner="selene" means that CQL3D will call the SELENE 
**      JAERI MHD equilibrium code to generate new geometry.
**      This code has not yet been ported to UNICOS. Typically
**      this option has been used with the "FR" routines for
**      ITER type current drive. 
**    partner="bramb" means that CQL3D will dump an eqdsk for
**      use by a ray tracing code.
**    partner="lsc" and lh="enabled",  then the  sign of n_parallel 
**      is changed when input ray data is read in.
**    partner="disabled" None of the above.
**      default="disabled"
**
*******************************************************
**
**  ANALYTIC SOURCE CONTROLS
**  namelist sousetup
*******************************************************
**
**  k is species index, m is source index for each species.
**
**  nso is the maximum number of sources per species k: it is the maximum
**    m used in the arrays below. Note: nso .le. nsoa (parameter)
**    default: nso=0 (set in code), in which case no
**                   variables with index m below will be referenced.
**
**  nsou is the number of cycles (time steps) between recomputing the source
**    The source is recomputed if mod(n,nsou)=0.
**    Except:  Freya is only called from tdinitl, and re-calc'd if
**             (n+1).eq.nonvphi or noffvphi.
**             Thus, presently, if nsou.lt.nstop and n doesn't pass 
**             through the nonvphi/noffvphi interval, the freya source
**             will be turned off at n=nsou+1 when the source is reset.
**             N.B.[YuP]: So, for NBI to be intact, MUST have nsou>nstop,
**             otherwise it will be reset to zero at n=nsou+1, 2*nsou+1, ...
**             in subr. sourcee.
**    N.B.:  If nonso.ne.0 (or 1), then NEED nso=1, to get proper
**           turn on (and off, if noffso.lt.nstop).
**    default: nsou=1
**
**  nonso(k,m) is the on cycle for the source
**    default; nonso(k,m)=100000
**
**  noffso(k,m) is the off cycle for the source [default: 100000].
**  [Also, appears to be used for parallel transport adjustment:
**         see tdchief.f, call wptramu/wptrafx, for cqlpmod='enabled']
**  [YuP/2021: Not anymore, changed to local hard-wired variable n_sub_wp]
**
**  For cql3d runs (i.e., lrzmax.ge.4), the following namelist variables
**    are expanded parabolicly in the radial direction, using
**    szm1z(k,m,0) to szm1z(k,m,1), etc.  The powers npwrsou(k=1,ngen) 
**    and mpwrsou(k=1,ngen) (real) are used in the manner of reden for all
**    of the source specification variables (see above discussion),
**    except that npwrsou(0) and mpwrsou(0) are used for asorz.
**  But, if asorz(k,m,0)=-1., then use direct input of asorz (see below).  
**    
**  The following rank=3 source variables are dimensioned (ngena,nsoa,0:lrza).
**
**  szm1z(k,m,l) is the (z/zmax)-peak of local source(k,m)
**    default:szm1z(k,m,0)=szm1z(k,m,1)=0.
**
**  szm2z(k,m,l) is the (z/zmax) width of local source(k,m)
**    default: szm2z(k,m,0)=szm2z(k,m,1)=.78
**    Note: for details, see the end of subr.soup.
**    Peak position and dispersion are set through facz=exp(-(zl-zm1)**2/zm2)
**    where zl==z(l), zm1=zmax*szm1z, zm2=(zmax*szm2z)^2, 
**    smaller szm2z means stronger localization around z=zmax*szm1z
**
**  soucoord="cart" means cartesian coordinates are used
**    for input; soucoord="polar" is the other option (see below)
**    default: soucoord="cart"
**
**  The next eight energy quantities are in units of Kev.
**  The next four input variable arrays are used only for soucoord="polar"
**
**  sem1z(k,m,l) is the energy peak of source(k,m)
**
**  sem2z(k,m,l) is the energy dispersion of source(k,m)
**
**  sthm1z(k,m,l) is the pitch angle (deg) peak of source(k,m)
**
**  scm2z(k,m,l) is the cos(pitch angle) dispersion of source(k,m)
**
**  The next four variable arrays are used only for soucoord="cart"
**
**  sellm1z(k,m,l) is the parallel energy of the peak of source(k,m)
**    default: sellm1z(k,m,0)=sellm1z(k,m,1)=0.
**
**  seppm1z(k,m,l) is the perpendicular energy of the peak of source(k,m)
**    default: seppm1z(k,m,0)=seppm1z(k,m,1)=0.
**
**  sellm2z(k,m,l) is the parallel energy dispersion of source(k,m)
**    default: sellm2z(k,m,0)=sellm2z(k,m,1)=1.
**
**  seppm2z(k,m,l) is the perpendicular energy dispersion of source(k,m)
**    default: seppm2z(k,m,0)=seppm2z(k,m,1)=1.
**
**  asorz(k,m,l) specified m'th current for general species k
**    in particles/cc/sec
**    default: asorz(k,m,0)=asorz(k,m,1)=0.
**    Alternatively, if asorz(k,m,0)=-1., then input the source profile as 
**       asorz(k,m,1:lrzmax)
**    Note: In the above arrays, suffix 'z' is important.
**          See subr. tdstin which makes a copy of sellm1z() to sellm1(),
**          and other arrays with 'z' suffix.
**
**  pltso="enabled" means plot contours of the source.
**    ="first" plot contours only if n=nonso(1,1).
**    Also, options to plot with colors are added [2018-02-07]
**    ="color" means plot contours of the source in color,
**    ="first_cl" plot contours in color only if n=nonso(1,1).
**    Also plot source integrated over theta, in manner of pltfofv (above)
**    default: pltso="first"
**
*******************************************************
**
**  SOURCE CONTROLS FOR KNOCK-ONS (electrons)
**  namelist sousetup
*******************************************************
**
**  knockon.ne."disabled", gives source of "knock-on" electrons
 
**    ="fpld_dsk" also writes disk file "fpld_dsk" with source*dtr 
**                for test of source (see sourceko.f and finit.f),
**                then exits.
**    ="fpld_ds1" also writes disk file "fpld_dsk1" with source*dtr
**                PLUS f(i,j,1,1) (no-renormalization), 
**                for test of source (see sourceko.f and finit.f),
**                and another file "fpld_dsk2" containing only f,
**                then exits.
**  komodel="mr", Marshall Rosenbluth mode, with primary electrons
**                approximated by pitch-angle averaged distribution,
**                and source directly from MR expression, averaged
**                over pitch angle and interpolated on to code u,theta grid,
**                as suggested by SCChiu.
**                Presently, this is the only choice, so komodel has
**                no effect (bh: 980501).
**  nonko=the on cycle for the source.
**  noffko=the off cycle for the source.
**  isoucof=0(mode for lower cutoff of source, per bh), 
**          .ne.0 (based on ucrit, per sc). (default=0)
**    Per bh:
**      soffvte=if positive, a factor (times v_thermal_e) for velocity below 
**          which the knockon source is set to zero (to avoid double
**          counting of FP and knockon collisions).
**         =if negative, factor is abs(soffvte)*sqrt(E_c/E)*v_thermal_e,
**            where E_c=2*E_Dreicer (see Fuchs et al, Phys. Fl. 29, 2931 
**            (1986).
**          (default: soffvte=3.0)
**  isoucof==1:
**    Per sc:
**      soffvte=(.gt.0.)a minimum value for source cutoff and defn. of runaways
**            given as a factor of vthe(lr_)/vnorm.
**            (.le.0) then faccof*ucrit is used for cutoff. See sourceko.f.
**      faccof= factor times ucrit for cutoff. (default=1.)  (sc)
**      (Should combine bh and sc methods for cutoff).
**  soffpr=a factor which is multiplied by the above source cutoff velocity 
**         (from soffvte or faccof) giving absolute value of the parallel  
**         momentum-per-mass below which the primary electron distribution 
**         for the knockon process is set equal to zero. 
**c  nkorfn=.ne.1, refined the theta grid calc of the knockon source.
**c         It is the number of theta mesh points between the theta
**c         grid points at which the knockon source is calculated.
**c         (default: nkorfn=1).  Evidently, no longer used [BH071012].
**  flemodel="pol" or "fsa" to invoke poloidally dependent, or flux
**           surface average calc of the parallel distribution, 
**           respectively, for plotting purposes (pltprppr).
**           (For knockon.ne."disabled", only a pitch angle averaged
**           reduced distn is used, and plotting of this is presently
**           "disabled" (bh: 980501).
**  jfl=number of grid points in parallel (reduced) distribution function
**      (Hard-wired to jx, for knockon.ne."disabled")
**      
**  xlfac=parallel velocity mesh spacing factor
**    (Not used if knockon="enabled", in which case the parallel
**     mesh is same as the x-mesh.)
**    xlfac=1. ==> uniform mesh;
**    xlfac<1. ==> geometric mesh with greater mesh resolution at x=0;
**    xlfac>1. ==> greater resolution at xmax.
**    xlfac<0. ==> xpctmdl,xpctlwr,xmdl,xlwr must be input.
**    xlfac<0. ==> the momentum (velocity) mesh contains three regions:
**    The region +/-[0,xllwr] will have jfl*xlpctlwr evenly spaced 
**    mesh points.
**    The region +/-[xllwr,xlmdl] will have xlpctmdl*jfl mesh points.
**    The region +/-[xlmdl,xmax] will have the balance of the jfl points.
**    default: xlfac=1.
**   For knock-on source:
**    Need to have nso.ne.0 to turn it on.  nsou is as for the
**    analytic sources. Need soucoord="disabled" to turn
**    off above analytic sources, when nso.ne.0. 
**    pltso controls plotting.  
**    Also need nonso(1,1) and noffso(1,1) set so source 
**      goes on.
**
**************************************************************
**
**  MAGNETIC GEOMETRY INPUT FOR "EQ" ROUTINES
**  namelist eqsetup
**************************************************************
**
**  eqmod="enabled" means the code assumes that the flux surface geometry
**    is general, through eqsource . Integrations are then performed
**    to determine the flux surfaces and the local magnetic field. The
**    input variable psimodel is overwritten and set equal to "spline".
**    default: ="disabled".   (lrzmax can be .ge.1).
**    ="disabled", then analytic model available, specified by
**                 rboxdst rbox zbox radmaj radmin bth and btor
**
**  Presently (June, '02) the input equilibria are up-down symmetrized
**  in the code by several methods, as defined by eqsym, below.
**  These are all kluges.
**  The code is not presently set up for non-up/down-symmetric equilibria.
**  (But, implemented non-up-down symmetric eqdsk, BH, Oct'09).
**
**  eqsym= Specifies the manner of up/down symmetrization.
**          "none", to be used for no-symmetrization, 
**                  when it becomes available (BH:Oct'09).
**                  Only valid for eqmod.eq."enabled",eqsource="eqdsk".
**                  Must have taunew='enabled'.
**                  Note:  Still require in eqdsk that ymideqd, up-down 
**                         mid-point of the z-mesh, is equal to 0.
**          "average", up-down-symmetrization by averaging psi-values
**             above and below z=0 plane. 
**             (Only possibility until June,'02).
**          "avg_zmag", up-down-symmetrization by averaging psi-values
**             above and below z=zmag plane.  zmag is given in the
**             equilibrium data, the z-value of the magnetic axis.
**             Then vertically shift data so zmag is zero, and similarly
**             adjust RF and diagnostics (X-ray) and NBI source.
**          "top", up-down-symmetrization by reflecting psi in the plane
**             above zmag to z.lt.zmag.  Then vertically shift 
**             data so zmag is zero, and similarly
**             adjust RF and diagnostics (X-ray) and NBI source. 
**          "bottom", up-down-symmetrization by reflecting psi in the 
**             plane below zmag to z.gt.zmag.  Then vertically shift 
**             data so zmag is zero, and similarly
**             adjust RF and diagnostics (X-ray) and NBI source.
**       default: eqsym="average" 
**
**  eqsource controls the equilibrium disk file source
**    default: eqsource="eqdsk"  (when eqmod="enabled")
**  eqsource="eqdsk" indicates an eqdsk file is to be read in. An eqdsk
**    file is a standard equilibrium code output file. See, for instance,
**    GA Tech report GA-A18036 by Helton and Bernard. Input variables
**    rboxdst rbox zbox radmaj radmin bth and btor are reset from eqdsk
**    data. This file, which must reside on the user's disk space, is
**    named eqdsk by default but can be reset using eqdskin below.
**  eqsource="topeol"  indicates the existence of a disk file called
**    topeol which is the output from a Culham equilibrium code. This
**    is not a standard choice for eqsource.
**  eqsource="ellipse"  while eqmod="enabled" generates elliptical cross-
**    section flux surfaces; this option requires, in addition, that
**    that input variables ellptcty, eqmodel, eqpower, fpsimodl, and rmag
**    be set. The limiter position is assumed to be located at (R,Z)=
**    (rmag,radmin). radmin is the minor radius.
**    This is a fragile option, used only for the original debug
**    and is not recommended for use in general.
**    default: eqsource="eqdsk"
**  eqsource="tsc" implies that the code will utilize a file generated
**    by the Princeton (S. Jardin) TSC transport code...this file
**    provides more than magnetic geometry information. Density and
**    temperature profiles for species are also provided by this
**    file. This file is called 'tscinp'. (Must have lrzmax.gt.1).  
**  eqsource="miller"    !New option added in 2020
**    Analytically defined equilibrium, D-shape. 
**    REF: R.L. Miller et al., "Noncircular, finite aspect ratio, local
**    equilibrium model", Phys. Plasmas, Vol. 5, No. 4, April 1998,
**    See Eqs.36-37.  
**    All related namelist variables start with "eq_miller_" prefix.
**    Set in cqlinput: 
**     eq_miller_rmag= ! Magnetic axis: major radius coord [cm]
**     eq_miller_zmag= ! Magnetic axis: vertical coord [cm]
**     eq_miller_btor= !Tor field at Geom. center of LCFS [Gauss], with sign.
**       !Note: the sign of eq_miller_btor will be saved 
**       ! as bsign variable (stored in comm.h)
**     eq_miller_radmin= ! Plasma minor radius [cm]
**     eq_miller_cursign=![+1. or -1.] Sign of Plasma Current in phi-direction,  
**       ! will be saved and used as cursign variable. 
**     eq_miller_psimag= ! Pol.flux at magn.axis [cgs] Set as positive.
**     eq_miller_psilim= ! Pol.flux at LCFS: Set smaller than psimag.
**     eq_miller_psi_n=  ! n and m powers for PSI(r) profile as in 
**     eq_miller_psi_m=  !  PSI(r)= psilim+(psimag-psilim)*(1-(r/a)^n)^m
**     eq_miller_deltaedge= ! Triangularity of LCFS (at r=radmin)[e.g., 0.4]
**     eq_miller_kappa=  ! Vertical elongation (const for all surfaces).
**     eq_miller_drr0= ! dR0/dr [e.g., -0.3]. We assume Shafr.shift=-drr0*r 
**    The equations for the (Rs,Zs) coordinates that define flux surfaces
**    can be found in subr. eq_miller(), see these lines:
**    Rs= rmag + drr0*r + r*dcos(tpbs)  !we assume Shafr.shift=-drr0*r
**    Zs= zmag + akappa*r*sint ! Vertical coord [cm] along magnetic surface
**    Other related values:  
**    Set eqsym='none' ! The only option for "miller" equilibrium. 
**    Also use  eqmod='enabled' for "miller" equilibrium case.
**
**  For file-based eqsource, set:
**  
**  eqdskin= character*512 name of input eqdsk file.  Default="eqdsk".
**
**  bsign=Multiplier of toroidal B-field terms read from eqdsk(btor,f,ff'),
**    and of the total B-field plus sign of npar in URFREAD_ input files.  
**    This is a  fix of neg. B-field problems in the URF routine
**    and possible elsewhere.  See urfread_.f (BobH, 980611). 
**    This is only for eqsource="eqdsk", btor.lt.0.  (default: bsign=1.0)
**
**  eqdskalt="enabled" means that once the current is computed
**    it is utilized to compute p' and ff'. These quantities are
**    then used to update the "eqdsk" file if eqmod="eqdsk".
**    eqdskalt="disabled"
**
**  povdelp determines the flux surface at which the calculation
**    takes place if rovera < 0. The value of psi (poloidal flux
**    coordinate) chosen is psi=psimag-(psimag-psilim)*povdelp.
**    psimag is psi at the magnetic axis; psilim is psi at the limiter.
**    In this code psimag.gt.psilim.
**    0. .le. povdelp .le. 1.
**    default:povdelp=.2
**
**  ellptcty is the ellipticity of the contours for eqsource="ellipse".
**    =0 gives circles.
**    ---> 1 give ellipses increasingly high in the Z direction.
**      default: =0
**
**  rmag gives the R position of the magnetic axis for eqsource="ellipse"
**    default: 100. (cm)
**
**  eqmodel chooses a model for the value of the poloidal flux function
**    on the flux surface generated by use of option eqsource="ellipse":
**    psi(R,Z)=factor*E**eqpower where
**    E=sqrt(Z**2+(R-rmag)**2/(1.-ellptcty**2)).
**    The factor is determined from the choice of bth.
**    default: eqmodel="power"
**    No other models currently exist (4/15/88).
**
**  eqpower is described just above.
**    default: eqpower=1.
**
**  fpsimodl="constant" means the f=R*btor is independent of psi.
**    This is the only model currently available for eqsource="ellipse".
**
**  rbox and zbox define the R length and the Z height of the
**    computational domain where the equilibrium is determined and
**    psi(R,Z) is defined.
**    It need be set only if eqsource="ellipse"
**    default:rbox=zbox=100
**
**  rboxdst is the value of R defining the inside edge of the box
**    They need be set only if eqsource="ellipse"
**    (above).
**    default: rboxdst=50.
**
**  atol and rtol are the absolute and relative accuracies demanded
**    of the LSODE orbit integrating routine.
**    A coupled pair of O.D.E's are solved to determine the
**    field line orbit.
**    default; both are 1.e-8
**
**  For the LSODE O.D.E. solver we specify:
**   methflag=10 for nonstiff Adams method (no Jacobian used).
**           =21 for stiff (BDF) method, user supplied Jacobian.
**               (Jacobian supplied in subroutine eqjac)
**           =22 for stiff method, internally generated Jacobian.
**    default:10
**
**  nconteqn(integer), and
**  nconteq(character*8)  play roles 
**    in interfacing between the input magnetic
**    configuration and the chosen flux surfaces where the 
**    Fokker-Plancking will take place. In order to select the flux
**    surface psi associated with a choice of rovera, an interpolation
**    is done with respect to an array of flux surface values eqpsi(l).
**    If nconteq.ne."psigrid", then
**    (integer) nconteqn .le. nconteqa generates an array 
**    eqpsi(l) evenly spaced in psi.
**  nconteq="psigrid" generates an array evenly spaced at R values in
**    the the midplane corresponding to the eqdsk from the magnetic
**    axis to the plasma edge.  The value of
**    nconteqn is calculated.  (For a 33x65 grid, a typical number
**    is found to be 12).
**    Old defaults: "psigrid" (previously highly recommended operating mode,
**                         although not sure why, BH990909),
**                        BH131001: NOT RECOMMENDED.  The
**                        nconteqn value generated by this method
**                        is roughly the number of R points in the
**                        eqdsk between the mag axis and psilim,
**                        which may be much too small for accurate
**                        binning in nfreya.
**    NEW Recommendation (as of 131001): 
**                        nconteq="disabled", nconteqn=40 or 50.
**            
**    NEW defaults: (BH, as of 171211)
**      nconteq="disabled"
**      nconteqn=50  [If nconteq="psigrid", code will stop and user will
**                    be forced to reset
**                    "psigrid" --> "disabled", or reset nconteqn=0
**                    in cqlinput.]
**
**  lfield  is the number of poloidal
**    points at which the O.D.E. integrator returns values for
**    R(s) and Z(s), where s parameterizes the contour.
**    default: lfield=250 (note: lfielda is not used anymore).
**
**   END "EQ" MODULE INPUT
**
**************************************************************
**
**  RF  HEATING AND CURRENT DRIVE INPUT FOR VLH MODULE
**  namelist rfsetup
**************************************************************
**
**  vlhmod="enabled" means apply a simple (phenomenological) RF
**    parallel(Landau) or perpendicular diffusion model in which the 
**    magnitude of the diffusion and the resonance regions 
**    (up to nmodsa) are determined through input.  
**    (The model is separate from the rigorous QL/ray tracing   
**    model to be found in the urf module (urfmod.ne."disabled")
**    or the more comprehensive single flux
**    model accessible with vlfmod="enabled".) 
**    The model will be relativistic
**    if relativ="enabled"; non-relativistic otherwise.
**    default:vlhmod="disabled"
**
**  nonrf(1:ngen) is the cycle at which RF excitation begins
**  noffrf(1:ngen) is the cycle at which RF excitation ends
**
**  nrf which was originally set through parameter nrfa
**    can be reset through input. nrf=0 means no RF calculation
**    is to be done using the "rf"(presently null),  "vlf",
**    or ""vlh" modules,
**    If nrf = 0, rf excitation may still be present through the
**    involvement of the "urf" and "rdc" modules.
**    default:nrf=0
**    nrf.ge.1 gives number of consecutive general species 
**    (starting with 1) for contributions by these phenomenological 
**    diffusion coefficients.
**    nrf should be .le.ngena.
**
**  vlhmodes=number of separate resonance regions specified.
**    There may be up to nmodsa (a parameter) regions.
**    default=1.
**
**  vparmin(1:nmodsa) and vparmax(1:nmodsa) define the lower
**    and upper edges of the region in v-parallel space 
**    where the resonances occur (at R=radmaj, see vprprop below). 
**    They are in units of clight.
**    These specifications are velocity/c (NOT momentum-per-mass/c!)
**    (vparmin(i).le.vparmax(i), and specifications may be negative.)
**    npar=1/vparmin and 1/vparmax.
**    Defaults are 0.0.
**    For vlhmod="enabled" only.
**
**  vprpmin(1:nmodsa) and vprpmax(1:nmodsa) define lower and
**    upper limits in perpendicular VELOCITY, in units of clight,
**    of the QL diffusion.  defaults= 0. and 1.e100, respectively.
**
**  vlhpolmn() and vlhpolmx() define lower and upper limits
**    in poloidal angle, of the QL diffusion. (degrees).
**    default: 0. and 180., respectively.
**
**  vdalp is the width of the transition region
**    between zero and maximum diffusion. It is in units of
**    fraction of total resonance region. .03 is reasonable.
**    If vdalp is not used, numerical inconsistencies are
**    magnified.
**    default: vdalp=.03
**
**  vprprop="enabled" means that k-parallel is assumed to scale
**    with R for purposes of computing the resonance region for
**    vlhmod="enabled". ="disabled" means k-parallel assumed
**    constant.
**    default:vprprop="disabled"
**
**  dlndau(1:nmodsa) is the magnitude of the RF diffusion coefficient
**    if vlhmod="enabled". It is in units of vth^2/tauei. Typical
**    values would be 1. or 2.
**    default:dlndau(1:nmodsa)=1.
**
**  vlhprprp(1:nmodsa)="parallel" gives diffusion in u-parallel direction.
**           "perp" gives diffusion in u-perp direction.
**    default:vlhprprp(1:nmodsa)="parallel"
**
**  vlh_karney=Karney-Fisch type cutoff of Dlh with momentum,
**            =1/(1+u/pth)**vlh_karney factor, where 
**             u=momentum-per-mass, pth=(fmass(k)*temp(k,lr_)),
**             k is first general species. 
**             Karney and Fisch, Phys. Fluids 28, 116 (1985)
**             arbitrarily use vlh_karney=1.0.
**    default=0.
**
**  vlhplse.ne."disabled", then turn on square wave oscillating 
**    vlh diffusion with frequency and pulse length specified 
**    thru vlhpon and vlhpoff.  The frequency is 1./(vlhpon+vlhpoff),
**    and vlhpon gives the on-time during the pulse.
**    If vlhplse.eq."tauee", these times are measured in units
**    of tauee, else in terms of seconds.
**    These parameters are also applicable for vlhmod="enabled".
**
**  vlhpon  is the time the repetative pulse is on.
**
**  vlhpoff is the time the repetative pulse is off.
**
**************************************************************
**
**  RF HEATING AND CURRENT DRIVE INPUT FOR VLF MODULE
**  namelist rfsetup
**************************************************************
**
**  
**  vlfmod="enabled", then use multi-flux surface, multi-
**    harmonic, cyclotron  model for rf quaslinear diffusion,
**    Diffusion coefficient is according to Stix book (1993).
**    Region of wave interaction on each flux surface, 
**    polarizations, wavenumbers, harmonics, and frequencies
**    are  specified for up to nmodsa modes or harmonics (but
**    not both simultaneously).  The code calculates the
**    resulting bounce-averaged coefficients at each surface.
**    Can do eqmod="enabled"-cases, or eqmod="disabled" with
**    psimodel="spline" (but not presently "axitorus").
**    ONLY RELIABLE for electron cyclotron cases with
**    taunew="disabled".  (See taunew description, BobH, 010320).
**    
**  nrf, nonrf, and noffrf are used as above (or can use
**    vlhplse, vlhpon, vlhpoff, as above).
**
**  vlfmodes = number of wave-types of quasilinear excitation. (default=1.)
**
**  vlfbes = "enabled", use full Bessel functions in QL coeff.
**    (only possiblity at moment.)
**
**  vlfnpvar= "1/R" gives 1/R variation to parallel refractive index vlfnp,
**            "constant" gives constant vlfnp as function of position.
**
**  Each of following vlf-namelist variables is dimensioned 1:vlfmodes:
**
**  vlfharms()= number of harmonics.  (default=1.)
**
**  vlfdnorm()= Strength of QL diffusion coefficient, normalized
**    to collisional diffusion coefficient v_th(k)**2/tau_coll(e or i).
**    (See vlf.f).
**
**  vlffreq()= excitation frequency (Hz).
**
**  vlfharm1()= designates harmonic number of first of 
**    vlfnharms harmonics, or of each mode in the vlfmodes.gt.1
*     case (default=0.)
**
**  vlfnp()= central parallel refractive index at minimum B on the flux
**    surface, of a spectrum defined with vlfdnp and vlfddnp, below.
**    (Presently must be .lt.1. for relativistic case, 
**    for other than 0'th harmonics, 
**    i.e., no relativistic anomolous doppler cases (hyperbolic
**    resonance) at this time.  (default=0.5)
**
**  vlfdnp()= full width of par. ref. index, at full power. (default=.2)
**
**  vlfddnp()= additional region of par. ref. index, in which the
**    QL coeff. tapers to 0, like 0.5*cos((npar-0.5*(vlfnp+vlfdnp))*pi/vlfdnnp)
**    (default=.1)
**
**  vlfnperp()= perpendicular refractive index.
**
**  vlfeplus()= complex E_plus/abs(E) polarization. (default=0.,0.)
**
**  vlfemin() = complex E_minus/abs(E) polarizaion. (default=0.,0.)
**    (The parallel rf polarization is 
**    sqrt(1.-cabs(vlfeplus**2)-cabs(vlfemin**2))
**
**  vlfpol() = center of poloidal region on flux surface with
**    non-zero QL coeff. (degrees). (default=0.)
**
**  vlfdpol() = full width of poloidal region with full QL coeff.
**    (360. gives full flux surface). (default=360.)
**
**  vlfddpol()= taper distance of poloidal QL coeff, analagous to
**    vlfddnp above (degrees). (default=20.)
**
**  vlfparmn()=An additional window is put on the quasilinear
**    diffusion coefficient, for testing which velocity region
**    of the diffusion coefficient causes what current.  This
**    is the minimum parallel momemtum, units of vnorm.
**    default=-1.e100
**  vlfparmx()=Max parallel momentum. default=+1.e100
**  vlfpermn()=Min perpendicular momentum. default=0.
**  vlfpermx()=Max perpendicular momentum. default=+1.e100.
**   
**   
**
**************************************************************
**
**  RF HEATING AND CURRENT DRIVE INPUT FROM INPUT DIFF COEFFS
**  namelist rfsetup
**************************************************************
**
**  
**  rdcmod="format1" or "aorsa", then read in externally calculated
**    velocity space diffusion coefficient D_u0u0 in cgs units.  This
**    is a bounce-averaged coeff giving diffusion in 
**    momentum-per-mass coordinates: u0=x*vnorm, where x is
**    the code normalized momentum-per-mass coordinate and pitch
**    angle theta0.  The u0,theta0 are at minimum-|B| on a flux 
**    surface.
**
**    "format1" is data formatted according to AORSA convention:
**    Data is given on a regular grid of u_par0 and u_perp0, 
**    with upar0: -unorm,+unorm, uperp0: 0,unorm.
**    (Velocity unorm, here, is specified in the file, and may
**    be different from cql3d vnorm (specified through enorm).
**    (for relativistic cases, unorm refers to momentum-per-rest-mass.)
**    If unorm in the file is less than cql3d vnorm, then the
**    diffusion coeffs above unorm will be set =0. on the cql3d grid.) 
**    Data is read in from file du0u0_input, with format shown
**    in subroutine rd_diff_coeff.
**
**    The read in data is interpolated onto the code
**    momentum-per-mass grid.
**    ===>  Presently, the rdcmod radial grid IS TAKEN TO   <====
**    ===>  BE THE SAME as specified rya(1:lrz) for cql3d.  <====
**    Except, if data is given on twice as many flux surfaces as
**    specified by lrz, then the data set is reduced by omitting
**    data for every second flux surface.
**
**   Multi-file "aorsa"/"format1" set up. Can input multiple diffusion
**     coeff files and ql diffuse separate or
**     the same general species depending on setting of
**     nrdcspecies(1:nrdc).  [BH110608]
**    
**   nrdc=number of "aorsa"/"format1" files to be read (must be .le.
**     parameter nrdca, presently set = 10). Default nrdc=1 (gives backwards 
**     compatibility); that is, the default input file is named du0u0_input.
**     It is assumed that the velocity and radial grids 
**     are the same in each of the input diff coeff files.
**   
**   rdcfile(1:nrdc)=file path for each file (up to 256 characters).
**             default: rdcfile(1)=du0u0_input, rdcfile(2:nrdca)="notset"
**   
**   nrdcspecies(1:nrdc)=general species index to which the respective
**     diffusion coefficient files are applied.
**     Default: nrdcspecies(1:nrdca)=1
**   
**   rdcscale(1:nrdc)=scaling of each of the diffusion coefficient files.
**     Default: rdcscale(1:nrdca)=1.
**           If nrdc=1 and namelist rdcscale(1).ne.1., then use this value.
**           However, if nrdc=1 and rdcscale(1).eq.1.(e.g., the default)
**           then set rdcscale(1)=pwrscale(1).
**           (This maintains backwards compatibility.)
**
**   rdcmod="format2" reads in velocity space diffusion coefficient 
**    from files du0u0_grid and du0u0_rnnn, where nnn=[001-n_psi], 
**    and n_psi is number of radii given in du0u0_grid.  Thus,
**    the set of velocity-space diffusion coeffs are given in a
**    separate file for each radius, as provided by the DC diffusion
**    coefficient code.   Formats are otherwise
**    similar to "format1".  Note restriction re rya()!
**    Presently not set up for multiple radial sets of diffusion coefficient
**    files (as in rdcmod="format1" case, but this feature can readily be 
**    added if there is a need,  BH180507).
**
**    default:  rdcmod="disabled"
**
**  plturfb (in namelist section setup, see plot controls above) gives
**    plots of urfb diff coeff.
**
**  rdc_netcdf=ne."disabled", then
**    rdcb,etc.,  written into a netcdf file mnemonic//'_nn_rdc.nc', 
**    where nn=number of diffusion coefficient input file.
**          ="one", or "enabled", write first set of diff coeff
**            to mnemonic//'_rdc.1.nc'
**          ="all", write mnemonic//'_rdc.n.nc', n=1:nrdc
**          (default="all")
**    The rdcb, rdcc, rdce, rdcf are on the cql3d vel space grid.
**    Amount of output affected according to netcdfshort='long_urf'.
**
**  pwrscale(1) [from urf section of code, presently used to
**               scale input value of diffusion coeffs, but
**               can be superceded by rdcscale(1),
**               (above) for rdcmod="format1"]
**
**  rdc_upar_sign=+1. to maintain u_par sign convention of du0u0_input,
**                -1., to reverse order.   This should be used for DC
**                when the toroidal mag field is negative(clockwise),
**                since upar in DC is pos in dirn of B-field,
**                but in cql3d, upar is positive when in toroidal dirn.
**                (default=+1.)
**
**  rdc_clipping='enabled' for clipping of diffusion coefficient spikes
**               which sometimes occur near the trapped-transiting bndry.
**               Clipping is according to "Running Median Filters and 
**               a General Despiker", John R. Evans, Bull. Seismological
**               Soc. Amer. (1982).  Efficacy of the method has been
**		 investigated with python plotting script
**               plot_DC_multiR_format_cqlf.py.
**               default='disabled', for backwards compatibility.
**
************************************************************************
**
**    BEGIN URF MODULE INPUT
**    namelist rfsetup
**
************************************************************************
**
**  urfmod.ne."disabled" means utilize the module for LH/EC/FW  excitation,
**    utilizing data ray tracing.  EC includes EBW. Must have iy.le.255.
**    default: urfmod="disabled"
**
**  Diffusion coefficents for up to nmodsa (parameter, presently = 155) 
**    wave types(rftype, defn'd below) may be simultaneously treated 
**    with damping from several harmonics from each wave type.  The sum
**    of the number of harmonics (including the 0th) of the wave types
**    must be .le.nmodsa.
**
**  In the cql3d source (not namelist, but for programming considerations):
**    mrf= number of rf "types" (and is determined from rftype(), below).
**         Each rf type is read in from one of mrf RF data files.
**         Different frequency waves must be read in as separate types. 
**    mrfn= number of rf "modes" (sum over mrf of the nharms())
**          It is required that mrfn.le.nmodsa, nmodsa being a 
**          parameter set in param.h.
**    krfn(1:mrfn)=  wave type index (in 1:mrf) for each rf mode
**          Thus nrfspecies(krfn(1:mrfn)) gives species index to which
**          a mode is applied.  See nrfspecies() nml input below.
**
**  The wave types (old, less general, deprecated input method) are 
**    indicated by:
**    ech="enabled" (default="disabled")
**    fw="enabled"   (default="disabled")
**    lh="enabled"   (default="disabled")
**    (If none of lh, ech, or fw are enabled, then urfmod will be
**     reset to "disabled", except as more recently specified in the
**     following paragraph.)  Alternatively, if urfmod="disabled" then
**     calculations as in this section are turned off.)
**     Each wave type can have a different wave frequency, however
**     can have only one wave frequency for each wave type.
**
**  Alternatively, wave types can be input as rftype(1:mrf).
**    [This is the more general, preferred input method for wave types.]
**    Wave ray tracing data for each wave type is specified through a
**    file format designator rfread, and file name rffile(1:mrf),
**    as given below.
**    The method can be used for inputing ray data for the same wave
**    mode but at a different frequency, using a separate rffile.
**    [Only put ray data for one wave frequency in each rffile.]
**    The rftype(1:mrf) namelist input is a generalization of the above
**    ech/fw/lh method;  the older method is retained for backward
**    compatibility of namelist input.
**  rftype(i,i=1:mrf): Input values for each i, are "ec",
**    "fw", or "lh"  [==>> NOTE: not "ech" <<==].
**    Only set consecutive values to be used, starting at rftype(1).
**    Default of rftype(1:nmodsa)="notset".
**    Note: different values of i=1,mrf can have same designator rftype.
**    This could be used for two fw sets of data at different frequencies.
**    Each frequency would have its own rffile() ray data set.
**
**  The code will not handle mixed use of ech/fw/lh and rftype() type of
**    input.  That is, either use the ech/fw/lh OR the rftype()
**    specification with ech/fw/lh='disabled'.
**
**  Ray data for the above wave types is input in files,
**    either in text or netcdf format, as specified by variable rfread:
**  rfread="text", then all input ray data from ascii files named 
**                rayech, if ech="enabled"
**                rayfw,  if fw="enabled"
**                raylh,  if lh="enabled"
**        ="netcdf", then there are more options for naming netcdf input
**                  files,a s given by rffile(1:mrf)
**        (default="text" for backward compatibility, but
**         more usual since about 2000 is "netcdf")
**  rffile(1:mrf)=names of input netcdf files (e.g., from genray or toray)
**    Default values for up to first three of rffile(1:nmodsa)=
**    "rayech.nc","rayfw.nc","raylh.nc".
**    If rrfile(1)="mnemonic", then the value of the &setup0 namelist 
**                 variable mnemonic is used to generate the rffile() names:
**                 rffile(1)=mnemonic//"_rf.nc"  [//=concatenate sign]
**                       (2)=mnemonic//"_rf.1.nc"
**                       (3)=mnemonic//"_rf.2.nc"
**                       ETC.
**                 This option is designed so that ray data will be
**                 output into the same file as it is read in from.
**                 The output file will contain updated data, and
**                 file proliferation will be reduced.
**    If rffile(1)=any name other than "mnemonic", then that name is
**                 used for the netcdf input file.  Set the corresponding
**                 rffile(2:mrf) if .gt.1 wave types are used.
**                 [User has to ensure that data files correspond
**                  to the order of the enabled subset of the rffile()=
**                  ech,fw,lh namelist variables above; for example,
**                  if only fw and lh are enabled, rffile(1) would
**                  give fw data, rffile(2) would give lh data.]
**
**  Note:  it is possible to apply the identical ray data files to
**         separate plasma species (for ngen.ge.2 cases), by entering
**         the identical path specs for sucessive rffile() entries.
**         The code detects the files have identical path specs and
**         directs the wave interaction to be with specified species
**         as given below in corresponding nrfspecies(1:mrf).
**         Wave damping for this repeated rffile data is computed by
**         summing both over the specified harmonics and over the species.
**
**  nharms(i=1,mrf).gt.0, then for each ray type calculate damping 
**            for harmonics nharm1(i) to nharm1(i)+(nharms(i)-1),
**            rather than according to nharm in the ray data file.
**            This is preferred input method.  
**        =0, then use harmonic number given in raylh, rayech, 
**            or rayfw, etc., files. The number of harmonics will be 1.
**            This input method is depracated.
**
**  nharm1(i=1,mrf)=lowest harmonic in the series for each wave type
**                  (default=0).
**
**  nrfspecies(i=1,mrf)=general species which each wave type is
**                      applied to.  default: nrfspecies(1:nmodsa)=1
**
**  Following call_ech/call_lh/call_fw not presently fully implemented.
**    [BH060314].
**  call_lh="enabled" means that CQL3D will call and utilize the Brambilla
**    ray tracing code (xbr).  ="disabled" means that the ray tracing code
**    will not be called, but that CQL3D will still expect the presence
**    of a raylh file (output from Brambilla code)
**    default: call_lh="disabled"
**
**  call_ech="enabled", call and utilize Toray code (toray). 
**          ="disabled", then still may get data from rayech file.
**
**  call_fw="enabled", call and utilize Brambilla code (xbr) for
**    fast waves. (Since polarizations come from the ray data file
**    this option can also be used for lh.
**
**  plturfb (described under plot controls above) gives plots 
**    of urfb diff coeff.
**
**  pwrscale(1:mrf)=scale factor to be applied to the power entering
**    the plasma as given in the ray data files, applied to each of
**    wave types used.
**
**  pwrscale1(),urftime(nurftime),nurftime:
**  Time-dependent scaling of power is achieved with an additional
**    multiplier of the pwrscale(1:mrf), above. This is implemented 
**    by specifying the array variable pwrscale1() at a sequence
**    of nurftime corresponding times, specified through in the 
**    variable urftime() (seconds).
**    Interpolation of pwrscale1 between time points is linear.   
**    If simulation time becomes greater than urftime(nurftime), 
**    then power scaling continues at the value pwrscale1(nurftime).
**    nurftime =0 gives no pwrscale1 scaling.
**    (nurftime must be .le. the parameter nbctimea, presently =101).
**
**    nurftime=0 is default.
**    pwrscale1(1:nbctimea) defaulted to 1.
**    urftime(1:nbctimea) defaulted to 0.
**
**
**  Setting central values =0. at the first time step
**    turns off profile generation through these namelist variable.
**    (default is all central and edge values =0.)
**
**
**  urfmult= multiplier of quasilinear diffusion coefficients
**           and damping calc in urf calcs.  (for test purposes).
**           default: urfmult=1.0
**
**  wdscale(1:nmodsa)=scale factor for n_parallel-width of rays,
**    applied to wdnpar, when read in from ray data file.
**
**  N.B.: Previous 2 variables, pwrscale and wdscale, will modify
**  ray data if it is output, as specified below by urfwrray="enabled". 
**
**  nbssltbl=the number of elements initially in the Bessel table.
**
**  nondamp= damping turned on in urf module, for n.ge.nondamp.
**           (The QL diffusion coeffs are calculated and used in the FP
**            equation, but the rays are not damped, until n.ge.nondamp.
**            This has been used for certain code diagnostic cases.)
**    Default: nondamp=0
**
**
**  The following namelist control variables are used to determine the
**    sequence of operations in the urf module:
**    nrfstep1, nrfstep2, nrfpwr, nrfitr1, nrfitr2, nrfitr3, urfncoef.
**    The purpose of these controls is to permit a "soft" turn-on
**    of the RF, avoiding formation of pathalogical distributions
**    during startup of the RF. In many situations, it is not
**    necessary to tune this process: just turn on full power.
**
**   nurf=A counter in the urf control subroutine urfchief,
**   nurf=the number of calls to urfchief which have resulted in
**        calculation or recalculation of the diffusion coefficients.
**        This variable will be incremented each time
**        n/integer(urfncoef*ncoef)*integer(urfncoef*ncoef).eq.n.
**        Generally, there is a solution of the FP eqn after each
**        call to urfchief.
**                                                  
**   The sequence of actions (as a function of 
**   increasing nurf.ge.0 at each call) is given by the following steps
**   (after each diffusion coeff calc, control returns to the calling
**    subroutine):
**   1. For nurf=0, calc or read ray data for nrfstep1 spatial steps 
**      along the ray [nrfstep1 presently is not used, i.e., read all data].
**   2. Calc. damping of ray data, and then resulting quasilinear
**      diffusion coeffs, using a fraction of the input power
**      = (1/2)**nrfpwr.
**      return.
**   3. Repeat step 2 for next nrfpwr calls, but with fractional input
**      power (1/2)**(nrfpwr-1), (1/2)**(nrfpwr-2),.... (1/2)**0. 
**      (This step is a no-op if nrfpwr=0). 
**   4. Iterate step 2 with full input power for next nrfitr1 calls.
**      (This step is a no-op if nrfitr1=0.  Usually want .gt.(nstop-nonrf)). 
**   5. Extend extendable rays by nrfstep2 steps.
**   6. Re-calc damping from ray data and then quasilinear diffusion 
**      coeffs.  Iterate this step nrfitr2 additional calls.
**   7. Steps 5 and 6 are carried out nrfitr3 times. 
** 
**   Thus choose 
**        nstop= 
**        (nrfpwr+1+nrfitr1+nrfitr3*(nrfitr2+1))*integer[urfncoef*ncoef]
**   if the above sequence is to be completed.
**   Default values:  [nrfstep1(:)=100 but presently not used], 
**                    nrfstep2=50
**                    nrfpwr=3
**                    nrfitr1=1, nrfitr2=1, nrfitr3=2
**                    urfncoef=1.0                             
**  
**
**  scaleurf="enabled" means rescale the contribution to urfb (the 
**      diffusion coefficient so that a particular ray does not 
**      "overdamp" on a  given flux surface volume. Note in the limit 
**      that the number of flux surfaces goes to infinity, 
**      overdamping would not happen. Here by invoking 
**      this option, we seek to override the possibility that
**      more power may be deposited by a ray than is actually in it,
**      due to the coarse grid. We recommend using this option.
**      default: scaleurf="enabled"
**
**  iurfcoll(1:mrf)="enabled" means add the collisional absorption of
**            the ray ....this information is passed in the ray data file
**            in variable salphac.
**       ="damp_out", then the damping coefficient along the ray 
**                    in the poloidal plane is written into the salphac
**                    (collisional absorption coeff) in the RF netcdf 
**                    file and (if urfwrray="enabled", 
**                    into the rayXXX ray data file).
**                    This is for code diagnostics (comparison of ray
**                    tracing damping coef with cql3d damping coeff).
**        NOTE:  In cases where the SAME wave is applied to
**               multiple general species (with nrfspecies, below),
**               usually will want this variable "enabled" only once.
**        default(1:mrf):"disabled"
**
**  iurfl(1:mrf)="enabled", means include the additional linear absorption
**                    from the ray data file [in ray data variable salphal]
**                    is added to the calculated damping or power
**                    flowing along the ray.  
**                    This may be used, for example, to add linear ion 
**                    damping calculated in the ray tracing code to 
**                    cql3d electron damping calculated in a cql3d
**                    electron simulation, or visa versa.
**        NOTE: In present GENRAY setup, it is always the electron linear
**              damping, see this line in GENRAY:
**              salphal(is)=2.d0*ckvipl_e ! electron damping coefficient.
**              It can be changed in future, to add linear damping on ions.
**        NOTE:  In cases where the SAME wave is applied to
**               multiple general species (with nrfspecies, below),
**               usually will want this variable "enabled" only once.
**        default(1:mrf):"disabled"
**
**
**   ieqbrurf designates the source of equilibrium data to be used by xbr.
**       Appropriate values are: 
**       ieqbrulh=0, to use Brambilla analytic "equilibria",
**               =3, to use standard eqdsk input.
**               =4, to use extended eqdsk, including density, and
**                   temperature profiles,...
**       If eqsource="tsc", ieqbrulh is reset to = 4.
**
**  urfdmp="secondd" means utilize the "second derivative" damping
**      diagnostic in computing the damping due to each ray as
**      it moves through the plasma. If urfdmp .ne. "secondd" the 
**      diagnostic uses an integration by parts technique to compute
**      the damping. We highly recommend "secondd" because convergence
**      of the FP solution for distn F and the QL modification of F
**      is determined by agreement between the sum of the damping of 
**      all rays and the absorbed power as computed from dF/dt. This 
**      latter diagnostic utilizes the "second derivative" approach so
**      consistency demands "secondd" for the rays.
**        default:"secondd"
**
**      If ndeltarho='enabled' then "secondd" option not implemented,
**        only urfdmp='firstd', calculating damping including an 
**        integration by parts removing the 2nd derivative in the
**        calculation with (df/dt)_ql, is implemented.  In a test
**        case with ndeltarho='disabled', the urfdmp='firstd' and
**        'second' results were quite close [BH110617]. 
**
**  urfrstrt="enabled"  => Option to "do not update delpwr power flowing 
**     along the ray".
**     For convergence studies with previously calculated ray data files...
**     (default="disabled")
**
**  urfwrray="enabled", then re-write ray data into ray data file
**            at end of run.   (default="disabled")
**            [BH120325:  Presently only able to write to the rayech,
**                        raylh, or rayfw text files.  Codes needs
**                        updating for write to .nc files and for new
**                        rftype()/rffile() data specification system.
**          ="krf_op", then output the rf mode ("krfx-xx") files.
**          ="enbl+krf", then both re-write ray data into ray data file
**                       and output the "krfx-xx" files.
**          BH191201: Updated ray data is written in to netcdf output file
**                    mnemonic_krfnnn.nc, where nnn= "wave type" number
**                    at the last time step.
**                    
**
***************************************************************
**
**  BEGIN TRANSPORT MODULE INPUT..... "TDTR" ROUTINES
**  namelist trsetup
**
***************************************************************
**
**  A radial diffusion and pinch operator is available as a means
**    for simulating energy and particle transport.  There are two
**    methods for integration of the resulting 3D 
**    Fokker-Planck equation:
**    (1)a splitting scheme which alternates between implicit in time
**    solutions of the velocity space equations and the radial equation,
**    and,
**    (2) an alternating- direction- implicit method which alternates
**    between  solving the velocity space equation implicitly with an
**    explicit radial term, and the radial equation implicitly
**    with an explicit velocity space term.
**    (3) soln_method='it3drv', as descibed above......
**
**  transp="enabled" allows for transport calculations,
**          (lrzdiff.ne."disabled", i.e. lrz=lrzmax, for cql3d)
**        ="disabled" excludes transport
**          default: "disabled"
**
**  nontran:  time-step at which transport is turned on (default=0)
**  nofftran: time-step at which transport is turned off(default=10000)
**
**
**  adimeth="disabled", with soln_method='direct', 
**           then use splitting scheme (default).
**         ="enabled", the use ADI method.
**           [BH_070525: As I presently understand, this ADI option has
**                       not worked for years.   It might be useful
**                       to go over it again.  I use the splitting
**                       algorithm, and expect to use fully-implicit,
**                       iterative 3D solution when it becomes available.
**                       BH100608: soln_method='it3drv' now available.]
**
**  nonadi=Value of time step at which radial transport starts.
**         (Applicable for adimeth="enabled")
**         default: 5
**
**  difus_type(k)="specify" (default), then specify according
**                       to difusr,difus_rshape,difus_vshape below.
**                       Superceded for difus_io=drrin or drrdrin (below).
**            ="neo_smpl", use ion neoclassical diffusion coeffs,
**                         constant in velocity space, as given
**                         in subroutine tdrmshst
**            ="neo_plus", use ion neoclassical diffusion coeffs,
**                         as above "new_smpl" choice,
**                         PLUS drr (below) as in "specify" case.
**            ="neo_trap", use ion neoclassical diffusion coeffs which
**                         are constant in trapped region of velocity
**                         space, but sqrt(R/r) larger than "neo_smpl"
**                         case.
**            ="neo_trpp", use ion neoclassical diffusion coeffs,
**                         as in "neo_trap" case,
**                         PLUS drr (below) as in "specify" case.
**
**  difus_io(k) option enables read/write of radial diffusion and/or pinch
**              velocity terms from/to a NetCDF file, for general species k.
**              Thus, can create
**              the diffusion coeffs in cql3d using the various
**              difus_type(k) options, output them to a .nc file, thus
**              creating a template for an input file, or 
**              test the io system by reading in the previously
**              output coefficients.  Alternatively, create the diffusion
**              coeff file with an external program, read it in.
**              Coefficients are on the native theta,vel,radial grids
**              in the code.
**              When reading in data (difus_io(k)="drrin" or "drrdrin", then
**              the read in coefficients supercede the above difus_type 
**              settings.
**  difus_io(k)="disabled" (default), then no read/write diffusion
**                          coeffs from/to a file.
**             ="drrout", then output the d_rr(i,j,k,l=1,lrz) and related
**                        data to a .nc netCDF file (units: cm**2/sec)
**             ="drrdrout", then also add d_r (i,j,k,l=1.lrz) pinch vel 
**                          term to the .nc file. (units: cm/sec)
**             ="drrin", then read d_rr(,,,) from a file, and apply at each
**                       time step.
**             ="drrdrin", then in addition to d_rr, read d_r(,,,).
**             Presently, there is no time-dependence enabled for the
**             read in or output coeff file.  In cql3d, d_rr is specied
**             independent of time step; however the pinch term d_r is 
**             calculated at each time step.  For "drrdrout" the
**             output pinch term is for the last time step.  For "drrdrin",
**             the same pinch term will be used for all time steps,
**             and pinch is set equal to "disabled".
**
**             For ngen.gt.1 cases, those species with 
**             difus_io(k).ne."disabled" must be the lower (first) 
**             k species, and have the same settin as difus_io(1).
**
**  difus_io_file=character*256 name/path of .nc file to read or write.
**            ="mnemonic", then the &setup0 namelist variable is used to
**              generate the file name mnemonic//"_difus_io.nc"
**                                             [//=concatenate sign]
**            =any name other than "mnemonic", then that name is used
**             for the NetCDF file. Include the .nc .
**             (default: "drrin.nc")    
**  
**  difus_io_drrscale(1:ndifus_io_t,k) and 
**  difus_io_drscale(1:ndifus_io_t,k) specify time-dependent scale
**      factors of drr and dr, respectively, for difus_io(k)="drrin"
**      and "drrdrin" cases (only difus_io_drrscale applied for drrin case).
**      k=1,ngen.
**  difus_io_t(1:ndifus_io_t)=time array  
**  ndifus_io_t=length of the relevant scale factor arrays (indep of k).
**      Defaults: ndifus_io_t=0, and difus_io_drrscale(,)/
**                difus_io_drscale(,)=1.,difus_io_t()=0.
**    The scaling is implemented by specifying the array
**    variables difus_io_drrscale() and difus_io_drscale() at a 
**    sequence of ndifus_io_t corresponding times, specified through
**    the variable difus_io_t() (seconds).
**    Interpolation of the scale factors between time points is linear.   
**    If simulation time becomes greater than 
**    difus_io_t(ndifus_io_t), 
**    then scaling continues at the corresponding last scaling.
**    ndifus_io_t =0 gives no scaling.
**    (ndifus_io_t must be .le. the parameter nbctimea, presently =101).
**
**  difusr is the central radial diffusion coefficient.
**     Units are: cm**2/sec 
**     default: 1.e4
**
**  difus_rshape(1:7) specifies scaling with radius:  
**     With cn denoting difus_rshape(n):n=1,7,
**     (c1 +c2*(rho/radmin)**c3)**c4 *(n_e/n_e0)**c5 
**                *(T_e/T_e0)**c6 *(Zeff/Zeff_0)**c7
**     The Maxwellian values are used for n_e and T_e, unless
**     colmodl=1, in which case self-consistent values are used.
**     defaults: c1=1.0,c2=3.0,c3=3.0,c4=1.0,c5=-1.0,c6=0.0,c7=0.0
**               (Corresponding to previous nominal spatial profile.).
**     NOTE: 050921, appears previous radial dependence specification
**                   was not properly implemented, giving constant
**                   radial dependence of the diffusion in cases
**                   where diffusion was independent of velocity.
**                   In this case, present comparable settings are
**                   c1=1.0, c2=0.0 ,c[3-7]=arbitrary [BH].
**  difus_rshape(8) specifies additional drr radial shape which
**     is constant to normalized radius difus_rshape(8), then
**     decreases to zero as 0.5*(1.+cos(2*pi*rho/(0.1*difus_rshape(8)))).
**     This could be used for added diffusion due to sawteeth.
**     With difus_rshape(1:2)=0., it will be the only added diffusion.
**     default: difus_rshape(8)=0.0    
**  
**  difus_vshape(1:4) specifies scaling with velocity:  
**     With cn denoting difus_vshape(n),
**     |vpar/vth(k,l=1)|**c1 /[1.+l_autocorr/lambda_mfp]**c2
**                 *(vprp/vth(k,l=1))**c3 /gamma**c4.
**     The l=1 "central" value of vth are used for normalization
**     of vpar and vprpr, for electrons or ions.
**     The autocorrelation length is taken to be pi*qsafety*radmaj.
**     The relativistic factor gives a heuristic cutoff of the
**     transport at high energy, e.g.,  due to drift orbit effects. 
**     defaults: c1=0.0,c2=0.0,c3=0.0,c4=0.0,
**               (gives previous constant in velocity space coeff).
**     [For Rechester and Rosenbluth-type diffusion due to small
**      stochastic magnetic field perturbations, c1=1.,c2=1.,c3=0.,c4=0.
**      The c2=1. term gives a collisional reduction of the 
**      collisionless(c1=1. term) diffusion, also in R&R.]
**
**  drr_scale= 1.d0(by default) !YuP[2021-08] added.
**             Scaling factor for Drr coefficient;
**             Used only in case of read_data='nimrod' 
**             (or could be other code, in future),
**             when time-dependent profiles of Drr are formed
**             from imported data on dB/B fluctuations.
**             [For other options described above, this scaling factor
**             is not needed, because the magnitude of Drr is determined 
**             by value of difusr, for example.]
**
**  difin(ryain): Defines the radial diffusion profile using a set of
**     points at r=ryain in a similar way as tein, etc. 
**     If sum(abs(difin))>0, then this option is used instead 
**     of the profile from difus_rshape(1-7).
**     The difin profile is also multiplied by difusr.
**     The additional contribution from difus_rshape(8) is still added.
**     It is used in the file tdtrdfus.f.
**     It is useful since a good pragmatic definition is 
**     D(r) = D0 * chie(r), with chie profile from power balance
**     and D0~0.2 (as tested on TCV for example).
**  
**  pinch: specifies options available to maintain the density
**     profile in the presence of radial transport by a radial
**     pinch velocity.  Various cases correspond to whether
**     the radial diffusion term d_rr and/or the pinch
**     velocity d_r are velocity dependent.
**     A Newton-Raphson iteration for d_r is available (soln_method
**     .ne.'it3drv') increasing the ability to maintain the density
**      profile for splitting method soln or radial transport eqns,
**     and is generally recommended (see "n" appended variables below).
**  pinch is only set up for ngen=1, at the present.  The target
**     density profile is the 1st Maxwellian species of electrons
**     or ions, depending on which species is being transported.
**     It is not difficult to extend it to multiple species.
**     Modifications in trtravct.f are required.
**     The specifics on pinch are:
**     
**  pinch="simple", velocity independent d_rr and d_r.
**        i.e., they depend only on plasma radius.
**        Non-cirular geometry effects are ignored.  
**        The pinch  velocity is computed to preserve the
**        initial density of the electron Maxwellian species
**        (but modified by advectr, as below).  
**        Relaxation or Newton-Raphson is not used.
**       ="case1", the advective coeff and the 
**        radial diffusion coeff are velocity independent.
**        Includes geometry+relaxation, not included in "simple".
**       ="case2", the advective coeff is indep of
**        velocity but the radial diffusion coeff can be
**        velocity dependent (as specified by difus_vshape).
**       ="case3",  we assume that the advective coeff has the same
**        velocity dependence as the radial diffusion coeff, i.e.,
**        d_r=adv(l)*d_rr. (d_rr vel shape specified by difus_vshape.)
**  "n" appended to the above options above (i.e., pinch=
**     "simplen", "case1n", or "case2n"), or, "case3n"
**     then,
**     the density is further maintained  by adjusting the radial
**     velocity at each time step using a Newton-Raphson interation
**     to obtain the velocity which keeps the density at the target,
**     i.e., at the initial density profile.
**     In general, the Newton techniques seems best, and adds
**     little computer time to the runs.
**    Default: pinch="disabled":  Radial advection term is 0.
**      (unless difus_io(1)="drrdrin").
**      (May also work with multiple species, but  
**       hasn't been tried: BH020304).
**       pinch="disabled" is also appropriate setting if difus_io(1)=
**       "drrdrin", but if not, it will be reset to "disabled".
**  
**
**  relaxden=relaxation coefficient for pinch term towards initial
**           electron density profile. 
**          =1., gives relaxation in one time step; smaller values
**              under-relax the density.  (default: 1.)
**           Prior to 091209, values of .le.1/dtreff were required.
**           Adding the 1/dtreff factor into the eqn==> use
**           relaxden as described above, i.e., .lt. but ~1. 
**
**  advectr=1., (default) is multiplier of the radial pinch velocity.
**         =0., gives zero radial pinch term. 
**
**
**  relaxtsp="enabled", then average distribution functions over 
**    previous time steps, before taking next velocity time step 
**    (adimeth="enabled").
**    default: "disabled"
**
**  {meshy, which is in namelist "setup" affects the radial diffusion.
**          "free" ==> theta meshes on flux surfaces independent
**                     of each other (not appropriate for 
**                     transp="enabled")
**           fixed_y=> theta meshes are same on each flux surface,
**                     giving radial diffusion at constant theta.
**           fixed_mu=>theta meshes on each flux surface chosen to
**                     give radial diffusion at constant magnetic moment
**                     mu.  (Most appropriate for transport due to
**                     low freq. perturbations).}
** 
**  PLOTTING OPTIONS FOR RADIAL TRANSP
**
**  pltdrr="enabled" means plot contours of Drr == d_rr(i,j,k,lr)
**      in momentum 2D space (v_par,v_perp) [normalized by vnorm or c]
**      at time steps selected by nplot=*,
**      and such that n.ge.nontran .and. n.le.nofftran
**    ="first" plot contours (black lines), only if n=nontran
**      (or first time the subr.pltmain is called at n.ge.nontran)
**    Also, options to plot with colors :
**    ="color" means plot contours of the Drr in color,
**      at time steps selected by nplot=*,
**      and such that n.ge.nontran .and. n.le.nofftran
**    ="first_cl" plot contours in color, only if n=nontran
**      (or first time the subr.pltmain is called at n.ge.nontran)
**    default: pltdrr="first"            Added: YuP[2021-08]
**         
**
**************************************************************
**
**  BEGIN "FR" MODULE INPUT
**  (THE CODING IS TAKEN FROM ONETWO - GA TECHNOLOGIES TRANSPORT CODE)
**  namelist frsetup
**************************************************************
**
**  This module can work for positive or negative ion based neutral
**  beams.  For negative ion beams, the full energy component is
**  generally very domoninant.   
**  For positive ion based NB sources such as TFTR, DIII-D, JET,
**  [quoting D.M. Thomas et al. (Rev. Sc. Instrum., 2012)]:
**  "As in all such arc discharge sources, the extracted positive ion 
**  current is distributed among the atomic and molecular ions D+, 
**  D_2+, D_3+, with additional small fractions due to such molecular 
**  ions as DO+, D_2 O+, and others. Subsequent collisions in a 
**  neutralizer cell ﬁlled with D_2 gas leads to the production of 
**  full, half, third, and other fractional energy components in the 
**  resulting neutral beam due to the various dissociation and charge
**  transfer processes.  While standard operation typically attempts
**  to maximize the full energy component, signiﬁcant half- and
**  third-energy components exist."
**
**
**  (Might want to look at subroutine freya comments for a little
**   more specific information.)
**
**  frmod           ="enabled" enables the NFREYA NB deposition.
**                      default:frmod="disabled"
**                  MUST also set 
**                   nso=max number of source per species, as above,
**                   nonso(kfrsou(1:nbeams),1:nbeams),
**                   noffso(kfrsou(1:nbeams),1:nbeams) in the above 
**                   sousetup namelist.
**                   (kfrsou(1:nbeams), below,  is set in setup namelist.)
**             N.B.[YuP]: For NBI to be intact, also MUST have nsou>nstop,
**             otherwise it will be reset to zero at n=nsou+1, 2*nsou+1, ...
**             in subr. sourcee.
**
**  nbeams          nbeams is the number of beams, indexed by ib, below.
**                  (max=kb=16). Used to set mb=nbeams within freya.f.
**
**  kfrsou(1:nbeams) [in namelist setup] = general specie index to which 
**                               the ib'th NBI source is applied to. 
**                               Default kfrsou(1:kb)=0, Maximum kb=16
**                               set in param.h. The number of distinct
                                 beams species must be .le. ngen.      
**                               The number of non-zero entries should be
**                               consistent with nbeams.
**                               Also, gives the species of the NB.
**                               Order kfrsou() so that beams of the
**                               same species are contigous 
**                               (beamlines involving the same species
**                               are assumed grouped together in kfrsou())
**                               Prior to 20171016, only one beam species
**                               could be specified, and kfrsou was a scalar.
**  frplt           ="enabled" allows a plot of a selected number of
**                    birth points on a cross-sectional tokamak 
**                    background.
**                  ="plotwrit" plots the nfrplt points and also
**                    outputs points to ascii file, freya_points.txt
**                  ="write_op" gives no plotting, but output points
**                    to ascii file, freya_points.txt
**                  default:frplt="enabled"
**  nfrplt          The approx number of birth points selected for 
**                    the plot (.le. npart, below).
**  smooth          smooth greater than .005 calls a source
**                    smoothing function (sub. frsmooth) which weights
**                    the FREYA generated source current by a Gaussian
**                    function with a standard deviation which is
**                    normally a small fraction of the radius, radmin.
**                    smooth would typically be about .1
**                    If it is less than .005, smoothing is disabled.
**                      default: smooth=.05
**                      [BH110314:
**                       Best probably to increase npart to 1000000.
**                       Had negligible effect on execution time.]

**  multiply        character*8 to enable multiplyn
**  multiplyn       If integer variable multiplyn is gt. 1, 
**                    and multiply.ne."disabled" we assume that each
**                    FREYA particle birth point represents the mean
**                    value of a normally distributed beamlet 
**                   (of multiplyn
**                    particles) with a standard deviation of bmsprd.
**                    multiply="disabled" means disengage this mechanism
**                    Use smooth first, if needed then try multiply.
**                    default: multiply="disabled",multiplyn=0
**  bmsprd          This standard deviation (see above) is normalized to
**                    to the minor radius, radmin. A typical choice 
**                    would be .1  [BH081204: Vestegial, as has no effect.]
**  kfrsou(1:nbeams)  (in namelist setup) is the gen species number of 
**                  the beam ions
**
**  nimp            nimp is the number of impurity species.
**                  (Only 0 or 1 is allowed, for now; see frinitz.f and frstup.f).
**                  YuP[2023] Revised version in subr.frstup to ID the impurity ion:
**                  Scan the Maxwellian_ion list, 
**                  skip those ions that are present among general_ion list
**                  (it happens in colmodl=3 runs), 
**                  take the first ion in the remaining list as the impurity ion.
**                  See subr.frstup for details.
**                  In screen/log file, verify this printout:
**                  'frstup: Impurity species index k=kimp_='
**                  where kimp_ points to the species index in cqlinput list.
**                  If no impurity (nimp=0), kimp_ is 0.
**
**  nprim           nprim is the number of primary species.  
**                  Apparently must be .ge.(number of beam species).
**                     (See frinitz.f, re ibion(). BH170807)
**                  YuP: In cqlinput, position those nprim ions in the beginning of list.
**                  E.g., if D and T beams are used, position D and T species
**                  in the beginning of general_species list, and use nprim=2.
**                  (Note that ngen can be larger than 2 in this example,
**                  if electrons are general species, too).
**
**  anglev(ib)      Vertical angle (degrees) between optical axis
**                    and horizontal plane; a positive value indicates
**                    particles move upward
**  angleh(ib)      Horizontal angle (degrees) between optical axis and
**                    vertical plane passing through pivot point and
**                    toroidal axis; a zero value denotes perpendicular
**                    injection, while a positive value indicates par-
**                    ticles move in the co-current direction.
**                    ==> Tangency Radius = sin(angleh) * rpivot
**  nsourc          Number of sources per beamline.
**                    If 1, source is centered on beamline axis.
**                    If nsourc=2, distinguish between the beamline
**                    axis and the source centerline (optical axis).
**                    The two sources are assumed to be mirror images
**                    through the beamline axis.
**                    In either case, the exit grid plane is 
**                    perpendicular to the beamline axis, and 
**                    contains the source exit grid center(s).
**                    If nsourc=2, the alignment of the sources w.r.t.
**                    the beamline axis is specified through bhofset,
**                    bvofset, and bleni (described further below).
**  bvofset(ib)     Vertical offset from beamline axis to center
**                    of each source (cm; used only for nsourc=2)
**  bhofset(ib)     Horizontal offset from beamline axis to center
**                    of each source (cm; used only for nsourc=2)
**  bleni(ib)       Length along source centerline (source optical axis)
**                   source to intersection point with the beamline axis
**  sfrac1(ib)      Fraction of source current per beamline coming
**                    from upper source (used only for nsourc=2)
**  bcur(ib)        Total current (amps) in ion beam (used only if bptor
**                    is zero)
**  bptor(ib)       Total power (watts) through aperture into torus; when
**                    nonzero, bptor takes precedence over bcur
**  bshape(ib)      Beam shape
**                    'circ' : circular
**                    'rect' : rectangular
**  bheigh(ib)      Height of source (cm)
**  bwidth(ib)      Width of source (cm); diameter for
**                     circular source.
**  bhfoc(ib)       Horizontal focal length of source (cm)
**  bvfoc(ib)       Vertical focal length of source (cm)
**  bhdiv(ib)       Horizontal divergence of source (degrees)
**  bvdiv(ib)       Vertical divergence of source (degrees)
**                  Checking in freya.f, the starting velocities
**                  at the source are distributed with angles
**                  to the normal given by 
**                  f(theta)=1/(theta_div*sqrt(2)) * exp(-(theta/theta_div)**2)
**                  where theta_div is bhdiv or bvdiv, as appropriate.
**                  [BH, 030611].
**  ebkev(ib)       Maximum particle energy in source (keV)
**  fbcur(ie,ib)    Fraction of current at energy ebkeV/ie
**  npart           Number of particles followed into plasma
**                    (suggest 100000)
**  npskip          Ratio of number of particles followed into plasma
**                    to number of source particles (suggest 1)
**  naptr           Total number of aperatures encountered by a particle
**                   as is moves from the source into the plasma chamber
**                    Maximum is specified by parameter nap (=10).
**                    First set of aperatures encountered by the particle
**                    are assumed centered on the source axis, and subseq
**                    aperatures are centered on the beamline axis; the
**                    distinction is made through ashape.
**  ashape(iap,ib)  Aperture shape.
**                   Prefix 's-' indicates source axis centered.
**                   Prefix 'b-' indicates beamline axis centered.
**                     's-circ'          'b-circ'
**                     's-rect'          'b-rect'
**                     's-vert'          'b-vert'
**                     's-horiz'         'b-horiz'
**                                       'b-d3d'
**                   (circ=circular aperature, rect=rectagular,
**                    vert=limits vertical height of source particles,
**                    horiz=limits horizontal height of source particles,
**                    d3d= special DIII-D polygonal aperature)
**                    'circ' : circular
**                    'rect' : rectangular
**  aheigh(iap,ib)  Height of aperture (cm)
**  awidth(iap,ib)  Width of aperture (cm); diameter for circular
**                    aperture
**  alen(iap,ib)    Length from source to aperature for 's-type' 
**                    aperatures and from exit grid plane along beamline 
**                    axis for 'b-type' aperatures.
**  blenp(ib)       Distance along beamline axis from source exit
**                    plane to the fiducial "pivot" point.
**  rpivot(ib)      Radial position of pivot (cm)
**  zpivot(ib)      Axial position of pivot (cm)
**  iexcit = -1     uses Stearn's average cross sections:otw. like =0
**  iexcit = 0      (default mode) Original Freeman-Jones cross section (crsecs)
**  iexcit = 1      (not active) use hexnb instead of crsecs to get cross 
**                  sections but neglect presence of excited states in beam
**                  2=use hexnb with excited beam states allowed.
**  iexcit = 5      use JET Atomic Data and Analysis Structure (ADAS) package
**                  Plasma ions are assumed to be fully stripped. Beam stopping 
**                  cross sections are returned in the array, sgxn, for each 
**                  beam and beam energy component, and flux zone. 
**                  NOTE: ascii cross section files are located 
**                        in the subdirectory ./data/adas/
**                        Can create a link in the working directory to a
**                        given location of the data/adas files.
**                  In iexcit=5 cases, set the following additional variables:
**                  fe_tk=factor (>1) to set upper limit of energy in 
**                     n*sigma array.  max(ebins) = max(ebkev)*fe_tk.
**                     Required for nonzero rotation cases.
**                     Default: fe_tk=0.d0
**                  ne_tk= number of equi-width energy bins used in forming
**                     n*sigma array.  required for nonzero rotation cases. Is
**                     internally reset to zero (used as flag to turn off 
**                     rotational effects on neutral stopping) if angular 
**                     rotation is not present (iangrot = 0).
**                  ds_tk= maximum trajectory increment (cm) used in subroutine
**                     INJECT1 to calculate psi(s), where s is the neutral 
**                     trajectory pathlength from the first closed flux surface
**                     it encounters.  Used for non-zero toroidal rotation,
**                     where mean free path as a function of path length
**                     is required.
**                     Default: ds_tk=0.d0
**  inubpat          switch controlling 2-d beam deposition calculation
**                   0, (default) bypass 2-d deposition
**                   MUST BE = 0 FOR CQL RUNS
**  NOTE:  following namelist variable, eleccomp, must be set in 2nd /setup/
**  eleccomp        If "enabled", then compute the analytic (Cordey-Start) 
**                  electron current degradation factor for NBI type 
**                  calculations.  This assumes that electrons are not
**                  a general species. 
**                  If they are, the electron contribution is computed
**                  directly from the excited electron current.
**                  Default="enabled".
**
**  ranseed=seed for random number generator, random.
**        Random() is used in Freya NB source modules,
**        but may be used elsewhere.   (default: ranseed=7**7)
**
**  beamplse .ne.'disabled', then turn on square wave oscillating 
**    NBI power with frequency and pulse length specified 
**    thru beampon and beampoff.  The frequency is 1./(beampon+beampoff),
**    and beampon gives the on-time interval during the pulse cycle,
**    beampoff gives the off-time interval during the pulse cycle.
**    If beamplse.eq."tauee", these times are measured in units
**    of tauee, else in terms of seconds.  Presently only set up for
**    one beam line (the first).
**  That is,
**  beampon  is the on time during each cycle.
**  beampoff is the off time during each cycle.
**  Note:  nonso and noffso are also used, to time-window the
**         beampon/beampoff, and the cycle will start with
**         beam on at time step nonso.
**         The user might want to coordinate the time-steps and the
**         beampon and beampoff times.  However, if dtr/dtr1 is
**         much less that beampon/beampoff, then one might be able
**	   to neglect exactly when the beam states change.
**  Note:  Applied to all particle sources, 1:nso. Presently,
**         all particles from the NBI are added into one source
**         array.  If needed, the code could readily be modified to
**         have the various beamlines contribute to separate
**         source arrays, dimension the beampxxx accordingly,
**         and apply the beamplse coding to each source [BH120606].
**
**  read_birth_pts="enabled", then read a set of birth point files
**         from NUBEAM (or similar format from other codes),
**         to compare with results from freya calculation [BH130401].
**         Beam geometry, as above will not be used.
**         Each subsequent call to freya will use the next file
**         in a sequence, coordinated in time to the cql3d run.
**         A single beam is assumed, with up to three beam energy 
**         components.
**         Make sure fr_gyro="disabled", unless want the gyro-shift.
**         default: "disabled"
**  nbirth_pts_files=number of birth point files, default=1
**  birth_pts_files=input nbirth_pts_file filenames into
**                 the namelist, in the order which they
**                 will be used in cql3d. character*256.
**                 Max list length=24.  default='notset'
**  nbirth_pts= should be same as in all the NUBEAM files. default=250000
**  The shine through, and total power is read from the files.
**  nbirth_points is checked.  Power and shine through numbers 
**  may vary from file to file.
**
**  fr_gyro="enabled" gives gyro-radius offset of calculated NBI
**          birth points, to the guiding center position birth point.
**          Should be disabled, in case this offset already included for
**          read_birth_pts="enabled", for example, from NUBEAM.
**          default = "disabled".  
** 
**
**      
**
** ---------------------------------------------------------------------
**   HOW TO UPDATE NBI/FREYA DEPOSITION AT EVERY TIME STEP
** ---------------------------------------------------------------------
** If plasma density and/or temperature change significantly,
** the NBI deposition should be updated.
** As of [2022-12] the procedure is developed as follows.
** Use the logic related to plasma rotation, but without toroidal speed,
** i.e., set (in &setup section of cqlinput):
**   nonvphi=2 (example)  and  noffvphi=10000 (example)
** The FREYA module will be called for time steps n according to 
**   if( (n+1).ge.nonvphi .and. (n+1).lt.noffvphi )then...
** See tdchief.f, line ~165. Note that in this logic condition
** the time step counter n starts with n=0.
** Leave default values of vphic()=0 and iprovphi="disabled",
** which means - no toroidal rotation of plasma in this case
** (if rotation is really needed, then - set those values).


** Description of nprim,nimp as in ONETWO TRANSPORT code, from which
** freya subroutines were originally extracted
** From onetwo, Version 3.0: cray102.f, Obtained from GA:onetwo_cvs_070823_genr
c ----------------------------------------------------------------------       @
c --- ION PARAMETERS                                                           @
c ----------------------------------------------------------------------       @
c% nprim       Number of      primary ion species (1 to 3)                     @
c% nimp        Number of     impurity ion species (0 to 2)                     @
c% namep(i)    Name   of ith  primary ion species                              @
c              'h' , protons                                                   @
c              'd' , deuterons                                                 @
c              't' , tritons                                                   @
c              'dt', mixture of d and t                                        @
c              'he', thermal alphas                                            @
c% namei(i)    Name   of ith impurity ion species                              @
c              'he', helium                                                    @
c              'c' , carbon                                                    @
c              'o' , oxygen                                                    @
c              'si', silicon                                                   @
c              'ar', argon                                                     @
c              'ti', titanium                                                  @
c              'cr', chromium                                                  @
c              'fe', iron                                                      @
c              'ni', nickel                                                    @
c              'kr', krypton                                                   @
c              'mo', molybdenum                                                @
c              'w' , tungsten                                                  @
c% fd          Number fraction of deuterons in d-t mixture                     @
c ----------------------------------------------------------------------       @
  


***************************************************************
**
**  BEGIN PARALLEL TRANSPORT MODULE INPUT..... "WP" ROUTINES
**   (CQLP CODE (CQL_PARALLEL))
**  
**  A new option, cqlpmod="enabled", has been introduced which in
**  effect makes up for a new code CQLP. This code solves the FP equations
**  keeping the variation along the magnetic field and thus without bounce
**  averaging. The dimensions are thus (u, theta, s), where s is distance
**  along the magnetic field line. In practice, one had only to change
**  the spatial dimension to be s instead of rho in the CQL3D case, and
**  to remove the bounce-averaging. This is why most of the arrays used
**  for CQL3D are used for CQLP, with the index l running over the s-mesh
**  instead of the radial-mesh. The following additional input parameters have
**  been introduced:
**  
**  In first namelist &setup0:
**  
**    cqlpmod:  "enabled" run CQLP (default: "disabled")
**    lsmax:    number of mesh points along s
**    ls:       number of s-mesh points for which the FP eq. is solved
**    lsdiff:   "enabled" allows for ls.ne.lsmax
**    lsindx(0:ls): list of the ls indices at which the FP eq. is solved
**    nlwritf:  ="enabled"  saves the distribution function at the end of 
**               the run (in file distrfunc).   default="disabled"
**    nlrestrt: ="enabled" reads the initial distribution function in 
**               file distrfunc.                default="disabled"
**    BH070414:  Variables nlwritf and nlrestrt have been changed from
**               logical to character*8.  Might need to adjust old NL.
**    
**  In namelist &setup:
**  
**    denpar:   density along s, defined in the same way as reden
**
**    temppar:  temperature along s, defined in the same way as temp
**
**    symtrap:  ="enabled" (default) assumes trapped region symmetric and solves
**              only for itl.le.i.le.iyh, and not for iyh+1.le.i.le.itu.
**              ="disabled" solves for all i in [1,iy]
**              For CQL3D, symtrap can be either but normally should be "enabled"
**              For CQLP, symtrap has to be "disabled" (set in ainsetva)
**    
**    tfac:     same as before, except add the option for tfac<0, with meshy="fixed_mu"
**    epsthet:  which constructs the theta mesh such that at each s position, we
**              have: y(iyh_(l),l)=pi/2-epsthet. Then distributes the rest of the
**              points. Default is epsthet=0.1
**    
**    ngauss,nlagran:
**              Determines how the integral over theta for computing the Legendre
**              coefficients is done:
**                 ngauss=0: old method: assumes f cst over [th(i),th(i+1)] and
**                                       integrates analytically the rest
**                 ngauss.ne.0: integrates using a ngauss-point formula on each
**                              sub-interval and interpolates the value of f at
**                              the Gauss-point using a nlagran-point Lagrange
**                              interpolation.
**    oldiag:   ="enabled" (default) uses old way of computing moments of f in
**              diaggnde subroutine (for reden, curr,..)
**              ="disabled" decompose f with Legendre polynomial as is done for
**              computing the collision operator.
**  
**    lbdry0:   ="enabled" changes matrix construction such that f is unique
**              at u=0. (default: enabled)
**
**    nummods:  Determines numerical scheme used to solve the spatial time-step
**              Following are the comments in subroutine wptramu,wptrafx:
**
**     0. The parameter nummods determines the numerical model used to solve
**        the equation along s. nummods is divided into classes of multiple
**        of 10: numclas=nummods/10 =>   nummods=[0,9]; [10,19]; [20,29] ...
**        numclas = 0: normal case
**                  1: assumes up/down symmetric case and periodic condition
**                     but compute equil. parameters on top half cross-section
**                     on ls/2+1 points and then copy the rest in sub. wploweq
**                     Otherwise, same as numclas=0 with sbdry="periodic"
**                  [YuP: it seems numclas>1 is not applicable in the code]
**        In each class, nummods is divided into two groups:
**           mod(nummods,10) < 5: 2-D FP velocity equ. solved on nodal s points
**                           > 4: velocity eq. solved at s+1/2
**        Then, in these two groups one has the following options, 
**        with numindx=mod(mod(nummods,10),5)=0,1,2,3 or 4, we have:
**          If numindx = 0,1: use backward diff. scheme + bound. cond at l=1
**                         2: "  back/forward depending on sign(cos(theta))
**                         3: " forward diff. + bound cond at l=ls
**                         4: " centered diff. + bound. cond. at l=1
**        Note: if sbdry="periodic" there might not be any extra boundary
**              conditions imposed (see wpbdry)
**        Default is nummods=1. Recommended 13 or 14, see comments by YuP below.
**        [YuP: From tests, nummods=10 gives exactly same results as nummods=11]
**        [YuP: For advection-dominated cases (low-collis; starnue<0.1) nummods=11,12,13
**         are the best options, but for higher collisionality (starnue>0.1),
**         those options give non-conservation of density and energy,
**         and visible dip or peak defect in distr. func. at v~0; 
**         so, for starnue>0.1, nummods=14 is best.]
**        [YuP: From theory of ADI methods,
**         the upwind differencing, in the direction of advection term,
**         helps to reduce odd-even (in s) oscillations, common for ADI when
**         advection term dominates over diffusion (collision) term.
**         For centered differencing (i.e. nummods=14),
**         in case of advection-dominated FPE, the refined s-grid 
**         helps to reduce odd-even oscillations; i.e., increase ls.] 
**              
**    lmidvel:  =0 (default) assumes velocity source term computed at s-nodes
**              =1 assumes velsou(l)=source term at l+1/2
**
**    laddbnd:  =1 (default) extra boundary conditions are imposed for the
**                 s-time-step, when sbdry="periodic"
**              =0 No extra boundary conditions (if sbdry="periodic")
**
**    sbdry:    ="bounded" (default) assumes bounded s interval
**              ="zeroslop" same as bounded, but impose zero df/ds at edges
**              ="periodic" assumes s is periodic (ex: on flux surface)
**
**    eseswtch: ="enabled", calculate nonlocal electrostatic electric
**               field to obtain constant electron flux with
**               sbdry="periodic", 
**               using Kupfer et al (PoP, 1995) method.
**               efswtch re-set "disabled", for now (until further
**               work/investigation of this functionality).
**              ="disabled" (default).
**
**    elpar0: coefficient of E_parallel(electrostatic)
**            default=0.
**
**    nkconro(): Gives indices of species which contribute to charge
**               density in Poisson equation
**
**    In namelist trsetup:
**
**    adimeth, nonadi: same as before, ADI method needed
**
**    nontran:  time-step at which transport is turned on (default=0)
**    nofftran: time-step at which transport is turned off(default=10000)
**
**    relaxtsp: ="enabled" define f_n+1 = 0.5*(f_n+1/2 + f_n+1) for
**                n in [nonavgf,nofavgf] 
**              ="disabled" keeps f_n+1 as determined after s-step
**                [YuP: "enabled" shows good stability]
**
**    nonelpr,noffelpr: time-steps in between which the parallel electric 
**                      field due to Poisson eq. is calculated.
**                      (This method not appreciably checked out, OS990315).
**
**    scheck:  Checking wp. See wpcheck.f
**             default="enabled"
**
**
**    To run CQLP, from a cqlinput file which was running CQL with a
**    DC electric field, here are the changes which should be made
**    (solving on flux surface => periodic)
**      
**    in first namelist setup:
**      cqlpmod="enabled",
**      lrz=1,lrzmax=(previous value of lrz),
**      lrzdiff="enabled", lrindx(1)=(desired radial position),
**      ls=20,lsmax=20
**
**    in second namelist setup:
**      symtrap="disabled",
**      nummods=14,
**      sbdry="periodic",
**      meshy="fixed_mu",
**      tfac=-0.5,
**      iy=80,jx=80,
**      denpar(1,0)=   ,
**      denpar(1,1)=   ,
**      denpar(2,0)=   ,
**      denpar(2,1)=   ,
**      denpar(3,0)=   ,
**      denpar(3,1)=   ,
**      temppar(1,0)=   ,
**      temppar(1,1)=   ,
**      temppar(2,0)=   ,
**      temppar(2,1)=   ,
**      temppar(3,0)=   ,
**      temppar(3,1)=   ,
**
**    in namelist trsetup:
**
**      transp="enabled",
**      adimeth="enabled",nonadi=0,
**      nontran=0,
**      relaxtsp="enabled",
**      nonavgf=5, nofavgf=10,
**      
**              
***************************************************************
**
**
*******************************************************************
**
**    BEGIN NAMELIST INPUT
**
********************************************************************
 $setup0
 ibox="box c08",
 iuser="mccoy",
 ioutput=6,
 lrz=1
 mnemonic="92_01_08"$
**
 $setup
 acoefne=-1.80,-7.83,+051.57,-353.68
 acoefte=8.01,-13.60,+08.69,-114.59
 bnumb(1)=-1.,
 bnumb(2)=5.,
 bnumb(3)=-1.,
 bootst="disabled",
 bth=2.80e+3,
 btor=2.8e+4,
 chang="noneg",
 colmodl=3,
 contrmin=1.e-12,
 dtr=1.
 eegy(1,1,1)=0.,
 eegy(1,2,1)=2.,
 eegy(2,1,1)=0.,
 eegy(2,2,1)=6.,
 elecfld(0)=20.e-15,
 elecfld(1)=20.e-15,
 enein=1.53e13,1.50e13,1.43e13,1.20e13,.5e13,.01e13,
 enloss(1)=200.,
 enmax=400.,
 enmin=40.,
 enorm=1200.,
 eoved=.00,
 ephicc=1.,
 fds=.2,
 fmass(1)=9.1095e-28,
 fmass(2)=3.3433e-24,
 fmass(3)=9.1095e-28,
 gamaset=0.,
 gsla=270.,gslb=35.,
 iactst="enabled",
 idskf="dskout",
 implct="enabled",
 ineg="disabled",
 iprone="parabola",
 iprote="parabola",
 iproti="parabola",
 irzplt(1)=1,
 irzplt(2)=2,
 irzplt(3)=3,
 irzplt(4)=4,
 irzplt(5)=0,
 irzplt(6)=0,
 iy=100,
 izeff="backgrnd",
 jx=100,
 kfrsou=0,
 kpress(2)="enabled",
 kspeci(1,1)="e",kspeci(2,1)="general",
 kspeci(1,2)="d",kspeci(2,2)="maxwell",
 kspeci(1,3)="e",kspeci(2,3)="maxwell",
 lbdry(1)="conserv",
 locquas="disabled",
 lossmode(1)="disabled",
 lz=30,
 machine="toroidal",
 manymat="disabled"
 meshy="free",
 mpwr=.5,5.,5.,5.
 mx=3,
 nchec=1,
 ncoef=1,
 ncont=3,
 nen=30,
 ngen=1,
 njene=6,
 njte=6,
 njti=6,
 nmax=2,
 noffel=10000,
 nonel=10000,
 nplot=4,
 nplt3d=15,
 npwr=2.,2.,2.,2.,
 nrskip=0,
 nrstrt=1,
 nstop=4,
 numby=30,
 nv=9,
 partner="disabled",
 plt3d="enabled",
 pltd="enabled",
 pltdn="disabled",
 pltend="enabled",
 pltfvs="disabled",
 pltinput="disabled",
 pltmag=1.,
 pltpowe="last",
 pltprpp="enabled",
 pltrst="enabled",
 pltstrm="disabled",
 pltvecal="disabled",
 pltvecc="disabled",
 pltvece="disabled",
 pltvecrf="disabled",
 pltvflu="disabled",
 pltvs="rho",
 profpsi="disabled",
 psimodel="axitorus",
 qsineut="disabled",
 radmaj=170.,
 radmin=140.,
 rd=45.,
 reden(1,0)=1.5e+13
 reden(1,1)=1.5e+13
 reden(2,0)=.5e+13
 reden(2,1)=.5e+13
 reden(3,0)=1.5e+13
 reden(3,1)=1.5e+13
 relativ="enabled",
 rfacz=1.,
 rmirror=7.5,
 rovera=.05,
 roveram=.00,
 rya(1)=.050,.100,.150,.175,.200,.225,.250,
 rya(8)=.300,.35, .4,.5,.6,.7,.8,.9,
 ryain=0.,.2,.4,.6,.8,1.0,
 rzset="disabled",
 softxry="disabled",
 syncrad="disabled",
 tauloss(1,1)=.3,
 tauloss(2,1)=0.,
 tauloss(3,1)=0.,
 tbnd=.002,
 tein=1.5,1.1,.7,.4,.2,.1,
 temp(1,0)=6.,
 temp(1,1)=6.,
 temp(2,0)=1.,
 temp(2,1)=1.,
 temp(3,0)=6.
 temp(3,1)=6.
 tfac=1.25,
 tfacz=1.,
 thet1=0.,0.,0.,0.,0.,0.,0.,0.,0.,
 thet2=-1.1932,-1.0814,-.9894,-.9087,0.,.9087,.9894,1.0814,1.1932,
 tiin=1.5,1.1,.7,.4,.2,.1,
 torloss(1)="disabled",
 veclnth=1.5,
 xfac=.6,
 xlwr=.085,
 xmdl=.25,
 xpctlwr=.1,
 xpctmdl=.4,
 ylower=1.22,
 yreset="disabled",
 yupper=1.275,
 zmax=408.$
**
 $trsetup
 difusr=2.e+17,
 pinch="simple",
 relaxden=1.,
 relaxtsp="enabled",
 transp="disabled",
 advectr=1.$
**
 $sousetup
 asor(1,1,1)=0.0e+13,asor(1,2,1)=0.0e+13,
 noffso(1,1)=21,noffso(1,2)=21,
 nonso(1,1)=0,nonso(1,2)=0,
 nso=0,
 nsou=1000,
 pltso="enabled",
 scm2(1,1)=.001,scm2(1,2)=10000.,
 sellm1(1,1)=1.,sellm1(1,2)=1.,
 sellm2(1,1)=1.,sellm2(1,2)=1.,
 sem1(1,1)=1600.,sem1(1,2)=0.,
 sem2(1,1)=.5,sem2(1,2)=25.,
 seppm1(1,1)=1.,seppm1(1,2)=1.,
 seppm2(1,1)=1.,seppm2(1,2)=1.,
 soucoord="disabled",
 sthm1(1,1)=5.,sthm1(1,2)=0.,
 szm1(1,1)=0.,szm1(1,2)=0.,
 szm2(1,1)=1.e+5,szm2(1,2)=1.e+5,
 knockon="disabled",
 nonko=10000,
 noffko=10000
 $
**
 $eqsetup
 atol=1.e-8,
 ellptcty=0.,
 eqmod="enabled",
 eqpower=2,
 eqsource="eqdsk",
 fpsimodl="constant",
 methflag=10,
 nconteq="psigrid",
 rbox=92.,
 rboxdst=120.,
 rmag=166.,
 rtol=1.e-8,
 zbox=92.$
**
 $rfsetup
 call_lh="disabled",
 call_ech="disabled",
 call_fw="disabled",
 dlndau=2.
 iurfcoll="disabled",
 lh="enabled",
 ech="disabled",
 fw="disabled",
 iurfl="disabled",
 nbssltbl=8000,
 nrfitr1=100,
 nrfitr2=0,
 nrfitr3=0,
 nrfpwr=0,
 nrfstep1(1)=400,
 nrfstep1(2)=400,
 nrfstep1(3)=400,
 nrfstep2=000,
 noffrf(1)=100000,
 nonrf(1)=0,
 nrf=1,
 scaleurf="enabled",
 urfdmp="secondd",
 urfmod="disabled",
 vlfmod="disabled",
 vlhmod="enabled",
 vparmax=.45,
 vparmin=.2,
 vprprop="disabled"$
**
 $frsetup
 aheigh(1,1)=250.,
 aheigh(1,2)=250.,
 aheigh(1,3)=250.,
 alen(1,1)=3650.,
 alen(1,2)=3650.,
 alen(1,3)=3650.,
 angleh(1)=39.3,
 angleh(2)=39.3,
 angleh(3)=39.3,
 anglev(1)=0.0,
 anglev(2)=0.0,
 anglev(3)=0.0,
 ashape(1,1)='s-rect',
 ashape(1,2)='s-rect',
 ashape(1,3)='s-rect',
 awidth(1,1)=60.,
 awidth(1,2)=60.,
 awidth(1,3)=60.,
 bcur(1)=110.,
 bcur(2)=110.,
 bcur(3)=110.,
 bhdiv(1)=.4,
 bhdiv(2)=.4,
 bhdiv(3)=.4,
 bheigh(1)=100.,
 bheigh(2)=100.,
 bheigh(3)=100.,
 bhfoc(1)=3650.,
 bhfoc(2)=3650.,
 bhfoc(3)=3650.,
 bhofset(1)=0.0,
 bhofset(2)=0.0,
 bhofset(3)=0.0,
 bleni(1)=3650.,
 bleni(2)=3650.,
 bleni(3)=3650.,
 blenp(1)=3650.,
 blenp(2)=3650.,
 blenp(3)=3650.,
 bmsprd=.1,
 bptor(1)=36.25e6,
 bptor(2)=36.25e6,
 bptor(3)=36.25e6,
 bshape(1)='rect',
 bshape(2)='rect',
 bshape(3)='rect',
 bvdiv(1)=.4,
 bvdiv(2)=.4,
 bvdiv(3)=.4,
 bvfoc(1)=3650.,
 bvfoc(2)=3650.,
 bvfoc(3)=3650.,
 bvofset(1)=0.0,
 bvofset(2)=0.0,
 bvofset(3)=0.0,
 bwidth(1)=20.,
 bwidth(2)=20.,
 bwidth(3)=20.,
 ebkev(1)=1300.,
 ebkev(2)=1300.,
 ebkev(3)=1300.,
 fbcur(1,1)=1.,
 fbcur(1,2)=1.,
 fbcur(1,3)=1.,
 fbcur(2,1)=00.,
 fbcur(2,2)=00.,
 fbcur(2,3)=00.,
 fbcur(3,1)=00.,
 fbcur(3,2)=00.,
 fbcur(3,3)=00.,
 frmod="disabled",
 frplt="enabled",
 ibcur=1
 iborb=0,
 iexcit=-1,
 inubpat=0,
 multiply="disabled",
 naptr=1,
 nbeams=3,
 nfrplt=400,
 nimp=1,
 npart=10000,
 nprim=2,
 npskip=1,
 nsourc=1,
 rpivot(1)=900.,
 rpivot(2)=900.,
 rpivot(3)=900.,
 sfrac1(1)=1.0,
 sfrac1(2)=1.0,
 sfrac1(3)=1.0,
 smooth=.125,
 zpivot(1)=0.0,
 zpivot(2)=100.,
 zpivot(3)=+200.$
end
end
end
LEAVE
THESE
HERE!


c***********************************************************************
c
c   Copyright R.W. Harvey and Yu. V Petrov
c   CompX company, Del Mar, California, USA
c   1995-2011
c   Below, "this program" refers to the CQL3D, alternatively designated
c   as CQLP when used in cqlpmod='enabled' mode (and alternative 
c   capitalizations), and associated source files and manuals.
c
c   The primary reference for the code is:
c   ``The CQL3D Fokker-Planck Code'',  R.W. Harvey and M.G. McCoy, 
c   Proc. of IAEA Technical Committee Meeting on Advances in Simulation 
c   and Modeling of Thermonuclear Plasmas, Montreal, 1992, p. 489-526, 
c   IAEA, Vienna (1993), available through NTIS/DOC (National Technical 
c   Information Service, U.S. Dept. of Commerce), Order No. DE93002962.
c   See also, http://www.compxco.com/cql3d.html, CQL3D Manual.
c
c                GNU Public License Distribution Only
c   This program is free software; you can redistribute it and/or modify
c   it under the terms of the GNU General Public License as published by
c   the Free Software Foundation; either version 3 of the License, or
c   any later version.
c
c   This program is distributed in the hope that it will be useful,
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c   GNU General Public License for more details.
c
c   You should have received a copy of the GNU General Public License
c   along with this program; if not, see <http://www.gnu.org/licenses/>.
c
c         --R.W. Harvey, CompX, Del Mar, CA, USA
c
c   E-mail:  rwharvey@compxco.com
c   Address: CompX company
c            P.O. Box 2672
c            Del Mar, CA 92014-5672
c
c   It will be appreciated by CompX if useful additions to CQL3D can be
c   transmitted to the authors, for inclusion in the source distribution.
c
c***********************************************************************

 