First, let’s explain briefly how the package works. The numerical models used in RPN are mapped on a handful of map projections: Gaussian grid (SEF), Polar Stereographic (MC2), and rotated lat-lon (GEM). This mapping is expressed as parameters such as the size of the grid in points (ni x nj), its resolution (in km or degrees of latitude/longitude), its orientation (relative to a reference meridian). Many grids are regularly spaced, but some are not: their grid length may vary along x and y direction.

For all sorts of use, we need to transform the output (or input) of these models to other map proJections. This mapping is actually done through lagrangian interpolation. Here is the basic series of steps one has to take to go through a grid interpolation.

• Define the input and output grids
• Compute the lat-lon coordinates of every point of the output grid
• From these lat-lon coordinates, compute the x-y position of these points on the input grid.
• Use these x-y positions to interpolate the values from the input grid using bi-cubic lagrangian interpolation.

Although this seems simple, the whole process can be complicated by things such as grid symmetry, irregular grid spacing, non-existence of data at the poles, and rotation of wind components along grids.

The cost of computation is split in 2 parts: the computation of the geographical mapping between the 2 grids, and the lagrangrian interpolation. The cost of the first is much higher than the cost of the second, although in theory it needs to be only done once. Once the mapping is established, it can be used over and over to interpolate many meteorological fields.

As mentionned, scint was the first implementation of the grid interpolators used in RPN utilities like PGSM. The grid types supported by scint were basically polar stereographic (‘N’ and ‘S’), latlon (‘L’, ‘A’ and ‘B’) and gaussian (‘G’). The package was scalar (in the vector computer sense) in its very essence: the interpolation routine scint2 had to be called for every point. Moreover, the computation of a position of a lat-lon point on a grid had to be computed for every point.

This was highly inefficient on a vector computer like the CRAY and the SX-4. This is why the package was re-written in 1991 to get optimized for vector computers. Some bookkeeping of the grid parameters was done to keep in memory the coordinates of the input and output grids, allowing it to reuse the transformed coordinates if no change in the input and output grids was detected. This, in addition of the vectorization of the interpolation process, led to a 6-time increase in speed of utilities like PGSM.

fscint solved most of the performance problems of scint, but some quirks, mostly having todo with the programming interface remained. Among those, mention:

• The user had to decide which routine to call: rgscint (if the grid type was "regular"), or igscint (if the grid type was "irregular", that is, a ‘Z’ grid).
• The user had to calculate the lat-lon coordinates of the output grid.
• The conversion routines between lat-lon and x-y coordinates were not very well documented.
• The user had to do all the management of the deformation axes used in the ‘Z’ grids (the famous ">>" and "^^" records, informally knowned in the building as tictic-tactacs), including how to look for them in an RPN standard file, how to read them while avoiding to get mixed up in the numerous IGs, and to understand thoroughly how this was holding together (not so obvious - many found there a steep learning curve).
• The wind interpolation had to be done in two parts: from the x-y components of the source grid to meteorological wind speed and direction, and from wind speed and direction to x-y components on the target grid.
• The package had destructive side-effects for some grid types. For example the source field was turned upside down for data fields stored in ‘A’, ‘B’, and ‘G’ grids stored north to south instead of south to north.
• The package was not very portable, relying heavily upon the POINTER statement and HPALLOC routine in the FORTRAN code.

As the use of the RPN standard file increased steadily in the canadian scientific community inthe 90’s, so did the need for a more user friendly package to be developed. The arriva1 of the GEM mode1 and its rotated lat-lon grids, very difficult to analyse visually on its native grids, spawned the need to extend the functionality of the package to be used in utilities like xrec and max. The arriva1 of coupled models, having to interpolate data steadily back and forth between atmospheric and oceanic models also fueled this need. And most of all, the complaints of many users saying the package and its mechanics were too complicated.

So here comes ezscint. Here is a list of the intended benefits of the new package.

• Black-box approach. The package will take care of the ‘Z’ grids by itself, if asked to. All the user needs to supply is the NI, NJ, GRTYP and IGl to IG4 parameters returned by FSTPRM.
• The package keeps up to 32 different grids and 64 grid sets (a grid set is composed of output grid) in memory, minimizing the cost of navigational computing. If this number is exceeded, the oldest grid parameters are overwritten.
• Once a grid set is defined (output/input grids), the interpolation calls are trivial:
  • ezsint(zout,zin) for scalar fields,
  • ezuvint(uuout,vvout,uuin,vvin) for vector fields.
• The interpolated winds are returned in grid or meteorological components.
• The lat-lon coordinates of any defined grid are returned by a single call
• The routines gdxyfll and gdllfxy are taking care of the conversion between lat-lon and grid coordinates.Values at given lat-lon or x-y points can be automatically interpolated.
• The package is entirely usable from C, in which most of the code is written. FORTRAN remains used where it is most efficient, in the core part of the interpolation. Memory allocation is all done in C, symplifying the FORTRAN part, and implying that a port to FORTRANless environment (Linux, Wintel) with f2c translation should be easier (if not very efficient).

General Use

To use the package, there is at least one routine that needs to be called, it is either ezqkdef or ezgdef_fmem . This defines at least one grid on which some operations are possible. When only this function is invoked, only 3 functions (albeit the most important) cannot be called. These are ezsint, ezuvint and ezwdint. This is typically used when one user wants to do interpolation on selected lat-lon points, or transformation of wind components on a grid.

When used to do regular grid interpolation, ezgdef or ezqkdef need to be called with the parameters of the target grid, and 2 more routines need to be called : ezdefset (to define the grid set) and either ezsint, ezuvint or ezwdint to do the interpolation.

Suppose we want to interpolate forecasted fields from the GEM grid to a Gaussian Grid.

Suppose all the parameters are declared somewhere above the following piece of code.

Suppose variables and fields ending with « gem » are defined on the GEM grid, and those with « gauss » are defined on the Gaussian Grid.

Here is the code that will do it.

gdgem   = ezqkdef(nigem, njgem, ‘Z’, ig1gem, ig2gem, ig3gem, ig4gem, iunit)
gdgauss = ezqkdef(nigauss, njgauss, ‘G’, ig1gauss, ig2gauss, ig3gauss, ig4gauss, -1)
iset    = ezdefset(gdgauss, gdgem)
ier     = ezsint(zgauss, zgem)
ier     = ezsint(ttgauss, ttgem)
ier     = ezuvint(uugauss, vvgauss, uugem, vvgem)

Now suppose that we want to interpolate output from the GEM model to selected lat-lon points, say, Montreal, Toronto and Vancouver. Suppose we the same assumptions as above. So here is the code

lat(1)  =   45.73
lon(1)  =  -73.75
lat(2)  =   43.40
lon(2)  =  -79.38
lat(3)  =   49.18 
lon(3)  = -123.18
gdgem   = ezqkdef(nigem, njgem, ‘Z’, ig1gem, ig2gem, ig3gem, ig4gem, iunit)
ier     = gdllsval(gdid, zout, zgem, lat, lon, 3)
ier     = gdllvval(gdid, uuout, vvout, uugem, vvgem, lat, lon, n)

There is a more complete example showing the difference between the fscint and ezscint packages.

Go to the list of functions, or the general index