NMC GrADS User-Defined Function #1
                                                      
			     REGRID
			 Version 2.0beta
                                                      

			  Mike Fiorino
                                                      
		    University of California
	     Lawrence Livermore National Laboratory
     Program for Climate Model Diagnosis and Intercomparison
		       P.O. Box 808 L-264
		       Livermore, CA 94551

		    fiorino@typhoon.llnl.gov

	     510-423-8505 (voice) 510-422-7675(fax)
                                                      
			17 January, 1995
				

PREFACE:

Regrid, the first in a series of Grid Analysis and Display System
(GrADS) user-defined functions developed originally for the
Development Division of the National Meteorological Center, has
been substantially improved and now features a simpler calling
syntax.  This is a BETA version, meaning I have tested it a lot
but have not exposed its warts to the world.  It is YOUR job to
BREAK it so I can correct it...

The original version can be found by

anon ftp to climate.gsfc.nasa.gov
cd /pub/fiorino
bin
mget *regrid*

The installation calls the routine regrid2beta.  This allows for
intercomparison with the original regrid and transition to the
new version.  When you are satified with its performance, then
you can edit the GrADS user-defined function table and replace
the regrid binary.

This document only assumes you have a basic knowledge of GrADS
and does not require experience with the original version regrid
v1.0.

regrid has NOT been tested on PCs; only Suns and SGIs...

GrADS USER DEFINED FUNCTIONS:

A user-defined function works like any other GrADS function
(e.g., hcurl(u,v) ), but is external to GrADS itself.  In
essence, the user-defined function capability is an interface
specification between GrADS and codes developed by users in
potentially other languages (e.g., C and FORTRAN).

The regrid user-defined function solves a common problem of
transforming HORIZONTAL 2-D gridded fields from/to different
resolutions/grid types for quantitative intercomparison.  For
example, a model monthly mean precipitation field on a T126
gaussian grid can be compared to an observed climatology on a
2.5x2.5 grid using regrid.  regrid offers many transform options
ranging from simple bilinear interpolation to box averaging with
"voting."  Additional methods can be added to regrid as needed.

Your comments on the format of this documentation and regrid
itself are greatly desired.  Furthermore, I will be maintaining a
list of users to contact when making bug fixes and upgrades.
Please direct your comments by email.

DESCRIPTION:

regrid transforms two-dimensional (2-D) lat/lon GrADS grids from
one grid type/resolution to another.  The input is any 2-D
lat/lon grid defined by the current GrADS lat/lon dimension
environment.  regrid handles input grids which are cyclically
continuous in longitude and excludes undefined input grid values
from participation in the transform.  If a valid transform cannot
be made (i.e., insufficient defined data points), the output grid
is set to undefined.  regrid supports two output grid types: 1)
lat/lon with fixed increments; and 2) gaussian.  Four transforms
are available: 1) box averaging for regridding fine to coarse
grids; 2) box averaging with "voting" for noncontinuous/index
data such, as soil type; 3) bilinear interpolation; and 4)
4-point bessel interpolation.

INSTALLATION AT NMC:

1)    Make sure you are running GrADS Version 1.4 or higher and that you start
      GrADS from a directory where you have write permission!

2)    Change your environment:

      For Suns:  setenv GAUDFT /export/nmc01/wd20mf/grads/udf/sun/udft
      For SGIs:  setenv GAUDFT /export/nmc01/wd20mf/grads/udf/sgi/udft   

3)    File Locations:

      docs (regrid.doc regrid.wp) and source (*.f) are in

      /export/nmc01/wd20mf/grads/udf/regrid

INSTALLATION OUTSIDE NMC:

Refer to the readme.regrid2 in regrid2beta.tar.

USAGE:

regrid(G1,X1,X2,C1,X3,X4,X5,X6)

