Schemes and grids for the PRISM interpolation library

A first list of  2D interpolation schemes to be implemented in the PRISM coupler are the following:

1. S1 - Nearest-neighbour interpolation function:

For each target grid point, the n nearest neighbours on the source grid, weighted or not by their distance, are averaged.
The weight of the different neighbours can be uniform, given by their distance from the target point, or given by the value of a gaussian function at their distance from the target point.
 
2. S2 - Nearest-neighbour gaussian weighted interpolation function:

For each target grid point, the n nearest neighbours on the source grid, weighted by the value of a gaussian function at their distance from the target point are averaged.
 
3. S3a - Bilinear (1st order) interpolation function:

Standard bilinear interpolation which involves 4 neighbours.
 
4. S3b - Bilinear (1st order) extrapolation function:

Standard bilinear extrapolation which involves 4 neighbours.
 
5. S4a - Bicubic (2nd order) interpolation function:

Standard bicubic interpolation which involves respectively 16 neighbours.
 
6. S4b- Bicubic (2nd order) extrapolation function:

Standard bicubic extrapolation which involves respectively 16 neighbours.
 
7. S5 - First order conservative remapping

This scheme will guarantee that the area-integrated field (e.g. water or heat flux) is conserved between the source and the target grid. For each target mesh, the weight assigned to each underlying source mesh is be proportional to the overlapped surface.
 
8. S6 - Second order conservative remapping

This scheme will also guarantee that the area-integrated field (e.g. water or heat flux) is conserved between the source and the target grid.  For each source grid point underlying the target mesh, a weight is given to the value of the source field at the source grid point and also to the value of the derivatives of the source field in both directions (see SCRIP documentation p.16). Second-order and higher order conservativ remapping are also proposed theoretically by Vintzileos and Sadourny, 1995.

The following grids should be supported for the above scheme. These grids have the following common characteristics:

• The grids may have masked grid points.
• The grids may have "holes" (i.e. they do not cover the whole sphere, e.g. regional grids).
• The grids may be global or regional (except the G3 reduced grid which is always global).
• The grids may have overlapping grid points (except the G3 reduced grid which is always global)..

1. G1 - Lat-lon grids:

The grid is given by the intersection of meridians and parallels. Each (i,j) grid point can be described by the value for the indices i and j of two 1-D arrays, latitude(j) and longitude(i).
• The latitudinal mesh sizes can be regular or irregular, i.e. latitude(j) is or not constant.
• The longitudinal mesh size can be regular or irregular,  i.e. longitude(i) is or not constant.
• The grid may overlap in longitude with N overlapping grid points.
• The grid may have a grid points at the pole and/or at the equator.
• The border of the meshes can be approximately placed half-way between two grid points, but the exact location of the border meshes should be given with the grid.
2.  G2 - Cartesian and streched and/or rotated grids:

Each (i,j) grid point can be described by the value for the indices i and j of two 2-D arrays, latitude(i,j) and longitude(i,j).
• The grid may overlap itself in the i and/or j direction.
• The exact location of the border meshes should be given with the grid.
3. G3 - Reduced grids

The grid is composed of a certain number of latitude circles, each one being divided into a varying number of longitudinal segments. The grid can be described by the number of latitude circles, N_lat, the latitudinal position of each circle npos(N_lat) and by an 1-D array giving the number of longitudinal segments for each latitude, n_seg(N_lat). The total number of grid points, N_tot, is the sum of all elements of n_seg. The grid can also be described by two 1-D arrays, latitude(N_tot) and longitude(N_tot).
• There is no overlap of the grid.
• There is no grid point at the equator nor at the poles.
• There are grid points on the Greenwich meridian.
• The exact location of the border meshes should be given with the grid.
4. G4 - Unstructured grids

The grid, with N_tot number of grid points, has no logical structure.  The grid must be described by two 1-D arrays, latitude(N_tot) and longitude(N_tot)
• The grid may overlap itself for some grid points.
• There may be a grid point at the equator or at the poles.
• There may be grid points on the Greenwich meridian.
• The mesh can have an arbitrary number of sides. The exact location of the border meshes should be given with the grid.

How to proceed
 

• Confirm that the optimal way of performing the interpolation is to calculate once at the beginning of the simulation the matrices of weights and addresses, keep these matrices in memory for the whole simulation, and use them for each interpolation in a matrix-vector product.
• Determine the best structure of the weight-and-address matrices for optimal efficiency in matrix-vector product.
• Determine the different interpolation schemes for the different grids. The sources to be evaluated are:
• Present oasis libraries (but not fscint)
• New fscint of RPN
• SCRIP library
• Renka's schemes (see articles in http://wwww.acm.org/calgo/)
• Optimize the different interpolation schemes for scalar fields the different grids:

Order of priority is:
• G1-S1, G1-S2, G1-S3a, G1-S4a, G1-S5,
• G2-S1, G2-S2, G2-S3a, G2-S4a, G2-S5,
• G3-S1, G3-S2, G3-S5,
• G4-S1, G4-S2, G4-S5,
• G1-S3b, G1-S4b, G2-S3b, G2-S4b,
• G3-S3a, G3-S4a, G3-S3b, G3-S4b,
• G1-S6, G2-S6, G3-S6, G4-S6
• G4-S3a, G3-S4a, G3-S3b, G3-S4b.
• Detremine and optimize treatment of vector fields.
• Parallelize the different interpolation schemes for the different grids.
• Keep contact with RPN and with SCRIP authors...

1. Should we be general and consider 1D and 3D interpolation at the same time?