1 trxplib
The `trxplib' library is used to load an `xplasma' from a TRANSP archive.
It is oriented to extraction of timeslices from the TRANSP data.  It
converts known standard TRANSP data to MKS units.

Related libraries:

  trread -- read TRANSP data (REAL precision, generally no MKS conversion).
            access to calculator & "multigraph" information.
            access to TRANSP namelist.

  xplasma-- store an MHD equilibrium (currently, axisymmetric equilibria
            only).  Return interpolated field and metric data.

Each of these libraries is an NTCC module in its own right.  These modules
have their own documentation (it won't be repeated here).

2 summary_of_calls
These routines are listed in a standard calling order:
------------------------------------------------------

trx_msgs    -- set lun and/or file for error/warning message output.

trx_connect -- connect to a TRANSP run via file or MDS+ interface.

  (now, trread routines can also be used to access the TRANSP archive).

trx_label   -- retrieve label for connected TRANSP run.

trx_tlims   -- return start and stop time of connected run.

trx_ntimes  -- return numbers of times in run (# of scalar times; # of
               profile times).

trx_time    -- set the current "time of interest" with averaging option.
               time can be reset repeatedly for connected run.

trx_gtime   -- get currently set "time of interest".

trx_init_xplasma -- re-initialize xplasma module.  should be done after
               each change of time of interest.

trx_mhd     -- acquire basic MHD equilibrium data at selected time.
               load xplasma geometry.

trx_wall    -- optional -- read TRANSP namelist for limiter / wall
               position; construct (R,Z) grid to cover plasma and
               scrape-off region up to the wall.

trx_bxtr    -- construct B field from MHD data; cartesian extrapolation 
               option if trx_wall is used.  load xplasma B field.

  (now, xplasma "eq_" routines can be used to interpolate field and
  metric quantities, compute flux surface averages, etc.)

trx_scal    -- fetch a scalar value from the TRANSP data archive.

trx_prof    -- fetch a profile from the TRANSP data archive.
               add it to xplasma.

trx_calc    -- use "RPLOT calculator" to create a new function for 
               access by trx_scal or trx_prof

trx_lintrans -- create a new profile g from an existing profile f, by
               the formula g = a*f + b, where a and b are given as
               input.

trx_nspec   -- determine the number of plasma species in the underlying
               TRANSP run

trx_spec_lbl -- obtain list of plasma species, charge, mass, and labels.
               This list can include fast ion species.

trx_spec_prof -- fetch specie densities or "temperatures" or
               <Eperp> or <Epll> profiles.

2 trx_msgs

 subroutine trx_msgs(lun,zfile)
 !
 !  redirect warning and error messages to the indicated lun & file.
 !  if the file on lun has already been opened, pass zfile = ' '
 !
   implicit NONE
 !
   integer lun           ! fortran lun for messages (input)
   character*(*) zfile   ! file for messages, or blank (input)

--------
This routine also redirects any message output from `trread' or
 `xplasma'.

2 trx_connect

 subroutine trx_connect(path,ierr)
 !
 !  ** fortran interface **
 !
 !
 !  connect to a TRANSP run at specified "path:"
 !
   character(*), intent(IN) :: path   ! path / runid
   integer, intent(OUT) :: ierr       ! completion code, 0=OK
 !
 !    MDS+ path syntax:
 !                 MDS+:<server-name>:<tree-name>(<shot-number>) or
 !                 MDS+:<server-name>:<tree-name>(<tok.yy>,<runid>)
 !
 !    Unix file path syntax:  <path>/<runid> or just <runid>
 !
 !    VMS file path syntax:   <disk>:[<dir>]<runid> or just <runid>

--------
trx_connect is always the first step in using trxplib to access
TRANSP data and load the xplasma module.  The caller should check
the status code for possible error.

Note:  for useful results:
The user needs to know (from a reliable source) which TRANSP run to
access.  The trxplib module does not have any facility for finding 
"the right" TRANSP run.

2 trx_label

 subroutine trx_label(zlbl)
 !
 !  return label for current TRANSP run
 !  if no run is connected, a blank label is returned
 !

   character*(*), intent(out) :: zlbl

--------
This routine can be successfully called after trx_connect has
been called.

2 trx_tlims
 subroutine trx_tlims(stime,ftime,ierr)
 !
 !  return the start time and stop time, seconds, of the current run.
 !  return ierr.ne.0 if no run is connected.
 !
   real*8, intent(out) :: stime      ! start time of run (seconds)
   real*8, intent(out) :: ftime      ! stop time of run (seconds)
 !
   integer ierr                      ! completion code, 0=OK

--------
This routine can be successfully called after trx_connect has
been called.

A non-zero completion code generally means, no run is connected.

2 trx_ntimes(insc,inpr,ierr)

 subroutine trx_ntimes(insc,inpr,ierr)
 !
 !   return the number of times in the run
 !
   integer, intent(OUT) :: insc     ! # of times, SCALAR f(t) data
   integer, intent(OUT) :: inpr     ! # of times, PROFILE f(x,t) data
   integer, intent(OUT) :: ierr     ! completion code: 0=OK

--------
This routine can be successfully called after trx_connect has
been called.

A non-zero completion code generally means, no run is connected.

2 trx_time
 subroutine trx_time(ztime,zdelta,iwarn,ierr)
 !
 !   set the time and delta_t for fetching/averaging data
 !
   real*8, intent(IN) :: ztime      ! time of interest, seconds
   real*8, intent(IN) :: zdelta     ! delta(time), +/- for avg, seconds
   integer, intent(OUT) :: iwarn    ! warn if ztime-zdelta or 
                                    !         ztime+zdelta
                                    ! ...is outside the limit
 
  integer, intent(OUT) :: ierr      ! completion code: 0=OK

--------

This routine sets the "time of interest".

2 trx_gtime
 subroutine trx_gtime(ztime,zdelta_t,ierr)
 !
 !   return currently set time of interest
 !   ierr is set if time of interest has not been set
 !
   real*8 ztime             ! time of interest (seconds) (returned)
   real*8 zdelta_t          ! +/- averaging time (seconds) (returned)
   integer ierr             ! completion coee, 0=OK (returned)
 !
--------

Use this routine to retrieve "time of interest" and +/- averaging
time information, as stored in a prior trx_time call.

2 trx_init_xplasma
 subroutine trx_init_xplasma(ierr)
 !
 !  re-initialize xplasma
 !  read rho-axis (TRANSP XB data) (time invariant).
 !
   integer, intent(OUT) :: ierr       ! completion code, 0=OK
 !
--------

Use this routine to clear and re-initialize a new xplasma, e.g. 
after a new "time of interest" is set with trx_time.

This routine would be used inside a loop over a sequence of times
after a single trx_connect call.  If there is no time loop, and only
a single time of interest per trx_connect call, then, this routine is
not strictly necessary (it is called by trx_connect itself).

2 trx_mhd
 subroutine trx_mhd(ntheta,ierr)
 !
 !  read TRANSP MHD equilibrium data; store in xplasma for use via
 !  interpolation routines
 !
   integer, intent(IN) :: ntheta    ! no. of theta grid pts [-pi,pi]
   integer, intent(OUT) :: ierr     ! completion code, 0=OK
 !
--------
This routine reads the TRANSP MHD equilibrium data at the selected
time of interest, and stores the data in the xplasma module.  The
data are:

  a)  an evenly spaced theta grid, ntheta points from -pi to pi.
  b)  a toroidal flux profile phi(rho), webers.  The number of 
      flux surfaces matches the number in the TRANSP run itself.
  c)  a poloidal flux profile psi(rho), webers.  The magnetic axis
      is at rho=0; by convention, phi(0)=psi(0)=0; also,
      phi(j+1).gt.phi(j) and psi(j+1).gt.psi(j).
  d)  a "g" profile:  R*Btoroidal, m*T.
  e)  a "pmhd" profile:  equilibrium pressure, Pascals.
  f)  a "q" profile (not strictly necessary since phi and psi are
      given, but, provided for convenience).
  g)  The flux surfaces themselves:  bicubic splines R(chi,rho) and
      Z(chi,rho).