G? = a GrADS grid expression (e.g., z or ave(z.3(t+0,t=120,1yr))
X? = a float parameter
C? = a character string parameters

REQUIRED params:  G1, X1

where:

   G1     =  any valid GrADS grid expression  

   X1     =  delta longitude (dlon) or # of gaussian longitudes on the
	     GLOBE

   X2     =  delta latitude (dlat) or # of gaussian latitudes on the GLOBE  

   C1     =  character option string
 
   X3     =  special case options #1

   X4     =  special case options #2

   X5 	  =  special case options #3

   X6     =  special case options #4

OPTION STRING:

C1 consists of a sequence of up to four options separated by the
underscore (_).  They are:

ba	-	box averaging (the default does NOT have to be set)
vt	-	"vote" interpolation or box averaging with voting
bl	- 	bi-linear interpolation
bs	-	3rd order Bessel interpolation

igXXX   -	input grid is gaussian with XXX the number of
		gaussian latitudes (e.g., ig92 for the NMC T62
		grid).  XXX must be >= 8 and a multiple of four.
		This param is used to invoke a more precise
		calculation of the boundary between gaussian grid boxes

maXXX  	- 

minimum % (XXX) area which must be covered with DEFINED
data boxes to be considered a valid interpolation.  Applies ONLY
to box averaging without voting when the input grid has undefined
points.  regrid v1.0 assumed am was 0 or that if ANY input grid
boxes contained defined data which intersected the output grid
produced a box average.  This was clearly too liberal and ma is
now set by default to 50% or that half the output grid box must
be covered with defined data to produced a defined output grid
point

gg - 

tells regrid that X1 and X2 are the number of GLOBAL gaussian
longitudes and latitudes on the output grid, even if the output
grid is a non global.  X1 and X2 define the gaussian grid.

p1 	-

special option which forces the (1,1) point of the output grid to
be anchored to a specific lon,lat defined by the pair (X3,X4) or
(X5,X6).  See special cases below:

SPECIAL CASES:

There are two pairs of special case params: (X3,X4) and (X5,X6).
The ordering depends on the order that the options are specified
in C1.  The example will make it clear how this works

Case 1:		Box averaging with voting (C1 = vt) 

	(X3,X4) or (X5,X6) set the fraction (called VT1 and VT2)
of an output grid box which must be covered by defined input grid
data for a "winner" to be chosen in the election.  The default is
VT1=VT@=0.5.

	VT1 = 0.0 - 1.0; minimum fraction of an output grid
point when there is only one candidate

   	VT2 = 0.0 - 1.0; same as VT1 except for three or more
candidates.  The fraction for two candidates is midway between X4
and X5.
   
   X4, X5 = 1.0 would require that 100% of the output grid box
must be covered by a single, unique value from the input grid
whereas X4,X5 = 0.0 would allow a winner to be chosen if ANY
candidates where running.  The default value of 0.5 means that a
simple majority is required to have a winner.
   
Case 2:      	Set the beginning lat/lon (the (1,1) point) of
		the output grid by overriding the automatic
		lat/lon bounds setting of regrid for LON-LAT
		output grids
   
   P1 = beginning longitude (lower left hand corner or the point
(1,1) on the grid) of the output grid

   P2 = beginning latitude of the output grid

          The number of points in longitude is then calculated
as:
 
   (eastward-most longitude of the input grid - X4)/abs(X1)+1
    
   and for latitude:

   (northward-most latitude of the input grid - X5)/abs(X2)+1
  
See the METHODOLOGY section for a description of the transforms  

EXAMPLES:

1) Regrid a global T62 gaussian grid (192x94) to a 2.5 deg lat/lon by box
   averaging,
 
   open /reanl1/pilot20/fluxgrb8508.ctl
   set x 1 192
   set y 1 94
   define p25=regrid(p,2.5,2.5,1)

Note: The lat/lon dimension environment is set using grid
coordinates (x,y) to make the input and output grids global. To
minimize future graphics/analysis calculations with the regridded
field, we use the GrADS define function to store the grid in
memory where it can be referenced as any other GrADS grid.

2) Regrid a 4x5 SiB vegetation type to a R15 (48x40) gaussian
grid using box averaging with "voting."  Require that at least
60% of the output grid box must be covered with a single
candidate value from the input grid for an election to occur.
Otherwise the output grid box is set to undefined.  Relax the
one-candidate election requirement to 20% when there are three or
more candidates,
 
   open /export/sgi18/fiorino/data/sib/sib.param.annual.ctl  
  set lon 0 360
  set lat -90 90
  define i21=regrid(index,48,40,12,0.60,0.20)
  set gxout grfill
  d index
  d i21

Note: During candidate selection, undefined input grid points do
not contribute to the fraction of the output grid box covered
with input grid boxes.  The best way to display index type data
in GrADS is to use the "grid fill" display option (set gxout
grfill).  GrADS will draw the grid box and color it according to
the grid value.

3) Regrid 1x1 Aviation run 500 mb z to 2.5x2.5 grid using bessel interpolation,
 
   open /export/sgi39/wd22sl/grads/avn/avn.93092800.fcst.ctl  
  set lev 500
  set lat 20 70
  set lon -140 -40
  d regrid(z,2.5,2.5,4)

Note: Box averaging would be more appropriate when regridding to
a coarser grid.

