What follows is a small set of examples. Note that many things are illustrated in each example, so there is some point in reading several of them.
This is an example which makes a contour map from SST data in the data catalog. This is intended as an input file to the program ingrid, though it would work with any model linked to ingrid as well.
All lines not between \begin{ingrid} and \end{ingrid} are treated as comments: text beginning with a '%' is treated as a comment as well.
#! /usr/local/bin/ingrid
Initialization: when FGFILE is called, it sets the stack so that it
contains only Ingrid: i.e. just the basic PostScript-like words. We
want to do laser printer plots using the data sources, so we start by
putting those objects on the stack.
\begin{ingrid}
PS SOURCES
\end{ingrid}
We now set the names of the plot files by using setplotname. This is
optional: if a name is not specified, one is created from the time and
date.
\begin{ingrid}
(test) setplotname
\end{ingrid}
We now do the actual plot
\begin{ingrid}
CACSSTM.ATL sstt % extracts the data from the data catalog
T 5.5 VALUE % restricts the stream to only month 5.5
X Y CONTOUR % makes a contour plot using X and Y as the axes
\end{ingrid}
Note that the range of neither X nor Y is explicitly altered. This means that the X range and the Y range are those of the original data set: i.e. data for the full XY range are plotted.
In this case, the plotfile generated is called ltest.001. It can be plotted by sendsumry or lp, previewed by gs. Note that a preamble file is required, but we normally have preloaded the definitions into the printer lwII, so that the preamble file is not required to print here. The print command is thus
lp ltest.001
or
sendsumry ltest.001
sendsumry allows many options (multiple plots per page, etc) and defaults to four per page. lp prints the file full size (one plot per page).
To preview
gs psplot.xps ltest.001
psplot.xps is the preamble file for Ingrid plots.
What follows is a short example which illustrates averaging, restricting the data range, and contouring.
It also makes the point that multiple plots can be generated with a single run, and input formatting is fairly flexible.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES % need plotting words and data catalog
CACSSTA % get SST anomalies
Y -5 5 RANGE % restrict data to near equatorial
Y AVERAGE % average over Y
X T CONTOUR % generate the contour plot
CACSSTA Y 10 20 RANGE Y AVERAGE
X T CONTOUR % same plot for a different latitude range
\end{ingrid}
While it is not strictly speaking necessary to understand the stack behavior of Ingrid in order to use it, eventually the understanding can be very useful. The stack behavior of the operators in this example is as follows
CACSSTA
RANGE
stream to have the
range [lo,hi] in coordinate grid. In the first
plot it restricts Y to [-5,5].AVERAGE
stream along
coordinate grid. In this case it averages along Y.CONTOUR
gridh
along the horizontal axis and gridv along the vertical.
Both PS and SOURCES simply put objects on the stack so
that the needed definitions would be available.
There are a number of parameters that can be set to control the appearance of plots. Plot files can also be edited fairly easily. Here we alter the aspect ratio, add the data grid, add an extra title, and specify the contour interval.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES % need plotting words and data catalog
OSUSFC.DATA sstt % get SST climatology
X -60 10 RANGE % restrict to Atlantic
/XOVY AUTO def % alter aspect ratio of the plot
% AUTO means same scale on both axes
/PLOTGRID 1 def % underlays data grid
/title (test plot) def % adds additional title to plot
DATA 1 STEP % contour in steps of 1 degree
X Y CONTOUR % makes series of (monthly) contour plots
\end{ingrid}
#! /usr/local/bin/ingrid
\begin{ingrid}
/ZPLOT { Z 1000 low RANGE Y Z CONTOUR } def
/XYPLOT { Z 250 VALUE /XOVY AUTO def X Y CONTOUR } def
/select { X 180 VALUE } def % word to pick out the data we want
PS (eos) setplotname % plotting words
SOURCES % data catalog
/S LEVITUSMEAN sal select % gets salinity
/fullname /$S$ def def % shortens label
/T LEVITUSMEAN temp select def % gets temperature
/P T S Z pressure def % calculates hydrostatic pressure
/TH T S P potemp % calculates potential temperature
/fullname /$theta$ def def % shortens label
/sigma TH S P densa 1000 mul def
TH ZPLOT % plots potential temperature
sigma ZPLOT % plots sigma theta
\end{ingrid}
Note that there is a currently a bug in Ingrid that tends to appear with Levitus data. Basically it is not smart enough to read the data file in non-sequential order, so it has trouble with functions of both temperature and salinity if more that one longitude is plotted. The fix is to force Ingrid to read all the data in at once, namely
/lev LEVITUSMEAN
X -60 20 RANGE Y -30 80 RANGE % restricts the data set
Z Y X D 4 REORDER % reorders the data
% (actually only reblocks)
toNaN % uses that reordered data
def
I hope this bug will go away eventually.
This is an example of how to prepare a file in "netCDF" format that can
be used as input for INGRID using readCDF. Full documentation on
netCDF is available.
(see section `Introduction' in netCDF User's Guide.)
This example uses standard netCDF library calls. It is much easier to use a locally developed library called ODB to generate netCDF files.
In this example we will convert a data file read from standard input. The file contains sst data defined on a cylindrical longitude-latitude grid, starting from 72W and 50S and increasing eastward and northward to 14E and 60N by 2 degree longitude and latitude increments. The file contains 12x39 monthly records from January of 1950 to December of 1988.
There are several steps in the process. First you create a short
description of the file, specifying the different dimensions of the data
and generally adding as much information about the dataset as possible.
You do not, however, put the actual data in this description. Then you
run ncgen, a utility that creates a FORTRAN (or C) program from
your dataset description. Then you modify the program so that the
actual data is written to the netCDF file. Finally, after compiling the
program, you run it, creating the netCDF file. This file can then be
read by Ingrid.
We start by creating a text version of the header of the netCDF file. In the text we define the dimensions, variables and one data record: the array "month" containing the months of the year. The text can be entered in the editor. We called the file sst.header (its path is /x/kushnir/sst.header).
netcdf fssta {
dimensions:
lon = 46 ;
lat = 56 ;
month = 12 ;
year = 39 ;
variables:
float lon(lon) ;
lon:units = "longitude" ;
lon:fullname = "Longitude" ;
float lat(lat) ;
lat:units = "latitude" ;
lat:fullname = "Latitude" ;
long month(month) ;
month:units = "month" ;
long year(year) ;
year:units = "year" ;
float ssta(year, month, lat, lon) ;
ssta:missing_value = 999.90002f ;
ssta:history = "COADS/GFDL - A . O" ;
ssta:long_name = "sea surface temperature anomaly" ;
ssta:units = "celsius" ;
// global attributes:
:source = "COADS/GFDL - A. Oort" ;
:base_time = "1960 1 15 12:00" ;
data: // asign data where possible
month = 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11;
}
Note that month varies from 0 to 11 consistent with INGRID. Having a gridtype of 1 would mean that Ingrid will treat that grid as being periodic. None of the grids in this example are periodic.
The next stage is to generate a fortran (or if you want c) program that uses the information in the header and will be used later to generate the netCDF file. To do this we execute the utility "ncgen", for example:
ncgen -f sst.header > gen_CDF.f
This will put the following program in the file gen_CDF.f:
ncgen -f sst.header
program ffssta
include 'netcdf.inc'
integer iret
* netCDF id
integer cdfid
* dimension ids
integer londim, latdim, monthdim, yeardim
* variable ids
integer lonid, latid, monthid, yearid, sstaid
* variable shapes
integer dims(4)
* corners and edge lengths
integer corner(4), edges(4)
* data variables
real lon(46)
real lat(56)
integer month(12)
integer year(39)
real ssta(46,56,12,39)
* attribute vectors
real floatval(1)
* enter define mode
cdfid = nccre ('ffssta.cdf', NCCLOB, iret)
* define dimensions
londim = ncddef(cdfid, 'lon', 46, iret)
latdim = ncddef(cdfid, 'lat', 56, iret)
monthdim = ncddef(cdfid, 'month', 12, iret)
yeardim = ncddef(cdfid, 'year', 39, iret)
* define variables
dims(1) = londim
lonid = ncvdef (cdfid, 'lon', NCFLOAT, 1, dims, iret)
dims(1) = latdim
latid = ncvdef (cdfid, 'lat', NCFLOAT, 1, dims, iret)
dims(1) = monthdim
monthid = ncvdef (cdfid, 'month', NCLONG, 1, dims, iret)
dims(1) = yeardim
yearid = ncvdef (cdfid, 'year', NCLONG, 1, dims, iret)
dims(4) = yeardim
dims(3) = monthdim
dims(2) = latdim
dims(1) = londim
sstaid = ncvdef (cdfid, 'ssta', NCFLOAT, 4, dims, iret)
* assign attributes
call ncaptc(cdfid, lonid, 'units', NCCHAR, 9, 'longitude', iret)
call ncaptc(cdfid,lonid,'fullname',NCCHAR,9, 'Longitude', iret)
call ncaptc(cdfid, latid, 'units', NCCHAR, 8, 'latitude', iret)
call ncaptc(cdfid, latid,'fullname',NCCHAR,8, 'Latitude', iret)
call ncaptc(cdfid, monthid, 'units', NCCHAR, 5, 'month', iret)
call ncaptc(cdfid, yearid, 'units', NCCHAR, 4, 'year', iret)
floatval(1) = 999.90002
call ncapt(cdfid, sstaid, 'missing_value', NCFLOAT, 1, floatval, i
1ret)
call ncaptc(cdfid, sstaid, 'history', NCCHAR, 20, 'COADS/GFDL - A
1 . O', iret)
call ncaptc(cdfid, sstaid, 'long_name', NCCHAR, 31, 'sea surface t
1emperature anomaly', iret)
call ncaptc(cdfid, sstaid, 'units', NCCHAR, 7, 'celsius', iret)
call ncaptc(cdfid, NCGLOBAL, 'source', NCCHAR, 20, 'COADS/GFDL - A
1. Oort', iret)
call ncaptc(cdfid, NCGLOBAL, 'base_time', NCCHAR, 15, '1960 1 15 1
12:00', iret)
* leave define mode
call ncendf(cdfid, iret)
* store month
corner(1) = 1
edges(1) = 12
data month /0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11/
call ncvpt(cdfid, monthid, corner, edges, month, iret)
call ncclos (cdfid, iret)
end
ncgenerated programThere are now a few things we have to add to this program so that we can use it to convert the data file. First we have to define the full path for the include file: "netcdf.inc". This is: /usr/local/include/netcdf.inc. Next, for reasons explained below, we redimension ssta(46,56). Finally, since we only specified the values of the array "month" we need to modify the program to generate or read in the other array values. For each array stored the subroutine "ncvpt" has to know the values of the arrays "corner" and "edges". The array "corner" specifies the first value of each of the data array indices (by the order they appear in the dimension statement). The array "edges" specifies the length of the segment stored in each call to "ncvpt" in terms of each of the indices of the array. Thus if the data array is dimensioned i1,i2, and i3 we can simply do:
corner(1)=1 corner(2)=1 corner(3)=1 edges(1)=i1 edges(2)=i2 edges(3)=i3 call ncvpt(cdfid,dataid,corner,edges,data,iret)
Alternatively, if we prefer to store the array with several calls to "ncvpt", each time putting out an i1xi2 field, we can do:
do j=1,i3
corner(1)=1
corner(2)=1
corner(3)=j
edges(1)=i1
edges(2)=i2
edges(3)=1
call ncvpt(cdfid,dataid,corner,edges,data(1,1,j),iret)
end do
In the example below, we use the latter form to output the data in "ssta". This array is defined in the program as having a dimension of 46x56. In the netCDF file it is recognized as having a dimension of 46x56x12x39. We read this array (from standard input) in segments of 46x56 and loop through the other two dimensions properly defining the values of "corner" and "edges". The arrays "lat", "lon", and "year" are calculated in the program and stored too. Note that year was defined as being =0 in 1960, to fit with the data catalog and INGRID's procedures. This is the final version of the program:
*-----------------------
* /x/kushnir/gen_CDF.f
*-----------------------
program fssta
include '/usr/local/include/netcdf.inc'
integer iret
* netCDF id
integer cdfid
* dimension ids
integer londim, latdim, monthdim, yeardim
* variable ids
integer sstaid, lonid, latid, monthid, yearid
* variable shapes
integer dims(4)
* corners and edge lengths
integer corner(4), edges(4)
* data variables
real ssta(46,56)
real lon(46)
real lat(56)
integer month(12)
integer year( 39)
character*4 varnam
* open output file
* attribute vectors
* enter define mode
cdfid = nccre ('fssta.cdf', NCCLOB, iret)
* define dimensions
londim = ncddef(cdfid, 'lon', 46, iret)
latdim = ncddef(cdfid, 'lat', 56, iret)
monthdim = ncddef(cdfid, 'month', 12, iret)
yeardim = ncddef(cdfid, 'year', 39, iret)
* define variables
dims(1) = londim
lonid = ncvdef (cdfid, 'lon', NCFLOAT, 1, dims, iret)
dims(1) = latdim
latid = ncvdef (cdfid, 'lat', NCFLOAT, 1, dims, iret)
dims(1) = monthdim
monthid = ncvdef (cdfid, 'month', NCLONG, 1, dims, iret)
dims(1) = yeardim
yearid = ncvdef (cdfid, 'year', NCLONG, 1, dims, iret)
dims(4) = yeardim
dims(3) = monthdim
dims(2) = latdim
dims(1) = londim
sstaid = ncvdef (cdfid, 'ssta', NCFLOAT, 4, dims, iret)
* assign attributes
call ncapt(cdfid, sstaid, 'missing_value', NCFLOAT, 1, 999.9,
1 iret)
call ncapt(cdfid, sstaid, 'history', NCCHAR, 20,
1'COADS/GFDL - A. Oort', iret)
call ncaptc(cdfid, sstaid, 'long_name', NCCHAR, 31, 'sea surface t
1emperature anomaly', iret)
call ncaptc(cdfid, sstaid, 'units', NCCHAR, 7, 'celsius', iret)
call ncaptc(cdfid, lonid, 'units', NCCHAR, 9, 'longitude', iret
1)
call ncaptc(cdfid, latid, 'units', NCCHAR, 8, 'latitude', iret)
call ncaptc(cdfid, monthid, 'units', NCCHAR, 5, 'month', iret)
call ncaptc(cdfid, yearid, 'units', NCCHAR, 4, 'year', iret)
call ncaptc(cdfid, NCGLOBAL, 'source', NCCHAR, 20, 'COADS/GFDL - A
1. Oort', iret)
call ncaptc(cdfid, NCGLOBAL, 'base_time', NCCHAR, 15, '1960 1 15 1
12:00', iret)
* leave define mode
call ncendf(cdfid, iret)
* store lon
corner(1) = 1
edges(1) = 46
do nl=1,46
lon(nl)=-76.+(nl-1)*2
end do
WRITE(6,*) lon
call ncvpt(cdfid, lonid, corner, edges, lon, iret)
* store lat
corner(1) = 1
edges(1) = 56
do nl=1,56
lat(nl)=-50+(nl-1)*2
end do
WRITE(6,*) lat
call ncvpt(cdfid, latid, corner, edges, lat, iret)
* store month
corner(1) = 1
edges(1) = 12
data month /0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11/
call ncvpt(cdfid, monthid, corner, edges, month, iret)
* store year
corner(1) = 1
edges(1) = 39
do ny=1, 39
year(ny)=(ny-11)*12
end do
call ncvpt(cdfid, yearid, corner, edges, year, iret)
* read data from standard input
do ny=1, 39
do nm=1,12
read (5,'(a4,i2,i4,/,(16f5.1))') varnam,imon,iyr,ssta
corner(1) = 1
corner(2) = 1
corner(3) = nm
corner(4) = ny
edges(1) = 46
edges(2) = 56
edges(3) = 1
edges(4) = 1
call ncvpt(cdfid, sstaid, corner, edges, ssta, iret)
end do
end do
call ncclos (cdfid, iret)
end
C Local Variables:
C compile-command: "f77 -o gen_CDF gen_CDF.f -lnetcdf -lsun"
C End:
In compiling the program you should link to the two libraries: "-lnetcdf" and "-lsun".
When using ingrid remember to call the variables by the same name you have used in the netCDF file. Thus if you want to contour the "ssta" field use:
(fssta.cdf) readCDF ssta lon lat CONTOUR
and not:
(fssta.cdf) readCDF ssta X Y CONTOUR
(Alternatively, you could call your grids X and Y).
Good luck!
This is an example of a simple color plot. It is mostly to demonstrate how to modify the color scale.
There is a default color scale, so that it is not necessary to specify one.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES % basic plotting words and data catalog
(ps) setplotname % set plot name
OSUSFC.DATA sstt T 1 VALUE % extract first month of sstt data
startcolormap % changing the color map: start
0. 30. RANGE % specify the data range
black % land color is first
darkgray % color for low values is second
blue 0. value % 0 is blue
green 15.value % 15. is green (smoothly interpolates colors)
255 255 0 RGB stripe % add a yellow strip at 15 (specifying as RGB)
green % go back to green
red 30. value % 30. is red (smoothly interpolates colors)
grey % color for high values
endcolormap % end of color map setting
X Y COLOR % makes color plot
\end{ingrid}
This example shows how to add your own FORTRAN functions to Ingrid. The first step is to create your own main program for Ingrid (this is unnecessary if you have already attached Ingrid to a model). The second step is to write a new "create" subroutine that adds the new words to Ingrid. The final step is to write and add the functions you want.
The main program for ingrid is rather simple. The only difference here is that we add the call to createMyWords
PROGRAM IngridMain
CHARACTER*80 PRGNAME
C creates various Ingrid modules, including your new one.
CALL createIngrid
CALL createHDF
CALL createEOS
CALL createMyWords
C Handles UNIX arguments
C It uses the first argument as an input file name if it exists,
C and gets its arguments from there. Otherwise, it reads from
C standard input
C FGFILE interprets the Ingrid command file
IF(IARGC().GE.1)THEN
CALL GETARG(1,PRGNAME)
OPEN(UNIT=1, STATUS='OLD', FORM='FORMATTED',
* FILE=PRGNAME)
CALL FGFILE(1)
ELSE
CALL FGFILE(5)
ENDIF
C DoTasks actually performs the tasks required by the Ingrid command file
CALL DoTasks
STOP
END
C Local Variables:
C compile-command: "f77 -o ingrid ingridmain.f -lingrid -lgrfps -lnetcdf
C -lsun -ldf -leos"
C End:
The following subroutine tells Ingrid what words correspond to the functions you wish to add.
All the definitions are added to the Ingrid: object, insuring
that they will always be available as long as the Ingrid:
object is on the stack.
SUBROUTINE createMyWords
EXTERNAL STREAMmyfn1, STREAMmyfn2
EXTERNAL sqrtsgn
CALL FGINTERP('Ingrid:')
CALL FGIoper('myfn1',STREAMmyfn1)
CALL FGIoper('myfn2',STREAMmyfn2)
CALL FGIscalarFilter('sqrtsgn',SQRTSGN,1,'/name /sqrtsgn cvx def')
CALL FGIpop
RETURN
END
Note that the EXTERNAL statements are very important.
They are required in FORTRAN in order to pass a function as an
argument to a subroutine, and if you leave it out, the program will
certainly core dump when you try to run it.
Now you need to write the FORTRAN code myfnstart1, etc. Note there are three sorts of functions you can add to Ingrid: an operator (which works on the stack directly), a scalarFilter (which operates on streams or real numbers, but it is limited to be a scalar function of its inputs), or a function (which operates in a more general way on streams.)
For an example of scalarFilter, I use a function already in Ingrid, sqrtsgn. Given a real FORTRAN function, FGIscalarFilter hooks it into Ingrid.
FUNCTION SQRTSGN(X) REAL SQRTSGN, X SQRTSGN = SIGN(SQRT(ABS(X)),X) RETURN END
See section FGIscalarFilter for more information on usingFGIscalarFilter.
See section Fortran Subroutines for more information on writing more general
stream filters.
While one usually filters one stream into another stream, it is important to be able to create a STREAM without an input stream. See section Fortran Subroutines for more information on writing general stream filters.
Given a subroutine
SUB(PARAM, T, OUT)
that outputs a block of data for each scalar T, we would like to connect it to Ingrid. We thus have the following FGIfunc call:
EXTERNAL SUB
...
CALL FGIfunc('sub',SUB,1,1)
i.e. the function has 1 input and 1 output.
This creates a word sub that has the following stack behavior:
( T PARAM -- )
i.e. it takes a grid T and a parameter block PARAM. We now need to make a
calls to NewBuffer for the output buffer. This stream
will correspond to the output array OUT of SUB.
So we have a FORTRAN program
SUBROUTINE SUB(p,T,OUT)
STRUCTURE /myS/
INTEGER NX
INTEGER NY
END STRUCTURE
RECORD /myS/ p
REAL T,OUT(p.nx,p.ny)
C we cleverly compute OUT for each T
...
RETURN
END
In our main program, after we call createIngrid, we call FGIfunc,
EXTERNAL SUB
...
CALL FGIfunc('mysub',SUB,1,1)
Now we create some Ingrid code to use our new function
/T /real 1. 1. 10. NewEvenGRID % we define the T grid
0 RECHUNK % we make it into a stream that feeds one
% element at a time into SUB
25 % gets the number of points from X
25 % gets the number of points from Y
2 integerarray astore % stores the two integers into an integer array
mysub % feeds the time and parameters to mysub [ nx ny ]
STREAM % all streams ultimately have the parent STREAM
% we are starting from scratch so we use STREAM
% directly .
25 25 mul % calculate the size of our output buffer
NewBuffer % create a new buffer which is filled by SUB
2 SetStreamIndex* % We now set the grids for the stream
% it is in 2D chunks
/X /real 1. 1. 25. NewEvenGRID % we define the X grid
/Y /real 1. 1. 25. NewEvenGRID % we define the Y grid
TaskStreams 0 get >T % gets the T grid from the input stream
* % ends SetStreamIndex
/name (mystream) def % a short name for our new stream
/fullname [ /MadeUP /Data ] def % a longer name for our new stream.
In practice, we present a clean interface to the outside world by writing a procedure, namely
mygreatsub
And we have the following definition
/mygreatsub
{
{ myX myY myT } inputs % creates object holding myX myY myT
myT 0 RECHUNK % gets time one element at a time
myX >npts myY >npts
2 integerarray astore % sets parameter block
mysub % feeds time and parameters to mysub
STREAM % parent of our output stream
TaskParameterBlock 0 get
TaskParameterBlock 1 get mul % calculates buffer size
NewBuffer
2 SetStreamIndex* myX myY myT * % sets the grids for the stream
/name (mystream) def % a short name for our new stream
/fullname [ /MadeUP /Data ] def % a longer name for our new stream.
(created by mygreatsub) addhistory
1output
} def
Then when we want to use it, we simply say
/X /real 1 1 25 NewEvenGRID /Y /real 1 1 25 NewEvenGRID /T /real 1 1 10 NewEvenGRID mygreatsub % created by mygreatsub X Y CONTOUR % a series of contour plots
If we want to use it to match the grids of the GLOBALHELLERMAN winds, for example, we could say
SOURCES GLOBALHELLERMAN taux X Y T mygreatsub X Y CONTOUR
We had a problem where we had two copies of a dataset -- SERVAIN and OLDSERVAIN -- and we wanted to make sure that OLDSERVAIN was contained in SERVAIN before we deleted it.
The first step was to take the difference between the two sets month by month. Since the data range is printed for each plot as it is made, we did not need to look at each plot.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES
SERVAIN OLDSERVAIN sub
Y low high RANGE % SERVAIN is stored high to low in Y
X Y COLOR % COLOR plots are the fastest
\end{ingrid}
As it turned out, the two data sets are not the same: starting in 1985, there are differences that are more than trivial. So now we plot the two different data sets and their difference.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES
SERVAIN
T 300 311 RANGE
Y low high RANGE
dup X Y CONTOUR
OLDSERVAIN
T 300 311 RANGE
Y low high RANGE
dup X Y CONTOUR
sub X Y CONTOUR
\end{ingrid}
Note that we use dup so that we can both plot the two streams
and use them in a subtraction.
Finally, we tried making a line plot at the spot in the map that had the largest difference in Jan 1985. It doesn't tell us much, however.
#! /usr/local/bin/ingrid
\begin{ingrid}
PS SOURCES
SERVAIN X -28 VALUE Y 18 VALUE T 276 311 RANGE dup T null LINE
OLDSERVAIN X -28 VALUE Y 18 VALUE T 276 311 RANGE dup T null LINE
sub T null LINE
\end{ingrid}
This is a FORTRAN code fragment that uses Ingrid to plot. It shows how to attach Ingrid to a FORTRAN program.
program coroort
c
c This program correlates the COADS data as reanalyzed by A. Oort
c of GFDL with a given input time series. The data are read from
c files of anomaly fields created by Steve Peng, residing on:
c /pahedra/poscratch/atlantic. They are piped from "zcat" since
c they are compressed, and the program reads them from standard input.
c Results are plotted out using INGRID
c
parameter (nx=46,ny=56,nxy=nx*ny)
parameter (alonw=-76, alone=14)
parameter (alats=-50, alatn=60)
parameter (nyrs=1987-1946+1)
parameter (zot=999.9)
c
real zdata(nxy) ! a single data grid
real xind(nyrs) ! index time series (one sample per year)
real absi(nyrs) ! time series of modulus of index
real zs(nxy) ! seasonal mean data
real zl(nxy) ! last season's data
real z1(nxy) ! first season's data
real za(nxy) ! ensemble mean data
real zr(nxy) ! rms height field
real zac1(nxy) ! 1-lag correlation field
real zsc(nxy) ! correlation <Z, index >
real zac(nxy) ! correlation <Z,|index|>
real alon(nx), alat(ny) ! arrays containing latitudes and latitudes
character*3 mname(12)
character*4 what,varnam
character*10 zname
character*24 string
character*48 period
logical ascii
logical BufferNeeded ! Ingrid-supplied function
data mname /'jan','feb','mar','apr','may','jun',
* 'jul','aug','sep','oct','nov','dec'/
...
c Plot data using INGRID subroutines
call createIngrid ! initialize ingrid
call SetIVAR1('X','longitude',nx,alonw,alone) ! define x-domain
call SetIVAR1('Y','latitude',ny,alats,alatn) ! define y-domain
call ModelDESC(period) ! set description
call FGINTERP('/missing_value 999.9 def') ! set missing value
call FGINTERP('X /fullname /longitude def pop')! fullname of X grid
call FGINTERP('Y /fullname /latitude def pop') ! fullname of Y grid
c Define the variables to be plotted
ID1=NewVAR('RMS','Rms value of '//varnam,'XY','XY')
ID2=NewVAR('R1AC','Lag-1 acr of '//varnam,'XY','XY')
ID3=NewVAR('SCOR','S X-corr '//varnam,'XY','XY')
ID4=NewVAR('ACOR','A X-corr '//varnam,'XY','XY')
ID5=NewVAR('SREG','S reg '//zname//varnam,'XY','XY')
ID6=NewVAR('AREG','A reg '//zname//varname,'XY','XY')
c Read INGRID instructions from input file
call FGFILE(7)
c Plot data
call FGPUT(ID1,zr)
call FGPUT(ID2,zac1)
if(BufferNeeded(ID3)then
do n=1,nxy
zdata(n)=zsc(n)
end do
what='COR'
call covxi(zdata,za,zr,nxy,aex,rex,ns,0.,what,.false.,zot)
endif
call FGPUT(ID3,zdata)
if(BufferNeeded(ID4)then
do n=1,nxy
zdata(n)=zac(n)
end do
call covxi(zdata,za,zr,nxy,aex,rex,ns,0.,what,.false.,zot)
endif
call FGPUT(ID4,zdata)
what='REG'
if(BufferNeeded(ID5)then
do n=1,nxy
zdata(n)=zsc(n)
end do
call covxi(zdata,za,zr,nxy,aex,rex,ns,0.,what,.true.,zot)
endif
call FGPUT(ID5,zdata)
if(BufferNeeded(ID6)then
do n=1,nxy
zdata(n)=zac(n)
end do
call covxi(zdata,za,zr,nxy,aex,rex,ns,0.,what,.true.,zot)
endif
call FGPUT(ID6,zdata)
c
stop
end
C Local Variables:
C compile-command: "f77 -o corOORTi corOORTi.f -lkushnir -lingrid -lnetcdf -lsun -lgrfps"
C End:
The BufferNeeded call is not required. BufferNeeded returns .FALSE.
if Ingrid has not been asked to do anything with the corresponding
buffer realization, allowing you to skip unnecessary calculations.
Note that FGPUT should be called even if BufferNeeded returns .FALSE.:
this lets Ingrid kept track of the current realization number for that
buffer (In this case each buffer has only one realization, so here it
makes no difference, but coding it this way is good practice).
This is an input file used with the program. Note that the program parameters are in the same file: Ingrid skips over everything until the \begin{ingrid} command.
47,86,1,4
WNTR SST1
1.2861868e+00, 7.1175910e-01, 1.3792630e+00, 1.3294256e+00,
1.1392760e+00, 1.7008809e+00, 1.5224471e+00, 1.2936962e+00,
1.2892774e+00, 1.2001539e+00, 9.9954061e-01, 3.0537025e-01,
6.2922901e-01, 1.4807101e+00, 1.5393012e+00, 7.6679767e-01,
3.3163543e-01, -4.8041282e-02, -5.8387122e-01, -8.8469503e-01,
-5.2401093e-01, -6.4082654e-01, -8.7956411e-01, -1.4008650e+00,
-7.4191162e-01, -1.1428326e+00, -8.0804998e-01, -5.8146370e-01,
-8.4327691e-01, -9.5064736e-01, -9.8119418e-01, -9.0490550e-01,
-1.0669201e+00, -7.7010233e-01, -7.1211606e-01, -7.5992757e-01,
-1.0772113e+00, -7.7296900e-01, -2.2595750e-01, -9.4640652e-01
\begin{ingrid}
PS SOURCES MODEL
(corOORT) setplotname
RMS
Y -40 60 RANGE
X Y CONTOUR
R1AC
Y -40 60 RANGE
X Y CONTOUR
SCOR
Y -40 60 RANGE
X Y CONTOUR
ACOR
Y -40 60 RANGE
X Y CONTOUR
SREG
Y -40 60 RANGE
X Y CONTOUR
AREG
Y -40 60 RANGE
X Y CONTOUR
\end{ingrid}
The plotting routines allow you to use a function of grid instead of the
grid itself. This allows you to make log plots, for example, by simply
taking the log of the grid before calling COLOR or
CONTOUR.
In the case of log plots, it is important not to feed log any negative numbers.
#! /usr/local/bin/ingrid
Tests of color map uneven gridding (log Z axis)
\begin{ingrid}
PS (color) setplotname
SOURCES LEVITUSMEAN
temp X 0 VALUE
Z high low RANGE
Y Z 100 add log COLOR
\end{ingrid}
Local Variables:
compile-command: "/usr/local/bin/ingrid < testcolor.in"
End:
This is an example of reading and regridding a dataset to use as model input. This was written by Zhao Yuechen, expanded comments by Benno.
First the total example,
PARAMETER(NX=20, NY=10)
real taux(NX,NY)
xs=110
xe=160
ys=-10
ye=10
NT = 24
ts = .5
te = 11.5
call createIngrid
call SetIVAR1('X','longitude',nx,xs,xe)
call SetIVAR1('Y', 'latitude',ny,ys,ye)
call SetIVAR1('T', 'monthtime',NT,ts,te)
idtaux=NewVar('taux', 'taux', 'XYT','XY')
call fginterp('SOURCES')
idtauxin=idstream('GLOBALHELLERMAN taux MODEL >taux gridtomatch ')
DO IT = 1 , NT
call fgget(idtauxin,taux)
write(*,*)'opo',taux
END DO
end
Now the same example with more explanation. At the top we declare our array and set some grid limits.
PARAMETER(NX=20, NY=10) real taux(NX,NY) xs=110 xe=160 ys=-10 ye=10 NT = 24 ts = .5 te = 11.5
The first step is to initialize Ingrid. This leaves three object open,
i.e. there are three sets of words available for use: Ingrid:
which contains the basic Ingrid words, SOURCES, which contains
all the entries in the data catalog, and MODEL, which will
contain any variables that you define. (It will contain those
definitions because it is the most recently opened object).
call createIngrid
The next step is to define the grids that you want to interpolate to. 'XYZT' are the best choices for grid names, because those are the grid names used in the data catalog (upper case, please).
call SetIVAR1('X','longitude',nx,xs,xe)
call SetIVAR1('Y', 'latitude',ny,ys,ye)
call SetIVAR1('T', 'monthtime',NT,ts,te)
The next step is to define variables. In particular, we are going to
define a variable taux which has the grids and blocking that we
want the data to have. We say that taux is a function of 'XYT'
with blocking 'XY', by which we mean that it is a function of X, Y, and
T with X the fastest varying coordinate, then Y, and finally T.
Furthermore, we want the data in blocks of 'XY', i.e. each time we
request data from Ingrid (using the FGGET call), we expect to get a set
of XY values, each call corresponding to a different time.
idtaux=NewVar('taux', 'taux', 'XYT','XY')
Ingrid has an object SOURCES which contains all the datasets in
the data catalog. It also contains MODEL, which contains all the
variable definitions made above; in this case, just taux. Here
we open the data catalog.
call fginterp('SOURCES')
Now we want to get an id for regridded taux. In particular, we want the GLOBALHELLERMAN taux data regridded to match the MODEL taux grids defined above.
idtauxin=idstream('GLOBALHELLERMAN taux MODEL >taux gridtomatch ')
having gotten our id, we can use it to extract the taux data
DO IT = 1 , NT call fgget(idtauxin,taux) write(*,*)'opo',taux END DO end
Go to the first, previous, next, last section, table of contents.