******************************************************************** ** ** 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!