Description of ONETWO's "inone" Input File Format updated 15 Sep 99
-------------------------------------------------
ONETWO Version ID: 15sep99m
c ----------------------------------------------------------------------
c --- SUMMARY OF * MAJOR * NEW ONETWO VERSION RELEASED LATE JUNE 1997
c ----------------------------------------------------------------------
c 1) TDEM mode of operations (i.e., time-dependent eqdsk mode)
c 2) IFS (Dorland-Kotschenreuther) confinement model
c (full 2d, 1d model implemented for circular plasmas only !!!!!!!)
c 3) New Weiland model (11 equations instead of 6)
c 4) Time-dependent beam power
c 5) Various changes (by Gary Staebler) to the Houlberg Bootstrap Model
c ----------------------------------------------------------------------
c --- SUMMARY OF MAJOR NEW ONETWO VERSION RELEASED MID-MARCH 1996
c ----------------------------------------------------------------------
c
c 1) SFAREA in summary page was changed to be the physical area.
c 2) All new input variables are described in file "cray102.f" as usual.
c 3) An example case is in file "inone.march" for shot 87977.
c 4) An example d-t fusion case is in file "inone.fusion" for shot 87977.
c 5) Use the ONETWO only with the new version of PREPLT and TRPLOT.
c 6) The 65x65 version of ONETWO will now process 33x65 EQDSKs as well.
c 7) The Weiland confinement model is installed but is not yet tested!!!
c 8) Beam-beam, beam-thermal, and thermonuclear d-d and d-t fusion rates
c are calculated with Bosch & Hale cross sections.
c
c ----------------------------------------------------------------------
c
*
c ----------------------------------------------------------------------
c --- GEOMETRY FLAG AND MACHINE NAME
c ----------------------------------------------------------------------
c The following flags appear in the FIRST LINE of the input file "inone",
c in the order given. Free-field format decoding is done to get the names.
c
c% codeid Geometry identification flag
c DEFAULT -> 'onedee' : 1-D run; circular or elliptical flux surfaces
c 'dee' : 1-1/2-D run; dee-shaped flux surfaces
c
c% machinei Machine identification flag
c 'doub-iii' : old Doublet III machine
c DEFAULT -> 'diii-d' : current DIII-D machine
c diffeq_methd =0 default,relaxation method solution
c =1 method of lines (not fully implemented)
c functions qq and qqq (in cray401.f) were renamed fqq and fqqq
c respectively. This was done to eliminate a compiler warning regarding
c the use of qq as both a functio and as an array.
c Appropriate changes were made in cray306.f (subroutine DELSOL)
c
c ----------------------------------------------------------------------
c
c ----------------------------------------------------------------------
c --- DEFINE THE THREE NAMELISTS (NAMELIS1, NAMELIS2, NAMELIS3)
c ----------------------------------------------------------------------
c
namelist /NAMELIS1/
. rmajor, rminor, elong, btor, nprim, nimp, namep, namei, fd,
. itenp, iteni, itte, itti, itxj, taupin, icenez, inenez, zfrac,
. adjzeff, w1saw, w2saw, w3saw, wneo, wneot, w3cla, w1typ, w12typ,
. w13typ, w2typ, w3typ, w4typ, typa, wstd, w1isl, w2isl, w3isl,
. nisl, w1mix, w2mix, bparzeff, w3mix, w4mix, w5mix, qmix, rsmixx,
. tsmix, tdmix, ipmix, w0mix, dtemix, dtimix, fusmix, mixpro,
. s3mix, s71mix, s18mix, trmix, jsxr, jco2, jzeff, jterow, nterow,
. imesh, nj, rnormin, njenp, njeni, njte, njti, njene, njzef,
. njcur, r, dr, rin, rout, zax, ene, enec, eneb, alpene, gamene,
. enp, enpb, enpc, alpenp, gamenp, eni, enib, enic, alpeni,
. gameni, te, teb, tec, alpte, gamte, ti, tib, tic, alpti, gamti,
. curden, xjb, xjc, alpxj, gamxj, zeff, zeffb, zeffc, alpzef,
. gamzef, xmtmdifs, rmtmdifs, bparenp, bparte, nbctim, bctime,
. totcur, iffar, enein, tein, tiin, zeffin, bpareni, bparti,
. time0, timmax, nmax, dt, dtmin, dtevmin, dtmax, relmin, relmax,
. relit, itmax, theta, timav, timprt, timplt, mprt, mplot, prtlst,
. pltlst, curdenin, jflux, jcoef, jsourc, jbal, jtfus, nneu,
. namen, iref, ipcons, gasflx, recyc, twall, wion, wrad, raneut,
. relneu, erneumax, idiagn, nengn, englstn, jprt, rtandn, rhdn,
. widths, nps, enes, tes, ttweak, relaxtyp,diffeq_methd,
. nw1pro, w1pro, nw2pro, w2pro, nw3pro, w3pro, nw4pro, w4pro,
. bl2in, fusnin, bparcur, ticin, voltin, qcin,
. zeflim, bparang, w33min, xdebug, ifred, ilimdt, vtyp, nvpro,
. vpro, ibaloo, voltav, betlim, betlim0, iangrot, itangrot,
. angrotin, rangrot, iwangrot, bparene, etaioff, renpin, reniin,
. rtein, rtiin, renein, rzeffin, rcurdein, splninpt, jgboot,
. jhirsh, implicit_fh, wrebut, relaxsol, relaxrebut, ddebug,
. resistive, ftcalc, qrebsmth, wshay, smult, scsmult, snexp,
. sbpexp, srexp, sbigrexp, stexp, sdtdrexp, srin, srout, suserho,
. skimult, srincev, sroutcev, relaxshay, tohmwant, ifsflag, aeh,
. ael, aih, ail, bh, bl, ch, cl, alfae, alfai, betah, sigma,
. gammah, lsfctr, tirlw, rlw_model, run_preplt, xrot, xeden,
. xsecder, ishayform, skimass, times_rgc, timee_rgc, irgc, rgc,
. rgca, rgcb, rgcc, rgcd, rgce_mult, rgci_mult, rgc_string,
. fusionvb, tportvb, namelistvb, zenvb, exptl_neutron_rate,
. neucgvb, wweiland, include_weiland, cer_ion, tdemvb, fiziksvb,
. squeeze, ikpol, kpolin, rkpol, bparkpol,wneo_elct_ifs,
. include_ifs, dorl_kotch, eps_adaptive, gridgenvb,
. speed_adaptive, freeze_adaptive, curve_eps_adaptive,
. spec_profiles, no_te_convection, no_ti_convection, rho_edge,
. timecrit, dtmaxcrit, fix_edge_te, fix_edge_ti, extend_seval,
. set_chie_chii,testing_NTCC,ce0_mgb,alpha_mgb,ci1_mgb,
. ci2_mgb,ce_bgb,cfe_mgb,cfe_bgb,cfi_mgb,cfi_bgbc1_g,c2_g,c_theta,
. include_itb, include_gf, dorl_kotche,dorl_kotchi,exbmult_ifs ,
. steps_per_plot
c
namelist /NAMELIS2/
. timbplt, beamon, btime, nameb, relnub, tfusbb, anglev, angleh,
. nashape, aheigh, awidth, bcur, bptor, blenp, nbshape, bleni,
. bheigh, bwidth, bhfoc, bvfoc, bhdiv, bvdiv, ebkev, fbcur,
. nbeams, naptr, alen, bvofset, bhofset, nsourc, sfrac1, mf,
. npart, npskip, rpivot, zpivot, ranseed, fionx, iddcal, fdbeam ,
. nbinject, rfmode, rfon, rftime, turnonp, turnonc, rfpow, idamp,
. xec, zec, irfplt, gafsep, freq, rfrad1, rfrad2, wrfo, wrfx,
. rnormin, njqin, rfrad1ic, rfrad2ic, necsmth, qine, qini, xdebug,
. a1rf, a2rf, wrfe, wrfi, nrfzon, rfzone, ichmod, betalm, relrf,
. nprf, iside, xkpar, nhigh, ykperp, navg, thetec, phaiec, hlwec,
. nampel, pelrad, vpel, nbgpel, timpel, inubpat, npat, ratwec,
. nray, ifus, iaslow, wtifus, ibcur, ibcx, ibslow, iborb, iyoka,
. ishot, itime, itrapech, itrapfi, wdelt, wgam, wohm, nqrad,
. qradr, qradin, refrad, ds_tk, fe_tk, ne_tk, iexcit, ilorent,
. mstate, ncont, izstrp, kdene, nbeamtcx, kdeni, kdenz, ksvi,
. ksvz, ksve, krad, ngh, ngl, nouthx, rnp, irfcur, ifb, rfcur,
. lmode, ifbprof, alphaf, hdepsmth, zrffw, nzrffw, pzrffw, htsfw,
. iswchfw, lifw, nihfw, timrfp, freqfw, rnpfw, impath, angrm2d,
. angrcple, nrayptrt, powersrt, nnkparrt, heightrt, maxrefrt,
. islofart, anzinfrt, anzsuprt, nthinrt, nfwsmth,
. nkfcd, xntor, rpant, gamloss, iterate_beam,
. iddfusrate, iddfusb, iddfusb_s, ddfusb_t, iddfusb_bulk,
. icalc_cxfactor, imaxbeamconvg, beam_beam_fusion,
. extcurrf, extqerf, extqirf, extcurrf_id, extqerf_id, extqirf_id,
. extcurrf_amps, extqerf_watts, extqirf_watts, extcurrf_rho,
. extqerf_rho, extqirf_rho, extcurrf_curr, extqerf_qe, extqirf_qi,
. extqerf_nj, extqirf_nj, extcurrf_nj, freyavb,
. beam_thermal_fusion, fast_ion_target, rtstcx, relaxden,
. relaxden_err,npart_mcgo,mcgo_output_file, use_Callen,
. neg_ion_source
c
namelist /NAMELIS3/
. ifixshap, mhdmode, xdim, ydim, redge, nlimiter, xlimiter,
. ylimiter, ifill, greentab, isym, mhdmethd, fixfcoil, ifitpsi,
. ifitprob, iecurr, ivessel, timeqbcd, flxeqbcd, curmax, toleq,
. tolcur, iteq, itcur, omeq, omcur, npsi, ieqmax, deqlst, delcap,
. delrho, minchisq, maxitr, ibypas, dtmine, ecurrt, fcoilcur,
. vesscur, irguess, eqdskin, pcurmhd, btormhd, vloopmhd, expmp2,
. volaray, voladj, volnudge, itvol, tolvol, limpos, dvadjmax,
. ieqprt, j2prt, iprtit, ieqdsk, rf_output, ispare, aspare,
. xdebug, ddebug, rmagax, zmagax, rscale, zscale, mhdonly,
. fitfcur, errpsilp, errmprbe, errfcoil, errecoil, errvescr,
. ivertsbl, rvloop, zvloop, ltest_code, optwi, optwe, optomegai,
. optomegae, tstark, sigstark, fwtstark, rstark, zstark, a1stark,
. a2stark, a3stark, a4stark, use_stark, timestark, xi_include,
. psifctr, tensionspl, mresidual, nfitpoints, splinefitbd,
. fitboundary, derwght, adotmult, f2d1mult, f2d2mult, f2d3mult,
. use_cnt1, use_cnt2, iterdb, iterdsc, itorfluse, irwflag,
. comp_methd_eqdsk, use_efit_cntr, curtype, pol_flux_lim,use_Bp0,
. cap_mult,deltat_fixed_boundary
*
c **********************************************************************
c *** FIRST NAMELIST (NAMELIS1) ****************************************
c **********************************************************************
c
c ----------------------------------------------------------------------
c --- MACHINE PARAMETERS
c ----------------------------------------------------------------------
c% rmajor Major radius (cm) to magnetic axis in 1-D runs.
c In 1-1/2-D runs this radius is used only to specify btor.
c Moreover, if irguess < 0 and rmajor = 0, then rmajor is
c read from the eqdsk file.
c% rminor Horizontal minor radius (cm) in 1-D runs (recalculated in
c 1-1/2-D runs)
c% elong(m) Elongation (height-to-width ratio) of elliptical flux
c surfaces in 1-D runs. If nbctim > 1, elong(2) .ne. 0
c specifies that this quantity is time-dependent (rminor
c is held fixed). See below on nbctim and bctime(m).
c% rin Inside major radius (cm) of vacuum vessel used to determine
c whether neutral beams or SXR chords are reentrant
c% rout Outside major radius (cm) of vacuum vessel used only by
c neutral beam plotting code NUBPLT
c% zax Vertical position (cm) of magnetic axis in 1-D runs
c (recalculated in 1-1/2-D runs)
c% btor Toroidal magnetic field (Gauss) at rmajor. In 1-1/2-D
c runs if irguess < 0 and btor = 0.0, then btor is read from
c the eqdsk file. if btor = -1.0e30 on input then btor is
c obtained from the third namelist of inone from the vector
c btormhd, in a fashion similar to the treatment of totcur
c (see totcur description). (btor is nominally time-independent
c for normal discharges. We allow time-dependent input
c for generality using btormhd).
c% tportvb integer: if set to 1, gives some detailed output to screen
c (primarily for diagnostic use)
c
c% namelistvb an integer, if set to 1, will print a message to the
c screen as each namelist is read. use this to isolate
c undefined variables in the namelists
c% zenvb an integer if set to 1 will print first few elements of
c subroutine ZEN (i.e., density) calculations to screen
c% tdemvb an integer if set to 1 will cause dump_data to be called
c (used for debugging)
c% fiziksvb an integer if set to 1 will cause FIZIKS to dump some info
c to be called (used for debugging)
c ----------------------------------------------------------------------
c
*
c ----------------------------------------------------------------------
c --- ION PARAMETERS
c ----------------------------------------------------------------------
c% nprim Number of primary ion species (1 to 3)
c% nimp Number of impurity ion species (0 to 2)
c% namep(i) Name of ith primary ion species
c 'h' , protons
c 'd' , deuterons
c 't' , tritons
c 'dt', mixture of d and t
c 'he', thermal alphas
c% namei(i) Name of ith impurity ion species
c 'he', helium
c 'c' , carbon
c 'o' , oxygen
c 'si', silicon
c 'ar', argon
c 'ti', titanium
c 'cr', chromium
c 'fe', iron
c 'ni', nickel
c 'kr', krypton
c 'mo', molybdenum
c 'w' , tungsten
c% fd Number fraction of deuterons in d-t mixture
c ----------------------------------------------------------------------
c
do i=1,kprim
namep(i) = ' '
end do
do i=1,kimp
namei(i) = ' '
end do
nprim = 1
nimp = 0
namep(1) = 'h'
fd = 0.5
tportvb = 0 ! controls diagnostic output
fiziksvb = 0
namelistvb = 0
zenvb = 0
neucgvb = 0
tdemvb = 0
*
c ----------------------------------------------------------------------
c --- TRANSPORT FLAGS
c ----------------------------------------------------------------------
c% testing_NTCC special flag used to compare solutios with NTCC code
c (intended for use by developers only)
c% itenp(i) 1, transport of primary ion species i is considered,
c i.e., profile is calculated
c 0, transport of primary ion species i is not considered,
c i.e., profile is specified
c% itte Same as above for electron temperature
c% no_te_convection Set to 1 to turn OFF convective ELECTRON energy flux
c
c% itti Same as above for ion temperature
c% no_ti_convection Set to 1 to turn OFF convective ION energy flux
c% itxj Same as above for current density
c% iangrot switch which includes ( = 1) or excludes (=0) the toroidal
c angular momentum. If iangrot = 0 (default), the code runs
c in its old mode of operation (i.e., toroidal angular
c momentum effects are totally absent). if iangrot = 1, then
c toroidal angular momentum effects are to be included,
c depending on the setting of switch itangrot. (the purpose
c of switch iangrot is to allow automatic setting of
c compatibility with older versions of the code which do not
c include toroidal rotation. As a consequence, two switches
c are required to include toroidal rotation).
c% itangrot switch which is used only if iangrot = 1. If iangrot=1 and
c itangrot = 1 the toroidal angular momentum equation is
c included in the set of transport equations to be solved.
c (i.e., the simulation mode).
c if iangrot = 1 and itangrot=0 the toroidal angular momentum
c is not transported (the initial angular rotation speed
c profile is independent of time,or input as a function of
c time, (i.e., the analysis mode).
c The source terms due to the presence of a nonzero rotation
c speed are incorporated in the remaining equations (the
c sources can be selectively turned off in both analysis and
c simulation runs-see below).Only first order corrections
c in rotation speed / ion thermal spee are included.
c default is itangrot = 1.0
c% ikpol If ikpol = 1, spline data for kpolin is expected
c (poloidal velocity of CER ion/Bp) (default = 0)
c% nbeamtcx switch related to torque calculation. should be set if
c angular rotation is used. see description under beam
c input (note nbeamtcx is input in the second namelist).
c ----------------------------------------------------------------------
c
do index=1,kprim
itenp(index) = 1
end do
c
testing_NTCC = 0 ! no tests are performed
itte = 1
itti = 1
itxj = 1
iangrot = 0
itangrot = 1
implicit_fh = .false.
no_te_convection = 0
no_ti_convection = 0
*
c ----------------------------------------------------------------------
c --- ANALYSIS MODE PARAMETERS
c ----------------------------------------------------------------------
c
c% taupin Particle confinement time (s) for each primary species
c with itenp(i) = 0 (same as global electron confinement time)
c% ttweak Time step (s) between adjustments (tweaks) to obtain
c a specified value for one of several parameters
c% fusnin Specified fusion neutron rate (s-1) obtained by
c adjusting wneo(3,3) or w3typ (w3typ>0 takes precedence
c over wneo(3,3); fusnin>0 takes precedence over ticin);
c see also the neutral beam parameters iddcal and fdbeam in
c the second NAMELIST
c% ticin Specified central ion temperature (keV) obtained by
c adjusting wneo(3,3) or w3typ (w3typ>0 takes precedence
c over wneo(3,3))
c% voltin Specified one-turn voltage (V) obtained by adjusting
c zeffc and zeffb
c% voltav Specified average voltage (V) across plasma obtained by
c adjusting zeffc and zeffb (voltav>0 takes precedence
c over voltin)
c% qcin Specified q on axis obtained by adjusting zeffc/zeffb
c% w33min Lower limit allowed for wneo(3,3) or w3typ
c% zeflim Upper limit allowed for zeff (if non-zero)
c% tohmwant Total desired ohmic current in steady state (can be negative)
c tohmwant is used only in subroutine SSCURDRV to determine
c required driven current in steady state
c ----------------------------------------------------------------------
c
tohmwant = -1.0e+100 ! turns off call to SSCURDRV
ttweak = 0.05
taupin = 1.0
fusnin = 0.0
ticin = 0.0
voltav = 0.0
voltin = 0.0
qcin = 0.0
w33min = 1.0
zeflim = 0.0
*
c ----------------------------------------------------------------------
c --- FLAG FOR TOGGLING THE RUNNING OF PREPLT
c ----------------------------------------------------------------------
c
c% run_preplt = .true. => DO call PREPLT postprocessor right before
c the end of ONETWO's execution, so that TRPLOT
c can be run subsequently. This is the DEFAULT.
c
c run_preplt = .false. => DO NOT call PREPLT postprocessor right before
c the end of ONETWO's execution.
c If you later decide that you want to run TRPLOT,
c remember to first run (standalone) PREPLT, since
c it would not have already been called by ONETWO.
c
run_preplt = .true. ! PREPLT will be run at the end of the ONETWO run
*
c ----------------------------------------------------------------------
c --- PARAMETERS FOR BASIC TRANSPORT MODELS
c ----------------------------------------------------------------------
c% wneo
c ((wneo(i,k), i=1,5), Weights for neoclassical transport
c k=1,5) if iangrot = 0 (i.e., toroidal momentum effects
c are absent), wneo is input exactly as
c before. the increase in dimension from
c 4 to 5 is accounted for in the code.
c if iangrot = 1 wneo must be input as a 5x5
c array! This means adding a trailing
c number to each row of the old 4 row wneo
c array and adding a 5th row of 5 elements.
c Note that this scheme will allow input of
c indicies to the wneo array by assigning
c each input value only if iangrot = 1. if
c iangrot = 0, statements of the type
c wneo(3,2) = 5, for example, will not work!
c for iangrot = 0, the entire 4x4 wneo array
c must be input. indicies assigned to the
c elements of the array must be obtained
c through the standard Fortran default.
c
c% jhirsh Selection Switch for Bootstrap Current Model
c = 0 (Default) The calculation is based on Hinton & Hazeltine
c (Rev. Mod. Phys. 48 (1976) 239).
c = 1 The bootstrap current is modeled using the formulae of
c Hirshman (Phys. Fluids 21 (1978) 1295).
c = 2 Similar to jhirsh = 1 but includes the fast ions in
c with the thermal ions.
c NOTE: jhirsh = 1,2 may modify the bootstrap
c current calculations near the magnetic
c axis by the poloidal field gyroradius.
c if this modification is not wanted use
c jhirsh = 11 instead of jhirsh = 1 and
c jhirsh = 22 instead of jhirsh = 2 ... HSJ
c = 88 The small collisionality, arbitrary aspect ratio,
c single ion fluid model of Hirshman 1988.
c (Phys. Fluids, vol 31 (10) 1988)
c No corrections for the behavior near the magnetic
c axis due to finite poloidal gyroradius are done.
c = 89 Similar to jhirsh = 88 but includes the fast-ion density
c gradient, assumed to be part of the single-ion fluid.
c = 95 The arbitrary collisionality, arbitrary aspect ratio,
c and multiple ion species model of Houlberg, 1995
c (Nuclear Fusion, to be published).
c Considered the most comprehensive model.
c No provision has been made to account for the fast ion
c bootstrap (Electrons from fast ions and alphas are
c included.)
c = 96 Houlberg model, including a fast ion pressure term for
c the beam and alpha particles. The beam and alpha
c pressures are derived from the averaged stored energy
c density (keV/cm3) of the beam and alphas, i.e.:
c press = nkT = 0.666667*W
c See subroutine NCLBOOT for more information.
c NOTE: Both jhirsh = 95 and 96 use the formula of Y. R. Lin-Liu
c R. L. Miller for the trapped fraction calculation
c unless ftcalc = 'exact' is set.
c
c% resistive character variable selects the resistivity
c to use according to:
c = 'hinton' old (small inverse aspect ratio)
c 'hinton' is the DEFAULT, valid for eps < 0.05
c = 'hirshman' (arbitrary aspect ratio)
c = 'kim' (arbitrary aspect ratio,
c valid in banana regime only)
c 'kim' is NOT CURRENTLY IMPLEMENTED
c
c% ftcalc character variable,
c ftcalc = 'exact' means evaluate the integrals
c along the psi contours to determine the
c non collisional trapped particle fraction.
c ftcalc = 'analytic' means use the (circular)
c inverse aspect ratio expansion:
c ft = 1.0-(1-delta)**2/(SQRT (1.0-delta**2)/
c (1.0 + 1.46 * SQRT (delta))
c where delta is horizontal inverse aspect ratio.
c The formula of Lin-Lui and Miller is used if
c jhirsh = 95 or 96
c ftcalc = 'analytic' is the DEFAULT
c NOTE: this affects bootstrap current calculations
c for models with jhirsh > 0 and also
c the plasma resistivity if resistive = 'hirshman'
c
c% cer_ion identifies the ion species which the charge
c exchange recombination (CER) diagnostic measured.
c The input angrot is then assumed to be Vtor/R along
c the Z=0 cord for the cer_ion species. This is used
c to compute the radial electric field if jhirsh > 95
c The poloidal and toroidal rotation of the main and
c cer ions is output to the netCDF file iterdb.nc .
c cer_ion can take the same values as namep, namei.
c The default = ' '
c
c% squeeze If .true. and cer_ion is set then ion orbit
c squeezing is applied to the NCLBOOT calculations.
c default is .false.
c
c% wneot Weight for neoclassical enhancement of
c resistivity due to trapped electrons
c% w1typ(nbctim), Weights for empirical ("take-your-pick")
c% w2typ(nbctim), transport. The diagonal elements may be
c% w3typ(nbctim), time-dependent if w*typ(2) .ne. 0 and
c% w4typ(nbctim) nbctim > 1 (see nbctim below)
c% w12typ, w13typ Off-diagonal empirical weights
c% vtyp(nbctim) Edge velocity (cm/s) for empirical pinch
c% typa(i), i=1,16 Exponents for empirical transport
c typa(i), i=9,16:
c typa(9) and typa(10) are exponents for the
c total or local integrated heating power,
c p, in Mw. Whether or not p is the local
c heating (in an annular volume region of
c width drho) or the global heating for the
c entire torus is determined by typa(16).
c if typa(16) .ne. 0, the global value is
c selected. if typa(16) = 0, the local value
c is selected. p is defined as,
c p = typa(11)*pohm +
c typa(12)*pbeame +
c typa(13)*pbeami +
c typa(14)*prfe +
c typa(15)*prfi
c the multiplier of the transport
c coefficient is p**typa(9) or p**typa(10),
c depending on the value of,
c pheat = pohm+pbeame+pbeami+prfe+prfi
c if pheat = pohm, typa( 9) is used.
c if pheat > pohm, typa(10) is used.
c no allowance for fusion currently exists.
c note that the factors pohm,pbeame,...,
c appearing above are local (typa(16) = 0.0)
c or global (typa(16) > 0.0). In effect,
c typa(9) and typa(10) are time-dependent
c exponents. This may cause some problems
c (oscillations?) though none are expected.
c since typa(9),typa(10) do not depend on
c the dependent variables (ni,te,..etc.),
c there will be a one time step lag in
c switching from typa(9) to typa(10) since
c the transport coefficients are calculated
c before the powers are known. Again, this
c is not expected to cause difficulties.
c This argument also applies at the initial
c time step (pohm and/or paux is zero at the
c initial time when the transport
c coefficients are set up).
c% nw1pro, w1pro Number of points (nw1pro) and scale factor
c profile (w1pro) for 'typ' particle diffusion,
c which may be time-dependent.
c% nvpro , vpro Number of points (nvpro) and scale factor
c profile (vpro) for 'typ' pinch velocity,
c which may be time-dependent.
c% nw2pro, w2pro Number of points (nw2pro) and scale factor
c profile (w2pro) for 'typ' electron thermal
c conductivity, which may be time-dependent,
c analogous to nqrad and qradin.
c% nw3pro, w3pro Number of points (nw3pro) and scale factor
c profile (w3pro) for 'typ' ion thermal
c conductivity, which may be time-dependent.
c% nw4pro, w4pro Number of points (nw4pro) and scale factor
c profile (w4pro) for 'typ' momentum diffusivity
c in units of 10**4 cm**2/sec
c% relaxtyp under-relaxation parameter for take-your-pick model
c during the iteration of the corrector steps the Rebut
c diffusivity will be under-relaxed if this parameter
c is < 1. set to 1 for no relaxation. under-relaxation
c may be necessary to get smooth solution.
c see also paramater relaxsol and ddebug(3).
c% wstd Weight for standard drift-ballooning model.
c% betlim, betlim0 An input limit on beta. When the total plasma
c beta passes beyond betlim, the 'typ'
c transport is increased by the factor
c (1.0 + (beta-betlim)/betlim*betlim0).
c betlim0 give the strength of the beta-limit
c turn-on.
c
c% set_chie_chii Switch added 6/25/98.
c This is a more general and accurate way
c of settig chie = const*xchi than using the wneo(2,2)
c value. If set to a non zero value then
c the electron chi will have added to it the term
c set_chie_chii*chi. The total electron diffusivity
c will be the sum of all active electron chi models as
c usual.
c
c% iwangrot angular momentum diffusivity model switch
c (ignored in analysis mode)
c = -2 (Mattor-Diamond model):
c Note this model predicts momentum
c diffusivity = ion thermal diffusivity
c and also gives a value for the electron
c thermal diffusivity. At present the ion
c and electron thermal diffusivities are
c calculated and printed out, but not used
c otherwise. Ref: Phys Fluids 31(1988)1180.
c The Mattor-Diamond model,as implemented here
c requires the specification of some multipliers
c as follows:
c% xmtmdifs(i) i = 1,nprim (nprim is number of primary ions)
c xmtmdifs(i) is an arbitrary multiplier of
c chi mometum (and chi ion thermal diffus.)
c for primary ion species i. That is,each
c primary ion species yields a chi given
c by eq. 43 (above ref.) A weighted sum of
c these individual diffusivites,with weight
c xmtmdifs(i),is used as the Mattor-
c Diamond momentum diffusivity.
c xmtmdifs(nprim+1) An impurity enhancement factor.
c According to ref
c This factor should multiply the
c expression
c to account for impurities. Normally
c this factor would be set to 1.0
c xmtmdifs(nprim+2) A threshold value for the etai modes.
c If etai is less than xmtmdifs(nprim+2)
c then the model gives zero diffusivity.
c Note: at present there is no soft turn-
c on of the etai modes in this model. This
c can cause oscillations in the angular
c momentum but does not affect the thermal
c ion and electron energy equations directly
c (there is some coupling through the ion
c energy equation but this should not
c cause any problems. If oscillations do
c show up a soft turn on model will have
c to be developed in subroutine DIFFUS.)
c xmtmdifs(nprim+3) A multiplier for the neoclassical
c momentum diffusivity which will be
c added to the Mattor-Diamond diffusivity.
c Thus for example we can set xmtmdifs(
c nprim+3) = 1.0 and xmtmdifs(i) = 0.0,i=
c 1,..nprim, to get the neoclassical
c diffusivity with the switch iwangrot
c set equal to 2. (of course we would
c normally achieve this more directly
c by setting iwangrot = 0 to begin with)
c Note that the neoclassical diffusivity
c is always multiplied by wneo(5,5) so
c that the actuall multiplier becomes
c wneo(5,5)*xmtmdifs(nprim+3).
c xmtmdifs(nprim+4) Multiplier for electron thermal
c diffusivity calculated from eq. (51)
c (above reference). Not significant at
c this time since the electron (and ion)
c thermal diffusivites are not used in
c any calculations.(printed out only)
c
c THE FOLLOWING MODELS APPLY ONLY TO MOMENTUM DIFFUSIVITY
c parabolic diffusivity:
c iwangrot = -1
c the model is defined by
c xkangrot(j) = (xmtmdifs(1)-xmtmdifs(2))*(1.0-
c roa(j)**xmtmdifs(3))**xmtmdifs(4)-xmtmdifs(2)
c
c iwangrot = 0 (default)
C the model is neoclassical (see subroutine DIFFUS)
c
c iwangrot = 1
c the angular momemtum diffusivity is constant
c xkangrot(j) = xmtmdifs(1) for all j
c
c iwangrot = 2
c the model is linear
c xkangrot(j) = xmtmdifs(1)*roa(j)+xmtmdifs(2)
c
c iwangrot> = 3
c the model is a spline
c in this case a minimum of three knots must be
c specified in the array rmtmdifs. The first knot
c must be at 0.0 and the last knot must be at 1.0
c
c NOTE the units of xkangrot are assumed to be
c cm**2/sec so set xmtmdifs accordingly.
c% xmtmdifs
c% rmtmdifs input arrays that define the angular momentum
c diffusivity model (see iwangrot above).
c
c --------------------- Rebut-Lallia-Watkins Model ---------------------
c
c% wrebut weight for Rebut-Lallia-Watkins diffusivities
c and diffusion coeffients. set wrebut = 0 (default)
c if this option is not wanted. set wrebut = 1.0 to
c get full Rebut model. wrebut multiplies the model
c so an arbitrary part can be added.
c% relaxrebut under-relaxation parameter for Rebut-Lallia-Watkins
c model. during the iteration of the
c corrector steps the Rebut diffusivity will be
c under relaxed if this parameter is less than 1.
c set to 1 for no relaxation. under relaxation
c may be necessary to get smooth solution.
c see also paramater relaxsol and ddebug(3).
c note that relaxrebut refers only to the
c anomalous part of the RLW diffusivity.
c if it is desired to relax (in addition to,
c or in place of, just the anomalous RLW value),
c then use ddebug(3).
c% qrebsmth = 0 no effect,
c = 1 smooth the q profile used in Rebut-Lallia model
c (note that q is not smoothed elsewhere in the code)
c
c% tirlw = 0 (default) means use the real ion temperature,
c ti, for the evaluation of the RLW model.
c = 1 use an effective ti, defined as
c
c ti*nth + 0.66*(wbeam+walp)
c ti = -----------------------------
c (nth+nbeam+nalp)
c
c where nth is the total thermal ion density
c (includes impurities)
c nbeam is the fast ion density due to beams
c nalp is the fast alpha particle density
c wbeam is the fast ion stored energy density
c due to beams
c walp is the stored fast alpha particle energy density
c
c% rlw_model = 'old' ==> original 1988 RLW model
c = 'new' ==> with corrections for ion rho-star scaling (DEFAULT)
c
c REFERENCES:
c
c Rebut, P.H., Lallia, P.P., Watkins, M.L., Plasma Physics
c and Controlled Nuclear Fusion Research, (Proc. 12th Int.
c Conf., Nice, 1988), Vol. 2, IAEA, Vienna (1989) 191
c
c Rebut, P.H., Watkins, M.L, Gambier, D.J., Boucher, D.,
c Physics of Fluids b 3 (1991) 2209
c
c Boucher, D., personal communication
c
c For the HSJ 8/22/95 update (rlw_model = 'new'),
c the reference is Seville, IAEA, 1994, Rosenbluth et al.
c
c -------------------- HSIEH ("SHAY") MODEL OF CONDUCTIVITY ------------
c
c% wshay multiplier for Hsieh model of ke. set wshay = 0.0
c to not use this model ,set wshay = 1.0 to use full
c Hsieh model (in addition to any other models that
c are turned on).
c% smult constant in Hsieh model (units of this constant
c depend on exponents choosen for Hsieh model.
c% scsmult logical,if true constant multiplier for Hsieh
c model, smult, is to be found from power balance
c analysis (using least squares method)
c NOTE that if scsmult is set then
c itte=0
c itti=0
c iterate_beam = .true.
c is internally set by the code.
c% snexp exponent of electron density in Hsieh model
c% sbpexp bp
c% srexp r
c% sbigrexp R
c% stexp Te
c% sdtdrexp d(Te)/dr
c% srin starting value of NORMALIZED rho at which Hsieh model
c becomes active
c% srout ending value of NORMALIZED rho at which Hsieh is active
c% srincev these parameters are similar to srin and srout
c% sroutcev they are used only if scsmult = true. in this
c case srincev is the starting value of NORMALIZED
c rho at which data for the least squares evaluation
c will be used. Similarly srout is the ending value
c of NORMALIZED rho at which data for the least
c squares evaluation will be used. It is required
c that srincev .ge. srin and sroutcev .le. srout.
c Also we must have srincev .le. sroutcev. If
c srincev = sroutcev then the single nearest point
c in NORMALIZED rho will be used. In this case
c a least squares anaylis is not possible because
c there are not enough degrees of freedom.
c% suserho logical, if true, then use flux surface average form
c of Hsieh model, including (grad-rho)**2 term.
c% skimult set ion k =skimult * shay (ke)
c (on top of any other ion ki that is active).
c of Hsieh model, including (grad-rho)**2 term.
c
c ================================
c 1995 Upgrade of the Hsieh Model:
c ================================
c
c% ishayform = 1 Use the dimensionally correct version of the model
c = 0 Use the old version with funky units
c NOTE: The dimensionally correct form is valid only for the
c default exponents
c% skimass Since the dimensionally correct version of the Hsieh
c (amu) model for electron thermal conductivity now has a
c 1/SQRT(elec mass) dependence, the ion multiplier
c SKIMULT should have a similar dependence on the ion
c mass. This is accomplished by using an effective
c atomic mass for the ions, in a.m.u. (The ratio
c sqrt(me/mp) is assumed to be included in the value
c of SKIMULT)
c For a deuterium plasma, use skimass= 2.0 (default),
c a hydrogen plasma use 1.0, etc. Its up to you to
c decide what to use in a multiple-species plasma!
c NOTE: skimass is auotmatically reset to 1.0 for
c ISHAYFORM = 0 to keep old form of model
c -Daniel Finkenthal
c
c ------------------------------------------------------------------------
c ---------------------------- WEILAND MODEL -----------------------------
c This model is described in:
c
c ... J. Weiland and H. Nordman, "Drift wave model for inward
c ... energy transport in tokamak plasmas," Institute for
c ... Electromagnetic Field Theory and Plasma Physics,
c ... Gothenburg, Sweden, (1992) CTH-IEFT/PP-1992-13 ISSN.
c
c ... J. Weiland, A.B. Jarm\'{e}n, and H. Nordman, "Diffusive
c ... particle and heat pinch effects in toroidal plasmas,"
c ... Nucl. Fusion {\bf 29} (1989) 1810--1814.
c
c ... and many more (see notes of Batemen)
c
c ONETWO INPUT SWITCES FOR THE WEILAND MODEL
c ------------------------------------------------------------------------
c% include_weiland = 1 must be set (in BOTH analysis and simulation
c mode) in order to activate this model.
c
c the Weiland particle and energy diffusivities will be used
c in simulation mode for any dependent variable for which the
c simulation mode flag is set (i.e., itenp, itte, itti) if
c% wweiland .ne. 0.0
c in this case the term
c wweiland*(d and/or chie and/or chi weiland)
c WILL BE ADDED TO any other models that are turned on.
c (So make sure other models are turned off )
c
c the switch relaxrebut is used for the Weiland model as well
c see comments under the RLW model for meaning of relaxrebut.
c
c
c In ANALYSIS MODE,or if wweiland =0.0, the d,chie,and chii for the
c Weiland model are calculated,printed and plotted for inspection
c (provided that include_weiland=1 of course) BUT ARE NOT USED
c IN ANY OTHER WAY.
c
c ----------------------------------------------------------HSJ/2/16/96-----
c
c OCTOBER 1995 MODIFICATION OF ELECTRON AND ION THERMAL CONDUCTIVITIES
c DUE TO A PHENOMENOLOGICAL CRITICAL GRADIENT IN THE TOROIDAL ROTATION
c SPEED PROFILE
c RGC (Rotational-Gradient-Critical) Model
c also known as the "DeBoo Phenomenological Model"
c
c This model multiplies the anomalous part of the electron and/or ion
c thermal conductivity by the factor
c rgca
c rgc_mult(rho)= ----------------------------------
c rgcb+rgcc*[abs(rg)/abs(rgc)]**rgcd
c
c Here rg and rgc are gradients (derivatives wrt rho, not normalized rho)
c and rgca, rgcb, rgcc, rgcd are input constants.
c abs(rg(rho)) is taken as the maximum of abs(rg(rho)) and abs(rgc(rho))
c so that their ratio is never less than 1. To reduce rgc_mult to 1
c you will typically want to take rgca=rgcb+rgcc, but that's up to you.
c
c The idea behind using this model is that rgc(rho) is determined a priori
c over some time interval (i.e., H mode). Then the model is applied over
c a different time interval (i.e., VH mode), or a different case with the
c now known rgc(rho) factor. Because we cannot switch from analysis to
c simulation mode during a onetwo run this typicaly means that you must do
c two runs, one to determine rgc, which would typically be done in
c analysis mode, and another run to use rgc in simulation mode.
c There is an exception to this rule because rgc is determined a priori
c (i.e., before we take any time steps away from time0).
c Hence you can set up your input
c file with kinetic profiles that are given for times before time0,
c and of course at least all the way up to timmax (use the bctimes
c vector to specify the times involved).
c Using times_rgc and timee_rgc you then select over what time interval
c you want rgc(rho) averaged. Typically you would not want the interval
c (times_rgc, timee_rgc) to overlap with (time0, timmax) but again that's
c up to you.
c
c% irgc turns on/off the model as follows:
c irgc=0 (default) model not used
c irgc=-1 don't apply the model,just determine
c the critical gradient factor rgc(i),i=1,2..nj
c and print it out in outone {times_rgc and
c timee_rgc must be set within bctime(1) to bctime(nbctim)}
c irgc=+1 Apply the critical gradient factor
c to the electron and ion thermal conductivity
c selected. If one or both times_rgc and timee_rgc are set to
c times outside the interval bctime(1) to bctime(nbctim)
c it is taken to mean that rgc(rho) is supplied in the inone
c file and should not be calculated.
c If times_rgc and timee_rgc are such
c that they form a subset of bctime(1),bctime(nbctim)
c then rgc will be calculated first, before the simulation
c is performed.
c (irgc is independent of ifsflag, so you could have
c both flow shear suppression effects on)
c
c% rgce_mult take electron conductivity to be
c xketotal=xkeano*rgce_mult*rgc_mult+xkeneo
c% rgci_mult take ion conductivity to be
c xkitotal=xkiano*rgci_mult*rgc_mult +xkineo
c here xkeano and xkiano are the anomalous part
c of the thermal conductivity (i.e., HSIEH, RLW, etc.)
c
c AT PRESENT RGC_MULT
c IS WIRED UP TO THE (ANOMALOUS PART)
c OF THE HSIEH, RLW, AND TAKE-YOUR-PICK MODELS!!
c% rgca
c% rgcb
c% rgcc
c% rgcd input factors see above
c% rgc(i) i=1,2..nj, see explanation above. NOTE if rgc is input
c in inone teh gradient should be wrt normalized rho.
c% times_rgc start time for determination of critical gradient
c factor rgc (sec)
c% timee_rgc end time for determination of critical gradient
c factor rgc (sec)
c YOU MUST SUPPLY APPROPRIATE EMPIRICAL TOROIDAL ROTATION
c SPEED PROFILES (SEE DESCRIPTION OF ANGROTIN...)
c THAT COVER THE RANGE OF TIMES SELECTED BY times_rgc
c and timee_rgc in order to determine rgc!
c
c% rgc_string character variable (up to 128 characters in length)
c will be echoed on output listing exactly as input
c
c ------------------------------------------------------------------------
c ------------------------------------------------------------------------
c ---------------------------- IFS MODEL ---------------------------------
c This model is described in:
c
c ... Dorland et. al. IFSR#677 "Comparison of Nonlinear
c Toroidal Turbulence Simulations with Experiments"
c
c Kotchenreuther et al. "Quantitative Predictions of Tokamak
c Energy Confinment From First Principles Simulations with
c Kinetic Effects." Phys. Plasmas 2,(6),june 1995
c
c ONETWO INPUT SWITCES FOR THE IFS MODEL
c
c (THIS MODEL IS INCOMPLETE AT THIS TIME, FOR EXAMPLE IT WILL
c NOT!!!! HANDLE NEGATIVE SHEAR . WE WILL HAVE TO AWAIT THE
c FUTURE DEVELOPMENTS FROM DORLAND AND KOTCHENREUTHER)
c Full 2d model is implemented. For 1d however we are currently
c limited to cirular plasmas ONLY!!
c
c ------------------------------------------------------------------------
c% include_ifs = 1 must be set (in BOTH analysis and simulation modes)
c in order to activate this model (default = 0)
c the IFS energy diffusivities will be used
c in simulation mode for any dependent variable for which the
c simulation mode flag is set (i.e., itte, itti) if
c flow shear suppression with the Staebler-Hinton model is possible in
c the IFS model, see switch ifsflg
c
c% dorl_kotch .ne. 0.0 (default = 0.0)
c in this case the term
c dorl_kotch*( chie and/or chi IFS )
c WILL BE ADDED TO any other models that are turned on.
c (So make sure other models are turned off )
c
c the switch relaxrebut is used for the IFS model as well
c see comments under the RLW model for meaning of relaxrebut.
c separate multipliers for the electron and ion channels are now (12/16/98)
c available :
c% dorl_kotche .ne. 0.0 (default = 0.0) electrons
c% dorl_kotchi .ne. 0.0 (default = 0.0) ions
c if dorl_kotch is set to a non zero value the code will set
c dorl_kotche = dorl_kotch
c dorl_kotchi = dorl_kotch
c hence to use dorl_kotche and dorl_kotchi separately make sure
c dorl_kotch = 0 ( which is the default)
c note that dorlkotchi is used as a multiplier for ion and momentum
c diffusivities as well.
c
c% exbmult_ifs E X B flow shear suppression multiplier(exclusive to
c the IFS model). This value is
c defaulted to 0.0 meaning flow shear sup. is turned off.
c set to 1.0 to get the full effect, etc.
c note that cer_ion must be set also for this to work
c use jhirsh =95,96 with this option
c
c
c% wneo_elct_ifs normally the IFS electron diffusivity is assumed
c to be sitting on top of wneo(2,2)*xchie_neo
c set wneo_elct_ifs to make the ELECTRON diffusivity
c sit on top of wneo_elct_ifs*xchii_neo
c In ANALYSIS MODE, the chie, and chii for the
c IFS model are calculated,printed and plotted for inspection
c (provided that include_ifs =1 of course) BUT ARE NOT USED
c IN ANY OTHER WAY.
c
c ------------------------------------------------- HSJ/2/10/97 --------
c
c% ddebug(i) i = 1,2,..50
c these are special switches used to turn on
c various models. We will eventually document
c all these switches as the need arises.
c The original programmer did not supply any
c info on these switches so the following is
c my interpretation of the ones I have used --- HSJ
c
c ddebug(1) unknown (not used?)
c ddebug(2) sets the value of jsteps (note floating
c to integer conversion). The diffusion
c coefficients d(i,k,j), i,k = 1,2..nion+2,
c i,k = nion+4,j=1,2,..nj,are averaged over
c 2*jsteps+1 spatial intervals. the rbp equation,
c nion+3, is not averaged. (the average can be
c just jsteps forwards or backwards, see ddebug(10))
c ddebug(3) a relaxation parameter for d(i,j,k). Normally
c you would underrelax d so 0 .lt. ddebug(3)
c .le. 1 makes sense here. default is 1.0
c (1.0 means no relaxation, 0.0 means take
c all of the old d and none of the new d).
c NOTE: ddebug(2) and ddebug(3) are also in effect
c for qexch, the anomalous energy exchange term
c between electrons and ions. qexch is zero
c unless wstd, the drift ballooning model, is
c turned on.
c ddebug(10) has an effect only if ddebug(2) is set.
c in this case the following is done:
c a) if ddebug(10) = 0.0 (default)
c then the diffusion coefficient is averaged
c symmetrically about each grid point
c b) if ddebug(10) .gt. 0.0
c then the diffusion coeffient is averaged
c only in the direction of increasing rho
c c) if ddebug(10) .lt. 0.0
c then the diffusion coeffient is averaged
c only in the direction of decreasing rho
c ----------------------------------------------------------------------
c ------------------------------------------------------------------------
c ------------------------------------------------------------------------
c --------- modified Bohm-gyro Bohm internal transport barrier model
c
c% include_itb set to 1 to turn this model on
c multipliers see the documentation for definitions.
c ce0_mgb,alpha_mgb,ci1_mgb,
c ci2_mgb,ce_bgb,cfe_mgb,cfe_bgb,
c cfi_mgb,cfi_bgbc1_g,c2_g,c_theta,
c ---------------------------- GLF23 MODEL ---------------------------------
c This model is described in:
c
c ONETWO INPUT SWITCES FOR THE glf23 MODEL
c
c ------------------------------------------------------------------------
c% include_gf = 1 must be set (in BOTH analysis and simulation modes)
c in order to activate this model.
c (default include_gf = 0)
c the glf energy diffusivities will be used
c in simulation mode for any dependent variable for which the
c simulation mode flag is set (i.e., itte, itti) if
c
c const values, not be changed by namelist input,are set here:
c
jshoot_gf = 0 ! shooting method not available in ONETWO .. don't change.
i_grad_gf = 0 ! default 0, D-V method i_grad=1 input gradients
c
c Namelist defaults for glf23 model:
c% jmm_gf = 0 ! grid number; jmm=0 does full grid jm=1 to jmaxm-1
c > angrotp_exp, ! 0:jmaxm exp plasma toroidal angular velocity 1/sec
c > ! if itport_pt(4)=0 itport_pt(5)=0
c > egamma_exp, ! 0:jmaxm exp exb shear rate in units of csda_exp
c ! if itport_pt(4)=-1 itport_pt(5)=0
c > gamma_p_exp, ! 0:jmaxm exp par.vel. shear rate in units of csda_exp
c > ! if itport_pt(4)=-1 itport_pt(5)=0
c > vphi_m, ! 0:jmaxm toroidal velocity m/sec
c > ! if itport_pt(4)=1 itport_pt(5)=0 otherwise output
c > vpar_m, ! 0:jmaxm parallel velocity m/sec
c > ! if itport_pt(4)=1 itport_pt(5)=1 otherwise output
c > vper_m, ! 0:jmaxm perp. velocity m/sec
c > ! if itport_pt(4)=1 itport_pt(5)=1 otherwise output
c > zeff_exp, ! 0:jmaxm ne in 10**19 1/m**3
c > bt_exp, ! vaccuum axis toroidal field in tesla
c > rho, ! 0:jmaxm 0 < rho < 1 toroidal flux surface label
c > arho_exp, ! unit length toroidal flux surface LCFS in meters
c > ! toroidal flux= B0*rho_phys**2/2,
c > ! B0=bt_exp, arho_exp=rho_phys_LCFS
c > gradrho_exp, ! 0:jmaxm dimensionless <|grad rho_phys |**2>
c > gradrhosq_exp, ! 0:jmaxm dimensionless <|grad rho_phys |>
c > !NOTE:can set arho_exp=1.,if gradrho_exp=<|grad rho |>
c > ! and gradrhosq_exp = <|grad rho |**2>
c > rmin_exp, ! 0:jmaxm minor radius in meters
c > rmaj_exp, ! 0:jmaxm major radius in meters
c > rmajor_exp, ! axis major radius
c > q_exp, ! 0:jmaxm saftey factor
c > shat_exp, ! 0:jmaxm rho d q_exp/ d rho
c > alpha_exp, ! 0:jmaxm MHD alpha from experiment
c > elong_exp, ! 0:jmaxm elongation
c > amassgas_exp, ! atomic number working hydrogen gas
c > alpha_e, ! 1 full (0 no) no ExB shear stab
c > x_alpha, ! 1 full (0 no) alpha stabilization with alpha_exp
c > !-1 full (0 no) self consistent alpha_m stab.
c ------------------------------------------------- HSJ/12/4/98 --------
c
set_chie_chii = 0.0
irgc = 0
rgca = 2.0
rgcb = 1.0
rgcc = 1.0
rgcd = 2.0
times_rgc = -1.0e35 ! checked in subroutine CDEBOO
timee_rgc = -1.0e35
do j=1,kj
rgc(j) = -1.0e35 ! inititalize to non-physical number
end do
c
do j=1,kddebug ! zero debug switches
ddebug(j) = 0.0
end do
ddebug(3) = 1.0 ! means all of new, none of old value
do i=1,4
do j=1,4
wneo(i,j) = 1.0
end do
end do
do i=1,4
wneo(i,5) = 0.0
wneo(5,i) = 0.0
end do
wneo(5,5) = 1.0
wneot = 1.0
w12typ = 0.0
w13typ = 0.0
wneo_elct_ifs = 0.0
do i=1,2
vtyp (i) = 0.0
w1typ(i) = 0.0
w2typ(i) = 0.0
w3typ(i) = 0.0
w4typ(i) = 0.0
end do
wrebut = 0.0 ! Rebut-Lallia-Watkins model off
relaxrebut = 1.0 ! relaxation of Rebut-Lallia diffusivity off
qrebsmth = 0 ! no smoothing of q in RLW model
tirlw = 0 ! use real ti in RLW model
rlw_model = 'new' ! use corrections for ion rho-star scaling
resistive = 'hinton' ! in tfact.i
ftcalc = 'analytic' ! in tfact.i
wweiland = 0.0
include_weiland = 0 ! exclude Weiland model calculations
include_ifs = 0 ! exclude IFS model
include_gf = 0 ! exclude glf23 model
dorl_kotch = 0.0
dorl_kotche = 0.0
dorl_kotchi = 0.0
c
c --- modified gyro Bohm defaults
c
ce0_mgb = 1.0
alpha_mgb = 1.0
ci1_mgb = 1.0
ci2_mgb = 1.0
ce_bgb = 1.0
cfe_mgb = 0.0 ! cfe_mgb or cfe_bgb must be non zero if include_itb=1
cfe_bgb = 0.0
cfi_mgb = 0.0 ! cfi_mgb or cfi_bgb must be non zero if include_itb=1
cfi_bgb = 0.0
c1_g = 100.0 ! units are kHz
c2_g = 0.0
c_theta = 0.0
include_itb = 0 ! do not include model by default
c
c --- Hsieh (aka "shay") model defaults
c
wshay = 0.0 ! turn off Hsieh model
scsmult = .false. ! don't calculate multiplier for Hsieh model
suserho = .false. ! use simple (rather than flux avg) form
smult = 0.0 ! constant in Hsieh model
snexp = 2.0 ! density exponent
sbpexp = 2.0 ! bp exponent
srexp = 3.0 ! rho exponent
sbigrexp = 1.0 ! R exponent
stexp = 1.0 ! te exponent
sdtdrexp = 2.0 ! dte/dr exponent
srin = -1.0 ! inside radius (normalized rho)
srout = -1.0 ! outside radius (normalized rho)
skimult = 1.0 ! ki = 1.0 * ke
sdenscale = 1.0e13 ! Hsieh model elec dens in units of sdenscale
c
c --- New 1995 upgrade of the Hsieh model
c
ishayform = 1 ! Use the dimensionally correct form of the model
skimass = 2.0 ! Effective ion mass (amu) for the ion model
c
c ----------------------------------------------------------------------
c
do i=1,16
typa(i) = 0.0
end do
nw1pro = 0
nvpro = 0
nw2pro = 0
nw3pro = 0
nw4pro = 0
do j=1,ksplin
do m=1,kbctim
w1pro(j,m) = 0.0
vpro(j,m) = 0.0
w2pro(j,m) = 0.0
w3pro(j,m) = 0.0
w4pro(j,m) = 0.0
end do
end do
relaxtyp = 1.0
wstd = 0.0
betlim = 100.0
betlim0 = 20.0
jgboot = 0
jhirsh = 0
cer_ion = ' '
squeeze = .false.
iwangrot = 0
c
c --- multipliers for Mattor-Diamond diffusivities (etai modes)
c
fimpurty = 1.0
wetaie = 1.0
etaioff = 1.0
do j=1,kion
wetai(j) = 1.0
end do
*
c ----------------------------------------------------------------------
c --- PARAMETER FOR CLASSICAL ION THERMAL CONDUCTIVITY
c ----------------------------------------------------------------------
c% w3cla Weight for classical ion thermal conductivity
c ----------------------------------------------------------------------
c
w3cla = 0.0
c
c ----------------------------------------------------------------------
c PARAMETERS FOR TRANSPORT MODELS ASSOCIATED WITH MHD ACTIVITY
c ----------------------------------------------------------------------
c% w1isl, Weights for island-induced transport
c% w2isl,
c% w3isl
c
c% w1saw, Weights for sawtooth-induced transport
c% w2saw,
c% w3saw (simple model)
c the simple model modifies
c a) particle diffusion coeffient by w1saw
c b) electron diffusivity by w2saw
c c) ion diffusivity by w3saw
c d) resistivity
c if the profiles are given (analysis mode) then
c the only effect is the change in resistivity which
c modfies the current diffusion (provided itxj = 1)
c
c% w1mix, Weights for sawtooth-induced transport
c% w2mix, (mixing model)
c% w3mix,
c% w4mix,
c% w5mix
c% mixpro 2, Generate mixed profiles consistent with helical
c flux conservation as suggested by Kadomtsev for
c a single-valued q profile and since extended by
c Parail and Pereverzev to the case of a double-
c valued q profile. Mix temperatures instead of
c energy densities; then renormalize to conserve
c energy.
c 1, Generate mixed profiles consistent with helical
c flux conservation as above, but mix energy densi-
c ties. This can lead to unphysical temperature
c profiles if the density profiles are peaked and
c not mixed.
c 0, Generate mixed profiles for densities and temper-
c atures that are flat.
c -1, Generate mixed profiles for densities and energies
c that are flat.
c% dtemix Change in electron temperature (keV) on axis during
c sawtooth oscillation (used when w2mix<0)
c% dtimix Change in ion temperature (keV) on axis during
c sawtooth oscillation (used when w3mix<0)
c% fusmix Relative change (dNdot/Ndot) in neutron rate during
c sawtooth oscillation (used when w3mix<0 and
c dtimix = 0)
c% qmix Critical value of the safety factor for
c triggering a sawtooth disruption
c% rsmixx Critical value of the normalized q = 1 radius (rs/a)
c for triggering a sawtooth disruption
c% s3mix Relative change (dS/S) in signal from SXR radial
c diode 3 in old Doublet III diode array during sawtooth
c oscillation (used when w2mix<0, dtemix = 0, and jsxr=1)
c% s71mix Relative change (dS/S) in signal from SXR vertical
c diode 71 in old Doublet III diode array during sawtooth
c oscillation (used when w2mix<0, dtemix = 0, s3mix=0, and
c jsxr = 1)
c% s18mix Relative change (dS/S) in signal from SXR side diode
c 18 in new Doublet III diode array during sawtooth
c oscillation (used when w2mix<0, dtemix = 0, and jsxr=2)
c% trmix Ratio of ion temperature change to electron tempera-
c ture change on axis during sawtooth oscillation
c (used when w3mix<0, dtimix = 0, and fusmix=0)
c% tsmix Critical value of the sawtooth period (s)
c for triggering a single sawtooth disruption
c% tdmix(1),(2) Critical values of the sawtooth period (s)
c for triggering double sawtooth disruptions
c% ipmix Flag for triggering a sawtooth disruption using
c the prescription of Parail and Pereverzev
c% w0mix Initial value for the width (cm) of the m = 1 island
c following a sawtooth disruption
c sawtooth disruption
c ----------------------------------------------------------------------
c
w1isl = 0.0
w2isl = 0.0
w3isl = 0.0
w1saw = 0.0
w2saw = 0.0
w3saw = 0.0
w1mix = 0.0
w2mix = 0.0
w3mix = 0.0
w4mix = 0.0
w5mix = 0.0
dtemix = 0.0
dtimix = 0.0
fusmix = 0.0
mixpro = 2
qmix = 0.0
rsmixx = 0.0
s3mix = 0.0
s71mix = 0.0
s18mix = 0.0
trmix = 0.0
tsmix = 0.0
tdmix(1) = 0.0
tdmix(2) = 0.0
ipmix = 0
w0mix = 0.0
*
c ----------------------------------------------------------------------
c --- FLAG FOR IDEAL BALLOONING MODE STABILITY CONSTRAINT
c ----------------------------------------------------------------------
c% ibaloo 1, calculate the maximum plasma pressure gradient allowed
c for stability to ideal ballooning modes, and increase
c energy transport, if necessary, to prevent this pressure
c gradient from being exceeded
c 0, do not perform any calculations related to ideal
c ballooning modes
c -1, calculate the maximum plasma pressure gradient allowed
c for stability to ideal ballooning modes, and print this
c out, but do not increase energy transport
c ----------------------------------------------------------------------
c
ibaloo = 0
*
c ----------------------------------------------------------------------
c --- MESH PARAMETERS
c ----------------------------------------------------------------------
c% imesh 0, uniform mesh
c USE ONLY UNIFORM MESH IN TIME-DEPENDENT EQDSK MODE, HSJ
c 1, nonuniform mesh with r(j) input
c 2, nonuniform mesh with dr(j) input
c 3, uniform mesh in r(j)**2
c% nj Number of mesh points (3 to kj)
c% r(j) Position (cm) of jth mesh point
c% dr(j) Width (cm) of jth mesh interval
c% rho_edge A normalized value of rho that defines the plasma boundary
c FOR PURPOSES OF SOLVING THE DIFFUSION EQUATIONS (ONLY)!
c This option is valid only in combination with the tdem mode.
c (The required input for tdem mode is given in NAMELIST 3).
c Normally rho_edge is left at its default value of 1.0
c which means the equations are solved on a grid which is
c determined by the last closed flux surface in the eqdsk.
c You may elect (for the IFS model for example) to solve
c over a smaller interval, expressed in normalized rho space,
c rho=[0., rho_edge]
c instead by specifying 0< rho_edge <= 1.
c The boundary codition at rho_edge for each of the profiles
c run in simulation mode is obtained by interpolation from
c the profiles given in the inone file. Note that this
c means that the input profiles given in inone MUST!!!!!
c be defined over rho=[0.,1.] even if rho_edge < 1.0 is input.
c
c IT IS ASSUMED THAT PROFILES OF THE QUANTITIES TO
c BE TRANSPORTED ARE INPUT AS A FUNCTION OF SPACE AND TIME
c JUST LIKE THEY WOULD BE IN ANALYSIS MODE.!!!!!!
c
c These profiles will
c be used only to generate the required boundary condition.
c (For Faraday's law we get the rbp profile as a function of
c space and time from the tdem mode)
c if you set rho_edge < 0.0 or rho_edge > 1.0 the code will
c ignore your request and use the default of rho_edge=1 !
c% fix_edge_te
c% fix_edge_ti
c set these values to a grid point less than nj to fix the te,ti
c tmeperatures at input values for j > fix_edge_te,ti.
c the values will be interpolated from the spline input of te,ti
c using this option disaples the rho_edge option !!!!!!!!!
c% eps_adaptive(i) ,i=1...nion+4
c selects the adaptive grid calculations.
c These calculations dynamically adjust the rho
c grid spacing as the solution evolves,placing more grid points
c near regions of large variation of the profiles selected by
c eps_adaptive(i).
c the index i corresponds to the
c profiles and hence depends on how the input is set up
c i=1 to nprim corresponds to primary ion densities
c i=nprim+1 to nion (where nion = nprim+nimp)
c corresponds to the impurity densities
c i=nion+1 corresponds to Te
c i=nion+2 corresponds to Ti
c i=nion+3 corresponds to RBP
c i=nion+4 corresponds to Toroidal rotation
c eps_adaptive(i) specifies the weight of each GRADIENT
c to be included in determining the adaptive grid.
c Note that you can include analysis mode profiles in this
c scheme.(I dont know why you would want to however).
c AN EXAMPLE:
c suppose we have an inone file with nprim=2,nimp=1
c eps_adaptive (1) =0.0 ! ignore 1'st primary ion density
c ! gradient in grid determination
c (2) =1.0 ! include 2'nd primary ion density
c ! gradient in grid determination
c (3) =0.0 ! ignore impurity ion density
c ! gradient in grid determination
c (4) =2.0 ! include te gradient with twice the
c ! density gradient importance
c (5) =1.0 ! include ti gradient
c (6) =0.0 ! exclude rbp gradient
c (7) =0.5 ! include rotation gradient at half
c ! the importance of density gradient
c ! this value would be included only
c ! if iangrot=1 in the input file
c Typically only the te gradient wpuld be selected with all
c other values set to 0.0 But the above generality is coded in
c if eps_adaptive(i)=0 for all i (default) then the
c adaptive grid calculations are not used.
c in addition to (or in place of) the gradient weighting
c given by eps_adaptive above there is identical weighting
c for curvature given by curve_eps_adaptive (1 to 7)
c
c% speed_adaptive a number between 0.0 and 1.0 that controls the
c rate at which the r grid moves. numbers near zero
c cause slow evolvement of the grid. speed_adaptive=1.0
c means switch to the newly calculated grid immediately.
c This generally is to be avoided since the term
c drho/dt at constant zeta would then introduce a large
c time dependent source perturbation.
c Default is 0.01. (we want the grid to evolve but
c generally it should evolve slowly enough that the
c code doesnt start changing time steps,etc due to
c sources caused by the changing grid).
c% gridgenvb set to 1 to get some output to the screen
c regarding the grid generation. (primarily for
c developers)
c% freeze_adaptive adaptive grid generation stops at time
c t .gt. freeze_adaptive (within a time step,no effort
c is made to stop exactly at the specified time)
c ----------------------------------------------------------------------
c
imesh = 0
nj = kj
do i=1,kk
eps_adaptive(i) = 0.0
curve_eps_adaptive(i) = 0.0
end do
speed_adaptive = 0.05
freeze_adaptive = -1.0e30 ! disable
gridgenvb = 0
rho_edge = 1.0
fix_edge_te = kj + 1
fix_edge_ti = kj + 1
*
c ----------------------------------------------------------------------
c --- INITIAL CONDITIONS (PROFILES) AND BOUNDARY CONDITIONS
c ----------------------------------------------------------------------
c The standard profiles for the particle densities, temperatures,
c current density, and Zeff have the following form:
c u(k,j) = ub(m,k) + (uc(m,k)-ub(m,k))*x(j)**alpu(m,k) ,
c with
c x(j) = 1.0 - (r(j)/r(nj))**gamu(m,k) .
c
c The central values, boundary values, and exponents at time point
c m (to be described shortly) are as follows:
c
c% enec(m) , eneb(m) , alpene(m) , gamene(m)
c% enpc(m,i) , enpb(m,i) , alpenp(m,i) , gamenp(m,i)
c% enic(m,i) , enib(m,i) , alpeni(m,i) , gameni(m,i)
c% tec(m) , teb(m) , alpte(m) , gamte(m)
c% tic(m) , tib(m) , alpti(m) , gamti(m)
c% xjc(m) , xjb(m) , alpxj(m) , gamxj(m)
c% zeffc(m) , zeffb(m) , alpzef(m) , gamzef(m)
c
c These produce profiles for the following variables:
c
c% ene(j) Density of electrons (1/cm**3)
c% enp(j,i) Density of ith primary species (1/cm**3)
c% eni(j,i) Density of ith impurity species (1/cm**3)
c% te(j) Electron temperature (keV)
c% ti(j) Ion temperature (keV)
c% curden(j) Current density (see below)
c% zeff(j) Effective charge number of ions
c
c The units for xjb and xjc are arbitrary, but curden is in A/cm**2,
c which is obtained by normalization to
c
c% totcur(m) Total plasma current (A)
c
c In 1-1/2-D runs if irguess < 0 and totcur(1) = 0, then
c totcur(1) is read from the eqdsk file.
c if totcur(1) = -1.0e30 then the current is obtained from the third
c namelist of inone,from vector pcurmhd,and if appropriate,pcurmhd
c is interpolated onto the bctime array to generate time-dependent
c totcur vector. bctime and timeqbcd must be commensurate for the
c interpolation to work. If not an error message will be generated.
c totcur has precedence over pcurmhd so pcurmhd will not be used if
c totcur(1) is set to something other than -1.0e30
c
c The particle densities may be specified in different ways, depending
c upon the value of the following parameter:
c
c% inenez 1, input electron density and zeff, and calculate the
c hydrogen and impurity densities. up to two
c hydrogen species and a single impurity species are
c allowed. for two hydrogen species the fraction
c of the first is specified by zfrac. (icenez = 0)
c Note: In simulation mode the electron density is held
c fixed in time. Hence as the primary ions diffuse
c the impurity density changes to keep zeff from
c changing.
c NOTE: if ifus = -1 ,up to three primary ions may be input
c see ifus=-1 for details
c 0, calculate electron density and zeff from primary and
c impurity ion densities. (icenez = 1)
c -1, input electron density and zeff, and calculate the
c hydrogen and impurity ion densities for initial time
c only. after initial time the electron density and
c zeff are calculated from hydrogen and impurity
c densities. (icenez = -1)
c% adjzeff if inenez = 1 or -1 zeff can lead to negative
c ion densitites. by setting adjzeff (which is an
c integer variable) = 1,the code will decrease zeff
c locally until a value for zeff is found that will
c lead to positive primary ion and impurity densitites.
c If zeff attains its minimum value of 1.0 no further
c adjustment is made.
c adjzeff = 0 (default ) does not make this adjustment.
c DO NOT USE ADJZEFF IN COMBINATION WITH TWEAKING OF ZEFF!
c
c Instead of the standard profiles, spline profiles may be specified.
c These profiles are controlled by the following parameters:
c
c% rnormin(j) normalized radii (from 0 to 1) at which input data
c are specified; these locations are called knots
c% enp(j,i) density of primary ion species i at knot j
c% njenp(i) >0, number of knots in spline profile
c =0, spline profile is not generated.
c% eni(j,i),njeni(i) - same for impurity ion species i.
c% te(j),njte - same for electron temperature.
c% ti(j),njti - same for ion temperature.
c% ene(j),njene - same for electron density.
c% zeff(j),njzef - same for zeff.
c%