Documentation of the FRANTIC NTCC Module (dmc: draft, Nov. 2001). ================================================================= By D. McCune (transp_support@pppl.gov) 27 Nov 2002. reviewed and edited DMC Jan 2009. Abstract: --------- The FRANTIC model performs a very fast 1d neutral gas transport calculation for tokamak core plasmas, taking into account charge exchange and impact ionization atomic reactions in a simplified nested cylindrical flux surface geometry. This model is a generalization of the original code by Steve Tamor, S. Tamor, J. Comput. Phys. 40 (1981) 104. with extensions to permit multiple ion and neutral species including Helium as well as Hydrogen isotopes. The calculation takes as input a neutral source, which can be an "edge" influx neutrals representing recycling or gas flow, or, a "volume source" due to recombination and/or the charge exchange halo around injected neutral beams. The main calculation is called once for each neutral source. All sources are assumed to be distributed evenly around the circularly shaped flux surfaces. The model also accepts as input data profiles of the rates of depletion of neutrals by collisions with fast ion species-- a dominant process in low density intensely heated tokamak experiments with signifant fast ion density fractions. FRANTIC then calculates neutral density and "temperature" and toroidal velocity profiles for each of the neutral species associated with the given neutral source. Particle, momentum and energy sources/sinks are calculated, as well as complete particle, momentum and energy balance information for each neutral specie. Atomic physics data comes from the PREACT module, which is an available component of the NTCC library, http://w3.pppl.gov/ntcc. Introduction: ------------- This document describes the FRANTIC NTCC module. This is a fortran-90 conversion of the traditional FRANTIC model that has been used in TRANSP for many years, and which was successfully benchmarked against older 1d Monte Carlo packages (Aurora models) from TRANSP and BALDUR in the early 1980s. The model gives a "qaulitatively correct" description of the transport effects of thermal neutral atom sources; correctness in detail would require a more accurate 2d or 3d geometry and a much more time consuming calculation. The fortran-90 FRANTIC allows the user to specify any number of ion and neutral species as well as the radial resolution of the model. Generally, Akima-Hermite or cubic spline interpolation is used for mapping profiles between the grid used by FRANTIC and the grid in the driver code; the interpolation software is based on the PSPLINE NTCC module (http://w3.pppl.gov/ntcc). The PSPLINE "ezspline" interface is used: "spline objects" for input and output profiles are passed down to FRANTIC from the driver program. In order to acquire the necessary definitions, the FRANTIC driver routine will need to contain the lines: use ezspline_obj use ezspline ... type(ezspline1_r8) :: ,... ! 1d spline objects. More information on the 1d profile splines is given in a separate section, below. The fortran-90 FRANTIC interface and code is implemented in 100% REAL*8 precision-- all floating point library subroutine arguments are in this precision. The units convention for the interface is MKS, except: KeV for temperatures. The input data consists of the list of plasma species, the desired resolution of the FRANTIC calculation, the driver code's grid for FRANTIC output, plasma geometry information, the temperatures, densities and rotation velocities of the ion species, and the fast ion depletion rates for each neutral specie. Neutral sources are specified as an influx rate with an average energy (edge sources) or as source density profiles (volume sources). All output associated with a neutral source scales linearly with the magnitude of the source; therefore, it is common practice to compute the "neutral transport response" for unitary sources for each source type; these can then be recombined later, algebraically, using the actual source magnitudes, with sums being taken as desired to compute total specie neutral densities, charge exchange power and momentum loss rates, etc. Such a renormalization method also allows the neutral gas model to be called less frequently, as renormalization rather than full recalculation is sufficient for many timesteps. However, the renormalization of the collected FRANTIC model results (and the decision when to redo the full calculation) is the responsibility of the driver code. The interface description is broken into parts as follows: a) input data specification (plasma geometry, species and parameters). b) input data specification (individual neutral sources). c) neutral transport calculation with debugging options. d) retrieval of output information associated with individual sources). The recommended general structure of the driver code is: [input data specifications (a)] do [loop over neutral sources] [source-specific input data specifications (b)] [neutral transport calculation (c)] [retrieve source-specific output information (d)]. enddo [renormalizations & sums as needed by driver code]. In order to use and test the FRANTIC NTCC modules, the following additional NTCC modules, upon which FRANTIC is dependent, are required: PREACT -- atomic physics database -- ionization and charge exchange rate coefficients for neutrals in a target plasma. R8SLATEC -- mathematical software used in some PREACT fits. PSPLINE -- spline interpolation software used for 1d profiles. EZCDF -- NetCDF portable binary file i/o package (for testing). PORTLIB -- portable access to operating system features. Because FRANTIC is implemented in fortran, it needs a LUN or logical unit number for output messages. This defaults to "6" which is mapped on all operating systems to standard output. To reroute FRANTIC output to an open file on a different LUN, make the call: integer ilun ! LUN of open file for FRANTIC messages call frantic_lun_set(ilun) The remaining sections give detailed information on the FRANTIC interface requirements. The final section gives a summary of the interface which may serve as a quick reference guide. A final note on the atomic physics framework of FRANTIC is in order. This model is intended for use in hot tokamak core plasma conditions only, where all light ion species are assumed to be fully stripped. For FRANTIC purposes, "charge exchange" refers only to atomic reactions taking the form: reagents: 1 neutral atom, 1 ion products: 1 ion, 1 neutral atom. Any reaction which ionizes the neutral atom but does not neutralize the ion (even if it is charge exchange in the sense of creating a not-fully-stripped light ion) is considered "impact ionization". The physical rationale for this is as follows. The importance of neutral transport comes from the fact that neutral atoms are not bound by the tokamak's confining magnetic field. Not-fully-stripped light ions are confined by the field and are fully stripped by ion- ion or ion-electron collisions in short order, in a high temperature tokamak plasma-- so the net effect is the same as if the initial collision had resulted in an impact ionization with a free electron source. Therefore, reactions taking the form: reagents: 1 neutral atom, 1 fully stripped ion products: 1 ion, 1 not-fully stripped ion i.e. "non-neutralizing charge exchange" reactions, are, for FRANTIC purposes, lumped together with classical impact ionization reactions: reagents: 1 neutral atom, 1 ion products: 2 ions, free electron(s) from the neutral atom. These considerations are reflected in the atomic physics database (PREACT) used by FRANTIC. More on 1d profile spline objects. ---------------------------------- The FRANTIC driver provides input data in the form of 1d profile splines (ezspline objects) defined angains a normalized flux coordinate label. The user may select any normalized flux surface label "x" provided it satisfies: x = 0 at the magnetic axis and x = 1 at the plasma boundary. Internally, FRANTIC will compute on bins that are equispaced in "x". When setting up 1d profile spline objects, the data needs to be defined on the x grid. The driver's x grid need not be evenly spaced but it must cover the entire range [0,1]. This often means, in practice e.g. for "zone centered" data, that 1/2 zone width extrapolations are needed to provide data at the end points. The "ezspline" objects in the NTCC PSPLINE module support two types of splines for representation and interpolation of 1d profiles: Cubic splines -- continuous and twice differentiable but subject to "ringing" in presence of noisy data. Akima Hermite splines -- continuous and once differentiable, much less prone to "ringing". 2nd order differentiability is not important to the FRANTIC model; therefore Akima Hermite splines are recommended for most profiles. Cubic splines should only be used for profiles which are known to be "smooth". Splines have boundary conditions at x=0 and x=1. These can usually be defaulted, in which case a boundary condition is chosen automatically that is "natural" for the data given. For some profiles f(x) it may be convenient or desirable to specify df/dx = 0 at x=0. To set up a spline to provide input data for FRANTIC, the following coding elements are needed (in a fortran-90 driver): use ezspline_obj ! import NTCC PSPLINE "ezspline" definitions use ezspline type(ezspline1_r8) :: f ! declare a 1d spline object for f(x). integer npts ! no. of pts in data used to define the spline. integer bcs(2) ! boundary condition options integer iok ! error check real*8, dimension(:), allocatable :: xdata, fdata ... npts = [...number of data points to cover x range [0,1]] ! in many models this would be the number of zones + 2 ! to allow for 1/2 zone width end point extrapolations. allocate(xdata(npts),fdata(npts)) ... set up xdata(:) and fdata(:) from input data ... ! boundary conditions: bcs = (/0, 0/) ! default boundary conditions at x=0 and x=1 ! or... bcs = (/1, 0/) ! specify df/dx at x=0, default at x=1 ! or... bcs = (/1, 1/) ! specify df/dx at x=0 and at x=1 call ezspline_init(f, npts, bcs(1:2), iok) if(iok.ne.0) [...handle error...] ! if non-default bcs are set f%bcval1min = [value for df/dx at x=0 (if needed).] f%bcval1max = [value for df/dx at x=1 (if needed).] f%isHermite = 1 ! for Akima Hermite. **Note** default is cubic spline. f%x1(1:npts) = xdata(1:npts) ! covering the range [0,1], this is the ! the x grid upon which the spline input data is defined. ! here the Spline coefficients are calculated: call ezspline_setup(f, fdata(1:npts), iok) if(iok.ne.0) [...handle error...] deallocate (xdata,fdata) ! f is still available, and can be ! passed to FRANTIC subroutines. ...[pass f spline to FRANTIC]... call ezspline_free(f, iok) if (iok.ne.0) [...handle error...] In most cases there will be a standard set of operations to map data from a driver code's grid into spline objects for passage on into FRANTIC. A subroutine should be written to encapsulate this standard set of operations conveniently. Input Data Specification (Plasma Geometry, Species, and Parameters). -------------------------------------------------------------------- This section gives a detailed description of how to set up input data for FRANTIC. A compressed reference description is given in the final section of this document. The first task is to assemble the plasma species list for FRANTIC. The plasma species list will be specified to FRANTIC, then the geometry, and then the profiles (densities, temperatures, etc., for each of the plasma species). Plasma species to be considered are: electrons, bulk plasma ions, impurity ions, and, fast ion species. The data needed for each plasma specie is described in what follows. Electrons: the density, temperature, and toroidal velocity profiles are needed. Bulk plasma species: Each isotope of Hydrogen, Helium, and (optionally) Lithium is to be counted as a separate ion specie, for which the following data will be required: scalar data: atomic number (Z) and atomic mass (A, AMU). profile data: temperature, density, and toroidal velocity. Impurity species: In terms of FRANTIC definitions, impurity ion species are all ions with atomic number Z > 3. Lithium, Z = 3, can be counted as either an impurity or as a separate non-impurity bulk plasma specie. Impurity ion species are important as they participate in collisions which efficiently ionize thermal neutral gas populations in tokamak plasmas. Internally, FRANTIC uses exactly one "fully stripped light impurity", which may, however, have zero density. If the driver code also has only one impurity, then, it need only provide to FRANTIC the Z and A of the impurity, and the impurity ion's density, temperature and toroidal velocity profiles. Although Z and A are usually thought of as scalars, if the model "fully stripped light impurity" actually represents a composite of multiple impurity species, a separate call allows Z and A to be defined as profiles. In the case where multiple impurities are present in the driver, These can also be fed individually to FRANTIC, which will then constructs its internal "fully stripped light impurity" using the formulae: nx = sum[all impurities i]{ n(i) } (density profile) Ax = 1/nx * sum[all impurities i]{ n(i)*A(i) } (average A profile) Zx = 1/nx * sum[all impurities i]{ n(i)*Z(i) } (average Z profile) Tx = 1/nx * sum[all impurities i]{ n(i)*T(i) } (average T profile) Vphix = 1/nx * sum[all impurities i]{ n(i)*Vphi(i) } (average Vphi profile) all of which are of course 1d profiles. Fast ion species: FRANTIC does not need detailed information at this point. For now, just the number of fast species and each of their atomic numbers "Z" need to be known. Later, when neutral-source-specific information is given, ionization and charge exchange rate data from the fast species will be needed. Before the data profiles can be specified, the species list must be assembled, as shown in the following pseudo-code: integer :: nions ! number of bulk plasma ion species integer :: nimp ! number of impurity ion species (can be zero) integer :: nfast ! number of fast ion species real*8, dimension(:), allocatable :: zplas,aplas,zfast [...get correct values of nions nimp and nfast...] allocate(zplas(nions+nimp),aplas(nions+nimp),zfast(nfast)) zplas(1:nions) = [Z of bulk ion species (atomic no.)] aplas(1:nions) = [A of bulk ion species (AMU)] zplas(nions+1:nions+nimp) = [Z of impurity species, (atomic no.)] aplas(nions+1:nions+nimp) = [A of impurity species, (AMU)] zfast(1:nfast) = [Z of fast ion species (atomic no.)] ! all arguments are input call FRANTIC_species_za(nions,nimp, & zplas(1:nions+nimp),aplas(1:nions+nimp), & nfast,zfast(1:nfast)) Note: if there is a single "composite" impurity and Z and A are to be made profiles, set nimp=1 and zplas(nions+1)=6 and aplas(nions+1)=12. There will be an opportunity to define the Z and A profile variation later. --------- Next, the resolution of the FRANTIC calculation must be chosen. The running time of FRANTIC goes up as n**2, where n is the number of plasma zones in the FRANTIC grid. Because of the approximate nature of the FRANTIC geometry, and the surface-averaged treatment of neutral sources, it would probably not be appropriate to go beyond n=50, at the most, in real applications. Even at n=50, FRANTIC is unlikely to carry a significant fraction of the computational load of an integrated tokamak plasma model. But, the driver code must make the choice. integer :: nzones_frantic ! number of FRANTIC plasma zones nzones_frantic = [the desired number of zones] ! all arguments are input call FRANTIC_profile_alloc(nzones_frantic) ! **** Note that this call must be made *after* the frantic_species_za call. It allocates all internal profiles needed by the FRANTIC calculation. Note also: interpolation effects can be avoided if the FRANTIC grid matches the grid of the main calculation. However, FRANTIC cpu time goes up as [grid size]**2. Recommendation: if the main calculation grid size is .le. 50 zones, have FRANTIC match the main calculation grid; otherwise set nzones_frantic to (at least) 50. --------- The next step is to specify the FRANTIC geometry, and to specify the driver code's grid, to be used for FRANTIC output profile subroutines. The following information is needed: integer :: ispline ! spline/Hermite option flag for Vol(x) and Rmaj(x) integer :: nbdys ! number of zone bdys in driver "x" [0,1] grid real*8 :: xbdys(nbdys) ! zone bdy values -- strict ascending order. ! xbdys(1)=0, xbdys(nbdys)=1 expected. real*8 :: vols(nbdys) ! volumes within each bdy -- m**3. ! vols(1)=0 expected. real*8 :: Rmaj(nbdys) ! major radii of centers of regions enclosed ! by each flux surface bdy. [...set up input arguments...] Note that xbdys(...) need not be evenly spaced. The number of plasma zones is nzones = nbdys-1. FRANTIC output routines will return (nzones) numbers representing (for example) a neutral density profile, as mapped to the driver code's grid. Set "ispline" as follows. If the driver's "x" coordinate is (r/a) or proportional to sqrt(magnetic flux), set ispline=1. If the driver's "x" coordinate is normalized magnetic flux (no sqrt), set ispline=0. When all is ready: call FRANTIC_radius(ispline, nbdys, & xbdy(1:nbdys), vols(1:nbdys), Rmaj(1:nbdys)) Internally FRANTIC treats the universe as a set of nested cylinders. The radius of the cylinder corresponding to flux surface x will be r(x) = sqrt( Vol(x)/(pi*2*pi*Rmaj(x)) ) Where Vol(x) and Rmaj(x) are splines constructed from the driver data provided in this call. --------- The final step in this section is to provide temperature, density, and toroidal velocity profiles for each thermal plasma species, and, if the model impurity represents a conglomerate, the impurity Z(x) and A(x) profiles as well. This is done by a series of calls of the form: ! all arguments are input call FRANTIC_profile_set(indx,n(x),T(x),Vphi(x)) with ezspline objects carrying the profile data, and where the integer indx specifies the specie: indx=0 for electrons indx=1:nions ... for thermal plasma ion species indx=nions+1:nions+nimp ... for thermal impurity ion specie where "nions" and "nimp" are as given in the preceding FRANTIC_species_za call. In terms of the arguments of the preceding FRANTIC_species_za call, if indx > 0, the profiles passed are for the ion species with atomic number Z=Zplas(indx) and atomic mass A=Aplas(indx). Note that profiles must be provided exactly once for *all* plasma species specified in FRANTIC_species_za; the main FRANTIC calculation will report an error if any are omitted or duplicated. If there is only one impurity specie, and its Z and A are to be set as profiles, this is done with the call: ! all arguments except IERR are input call FRANTIC_impurity_ZAprof(Z(x),A(x),IERR) with spline arguments carrying the Z and A profiles, and an output completion code, integer IERR, equal to zero in case of normal completion. If IERR is set to a non-zero value, this means the relation A = approximately 2*Z is not satisfied, indicating a probable error in the passed profile data. Input Data Specification (individual neutral sources). ------------------------------------------------------ Once the basic profile data and geometry have been specified, the FRANTIC calculation is usually called in a loop over neutral sources. This section covers input data that should be specified which is specific to each neutral source. ------- Neutral specie selection: The first step is to identify the neutral source specie of the upcoming source. However, it should be borne in mind that each neutral source causes N neutral populations to be created, because of charge exchange. Thus, an H neutral source in an H,D,T plasma will result in the presence H, D, and T neutral populations. On the other hand, an H neutral source in an He++ plasma will not create neutral He, since H atoms do not carry the requisite two electrons. In defining the neutral source specie, it will be convenient to refer to the Zplas and Aplas arrays used in the preceding FRANTIC_species_za call. Let ig be in the range (1:nions), indexing a non-impurity ion specie, the specie of the neutral source under consideration. Also, let there be an integer array i0map(1:nions) available. Then, call FRANTIC_neutral_specie(Zplas(ig),Aplas(ig),I0MAP(1:nions)) I-zero-map, I0MAP(1:nions) is returned: for igi in the range (1:nions): I0MAP(igi)=0 if there are no charge exchange neutrals for plasma specie igi arising from this source I0MAP(igi)=k if charge exchange does create a neutral source for this specie-- neutral species #k. I0MAP is important, because it gives the mechanism for referring to specific neutral populations, when specifying further data or retrieving FRANTIC calculation outputs later on. ------- Neutral depletion (or "sink") rates due to presence of fast ion species: In many tokamak experiments, the presence of fast ions has a dominant effect on core plasma neutral transport. For each neutral specie "k" created by the current neutral source, the following rate profiles should be provided: snbii(x) = n*sigma*v (1/sec rate) of impact ionization of specie "k" neutrals by all fast ion species (sum over all fast species). and also the charge exchange rate due to each individual fast ion specie "j", "j" ranging over 1:nfast: snbcx[j](x) = n*sigma*v (1/sec rate) of neutralizing charge exchange of neutral specie "k" with fast ion specie "j". The charge exchange rates have to be separated by fast species, in order to allow a correct calculation of the free electron source, given that the Z's of the fast ion species and the neutral species need not be the same. The code to transfer the requisite rate profiles will look something like this: do i=1,nions k=i0map(i) if(k.gt.0) then [form snbii(x) spline for species i/k neutrals] call FRANTIC_bsink_ii(k,snbii(x)) do j=1,nsfast [form snbcx(x) spline for species i/k neutrals by fast ions "j"] call FRANTIC_bsink_cx(k,j,snbcx(x)) enddo endif enddo ------- Neutral source parameters: If this source is an "edge" source, coming in through the plasma boundary: zflux = [source magnitude, in ptcls/m**2/sec entering through the plasma boundary, or the equivalent for volume source] zt0 = [edge source "temperature", KeV] call FRANTIC_edge_source(zflux,zt0) (average energy of an edge source neutral is 1.5*zt0 KeV). If this source is a "volume" source, there is an optional normalization: The total number of neutral atoms per second divided by the surface area of the plasma boundary "zflux" expresses the total source magnitude in terms of an equivalent edge source. The shape of the source profile is a spline, proportional to the source density profile in ptcls/m**3/sec would then be: Svol(x) = spline of neutral source profile inside plasma call FRANTIC_volume_source(zflux,Svol(x)) Alternatively, Svol(x) can already be normalized in units of #/m^3/sec, in which case this call should be made: call FRANTIC_volume_source(0.0_r8,Svol(x)) Here 0.0_r8, a constant of real*8 precision, indicates no renormalization. The source temperature is taken as the ion temperature of the plasma specie which is contributing this volume neutral source. (In terms of the physics the volume neutral sources come either from collisional recombination of ions and electrons to form neutral atoms-- usually insignificant in a tokamak core plasma-- or, more importantly, the "beam halo" thermal neutral source associated with neutral beam deposition). The input data specification is now complete; a neutral transport calculation can now be performed. Neutral Transport Calculation with Debugging Options. ----------------------------------------------------- The neutral transport can be calculated immediately, by means of the call: call FRANTIC_NEUTS(IERR) This will carry out the calculation and prepare output data which can then be retrieved by the driver program. The integer status code IERR is returned. The normal status code is IERR=0. IERR will generally be set, and an error message written, if input data expected by FRANTIC is missing. Some extra debugging calls are available. If it is desired to save the input data for the upcoming FRANTIC_NEUTS call, this can be done. Let filename (character string) = name of NetCDF file to create. call FRANTIC_input_save(filename, IERR) then, on output, the status code IERR is set to a non-zero value, if the NetCDF file (filename) cannot be written. This call is strictly optional, but could be used if there was a debugging problem inside FRANTIC itself which required appraisal by the authors. The NetCDF file provides input data for a standalone FRANTIC calculation, the "frantic_test" test driver that is provided with the FRANTIC NTCC module. FRANTIC_input_save should be called just before FRANTIC_NEUTS. An additional debugging call is available *after* the FRANTIC_NEUTS call. This call produces some text output containing a partial summary of the results of the neutral transport calculation: call FRANTIC_prinsum(ilun) where "ilun" is a fortran i/o unit number on an open file, or, "6" for fortran standard output. Output Information for Each Neutral Source. ------------------------------------------- After the FRANTIC neutral transport calculation is complete, the following profile information is available: o for each neutral specie involved (directly or via charge-exchange): -- n0,T0,vphi0 (neutral density m**-3, temperature KeV, and toroidal velocity m/sec) -- Si (ionization source profile, m**-3/sec). o total for this neutral source: -- Se (electron source profile, m**-3/sec) -- Scx (charge exchange reaction rate, m**-3/sec) -- Pcx (charge exchange, loss to ions, W/m**3) -- Pii (ionization, source to ions, W/m**3) -- Tqcx (charge exchange, loss to ions, Nt-m/m**3/sec) -- Tqii (ionization, source to ions, W/m**3) -- Pcxsf (charge exchange "source friction", W/m**3) -- Piisf (ionization "source friction", W/m**3) -- Seb (electrons -> charge-exchange with fast ions) -- Sflx01 -(div(electron-influx) carried by neutrals, m**-3/sec) -- Sflx02 (div(electron-outflux) carried by neutrals, m**-3/sec) -- Pflx01 -(div(power-influx) carried by neutrals, W/m**3) -- Pflx02 (div(power-outflux) carried by neutrals, W/m**3) -- Tqflx01 (div(momentum-influx) carried by neutrals, Nt-m/m**3/sec) -- Tqflx02 (div(momentum-outflux) carried by neutrals, Nt-m/m**3/sec) The "source friction" terms are associated with the difference in the average angular velocity of an ion source relative to the average angular velocity of the ion population itself. In the context of FRANTIC output, the term "momentum" refers to "angular momentum". Because d/dt[neutral density] is negligible, the above output profiles allow the following balance checks to be computed: [particle-balance] = Z*[volume-source] + Sflx01 - Se - Seb - Sflx02 (evaluated in terms of the electrons carried in by the neutral source). [power-balance] = [volume-source] + Pflx01 + Pcx -Pii -Pcxsf -Piisf -Pflx02 and [angular-momentum-balance] = [volume-source] + Tqflx01 + Tqcx -Tqii -Tqflx02 Note that (a) all of the above profiles scale linearly with the magnitude of the neutral source, and (b) [volume-source] carried electrons, power, and torque density profiles are not calculated by FRANTIC, but are readily constructed from data proved to FRANTIC as input. The quantities needed for the plasma particle, power, and momentum balance equations are: {Si, Se, Pcx, Pii, Tqcx, Tqii, Pcxsf, Piisf}. The profiles {n0, T0, vphi0} may be useful e.g. to enable a fast ion model to estimate its own losses due to charge exchange. FRANTIC provides a separate call to fetch each output profile. The fetch operation includes an interpolation or mapping to the driver code's grid, using the information provided in the "FRANTIC_radius" call. The user provides a buffer: integer :: nbdys ! no. of driver grid bdys, see FRANTIC_radius. integer :: nzones ! = nbdys - 1 real*8 :: outbuf(nzones) If the quantity is specific to a neutral species k=I0MAP(i) corresponding to ion specie i, then: call FRANTIC_get_(k,outbuf(1:nzones)) where I0MAP(1:nions) contains the indexing information for neutral species, as discussed above, cf the "FRANTIC_neutral_specie" call. For example: real*8 dens0(1:nzones,1:nions). [...after FRANTIC_neuts call...] do ig=1,nions k=I0MAP(ig) if(k.gt.0) then call FRANTIC_get_n0(k,outbuf(1:nzones)) dens0(1:nzones,ig)=dens0(1:nzones,ig)+outbuf(1:nzones) else dens0(1:nzones,ig)=0 ! no "ig" charge-exchange neutrals this time. endif enddo If the quantity is not specific to a neutral species, then simply: call FRANTIC_get_(outbuf(1:nzones)). Thus for example call FRANTIC_get_pcx(outbuf(1:nzones)) to retrieve the charge exchange power loss (from the thermal ions) due to the current neutral source, and call FRANTIC_get_pii(outbuf(1:nzones)) to retrieve the ionization (and charge exchange recapture) power profile. ------- In addition to retrieval of profile information, the following call retrieves the neutral outfluxes (across the plasma boundary) associated with the current neutral source: The outflux information can be useful when setting up renormalization schemes-- i.e. mechanisms for reusing FRANTIC outputs without doing the full neutral transport calculation, when the magnitudes of neutral sources change a little bit. real*8 :: outflux ! outflux for this neutral specie (#/m**2/sec) real*8 :: avg_energy ! average energy of escaping neutrals (KeV) call frantic_get_outflux(k,outflux,avg_energy) where integer "k" comes from I0MAP(1:nions) as already shown. Note: any output routine called with k=0 will return zeroes in its output arguments. An alternative routine that can be used for normalizations in terms of #atoms/sec (for edge sources) can be had as follows: real*8 :: influx_aps ! influx, #atoms/sec call frantic_get_inps(influx_nps) Multiplying FRANTIC edge source outputs which scale linearly (i.e. sources, fluxes, and neutral densities) by 1/(influx_nps) yields quantities "per incident atom". Note: subroutine frantic_get_inps(...) returns ZERO for a volume source. Interface summary. ------------------ In what follows: 1d ezspline objects are indicated by the notation (x), e.g. "Vol(x)". Scalar variables with names beginning with "i", "j", "k", or "n" are integers; all other arguments are real*8 floating point numbers. Array arguments are indicated by subscripting notation, e.g. "Zfast(1:nfast)". a) basic setup: call FRANTIC_lun_set(ilun) ! set LUN for FRANTIC messages nions = [number of non-impurity thermal ion species] nimp = [number of impurity thermal ion species] nfast = [number of fast species] ! give Z's and A's of ion species. call FRANTIC_species_za(nions,nimp, & zplas(1:nions+nimp),aplas(1:nions+nimp), & nfast,zfast(1:nfast)) ! in the zplas and aplas arrays, the first "nions" elements are ! for the non-impurity ions. nzones_frantic = [number of zones to use in FRANTIC calculation]. call FRANTIC_profile_alloc(nzones_frantic) ! the above call ordering is recommended. At this point, all of ! FRANTIC's internal memory has been allocated. b) define the driver's grid (to be used for FRANTIC output) and provide the data needed for FRANTIC "shifted nested cylinders" geometry. The following data are required: ispline = 1 for driver's x = (r/a) or normalized sqrt(magnetic-flux) = 0 for driver's x = normalized-magnetic-flux (no sqrt). nbdys = no. of bdys in driver's x grid xbdys(1:nbdys) = the grid, strict ascending, xbdys(1)=0, xbdys(nbdys)=1 vols(1:nbdys) = enclosed volumes, strict ascending, vols(1)=0 Rmaj(1:nbdys) = major radius of flux zond bdys (midplane center) call FRANTIC_radius(ispline, nbdys, & xbdy(1:nbdys), vols(1:nbdys), Rmaj(1:nbdys)) c) specify plasma parameters: n(x) = electron density profile T(x) = electron temperature profile Vphi(x) = electron toroidal velocity profile call FRANTIC_profile_set(0,n(x),T(x),Vphi(x)) For each ion species i (1:nions+nimp) n(x) = specie i density profile T(x) = specie i temperature profile Vphi(x) = specie i toroidal velocity profile call FRANTIC_profile_set(i,n(x),T(x),Vphi(x)) (i=1:nions -- non-impurities) (i=nions+1:nions+nimp -- impurities) If only one impurity was specified (nimp=1), and, radially varying Z and A profiles are wanted for the impurity, call FRANTIC_impurity_zaprof(Z(x),A(x),IERR) where integer IERR is set to a non-zero value if an error is detected. ------- The following steps are performed in a loop over all neutral sources; this is represented by indentation: d) specify neutral source species; get list of neutral populations that will be created (by charge exchange). (integer ig in range 1:nions -- indexing a non-impurity plasma specie): call FRANTIC_neutral_specie(Zplas(ig),Aplas(ig),I0MAP(1:nions)) (Zplas and Aplas as used in FRANTIC_species_za). I-zero-map, I0MAP(1:nions) is returned: I0MAP(igi)=0 if there are no charge exchange neutrals for plasma specie igi arising from this source I0MAP(igi)=k if charge exchange does create a neutral source for this specie-- neutral species #k. e) specify rate profiles for ionization of thermal neutrals by fast ion species, and for charge exchange of thermal neutrals with fast ion species: do i=1,ng k=i0map(i) if(k.gt.0) then ! ionization rate (1/sec) summed over ALL fast species [form snbii(x) spline for species i/k neutrals] call FRANTIC_bsink_ii(k,snbii(x)) do j=1,nsfast ! charge exchange rate (1/sec)... [form snbcx(x) spline, rate of charge exchange of species i/k neutrals with species j fast ions] call FRANTIC_bsink_cx(k,j,snbcx(x)) enddo endif enddo f) specify edge or volume source: zflux = [source magnitude, in ptcls/m**2/sec entering through the plasma boundary, or the equivalent for volume source] zt0 = [edge source "temperature", KeV] call FRANTIC_edge_source(zflux,zt0) --or-- Svol(x) = spline of neutral source profile inside plasma call FRANTIC_volume_source(zflux,Svol(x)) g) compute the neutral transport for this neutral source: call FRANTIC_NEUTS(IERR) Status code is set to zero for normal successful completion. IERR <> 0 generally indicates an error in the input data; a message will have been written. (some additional debugging calls are available, and are described in detail above). h) fetch output profiles: real*8 :: outbuf(nzones) ! nzones = nbdys-1, nbys as in FRANTIC_radius. For profiles associated with a neutral species generated by the neutral source just processed by FRANTIC_NEUTS: k=I0MAP(ig) ! map from ion specie to neutral specie index call FRANTIC_get_(k,OUTBUF) example: call FRANTIC_get_n0(k,OUTBUF) ! neutral density, m**-3 For profiles associated with the current neutral source but not a particular neutral specie: call FRANTIC_get_(OUTBUF) example: call FRANTIC_get_pcx(OUTBUF) ! charge exchange Pwr W/m**3 j) fetch output scalars: real*8 :: outflux ! outflux, #/m**2/sec, returned. real*8 :: avg_energy ! average energy/ptcl, KeV, returned. k=I0MAP(ig) ! map from ion specie to neutral specie index call FRANTIC_get_outflux(k,OUTFLUX,AVG_ENERGY) Here the loop over neutral sources ends. ------- y) free up FRANTIC memory after all calculations are completed. call FRANTIC_free [it is possible to use FRANTIC_free_profiles to delete profile arrays only, while retaining species information. This might be done e.g. to compute different sources on different grids]. z) the driver may want to compute renormalizations and/or summations of the various FRANTIC results.