2 trx_wall
 subroutine trx_wall(ilun,inumR,inumZ,ierr)
 !
 !  get the limiter locations; form (R,Z) grid to cover just 
 !  this space...
 !
   integer ilun                ! lun for file read (input)
   integer inumR,inumZ         ! cartesian overlay grid sizes (input)
   integer ierr                ! completion code:  0=OK
 !

-- or -- 
 subroutine trx_wall_RZ(ilun,Rmin,Rmax,inumR,Zmin,Zmax,inumZ,ierr)
 !
 !  read from the TRANSP namelist the TRANPS circles & lines 
 !  limiter specifications, then create this limiter in xplasma, 
 !  and create an inumR x inumZ cartesian box that covers a 
 !  rectangular region enclosing both the core plasma and the 
 !  scrape-off plasma (btw limiter and core plasma)
 !
   integer ilun                ! lun for file read (input)
   real*8 Rmin,Rmax            ! Rmin,Rmax of grid (0.0 for autogrid)
   real*8 Zmin,Zmax            ! Zmin,Zmax of grid (0.0 for autogrid)
   integer inumR,inumZ         ! cartesian overlay grid sizes (input)
   integer ierr                ! completion code:  0=OK
 !
 ! **NOTE** Rmin,Rmax,Zmin,Zmax modified on output
 !          to give the actual limits used, i.e. taking limiter 
 !          into account.  If these quantities are non-zero on
 !          input, they define a minimum range to be covered,
 !          which may be expanded to reach the limiter given
 !          in the TRANSP namelist, if necessary.
 !
