Documentation of the FRANTIC NTCC Module (dmc: draft, Nov. 2001).
=================================================================
By D. McCune (dmccune@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) :: <spline_objects>,...   ! 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) :: 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_<item>(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_<item>(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
<function-name>(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_<item-name>(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_<item-name>(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.

