This is the help file for TRREAD -- a fortran library for reading
TRANSP output and namelist data.
This fortran library gives access to TRANSP output data. The TRANSP
output data can be stored in TRANSP legacy formats, in NetCDF format,
or in MDS+.
The routines in this library give access to the core functionality
of the legacy TRANSP `rplot' output data browser, that is, access to:
-- establishment of connection to a TRANSP run.
-- table of contents
-- lists of scalar functions f(t), with labels & physical units.
-- lists of profile functions f(x[j],t), with labels & units,
and associated x axis data which can also vary with time t.
-- lists of "scalar multigraphs" -- groupings of related scalar
functions all with the same physical units and each with
associated sign factor +/-1
-- lists of "profile multigraphs" -- groupings of related profile
functions all with the same phsical units, x axis, and each
with associated sign factor +/-1.
-- (multigraphs are typically used to construct particle-,
momentum-, or power-balance plots).
-- special lists-- e.g. lists of all plasma particle species.
-- the data itself
-- timebase for scalar functions f(t).
-- timebase for profile functions f(x,t).
-- the scalar functions {f(t)}.
-- the profile functions {f(x,t)} (the time evolving x axes are
also considered profile functions).
-- the calculator
-- string based, algebraic syntax.
-- standard arithmetic operators and mathematical functions.
-- numerical integration and differentiation.
-- geometry-dependent operators such as volume integrators and
flux-surface-averaged gradient operators.
-- constants, scalar functions, and categories of profile functions
as supported data types.
-- ability to create user defined scalar and profile functions.
-- ability to create and modify multigraph definitions.
-- ability to access data from multiple runs for comparison.
-- extraction of MHD equilibrium data at specified time or averaged
over a time range.
-- extraction of scalar or profile information at a specified time or
averaged over a time range.
-- extraction of TRANSP namelist data
-- limiter locations
-- TRANSP run control options
-- etc., etc.
The library routines can be used to drive time-dependent or time-slice
based post-processors, acquire data for loading of relational databases,
or to build a modern (e.g. IDL or AVS/express based) GUI-oriented data
browser tool to replace the ancient command line `rplot' data browser.
Consistent with the RPLOT legacy software on which it is built, the
trread software is optimized for the detailed examination of data
from a single TRANSP run-- i.e. one run at a time.
When connected to a particular run, data from other runs can be
fetched e.g. for comparison, by using RPLOT calculator commands.
This works reasonably well, but, the software is not optimal for
intensive cross-run comparisons.
For studying trends of select data extracted from many different
runs, the usual approach is to load a relational database. Relational
database software is better optimized for these types of queries.
TRREAD was updated ca. September 2001, to support labeling of
members of multigraphs based on a single data item from multiple
distinct runs.
The maximum length of data identifiers or "abbreviations" was
increased from 10 characters to 21 characters length.
This means that, although individual TRANSP run databases still
have a limit of 10 characters for names of data items stored on
disk, longer names can be generated internally -- specifically,
calls to `run_mg' can result in the creation of longer identifiers
with imbedded runid aliases-- for example:
NEUTT$37065K17 -- total neutrons vs. time, run 37065K17
NEUTT$37065K18 -- total neutrons vs. time, run 37065K18
TRREAD interfaces (including the calculator) now support the longer
names, but, applications built over TRREAD may need uprading to
access the new functionality; this could affect user interfaces and
the labeling of graphical output.
Also new in September 2001: The calculator command %GS2FETCH can
acquire tabulated GS2 data (cf R. Budny) as profile functions in
a TRREAD session. For more information contact dmccune@pppl.gov
or rbudny@pppl.gov.
Routines in the trread library just read the single precision
TRANSP data. Separate NTCC modules contain related software.
For example, one can use `trxplib' to extract a data timeslice from
a TRANSP run in order to setup a standalone test driver for a
physics module:
test driver
(1) acquire data
-->trxplib
-->trread (connect to TRANSP run, read data)
-->xplasma (setup interpolation on data)
-->pspline (interpolation primitives: set up coefficients)
(2) call physics module which uses the data
-->[physics module]
-->xplasma (retrieve plasma parameters, field & metric
information in coordinates chosen by the physics module)
-->pspline (interpolation primitives)
The `trxplib' library uses `trread' to extract timeslices from the data
and then uses `xplasma' to set up continuous differentiable interpolation
and mapping between cartesian/cylindric coordinates and the original
magnetic flux coordinates of the plasma model. `trxplib' also knows
most TRANSP physical units conventions and can usually do an MKS
conversion of the data. It also convers the data to REAL*8 precision.
With `trxplib' one can access TRANSP data via interpolation, without
dealing explicitly with TRANSP numerical grids.
Access via the `trread' library is at a "lower level" and does expose
the TRANSP numerical grids. On the other hand `trread' allows one
to access time evolving data, whereas `trxplib' and `xplasma' are
oriented to time slices.
The `xplasma' library presents the interface to the interpolation and
coordinate mapping software, with extensions to reconstruct plasma
parameters and magnetic field vectors and their gradients in user
selected coordinate systems. Interpolation, mapping, and field
reconstruction software is all fully vectorized. While the plasma
state is stored inside xplasma's internal fortran-90 data module,
argument vectors (e.g. locations of particles or ray bundles) are
all handled through subroutine arguments and local variables, to
yield an architecture naturally consistent with common parallelization
strategies.
The `pspline' library provides the basic interpolation software--
cubic/bicubic/tricubic splines and hermite piecewise cubic as well
as piecewise linear/bilinear/trilinear interpolation.
On most systems the `trread' library will be stored as `libtrread.a'
in an area reserved for ntcc binary object libraries. For example,
on the PPPL unix cluster these libraries are found in $NTCCHOME/lib.
In traditional TRANSP development environments the library filename is
shortened to `trread.a' and is found in $LOCAL/lib.
Programs which thoroughly exercise trread functionality will need to
link several other ntcc libraries as well. Specifically:
(Last update: clf 17 Feb 2009)
UNIX libraries (for loading) -- NOTE the order is important!
libtrread.a ! trread
libtr_getnl.a
librp_kernel.a ! "rplot kernel" library
librplot_io.a ! "rplot i/o" library
libxdatmgr.a ! rplot COMMON buffer memory management
libmdstransp.a ! rplot/TRANSP/MDS+ interface library
libinterp_sub.a ! legacy interpolation primitives
libsmlib.a ! legacy smoothing primitives
libcomput.a ! miscellaneous computational routines
libvaxonly.a ! legacy portability routines
libportlib.a ! standard portability routines
(MDS+, NetCDF -- these libraries are optional but recommended).
libsgdummy.a ! dummy routines for sglib.
libmds_dummy.a ! dummy routines for mdsplus
It is important that the libraries be specified in the order shown.
Note: if you need "Ufiles", link with:
-ltrread -ltr_getnl -lrp_kernel -lrplot_io \
-lxdatmgr -linterp_sub -lsmlib -luflib \
-lmds_sub -lufhdf -lmdstransp -lsgdummy \
-lcomput -lvaxonly -lportlib
Applications loaded with MDS+ can access TRANSP data directly over
the internet via an MDS+ client/server connection. CAUTION: the
user must arrange for permission to access server data. This is
usually administered on a tokamak-by-tokamak basis; typically the
user will need to establish an account on the server machine and
arrange for that account to have access permission to the desired
data. The MDS+ server will then create a server process in the
user's account on the server machine to mediate data requests.
MDS+ server access permission / security requirements are in
addition to any general firewall requirements set at the server
site.
Since networked access can be slow, but disk space is cheap and
TRANSP data once written never changes, trread supports a mechanism
for caching data to local disk. Users enable this by defining
the RPLOT_CACHE environment variable. For example (csh):
> setenv RPLOT_CACHE TRUE
the directory ~/.rplot_cache will be created, and cache files
will be written in subdirectories of this directory, on a run
by run basis.
> setenv RPLOT_CACHE /scratch/joe/rplot
the directory /scratch/joe/rplot/.rplot_cache will be created,
and cache files will be written in subdirectories of this
directory, on a run by run basis.
Even if cached data files are present, the cache is "validated" by
checking the TRANSP run date against a cached value. If this does
not match, the cache data is invalidated, removed, and MDS+ networked
access is employed.
Generally, if data from a given TRANSP run is to be accessed more than
once, it is worthwhile from a performance perspective to enable the
cache. Note, however, that the trread software does not actively
manage the cache in the sense of removing old data to make room for
new data. Usually, therefore, it will be best to set RPLOT_CACHE to
use a "scratch disk" which gets skimmed by the system from time to
time.
No startup calls are required. However, optional calls are available
to control message output. By default, errors and warning messages
from trread are written on fortran unit 6, i.e. stdout. Use the
following calls to change this behavior.
(optional:)
to change the LUN and/or file where messages are written:
INTEGER ilun ! FORTRAN i/o logical unit number
CALL PLC_MSGS(ilun,' ')
...directs messages to unit ILUN (file opened by caller)
CALL PLC_MSGS(ilun,'msg.file')
...directs messages to "msg.file" which will be opened
on unit (ilun) by the PLC_MSGS subroutine
There are two alternative methods:
(a) "integrated method" -- connect to the run and return standard
information such as a status code, a run label, the size of the
run's timebase, and size of largest profile function.
(b) "call sequence method" -- separate calls are made for making
the connection, for retrieving the label, and for retrieving sizes
of the run's timebase and largest profile function.
First the call sequence method is shown. Then, the call for the
integrated method is shown.
****** call sequence method ******
** to connect to a TRANSP run **
CHARACTER*(*) path ! (IN) path to run, path/runid
INTEGER ier ! (OUT) completion code, 0 = normal
CALL RCONNECT(path,ier)
The caller should check the status code to verify that the connect
was successful.
Notes on PATH syntax:
! MDS+ path syntax: two variants:
! MDS+:<server-name>:<tree-name>(<shot-number>) or
! MDS+:<server-name>:<tree-name>(<tok.yy>,<runid>)
!
! (the appropriate choice depends on the server).
!
! Unix file path syntax: <path>/<runid> or just <runid>
!
! VMS file path syntax: <disk>:[<dir>]<runid> or just <runid>
!
! examples:
! CMOD MDS+ "MDS+:CMODA.PSFC.MIT.EDU:TRANSP03(12345678)"
! OLD PPPL MDS+ "MDS+:TRANSPGRID.PPPL.GOV:TRANSP(TFTR.88,37065Y82)"
! unix, local files "$ARCDIR/TFTR/88/37065Z19"
!
** to fetch a generic label for a run, after RCONNECT call **
CHARACTER*(*) rlabel ! about 50 characters
CALL GRUNLB2(rlabel)
... the label might be useful for messages or plots.
** to fetch basic run data size statistics
It might be useful to know the largest number of words required by any
scalar or profile function. This can be found with:
INTEGER nsctime,nprtime,nxmax,nmax
CALL RPSTATS(nsctime,nprtime,nxmax,nmax)
which returns
nsctime = number of timepoints in scalar functions timebase
nprtime = number of timepoints in profile functions timebase
nxmax = largest number of X axis points occurring in run
nmax = largest total number of words for any function in
the run database (generally nmax = nprtime*nxmax).
****** integrated method ******
CALL KCONNECT(path,rlabel,nsctime,nprtime,nxmax,nmax,ier)
(individual arguments are as in the call sequence method;
PATH is input, all others are output).
TRANSP data consists of the following:
* a SCALAR FUNCTIONS timebase, and a set of named scalar functions
of time f(t). Each function has a name (up to 10 characters long,
alphanumeric, no imbedded blanks), a label (up to 32 characters
long) and a physical units label (up to 16 characters long).
* a PROFILE FUNCTIONS timebase, and a set of named profile functions
of time and the additional coordinate, f(x,t). As with the
scalars, each function has a name, label, and physical units.
In addition, each profile function has an "x axis type" attribute
(an integer) which identifies x. This x is itself a profile
function of time but generally satisfies the condition of being
monotonic increasing at all times.
note that the SCALAR FUNCTIONS and PROFILE FUNCTIONS do not
always have the same timebase. The dimensioning of profile
function data is always (nx,nt) where nx is the number of
points associated with x, and nt is the number of time points
in the profile function timebase. nx is the dimension of
contiguous storage (the first dimension in the fortran
convention).
* a set of MULTIGRAPHS, each of which consists of a label and a
"signed" list of function names. The functions in the list are
all of the same type, i.e. scalar or profile all with the same
x axis type. Also, generally all multigraph function members
share a common physical units label, which becomes the physical
units label of the multigraph. Each function in the list has
an associated sign, +/-1, which is meant to be applied prior to
display of the data. The idea here is to support "power balance"
type plots, which generally show a collection of terms related
to eachother in an additive equation, e.g. "beam heating",
"RF heating" balanced against "convection loss", "conduction
loss", and "d/dt(energy content)". Because of the way these
quantities are defined inside TRANSP, sometimes a factor of -1
needs to be applied before display. Multigraphs contain a
maximum of 15 functions.
* via the RPLOT CALCULATOR, the user can define new profile or
scalar functions, and multigraphs, and, members can be added to
or removed from multigraphs. Function and multigraph labels
can be modified. Data from other runs can be accessed.
Functionality not available through direct trread calls can
often be found by using the calculator.
* the software enforces uniqueness of all names across all lists:
if "XYZ" names a multigraph, it will not also name a scalar or
profile function.
* a general guide to the use of the RPLOT CALCULATOR is provided
ftp://ftp.pppl.gov:/pub/transp/codata/doc/rpcalc.doc
The calculator controls evaluation of arithmetic expressions based
on TRANSP data. (In the TRREAD library, the calculator is accessed
by calling `rpcalc' with a string argument containing the desired
calculator expression or command). The following examples illustrate
some calculator capabilities:
PCUR*VSUR # plasma current * surface voltage, vs. time
UETMP = VOLINT(1.5*1.601E-19*NE*TE) # integrated electron energy
# vs time and flux coordinate, stored as "UETMP"
UITMP = VOLINT(1.5*1.601E-19*NI*TI) # ditto the ions, as "UITMP".
USCALAR = %XINTERP(X=1.0,EXPR=(UETMP+UITMP))
# ion energy + electron energy, integrated out to
# the plasma boundary (X=1.0), vs. time only;
# saved as "USCALAR".
Calculator expressions can be passed as character strings to the
appropriate `trread' subroutines; results are returned in user
supplied data buffers.
The calculator can also be used to fetch data from a different run
and map it to the current run's timebase, as in:
%RFETCH("MDS+:TRANSPGRID.PPPL.GOV:TRANSP:JET.00",50632C02,TE,TE2)
which fetches the TRANSP run JET.00 50632C02 from the PPPL MDS+
server, reads that run's TE data, maps that to the current run's
database, and stores the result as TE2.
This can be the basis of a quite general method for comparing data
from multiple runs.
For detailed information on the calculator see:
ftp://ftp.pppl.gov:/pub/transp/codata/doc/rpcalc.doc
In much of the software, a data "subtype code", typically referred to
as "istype", specifies the type of a given data item:
INTEGER istype
istype = -1 -- scalar function of time
istype = +1 -- profile function of time and X axis #1
istype = +2 -- profile function of time and X axis #2
...
istype = +N -- profile function of time and X axis #N.
Generally in TRANSP:
istype=1 ... identifies f(x,t) where x is the flux surface label
and the profile data is zone-center oriented.
examples: temperatures, densities, Zeff...
istype=2 ... identifies f(x,t) where x is the flux surface label
and the profile data is zone-boundary oriented.
examples: transport coefficients, radial velocities...
the flux surface label is
x = sqrt((toroidal_flux)/(toroidal_flux_at_boundary))
other istype codes will correspond to "midplane major radius",
energy grids or chord indices for diagnostic simulations, etc.
The trread library is probably best used in conjunction with
dynamic memory allocation. When there is a data item or list
that can have variable length, there will be a corresponding
pair of routines, the first of which returns the length, and
the second of which returns the actual data or list. Between
the two calls, dynamic allocation of the necessary buffer can
take place.
For data access, another approach will be to look at the maximum
data object size returned at run connect time by KCONNECT (or
later by a call to RPSTATS). This could be used to create buffers
that would persist for the duration of connection to a particular
run, which could contain e.g. the complete data set for any
multigraph, and which could be used repeatedly.
It may be useful to know how many profile X axis types exists in the
TRANSP run, and what their names are:
** getting total number of x axis types in the run **
** getting a list of known profile x axis type names **
Summary:
Subroutine RPNUMX returns the number of defined profile X axes.
Subroutine RPXNAME returns the X axis name, given the number.
Details:
INTEGER istype,naxes,nxmax
PARAMETER (nxmax=100) ! should be big enough...
CHARACTER*21 xnames(nxmax)
CALL RPNUMX(naxes) ! get number of known x axes in run
...and assuming nxmax was big enough...
DO istype=1,naxes
CALL RPXNAME(istype,xnames(istype))
! get function name for each x axis
ENDDO
Profile function "subtype codes", generally referred to by the name
"ISTYPE", identify the function's non-temporal x axis.
** get current length of funtion or multigraph list **
** acquiring lists of function and multigraph names **
** applications bear in mind these lists can change dynamically **
Summary:
Subroutine RPNLIST gives the length of a selected list;
Subroutine RPLIST returns the actual list.
List elements are CHARACTER*21 names.
Use these routines to build a "table of contents" of the TRANSP run.
In a dynamic memory context, use RPNLIST to figure out the minimum
size of an array that can hold a particular list of names.
RPNLIST Details:
CHARACTER*(*) ztype ! type of list (input)
INTEGER istype ! subtype code (input)
INTEGER isize ! length of list (output)
CALL RPNLIST(ztype,istype,isize)
...example: CALL RPNLIST('PROFILE',0,inum)
...more details:
subroutine RPNLIST(ztype,istype,isize)
C
character*(*) ztype ! input: list type
C
C choose ztype = "scalar" or "profile" or "multigraph"
C
integer istype ! input: list subtype
C
C if ztype = "scalar", istype is ignored. for "multi" or "profile":
C
C set istype = 0 to not restrict the list by subtype.
C set istype = -1 for scalars only (e.g. if zytpe.eq.'multi')
C set istype = +N for type N profiles only
C
integer isize ! output: list size
C
C ztype can be 'profile', 'scalar', or 'multigraph'
C actually only the first 5 characters are tested with a case blind
C compare.
C
C see subroutine rplist ... to get the actual names.
C see subroutine rptype ... to get subtype data associated with names.
** EXAMPLES **
CALL RPNLIST('scalar',0,isize) ! returns total no. of scalar fcns
CALL RPNLIST('profile',0,isize) ! total no. of profile functions
CALL RPNLIST('multi',0,isize) ! total no. of multigraphs
CALL RPNLIST('profile',2,isize) ! no. of profile fcns, x axis type 2
CALL RPNLIST('multi',-1,isize) ! no. of scalar function multigraphs
CALL RPNLIST('multi',3,isize) ! no. of multigraphs, x axis type 3
C ztype can also be qualified with a "substring", using the syntax
C <list-subtype>:<substring>. If a substring is given, only names
C which contain the specified substring will be counted.
** EXAMPLE **
CALL RPNLIST('profile:NIMP_',0,isize)
! get number of profile functions whose names contain
! "NIMP_" as a substring. The substring test is based on
! a case blind comparison.
Once RPNLIST has been used to determine the length of the list...
...use RPLIST to get the actual list of names; the first two arguments
of RPLIST work the same as with RPNLIST.
RPLIST Details:
CHARACTER*(*) ztype ! type of list (input)
INTEGER istype ! subtype code (input)
INTEGER namax ! size of alist array (input)
CHARACTER*21 alist(namax) ! list returned (output)
INTEGER isize ! length of list returned (output)
INTEGER ier ! completion code (output) 0=OK
CALL RPLIST(ztype,istype,alist,namax,isize,ier)
...example: CALL RPLIST('PROFILE',0,alist,namax,isize,ier)
where plist is an array of character*21 which could have
been allocated to length namax based on prion RPNLIST call.
...more details:
subroutine rplist(ztype,istype,alist,namax,isize,ier)
C
C return sorted list of RPLOT names of requested type.
C
character*(*) ztype ! input: list type
C
C choose ztype = "scalar" or "profile" or "multigraph"
C
integer istype ! input: list subtype control
C
C if ztype = "scalar", istype is ignored. for "multi" or "profile":
C
C set istype = 0 to not restrict the list by subtype.
C set istype = -1 for scalars only (e.g. if zytpe.eq.'multi')
C set istype = +N for type N profiles only
C
integer namax ! input: size of list array
C
character*(*) alist(namax) ! output: the list
C
integer isize ! actual number of items in list
C
integer ier ! error code, 0=OK
C
C ztype can be 'profile', 'scalar', or 'multigraph'
C actually only the first 5 characters are tested with a case blind
C compare.
C
C the substring qualifier for ztype, valid for RPNLIST and described
C above, is valid for RPLIST as well.
C
C note that because of RPLOT calculator action, these lists can grow
C dynamically and may need to be refetched...
C
C see subroutine rpnlist ... to get current list size
C see subroutine rptype ... to get subtype data associated with names.
C
C error codes (no message written):
C
C ier=1 invalid "ztype" input
C ier=2 width of "alist" elements less than what is required
C (n.b. character*10 is standard)
C ier=3 length of "alist" does not accomodate actual list.
** getting information on a particular named entity **
Summary:
-- use RPLABEL to get the labeling and type information
on a particular named item in the run database. See also RPDIMS
(below).
-- use RPDIMS to get dimensioning information on a specific named
item, by using its subtype code returned from RPLABEL.
Details:
CHARACTER*21 name ! item name (input)
CHARACTER*32 label ! item label (output)
CHARACTER*16 units ! physical units (output)
INTEGER imulti ! multigraph flag (output =1 if TRUE)
INTEGER istype ! subtype code...
CALL RPLABEL(name,label,units,imulti,istype)
...example: CALL RPLABEL('TE',label,units,imulti,istype)
...more details:
subroutine rplabel(abbrev,zlabel,zunits,imulti,istype)
C
C given an RPLOT name (abbrev) return the label (zlabel) and units
C (zunits)
C
character*(*) abbrev ! input name (C*21)
character*(*) zlabel ! output label (C*32)
character*(*) zunits ! output units label (C*16)
C
integer imulti ! output =1 if a multigraph
integer istype ! output data subtype code
C
C istype = -1 -- scalar f(t) function or multigraph
C istype = 1,2,... -- various f(x,t) profile functions /
C multigraphs
C
C if the name is not recognized,
C zlabel = zunits = 'error' are returned, and
C imulti = istype = -1000 are returned.
** getting dimensioning information **
** procedure: use RPLABEL to get the function subtype code (istype);
then use RPDIMS to get dimensioning information associated with
this subtype.
INTEGER istype ! subtype code (input)...
INTEGER IRANK ! dimensionality or rank (output)
INTEGER IDIMS(8) ! dimensioning information (output)
CHARACTER*21 ZXABB ! x-axis abbreviation (output)
INTEGER IER ! completion code (output) 0=OK
CALL RPDIMS(istype,irank,idims,zxabb,ier)
...examples:
CALL RPDIMS(-1,irank,idims,zxabb,ier)
...returns dimensioning information on a scalar function of time:
irank = 1, idims(1)=[# of points in scalar functions timebase],
zxabb=' ' (no x axis), ier=0
CALL RPDIMS(1,irank,idims,zxabb,ier)
...returns dimensioning information on any function having istype=1:
irank = 2, idims(1)= [no. of points in x axis],
idims(2)= [no. of points in profile functions timebase],
zxabb='X' [or some other X axis function id], ier=0
...more details:
subroutine rpdims(istype,irank,idims,zxabb,ier)
C
C return dimensioning information for given function subtype code
C
C subroutine "rplabel" returns subtype codes given the function name.
C
integer istype ! input: function subtype code
C
C istype = -1 -- scalar f(t) function or multigraph
C istype = 1,2,... -- various f(x,t) profile functions /
C multigraphs
C
integer irank ! output: rank or dimensionality
integer idims(*) ! output: sizes of dimensions
C
C idims(1) = 1st (contiguous storage) dimension
C idims(2) = 2nd dimension -- if irank.ge.2
C ...etc...
C
C idims(j) is not referenced if j.gt.irank.
C nothing above irank=2 exists in TRANSP databases as of 1999
C (dmc).
C
character*(*) zxabb ! for profiles: X axis function
C
C the non-temporal coordinate is represented by the named function
C this is a default specification; other X axis functions can be
C chosen.
C
integer ier ! exit code, 0=OK, 1= invalid istype.
------------------------------
** getting contents of an individual multigraph **
Summary: return the labeling information and contents of a
named multigraph.
CHARACTER*21 name ! name of multigraph (mg) (input)
INTEGER istype ! subtype code of mg contents (output)
CHARACTER*32 label ! multigraph label (output)
CHARACTER*16 units ! multigraph units (output)
INTEGER infuns ! no. of functions in mg (output)
INTEGER isigns(15) ! member sign factors +/-1 (output)
CHARACTER*21 zmembs(15) ! names of member functions (output)
INTEGER ier ! completion code, 0=OK (output)
CALL RPMULTI(name,istype,label,units,infuns,isigns,zmembs,ier)
...example: CALL RPMULTI('PTEMP',istype,...)
...more details:
subroutine rpmulti(zname,istype,zlbl,zuns,infuns,isigns,zmembs,
> ier)
C
C fetch information on a named multigraph
C
character*(*) zname ! multigraph name, input
C
integer istype ! x axis type code, output
character*(*) zlbl ! multigraph label (C*32), output
character*(*) zuns ! multigraph units (C*16), output
integer infuns ! number of function members, output
integer isigns(*) ! sign code for each member, output
character*(*) zmembs(*) ! name of each member, output
C
integer ier ! completion code, output, 0=OK
C
C abnormal completion means: multigraph name invalid.
C
...notes: isigns(...) is an ARRAY of integers
zmembs(...) is an ARRAY of CHARACTER*21
...since time immemorial the array size for multigraph
member lists has been... 15.
** ROUTINES FOR FETCHING ACTUAL DATA **
** applications may want to use the "table of contents" routines
RPLABEL, RPDIMS et al, to get dimensioning information first,
thus allowing dynamic allocation of properly sized buffer arrays
for receiving the data itself.
Summary:
RPTIME_S, RPTIME_P, RPSCALAR and RPROFILE return actual
numerical data associated with specific items in the run database.
RPCALC returns data formed by evaluating an arithmetic expression
which uses TRANSP data items as operands.
(for individual functions, not multigraphs)
CHARACTER*21 zname ! name of desired function (input)
INTEGER ibufsize ! buffer size (input)
REAL zbuf(ibufsize) ! array to contain data returned
INTEGER iret ! actual number of data points returned.
INTEGER ier ! completion code returned, 0=normal
CALL RPTIME_S(zbuf,ibufsize, iret) ! scalar timebase
CALL RPTIME_P(zbuf,ibufsize, iret) ! profile timebase
CALL RPSCALAR(zname, zbuf,ibufsize, iret, ier) ! f(t) data
CALL RPROFILE(zname, zbuf,ibufsize, iret, ier) ! f(x,t) data
Summary: to calculate an item and return the result in a memory
buffer, use RPCALC.
to calculate an item and not return it, e.g. as an
intermediate step in a longer calculation, use RPCAL0.
CHARACTER*512 expr ! calculator expression string (input)
INTEGER ibufsize ! buffer size (input)
REAL zbuf(ibufsize) ! array to contain data returned
INTEGER iret ! actual number of data points returned.
INTEGER ier ! completion code returned, 0=normal
! more output:
INTEGER istype ! x axis type of result of expression
INTEGER iwarn ! non-zero if arithmetic errors were detected
...
CALL RPCAL0(expr, iwarn, ier)
CALL RPCALC(expr, zbuf,ibufsize, iret, istype, iwarn, ier)
example:
CALL RPCAL0('TMP=NE*TE',iwarn,ier)
(check for errors)
CALL RPCALC('NI*TI + TMP', zbuf,ibufsize, iret, istype,
> iwarn, ier)
(check for errors)
...more details
subroutine RPCAL0 -- expr, iwarn, ier as in RPCALC.
subroutine RPCALC(zexpr,zbuf,ibufsize,iret,istype,iwarn,ier)
C
C fetch an rplot calculator result.
C data only, use rptime_p or rptime_s for timebase.
C use rpdims for dimensioning information.
C
character*(*) zexpr ! calculator expression (input)
integer ibufsize ! size of buffer for data (input)
C
C output arguments:
C
real zbuf(ibufsize) ! buffer into which to write result
C
integer iret ! number of data points written
integer istype ! data type code of result
C
C istype=-1 -- scalar f(t)
C istype= 1 -- f(x,t), x is TRANSP zone ctrs if this is TRANSP data
C istype= 2 -- f(x,t), x is TRANSP zone bdys if this is TRANSP data
C ..etc..
C istype = 0 indicates an error
C
integer iwarn ! arithmetic warning, 0 = OK
integer ier ! completion code, 0 = OK
C
C ier=1 means: error in the calculator, messages were written.
C ier=2 means: insufficient buffer space
C
C iwarn.ne.0 means: arithmetic errors were trapped during evaluation.
C example: expr is "tmp=1/pbe" and pbe (beam heating) is zero at
C some times during the experiment. Then iwarn=1 is set; the data
C points for which an error occurred are zeroed. The caller does
C not have access to detailed information about which points had
C the error.
C
C a side effect of errors:
C iret=0 if an error occurs.
C
Summary:
First use RPMULTI to get the multigraph member names (mgnames), integer
sign factors (mgsigns, +/- 1 for each member), and number of members
(nnames) and generic labeling information. (RPMULTI is described in
the section on routines describing the contents of the TRANSP output
database). Then, the multigraph data, or the data transformed e.g.
by a volume integral, is fetched with one of the following
calls. The data will be returned in a 2d array
The following arguments in all multigraph data access subroutine calls.
Details: (mg = short for multigraph):
PARAMETER (mgmax=15) ! max no. of members in mg's
! obtained from RPMULTI, input to multigraph read routine:
CHARACTER*21 mgnames(mgmax) ! multigraph names
INTEGER mgsigns(mgmax) ! +/-1 sign factors
INTEGER nnames ! actual no. of names, this mg.
! buffer for output data
REAL zmgbuf(<item size>,mgmax) ! where <item size> is greater
! than or equal to the number
! of data points in each member
! function.
INTEGER ibufsize ! ibufsize = <item size>
! output other than the data itself...
INTEGER iret ! actual no. data pts *per item*
INTEGER istype ! item subtype code
INTEGER iwarn ! arithmetic warning flag
INTEGER ier ! error flag
Just read the multigraph data...
...
call RPMG0CAL(mgnames,mgsigns,nnames,zmgbuf,ibufsize,
> iret,istype,iwarn,ier)
returns the multigraph data in zmgbuf(*,*), and:
iret = actual multigraph member data item size, .le.
ibufsize. note iret is the item size, not the
total amount retrieved!
istype = x axis type of data
iwarn = warning flag (0 = no warning)
ier = error flag (0 = no error)
Read the multigraph data, running each item through a one step
calculator program
RPLOT calculator expressions are used to return standard transforms of
multigraph data. For example,
character*50 zxpr
zxpr='volint(@)' ! multigraph names substituted at "@" symbol.
call RPMG1CAL(mgnames,mgsigns,nnames,zxpr,zmgbuf,ibufsize,
> iret,istype,iwarn,ier)
Will cause zmgbuf(*,*) to be filled with a set of profiles of volume
integral transforms of the original data. An error flag would be set if
the named items were of the wrong type to be volume integratable.
On output,
iret = actual multigraph member data item size, .le.
ibufsize. note iret is the item size, not the
total amount retrieved!
istype = x axis type of data
iwarn = arithmetic warning flag (0 = no warning)
ier = error flag (0 = no error)
If iwarn is set, it means arithmetic errors occurred during
evaluation of the calculator expression.
Note that the calculator expression can change `istype' to something
other than the `istype' of the original multigraph contents.
Multi-step calculations can also be performed:
parameter (maxprog=10)
character*512 cprog(maxprog)
C 2 step calculation on each member:
prog(1)= <whatever#1>
prog(2)= <whatever#2>
call RPMG2CAL(mgnames,mgsigns,nnames,prog(1),prog(2),zmgbuf
> ibufsize, iret,istype,iwarn,ier)
C 3 step calculation on each member:
prog(1)= <whatever#1>
prog(2)= <whatever#2>
prog(3)= <whatever#3>
call RPMG3CAL(mgnames,mgsigns,nnames,
> prog(1),prog(2),prog(3),zmgbuf,
> ibufsize, iret,istype,iwarn,ier)
C arbitrary sized calculation
prog(1) = <whatever#1>
...
prog(nsteps) = <last-whatever>
call RPMGCALC(mgnames,mgsigns,nnames,prog,nsteps,zmgbuf,
> ibufsize, iret,istype,iwarn,ier)
C an example of a multistep calculation would be:
C ("@" symbol shows where to substitute multigraph data in expressions).
call RPMG2CAL(mgnames,mgsigns,nnames,
> '_tmp_norm = %TIME_TRACE(I1.0, VOLAVG(@))',
> '@/_tmp_norm'
> zmgbuf,ibufsize, iret,istype,iwarn,ier)
C ...the first step extracts as a function of time the last point
C of the volume average transform of multigraph members (i.e. the
C volume average over the entire plasma), and then, the data returned
C is the original data "@" renormalized by dividing out _tmp_norm.
C
C **history note**
C This is the "H(r)" transform traditionally used by experimental physicists
C to characterize the peakedness of neutral beam deposition or heating
C profiles.
On output,
iret = actual multigraph member data item size, .le.
ibufsize. note iret is the item size, not the
total amount retrieved!
istype = x axis type of data
iwarn = arithmetic warning flag (0 = no warning)
ier = error flag (0 = no error)
If iwarn is set, it means arithmetic errors occurred during
evaluation of the calculator expression.
Note that the calculator expression can change `istype' to something
other than the `istype' of the original multigraph contents.
T1SCALAR -- extract a scalar number from f(t) at the selected time
T1PROFIL -- extract a profile from f(x,t) at the selected time
T1MHDEQ -- get the TRANSP MHD equilibrium at the selected time
RPTIMAV -- extract time slice(s) from time dependent data already
read into memory.
fortran-90 code for extracting scalar data from the currently
connected TRANSP run.
subroutine t1scalar(zname,zlabel,zunits,ztime,zdelta,zvalue,ierr)
!
! fetch a scalar function of time and interpolate to the indicated
! time value (ztime), optionally with averaging (+/- zdelta)
!
character(*), intent(in) :: zname ! name of scalar function
character(*), intent(out) :: zlabel ! function label
character(*), intent(out) :: zunits ! function units label
!
real, intent(in) :: ztime ! time (seconds) to which to interpolate
real, intent(in) :: zdelta ! +/- averaging time (seconds)
!
real, intent(out) :: zvalue ! the interpolated or averaged value.
!
integer, intent(out) :: ierr ! completion code, 0=OK
fortran-90 code for extracting profile data from the currently
connected TRANSP run.
subroutine t1profil(zname,zlabel,zunits,ztime,zdelta, &
istype,zprofil,nmax,ngot,ierr)
!
! fetch a profile function of time and interpolate to the indicated
! time value (ztime), optionally with averaging (+/- zdelta)
!
character(*), intent(in) :: zname ! name of profile function
character(*), intent(out) :: zlabel ! function label
character(*), intent(out) :: zunits ! function units label
!
real, intent(in) :: ztime ! time (seconds) to which to interpolate
real, intent(in) :: zdelta ! +/- averaging time (seconds)
!
integer, intent(out) :: istype ! type of x-axis associated with profile
integer, intent(in) :: nmax ! dimension of zprofil array
real, dimension(nmax), intent(out) :: zprofil ! the output profile
integer, intent(out) :: ngot ! actual size of profile returned
!
integer, intent(out) :: ierr ! completion code, 0=OK
fortran-90 code for extracting an MHD equilibrium dataset from a
TRANSP run. There is an MKS units conversion option.
subroutine t1mhdeq(ztime,zdelta,nsmax,ntheta,nsgot,units_option, &
theta, &
Rarr,Zarr, &
rho,psi,pmhd,qmhd,gmhd, &
tflux, pcur, ierr )
!
implicit none
!
! read in the TRANSP MHD equilibrium at a particular time
! or averaged over a range of times
!
! in seconds:
!
real, intent(in) :: ztime ! time at which equilibrium is wanted
real, intent(in) :: zdelta ! +/- zdelta for time average, or 0.0
!
! dimensioning info:
!
integer, intent(in) :: nsmax ! max no. of surfaces, counting axis
integer, intent(in) :: ntheta ! no. of theta points in Rarr,Zarr
integer, intent(out) :: nsgot ! actual no. of surfaces counting axis
!
! output physical units option
!
character(*), intent(in) :: units_option
! "TRANSP" for TRANSP traditional units
! "MKS" for MKS units:
! TRANSP MKS
! Rarr,Zarr cm m
! rho (dimensionless)
! psi Webers/rad (same for both)
! Pmhd Pascals (same for both)
! Qmhd (dimensionless)
! Gmhd T*cm T*m
! Tflux Webers (same for both)
! Pcura Amps (same for both)
!
real, intent(in), dimension(ntheta) :: theta ! theta grid, input
!---------------------------
! equilibrium geometry:
real, intent(out), dimension(ntheta,nsmax) :: Rarr ! R(theta,rho)
real, intent(out), dimension(ntheta,nsmax) :: Zarr ! Z(theta,rho)
!
! normalized sqrt(toroidal flux):
real, intent(out), dimension(nsmax) :: rho ! sqrt(phi)/sqrt(philim)
!
! poloidal flux profile
real, intent(out), dimension(nsmax) :: psi ! Webers/radian
!
! MHD pressure
real, intent(out), dimension(nsmax) :: pmhd ! (see units_option)
!
! q profil
real, intent(out), dimension(nsmax) :: qmhd
!
! g profile
real, intent(out), dimension(nsmax) :: gmhd ! (see units_option)
!
! total enclosed flux
real, intent(out) :: tflux ! Webers
!
! total toroidal plasma current
real, intent(out) :: pcur ! Amps
!
!----> completion code
!
integer, intent(out) :: ierr ! completion code, 0=OK
If time dependent data has been read in, RPTIMAV can be used to
extract a timeslice.
RPTIMAV time interpolates individual function or multigraph
data to a specified time, or, optionally, time averages at the
specified time +/- delta_t (for interpolation only, use
delta_t=0.0):
input:
INTEGER ibufsize,nfcns ! buffer dimension information
REAL zmgbuf(ibufsize,nfcns) ! data buffer
INTEGER istype ! data subtype code
REAL ttarg ! time-of-interest (seconds)
REAL delta_t ! +/- averaging time (seconds)
INTEGER ioutsize ! output buffer size / result
output:
REAL outbuf(ioutsize,nfcns) ! output data buffer
INTEGER ier ! completion code, 0=OK
...
call RPTIMAV(zmgbuf,ibufsize,nfcns,istype,ttarg,delta_t,
> outbuf,ioutsize, ier)
C
C where:
C zmgbuf = the data to be interpolated or averaged
C from a prior RPMGCALC or RPMG*CAL call, or an
C individual function or calculator result (set nfcns=1).
C
C ibufsize = 1st array dimension of zmgbuf, (max) size per item.
C
C nfcns = number of functions to time average
C
C istype = function subtype code (will determine actual item size)
C
C ttarg = time to interpolate to or average around
C delta_t = average +/- delta_t around ttarg
C if delta_t is NEGATIVE then 1/[time average of (1/data)]
C is computed. ("double inverse averaging").
C
C output:
C outbuf = array into which to write results of interpolation
C ioutsize = (max) size per result
C
C ier = completion code, 0 = normal
C
C RPTIMAV declares
C real zmgbuf(ibufsize,nfcns)
C real outbuf(ioutsize,nfcns)
C
** GUARD character **
The RPLOT calculator supports "expressions" and "commands" (for details
see ftp://ftp.pppl.gov:/pub/transp/codata/doc/rpcalc.doc). The parser
uses a leading "guard character" to determine that the following
string is a command rather than an expression. This guard character
can be reset, although care should be taken not to reuse a character
that has other meaning.
The default guard character is "%". To get the current guard character,
CHARACTER*1 cguard
CALL RPGETGC(cguard)
or, to SET the guard character,
cguard='!'
CALL RPSETGC(cguard)
(The guard character should not be reset unless the programmer feels
strongly that "%" is just too ugly... the should first carefully read
the calculator documentation and help messages before selecting a
new guard character).
** HELP messages **
The following routines generate ascii output, in the file allocatd
for messages (see initialization). They describe aspects of the
RPLOT calculator:
call rpgetgc(cguard)
call plcintro(cguard) -- introductory information; "%" as guard
character for RPLOT calculator commands.
call plclis -- list of arithmetic expression operators
call plcmds(cguard) -- list of calculator commands.
** UNITS conversion **
Generally, the calculator cannot automatically process units labels
of inputs to determine the physical units of the output of a calculator
expression. However, for a few simple operators, units conversions
are supported. These operators are:
VOLINT(...) -- volume integral (#/cm3 --> #)
ARINT(...) -- area integral (A/cm2 --> A)
FLXINT(...) -- "flux" integral (#/cm3/sec --> #/cm2/sec)
GRAD(...) -- gradient (anything --> anything/cm)
LOGDERIV(...) -- inv. scale length (anything --> cm**-1)
SCALEN(...) -- scale length (anything --> cm)
TIME_DERIV(...) -- time derivative (anything --> anything/sec)
It may be useful to access these transformations e.g. to modify
the units string from a multigraph when plotting a volume
integral transform. To do so:
CHARACTER*16 units
CHARACTER*10 transform
INTEGER iok
CALL RPFIXUNS(units, transform, iok)
details:
subroutine rpfixuns(units,transform,iok)
C
C call up an RPLOT units transformation.
C
C does stuff like #/cm3/sec -> #/sec (e.g. for "VOLINT")
C
C in/out
character*(*) units ! units to be transformed
C
C in
character*(*) transform ! type of transform
C
C out
integer iok ! =1 if transform type is known.
C
C if iok.ne.1 on exit, units are unchanged. no message is generated.
C
C iok=0 means transform not recognized.
C iok=2 means transform recognized but units change call did not
C change the units string
C
example of a subroutine which uses several of the calls described
above: trread/rpcalcd.for:
subroutine rpcalcd(zexpr, rlabel,
> zdata, zxdata, ndmax, ztdata, ntmax,
> nxgot, ntgot, istype, iwarn, ier)
c
c call the rplot calculator routine rpcalc.
c return a copy of the data computed, the corresponding x axis data
c if applicable, and the corresponding timebase data.
c
c input:
character*(*) zexpr ! calculator expression (cf rpcalc).
character*(*) rlabel ! run label (for error message)
c
c storage buffers filled on output:
integer ndmax ! max size of data item (INPUT)
integer ntmax ! max size of timebase (INPUT)
c
real zdata(ndmax) ! buffer for calculator results data
real zxdata(ndmax) ! buffer for x axis data
real ztdata(ntmax) ! buffer for timebase
c
c scalars set on output:
c
integer nxgot ! no. of x pts, 1 for scalar f(t) data
integer ntgot ! no. of time pts in timebase
c
c ntgot words are written in ztdata;
c nxgot*ntgot words are written in zdata
c if (istype.ge.1) nxgot*ntgot words are written in zxdata
c ...the x variation is stored contiguously, i.e. the fortran 2d
c array declaration would be array(nxgot,ntgot)
c
integer istype ! data item type code
c
c istype=-1 -- scalar f(t)
C istype= 1 -- f(x,t), x is TRANSP zone ctrs if this is TRANSP data
C istype= 2 -- f(x,t), x is TRANSP zone bdys if this is TRANSP data
C ..etc..
C istype = 0 indicates an error
C
integer iwarn ! arithmetic warning, 0 = OK
integer ier ! completion code, 0 = OK
c
c--------------------------------------------------------------------
c
integer str_length,lunzer
c
integer irank,idims(10)
character*21 zxabb
c
c--------------------------------------------------------------------
c
c OK... call the calculator
c
ilunz=lunzer(0)
c
call rpcalc(zexpr,zdata,ndmax,igot,istype,iwarn,ier)
c
c message on warning or error
c
if(max(ier,iwarn).ne.0) then
ile=str_length(zexpr)
ilb=str_length(rlabel)
write(ilunz,1001) rlabel(1:ilb),zexpr(1:ile),iwarn,ier
1001 format(/
> '%rpcalcd: runid: ',a/
> ' expression: ',a/
> ' iwarn=',i6,' ier=',i6)
endif
if(ier.ne.0) then
nxgot=0
ntgot=0
return
endif
c
if(istype.eq.-1) then
c
c if a scalar function result: have data already, just get timebase
c
nxgot=1
ntgot=igot
call rptime_s(ztdata,ntmax,igot)
return
else
c
c if a profile function result: have data already, get profile
c timebase and x axis function
c
call rpdims(istype,irank,idims,zxabb,ier)
if(irank.gt.2) then
write(ilunz,*)
> ' ??rpcalcd: irank.gt.2 object not supported!'
nxgot=0
ntgot=0
ier=99
else
nxgot=idims(1)
ntgot=idims(2)
endif
c check time dimension
if(ntgot.gt.ntmax) then
write(ilunz,*) ' ??rpcalcd: ntgot.gt.ntmax, ntgot=',ntgot,
> ' ntmax=',ntmax
nxgot=0
ntgot=0
ier=88
else
c get timebase
call rptime_p(ztdata,ntmax,igot)
c get x axis
call rprofile(zxabb,zxdata,ndmax,iret, ier)
endif
c
endif
c
return
end
The following routines are available:
rd_nspecies -- return information on total number of distinct
plasma species, and numbers in each of several
categories of plasma species
rd_species -- get list of all plasma species
rd_thspec -- get list of non-impurity thermal species
rd_thxspec -- get list of impurity thermal species
rd_bmspec -- get list of beam species
rd_rfspec -- get list of RF tail species
rd_fuspec -- get list of fusion product species
subroutine rd_nspecies(n_species,n_thi,n_thx,n_bi,n_rfi,n_fusi)
!
! return the number of species in the (currently open) TRANSP run
! return n_species=0 if no run is currently open.
! return n_species=-1 if some other error occurred (rare).
!
implicit NONE
!
integer, intent(out) :: n_species ! total no. of species including electrons
integer, intent(out) :: n_thi ! no. of "non-impurity" therm. ion species
integer, intent(out) :: n_thx ! no. of "impurity" thermal ion species
integer, intent(out) :: n_bi ! no. of "beam" ion species
integer, intent(out) :: n_rfi ! no. of "RF tail" ion species
integer, intent(out) :: n_fusi ! no. of "fusion product" ion species.
subroutine rd_species(maxn,n_species,zlbla,abray,itype,ifast,zz,aa,izc)
!
! find all plasma species -- return information in passed arrays
!
implicit NONE
!
integer, intent(in) :: maxn ! max no. of species (array sizes)
integer, intent(out) :: n_species ! actual no. of species (-1 if error)
!
character*20, intent(out) :: zlbla(maxn) ! descr. label for each specie
!
character*21, intent(out) :: abray(4,maxn) ! specie profile names
!
! for specie j:
! abray(1,j) = name of density profile-- for all species (n/cm**3).
! abray(2,j) = name of temperature or "average energy" profile.
! abray(3,j) = name of perpendicular energy density profile or
! average perp energy per particle or " ".
! abray(4,j) = name of parallel energy density profile or
! average pll energy per particle or " ".
! (generally abray(3:4,j) will be set for fast species only)
!
integer, intent(out) :: itype(maxn) ! species type code
!
! 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
!
! itype(j)=+4 -- beam ion specie
! itype(j)=+5 -- rf tail ion specie
! itype(j)=+6 -- fusion product ion specie
!
integer, intent(out) :: ifast(maxn) ! thermal/fast flag
!
! ifast(j)=0 -- thermal electron or ion specie,
! abray(2,j) contains temperature, eV, abray(3:4,j)=" ".
! ifast(j)=1 -- fast specie, abray(2,j) is avg energy / ptcl <E>, eV
! abray(3:4,j) contain perp, pll energy density, J/cm**3
! ifast(j)=2 -- fast specie, abray(2,j) is "temperature" = 2/3 <E>, eV
! abray(3:4,j) contain perp, pll <E>/ptcl, eV
!
! (other fast specie configurations may be added later)
!
real, intent(out) :: zz(maxn) ! charge of specie
real, intent(out) :: aa(maxn) ! mass of specie (amu)
integer, intent(out) :: izc(maxn) ! atomic number of species
!
!---------------------------------
The following routines have a subset of the arguments of rd_species.
For argument details see rd_species.
The argument "ngot" (integer, intent(out)) specifies the number of
species found, which can be zero for routines other than rd_thspec.
subroutine rd_thspec(maxn,abray,zz,aa,izc,ngot) ! non-impurity, thermal
subroutine rd_thxspec(maxn,abray,zz,aa,izc,ngot) ! impurity, thermal
note: if rd_thxspec returns ngot=1, zz=aa=0, this means that the
impurity Z and A are functions of time: read the scalar
functions "XZIMP" and "AIMP".
subroutine rd_bmspec(maxn,abray,zz,aa,izc,ngot) ! beam species (ifast=1)
subroutine rd_rfspec(maxn,abray,zz,aa,izc,ngot) ! RF tail species (ifast=2)
subroutine rd_fuspec(maxn,abray,zz,aa,izc,ngot) ! fusion products (ifast=1)
We are open to suggestions for other useful aggregate lists.
To get in touch with the `trread' authors send email to
dmccune@pppl.gov
The trread library includes a facility for reading and parsing
TRANSP namelists.
Summary: first read the namelist. Then, named items can be
extracted.
---> To read the namelist:
INTEGER ilun ! fortran i/o unit number (input)
INTEGER ier ! status code returned (0=OK)
call tr_getnl_text(ilun,ier) ! read in the namelist
---> To extract named items
CHARACTER*32 ZNAME ! name of desired item (input)
INTEGER MAXVAL ! array dimension (input)...
INTEGER intvec(MAXVAL) ! INTEGER vector (output)
LOGICAL logvec(MAXVAL) ! LOGICAL vector (output)
REAL r4vec(MAXVAL) ! standard REAL vector (output)
REAL*8 r8vec(MAXVAL) ! REAL*8 vector (output)
INTEGER istat ! status code (output)
!
! on output:
!
! normal returns:
! istat = 0 -- no values found
! istat = N, 1.le.N.le.maxval -- N values found, stored in intvec(1:N)
!
! error returns:
! istat = -1 -- error (e.g. namelist was never read)
! istat = -N, N.gt.maxval -- N values were found, intvec(1:maxval) set
! to the first maxval of them, there are too many values to
! return them all
!
CALL TR_GETNL_INTVEC(zname,intvec,maxval,istat)
CALL TR_GETNL_LOGVEC(zname,logvec,maxval,istat)
CALL TR_GETNL_R4VEC(zname,r4vec,maxval,istat)
CALL TR_GETNL_R8VEC(zname,r8vec,maxval,istat)
---> Artificial Example:
! declarations
integer ilun ! i/o unit
integer ier ! status code (namelist file read)
integer istat ! status code (namelist item lookup)
integer ngmax ! #of plasma species (namelist)
real*8 aplasm(5) ! A of plasma species (namelist)
real*8 backz(5) ! Z of plasma species (namelist)
! read namelist file
ilun=99
call tr_getnl_text(ilun,ier)
if(ier.ne.0) then
...handle error...
endif
! get number of plasma species
call tr_getnl_intvec('NGMAX',ngmax,1,istat)
if(istat.ne.1) then
...treat as error...
else if(ngmax.gt.5) then
...treat as error...
endif
! get A of plasma species
call tr_getnl_r8vec('APLASM',aplasm,ngmax,istat)
if(istat.ne.ngmax) then
...treat as error...
endif
! get Z of plasma species
call tr_getnl_r8vec('BACKZ',backz,ngmax,istat)
if(istat.ne.ngmax) then
...treat as error...
endif
(A more substantial example is in the trxplib source code: see
routines trx_wall.f90, trx_getnlims.f90, trx_getcrlim.f90, and
trx_getlnlim.f90)
---> Cautionary Notes...
The TRANSP code maintains default values for namelist variables.
This means that some runs may rely on default values for some items,
in which case there will be no specific mention of the item in the
namelist file. **The software described here does not have access
to the default values at this time**. Any attempt to read a defaulted
variable will yield istat=0, not found. If there is a requirement
to change this, something can be done. However, reliability will be
less than 100% over time, because, the TRANSP code evolves, the namelist
evolves, and the code's default values for namelist controls can change.
On the other hand, certain information cannot be defaulted, or the
default has a well known and stable interpretation, so that a namelist
read can be reliably employed. Bottom line: seek advice from a TRANSP
expert when contemplating a namelist read to fetch TRANSP data.
A final note: the tr_getnl_*vec(...) routines can fetch scalar or
singly dimensioned vector integer, logical, and floating point data.
At present there is no support for fetching character data or data
elements from multiply dimensioned arrays. Thus there are a small
number of items in the TRANSP namelist that cannot be accessed with
the software as currently implemented.
Most of the trread software is organized around the idea of
"connecting" to a single "primary" run, and then reading data
out of that run exclusively.
However, the following routines allow access to data from "secondary"
runs without disrupting the connection to the primary run.
Indeed, these routines are the legacy callable interface to TRANSP
data. The earlist versions of TRPROFIL and TRSCALAR are over ten
years old.
Summary:
TRSCALAR -- read a named function f(t) from a TRANSP run
TRPROFIL -- read a named function f(x,t) from a TRANSP run
These routines have similar arguments:
! input:
CHARACTER*(*) DISK ! VMS DISK **or** MDS+ SERVER & TREE NAME
CHARACTER*(*) DIR ! UNIX PATH **or** PPPL MDS+ tok.yy ID
CHARACTER*(*) RUNID ! TRANSP RUND **or** non-PPPL MDS+ shot no.
CHARACTER*(*) ABBREV ! TRANSP id for data item sought
INTEGER MAXTIMES ! max. no. of time points in data returned.
! input, TRPROFIL ONLY:
INTEGER MAXDATA ! data buffer size limit
! output:
CHARACTER*32 LABEL ! TRANSP 32 character label for item
CHARACTER*16 UNITS ! TRANSP 16 character units label for item
! output, TRPROFIL ONLY:
INTEGER ITYPE ! X-axis type code for data returned
INTEGER INX ! no. of X points in data returned.
! output:
INTEGER NTIMES ! number of time points returned.
REAL TIMES(MAXTIMES) ! TIMES(1:NTIMES) = the time points returned.
! output, TRSCALAR ONLY:
REAL SDATA(MAXTIMES) ! SDATA(1:NTIMES) = the data points returned.
! output, TRPROFIL ONLY:
REAL DATA(MAXDATA) ! DATA(1:NTIMES*INX) = the data points
! returned. ordering: INX points at
! TIME(1), then INX points at TIME(2),
! etc., contiguously stored.
! output:
INTEGER IERR ! completion code, 0=OK (error codes see below)
... ...
CALL TRSCALAR( DISK,
1 DIR, RUNIDIN, ZABBREV,MAXTIMES,
2 LABEL, UNITS,
3 NTIMES, TIMES, SDATA, IERR )
... ...
CALL TRPROFIL( DISK,
1 DIR, RUNIDIN, ZABBREV,MAXTIMES, MAXDATA,
2 LABEL, UNITS, ITYPE, INX,
3 NTIMES, TIMES, DATA, IERR )
! error codes...
C IERR - INTEGER ERROR MESSAGE. RETURNED = 0 IF OK.
C 0 = OK
C 2 = ERROR READING MF FILE.
C 3 = UNEXPECTED EOF READING MF FILE.
C 5 = TRIED TO READ OLD FORMAT TF FILE.
C 15 = PROFILE ABBREVIATION NOT FOUND
C 16 = COULDN'T OPEN MF FILE
C 17 = COULDN'T OPEN TF FILE
C 26 = RUNID IN FILE IS DIFFERENT THAN ASKED FOR
C 101 = # OF TIMES IN FILE > MAXTIMES.
C MAXTIMES TIMES WERE READ IN.
C 102 = # OF PROFILES > DIMENSION OF BUFFER.
C NO DATA READ IN.
c 103 = TIMES RETURNED ARE NOT MONOTONICALLY
C INCREASING. DATA IS ALSO RETURNED.
C 104 = # OF TIMES * # IN X DIRECTIONS > MAXDATA
C 105 = # OF POINTS IN X DIRECTION = 0
C
-----------------
MDS+ examples:
PPPL TRANSP archives -- runs are identified by tok.yy & runid string;
treename = TRANSP
DISK = "MDS+-<server_name>@<treename>"
DIR = "tok.yy"
RUNID = "<TRANSP_runid>"
! fetch plasma current (PCUR) vs. time...
CALL TRSCALAR( 'MDS+-TRANSPGRID.PPPL.GOV@TRANSP',
> 'JET.00', '50632C02', 'PCUR', MAXTIMES,
2 LABEL, UNITS,
3 NTIMES, TIMES, SDATA, IERR )
non-PPPL servers -- runs are identified by treename & shot number
DISK = "MDS++<server_name>@<treename>"
DIR = " "
RUNID = "<shot_number>" ! in ascii
! fetch electron temperature profile evolution data...
CALL TRPROFIL( 'MDS++CMODA.PSFC.MIT.EDU@TRANSP01',
> ' ', '123456789', 'TE', MAXTIMES, MAXDATA,
2 LABEL, UNITS, ITYPE, INX,
3 NTIMES, TIMES, DATA, IERR )
-----------------
non MDS+ examples
NetCDF or legacy format TRANSP run on UNIX filesystem:
DISK = " "
DIR = "<unix path>"
RUNID = "<TRANSP_runid>"
! fetch electron temperature profile evolution data...
CALL TRPROFIL( ' ', '/u/transp/results/JET.98',
> '46243C02', 'TE', MAXTIMES, MAXDATA,
2 LABEL, UNITS, ITYPE, INX,
3 NTIMES, TIMES, DATA, IERR )
(files are at /u/transp/results/JET.98/46243C02* ... the DIR string
can also have a leading environment variable, which will be
translated, e.g. "RUNDATA/JET.98").
NetCDF or legacy format TRANSP run on VMS filesystem:
DISK = "<disk name>", trailing ":" optional
DIR = "<directory>", no brackets
RUNID = "<TRANSP_runid>"
! fetch electron temperature profile evolution data...
CALL TRPROFIL( 'RUNDATA', 'TRANSP.JET.98',
> '46243C02', 'TE', MAXTIMES, MAXDATA,
2 LABEL, UNITS, ITYPE, INX,
3 NTIMES, TIMES, DATA, IERR )
(files are at RUNDATA:[TRANSP.JET.98]46243C02*.*)
There is a special routine for the situation where it is desired to
read a single data item (scalar or profile function) from several
distinct runs. A call to this routine causes a multigraph to be
created containing items of the form
<item>$<runid_1>
<item>$<runid_2>
<item>$<runid_3>
...etc...
(up to 15 runs)
--> note that this causes data ids with name lengths greater than 10
to be created; indeed, throughout TRREAD the maximum length of a
data identifier or "abbreviation" has been increased from 10 to 21
characters, to accomodate the $<runid> suffix.
For example:
character*10 mgname
character*10 item
integer :: numruns=3
character*80 paths(numruns)
character*10 aliases(numruns)
integer inum_cur ! set inum_cur to indicate "current run"
! i.e. the run accessed via a prior
! KCONNECT or RCONNECT call
integer iwarn(numruns),ierr
mgname = 'neut_mg' ! multigraph name
item = 'neutt' ! function id -- total neutrons, f(t)
inum_cur = 0 ! (not required to be set)
paths(1) = '/transp/results/dmccune/output/TFTR.88/37065S10'
paths(2) = 'MDS+:TRANSPGRID.PPPL.GOV:TRANSP(TFTR.88,37065K12)'
paths(3) = 'MDS+:TRANSPGRID.PPPL.GOV:TRANSP(TFTR.88,37065K13)'
aliases(1) = '37065S10'
aliases(2) = '37065K12'
aliases(3) = '37065K13'
(Here runids are used for aliases but in another context, something
else e.g. MDS+ tree ids, or, shot numbers, might want to be used).
call run_mg(mgname,rpitem,numruns,paths,aliases,inum_cur, &
& iwarn,ierr)
ierr is set if there is an error such that the multigraph
"neut_mg" cannot be created; iwarn is set for any run where
either an item could not be accessed or an existing item in
the TRREAD session was replaced due to the name already being
in use.
if the call is successful, scalar functions
NEUTT$37065S10
NEUTT$37065K12
NEUTT$37065K13
are created and grouped in the multigraph
NEUT_MG
calculator expressions such as
DNEUTT,"Neutron Difference","N/SEC" = NEUTT$37065K13 - NEUTT$37065K12
can be formed, and the multigraph contents can be manipulated by
calculator commands e.g. %MG_ADDFUN and %MG_DELFUN: for standard
information on the calculator see
ftp://ftp.pppl.gov:/pub/transp/codata/doc/rpcalc.doc
The subroutine header follows...
subroutine run_mg(mgname,rpitem,numruns,runPath,runAlias,inum_cur,iwarn,ier)
! create a multigraph of a single item (rpitem, e.g. "TE" or "PCUR")
! read from several different runs.
! create the member functions for the multigraph, forming labels of the
! form <rpitem>$<runAlias(i)> for i = 1 to numruns
! inum_cur.ne.0 indicates that this runPath corresponds to the "main"
! run to which the session is connected (via prior call to rconnect or
! kconnect), hence we just copy the data item and do not have to read
! from a separate run database (no trprofil or trscalar call needed)
implicit NONE
character*(*), intent(in) :: mgname ! name of multigraph (mg)
character*(*), intent(in) :: rpitem ! item out of which to build mg
! the multigraph label and phys. units will be taken from rpitem in
! the current run (which must exist)
integer, intent(in) :: numruns ! no. of runs (.le. 15)
character*(*), intent(in) :: runPath(numruns) ! run paths (rconnect syntax)
character*(*), intent(in) :: runAlias(numruns) ! run aliases (max 10 chars)
! for description of runPath syntax see rconnect subroutine
! note each runAlias must be unique
integer, intent(in) :: inum_cur ! if non-zero -- index to current run
integer, intent(out) :: iwarn(numruns) ! warn on data access failures
integer, intent(out) :: ier ! completion code, 0=OK
! ier.ne.0 means an error with the arguments or that *all* data access
! attempts for *all* runs failed...
! ier=1 mg <mgname> already exists
! ier=2 item <rpitem> does *not* exist in current run
! ier=3 numruns.le.0 .or. numruns.gt.15
! ier=4 inum_cur.lt.0 .or. inum_cur.gt.numruns
! ier=5 alias error:
! runAlias(i).eq.runAlias(j) for some i.ne.j
! ier=6 **all** data access attempts failed
! ier=7 (... some other error ...)
! iwarn(j)=0 is the normal return
! iwarn(j)=1 -- data access failure on runPath(j)
! iwarn(j)=2 -- <rpitem>$<runAlias(j)> already exists -- replaced.
This Document was created by hlptohtml
Written By:Manish Vachharajani(mvachhar@pppl.gov)