--------

This routine reads the TRANSP namelist to determine the limiter or
wall position as defined for the run, and creates a cartesian grid
which covers both the core plasma and the limiters / walls all
around it.

If an extrapolation of the B field and/or plasma parameters onto
an (R,Z) grid extending beyond the core plasma is desired, this 
routine should be used.  Otherwise it should be omitted.

2 trx_bxtr
 subroutine trx_bxtr(ibccw,ipccw,irzflag,ierr)
 !
 !  form the B field bicubic splines in the core plasma, and,
 !  optionally, extrapolate field to cover (R,Z) rectangular grid.
 !
   integer, intent(in) :: ibccw ! Btoroidal:  1 for ccw, -1 for cw
   integer, intent(in) :: ipccw ! Itoroidal:  1 for ccw, -1 for cw
   integer, intent(in) :: irzflag ! =1:  extrapolate field to (R,Z) region
 !
   integer, intent(out) :: ierr ! completion code, 0=OK
 !
--------
(after trx_mhd) form a representation of the B field components, based
on the data loaded into xplasma by trx_mhd.

"ccw" means "pointed counter-clockwise looking down on machine from
             above".

If irzflag=1 is set, a prior call to trx_wall shall have defined the
(R,Z) grid to which to extrapolate.

If irzflag=1, options are available: see subtopic
3 extrapolation_options

TRANSP data defines a plasma equilibrium inside a prescribed plasma
boundary that was input to the code.  Generally, TRANSP does not have
reliable data on the field beyond the plasma boundary; nor do many
TRANSP runs have information e.g. on coil positions and currents which 
would be necessary to reconstruct this properly.

Nevertheless, it is handy at times to be able to extrapolate a field
into an [R,Z] box containing the plasma boundary, for various purposes.

