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% curden(j),njcur - same for current density.
c Spline profiles specified by the preceding controls are
c time-independent. Time-dependent spline profiles may be
c specified for the following quantities. (The spline knot
c values are interpolated linearly in time.):
c
c% tein(j,m) - electron temperature at knot j at time point m.
c% njte - same as defined above.
c% tiin(j,m),njti - same for ion temperature.
c% enein(j,m),njene - same for electron density.
c% zeffin(j,m),njzef- same for zeff
c CAUTION: all non-zero njxx MUST have the same value if splninpt = 'old'
c if splninpt = 'new' this restriction is removed,and curden is also
c allowed to be time-dependent. splninpt = 'old' is retained to allow full
c backwards compatibility. However it is recommended that splninpt = 'new'
c is used for time-dependent spline profiles from here on out.
c default is splninpt = 'old' to be consistent with old input files
c **********************************************************************
c *** SPLNINPT = OLD IS NO LONGER ACCEPTED *****************************
c **********************************************************************
c specifically the following input can be set using different knot sets:
c% enp(j,i) renpin(j,i) njenp(i) bparenp(4,i) (no time dep.)
c% eni(j,i) reniin(j,i) njeni(i) bpareni(4,i) (no time dep.)
c% curdenin(j,m) rcurdein(j,m) njcur bparcur(4,m)
c% tein(j,m) rtein(j,m) njte bparte(4,m)
c% tiin(j,m) rtiin(j,m) njti bparti(4,m)
c% enein(j,m) renein(j,m) njene bparene(4,m)
c% zeffin(j,m) rzeffin(j,m) njzef bparzeff(4,m)
c% angrotin(j,m) rangrot(j,m) bparang(4,m) (see below)
c% kpolin(j,m) rkpol(j,m) bparkpol(4,m) (set ikpol=1)
c tein, tiin, enein, zeffin and curdenin may be input as
c te , ti , ene , zeff and curden respectively,
c if there is no time dependence.
c here index j is the knot and i is the species as usual
c index m is the time index (see bctime).
c rxx specifies the knots (0.0 to 1.0) as rnormin but with
c possibly different knots for each input profile.
c If njxx .ne. 0 but the corresponding rjxx is not set,
c then rnormin will be used for that profile.
c the njxx do not have a time index m because they are not used to
c determine the number of knots in rxx(j,m) at time point m. Rather,since
c the first knot must be at zero and the last knot must be at 1.0,the code
c discovers for itself what the number of knots is at any time m > 1.
c however,in order to let the code know that spline (rather than parabolic)
c input profiles are given,each of the njcur,njte,njti,njene,and njzef must
c be set to the initial number of knots at time m = 1.0
c for the rotation speed profile parabolic input is not allowed so no njxx
c switch is present in the input.
c The bparxx arrays are used to set boundary conditions.
c (on the spline input,-has nothing to do with boundary conditions on
c transport equations)
c normally the b.c. at rho = 0 is zero gradient so bparxx(1) and bparxx(2)
c do not have to be set since zero gradient at rho=0 is the default.
c however to set the gradient at rho = 0 to the value x,specify
c bparxx(1) = -1.0e30,bparxx(2) = x
c at rho = 1.0 the default boundary condition is the natural spline
c to change it to a gradient condition,say gradient = y, proceed as above:
c bparxx(3) = -1.0e30,bparxx(4) = y
c NOTE: at present only a gradient b.c. is allowed as input.
c when using the splninpt='new' option,you may,but do not have to
c specify profiles for all times nbctim. Any profile and knot set
c which is required but is not given in the input will be copied
c from the last previous time point at which it is available. Hence
c profiles which do not change (and the corresponding knot set ) need
c only be given onece. extensive error checking is done
c and a summary print of the input spline profile parameters is given
c in the outone file.
c because there is no easy way to default or check the bparxx arrays,
c these arrays must be given explicitly at all times m! The only
c exception is if bparxx(j,m) is zero for all m and j = 1...4 In this
c case the bparxx array need not be specified at all.
c
c% angrotin(j,m) angular rotation speed,rad/sec,at knot rangrot(j,m)
c at time bctime(m). if angrotin is time-independent but
c nbctim .gt. 1 then the values of angrotin for other times
c do not have to be supplied(see rangrot for further info)
c% rangrot(j,m) similar to rnormin but for angular rotation speed
c profile. rangrot(j,m) gives the j'th knot
c at time bctime(m). Since the first knot must
c always be at zero and the last knot must always
c be at one the code automatically determines the number
c of knots input so no variable specifying the number
c of knots at time bctime(m) is required.
c the range on j is [3,ksplin]
c if the angular rotation speed is to be time-independent
c then simply specify angrotin(j,m) with m = 1. for m = 2..
c nbctim leave angrotin(j,m) unspecified (i.e., 0.0).
c the knot locations and number of knots at each
c time point m are not required to be the same. if they
c are the same however the user need only give the values
c for m = 1 in rangrot(j,m). More generally if angrotin
c and/or rangrot is the same for all times greater than
c the last specified time the values do not have to be
c repeated.
c there is no plan at present to allow input of rotation
c speed profiles in terms of the parabolic parameter set.
c% bparang(4,kbctim) boundary condition array for angrotin.
c note that angrotin, bparang and rangrot are independent of the setting
c of splninpt (i.e., works with either splninpt = 'old' or 'new')
c (there is no backward compatibility problem here).
c
c% kpolin(j,m) poloidal velocity of cer ion/Bp at knots. This is only
c used to compute the radial electric field from the CER
c data if cer_ion is set. Output of Epsi_exp and Kpol_exp
c is to the netCDF file iterdb.nc (must set iterdb = -1)
c
c% rkpol(j,m) knot location for kpolin similar to rangrot above
c
c% bparkpol(4,kbctim) boundary conditions for splines similar to bparang above
c
c Time-dependent profiles and boundary conditions are controlled by
c two parameters:
c
c% nbctim Number of times at which boundary conditions
c are input (up to kbctim). See note below.
c% bctime(m) Times (sec) at which boundary conditions are input
c
c All preceding quantities with subscript m,
c plus elong(m), gasflx(m,i), w*typ(m), w*pro(j,m) and qradin(j,m),
c are given at the time values in bctime.
c If nbctim = 1, all quantities are time-independent.
c If nbctim>1, the second value of each quantity specifies whether
c it is time-dependent. If the second value is 0 (default), the
c quantity is time-independent, and the first value is used
c at all times. If the second value is not zero, it is assumed
c that nbctim values have been given for the quantity.
c (Use a very small number if the second time-dependent quantity is 0)
c
adjzeff = 0
pfact = 1.0
c
do i=1,kbctim
if (i .gt. 1) pfact = 0.0 ! else logic in SPECIFY won't work
do k=1,kprim
enpc(i,k) = pfact * 5.0e13
enpb(i,k) = pfact * 1.0e13
alpenp(i,k) = pfact * 1.0
gamenp(i,k) = pfact * 2.0
end do
do k=1,kimp
enic(i,k) = pfact * 5.0e10
enib(i,k) = 0.0
alpeni(i,k) = pfact * 1.0
gameni(i,k) = pfact * 2.0
end do
tec(i) = pfact * 1.0
teb(i) = pfact * 0.2
alpte(i) = pfact * 1.0
gamte(i) = pfact * 2.0
tic(i) = pfact * 0.4
tib(i) = pfact * 0.2
alpti(i) = pfact * 1.0
gamti(i) = pfact * 2.0
xjc(i) = pfact * 1.0
xjb(i) = pfact * 0.1
alpxj(i) = pfact * 1.0
gamxj(i) = pfact * 2.0
enec(i) = pfact * 5.0e13
eneb(i) = pfact * 1.0e13
alpene(i) = pfact * 1.0
gamene(i) = pfact * 2.0
zeffc(i) = pfact * 1.0
zeffb(i) = pfact * 1.0
alpzef(i) = pfact * 1.0
gamzef(i) = pfact * 2.0
totcur(i) = pfact * 7.0e5
end do
c
nbctim = 1
inenez = 1
do 1830 i=1,kprim
1830 njenp(i) = 0
do 1840 i=1,kimp
1840 njeni(i) = 0
njte = 0
njti = 0
njene = 0
njzef = 0
njcur = 0
splninpt = 'old' ! don't set to new (will defeat error checking)
call zeroa (tein,ksplin*kbctim)
call zeroa (tiin,ksplin*kbctim)
call zeroa (enein,ksplin*kbctim)
call zeroa (zeffin,ksplin*kbctim)
call zeroa (angrotin,ksplin*kbctim)
call zeroa (kpolin,ksplin*kbctim)
call zeroa (bparkpol,4*kbctim)
call zeroa (bparang,4*kbctim)
call zeroa (bparenp,4*kprim)
call zeroa (bpareni,4*kimp)
call zeroa (bparte,4*kbctim)
call zeroa (bparti,4*kbctim)
call zeroa (bparcur,4*kbctim)
call zeroa (bparzeff,4*kbctim)
call zeroa (bparene,4*kbctim)
*
c ----------------------------------------------------------------------
c --- SOLUTION CONTROL PARAMETERS
c ----------------------------------------------------------------------
c% time0 Initial time (sec)
c% timmax Maximum time (sec) set timmax .le. time0 for snapshot
c mode. note that in snapshot mode all dependent variables
c must be run in analysis mode and it is implicitly
c assumed that all time derivatives are zero. At present
c specifying profiles as a function of time will not yield
c the desired time derivatives in the snapshot mode.
c also note that energy source terms due to the rate of change
c of beam or fusion ion densities must be neglected since
c these rates are not available at the initial time.
c Due to roundoff it is best if timmax is set slightly less
c (say one millisec ) than bctime(nbctim) if multiple
c boundary values are specified on input.
c% nmax Maximum number of time steps
c% dt Initial time step (sec)
c% dtmin Minimum time step (sec) (suggest 1.0e-6)
c% dtevmin Minimum time step for plotting, printing, beam on (sec)
c% dtmax Maximum time step (sec) (suggest timmax/200)
c% relmin Minimum delmax, where delmax is the maximum relative
c change from one time point to the next in any dependent
c variable. If delmax < relmin, the time step is doubled.
c (suggest 0.10)
c% relmax Maximum delmax. If delmax > relmax, the time step is
c halved. (suggest 0.30)
c Note: both simulation and analysis mode profiles are
c checked. As a consequence the time step will be cut
c if the change in a profile is greater than relmax,even if
c that profile is not run in simulation mode. In particular if
c you run in total analysis mode (i.e., all itran(i) =0) then
c no profiles are transported but the change in
c input profiles and beams may still cut doewn the time step.
c The neutral density equation has to be solved even in this
c case and may also result in a decreased time step.
c Set relmax to a large number to avoid this as appropriate.
c% relit Maximum delit, where delit is the maximum relative
c difference between the predictor and corrector in any
c dependent variable. If delit > relit, the corrector
c calculation is iterated. (suggest 0.01)
c set relit = a large number and dtmax a small time step
c to mimick Euler type solution method HSJ
c% itmax Maximum number of corrector iterations (suggest 10)
c (Note: includes particle conservation iterations.)
c% theta Time weighting of diffusion term in difference
c equations:
c 0.0, fully explicit (unstable)
c 0.5, Crank-Nicolson (marginally stable)
c 1.0, fully implicit (stable)
c (suggest 0.8)
c% ilimdt 1, allow restriction of dtmax to numerically stable
c regime, as determined by the routine DTLIM (default).
c 0, disable this: dtmax is not touched.
c The standard Euler approach can be mimicked by setting
c ilimdt=1, relit = a large number,dtmax a small time step
c (small means relative to the rate at which the sources are
c forcing the solution to change,typically .1 to 10 ms) HSJ
c% timav >0, time average solution over the interval timav;
c output time-average results at last time point
c 0, no effect (default)
c <0, time average solution over the interval timav;
c output last time results in addition to time-average
c results
c% iffar 1, Fast Faraday's law: speeds up the current diffusion in
c Faraday's law. Useful when steady state is desired.
c This option is invoked automatically when 'tweaking'
c with voltin or qcin non-zero.
c 0 no effect (default)
c
c% relaxsol relaxation parameter for the solution during
c the corrector iterations. set relaxsol = 1.0 to
c disable the relaxation. otherwise relaxsol should
c be set to a number between 0 and 1. to underelax
c the solution,if it oscillates due to a nonlinearity
c in the transport coefficients. relaxsol only
c needs to be used for models such as Rebut-Lallia
c when running te in simulation mode to prevent
c extremely short time steps.
c ----------------------------------------------------------------------
c
relaxsol = 1.0
time0 = 0.0
timmax = 200.0e-3
timecrit = 1.0e+39
nmax = 200
dt = 1.0e-3
dtmin = 1.0e-6
dtevmin = dtmin
dtmax = 100.0e-3
dtmaxcrit = dtmax
relmin = 0.10
relmax = 0.30
relit = 0.01
theta = 0.8
itmax = 10
ilimdt = 1
timav = 0.0
iffar = 0
idterr = 0
ineucg = 0
ifreya = 0
ifreya_old = 0
do k=1,krf
irfcnt(k) = 0
end do
diffeq_methd = 0
*
c ----------------------------------------------------------------------
c --- OUTPUT PARAMETERS
c ----------------------------------------------------------------------
c% timprt Time interval (sec) between printout
c% timplt Time interval (sec) between 3-D plots
c% mprt Number of time steps between printout
c% mplot Number of time steps between 3-D plots
c% jterow List of grid points for te plot -
c plot te(jterow(i)) versus time
c for i = 1 to nterow
c% nterow See above comment
c% prtlst List of times (sec) for printout (max = 10)
c NOTE: the 10 is hardwired in the code@##! if you change it
c you need to be verrry careful
c prtlst also controls times at which the "analysis mode" variable
c parameters are plotted
c These are the sequence of plots thermal conductivites,
c temperatures,densities,confinement times ...etc that follow the
c 3D plots in the trplot.cgm file output. THE ACTUAL TIME
C PLOTTED IS THE CORRECTOR PIVOT POINT TIME
c GIVEN BY PRTLST(I)-(1-THETA)*(CURRENT TIME STEP)
c (see desription of theta ) HENCE THE TIMES IN THE PLOTS WON'T
c MATCH THE PRTLST TIMES EXACTLY. (LIMIT THE MAXIMUM CURRENT TIME
c STEP WITH DTMAX IF THIS IS A PROBLEM FOR YOU)
c
c% pltlst List of times (sec) for 3-D plots (max = 30)
c Also see above prtlst list!
c ONE and only ONE 3d plot will be made for each
c of the parameters which is designated as a 3D quantity.
c pltlst guarantees that the profile as a function of rho at the
c times pltlst(i) will be part of this single 3d plot, however
c the dimensions of the arrays may conspire to defeat this command)
c Since info for the 3d plots is generated at the end of each
c time step this option really only makes sense if there are
c too few time steps in your analysis. However the number of time
c steps is best controlled using dtmax so, in my (HSJ) opinion,
c this switch is rather useless. We may remove it in the future.
c it is here only for backward compatibility. -- HSJ)
c
c% jprt Number of mesh intervals between printout
c% jflux 0, do not print out components of fluxes
c 1, do print out components of fluxes
c% jcoef 0, do not print out transport coefficients
c 1, do print out transport coefficients
c% jsourc 0, do not print out comps. of particle and energy sources
c 1, do print out comps. of particle and energy sources
c% jbal 0, do not print out balance tables
c 1, do print out balance tables
c% jtfus 0, do not print out detailed fusion information
c 1, do print out detailed fusion information
c if ifus=1,jtfus=1 is forced.
c JSXR IS OUT OF DATE (OBSOLETE, USES DOUBLET III STUFF
c don't use this option blindly) HSJ
c% jsxr 0, do not calculate and print SXR info
c 1, do calculate and print SXR info for old Doublet III
c diode arrays
c 2, do calculate and print SXR info for new Doublet III
c diode arrays
c LEAVE JSXR SET AT 0 UNTIL FURTHER NOTICE. HSJ
c% jco2 0, do not calculate and print CO2 info
c 1, do calculate and print CO2 info for Doublet III
c (line-averaged electron density for various
c tangency radii)
c% jzeff 0, do not calculate and print ZEFF info
c 1, do calculate and print ZEFF info for Doublet III
c (visible continuum light emission for various
c tangency radii)
c LEAVE JZEFF SET AT 0 UNTIL FURTHER NOTICE. HSJ
c% spec_profiles(j) j=1,...up to nj
c spec_profiles(j)=k is used to print out the
c values te(k) and ti(k) at radial grid point k
c as a function of time. (note that this
c is at constant zeta, not necessarily at constant rho)
c% steps_per_plot integer switch,default value of 1,that determines the
c number of time steps between calls to the plotting
c routine trplot. By setting steps_per_plot to an
c integer greater than 1 the output plot file can be
c reduced in size if coarser plots are acceptable.
c ----------------------------------------------------------------------
c
call zero_intg (spec_profiles, kj)
timprt = 1.0e6
mprt = 0
jprt = 5
timplt = 0.050
mplot = 0
nterow = 0
jflux = 0
jcoef = 0
jsourc = 1
jbal = 0
jtfus = 0
jsxr = 0
jco2 = 0
jzeff = 0
steps_per_plot = 1
do i=1,30
pltlst(i) = 1.0e6
if (i .le. 10)
. prtlst(i) = 1.0e6
end do
*
c ----------------------------------------------------------------------
c --- NEUTRAL TRANSPORT PARAMETERS
c ----------------------------------------------------------------------
c% nneu Number of neutral species (0, 1 or 2)
c% namen Names of neutral species: For these primary species,
c neutral transport and radiative recombination of ions
c will be modeled in the code. (Note that He is treated
c as if there were no singly-ionized state.)
c Permitted for primary species #1 and #2, only.
c FOR TWO NEUTRAL SPECIES INPUT THE NEUTRALS IN THE SAME ORDER
c AS THE PRIMARY IONS HSJ
c NOTE The following parameters that control neutral profiles
c take effect only if the primary ions are run in
c simulation mode. For analysis mode use taupin to
c get some control.
c% gasflx(m,i) Flux of injected neutral gas for primary species i.
c (1/cm**2-sec). If time-dependent, specify nbctim and
c bctime(m).
c% ipcons(i) 0, recycling of species i (default) with coeff recyc (iref = 1)
c 1, adjust recycling of species i to maintain a constant
c number of ions in the plasma.
c% recyc(i) recycling coefficient of species i. (default = 1.0)
c% twall Temperature (keV) of neutrals emitted from wall
c% wion Electron energy loss (keV) per ionization due to
c ionization
c% wrad Electron energy loss (keV) per ionization due to
c radiation
c% raneut Effective plasma radius (cm) seen by neutrals;
c if input as zero, raneut is calculated internally as
c SQRT (elong(1))*rminor (elong may be time-dependent).
c% relneu Minimum relative change in ion density or 1/4 of ion
c temperature between neutral transport calculations
c (suggest 0.10)
c% erneumax Maximum epsneu, where epsneu is the relative error in
c the total number of primary ions in the plasma. If
c epsneu > erneumax, the neutral density profile is
c adjusted and the corrector calculation is iterated.
c (suggest 0.005)
c
c% widths Width of plasma scrape-off layer surrounding the
c main plasma (cm). Default = 0.0, no scrape-off layer.
c Transport is not modeled in this layer: it is used
c to obtain more realistic neutral densities within the
c main plasma.
c% nps Number of points in the scrape-off layer (suggest 4).
c These are added to the nj points in the main plasma.
c% enes Electron density in the scrape-off layer (1/cm**3).
c Default = 0.0 causes eneb to be used. The charge
c state of impurities and proportions of primary ion
c species used in the layer are the same as those at
c the edge of the main plasma.
c% tes Electron and ion temperature in the scrape-off
c layer (keV). Default = 0.0 causes teb and tib
c to be used.
c
c% idiagn Neutral diagnostic flag
c 0: do not calculate neutral diagnostic quantities
c 1: calculate specific flux of outgoing neutrals -
c rtandn,rhdn define the detector line-of-sight chord
c% rtandn Tangency chord of detector chord in major plane (cm).
c Default = 0.0, chord is radial.
c% rhdn Height above midplane of detector chord (cm).
c Default = 0.0, chord is in midplane.
c% nengn Number of energies at which spectrum is calculated
c (max = 50)
c% englstn List of energies (keV) for neutral spectrum
c% neucgvb verbose switch for neucg. set to 1 to get messages to screen
c ----------------------------------------------------------------------
c
nneu = 0
namen(1) = 'h'
namen(2) = ' '
erneumax = 0.005
do i=1,2
ipcons(i) = 0
recyc (i) = 1.0
do m=1,kbctim
gasflx(m,i) = 0.0
end do
end do
twall = 0.010
wion = 0.014
wrad = 0.006
raneut = 0.0
relneu = 0.10
widths = 0.0
nps = 4
enes = 0.0
tes = 0.0
idiagn = 0
nengn = 0
rtandn = 0.0
rhdn = 0.0
c
c ----------------------------------------------------------------------
c *** debug information ***
c ----------------------------------------------------------------------
c
do i=1,20
xdebug(i) = 0.0
end do
c
c ----------------------------------------------------------------------
c% ifred: 1 turn on Fred Marcus output; 0 turn off Fred Marcus output
c ----------------------------------------------------------------------
c
ifred = 0
*
c **********************************************************************
c *** SECOND NAMELIST (NAMELIS2) ***************************************
c **********************************************************************
c
c ----------------------------------------------------------------------
c YOKA ANALYSIS RUN PARAMETERS
c ----------------------------------------------------------------------
c% iyoka 0, Not a Yoka analysis run (default)
c 1, Yoka analysis run: compute and print out quantities
c to be stored in the Yoka Neutral Beam data bank.
c (Sets jsourc = 1 to compute needed quantities.)
c 2, Same as 1, but print an extra copy of the transport
c summary page at the end of the qikone file.
c% ishot Doublet III or DIII-D shot number for this analysis run.
c NOTE: ishot is used as part of the identifier
c for the eqdsk name for eqdsks created by
c ONETWO (this has nothing to do with input
c eqdsks which are required if irguess .ne. 0)
c Hence you may wish to assign some sort of
c identifying number to ishot even if this is
c not an official Yoka-type run.
c In tdem mode ishot should be set to the
c shot number in the netCDF file.
c% itime Time at which the shot is being analyzed (msec).
c ----------------------------------------------------------------------
c
iyoka = 0
ishot = 0
itime = 0
*
c ----------------------------------------------------------------------
c --- SOURCE MULTIPLIERS
c ----------------------------------------------------------------------
c% wdelt: multiply qdelt by wdelt
c% wgam : multiply qgam by wgam
c% wohm : multiply qohm by wohm
c% angrm2d(i) i = 1,4 multipliers for 2d angular rotation source terms
c see documentation for explanation.
c% angrcple multiplier for energy flux due to angular rotation
c in ion energy equation.
c ----------------------------------------------------------------------
c
wdelt = 1.0
wgam = 1.0
wohm = 1.0
angrm2d(1) = 1.0
angrm2d(2) = 1.0
angrm2d(3) = 1.0
angrm2d(4) = 1.0
angrcple = 1.0
*
c ----------------------------------------------------------------------
c --- RADIATION PARAMETERS
c ----------------------------------------------------------------------
c% nqrad = 0, calculate radiation power density
c > 0, number of values in an input table of radiation
c power density (watts/cm**3) as a function of
c normalized radius
c% qradr Normalized r values for table.
c% qradin Power density values (watts/cm**3).
c These may be time-dependent:
c qradin(1 - nqrad,m) are for bctime(m).
c A linear interpolation in time is used.
c% refrad Wall reflection coefficient for synchrotron radiation.
c This parameter lies between 0 and 1, with the latter
c value turning off the radiation.
c ----------------------------------------------------------------------
c
nqrad = 0
do j=1,ksplin
do m=1,kbctim
qradin(j,m) = 0.0
end do
end do
refrad = 1.0
*
c ----------------------------------------------------------------------
c --- NEUTRAL BEAM HEATING PARAMETERS
c ----------------------------------------------------------------------
c The default is for co-injection. To get ctr injection
c set angleh to a negative value (usually -19.5 degrees)
c and switch the sources (sfrac(1) becomes sfrac(2), etc.)
c Be sure to check the graphic output from program NUBPLT if you
c change the default parameters!
c
c% iterate_beam logical, if .true., allow up to imaxbeamconvg iterations
c default: iterate_beam = .false.
c% imaxbeamconvg iterations are done to get consistency
c in the thermal and fast ion densities (ONLY).
c default: imaxbeamconvg = 5
c% relaxden relaxation parameter for beam iteration,must be in (0,1]
c with 1 meaning no relaxation.
c% relaxden_err relative error for beam convergence (The relative error is
c measured in the change in electron density if icenez=0
c and as a relative change in the thermal ion species
c corresponding to the beam if icenez=1)
c% timbplt Times (up to 5) to produce data for Freya-like plots
c of beam deposition. output is processed by the NUBPLT code.
c Defaulted to OFF.
c timbplt(1) .le. time0 .and. beamon .lt. time0 gives o/p
c for initial time.
c% beamon Time (sec) at which beam heating is turned on
c% btime Time interval (sec) during which beam heating is on
c% ibcur Flag for neutral beam driven current
c 1, include beam driven current
c 0, neglect beam driven current
c% use_Callen integer,has meaning only if mcgo run is done.
c Intended for checking out some model details.
c Users should leave it set at the default value of 0.
c% ibcx Flag for secondary charge exchange between fast ions
c and thermal neutrals
c 1, include secondary charge exchange (default)
c 0, neglect secondary charge exchange
c NOTE: thermal ion particle source due to secondary charge
c exchange is not included in the ion density equation,
c independent of the setting of ibcx.
c
c% nbeamtcx switch used to determine the form of the beam torque.
c if nbeamtcx = 0 (default) then the loss of beam torque due
c to secondary charge exchange of fast ions with thermal
c neutrals is neglected. if nbeamtcx = 1 the transfer of angular
c momentum from fast ions to thermal ions and electrons is
c modified to account for the cx losses (done with array ssprcxl)
c% ibslow Flag for neutral beam slowing down
c 1, include neutral beam slowing down
c 0, neglect neutral beam slowing down
c% fionx Allows testing for non-classical slowing down
c (see subroutine SLOW1).
c
c% fast_ion_target integer,if set to 1 will let the incoming neutral beam
c see the stored fast ion density
c no attempt to modify the stopping rates due to the
c fact that the fast ion distribution is not Maxwellian
c has been made!!!!!!!!!!!!!!!!!!!!!!!
c Note that if the stored fast ion density is large enough
c that it is of concern then it is also large enough so that
c a simple one pass linearization doesnt make sense. So you
c should use this option in conjunction with the iterate
c beam option.
c For Mcgo coupled runs the birth points of the fast ions
c are determined by the Freya in Onetwo. Hence this switch
c will affect the Mcgo results as well. Note that Mcgo
c itself does not let the Monte Carlo particles see the
c stored fast ion density. As an approximate remedy for
c this situation you can set fast_ion_target = -1 when
c running a Mcgo coupled case. This will add the stored fast
c ion density determined in Freya to the thermal ion density
c used in Mcgo, so that the Mcgo particles see a correct
c total density.
c (If Mcgo is run with fast_ion_target = 1 then
c the birth points feed to Mcgo will account for the stored
c fast ions but subsequent slowing down of the Monte Carlo
c particles will see only the thermal ions)
c% rtstcx A factor between zero and 1 used to fudge secondary charge
c of fast ions. ie the original fast ion charge exchanges
c with a thermal neutral. The resulting fast neutral
c will however most likely be reionized. To model this
c reionization in Onetwo we don't do the secondary charge
c exchange calculation, which requires a Monte Carlo approach.
c Instead we assume that the probability against charge exchange
c of a fast ion, normally taken as
c {(vbeam**3+vcrit**3)/(v**3+vcrit**3)}**A
c where A=(-taus/(3.*taucx)) is the ratio of slowing down
c to charge exchange lifetime is given instead by
c {(vbeam**3+vcrit**3)/(v**3+vcrit**3)}**B
c where B=rtstcx*A. To decrease charge exchange losses
c set rtstcx less than 1. Note that this can be interpreted
c as assuming that the neutral density, enn, is replaced by
c rtstcx*enn for the purposes of charge exchange probability
c calculations ONLY. The actual neutral density is not changed!
c
c Set rtstcx to a negative number to use sigma*v, rather
c than (i.e., Maxwellian average), to determine
c the mean lifetime of a fast ion against charge exchange.
c The energy and velocity used will be that of the beam
c corrected for bulk rotation of ions (note the assumption
c that neutrals rotate as ions do).
c If rtstcx is set to a number less than -5. then
c the average fast ion energy (again corrected for rotation)
c for each beam slowing down distribution is used.
c
c% nameb Name of neutral species in beam
c 'h', 'd', 't', 'dt'
c MUST BE PRIMARY ION SPECIES #1 or #2.
c if nameb = 'h' then fdbeam defaults to natural isotopic
c content of d in 'h'.
c if nameb = 't' fdbeam is explicitly set to 0
c if nameb = 'dt' fdbeam is set to fd, (the fraction
c of d in thermal dt mixture) USER HAS NO CHOICE
c UNLESS IFUS = -1 is selected
c NOTE that nameb=dt selects a SINGLE fast ion fluid
c with a fictitious mass !!!!
c% fdbeam Fraction of deuterium atoms in neutral beam.
c% relnub Minimum relative change in ion density or electron
c temperature between neutral beam calculations
c (0.10 is suggested)
c% tfusbb 'thermal fusion balance for beams', fraction by which
c the net energy gain from thermal fusion must exceed
c the net energy loss for automatic beam turnoff.
c If tfusbb = 0, automatic beam turnoff is not done.
c% iddcal Flag controlling treatment of beam effects on fuscal,
c the calculated fusion neutron rate:
c 0 = do not include knock-on or beam-d neutrons in fuscal
c 1 = include only knock-on neutrons in fuscal
c 2 = include only beam-d neutrons in fuscal
c 3 = include both knock-on and beam-d neutrons in fuscal
c% ranseed Starting seed for random number generator used in the
c FREYA determination of the beam deposition
c% npart Number of particles followed into plasma (suggest 10000)
c% npskip Ratio of number of particles followed into plasma
c to number of source particles (suggest 1)
c% iborb Flag for modeling orbit effects on beam-generated fast ions
c 3 = use the Monte Carlo Code MCGO to model the fast
c ion slowing down. (The ionization of the injected neutrals
c is still done by FREYA as usuall) This is a very expensive
c option and shouldn't be used indiscriminately. It will take
c about 60 mins of CPU time to follow 2500 ions and
c generate the fast ion density and deposition profiles.
c If iborb=3 is used the following input is required.
c npart_mcgo the number of fast ions to be followed
c during the slowing down process. MCGO
c will follow the guiding center of each
c ion until it is thermalized or lost
c (by bad orbit or charge exchange)
c The maximum you can set here is 10000
c Typically 2000-5000 ions are sufficient
c to yield a reasonable sampling.
c DO NOT CONFUSE NPART_MCGO with the
c FREYA NPART SETTING. Npart and
c npart_mcgo do not have to be the same
c but npart MUST be >= npart_mcgo.
c MCGO will take a random sample of
c npart_mcgo initial fast ion launch
c parameters from the list of length
c npart generated by FREYA. Proper
c consideration is given to the full
c half and third energy components of the
c beam. (Npart_mcgo is the total for all
c three components)
c
c -3 iborb= -3 means that all of the mcgo preparatory
c calcs will be done and the input files required to run
c mcgo will be generated but mcgo itself will not be run
c Using this option you can run mcgo manually,
c (at a later time) rather than
c having onetwo control the execution. Note that iborb =-3
c disables most of the freya calculatiosn in onetwo so it
c should not be used for normal transport runs.
c
c In order to use this option the files
c mcgo_output_file(1,2)
c mentioned below must not be specified in inone
c (or specify them as zero-length files:
c mcgo_output_file(1) = '' )
c
c Another way to run onetwo is to use pre-existing mcgo
c results rather than have onetwo spawn mcgo during a
c transport run. To use this option set iborb =-3 and
c specify the name of the file(s) that onetwo is to read
c instead of running mcgo. One such file will be read for
c each beam so the number of files you specify has to
c equal the value of nbeams in inone. The files are
c specified in namelist 2 according to
c mcgo_output_file(1) = 'XXXXX'
c mcgo_output_file(2) = 'XXXXX'
c (substitute actual file names for the XXXXX )
c A time-dependent version of this has NOT been programed!
c So the same files will be read over and over again in a
c time-dependent ONETWO run. Obviously the propensity to
c get out of sync with the inone file is great here so use
c this option cautiously !!!!!!
c mcgo_12_output MCGO will create an output file that
c can be read by Onetwo.
c You can use
c this file on subsequent Onetwo runs
c as well, without running MCGO,
c by specifying the file name of the
c MCGO data in the variable
c mcgo_results_12
c if mcgo_results_12.nc is set to
c a valid file
c then this file will be used even if
c iborb=3 (and MCGO will not be executed)
c If mcgo_results_12 is not a valid
c file name
c then MCGO will be executed only if
c iborb=3. In this case a file with
c the name mcgo_onetwo.dat will be
c created.
c We would like to do this using PVM (parallel
c virtual machine) instead of a file
c interface eventually so the above may change.
c Other information that MCGO needs is part of the
c standard Onetwo input so no further user intervention
c is required.
c 2 = new prompt loss model of Stambaugh (implemented only for
c codeid = dee i.e., two d case)
c (REF: "Calculating Orbit Loss in a Separatrix Bounded Tokamak",
c DIII-D Physics Memo No. 9303, 16 March 1993)
c 1 => do model orbit effects (this is the original method)
c 0 => do not model orbit effects
c% itrapfi Flag for trapped fast ions. If iborb = 1 then
c setting itrapfi = 1 will modify the beam driven
c current by the initial trapped ion fraction.
c (Subsequent pitch angle diffusion is NOT taken into account).
c itrapfi = 0 neglects this effect. If iborb = 0,
c itrapfi has no effect. itrapfi = 0 is default.
c% ds_tk maximum trajectory increment (cm) used in subroutine INJECT to
c calculate psi(s), where s is the neutral trajectory
c pathlength from the first closed flux surface it
c encounters. used for non-zero toroidal rotation, where
c mean free path as a function of path length is required.
c% fe_tk factor (>1) to set upper limit of energy in n*sigma array.
c max(ebins) = max(ebkev)*fe_tk. required for nonzero
c rotation cases.
c% ne_tk number of equi-width energy bins used in forming n*sigma
c array. required for nonzero rotation cases. is internally
c reset to zero (used as flag to turn off rotational effects
c on neutral stopping) if angular rotation is not present
c (iangrot = 0).
c% iexcit Selection switch for atomic cross section data and model:
c 0, Use the fundamental atomic data of Freeman & Jones (1972).
c This option is not advised since this atomic data is
c considered outdated and excited state effects are
c not considered.
c 1, Use hexnb package but do not include excitations in
c its calculation of cross sections.
c 2, Use hexnb package, include excitations in calculation
c of cross sections.
c NOTE: Options 1 and 2 are not advised because the hexnb
c model has been found to be flawed and is based on
c atomic data that is considered outdated.
c (Boley et al, Nuclear Fusion 29, 1984)
c 5, Use the JET-ADAS effective stopping cross sections.
c This model is preferred since it is based on the most
c recent atomic data available and includes multi-step
c ionization processes due to excited states. This is an
c important effect for considering neutral beam penetration
c into fusion-grade plasmas. For a detailed description
c of the ADAS Atomic Data and Analysis Structure (ADAS)
c developed by JET, see Finkenthal, 1994.
c Note: All ions are considered fully stripped.
c (Daniel Finkenthal Ph.D. Thesis, UC-Berkeley, 1994)
c
c% neg_ion_source(i) i=1,...,nbeams
c an integer array, set to 1 to indicate
c that the a NEGATIVE ion source is used for
c neutral beam line i. At this time, the only
c effect of this switch is to set the neutralization
c efficiency, arbitrarily, to 98%, independent of the
c negative ions energy and to eliminate the second and
c third energy components from the beam.
c
c% izstrp(i) i = 1,2..nimp set to 0 for coronal equilibrium
c values of and . set to 1 for fully stripped impurities
c% mstate principal quantum number n above which excitations
c count as ionizations.
c% inubpat two-dimensional beam deposition option
c 0, do not calculate beam deposition in (r,z) coordinates
c 1, calculate beam deposition on (r,z) grid (default = eqdsk
c grid). write deposition array and n = 3 excited state
c fraction to file 'beamdep' for standalone analysis.
c% npat neutral beam deposition (r,z) grid dimensions, used if
c inubpat = 1.0
c npat(1) = number of elements in 'r' (<=2*nw)
c npat(2) = number of elements in 'z' (<=2*nh)
c defaults, npat(1) = nw, npat(2) = nh
c% mf Number of flux zones plus 1 (max = 81)
c
c In the following list the index ib designates the beam injector,
c while ie refers to one of the three energy components.
c iap refers to one of the apertures.
c
c% nbeams Number of neutral beam injectors (1 or 2)
c% nsourc Number of sources per beamline.
c If 1, source is centered on beamline axis.
c If nsourc = 2, distinguish between the beamline
c axis and the source centerline (optical axis).
c The two sources are assumed to be mirror images
c through the beamline axis.
c In either case, the exit grid plane is perpendicular
c to the beamline axis, and contains the source
c exit grid center(s).
c If nsourc = 2, the alignment of the sources w.r.t.
c the beamline axis is specified through bhofset,
c bvofset, and bleni (described further below).
c% naptr Total number of apertures encountered by a particle
c as is moves from the source into the plasma chamber.
c Maximum is specified by parameter nap ( = 10).
c First set of apertures encountered by the particles
c are assumed centered on the source axis, and subsequent
c apertures are centered on the beamline axis;
c the distinction is made through ashape.
c% anglev(ib) Vertical angle (degrees) between optical axis
c and horizontal plane; a positive value indicates
c particles move upward
c% angleh(ib) Horizontal angle (degrees) between optical axis and
c vertical plane passing through pivot point and
c toroidal axis; a zero value denotes perpendicular
c injection, while a positive value indicates par-
c ticles move in the co-current direction
c% bvofset(ib) Vertical offset from beamline axis to center
c of each source (cm; used only for nsourc = 2)
c% bhofset(ib) Horizontal offset from beamline axis to center
c of each source (cm; used only for nsourc = 2)
c% bleni(ib) Length along source centerline (source optical axis) from
c source to intersection point with the beamline axis.
c% sfrac1(ib) Fraction of source current per beamline coming
c from upper source (used only for nsourc = 2)
c (the upper source is the more perpendicular,
c or right source normally)
c% bcur(ib) Total current (a) in ion beam (used only if bptor
c is zero)
c% bptor(ib) Total power (w) through aperture into torus; when
c nonzero, bptor takes precedence over bcur
c% bshape(ib) Beam shape
c 'circ': circular
c 'rect': rectangular
c 'rect-lps': rect. long pulse source (DIII-D only)
c a choice of short or long pulse sources is
c available by injector (not by source). one
c or both injectors may be long pulse by
c setting one or both to 'rect-lps'
c CAUTION: DIII-D sources are defaulted to lps specs.
c It is the user's responsibility to overide these for
c sps configuration(s).
c% bheigh(ib) Height of source (cm)
c Default is bshape(1) = bshape(2) = 'rect-lps'
c% bwidth(ib) Width of source (cm); diameter for circular source.
c% bhfoc (ib) Horizontal focal length of source (cm)
c% bvfoc (ib) Vertical focal length of source (cm)
c% bhdiv (ib) Horizontal divergence of source (degrees)
c% bvdiv (ib) Vertical divergence of source (degrees)
c% ebkev (ib) Maximum particle energy in source (keV)
c% fbcur (ie,ib) Fraction of current at energy ebkev/ie
c Note that this is the current fraction at the source,
c before it enters the neutralizer.
c% ashape(iap,ib) Aperture shape.
c Prefix 's-' indicates source axis centered.
c Prefix 'b-' indicates beamline axis centered.
c 's-circ' 'b-circ'
c 's-rect' 'b-rect'
c 's-vert' 'b-vert'
c 's-horiz' 'b-horiz'
c 'b-d3d'
c circ = circular aperture,
c rect = rectangular,
c vert = limits vertical height of source particles,
c horiz = limits horizontal height of source particles,
c d3d = special DIII-D polygonal aperture
c% aheigh(iap,ib) Height of aperture (cm)
c% awidth(iap,ib) Width of aperture (cm); diameter for circular
c aperture
c% alen(iap,ib) Length from source to aperture for 's-type' apertures,
c and from exit grid plane along beamline axis for
c 'b-type' apertures.
c% blenp(ib) Distance along beamline axis from source exit
c plane to the fiducial "pivot" point.
c% rpivot(ib) Radial position of pivot (cm)
c% zpivot(ib) Axial position of pivot (cm)
c
c% hdepsmth set this param to a positive value (gt.0.0) to turn off
c the smoothing of hibrz and hdepz in subroutine POSTNUB.
c if this option is used then enough zones must be specified
c for adequate resolution (zones = number of radial grid points)
c and enough injected neutrals must be followed to minimize
c the statistical noise enough so that no greatly uneven
c profiles result! this option was added because the smoothing
c of the profiles by subroutine SMOOTH can lead to unphysical
c peaking of the birth and deposition profiles.
c
c% freyavb integer switch used for debug purposes
c if freyavb .ne. 0 then some FREYA-related diagnostic
c output will be written to the screen
c ----------------------------------------------------------------------
c
do i=1,5
timbplt(i) = 1.0e6
end do
do i=1,kb
neg_ion_source(i) = 0 ! default all beam lines to std sources
end do
iterate_beam = .false.
rtstcx = 1.0
relaxden = 1.0 ! no relaxation in beam iterations
relaxden_err = 0.1 ! default is large,depends on npart
imaxbeamconvg = 5
fast_ion_target = 0 ! neglect fast ions as a beam target
beamon = 1.0e6
btime = 0.0
ibcur = 1
itrapfi = 0
use_Callen = 0
ibcx = 1
ibslow = 1
nbeamtcx = 0
fionx = 0.0
nameb = 'h'
relnub = 0.1
tfusbb = 0.0
iddcal = 3
fdbeam = 0.150e-3 ! isotopic content of d in h
c this default is for h beams
ranseed = 7**7
npart = 10000
npart_mcgo = 3000
hdepsmth = -1.0
freyavb = 0
c
c --- smoothing is normally on by above line
c
npskip = 1
iborb = 1
inubpat = 0
npat(1) = nw
npat(2) = nh
mf = 41
nbeams = 1
nsourc = 2
c
c DIII beam input
c
if (machinei .ne. 'doub-iii') go to 9020
naptr = 2
do i=1,kb
anglev(i) = 0.0
angleh(i) = 13.5
bvofset(i) = 39.75
bhofset(i) = 0.0
bleni(i) = 553.88
bcur(i) = 110.0
bptor(i) = 3.5e6
nbshape(i) = 'rect'
bheigh(i) = 10.0
bwidth(i) = 40.0
bhfoc(i) = 480.0
bvfoc(i) = 550.0
bhdiv(i) = 1.4
bvdiv(i) = 0.45
ebkev(i) = 80.0
fbcur(1,i) = 0.6
fbcur(2,i) = 0.3
fbcur(3,i) = 0.1
sfrac1(i) = 0.5
blenp(i) = 486.13
rpivot(i) = 270.0
zpivot(i) = 89.0
nashape(1,i) = 's-vert'
nashape(2,i) = 'b-horiz'
aheigh(1,i) = 10.1
awidth(1,i) = 0.0
aheigh(2,1) = 0.0
awidth(2,i) = 32.0
alen(1,i) = 442.0 * 1.0026
alen(2,i) = 456.0 * 1.0026
end do
go to 9030
c
c DIII-D beam input. DEFAULT IS LONG-PULSE-SOURCE SPECIFICATIONS.
c
9020 naptr = 4
c
do i=1,kb
anglev(i) = 0.0
angleh(i) = 19.5
bvofset(i) = 0.0
bhofset(i) = 42.074
bleni(i) = 556.808
bcur(i) = 110.0
**** bptor(i) = 10.0e6
bptor(i) = 0.0e6
nbshape(i) = 'rect-lps'
bheigh(i) = 48.0
bwidth(i) = 12.0
bhdiv(i) = 0.50
bvdiv(i) = 1.3
fbcur(1,i) = 0.7
fbcur(2,i) = 0.2
fbcur(3,i) = 0.1
bhfoc(i) = 1.0e100
bvfoc(i) = 1.0e3
ebkev(i) = 75.0
sfrac1(i) = 0.5
nashape(1,i) = 's-rect'
nashape(2,i) = 's-rect'
nashape(3,i) = 'b-d3d'
nashape(4,i) = 'b-circ'
aheigh(1,i) = 47.8
awidth(1,i) = 13.8
alen(1,i) = 186.1
aheigh(2,i) = 48.0
awidth(2,i) = 17.7
alen(2,i) = 346.0
alen(3,i) = 449.0
awidth(4,i) = 50.9
alen(4,i) = 500.0
blenp(i) = 539.0
rpivot(i) = 286.6
zpivot(i) = 0.0
end do
c
c parameters for including rotation in neutral stopping
c
9030 ds_tk = 5.0
fe_tk = 1.1 !passed to nbsgxn where it is used as ebfac
ne_tk = 20 !passed to nbsgxn where it is used as nebin
c
c --- following parameters are used in subroutine HEXNB
c
kdene = 1
kdeni = 1
kdenz = 1
ksvi = 0
ksvz = 0
ksve = 0
krad = 1
ngh = 10
ngl = 10
iexcit = 0
ilorent = 0
mstate = 4
ncont = 30
do 1710 j=1,kprim
znipm(j) = 0.0
1710 atwpm(j) = 0.0
do 1720 j=1,kimp
iz(j) = 0
izstrp(j) = 0
c
c note izstrp = 0 implies coronal equilibrium for impurity j in hexnb
c
atwim(j) = 0.0
1720 zniim(j) = 0.0
*
c ----------------------------------------------------------------------
c --- RF HEATING PARAMETERS
c ----------------------------------------------------------------------
c PARAMETERS USED BY ALL OR SEVERAL RF MODELS:
c Parameter krf (=9) gives possible number of RF models to be used.
c nth elements of quantities dimensioned by krf are to be associated
c with each other.
c% rfmode(krf) Flag identifying RF heating model
c 'ech' : electron cyclotron heating - ray tracing code
c 'ich' : ion cyclotron heating - T.K. Mau model
c 'input' : input source
c 'wedge' : input wedge-shaped source
c 'fb' : give profiles of RF current, or RF power, as input.
c current drive efficiency according to Fisch-Boozer.
c (Harvey-Freije)
c 'fastwave' : 1D or 2D fastwave heating and current drive.
c (S.C. Chiu and Bob Harvey)
c 'raytrace' : LH and fast wave raytracing model
c (adapted from Brambilla, S.C. Chiu, T.K. Mau, Bob Harvey)
c 'fastcd' : Kupfer's Ergodic Fast Wave Heating and Current Drive model
c% rfon(krf) Time (sec) at which RF heating is turned on
c% rftime(krf) Time interval (sec) during which RF heating is on
c% turnonp(krf) Time during which RF power is ramped up to full power,
c i.e., power is multiplied by
c 0.5*(1.0-COS ((time-rfon)/turnonp)). (default = 0.0)
c% turnonc(krf) Same as turnonp, except for RF current.
c% rfpow(krf) Total RF power (W)
c% rnp(krf) Parallel refractive index
c% freq(krf) Wave frequency (Hz)
c% irfcur(krf) current drive switch
c% relrf Minimum relative change in electron density and
c temperature profiles between ech calculations
c
c PARAMETERS USED BY INPUT MODEL
c% rnormin Normalized radii at which power densities are input;
c may be input in the first namelist instead
c Set following array elements corresponding to rfmode(k) = 'input'
c% njqin(krf) Number of radial points in input profiles
c% qine(kj,krf) Power density profile (W/cm**3) input to electrons
c% qini(kj,krf) Power density profile (W/cm**3) input to ions
c Note: If rfpow > 0, then qine and qini are normalized such
c that the total source power equals rfpow.
c
c PARAMETERS USED BY WEDGE MODEL
c% rfrad1(krf) Inner minor radius of wedge (cm).
c% rfrad2(krf) Outer minor radius of wedge (cm).
c rfrad1,2 are used to set up qini and qine in this subroutine.
c Unfortunately the point at which this is done uses the 1D
c rminor grid to establish the support set for qine,qine.
c The eqdsk r grid can be approximated by setting
c rminor = eqdsk rminor value
c (this value isnt know until you run the code at least once).
c% a1rf(krf) Source magnitude at rfrad1 (arb. units)
c% a2rf(krf) Source magnitude at rfrad2 (arb. units)
c% wrfe(krf) Fraction of power deposited in electrons.
c% wrfi(krf) Fraction of power deposited in ions.
c Note: the source profile is normalized such that the total
c input power equals rfpow.
c Note: The wedge model is an alternate way of setting the
c qine(i,k) and qini(i,k), variables, i = 1,nj.
c Set rfmode(k) = 'wedge'.
c Note: If irfcur(model) .ne. 0 then subroutine RFCUR_MODEL
c will be called to calculate an approximate RF
c driven current. The current drive model in this
c subroutine is only a crude approximation so
c do not use it indiscriminately HSJ!!!!!!!!!!!!!!!!!!
c Note: more than one wedge model can be active simultaneously
c thus a more realistic energy deposition can be achieved
c by combining multiple wedges HSJ
c
c PARAMETERS USED BY FISCH-BOOZER MODULE.
c% rfon,rftime,turnonp,turnonc,rfpow,rfmode,relrf, basically as above.
c% rfmode = 'fb' for fisch-boozer model, with cordey et al efficiencies.
c% lmode = 0, for lower hybrid efficiency; harmonic number (1 to 3) for ech.
c% irfcur = 1/0/-1 for co/no/counter
c% itrapech switch for modifying ech current drive by parabola
c suggested by v. chan. if itrapech = 1 and RF cur drive is
c on then the calculated RF current is multiplied by
c (1.0-eps**0.7)**0.73 where eps is inverse aspect ratio.
c the powers .7 and .73 are a crude fit to the curve
c fig.2 in cohen 'effects of trapped electrons on
c current drive',to be published in phys of fluids.
c% rnp(krf) = positive:
c parallel refractive index. used with te and ene to find efficiency
c according to the Fisch-Boozer v-parallel current drive formula.
c negative:
c absolute value is a multiplier of SQRT (2*te/me) giving parallel
c phase velocity of wave in terms of the local thermal velocity.
c% ifb= 1, specify total RF power and profile factors, and calculate
c RF power density and current profiles and total RF current.
c 2, specify total RF current and profile factors, and calculate
c RF current and power density profiles, and total RF power.
c% ifbprof= 1,parabolic profile.
c (1.0-(r**2/a**2)**alphaf(2))**alphaf(1)
c 2,guassian profile.
c EXP (-(r/(alphaf(1)*a))**2 - boundary value.
c 3,shifted guassian profile.
c EXP (-((r-alphaf(2)*a)**2/(alphaf(1)*a)**2) - boundary value,
c except, constant inside r = alphaf(3)*a at value of above
c function.
c% alphaf= as above.
c% rfcur= total current, in case current to be specified.
c% rfpow= total power, in case power profile to be specified.
c
c PARAMETERS USED BY ICH MODELS:
c% freq(krf) Frequency (1/sec) of applied wave
c% iside +1 for high field launch; -1 for low field launch
c% navg apply moving average to output radial power absorbtion
c profiles by averaging over navg adjacent points
c (total of 2*navg+1 points in average).
c
c PARAMETERS USED BY ECH (RAY-TRACING) MODEL ONLY:
c% irfplt If irfplt .ne. 0, then plots are made each time TORAY is called.
c Otherwise, the output files from the last call to TORAY can
c be processed by manually by running XPLOT after the completion
c of ONETWO. Note: if multiple TORAYs are run only the last
c TORAY run remains, not all of the TORAYs at the last time.
c The logic implementing the following definitions of IRFPLT and TIMRFP
c was found to have a bug (TCL 6/17/91). It has been disabled (CRAY40)
c and the above definition is implemented in versions later than 6/10/91.
c
c% irfplt If .ne. 0, plotted o/p from the ray tracing code is enabled.
c Plots are made at the first call to ech, and thereafter
c whenever the central temperature changes by more that
c 20%, unless timrfp(1) is < 1.0e-6.
c Plots are only made at times timrfp(i),i = 1,5, for each
c time which is less than 1.0e6.
c% timrfp(5) plotted output from ray tracing is given at these times,
c when set to less than the default value (1.0e6).
c% gafsep For codeid .ne. 'onedee', then this is the fraction
c of poloidal flux smoothed over at the plasma
c magnetic axis, and sets the mesh spacing for passage
c of output profiles from TORAY to ONETWO, i.e. 1/gafsep points.
c If = 0.0 and/or codeid = 'onedee', then the anayltic,
c circular, concentric flux surface model is used for the
c ray tracing.
c% freq(krf) Frequency (1/sec) of applied wave
c% idamp = 0 no damping calculation
c idamp = 1 near first harmonics, using Matsuda & Hsu's routine (dampga)
c idamp = 2 near second harmonics, using dampga
c idamp = 3 near third harmonics, using damprm
c idamp = 7 good for between 2nd and 3rd, using damprm
c idamp = 8 near first harmonics, using Mazzucato's routine
c for fully relativistic damping.
c Fidone, and Granata code.(Physics of Fluids, p3745, 1988)).
c (Added by Bob H. 7-26-88; see helpme in l??ech??
c idamp = negative of above numbers, then damping is as for the
c corresponding positive number, but the EQDSK is not
c reprocessed with GAFIT to form the PSICOF coefficient file.
c This saves computer time and is valid if the equilibrium has
c not changed. (the file involved is "mhddat".)
c within the ONETWO libfile).
c% xec(krf) Major radius of source.
c% zec(krf) Height above midplane (cm) of source.
c% wrfo(krf) Fraction of RF power launched as O-mode.
c Note that if wrfo > 0.95, X-mode power will be neglected,
c and if wrfo < 0.05, O-mode power will be neglected.
c
c% necsmth smoothing control parameter:
c necsmth = 0 (the default)
c = => no smoothing
c necsmth = j (a positive integer)
c = => apply j passes of a boxcar smoother
c to the eccd current and ech heating profile.
c Note that a large j will flatten out the profiles.
c We probably don't want to use j > 3
c
c% hlwec angular dispersion of rays ??? degrees ??
c% ratwec specifies antenna shape:
c ratwec = 1 => circular antenna shape
c ratwec > 1 => more elongated shape ?
c
c PARAMETERS USED BY ICH MODEL ONLY:
c% xkpar k-parallel (1/cm)
c% nhigh number of terms to retain in K
c Following two varialbes dimensioned, in view of WEDGE model, above.
c Set values appropriate to rfmode(k).
c% rfrad1ic Inner major radius (cm) for RF grid points
c% rfrad2ic Outer major radius (cm) for RF grid points
c If rfrad1ic,rfrad2ic = 0 (default) they are calculated as the
c intersection of the horizontal chord at ylaunch with:
c for onedee case- the outer flux surface;
c for other cases- the limiter surface.
c% nrfzon Number of zones in the RF grid. (Each zone may have
c a different density of grid points.)
c% nprf(1) Should always = 1
c nprf(i) (i = 2,nrfzon+1) Point number corresponding to the
c location of rfzone(i). Thus nprf(nrfzon+1) = total number
c of RF grid points, and there are nprf(i+1)-nprf(i) steps
c in the zone between rfzone(i) and rfzone(i+1).
c% rfzone(i) Location of zone boundaries for RF grid, normalized so
c that -1 corresponds to rfrad1ic, and +1 corresponds to rfrad2ic.
c Note that the RF grid does not have to stretch from rfrad1ic
c to rfrad2ic (e.g. rfzone(1) may be >-1.0).
c% relrf Minimum relative change in electron density and
c temperature and ion temperature profiles between
c ich calculations.
c% ichmod 1 to use t.k.mau code only
c 2 to use s.c.chiu code only
c 3 to use t.k.mau code when beta-ion(0) < betalm,
c and use s.c.chiu code when beta-ion(0) >= betalm
c% betalm Limit for beta-ion(0) when ichmod = 3
c% ykperp k-poloidal (not implemented)
c% ylaunch Vertical launch point for wave (cm above midplane)
c (Not implemented.)
c (Note that midplane differs for onedee and 2d cases.)
c
c PARAMETERS USED BY FASTWAVE MODULE
c% codeid, rfmode, rfon, rftime, turnonp, turnonc, rfpow, irfcur, relrf
c% nzrffw number of fw energy channels ( .le. kzrf = 11)
c% zrffw(kzrf) centerline height of the fw energy channel (cm)
c% freqfw(kzrf) frequency of each energy channel
c% rnpfw(kzrf) parallel refractive index of each energy channel
c% pzrffw(kzrf) fraction of rfpow into fw energy channel
c% htsfw(kzrf) vertical extent of RF channel at outboard, rmajor, and
c inboard positions (cm). (Same values used for each energy
c channel)
c% iswchfw 1, ttmp+landau damping (use for ICRF).
c 2, landau damping (use for LHRF).
c% nspfw number of primary ion species on which there is ion damping.
c (nspfw was eliminated by s.c.chiu on 21 feb 92)
c% lifw(kprim) starting harmonic number for damping calculations.
c% nihfw(kprim) number of harmonics calculated.
c% impath multiple absorption switch, =1 for single path absorption.
c =2 for multiple path
c
c PARAMETERS USED BY LH AND FASTWAVE, RAYTRACE MODULE:
c (rfmode = 'raytrace')
c codeid, rfmode, rfon, rftime, freq, irfcur, relrf
c% nrayptrt maximum number of steps along each ray (step size
c set elsewhere (in a raytracing input file).
c% powersrt(krt) power (Mwatt) of each n parallel interval.
c% nnkparrt(krt) number of rays representing an n_parallel interval.
c% anzinfrt(krt) lower (in absolute value) n_parallel of interval.
c% anzsuprt(krt) upper (in absolute value) n_parallel of interval.
c% nthinrt number of poloidal locations along antenna,
c represented by rays.
c (Thus number of rays is (sum of nnkparrt(krt))*nthinrt).
c% nfwsmth smoothing parameter for fast waves; meaning the same
c as necsmth for ech (0 means no smoothing,
c n .ne. 0 means n-path smoothing)
c% heightrt full height of the antenna (cms).
c% maxrefrt Maximum number of reflections from plasma edge.
c (1 for single pass for LH, 2 for single pass for FW).
c% islofart 1, for launch in LH mode, -1 for lauch in FW mode.
c ....... There are other seldom changed input variables in
c the input (command) file to the raytracing code.
c
c PARAMETERS USED BY FASTCD MODULE
c (rfmode = 'fastcd')
c codeid, rfmode, rfon, rftime, rfpow, turnonp, turnonc, irfcur, relrf,
c% nkfcd(krf) number of toroidal modes used in the calculation
c% xntor(kzrf,krf) toroidal mode numbers of the antenna spectrum
c% rpant(kzrf,krf) relative power on each toroidal mode number,
c in arbitrary units
c% gamloss(krf) phenomenological parasitic damping rate
c
c PARAMETERS USED TO SUPPLY EXTERNAL RF CURRENT FROM INONE FILE
c% extcurrf(l) master switch, if = 0.0 then model l is turned off and
c following parameters are not used. If .ne. 0.0 then
c extcurrf*extcurrf_curr is added to the currf array.
c (The currf array contains the sum total of all RF current
c sources. No time dependence of extcurrf_curr
c current is available. It is always on or always off.)
c% extcurrf_amps(l) if equal to zero then you get a total RF current
c given by the integral of extcurrf*extcurrf_curr
c if extcurrf_amps is .ne. 0.0 then the product extcurrf*
c extcurrf_curr is renormalized to give extcurrf_amps total
c current. (in this case the value of extcurrf is irrelevant
c so long as it is not 0.0)
c% extcurrf_id(l) a character variable (must be less than 80 characters)
c used to identify extcurrf_curr in output
c% extcurrf_nj(l) an integer which gives number of values in extcurrf_curr.
c if extcurrf_nj equals nj then extcurrf_rho is assumed equal
c to the internal rho (so extcurrf_rho need not be given).
c If extcurrf_nj is less than nj then it is assumed that
c extcurrf_curr
c is a spline representation and exactly extcurrf_nj values
c of NORMALIZED extcurrf_rho must be given in inone,with the
c first value exactly 0.0 and the last value exactly 1.0
c The number of values in the spline representation must be
c .le. kpsi (values of extcurrf_nj less than 33 are always ok,
c see include file mhdpar.i for exact value of kpsi )
c% extcurrf_rho(j,l) j=1,2..extcurrf_nj NORMALIZED rho values for spline
c representation of extcurrf_curr.(zero gradient at rho=0
c and natural spline at plasma edge are used)
c% extcurrf_curr(j,l) j=1,2...extcurrf_nj The RF current in amps/cm**2
c that is to be added to currf. (this will be the total
c RF current if no RF models are turned on, or will be an
c additional RF current if other RF models are in use as well)
c See subroutine GET_EXTERNAL_RFCUR for further details.
c NOTE the subscript l ranges from 1 to krf and corresponds to
c the rf models given above. Thus an external current and
c electron, iuon heating components can be input for each model
c
c AN IDENTICAL SET OF INPUTS FOR (RF) HEATING OF ELECTRONS AND IONS IS
c ALSO AVAILABLE. THE RELEVANT INPUT PARAMETERS ARE
c for electrons:
c variable analogous meaning as variable above for currf
c ------------ ---------------------------------------------
c% extqerf extcurrf
c% extqerf_rho extcurrf_rho
c% extqerf_qe extcurrf_curr
c% extqerf_nj extcurrf_nj
c% extqerf_id extcurrf_id
c% extqirf_watts extcurrf_amps
c
c extqerf_qe is input in watts/(cm**3 ) (you can input this in keV/cm**3 sec
c by setting the multiplier extqerf to 1.602e-16 joules/keV)
c extqerf_watts is input in watts
c
c for ions:
c variable analogous meaning as variable above for currf
c ------------ ---------------------------------------------
c extqirf extcurrf
c extqirf_rho extcurrf_rho
c extqirf_qe extcurrf_curr
c extqirf_nj extcurrf_nj
c extqirf_id extcurrf_id
c extqirf_watts extcurrf_amps
c
c extqirf_qe is input in watts/(cm**3 ) (you can input this in keV/cm**3 sec
c by setting the multiplier extqirf to 1.602e-16 joules/keV)
c extqirf_watts is input in watts
c
c ----------------------------------------------------------------------
c
do k=1,krf
rfmode(k) = 'no rf'
rfon(k) = 1.0e6
rftime(k) = 0.0
turnonp(k) = 0.0
turnonc(k) = 0.0
rfpow(k) = 0.0
rnp(k) = 1.0
freq(k) = 60.0e9
irfcur(k) = 0
xec(k) = rmajor - rminor
zec(k) = 0.0
wrfo(k) = 0.0
njqin(k) = 0
a1rf(k) = 0.0
a2rf(k) = 0.0
wrfe(k) = 0.0
wrfi(k) = 0.0
rfrad1(k) = 0.0
rfrad2(k) = 0.0
end do
c
necsmth = 0 ! 0 means no smoothing of qrfe and currf
irfplt = 0
iside = 1
navg = 0
gafsep = 0.02
ylaunch = 0.0
xkpar = 0.0
ykperp = 0.0
ichmod = 1
betalm = 0.0
nhigh = 4
rfrad1ic = 0.0
rfrad2ic = 0.0
relrf = 0.1
nrfzon = 1
nprf(1) = 1
nprf(2) = 101
nprf(2) = kj ! I think this should be kj ... HSJ
rfzone(1) = -1.0
rfzone(2) = +1.0
do 1750 i=1,5
1750 timrfp(i) = 1.0e6
itrapech = 0
ifb = 1
rfcur = 0.0
lmode = 0
ifbprof = 2
alphaf(1) = 0.2
alphaf(2) = 0.3
alphaf(3) = 0.2
c
nzrffw = 1
do k=1,kzrf
zrffw(k) = 0.0
pzrffw(k) = 1.0
freqfw(k) = 900.0e6
rnpfw(k) = 2.0
htsfw(k) = 20.0
end do
iswchfw = 1
**** nspfw = 1
lifw(1) = 1
nihfw(1) = 20
impath = 1
nrayptrt = 600
c
do k=1,krt
powersrt(k) = 0.0
nnkparrt(k) = 0
anzinfrt(k) = 0.0
anzsuprt(k) = 0.0
end do
c
nthinrt = 1
nfwsmth = 0 ! no smoothing
heightrt = 20.0
maxrefrt = 2
islofart = -1
nrfrad = 101
c
c initialize input parameters for FASTCD model (Y. R. Lin-liu)
c
do l=1,krf
nkfcd(l) = 1
do k=1,kzrf
xntor(k,l) = 0.0
rpant(k,l) = 0.0
end do
gamloss(l) = 0.0
end do
c
do l=1,krf
extcurrf(l) = 0.0 ! no external RF current from inone
extqerf(l) = 0.0 ! no external RF heating of electrons
extqirf(l) = 0.0 ! no external RF heating of ions
extcurrf_id(l) = 'NONE' ! external RF current
extqerf_id(l) = 'NONE' ! external RF electron heating
extqirf_id(l) = 'NONE' ! external RF ion heating
extcurrf_amps(l) = 0.0
extqerf_watts(l) = 0.0
extqirf_watts(l) = 0.0
rf_ext_curtot(l) = 0.0
rf_ext_qetot(l) = 0.0
rf_ext_qitot(l) = 0.0
cconst = 0.0
call multpl1 (extcurrf_curr(1,l),kj,cconst) ! zero for starters
call multpl1 (extqerf_qe(1,l) ,kj,cconst) ! zero for starters
call multpl1 (extqirf_qi(1,l) ,kj,cconst) ! zero for starters
end do
*
c ----------------------------------------------------------------------
c --- PELLET FUELING PARAMETERS
c ----------------------------------------------------------------------
c
c PARAMETERS USED BY OAK RIDGE PELLET MODEL (HOULBERG ET AL.)
c
c Note that this model is discontinuous with respect to plasma
c transport time: the plasma profiles are changed between time steps.
c The pellet is injected perpendicularly (along a major radius).
c
c% nampel Name of pellet species
c 'h', 'd', 't', 'dt'
c Must be primary species #1 or #2.
c% pelrad Initial pellet radius (cm)
c% vpel Pellet velocity (cm/sec)
c% timpel(i) Pellet fueling times (sec, up to 10 by increasing time)
c% nbgpel Number of fast beam ion energy groups (suggest 18)
c
c ----------------------------------------------------------------------
c
nampel = 'none'
pelrad = 0.07
vpel = 800.0e2
nbgpel = 18
do 1760 i=1,10
1760 timpel(i) = 1.0e6
*
c ----------------------------------------------------------------------
c --- FUSION FLAGS
c ----------------------------------------------------------------------
c% ifus 1, calculate sources and sinks associated with fusion
c 0, do not calculate sources and sinks associated with fusion
c NOTE THAT ifus=1 IS DEFAULT so use this switch to turn
c fusion calcs off.
c -1 allows input of 'dt' beam and 'd','t','he' thermal species.
c NOTE: USING HE AS THE THIRD ION IS NOT YET POSSIBLE.
c SELECT ONLY D AND T AS PRIMARY IONS FOR NOW:
c ifus =-1 requires that you input the following
c 1) primary ions are individual d and t fluids
c he is an optional third primary ion fluid
c 2) the two neutral species are d and t,no exceptions
c (it is not possible to model neutral he at present,
c probably never will be)
c 3) a single impurity species is present (could be he
c if he is not selected as a primary species)
c 4) electron density and zeff profiles (time dependent)
c MUST be given with inenez set to 1
c 5) if he is not a primary ion you must evolve either
c the d or t density in simulation mode.The primary
c species which is run in analysis mode is then obtained
c by charge balance
c 6) if he is a primary ion you must evolve he and one
c of the two hydrogenic species d or t.
c The remaining one will again be abtained
c by charge balance.
c 7) the code will automatically adjust zeff if necessary
c to get charge balance (if possible). It is however not
c reasonable to account for the implied dynamics that is
c associated with such a change in zeff simultaneously.
c Instead a file "zeff.mod" is created which details the
c changes that were made. You can incorporate these
c changes as appropriate in the inone file
c and do the run again.
c (non-sequential time values can occur in this file due
c to predictor/corrector iterations. Simply use the
c sequential values)
c You are responsible for managing this file between runs.
c The code will destroy any existing "zeff.mod" file and
c create a new one for the current run only if necessary.
c Hence an old "zeff.mod" file might erroneously be
c interpreted as being associated with the current run
c Note that the file begins with a date/time stamp.
c 8) the beam species must be (single fluid) dt
c with fraction of d given by fdbeam
c (if the beam is not on then ifus = -1 reproduces one
c of the other modes of operation)
c 9) at the initital time you must give the initial condition
c (i.e., initital density) of he if it is a primary ion
c 10) you must either give the initital condition for
c the hydrogenic species which will be run in simulation
c or you can specify that the ratio of the two hydrogenic
c species is zfrac (AT THE INITIAL TIME ONLY).
c If zfrac is set between 0.0 and 1.0 then it is assumed
c that you want to use zfrac.(set inenez=+1 NOT -1)
c The boundary conditions for the hydrogenic
c densities at r/a=1.0 will be determined by zfrac if
c selected or by the value of the initial profile at
c r/a=1.0 if zfrac is not used. The boundary condition
c for he is given by the initial he density profile at
c r/a=1. These boundary conditions will be assumed time
c independent unless you give non zero boundary values in
c enpb(n,k) where n corresponds to bctime(n) and k is the
c index for the hydrogenic and he densities to be
c run in simulation mode.
c (at r/a=0.0 zero gradient in rho is enforced as usual)
c taupin controls the edge neutral flux for the
c hydrogenic species run in analysis mode.
c
c% exptl_neutron_rate(j) j=1,2..nbctim,measured neutron rate, #/sec,
c as a function of bctime(j). This rate will be transferred to
c the output routines for comparison purposes (only).
c
c% beam_beam_fusion (integer)
c if ifus ne 0, then setting the integer beam_beam_fusion = 1
c will calculate the beam-beam interactions. These calculations
c take some time (about 1 min per beamline, each time NFREYA
c is called).
c Typically the beam-beam contribution will be less than a 10
c percent of the thermonuclear contribution. For this reason
c you will most likely want an indication of what the beam-beam
c rate is only at the initital and final time points. set:
c
c beam_beam_fusion = i , where i=1,2,...
c to calculate the beam beam
c reaction every i'th time that the
c beam package is called. In between
c these time the approximation
c described below is used. If you do
c not want to use the approximation
c simply set beam_beam_fusion=1 which
c will force a call to the package
c each time the beam is called.
c The exact calculations will always
c be done at t=time0.
c PROVIDED that the run makes it
c through to t=timmax the exact calc.
c will also always be done at that
c time.
c If the run terminates prematurely
c the last beam-beam (and also
c beam-thermal,see below) may be an
c exact calculation or an approximate
c result,there is no way of knowing
c a priori.
c
c The beam_beam rates are
c APPROXIMATED (at times that they are not actually computed)
c by assuming that the dependence of the beam integrated
c reaction rate has not changed (due to te evolvement for
c example- which affects the critical fast ion speed)
c from the last time the calculations were done. (Scaling with
c the appropriate fast ion densities will be done however)
c
c beam_beam_fusion = 0 will skip the calculation altogether
c (the results are stored in beam_beamddn, beam_beamdtn
c and beam_beam_tt2n)
c% beam_thermal_fusion (integer)
c same as beam_beam_fusion except for beam-thermal interactions
c iddfusb_bulk = 1 MUST BE SET (see description of
c iddfusb_bulk elsewhere)
c% iaslow 1, account for alpha particle slowing down
c 0, do not account for alpha particle slowing down
c% wtifus if wtifus > 0.0 then wtifus fraction of the alpha energy
c will go to the ions
c% fusionvb integer, if set to 1 gives some detailed output to screen
c (primarly for diagnostic use)
c
c% iddfusrate 1 use new (thermal) d(d,n)he3 rate coefficient
c (highly recommended but not the default since
c we wish to remain comensurate with other input)
c (see Bosch & Hale, Nuc. Fus. 32, No.4 (1992) 611
c
c 0 (default) use old NRL Formulary rate coefficient
c A graphical comparison of the
c two rate coefficients is available in
c /u/stjohn/ONETWOPAGES/neutron_rate/
c thermal_rate/ddnhe.?
c where ? is tex and/or ps
c
c BEAM-THERMAL ION FUSION IS MODELED BY ASSUMING THAT THE FAST ION
c DISTRIBUTION FUNCTION CAN BE REPRESENTED BY THE ANALYTICAL
c FORM, GIVEN FOR EXAMPLE IN GAFFEY, J., JOUR. PLASMA PHYS. (1976)
c VOL 16, PART 2, PG. 149
c
c iddfusb = 0 default, use the old (as presented in
c S. SCOTT's thesis),
c approximation to the neutron rate for the
c d(df,n)he3 reaction (df = fast deutron)
c Thermal motion of ions is neglected
c in this model, the effect of bulk rotation
c of thermal ions is approximated by shifting
c the injected beam energy.
c
c = 1 user selected, do a more detailed calc
c of the d(df,n)he3 rate as selected by the
c switches described below.
c (the Bosch & Hale
c cross sections are used throughout)
c ibslow NOTE if ibslow=0 is selected( see beam input
c section) the code forces iddfusb=0
c
c the following are required if iddfusb = 1:
c
c +------- iddfusb_s = 0,
c | use model which neglects bulk and thermal
c | motion of ions.
c | This model is similar to the S. SCOTT model
c | (called when iddfusb = 0, see above
c | but uses Bosch & Hale cross
c | sections and does the integral over
c | the energy of the fast ion slowing down
c | distribution numerically (i.e., essentially
c | exactly,without the approximations involved
c | in order to get an analytically integrable
c | result). The result obtained with this
c | method is correct only if the thermal ions
c | have negligible motion and hence the
c | actual neutron rate will be
c | underestimated at all temperatures. This model
c | does not attempt to account for bulk rotation of
c | thermal ions since such a result is misleading
c | unless the thermal motion is also included. The
c | intended primary use of this option is to
c | compare ONETWO results with other codes.
c |
c | ibcx ibcx is an input in the beam section used
c | to control inclusion/exclusion of fast ion
c | charge exchange. This switch affects the neutron
c | rate.
c | Note that if ibcx=0 is selected the code will
c | force iddfusb_s=0 for consistency
c |
c |--------- iddfusb_s = 1 use more refined model as determined by
c | following set of switches:
c |
c | iddfusb_bulk = 0 neglect bulk motion (NEUTRON RATE CALCS ONLY)
c | = 1 include bulk motion (default)
c | Note that this only refers to the model
c | obtained when iddfusb_s=1.
c | bulk motion is always negelected in the model
c | used when iddfusb_s=0
c | THIS SWITCH ONLY CONTROLS THE BEAM ENERGY
c | USED IN THE NEUTRON RATE CALCS. THE
c | MODIFICATION OF THE FAST ION DISTRIBUTION
c | FUNCTION {DUE TO AN ALTERED DEPOSITION
c | PROFILE OF FAST IONS CAUSED BY CHANGES IN
c | CHARGE EXCHANGE CROSS SECTIONS} IS NOT
c | AFFECTED BY THIS SWITCH. TO GET TOTAL
c | SELF CONSISTENCEY WITH NO BULK MOTION YOU MUST
c | RUN WITH NO TOROIDAL ROTATION !!!!!!!!!
c |
c | icalc_cxfactor = 0 Use this switch to control charge exchange
c | treatment. icalc_cxfactor=0 (default ) means
c | that the charge exchange probability is taken
c | as independent of the fast ion speed.
c | This assumption is made in the ONETWO fast ion
c | model and thus should also be used when the
c | neutron rate is calculated.
c |
c | = 1 Use icalc_cxfactor=1 if you want to do the
c | charge exchange probability integral instead
c | of assuming it to be constant. This method
c | is of course slower than the above method
c | and is not fully consistent with other fast ion
c | "stuff" in the code . It will give you an
c | indication how important the assumption of
c | constant cx probability is. See the ONETWOPAGES
c | docs for an exact description of this switch.
c |
c | THE FOLLOWING IS CURRENTLY NOT AVAILABLE
c |
c | ddfusb_t = le 1.0, (default = 1.0)
c | neglect neutrons produced by part of fast
c | ion distribution function that is above the
c | energy ddfusb_t*e0,where e0 is the
c | injected energy normally ddfusb_t =1.0
c | should be used. Values less than 1 are not
c | useful generally.
c | (note : the model assumes the fast ion high-energy
c | tail is due to collisions with electrons only.
c | The collision of fast ions amongst themselves,
c | indpendent of the three injection energies,
c | is neglected (i.e., the density of fast ions
c | is supposed to be small for the analytic
c | distributions used here to make sense).
c | ddfusb_t > 1.0 (user-specified)
c | account for neutrons produced by
c | part of distribution function between
c | the injected energy e0 and ddfusb_t*e0.
c | (normally,due to the rapid decay of the
c | fast ion distribution function above the injected
c | energy, ddfusb_t=1.25 should suffice)
c |
c some info and graphs are available in ONETWOPAGES
c
c ----------------------------------------------------------------------
c
ifus = 1
iaslow = 1
wtifus = 0.0
iddfusrate = 0 ! governs thermal d(d,n)he3 only
iddfusb = 0
iddfusb_s = 0
ddfusb_t = 1.0 ! neglect fast ions above inject energy
iddfusb_bulk = 1 ! include bulk motion of ions
icalc_cxfactor = 0 ! use simple cx model
beam_beam_fusion = 0 ! no beam_beam fusion
beam_thermal_fusion = 1 ! include beam_thermal fusion
fusionvb = 0
do j=1,kbctim
exptl_neutron_rate(j) = 0.0
end do
*
c **********************************************************************
c *** THIRD NAMELIST (NAMELIS3) ****************************************
c **********************************************************************
c
c ----------------------------------------------------------------------
c RUN TYPE FLAGS
c ----------------------------------------------------------------------
c% ltest_code set to 1 to do some mhd/transport coupled testing.
c the user not involved in development of the code
c should leave ltest_code = 0 (its default value).
c% use_Bp0 differentiation of the pressure can be done with
c finite differences (use_Bp0= 0, default) or by using
c the chain rule (use_Bp0=1):
c d Press /d psi = d Press /d rho * d rho/d psi
c where dhro/d psi = 1/(R0 Bp0)
c for information on what this switch does,see HSJ
c% curtype set to 0 (default) to use in determination
c of ohmic current. ( this is the old onetwo way)
c set curtype to 1 to use < J*Bt/Bt0>.
c setting curtype=1 also modifies the calculation
c of the total current contributions by accounting for
c the diamagnetic effect.
c% ifixshap 1: Fix flux surfaces in time. Perform only a single,
c initial equilibrium calculation or read the equilibrium
c from an eqdsk,see irguess below.
c 0: Evolve flux surfaces dynamically. Perform several
c equilibrium calculations.
c use ifixshap = 0 if tdem mode is selected, see below.
c ( Tdem mode is the time dependent eqdisk mode)
c% mhdmode 'no coils': The flux (in volt-sec/rad) is specified
c on the boundary of the MHD grid (using array flxeqbcd).
c the coils are not explicitly modeled in this case.
c You must simulate effects of coils outside the
c MHD grid indirectly by using appropriate boundary
c conditions. Toroidal currents (other than that due to
c the plasma) are not allowed inside the mhdgrid.
c Consequently it is usually necessary to modify
c the mhdgrid. (You MUST set xdim,ydim and redge!)
c Use 'no coils' option if tdem mode is selected (see below)
c 'coils': Coil option -- model field-shaping coils,using
c psi loop and/or magnetic probe values,see below.
c this option requires a full coil set (through the
c associated Green's table) and is implemented only
c for DIII-D type equilibria. non DIII-D runs
c must use the 'no coils' option.
c% mhdmethd set mhdmethd = 'green' to get solution of Grad-Shafranov
c equation using Green's function for entire
c MHD grid (this is the slowww method.It should not
c be used routinely. It was included here primarily
c as an aid to verify solutions obtained with the
c new cyclic reduction solver used in subroutine CMPLTCYR.)
c set mhdmethd = 'cycred' to use cyclic reduction
c
c method of inverting del-star(default).
c This method is much faster and should be used routinely.
c mhdmode = 'no coils' for fixed boundary cases
c mhdmethd = 'cycred' yields the cyclic reduction
c solution of the fixed boundary problem.
c mhdmethd = 'sorpicrd' yields the successive
c overrelaxation method.
c mhdmethd ='tdem' for time-dependent eqdsk mode
c if mhdmethd ='tdem' (i.e., the multiple eqdsk case)
c then this list of times is not used
c (the netCDF file shot_12.cdf supplies the data instead)
c The tdem mode also requires that you set the following
c switches (these are set automatically for you even if
c you set them differently or didn't set them at all)
c if mhdmethd = 'tdem' then the following settings are FORCED
c ieqdsk = 0
c mhdmode = 'no coils'
c ifixshap = 0
c also eqdskin must be set to an appropriate netCDF
c file (shot_12.cdf) which we CANNOT set for the user
c% ishot ishot should be set to shot number in netCDF file
c The tdem mode NEVER solves the MHD (i.e., Grad-Shafranov)
c equation. Instead it reads a precomputed file of MHD
c data (in netCDF format) as a function of SPACE AND TIME
c and interpolates this data for the values required at each
c time step.
c
c if mhdmode = 'no coils' and mhdmethd = 'sorpicrd' then
c% optwe optwe = .true. is used to tell the code that optimal
c external relaxation parameter is to be found. if
c if optwe = .false. then it is assumed that an optimal
c% optomegae value is input in optomegae. failing this the code will
c use a crude estimate to get optomege.
c% optwi optwi = .true. is used to tell the code that optimal
c internal relaxation parameter is to be found. if
c optwi = .false. then it is assumed that an optimal
c optomegai value is input in optomegai. failing this the code
c will use a crude estimate to get optomegai
c
c --- THE FOLLOWING IS IMPLEMENTED FOR SYMMETRIC SOR ONLY:
c
c mresidual
c nfitpoints
c splinefitbd
c fitboundary
c
c NOTE irguess specifies the source of the initial
c equilibrium. Thereafter ifixshap determines if
c further equilibria are to be calculated. Both
c the initial and subsequent equilibria depend
c on how the boundary condition for the poloidal
c flux in the Grad-Shafranov equation is specified
c using the switch mhdmode.
c if irguess .lt. 0 and ifixshap = 1 then the Grad-Shafranov
c equation will never have to be solved and boundary
c conditions are irrelevant.
c however if ifixshap = 0 AND/OR irguess .ge. 0
c then the Grad-Shafranov equation will have to be solved
c at least once and consequently you must supply
c all information required to arrive at a solution
c (in particular since this is a free boundary code
c the limiter points as well as the boundary conditions
c determine the solution)
c
c% fixfcoil INTEGER input number, pertains to treatment of
c fcoil currents in the mhdmode = 'coils' option.
c not used if mhdmode = 'no coils'.
c (also not used if ifixshap = 1 and irguess .ne. 0)
c set fixfcoil = 0 if the fcoil currents
c are to be calculated as part of the MHD solution
c procedure. This is done by simultaneously adjusting
c all fcoil currents in such a way that the value of
c psi on the psi loops, magnetic probes
c and experimental fcoil currents
c values are reproduced in a least squares sense.
c the input
c fcoil currents will be compared to the calculated
c fcoil currents in the output file and plot and
c a figure of merit describing the discrepancy
c between the calculated and input values will be
c generated if the f coil current values are input
c (in the third namelist of inone, as fcoilcur).
c
c set fixfcoil = 1 if the f coil currents are to be
c taken from the input file and therefore NOT
c CALCULATED during the MHD solution process.
c in this case the calculated psi loop and probe
c values will generally not match the input values
c The calculated and input psi loop and probe values
c will be compared in the output file and plots.
c A figure of merit describing the discrepancy
c between the calculated and input values
c will be generated for this case if the input
c values are present in the third namelist of inone.
c
c No effort other than the adjustment of fcoil
c currents described above is made to minimize
c the discrepcancy between the calculated and
c measured f coil currents or between the calculated
c and measured psi loop and probe values. We do not
c have a feedback mechanism in place to adjust
c resistivity or other parameters to force
c agreement.
c YOU SHOULD NORMALLY USE FIXFCOIL = 0
c
c if fixfcoil = 0 then the user may select the
c experimental data used to determine the fcoil
c currents as follows:
c% ifitpsi =1 means use experimental psi loop values
c (obtained by interpolation in time from the input
c values stored in flxeqbcd) in determination
c of the fcoil currents.
c ifitpsi =0 means ignore psi loop values in determination
c of fcoil currents.
c% ifitprob =1 means use experimental magnetic probe values
c (obtained by interpolation in time from the input
c values stored in expmp2) in determination
c of the fcoil currents.
c ifitprob =0 means ignore magnetic probe values in determination
c of fcoil currents.
c% fitfcur =1 means use experimental fcoil currents
c (obtained by interpolation in time from the input
c values stored in fcoilcur) in determination
c of the fcoil currents.
c fitfcur =0 means ignore experimental fcoil currents
c
c NOTE: the difference between fixfcoil and fitfcur
c is as follows. Fitfcur only has meaning if
c fixfcoil=0. For only then do we calculate fcoil
c currents. There are three sources of information
c available for determination of f coil currents:
c a) the psi loops
c b) the magnetic probes
c c) the experimentally measured f coil currents
c d) Rogowski coils (not implemented yet)
c It is generally not a good idea to set the calculated
c fcoil currents equal to the measured ones (by setting
c fixfcoil = 0). Rather the best results are abtained
c by determining the fcoil currents through a least
c squares determination. fitfcur is the switch that
c determines wheter or not the experimental f coil
c currents will be included as part of the data in
c the least squares fit.
c YOU SHOULD NORMALLY USE FITFCUR = 0 (the experimental
c fcoil currents apparently are responding to eddy
c effects not modeled in the MHD calculations so
c it is best not to use them. Quite often using
c the experimental f coil currents causses the
c MHD equilibrium not to converge).
c
c at least one of ifitpsi,ifitprob must be =1
c otherwise there is no information with which
c the fcoil currents can be determined.
c feedback from plasma resistivity or (are we
c crazy enough?) transport parameters is miles
c or years down the road.
c% iecurr =1 means include ecoil currents in MHD solution.
c iecurr =0 means exclude ecoil currents in MHD solution.
c% ivessel =1 means include vessel currents in MHD solution.
c ivessel =0 means exclude vessel currents in MHD solution.
c if the ecoil and or vessel currents are included
c in the MHD model they must be given explicitly as
c input functions of time (in arrays ecurrt and
c vesscur respectively). No models exist in ONETWO
c to allow calculation of these quantities at present.
c furthermore the vessel and ecoil currents are treated
c as exact (i.e., not measured) values.
c if the fcoil currents are fit (i.e., fixfcoil = 0) then the following
c least squares weight factors must also be supplied. These are
c necessary so that psiloop,probe and fcoil current data can be
c handeled in a way that gives reasonable weight to each. Without
c weighting the experimental fcoil currents would force the solution
c to ignore the psiloop and probe values.
c
c% errpsilp weight for psi loops
c% errmprge magnetic probes
c% errfcoil experimental f coil currents
c% errvescr experimental vessel currents
c% errecoil ecoil
c each of the above five weights is to be input
c as a fraction of the measured value. for example
c if errpsilp = 0.01 then it will be assumed that all
c psiloop values have a standard deviation of 1% of
c their input values.
c% minchisq set minchisq = 1 to roughly minimize
c chisq each time the MHD package is called.
c set minchisq = 0 to turn off this calculation.
c NOTE: the only degree of freedom we have to
c minimize chisq is the total toroidal current,
c (which enters the transport calcs only as a
c boundary condition). Hence using minchisq = 1
c has the effect of changing the boundary
c condition (on Faraday's law). Using this
c option is expensive,but may lead to better
c statistical agreement with the measured
c magnetics. However if chisq is too
c large then just changing the current
c normalization may not reduce chisq enough.
c The only recourse we have then is to modify
c the transport model (particularly Faraday's law)
c
c MOTIONAL STARK EFFECT:
c --- motional stark effect vectors are as follows:
c --- nstark max number of channels (parameter,defined in mhdpar.i)
c --- Tstark(nstark,mxtbcmhd) tangent of mag. pitch angle
c --- sigstark(nstark,mxtbcmhd) error in tstark
c --- fwtstark(nstark,mxtbcmhd) fit weight of each channel
c --- rstark(nstark,mxtbcmhd) major rad. of measurement
c --- zstark(nstark,mxtbcmhd) elevation of measurement
c --- coefficient vectors for calculation of tstark.
c --- a1-4stark(nstark,mxtbcmhd)
c% timestark(i) i = 1,2..mx of mxtbcmhd times at which MSE data are given
c% use_stark logical, if true stark effect measurements
c will be used, otherwise they will be ignored.
c use_stark = .false. is the default
c ----------------------------------------------------------------------
c
optwe = .false.
optwi = .false.
optomegai = 0.0
optomegae = 0.0
splinefitbd = .false.
fitboundary = .false.
ltest_code = 0
ifixshap = 1
curtype = 0
fixfcoil = 0
fitfcur = 0
errpsilp = 0.03
errmprbe = 0.05
errfcoil = 0.10
errvescr = 0.10
errecoil = 0.10
ifitpsi = 1
ifitprob = 1
iecurr = 1
ivessel = 0
use_Bp0 = 0
call zeroa (fwtstark,mxtbcmhd*nstark)
call zeroa (rstark,mxtbcmhd*nstark)
call zeroa (zstark,mxtbcmhd*nstark)
call zeroa (tstark,mxtbcmhd*nstark)
call zeroa (timestark,nstark)
call zeroa (sigstark,mxtbcmhd*nstark)
call zeroa (a1stark,mxtbcmhd*nstark)
call zeroa (a2stark,mxtbcmhd*nstark)
call zeroa (a3stark,mxtbcmhd*nstark)
call zeroa (a4stark,mxtbcmhd*nstark)
use_stark = .false.
mhdmode = 'coils'
mhdmethd = 'cycred'
*
c ----------------------------------------------------------------------
c --- GEOMETRIC PARAMETERS
c ----------------------------------------------------------------------
c% xdim Width of box (cm) when MHDMODE = 'no coils'
c% ydim Height of box (cm) when MHDMODE = 'no coils'
c% redge Inside radius of box (cm) when MHDMODE = 'no coils'
c Note: when using the 'no coils' option, xdim, ydim
c and redge define the computational MHD grid
c vectors rmhdgrid(i), zmhdgrid(j), i = 1,..nw,
c j = 1,..nh. Since there are no coil currents
c in this option it usually does not make sense
c to define xdim,ydim,redge in the same way as
c is done with the coils option. Instead you will
c probably find that the boundary of the grid
c must be very close to the plasma,in order that
c the interpolation between boundary values makes
c sense and that shape control is possible.
c don't forget that the limiter still limits the
c radial and vertical extent of the plasma. In
c particular setting up an mhdgrid with xdim,ydim,
c redge means that you will probably also have to
c define a consistent limiter point set as well.
c if mhdmode = 'no coils' xdim,ydim,redge MUST be
c input (there is no default here).
c if mhdmode = 'coils' xdim,ydim,redge will not
c be used(the grid will be generated from the
c eqdsk or the Green's table as appropriate).
c
c% nlimiter Number of points defining limiter (max = maxlimpt)
c the limiter can be specified in one of three ways:
c a)in the third namelist of file inone. to use
c this option set nlimiter equal to the number
c of (xlimiter,ylimiter) points given in the
c third namelist of inone (and specify nlimiter
c values of xlimiter,ylimiter).
c b)the default values given below may be used.
c leave nlimiter UNSPECIFIED and DO NOT include
c the (xlimter,ylimiter) points in file inone to get
c this option.
c c)the limiter points are read from the eqdsk.
c set nlimiter equal to a negative number
c in the third namelist of inone to get this option.
c you must of course supply an eqdsk for this option.
c (most eqdsks have limiter points in them. If it
c does not a warning will be printed and the code
c will stop). Use this option for 'tdem' mode also
c (see RESTART PARAMETERS below for eqdsk specification)
c Note: When running with an eqdsk (i.e., irguess .ne. 0) it is recommended
c that option c be used since this will presumably be the
c correct limiter position for the equilibrium to be
c analyzed. If MHD calculations are done in ONETWO
c the results WILL DEPEND on the limiter position
c if the plasma touches the limiter at any time during
c the analysis. It is up to you to make sure that
c the limiter points used are correct! Consequently
c option b,above, should be used only if a fixed
c equilibrium is used,since in that case,the limiter
c position is irrelevant.
c It is now possible to generate the initial
c equilibrium in ONETWO so no eqdsk exists a priori.
c (this is done using irguess = 0,see below).
c in this case option a,or c can be used.
c if option c is used then the eqdsk can be any
c eqdsk which has the appropriate limiter points.
c all other information in the eqdsk will be ignored
c provided that irguess = 0.
c [IRGUESS = 0 IS CURRENTLY NOT IMPLEMENTED]
c
c% xlimiter(i) List of x coordinates for limiter (m)
c% ylimiter(i) List of y coordinates for limiter (m)
c% ifill 1, Fill in additional limiter points
c 0, Do not fill in additional limiter points
c ifill is not used
c ifill is ignored in input file, no filling of
c limiter will be done
c% greentab Name of Green's function table
c unless the Green's table name was changed from
c the standard defaults, the code will determine
c the proper name if greentab is not specified
c in the input. The default name will change with
c the equilibrium grid size. It consists of 8
c characters, starting with a t, followed by number
c of grid points in r direction, nw, followed by number
c of grid points in z direction, nh, followed by
c however many characters of 'd3d' required to
c complete an 8-character name. example:
c 't65129d3' would be the default name for
c for a 65 by 129 Green's table for d3d. No option
c for machines other than d3d exists.
c NOTE YOU CAN ONLY USE A GREEN'S TABLE WHICH
c MATCHES IN SIZE (NW,NH,NCOILS,ETC)
c THE PARAMETERS SET IN ONETWO. A CHECK FOR THIS
c IS MADE IN SUBROUTINE READGREN AND ERROR MESSAGE
c WILL BE ISSUED IF INCONSISTENCIES ARE DETECTED.
c A GREEN'S TABLE IS REQUIRED WHEN
c mhdmode ="coils" and irguess .ge. 0
c (independent of ifixshap)
c if irguess < 0 and ifixshap = 1, then a Green's table is not required
c NO GREEN'S TABLE IS REQUIRED IF MHDMODE ="NO COILS".
c HOWEVER YOU MUST SUPPLY BOUNDARY CONDITIONS
c FOR PSI IN THIS CASE (SEE MHDMODE)
c ----------------------------------------------------------------------
c
if (machinei .ne. 'doub-iii') go to 9080
xdim = 0.0
ydim = 0.0
redge = 0.0
greentab = 'default'
ifill = 1
nlimiter = 69
c
xlimiter( 1) = 0.961
xlimiter( 2) = 0.961
xlimiter( 3) = 0.961
xlimiter( 4) = 0.961
xlimiter( 5) = 0.961
xlimiter( 6) = 0.961
xlimiter( 7) = 0.961
xlimiter( 8) = 0.961
xlimiter( 9) = 0.961
xlimiter(10) = 0.961
xlimiter(11) = 0.961
xlimiter(12) = 0.961
xlimiter(13) = 0.961
xlimiter(14) = 0.961
xlimiter(15) = 0.961
xlimiter(16) = 1.134
xlimiter(17) = 1.300
xlimiter(18) = 1.400
xlimiter(19) = 1.500
xlimiter(20) = 1.600
xlimiter(21) = 1.680
xlimiter(22) = 1.908
xlimiter(23) = 1.850
xlimiter(24) = 1.855
xlimiter(25) = 1.859
xlimiter(26) = 1.861
xlimiter(27) = 1.861
xlimiter(28) = 1.860
xlimiter(29) = 1.858
xlimiter(30) = 1.854
xlimiter(31) = 1.849
xlimiter(32) = 1.842
xlimiter(33) = 1.833
xlimiter(34) = 1.823
xlimiter(35) = 1.812
xlimiter(36) = 1.799
xlimiter(37) = 1.784
xlimiter(38) = 1.767
xlimiter(39) = 1.749
xlimiter(40) = 1.729
xlimiter(41) = 1.707
xlimiter(42) = 1.724
xlimiter(43) = 1.666
xlimiter(44) = 1.624
xlimiter(45) = 1.624
xlimiter(46) = 1.666
xlimiter(47) = 1.724
xlimiter(48) = 1.689
xlimiter(49) = 1.712
xlimiter(50) = 1.733
xlimiter(51) = 1.752
xlimiter(52) = 1.770
xlimiter(53) = 1.785
xlimiter(54) = 1.799
xlimiter(55) = 1.812
xlimiter(56) = 1.822
xlimiter(57) = 1.832
xlimiter(58) = 1.839
xlimiter(59) = 1.845
xlimiter(60) = 1.850
xlimiter(61) = 1.853
xlimiter(62) = 1.855
xlimiter(63) = 1.855
xlimiter(64) = 1.854
xlimiter(65) = 1.851
xlimiter(66) = 1.847
xlimiter(67) = 1.908
xlimiter(68) = 1.680
xlimiter(69) = 1.134
c
ylimiter( 1) = -1.300
ylimiter( 2) = -1.000
ylimiter( 3) = -0.900
ylimiter( 4) = -0.800
ylimiter( 5) = -0.600
ylimiter( 6) = -0.400
ylimiter( 7) = -0.200
ylimiter( 8) = 0.000
ylimiter( 9) = 0.200
ylimiter(10) = 0.400
ylimiter(11) = 0.600
ylimiter(12) = 0.800
ylimiter(13) = 0.900
ylimiter(14) = 1.000
ylimiter(15) = 1.300
ylimiter(16) = 1.455
ylimiter(17) = 1.455
ylimiter(18) = 1.455
ylimiter(19) = 1.455
ylimiter(20) = 1.455
ylimiter(21) = 1.455
ylimiter(22) = 1.220
ylimiter(23) = 1.220
ylimiter(24) = 1.170
ylimiter(25) = 1.120
ylimiter(26) = 1.070
ylimiter(27) = 1.020
ylimiter(28) = 0.970
ylimiter(29) = 0.920
ylimiter(30) = 0.870
ylimiter(31) = 0.820
ylimiter(32) = 0.770
ylimiter(33) = 0.720
ylimiter(34) = 0.670
ylimiter(35) = 0.620
ylimiter(36) = 0.570
ylimiter(37) = 0.520
ylimiter(38) = 0.470
ylimiter(39) = 0.420
ylimiter(40) = 0.370
ylimiter(41) = 0.320
ylimiter(42) = 0.320
ylimiter(43) = 0.230
ylimiter(44) = 0.064
ylimiter(45) = -0.064
ylimiter(46) = -0.230
ylimiter(47) = -0.320
ylimiter(48) = -0.320
ylimiter(49) = -0.370
ylimiter(50) = -0.420
ylimiter(51) = -0.470
ylimiter(52) = -0.520
ylimiter(53) = -0.570
ylimiter(54) = -0.620
ylimiter(55) = -0.670
ylimiter(56) = -0.720
ylimiter(57) = -0.770
ylimiter(58) = -0.820
ylimiter(59) = -0.870
ylimiter(60) = -0.920
ylimiter(61) = -0.970
ylimiter(62) = -1.020
ylimiter(63) = -1.070
ylimiter(64) = -1.120
ylimiter(65) = -1.170
ylimiter(66) = -1.220
ylimiter(67) = -1.220
ylimiter(68) = -1.455
ylimiter(69) = -1.455
go to 9090
c
c --- d3d limiter points
c
9080 xdim = 1.70
ydim = 3.20
redge = 0.84
greentab = 'default'
ifill = 1
nlimiter = 33
c
xlimiter( 1) = 1.0139
xlimiter( 2) = 1.0139
xlimiter( 3) = 0.9989
xlimiter( 4) = 0.9989
xlimiter( 5) = 1.1378
xlimiter( 6) = 1.6703
xlimiter( 7) = 1.9055
xlimiter( 8) = 2.1406
xlimiter( 9) = 2.2035
xlimiter(10) = 2.3293
xlimiter(11) = 2.4001
xlimiter(12) = 2.3480
xlimiter(13) = 2.3590
xlimiter(14) = 2.3800
xlimiter(15) = 2.3820
xlimiter(16) = 2.3800
xlimiter(17) = 2.3590
xlimiter(18) = 2.3480
xlimiter(19) = 2.4001
xlimiter(20) = 2.3289
xlimiter(21) = 2.2024
xlimiter(22) = 2.1391
xlimiter(23) = 1.8221
xlimiter(24) = 1.8221
xlimiter(25) = 1.5630
xlimiter(26) = 1.5630
xlimiter(27) = 1.2750
xlimiter(28) = 1.2750
xlimiter(29) = 1.1440
xlimiter(30) = 1.0201
xlimiter(31) = 0.9989
xlimiter(32) = 0.9989
xlimiter(33) = 1.0139
c
ylimiter( 1) = 0.0000
ylimiter( 2) = 0.3998
ylimiter( 3) = 0.4148
ylimiter( 4) = 1.2443
ylimiter( 5) = 1.3832
ylimiter( 6) = 1.3832
ylimiter( 7) = 1.1930
ylimiter( 8) = 1.0027
ylimiter( 9) = 0.8501
ylimiter(10) = 0.5449
ylimiter(11) = 0.3922
ylimiter(12) = 0.3260
ylimiter(13) = 0.2220
ylimiter(14) = 0.1120
ylimiter(15) = 0.0000
ylimiter(16) = -0.1120
ylimiter(17) = -0.2220
ylimiter(18) = -0.3260
ylimiter(19) = -0.3922
ylimiter(20) = -0.5458
ylimiter(21) = -0.8529
ylimiter(22) = -1.0064
ylimiter(23) = -1.3832
ylimiter(24) = -1.3682
ylimiter(25) = -1.3682
ylimiter(26) = -1.3832
ylimiter(27) = -1.3832
ylimiter(28) = -1.3682
ylimiter(29) = -1.3682
ylimiter(30) = -1.2443
ylimiter(31) = -1.2443
ylimiter(32) = -0.4148
ylimiter(33) = -0.3998
*
c ----------------------------------------------------------------------
c --- BOUNDARY CONDITIONS
c ----------------------------------------------------------------------
c% deltat_fixed_boundary(mxtbcmhd)
c time intervals (in sec) between fixed boundary
c equilibrium calculations. Note that this input is valid
c only for the fixed boundary case( free boundary and
c tdem cases must input the times explicitely in timeqbcd)
c the first equilibrium is at time0 (but see irguess which
c could eliminate the first equilibrium time) the last
c equilibrium is at time timmax PROVIDED that mxtbcmhd
c times deltat_fixed_boundary plus time0 is large enough
c to span this time interval. (if it is not an error exit is
c taken) . If some values are set in timeqbcd then
c deltat_fixed_boundary(1) becomes the first interval
c after the last time given in timeqbcd.
c default deltat_fixed_boundary(j)=0.0,j=1...mxtbcmhd
c
c% timeqbcd(m) the times for equilibrium
c calculations can be (must be if deltat_fixed_boundary
c is not set ,see above) input in this vector.
c timeqbcd is a
c list of times (sec) for time-dependent boundary
c conditions on the magnetic flux. The maximum
c number of times allowed is mxtbcmhd ( = 1000)
c if mhdmode ='no coils' then this is a list of times
c at which equilibrium calculations are to be done
c (and possibly eqdsks written). the code may do
c additional MHD calculations, if convergence is
c not achieved.
c NOTE: if mhdmethd ='tdem' (i.e., the multiple eqdsk case)
c then this list of times is not used
c (the netCDF file shot_12.cdf supplies the data instead)
c
c% flxeqbcd(i,m) Flux boundary conditions for MHD equilibrium.
c index m corresponds to time timeqbcd(m) and index i
c corresponds to the psi loop values if mhdmode = coils
c For mhdmode = 'no coils', flxeqbcd is not used.
c
c The no coils option solves a fixed boundary equilibrium
c problem. the plasma boundary is read from the eqdsk,the
c mhdgrid (rmhdgrid,zmhdgrid) is calculated from the eqdsk
c and initially psincbcd is set to the boundary values of
c psi on the MHD grid as follows:
c
c nh > > > > > nh+nw-1
c ^ v
c ^ v
c ^ v
c ^ v
c 1 < < < < < nh+nh+nw-1
c
c That is, psincbcd has the flux (in v-sec/rad) stored
c (for time point 1) starting at the lower left corner
c and progressing clockwise. there are
c 2*(nh+nw-2) values of psi stored in psincbcd in all.
c after the initial time, psincbcd is internally set to
c whatever values are required in order to force the given
c plasma boundary to be a flux surface. thus no boundary
c condition input is required from the user.
c Point 1 above corresponds to (rmhdgrid(1),zmhdgrid(1)),
c point nh+nw-1 corresponds to (rmhdgrid(nw),zmhdgrid(nh)) etc.
c (if you want to change the number of grid points from nw,nh
c to something else you will have to recompile the code. BUT
c nh must still be 2**m+1 for some integer m, nw is arbitrary).
c% curmax(i) maximum current (amps) that coil number i can live with.
c (ignored in this version of the code.)
c ----------------------------------------------------------------------
c
9090 icoil = 0
nfcoilmax = nfcoil
if(nflxeqbcd .lt. nfcoil)nfcoilmax=nflxeqbcd
do 2610 i=1,nfcoilmax
curmax(i) = 1.0e30
do 2610 j=1,mxtbcmhd
2610 flxeqbcd(i,j) = 0.0
do 2620 i=1,mxtbcmhd
deltat_fixed_boundary(i)=0.0
2620 timeqbcd(i) = -1.e30
*
c ----------------------------------------------------------------------
c --- SOLUTION CONTROL PARAMETERS
c ----------------------------------------------------------------------
c% iteq Maximum number of inner iterations to obtain a free
c boundary equilibrium for a given current density
c (suggest 55)
c% itcur Maximum number of outer iterations to obtain a self-
c consistent current density, i.e., to find an ff' consistent
c with the p' and obtained from transport
c (suggest 3)
c% toleq Tolerance for inner iteration loop (suggest 1.0e-5)
c% tolcur Tolerance for outer, current iteration loop (suggest 1.0e-3)
c% omeq Relaxation parameter for inner iterations (suggest 0.8)
c set omeq = 0.0 to get dynamic relaxation (recommended)
c% omcur Relaxation parameter for outer, current iterations
c omcur is used in subroutine PFPRIM (cray209.f) to relax
c pprim. pprim is defined over the psir(1...nj) grid
c (suggest 0.5)
c% isym 1: Symmetrize solution
c 0: Do not symmetrize solution
c note: at the present time the code makes
c the solution up/down symmetric about the grid point
c nh/2+1 by first calculating the total solution and
c then overwriting the lower half of the solution with
c a copy of the upper half.
c% npsi Number of flux surfaces on which the averages are done
c (suggest 51, the same as standard transport grid, max=kpsi
c% ivertsbl switch for vertical stability
c set ivertsbl = 0 to disable the internal coding
c that tries to stop vertical drift of the plasma
c if it is sensed to occur. set ivertsbl = 1 if vertical
c drift is a problem and stabilization may help.
c (this is a real problem for low current highly
c elongated plasmas. EFIT also sees this drift)
c The method used can be found in subroutine FREEBDRY.
c IVERTSBL = 1 HAS AN EFFECT ONLY IF FIXFCOIL=0
c USE ONLY IVERTSBL = 0 UNTIL FURTHER NOTICE. HSJ
c% tensionspl tension spline parameter. if not equal to 0.0
c then a tension spline will be used to interpolate
c q from the MHD to the transport grid. to get
c infinite tension set tensionspl = 1.0e30. otherwise
c reasonable values range from 1 to 2000.
c Set tensionspl to a negative number, to let the code
c decide what tension to use.
c NOTE: normally you want to use tensionspl = 0
c tensionspl is intended only for problem cases where
c q should be monotonic but isnt due to rapid variaton
c of the spline representation of q near the boundary.
c% extend_seval integer,set to 1 to allow spline interpolation outside
c the interpolating table in subroutine seval (not a good
c idea in general) At this time seval is used only to
c interpolate beam stopping cross sections for the ADAS
c cross section set !!!! Any attempt to go outside the
c range of the table will cause the limiting value in
c the table to be returned. THERE IS NO EXTRAPOLATION OF THE
c TABLE !!!!!
c Default is 0 which stops the code when an out of
c bounds interpolation is performed.
c
c ----------------------------------------------------------------------
c
extend_seval = 0
iteq = 60
itcur = 3
ivertsbl = 0
toleq = 1.0e-6
tolcur = 1.0e-4
omeq = 0.75
omcur = 0.5
isym = 0
npsi = kpsi ! parameter kpsi is set in MHDPAR
tensionspl = 0.0 ! don't use tension spline interpolation
*
c ----------------------------------------------------------------------
c --- DYNAMIC CONTROL PARAMETERS
c ----------------------------------------------------------------------
c% ieqmax Number of additional equilibria to be calculated, at times
c given by deqlst(j),j = 1,2..ieqmax . for the fixed
c boundary case (mhdmode = 'no coils') ieqmax should
c be set to zero. for the free boundary case
c (i.e., mhdmode ='coils') ieqmax may be zero or not
c (see below).
c% deqlst(i) i = 1,2...ieqmax. list of times(sec) for additional equilibria
c calculations.
c If no additional equilibria are required then
c deqlst need not be set and ieqmax should be set to 0.
c NOTE
c equilibria are always calculated at times
c given by timeqbcd(i), i = 1,2..itbcmhd.
c deqlst provides a way to get additional equilibria
c at times other than those at which the MHD
c boundary conditions are specified. If ibypass
c is set to 0 then the code may perform additional
c equilibrium calculations because the time step
c between equilibria may be reduced due to lack
c of convergence. an absolute maximum of keqmax
c equilibria may be calculated. if keqmax is exceeded
c the code will stop.
c (keqmax is not an input, it is a fixed parameter).
c for the fixed boundary case (i.e., mhdmode = 'no coils')
c deqlst is redundant and hence is not used! instead
c timeqbcd is used to determine when equilibria are
c to be calculated (again the code may calculate
c equilibria at additional times if necessary).
c Note that in the fixed boundary case timeqbcd is
c the time at which equilibria are calculated rather
c than the time at which boundary conditions are given
c and (consequently) equilibria are calculated.
c An eqdsk with name g.time may be written for each
c of the times given in timeqbcd and in deqlst.
c These eqdsks may be used
c to restart the transport code, see restart
c parameters, and setting of ieqdsk.
c The eqdsks generated are the same as EFIT eqdsks
c but contain additional information appended to the
c the normal EFIT eqdsk.
c% delcap Maximum relative change allowed in geometric factors
c before and after an equilibrium calculation
c (suggest 0.10)
c% xi_include logical, set to .true. if eps, xhm2, xi11, xi33,
c and xips are to be included under delcap
c convergence criterion.
c% delrho Maximum relative change allowed in the value of rho
c at the plasma boundary before and after an equilibrium
c calculation (suggest 0.005)
c% maxitr Maximum number of passes made over a single transport/mhd
c time step, to converge rho on plasma
c boundary and geometric parameters used in transport, suggest 3
c% ibypas Bypass flag
c 0: equilibrium time step is halved down if convergence
c criteria on delcap and delrho are not met after maxitr
c passes
c 1: convergence test is bypassed after maxitr passes
c NOTE: RUNNING WITH IBYPAS = 1 SHOULD BE DONE ONLY
c WITH THE INTENTION OF "LETS SEE WHAT HAPPENS".
c IT DOES NOT IMPLY THAT THE RESULTS WILL BE CORRECT !
c% dtmine Minimum time step between equilibrium calculations
c% psifctr controls the outermost psi surface value at which
c FLUXAV does the contour integral calculations.
c set psifctr = 0.0 to use the boundary value of psi.
c {Note that this boundary value of psi
c could be different than
c the original value of psilim
c on the eqdsk because pol_flux_lim,
c see below, was used.}
c psifctr = 0.01 means boundary will be pulled in
c 0.01*(delta psi). where delta psi is (psimag-psilim)/(npsi-1)
c since npsi is an input param you can use it in conjunction
c with psifctr to determine the outermost flux surface to use.
c reasonable values for psifctr are 0.0 to about 0.9 (if
c you set psifctr = 1 then the last two values of psi on
c the grid become identical so the code won't let you
c set psifctr greater than 0.99)
c
c% pol_flux_lim A switch introduced so that onetwo could match other codes
c in the treatment of the plasma boundary. The value of psi
c at the plasma boundary is taken as
c psilim = psiaxis + pol_flux_lim*(psilim_eqdsk-psiaxis)
c where psiaxis is the mag axis psi value, psilim_eqdsk is
c the value of the plasma boundary as specified in the eqdsk.
c (This switch is not availble in tdem mode.)
c Default value is pol_flux_lim = 1.0
c% derwght a relaxation parameter for the cap (i.e., fcap, gcap,
c hcap) parameters. derwght is defaulted to -0.5.
c which means 1/2 old plus 1/2 new estimate is
c taken. to get mostly old estimate use
c derwght = -0.95, to get mostly new use derwght = -0.05
c any value between -1.0 and 0.0 is valid.
c% cap_mult a multiplier for the time derivatives
c dgdt,dhdt,dfdt,dr2dt,drcapdt,dr2idt,depsdt,
c dxhm2dt,dxi11dt,dxi33dt,dxipsdt,drhoa_dt_geom
c most users should leave cap_mult set at its default
c value of 1.0 (it is used for code verification purposes only)
c
c% adotmult multiplier for adot. defaulted to 1.0
c
c% f2d1mult
c% f2d2mult
c% f2d3mult multipliers for 2d sources in Faraday's Law
c
c% use_cnt1 logical
c% use_cnt2 logical: set only one of use_cnt1 or use_cnt2 to true.
c set the other to false.
c if use_cnt1 is true then the modified CNTOUR routine,
c which finds boundary points using steps perpendicular
c to the gradient of psi, is used.
c if use_cnt2 is true then the ray search method is used.
c ----------------------------------------------------------------------
c
use_cnt1 = .true.
use_cnt2 = .false.
timlong = 1.0e100
ieqmax = 1
do 2520 i=1,keqmax
2520 deqlst(i) = timlong
maxitr = 2
delcap = 0.01
delrho = 0.001
derwght = -0.5
cap_mult = 1.0
adotmult = 1.0
f2d1mult = 1.0
f2d2mult = 1.0
f2d3mult = 1.0
ibypas = 0
xi_include = .false.
psifctr = 0.0
pol_flux_lim = 1.0
*
c ----------------------------------------------------------------------
c --- STARTUP/RESTART PARAMETERS
c ----------------------------------------------------------------------
c% irguess specifies the source of the initial equilibrium
c AND SOURCE OF INITIAL CONDITIONS FOR TRANSPORT EQUATIONS.
c 0: initial equilibrium is not available. Therefore
c construct this equilibrium from scratch,
c using current and pressure profiles in the
c first namelist of inone.(An internal guess routine
c is used to generate the initial equilibrium. The
c starting shape for psi is controlled with rmagax,
c zmagax,rcsale,zscale,see below. The initial conditions
c for the transport equations are obtained from file inone
c and are used to generate the current profile required
c to solve the Grad-Shafranov equation. An iterative solution
c of the Grad-Shafranov equation is required because the initial
c transport profiles are given as a function of rho
c which is not known until psi(r,z) is established. Psi
c also depends on boundary conditions supplied for the
c Grad-Shafranov equation (in the third namelist of inone) so in
c effect we are trying to force the psi to rho mapping
c to be consistent with the current profile and the MHD
c boundary conditions (by forcing the poloidal current
c function f(psi) to the required values). This may or
c may not lead to a converged solution. If convergence
c is not achieved you may have to alter the current profile
c in inone and/or the boundary conditions for the Grad-Shafranov
c equation. It is also possible (but not likely) that
c a different initial guess for psi will lead to convergence,
c (see rmagax,zmagax,rscale,zscale below).
c DO NOT USE IRGUESS = 0 UNTIL FURTHER NOTICE! HSJ
c = 1: read initial equilibrium data from file eqdskin
c and do initial equilibrium calculation.
c (this should reproduce the equilibrium read in,
c and serves as a check on the input equilibrium)
c = 2: read p prime
c and f-fprime from eqdsk. construct curden from
c this pprime and ffprime. then use kinetic data
c from first namelist of inone to get new p-prime.
c use new p-prime together with curden to generate
c an f-fprime that is consistent with curden
c calculated from eqdsk and kinetic pressure profiles.
c = -1: read initial equilibrium data from file eqdskin;
c use input equilibrium for initial time point,
c and read the parameters rmajor, btor, and totcur from
c eqdskin if their values are zero in the first
c namelist . If the values are not zero in the
c first namelist then rescale the initial equilibrium
c accordingly. Use this initial equilibrium to generate
c the MHD dependent geometric coefficients for the
c transport equations directly.
c The initial conditions required for transport equations
c are obtained from the first namelist of inone.
c (this is the method which was used in the
c past to make time-independent 1.5 d runs,WITH THE
c FOLLOWING IMPORTANT EXCEPTION:
c The current profile,curden,is no longer an input
c quantity. The current profile is determined from
c the eqdsk so that it is consistent with the geometric
c parameters (mostly gcap). In the old version the eqdsk
c current profile was used only if you specified the
c input current profile as a parabola. If you used
c spline input for curden,the eqdsk current profile
c given by pprim and ffprim was not used.
c
c If irguess = 0 you have some limited control over the
c intial guess used for psi,by specifying the following
c parameters:
c% rmagax rmagax and zmagax are guesses for the location of
c% zmagax the magnetic axis (in cm,must be consistent with
c definition of mhdgrid.)
c% rscale rscale and zscale are scaling factors for the
c% zscale intial psi guess.
c The initial shape for psi is
c
c psi(R,z) = a * EXP (-(R-RMAGAX)**2/RSCALE-(Z-ZMAGAX)**2/ZSCALE)
c
c the constant a is calculated by the code (so that
c it is consistent with totcur)
c if rmagax,zmagax,rscale,zscale are not set but are
c required because irguess = 0 the code will pick
c appropriate? values.
c
c% mhdonly =0 no effect
c =1 run only a single equilibrium calculation, do not
c do any transport! This option was created primarily
c to check out MHD calculations in the context of the
c transport code. It may be useful in certain
c circumstances since it allows more general
c specification of current profiles than does EFIT.
c if a large number of cases with mhdonly = 1
c are to be run you may wish to make a
c special version of ONETWO that does not
c include the transport section. to do this
c simply compile and link the code using only
c cray10,cray20 and cray40. the resulting
c executable will have some unsatisfied external
c references (which are not used) but will
c nevertheless run correctly.
c
c% eqdskin Name of input eqdsk file (default = 'eqdskin')
c NOTE: If this is a multiple eqdsk run (mhdmethd = 'tdem') then
c eqdskin is the name of the data file to be read, NOT the
c name of any particluar eqdsk. The data file is created
c by the standalone code MEPC. the file name is usually
c "shot_12.cdf"
c% use_efit_cntr set to 0 to find the plasma boundary from the information
c on the eqdsk.
c set to 1 to use the cntour given in the eqdsk rather than
c finding it again
c ----------------------------------------------------------------------
c
irguess = 0
mhdonly = 0
eqdskin = 'eqdskin'
rmagax = 0.0
zmagax = 0.0
rscale = 0.0
zscale = 0.0
use_efit_cntr =0
*
c ----------------------------------------------------------------------
c --- STAEBLER-HINTON MODEL
c ----------------------------------------------------------------------
c
c Input variables aeh, ael, aih, ail, bh, bl, ch, cl, alfae, alfai,
c betah,sigma and gammah are defined by the following model
c developed by G.M. Staebler and F.L. Hinton*:
c
c (xe is electron energy difusivity)
c Xe = A * Xe where A = aeh + ael
c ---------------------------
c 1 + (alfae * sperp)**gammah
c
c (xi is ion energy diffusivity)
c Xi = B * Xi where B = aih + ail
c ---------------------------
c 1 + (alfai * sperp)**gammah
c
c (Di is particle diffusivity)
c Di = C * Di where C = bh + bl
c ---------------------------
c 1 + (betah * sperp)**gammah
c
c ( ui is momentum diffusivity)
c ui = D * Xi where D = ch + cl
c ---------------------------
c 1 + (sigma * sperp)**gammah
c
c If the flow shear suppression flag, ifsflag, is set to 0 the
c model will not be used; if it is set to 1, the coefficients
c A,B,C,D, will be applied . At the present time
c this model acts only on the RLW and HSIU models of diffusivity.
c Neoclassical values of xe,xi,di,ui are added to the above turbulent
c values after the above modification is made. Thus for example the
c total momentum diffusivity is
c ui = D*xi +uneo
c where uneo is the neoclassical momentum diffusivity ,D is the
c multiplier defined above and xi is either the RLW or HSIU model
c of ion energy diffusivity.
c
c The constant variables are saved in array fs in the following
c order:
c
c% fs( 1) = ifsflag
c fs( 2) = aeh
c fs( 3) = ael
c fs( 4) = aih
c fs( 5) = ail
c fs( 6) = bh
c fs( 7) = bl
c fs( 8) = ch
c fs( 9) = cl
c fs(10) = alfae
c fs(11) = alfai
c fs(12) = betah
c fs(13) = sigma
c fs(14) = gammah
c fs(15) = lsfctr
c% lsfctr is a switch used in determination of sperp
c (see subroutine csperp for actual defn of sperp)
c if lsfctr = 0 then the magnetic shear scale length
c used in sperp is given by
c ls = (R*q**2)/( r*dq/dr)
c if lsfctr = 1 then we take the constant length
c ls = Rmajor
c
c fs(16) = xrot * Multiplier of toroidal frequency derivative
c component in Sperp term. Set to 1.0 by default.
c fs(17) = xeden * Multiplier of electron density derivative
c component in Sperp term. Set to 1.0 by default.
c fs(18) = xsecder * Multiplier of second derivative
c component in Sperp term. Set to 1.0 by default.
c
c *Reference: G.M. Staebler and F.L. Hinton, Particle and Energy
c Confinement Bifurcation in Tokamaks, Phys. Fluids B 5
c (4), p. 1281-1288, April 1993.
c ----------------------------------------------------------------------
c
ifsflag = 0
aeh = 0.0
ael = 1.0
aih = 0.0
ail = 1.0
bh = 0.0
bl = 1.0
ch = 0.0
cl = 1.0
alfae = 1.0
alfai = 1.0
betah = 1.0
sigma = 1.0
gammah = 2.0
lsfctr = 1
xrot = 1.0
xeden = 1.0
xsecder = 1.0
*
c ----------------------------------------------------------------------
c --- VOLUME CONTROL PARAMETERS
c ----------------------------------------------------------------------
c
c Turn volume control OFF. (LEAVE IT OFF UNTIL FURTHER NOTICE - HSJ)
c% volaray(m) Volume (cm**3) desired at time timeqbcd(m)
c% voladj Volume adjustment parameter. A value of 0. is
c suggested initially. After a restart the last
c value calculated for voladj should be input.
c% volnudge Change in volume adjustment parameter (suggest 0.01)
c% dvadjmax ???
c% tolvol Maximum relative difference allowed between actual
c and desired volumes. A value of 0.005 is suggested
c if volume control is desired. A large value will
c% itvol Maximum number of attempts made to obtain desired
c volume (suggest 6)
c% limpos Side on which plasma touches limiter. Specify 'left'
c or 'right' if volume control is desired. A value of
c 'none' will turn volume control off.
c% itorfluse set to 1 to calculate the toroidal flux inside
c the outermost flux surface using rectangular area
c elements. This method avoids problems in integration
c near the x point but may introduce some inaccuracy
c due to the coarse [ 33 by 65] grid normally used
c for DIII-D. The usual flux surface average method
c of calculating the total toroidal flux is (essentially)
c independent of the grid size but may be somewhat inaccurate
c due to integration near the x point (the actual contour
c used is pulled away from the x point some small amount)
c For plasmas without an x point itorfluse = 0 (which is
c the default) is recommended. If there is/are x point(s)
c then itorfluse = 1 may be necessary. if in doubt try both!
c (the max value of rho is what is affected)
c
c ----------------------------------------------------------------------
c
limpos = 'none'
itorfluse = 0
do i=1,10
volaray(i) = 1.0e30
end do
tolvol = 1.0e30
voladj = 0.0
volnudge = 0.01
itvol = 6
dvadjmax = 0.0001
*
c ----------------------------------------------------------------------
c --- OUTPUT PARAMETERS
c ----------------------------------------------------------------------
c% ieqprt 1: produce copious output concerning equilibrium and
c flux surface averages
c 0: (default) produce limited output
c% iterdb 1: means create a file of ITER database quantities
c A file named "iterdb" will be created with the results at
c all times in it. Note that any previously existing "iterdb"
c file will be destroyed when the new one is created. The file
c will be written at the end of the run if MHD calcs are not
c done, so that only a single time slice is created. If MHD
c calcs are done, then the file will be written at the end of
c each transport/equilibrium cycle.
c iterdb = 1 is implemented only for codeid = dee
c 0: (default) do not create the database file
c -1: same as 1, but the file is instead named "iterdb.nc", and is
c in portable, self-describing "netCDF" format instead of ASCII
c% iterdsc 1: print out a short description of each ITER database
c variable before printing out the values
c (iterdb must be set of course)
c 0: (default) don't print the description
c iterdsc is currently not an option.
c the code sets it to 1 (to be compatible with
c read operation; see irwflag)
c% irwflag flag used to indicate if read (irwflag=0)
c or write (irwflag=1) operation is to be performed
c on the file "iterdb" (this option is
c included so other codes can easily read the file
c using the same subroutine that originally wrote
c the file. At present there is no way to use
c the information read from file "iterdb"
c in ONETWO so this flag should be left at 0
c -which is the default value- ).
c% j2prt Number of flux surface intervals between output
c% iprtit 1: print output at every iteration in the
c equilibrium/transport cycle
c 0: print output only for the last iteration in each
c equilibrium/transport cycle
c% ieqdsk 1: write eqdsk files at times given by timeqbcd and deqlst
c (this could lead to a large number of files)
c 0: do not write eqdsk file
c -1: generate only one eqdsk file at the final time
c (or at the last time a successful step was taken)
c% comp_methd_eqdsk 0: use standard method to generate p and f (for print out
c to eqdsk only)
c 1: use modified method SEE SUBROUTINE WRTEQDSK FOR
c DEFINITIONS OF THE TWO THEORETICALLY EQUIVALENT,
c BUT COMPUTATIONALLY DIFFERENT, METHODS
c% rf_output 0: discard standard output of TORAY and GAFIT
c 1: allow standard output of TORAY and GAFIT to flow <- DEFAULT
c 2: redirect standard output of TORAY and GAFIT to
c "toray.log" and "gafit.log" respectively ( overwriting)
c 3: redirect standard output of TORAY and GAFIT to
c "toray.log" and "gafit.log" respectively (concatenating)
c 4: redirect standard output of TORAY and GAFIT to
c "toray_NN.log" and "gafit_NN.log" respectively, NN = 01,02,...
c
c ----------------------------------------------------------------------
c
ieqprt = 0
j2prt = 5
iprtit = 0
ieqdsk = 0 ! changed from 1 to 0 by HSJ on 11/29/94
rf_output = 1
iterdb = 0
iterdsc = 0
irwflag = 0 ! WRITE (rather than READ) ITER database file
comp_methd_eqdsk = 0
c