Go to the previous, next section.

Examples with discussion

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.

Simple plot of SST data from the data catalog

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.

Average equatorial SST anomaly plot

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
( -- stream ) returns stream of CAC SST anomalies.
RANGE
( stream grid lo hi -- stream' ) modifies stream to have the range [lo,hi] in coordinate grid. In the first plot it restricts Y to [-5,5].
AVERAGE
( stream grid -- stream' ) averages the data in stream along coordinate grid. In this case it averages along Y.
CONTOUR
( stream gridh gridv -- ) generates contour plots, using 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.

Plot Appearance

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}

Oceanic equation of state

#! /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.

Preparing a netCDF File

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.

Describing the dataset in CDF

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.

Running ncgen

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

Modifying the ncgenerated program

There 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:

Linking

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!

Color Plot

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}

Adding FORTRAN fns

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.

Writing a new main program for Ingrid

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:

createMyWords

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.

Creating a STREAM

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
( X Y T -- stream' ) creates a new stream from three grids.
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

Comparing datasets

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}

Attaching Ingrid to a FORTRAN program

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}

Log axes

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:

Reading Data into a model

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 previous, next section.