The xplasma module provides two extrapolators.  Both are based on a
pseudo-flux-surface coordinate grid extrapolation beyond the plasma
boundary, with radial coordinate rho > rho_bdy, and poloidal coordinate
theta.

  a) an ad-hoc extrapolator:

     Bphi = (R*Bphi)|vacuum / R

     let Lpol(rho) = poloidal path length of extrapolated rho surface
     (numerically evaluated by xplasma).

     mod(Bpol(rho,theta)) = 
             mod(Bpol(rho_bdy,theta))*Lpol(rho_bdy)/Lpol(rho)

     direction of Bpol(rho,theta) is tangent to the rho surface, with
     sign consistent with the field at the boundary.

     a "roughly consistent" psi(rho,theta) is also provided, but...

     -> del.B = 0 is not satisfied
     -> del x B = 0 is not satisfied even where there is no current

     -> Bpol = grad(psi)/R is not satisfied.

     on the other hand, the extrapolated fields have reasonable 
     magnitudes and smooth behavior.  The behavior near the plasma
     boundary has roughly matched EFIT calculations for several test
     cases.

     ** this method is adequate for most purposes **

     ** this is the default **

  b) a del.B = 0 extrapolator.

     ** NO LONGER AVAILABLE ** (dmc Apr 2006)

     If one imposes del.B = del x B = 0 and extrapolates outward from
     the boundary numerically, on any grid, the result is exponentially 
     unstable.  The physical reason for this is that there are coils
     with strong currents out there beyond the plasma boundary, which
     determine the plasma position and shape.  An extrapolation which
     ignores this misses something crucial.  Unfortunately, the coil
     positions and currents are unknown to TRANSP.

     However, it is possible to have

     del.B = 0, del x B = 0 very near the boundary, and then 
     introduce current, del X B <> 0, to force a "smooth stabilization" 
     of the extrapolated field -- in place of the actual coils.

     This works, but the details of the form of imposition of current
     are somewhat arbitrary, and can give rough results for some 
     combination of input data.

     to request this extrapolator:
       call trx_bxtr_method_set(2)

     to restore the default extrapolator (a) above:
       call trx_bxtr_method_set(1)

     to see the current setting (1 or 2):
       call trx_bxtr_method_get(imethod)   ! imethod integer, intent(out).

     NOTE: setting imethod=2 no longer has any effect!

  c) a future option -- get the external B field from EFIT.  This would
     be a better choice in the long run ... and should be done if the
     external field is ever really needed in the TRANSP model or post-
     processors.

Edge smoothing option:  

Even when method (a) is used, the transition of field variation in
[R,Z] space can be sharp at the boundary.  This is because even if
the extrapolation itself is smooth, the extrapolated (rho,theta)
coordinate system's behavior is continuous but not differentiable
across rho_bdy.  (It is really hard to avoid singularities for all
input equilibria, if one tries to force a differentiable coordinate
system extrapolation).

As an option to ameliorate this, an edge smoothing of the B(R,Z) 
field representations can be specified:

   real*8 zsmooth  ! meters

   call trx_bxtr_smooth_set(zsmooth)  ! set smoothing (default: NONE)

   call trx_bxtr_smooth_get(zsmooth)  ! retrieve currently set value.

The data at the edge is smoothed by averaging +/-zsmooth m in both the
R and Z directions.  Data of distance d >= zsmooth from the edge is not
smoothed.  Data of distance d < zsmooth from the edge is smoothed by
averaging +/- d in both the R and Z directions.  The averaging is done
with a pyramidal weight function which tends to return numerically
smooth results.

The xplasma software imposes an upper limit of 0.05*[Rmax-Rmin] where
Rmax and Rmin are the radial grid extrema.  If zsmooth exceeds this 
value, the limit value is used instead.

2 trx_scal
subroutine trx_scal(zname,zuns_mks,zval,ierr)
 ! ...read a TRANSP scalar function value 
 !    at current time of interest; convert to MKS units. 
 ! ...arguments
 
   character(*), intent(in)  :: zname    ! name of desired scalar
   character(*), intent(out) :: zuns_mks ! mks units label of profile
   real*8, intent(out)       :: zval     ! scalar value returned
   integer, intent(out)      :: ierr     ! completion code, 0=OK
 !
 ! if no MKS conversion is possible, a warning will be issued and
 ! data will be supplied without units conversion.  However, units
 ! conversions are available for most common TRANSP scalar functions.
