IZOGRID - A new tool for setting up orthogonal curvilinear grids for oceanic modeling

Virtually all modern structured-grid ocean modeling codes are written in orthogonal curvilinear coordinates in horizontal directions, yet the overwhelming majority of modeling studies are done using very simple grid setups - mostly rectangular patches of Mercator grids rotated to proper orientation.  Furthermore, in communities like ROMS, we even observe decline in both interest and skill of creating curvilinear grids over long term.  This is caused primarily by the dissatisfaction with the existing tools and procedures for grid generation due to inability to achieve acceptable level of orthogonality errors.  Clearly, this causes underutilization of full potential of the modeling codes.

To address these issues, a new algorithm for constructing orthogonal curvilinear grids on a sphere for a fairly general geometric shape of the modeling region is implemented as a compile-once - use forever software package.  Theoretically one can use Schwartz-Christoffel conformal transform to project a curvilinear contour onto rectangle, then draw a Cartesian grid on it, and, finally, apply the inverse transform (the one which maps the rectangle back to the original contour) to the Cartesian grid in order to obtain the orthogonal curvilinear grid which fits the contour.  However, in the general case, the forward transform is an iterative algorithm of Ives and Zacharias (1989), and it is not easily invertible, nor it is feasible to apply it to a two-dimensional object (grid) as opposite to just one-dimensional (contour) because of very large number of operations.  To circumvent this, the core of the new algorithm is essentially based on the numerical solution of the inverse problem by an iterative procedure - finding such distribution of grid points along the sides of curvilinear contour, that the direct conformal mapping of it onto rectangle turns this distribution into uniform one along each side of the rectangle.  Along its way, this procedure also finds the correct aspect ratio, which makes it possible to automatically chose the numbers of grid points in each direction to yield locally the same grid spacing in both horizontal directions.  The iterative procedure itself turns out to be multilevel - i.e. an iterative loop built around another, internal iterative procedure.  Thereafter, knowing this distribution, the grid nodes inside the region are obtained solving a Dirichlet elliptic problem.  The latter is fairly standard, except that we use "mehrstellenverfahren" discretization, which yields fourth-order accuracy in the case of equal grid spacing in both directions.  The curvilinear contour is generated using splines (cubic or quintic) passing through the user-specified reference points, and, unlike all previous tools designed for the same purpose, it guarantees by the construction to yield the exact 90-degree angles at the corners of the curvilinear perimeter of grid.

Overall, with the combination of all the new features, it is shown that it is possible to achieve very small, previously unattainable level of orthogonality errors, as well as make it isotropic -- local distances between grid nodes in both directions are equal to each other.