4) Regrid 1x1 Aviation run 500 mb z to 2.5x2.5 grid using box averaging and
   setting the    grid to start at a specific lat/lon,
   
   open /export/sgi39/wd22sl/grads/avn/avn.93092800.fcst.ctl  
  set lev 500
  set lat 20 70
  set lon -140 -40
  d regrid(z,-2.5,-2.5,1,-141.0,21.0)

Note: the -2.5 for dlat/dlon tells regrid that the fourth and
fifth parameters will be setting the beginning lat/lon of the
output grid.  The grid spacing of the output grid is abs(-2.5)
deg.

RESTRICTIONS:

1) regrid currently is limited to input/output grids of dimension
730x380 (~T225), but this can be easily changed by modifying the
nimax,njmax parameters in FORTRAN source in regrid2beta.f

2) Any valid GrADS grid can be regridded.  However, GrADS (V1.5)
currently only supports lat/lon grids where the mapping between
grid and world coordinates is one dimensional, i.e.,
longitude(i,j)=f(i) vice longitude(i,j)=f(i,j).
 
3) Only two output grid types have been implemented: 1) lat/lon
(dlat does NOT have to equal dlon); and 2) gaussian grids.
Mercator output grids could be added as lon(i,j)=f(i) and
lat(i,j)=f(j) in this projection.

4) The lat/lon bounds of the output grid are calculated within
regrid to encompass as much of the input grid as possible.  For
uniform lat/lon output grids, the number of grid points and the
starting lat/lon are based on fitting the input grid within a
"template" global grid with the following properties: 1) a grid
spacing specified by X1 (dlon) and X2 (dlat); 2) starts at 0W
with a grid point on the equator; and 3) is cyclically continuous
in longitude.  The lat/lon bound setting algorithm can be
overridden by setting X1 and/or X2 to the negative of the
required output grid spacing.

SUPPORT:

The FORTRAN source code is available upon request and you are
free to use the regrid routines in other programs or to modify
the algorithms.  However, I would prefer that you send bug
reports, comments and suggestions for improvements directly to me
(email is best) so that I can maintain a single "standard"
version.

METHODOLOGY:

The first step in the regrid transform is to define a
relationship between the input and output grids within a common
frame of reference or coordinate system.  regrid bases the
inter-grid relationship on "world" coordinates, and the GrADS map
between grid coordinates (i,j) and world coordinates (lon, lat).
As noted above, the world-grid mapping in GrADS is
one-dimensional.  Thus, the world-grid map of an input GrADS grid
takes the form,

   lat(i,j)=f(j) and lon(i,j)=g(i).

By specifying a similar mapping for an output GrADS grid of the
form
 
   LAT(I,J)=F(J) and LON(I,J)=G(I), 

as defined by the input parameters X1, X2 and X3-6, regrid
calculates,

   X(I)=i(G(I)) and Y(J)=j(F(J)),

where i(G(I)) is the location of the output grid with respect to
the input grid dimension i and j(F(J)) for j.

For simplicity, and greater generality, regrid assumes that the
grid point is at the center of a rectangular grid box and maps
the location of the boundaries of the output grid box to that of
the input grid box.  By default the boundaries are assumed to lie
midway between grid points and while this is not strictly true
for a gaussian grid near the poles, it is close nonetheless.  The
boundaries for gaussian grids can be calculated by specifying
igXXX in C1.  The reason why this cannot be automatic is that
GrADS does not directly support gaussian grids (i.e., there is no
ydef gauss 40 option in the data descriptor .ctl file, just
linear and levels).

Given the inter-grid map X(I) and Y(J), regrid uses two basic
methods for doing the transform: 1) box averaging; or 2)
interpolation.  Box averaging is simply the area-weighted
integral of all input grid boxes which intersect an output grid
box, divided by the area of the output grid box.  This approach
is most appropriate: 1) for transforming from fine (e.g., dlon =
1 deg) to coarse grids (e.g., dlon = 2.5 deg); and 2) when
approximate conservation of an integral quantity (e.g., global
average) is desired.

Box averaging is also useful for regridding noncontinuous,
parametric or "index" data.  For example, suppose you have been
given a 0.5x0.5 deg global grid of vegetation type and want to
use these data in an R43 global model.  The intuitive solution is
to assign the output grid the value of the intersecting input
grid box(es) which account(s) for the greatest percentage of the
output grid box surface area.  In the example of vegetation data,
if 70% of the output grid box is covered by deciduous forest,
then it might be reasonable to call the output grid deciduous
forest.  However, if there were 5 distinct vegetation types or
"candidates" available, then regrid, being an American function,
holds an "election" and select a "winner" based on the candidate
covering the greatest percentage of the total surface area in the
output grid box.  Of course, coming from an imperfect democracy,
the election can be "rigged" for a desired effect....