--------
The scalar f(t) data will be interpolated to or time averaged over
the currently set time of interest.  The name `zname' is the name
used in the TRANSP data archive (case insensitive).
2 trx_prof
 subroutine trx_prof(zname,zuns_mks,iordr,ident,ierr)
 !---------------------------------------------------------
 !  get units conversion factor, then, read & scale profile
 !
 ! ...arguments
   character(*), intent(in) :: zname     ! name of desired profile
   character(*), intent(out) :: zuns_mks ! mks units label of profile
   integer, intent(in)      :: iordr     ! desired order of fit
   integer, intent(out)     :: ident     ! xplasma object id#
   integer, intent(out)     :: ierr      ! completion code, 0=OK
--------
The profile f(x,t) data will be interpolated to or time averaged
over the currently set time of interest.  The name `zname' is the 
name used in the TRANSP data archive (case insensitive).  Only
TRANSP profile functions of flux surface or flux zone, and time,
can be accessed.

The "fit order" determines the type of interpolation vs. "rho"
to be used by xplasma:

  iordr = 2  --  spline          f'' continuous
  iordr = 1  --  Akima Hermite   f' continuous, less "ringing"
  iordr = 0  --  piecewise linear

2 trx_calc
 subroutine trx_calc(zname,zdescr,zunits,zexpr,iwarn,ierr)
 !
 !  invoke the rplot calculator -- return no data; define no
 !  xplasma object.
 !
 !  the calculator can create new named items that can be 
 !  loaded into xplasma by a subsequent trx_prof call.
 !
   implicit NONE
   INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100)
 !
 !  input:
   character*(*) zname  ! name of item to define in computation
   character*(*) zdescr ! description of item
   character*(*) zunits ! physical units
   character*(*) zexpr  ! calculator expression
 !
 !  output:
   integer iwarn,ierr   ! rpcalc warning and error flags
--------
For more information on the "rplot calculator" see the trread
documentation.  Note that effective use of the calculator requires
knowledge of the contents of the TRANSP archive, and, no automatic
MKS conversion is available.  The user is responsible for 
determining the physical units of the calculated data, as
appropriate to the formula provided to the calculator expression
parser.

trx_calc creates a new named function, known to trread.  If this is
to be added to xplasma, use trx_prof, giving the new name.  An MKS
units warning may occur.

2 trx_lintrans
subroutine trx_lintrans(id_orig,factor,offset,new_name,id_new,ierr)
 !
 !  form a new function of the form factor*[old function]+offset
 !    i.e. a linear transformation of an old function
 !
 !  at present only items of form f(rho) are supported, and the
 !  interpolation order must be at least 0 (piecewise linear).
 !
 !  these restrictions can be moved later, if there is a reason to do so...
 !  dmc 1 May 2001
 !
  implicit NONE

  integer, intent(in) :: id_orig        ! xplasma id of original function
  real*8, intent(in) :: factor,offset   ! defining the linear transformation
  character*(*), intent(in) :: new_name ! desired name of new function
 
  integer, intent(out) :: id_new        ! xplasma id of resulting new function
  integer, intent(out) :: ierr          ! completion code, 0=OK

 !-----------------------------
2 trx_nspec

