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).
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.
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'.
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.
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.
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.
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.
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".
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.
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).
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).
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.
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
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.
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).
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
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.
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
!-----------------------------
(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
** 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
!
!---------------------------------------
** 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
!
!----------------------------------
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.
This Document was created by hlptohtml
Written By:Manish Vachharajani(mvachhar@pppl.gov)