This grid transform strategy is invoked using the "vote" option
in box averaging (vt in C1).  Conditions on the percentage of the
output grid box (number of candidates and what it takes to get
elected) can be finely controlled by the X4 and X5 parameters.

Perhaps the most conventional way of regridding meteorological
data (e.g., 500 mb geopotential heights) is interpolation because
weather data tend to be continuous .  regrid features a 4x4 point
bessel interpolation routine developed at Fleet Numerical
Meteorology and Oceanography Center (courtesy of D. Hensen,
FNMOC).  While this routine is in wide use at FNMOC, the regrid
implementation has been substantially improved to handle more
general input grids.

First, bilinear interpolation is performed at all points to
produce a "first guess."  Improvements to the bilinear "first
guess" are made using the higher-order terms in the bessel
interpolator, but only if the bessel option is set (i.e., bs in
C1).  Second, an undefined value in the 2x2 bilinear stencil
causes the output grid to be undefined.  If the bilinear value is
defined, but any of the points in the larger 4x4 bessel stencil
are undefined, the output grid is assigned the bilinear value.
The third improvement is that the routine handles grids which are
cyclically continuous in longitude.

It is generally recommended that you use the bessel option when
interpolating because the higher-order terms in the polynomial
interpolation function produce a closer fit to the original data
in regions of rapid changes (e.g., large gradients of the
gradient around low pressure centers).

ADDITIONAL NOTES:

1) Where to run GrADS

There is normally no problem firing up GrADS from ANY directory
because GrADS typically only READS data, but with user-defined
functions and regrid, you'll need to start GrADS from a directory
where you have write permission.  The reason is that when a
user-defined function is invoked within GrADS, a file with the
grid data is written out by GrADS and then a file written by your
function is read back into GrADS for further processing.
Although this restriction could be modified by changing the code
and the user-defined function table, the default is for regrid to
read/write from the directory where GrADS was first invoked.  You
will notice when exiting GrADS that regrid leaves behind three
files: 1) udf.regrid.gfi (input to regrid from GrADS); 2)
udf.regrid.gfo (output from regrid returned to GrADS); and 3)
udf.regrid.out (diagnostic print file).  Use the GrADS shell
escape in a script to remove these regrid "droppings" (e.g., !rm
udf.regrid.*).  Future user-defined functions will the prefix
"udf." so that the rm command will take out all droppings from
other user-defined functions.

2) Using regridded fields in other GrADS functions

The only down side to a regridded field is that its dimension
environment cannot be controlled by its "grid" coordinate system.
The best way to understand this is by an example. Suppose you
regrid a T62 global Gaussian grid (192x94) to a uniform 2.5 deg
grid (144x73) using box averaging and the GrADS define
capability, e.g.,

define p25=regrid(p,2.5,2.5,ba) 

You now want to calculate the global average of the original
field p and the defined regridded field p25.  The only
unambiguous way (using all grid boxes) of doing this calculation
for p would be,

d aave(p,x=1,x=192,y=1,y=94)

and not,

d aave(p,lon=0,lon=360,lat=-90,lat=90)
 
This is because the cyclic continuity feature in GrADS would
count grid boxes at the prime meridian twice, i.e., GrADS would
really be doing,

d aave(p,x=1,x=193,y=1,y=94)

Trying to find the global average of the 2.5 deg regridded field
p25 using,

d aave(p25,x=1,x=144,y=1,y=73)

would not yield a global average even though p25 IS 144x73!
However,

d aave(p25,x=1,x=192,y=1,y=94)

would because GrADS converts the grid coordinate range to
(x=1,x=192) to world coordinates (lon=0,lon=360-1.875) and grabs
all grid boxes in p25 within that range when putting together the
data for the areal averaging calculation.  Despite this
restriction on grid coordinates, you can access specific chunks
of a regridded defined variable using world coordinates.  Suppose
you want to look at the latitudinal variation of the u wind
component at 180 deg and 500 mb on a 1 deg grid, e.g.,

set lev 500
set lon 180
set lat -90 90
d u

if the you had,

define u5=regrid(u,5)

you could then,

d u5(lon=175)

but not,

d u5(x=1)

3) Diagnostic messages from regrid

regrid sends information on the transform process
(starting/ending lat/lon, number of grid points and the
regridding method) to the terminal window where you are running
GrADS.  Additionally, errors during the call to regrid (e.g.,
trying to regrid a two-dimensional longitude-pressure cross
section) will be displayed, the process terminated, and control
returned back to GrADS.

Future releases of regrid will allow this diagnostic information
to be passed back to a script for additional control of GrADS
processing.  However, this requires a change to the specification
of the user-defined function interface and GrADS itself.