Integration of Pattyn 2003 Model
Because Pattyn's diagnostic model is written with slightly different data formats than CISM typically uses, certain data adaptation procedures need to be performed to get the models to work together. These adaptations fall generally into two categories: staggering and transposition.
CISM places certain quantities related to velocities at the centroids of the grid on which the ice geometry is stored (usually called the "velocity grid" or the "staggered grid"). However, Pattyn's model assumes that ice geometry and velocities are co-located. This leaves us with the option of which grid to run Pattyn's model on. If it is run on the ice grid, then quantities such as initial velocity estimates need to be moved from the velocity grid to the ice grid, and the resulting velocities need to be staggered at the end of the diagnostic computation. If it is run on the velocity grid, then the ice geometry needs to be placed onto the staggered grid before the velocities are computed. In each case, moving from one grid to the other involves a linear interpolation based on the surrounding data points. More simply, to find the value of a point on the ice grid from the four surrounding points on the velocity grid or vice versa, the values of the known points are simply averaged.
Pattyn's model can be run on either grid; this option was implemented out of concern that the averaging introduced by staggering the ice geometry was leading to spuriously high velocity values is ISMIP-HOM experiments. It has since been shown that whether the velocities are computed on a staggered or unstaggered grid makes little difference in the results of these experiments. However, the options still have some utility, as the lateral boundary condition of an ice shelf is easier to apply when computing velocities on the ice grid.
The second data transformation concern is the fact that Pattyn uses a different coordinate ordering in both two-dimensional and three-dimensional arrays than does CISM. Fields in CISM are indexed using the x coordinate first, followed by the y coordinate. Pattyn reverses this, requiring all fields to be transposed. Additionally, CISM places the vertical coordinate before the horizontal coordinates, while Pattyn places them after. Thus, whereas one would index a three-dimensional CISM field by specifying (z,x,y), a three-dimensional field in Pattyn's model would be indexed by specifying (y,x,z). These alternate coordinate orderings require all field inputs and outputs to Pattyn's model to undergo a transposition.
In addition to transposing the inputs and outputs, the transposition with respect to CISM's coordinate ordering requires care when calling back out to CISM libraries. In order to avoid code duplication, Pattyn's code calls the derivative machinery in the GLIDE Derivs module. However, because the Glide module assumes that the dimension in an array is x and the second dimension is y, calls from withing Pattyn's model must call dfdx functions for a derivative with respect to y and vice versa!
These data transformation issues in general and the transposition in particular can introduce subtle errors have already caused a number of costly bugs. They must be kept under close consideration when contributing to the Pattyn diagnostic module.