Go to the previous, next section.
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.001or
sendsumry ltest.001sendsumry 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.001psplot.xps is the preamble file for Ingrid plots.
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.PS and SOURCES simply put objects on the stack so
that the needed definitions would be available.
#! /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.
(@xref{top,netCDF,Introduction,/usr/local/emacs/info/netcdf,netCDF Use
}'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.
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.
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 programcorner(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:
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 CONTOURand not:
(fssta.cdf) readCDF ssta X Y CONTOUR(Alternatively, you could call your grids X and Y).
Good luck!
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}
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:
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
/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
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}
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}
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:
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 previous, next section.