(note that the TRREAD routine rd_nspecies can be used directly, if
separate information is desired on:

  n_thi -- number of "non-impurity" thermal ion species
  n_thx -- number of "impurity" thermal ion species
  n_bi  -- number of beam ion species
  n_rfi -- number of rf tail (minority) ion species
  n_fusi -- number of fusion product ion species.

This routine just returns the total number, which includes a slot for
the electrons.

subroutine trx_nspec(n_species)
 !
 !  return the number of plasma species (including electrons) present
 !  in the currently open run.
 !
 !  if no run is open, return n_species=0.
 !
  use trx_module
 !
  implicit NONE
 !
  integer, intent(out) :: n_species
 !
 !------------------------------------------------
 !
  integer n_thi,n_thx,n_bi,n_rfi,n_fusi
  integer lunzer
 !
 !------------------------------------------------
 !
  if(nsurf.eq.0) then
     write(lunzer(0),*) '?trx_nspec:  use trx_connect to open a run, first!'
     n_species=0
  else
     call rd_nspecies(n_species,n_thi,n_thx,n_bi,n_rfi,n_fusi)
  endif
end subroutine trx_nspec

2 trx_spec_lbl

 ** it is recommended to use trx_nspecies to get the list size first,
 then allocate list arrays of the correct size to hold the output of
 this routine.  However, the argument "maxn" supports a call from a
 statically allocated f77-style program as well.

 ** note that the Z & A values returned are REAL*8 -- not integer.
 A TRANSP model impurity can have non-integer values of Z and A.

subroutine trx_spec_lbl(maxn,n_species,slbl,itype,Zr8,Ar8,IZc)
 !
 !  return the number of plasma species, and, 
 !  return a character*20 label for each plasma species present in the run, 
 !  and, return a (real*8) Z and A for each specie and IZc atomic number
 !
 !  for some model species, Z and A vary in time; therefore the time must
 !  be selected before this routine is called.
 !
  use trx_module
 !
  implicit NONE
  integer, intent(in) :: maxn   ! max no. of species (dimension of slbl)
  integer, intent(out) :: n_species   ! actual no. of species found
  character*20, intent(out) :: slbl(maxn)   ! the species labels
  integer, intent(out) :: itype(maxn) ! type code of species (see below)
  real*8, intent(out) :: Zr8(maxn)    ! Z of species
  real*8, intent(out) :: Ar8(maxn)    ! A of species
  integer,intent(out) :: IZc(maxn)    ! atomic number of species
 !
 ! return n_species=0 if an error occurs
 !
 ! itype codes:  for j'th specie:
 !
 !   itype(j)=-1  -- electrons
 !
 !   itype(j)=+1  -- non-impurity thermal specie, usually H or He isotope
 !                   can be Li
 !   itype(j)=+2  -- impurity thermal specie:  Z and A are known, constant
 !
 !   itype(j)=+3  -- impurity thermal specie:  Z and A are functions of time
 !                   "model impurity" -- could represent a hybrid; non-integer
 !                   Z and A values possible.
 !
 !   itype(j)=+4  -- beam ion specie
 !   itype(j)=+5  -- rf tail ion specie
 !   itype(j)=+6  -- fusion product ion specie
 !
 !---------------------------------------

2 trx_spec_prof

 ** it is recommended to use trx_nspecies to get the list size first,
 then allocate list arrays of the correct size to hold the output of
 this routine.  However, the argument "maxn" supports a call from a
 statically allocated f77-style program as well.

 ** note that for fast ion species "T" really means 2/3 <E>, where <E>
 is the average energy per particle; <E> = <Eperp> + <Epll> can be 
 verified.

 ** the units of the profile splines are: 

    densities:  #/m**3
    T's and E's:  KeV

 use xplasma routines to evaluate the profiles and/or their derivatives
 along given position vectors.

subroutine trx_spec_prof(maxn,n_species,iorder,zident,id_array,ierr)
 !
 !  set up the indicated profiles:
 !    zident.eq.'N'  --  densities, for all plasma species
 !    zident.eq.'T'  --  temperatures, for all plasma species
 !    zident.eq.'<EPERP>' -- average perpendicular energy / particle, all...
 !    zident.eq.'<EPLL>' -- average parallel energy / particle, all...
 !
 !  it is allowed for some or all of the profiles to be already defined.
 !
  implicit NONE
 !
  integer, intent(in) :: maxn           ! max no. of species expected
  integer, intent(out) :: n_species     ! actual no. of species found
  integer, intent(in) :: iorder         ! fit order (e.g. 1 for Hermite)
  character*(*), intent(in) :: zident   ! "N" or "T" or "<EPERP>" or "<EPLL>"
 !    (case insensitive test is used)
 !
  integer, intent(out) :: id_array(maxn)  ! ids of indicated profiles
  integer, intent(out) :: ierr          ! completion code, 0=OK
 !
 !----------------------------------


2 trxpl_test_program

Along with the trxplib program comes a test driver, `trxpl',
which is actually useful in its own right (for example, it
can convert a TRANSP equilibrium timeslice into G-EQDSK format).

The `trxpl' source can be used as an example of programming 
with the trxplib library.