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 defI 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
ncgen
erated 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 doIn 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 ENDIn 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 } defThen 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.