********************************************************************
**
** TO RUN CQL3D SET NAMELIST AS DIRECTED BELOW
** AND
** NAME THIS FILE cqlinput
** 2011-2-18
**
** 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.
** 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.]
**
** 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 determines graphical o/p (PS only at moment, so no effect.)
**
** 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="enabled" inhibits all plots.
** default: noplots="disabled"
** [Presently inhibits old unimplemented graflib plots,
** necessary for code to run.]
**
** 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
** and external script containing
** uname -a > uname_output
** pwd > pwd_output
** ="disabled", then no system calls and not reference
** to uname_output and pwd_output files.
** This gets around lack of some special functions,
* e.g., second(), with some compilers.
**
** 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.
** Text distrfunc will contain only namelist and profiles
** at the end of the run.
** nlrestrt="enabled" reads the initial distribution function in
** text file 'distrfunc'. default="disabled"
** ="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 for reading
** namelist from &setup and others (except %setup0),
** ensuring compatibility of grids, etc. (Can change
** mnemonic in &setup0.)
** ='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) 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.
**
************************************************************************
** 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" cases will require dtr values less than or
** of order the collision time of the nonthermal particles,
** to get accurate results.
** 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 for plotting of data for
** individual flux surfaces.
** default: nplot(1:nplota)=0
**
** 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="enabled", 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).
**
** lz is the number of poloidal (or field line mesh points). Bounce-
** averages are computed over this mesh. lz must be .le. lza
** (a parameter).
** default: lz=lza (lza is a parameter)
**
** 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).
**
** 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)
** 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.
**
** 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=.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,
** ="sqvolume" 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.
** default: ineg="enabled"
** ="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.
**
** implct= "enabled" forces implicit differencing and
** Gaussian elimination. implct = "disabled" forces operator
** splitting (used for time dependent calculations).
** If "enabled", use LAPACK linear algebra subroutines.
** 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="enabled". 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.
** 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 after the 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.
** lbdry(k)="conserv" is the ONLY option for implct="disabled"
** default: lbdry()="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 is applied if lbdry(i)="conserv".
** 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 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 Kev 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"
** 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 three possible
** values (below) governing the electron temperature, ion temperature
** and densities of the initial profiles for the species.
** We have the following possibilities:
** iprote="parabola","spline" and "asdex" (similarly for iproti and iprone)
** Each option is described below.
** default: "parabola" for everything
** 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(k)=2. for all k; k,....,ntotal, npwr(0)=2.
** default: mpwr(0)=1.; otherwise mpwr(k)=3.
**
** 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.
** 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 x coordinate)
** 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 (= 101 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 at the same ryain(1:njene).)
**
** 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", or "spline":
**
** = "disabled", then this option has no effect.
** = "parabola" or "spline", 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.
**
** zeffscal=scale factor for zeffin spline input (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...)
** I recommend 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 by specifying
** central and edge values for parabolic-type profiles
** as a function of time. The exponents of the profiles
** are fixed in time, and are those used elsewhere in the
** code.
**
** Two methods are implemented for specification of the time-
** dependent 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.
** (default=0)
** bctime(1:nbctime)=times at which data given (seconds).
** (default= 0, 1., 2., ..., float(nbctimea-1))
**
** 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.)
**
** 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).
** 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 npwrzeff,mpwrzeff,
** specify zeff profiles as described above.
** NOTE: if using this option (i.e., zeffc(1).ne.0.), need to have
** iprozeff="parabolic" (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 npwrelec(real),mpwrelec(real)
** specify electric field profiles, unless current profile
** control is enabled (volts/cm).
** 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.
**
**************************************************************
**
** 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 collisional coefficients (general
** species and background species included)
** 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 would keeps accurate track of momentum
** transfer but would allow the distribution to run away in a tau_e
** time or so, this option would keep the bulk at the background
** electron temperature, but still calculate 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
** 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
** species in the calculation of the collision integral coefficients.
** kfield(k)="disabled" means it is excluded.
** default: kfield(k)="enabled"
**
** gamaset .ne. 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 ions.
** **This needs attention.
**
** gamafac = "enabled", then for electrons as (only) general species,
** multiply the Coulomb logrithms (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 become significant at electron energies
** of order of greater than me*c**2.
**
** 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 F coeff.
** (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 (input) are loss orbits and
** are lost from the system with characteristic time equal to
** the bounce time.
** "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 are lost.
** "mirrsnk"==> means both mirrorcc and snk0 options are used.
** default: lossmode(k)="disabled"
** "simplban"==> then trapped general species ions (or electrons) with
** (banana width (poloidal gyroradius) plus gyroradius) greater
** distance to the outer flux surface are lost in a bounce time.
** See subroutine losscone for a few more details.
** "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*512 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 and eperc: 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 an loss term -f/taulos,
** treated implicitly.
** 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 "flutterx"
** model.
** 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= (dist. 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.
** 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
** parabolic profiles (with iproelec="parabola"),
** 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.)
**
** 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 iterations for electric
** field used AT EACH TIME STEP, to achieve the target
** current within fractional accuracy currerr. For
** each iteration, the electric field is reset according 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)/(amax(curr,currxj),
** where currxj is the parallel flux surface area average
** target current density).
** (See efield.f file).
** (default=0.5)
**
** 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)
** for machine="toroidal". Electric field values
** are evaluated at the major radius of the
** magnetic axis, and field varies on a flux
** surface like 1/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 elecfdl(0),elecfld(1),
** npwrelec,mpwrelec. default: iproelec="parabola".
**
** elecin(1:njena)=spline input at the ryain, of electric field (volts/cm),
** evaluated at major radius of the magnetic axis.
**
** 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.
**
** 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. Positive eoved input is allowed only for lrz=1.
** 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
**
***************************************************************
**
** 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)
** (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)
** pltd="enabled" plot f(n+1)
** 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
**
** 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)).
** jpxy=number of parallel velocities for which vectors are plotted.
** (default: min(101,jx+1)).
**
**
** 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 (_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"
**
** plturfb.ne."disabled" give plotting of QL B coefficient when
** urfmod="enabled"
** default: plturfb="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)=0
** 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 used longer filename.
** Note: file name idskfrf should be quoted: "...".
**
**
** netcdfnm.ne."disabled", then create a netCDF output file
** named according to the character*40 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.eq."enabled", 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.eq."enabled")
** netcdf files, omitting standard f(i,j,k,l) and
** urfb QL coeff at last time step.
** ="longer_f", i.e., not short, gives f at every
** time step (and urfb at last time step).
** ="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 currf(u,r)
** at each step (else only a last time step).
** ='long_urf", urfa,urfb,urfc,urfd,urfe,urff into
** _rf.nc file at last time step, 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.
**
***************************************************************
**
** 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 are output to the
** mnemonic.nc netCDF file. Multiplying by v=v0 gives
** the shift in the normalized minor radius coordinate.
** 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".
**
**
** 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)
**
***************************************************************
**
** 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(thet2)
**
** 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).
**
** 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.
** 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).
** default: npa_diag="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_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(thet2_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:njenn,:) 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.
** 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) = density 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.
**
** 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).
**
** isigsgv1 = 0, Compute reactions of test (general) distribution with
** self and others.
** 1, Also do calc. for Max. with same energy as test distn.
** (default=0).
**
** isigsgv2 = 0, Background Maxwellian distribution not included
** in calculation of reactions,
** = 1, Include background Maxwellian distn in calc of reactions.
** (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.
**
**
**
**
***************************************************************
**
** 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)
** Only use lrzdiff.eq."disabled" with this option.
** 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)
** in banana limit.
**
** 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="eneabled", 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)=1.
** 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 when the source is reset.
** 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']
**
** 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
**
** soucoord="cart" means cartesian coordinates are used
** for input; soucoord="polar" is the other option (see below)
** default: soucoord="cart"
**
** The next eight quantities are in units of Kev.
** The next four input variables 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 variables 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)
**
** pltso="enabled" means plot contours of the source.
** ="first" plot coutours only if n=0.
** 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 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(cutoff of source, per bh), 1(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 time 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).
**
** 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".
** 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
** 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).
**
** 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).
** defaults: "psigrid" (highly recommended operating mode,
** although not sure why, BH990909),
** nconteqn=0
**
** lfield should be .le. lfielda and 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=lfielda and lfielda should be around 250.
**
** 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="enabled")
** 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 single flux surface, multi-
** harmonic, cyclotron model for rf quaslinear diffusion,
** Diffusion coefficient is according to Stix book (1993).
** Region of wave interaction on the 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.
** 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_te**2/tau_coll. (See code).
**
** 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 relativistice 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.
**
** "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()!
**
**
** rdcb diffusion coeffs are plotted for all flux surfaces.
** rdcb written into a netcdf file mnemonic//'_rf.nc'.
**
** default: rdcmod="disabled"
**
** pwrscale(1) [from urf section of code, presently used to
** scale input value of diffusion coeffs.]
**
** 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.)
**
**
**
************************************************************************
**
** BEGIN URF MODULE INPUT
** namelist rfsetup
**
************************************************************************
**
** urfmod="enabled" 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 from up to nmodsa (parameter, typically = 6)
** wave types 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:
** mrf= number of rf "types" (and is determined from rftypes())
** mrfn= number of rf "modes" (sum over mrf of the nharms())
** Thus, it is required that mrfn.le.nmodsa, nmodsa being a
** parameter set in param.h.
**
** The wave types 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.)
**
** Alternatively, wave types can be input as rftype(1:mrf).
** [This is more general, preferred input method of wave types.]
** Wave ray tracing data for each wave type is specified through
** rfread and rffile(1:mrf), as given below.
** 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".
**
** The code will not handle mixed use of ech/fw/lh and rftype().
** 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 input ray data from files
** rayech, if ech="enabled"
** rayfw, if fw="enabled"
** raylh, if lh="enabled"
** ="netcdf", then there are more options for naming the input files,
** as given by rffile(1:mrf)
** (default="text" for backward compatibility,
** 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 &setup 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:nmodsa) if .gt.1 wave types are used.
** [User has to see that data files correspond
** to the order of the enabled subset of the
** 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.]
**
** 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.
**
** 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.
**
** Time-dependent scaling of power is achieved with an addtional
** 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.
** 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).
** 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, nrfstep2=50
** nrfpowr=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 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"
**
** urfrstrt="enabled" => "do not update delpwr" option
** 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")
**
***************************************************************
**
** 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,
** ="disabled" excludes transport
** default: "disabled"
**
** adimeth="disabled", then use splitting scheme (default)
** ="enabled", the use ADI method.
** [BH_070525: As I presently understand, this 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.
**
** 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.
** ="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.
**
** 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.]
**
** 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_rr 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.
** 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.
** (May also work with multiple species, but
** hasn't been tried: BH020304)
**
**
** 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).}
**
**************************************************************
**
** BEGIN "FR" MODULE INPUT
** (THE CODING IS TAKEN FROM ONETWO - GA TECHNOLOGIES TRANSPORT CODE)
** namelist frsetup
**************************************************************
**
** (Might want to look at subroutine freya comments for a little
** more specific information.)
**
** frmod ="enabled" enables the NFREYA NB deposition.
** default:frmod=frmoda (parameter)
** MUST also set nso=1, nonso(kfrsou,1),
** noffso(kfrsou,1) in the above sousetup namelist.
** (kfrsou, below, is set in setup namelist.)
**
** kfrsou [in namelist setup] = general specie which NBI source
** is applied to. Default=0, so must
** set it in the namelist.
**
** frplt ="enabled" allows a plot of a selected number of
** birth points on a cross-sectional tokamak
** background.
** default:frplt="enabled"
** nfrplt The approx number of birth points selected for
** the plot.
** 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, no smoothing is allowed.
** default: smooth=.1
** [BH081204: Vestegial, as has no effect.
** 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 (in namelist setup) is the species number of
** the beam ions
**
** nimp nimp is the number of impurity species.
** nprim nprim is the number of primary species.
** nbeams nbeams is the number of beams, indexed by ib, below.
** 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 (a) in ion beam (used only if bptor
** is zero)
** bptor(ib) Total power (w) 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 10000)
** 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) normal freya run
** iexcit = 1 use hexnb instead of crsecs to get cross sections
** but neglect presence of excited sates in beam
** iexcit = 2 use hexnb with excited beam states allowed.
** inubpat switch controlling 2-d beam deposition calculation
** 0, (default) bypass 2-d deposition
** MUST BE = 0 FOR CQL RUNS
** NOTE: follow namelist variable, eleccomp, must be set in 2nd /setup/
** eleccomp If "enabled", the 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)
**
***************************************************************
**
** 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.
**
** 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 => [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"
** 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)
**
** 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 investigating this case).
** ="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
** nofftran: time-step at which transport is turned off(?)
**
** 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 determnined after s-step
**
** 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
**
********************************************************************
$setup
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!