Theoretical details

Horizontal interpolation methods

Consider an observation point ${\rm P}$ with with longitude and latitude $({\lambda_{}}_{\rm P}, \phi_{\rm P})$ and the four nearest neighbouring model grid points ${\rm A}$, ${\rm B}$, ${\rm C}$ and ${\rm D}$ with longitude and latitude ( $\lambda_{\rm A}$, $\phi_{\rm A}$), ( $\lambda_{\rm B}$, $\phi_{\rm B}$) etc. All horizontal interpolation methods implemented in NEMO estimate the value of a model variable $x$ at point $P$ as a weighted linear combination of the values of the model variables at the grid points ${\rm A}$, ${\rm B}$ etc.:

$${x_{}}_{\rm P} \hspace{-2mm} = \hspace{-2mm} \frac{1}{w} \left( {w_{}}_{\rm A} {x_{}}_{\rm A} + {w_{}}_{\rm B}... ...{\rm B} + {w_{}}_{\rm C} {x_{}}_{\rm C} + {w_{}}_{\rm D} {x_{}}_{\rm D} \right)$$ (12.1)

where ${w_{}}_{\rm A}$, ${w_{}}_{\rm B}$ etc. are the respective weights for the model field at points ${\rm A}$, ${\rm B}$ etc., and $w = {w_{}}_{\rm A} + {w_{}}_{\rm B} + {w_{}}_{\rm C} + {w_{}}_{\rm D}$.

Four different possibilities are available for computing the weights.

1.

Great-Circle distance-weighted interpolation. The weights are computed as a function of the great-circle distance $s(P, \cdot)$ between $P$ and the model grid points $A$, $B$ etc. For example, the weight given to the field ${x_{}}_{\rm A}$ is specified as the product of the distances from ${\rm P}$ to the other points:

$${w_{}}_{\rm A} = s({\rm P}, {\rm B}) s({\rm P}, {\rm C}) s({\rm P}, {\rm D})$$

where

$$s\left( {\rm P}, {\rm M} \right) \hspace{-2mm} = \hspace{-2mm} \cos^{-1} \! \left\{ \sin {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm M... ...s {\phi_{}}_{\rm M} \cos ({\lambda_{}}_{\rm M} - {\lambda_{}}_{\rm P}) \right\}$$ (12.2)

and $M$ corresponds to $B$, $C$ or $D$. A more stable form of the great-circle distance formula for small distances ($x$ near 1) involves the arcsine function ($e.g.$ see p. 101 of Daley and Barker [2001]:

$$s\left( {\rm P}, {\rm M} \right) \hspace{-2mm} = \hspace{-2mm} \sin^{-1} \! \left\{ \sqrt{ 1 - x^2 } \right\}$$

where

$$x \hspace{-2mm} = \hspace{-2mm} {a_{}}_{\rm M} {a_{}}_{\rm P} + {b_{}}_{\rm M} {b_{}}_{\rm P} + {c_{}}_{\rm M} {c_{}}_{\rm P}$$

and

$${a_{}}_{\rm M} \hspace{-2mm} = \hspace{-2mm} \sin {\phi_{}}_{\rm M},$$

$${a_{}}_{\rm P} \hspace{-2mm} = \hspace{-2mm} \sin {\phi_{}}_{\rm P},$$

$${b_{}}_{\rm M} \hspace{-2mm} = \hspace{-2mm} \cos {\phi_{}}_{\rm M} \cos {\phi_{}}_{\rm M},$$

$${b_{}}_{\rm P} \hspace{-2mm} = \hspace{-2mm} \cos {\phi_{}}_{\rm P} \cos {\phi_{}}_{\rm P},$$

$${c_{}}_{\rm M} \hspace{-2mm} = \hspace{-2mm} \cos {\phi_{}}_{\rm M} \sin {\phi_{}}_{\rm M},$$

$${c_{}}_{\rm P} \hspace{-2mm} = \hspace{-2mm} \cos {\phi_{}}_{\rm P} \sin {\phi_{}}_{\rm P}. \nonumber$$

2.

Great-Circle distance-weighted interpolation with small angle approximation. Similar to the previous interpolation but with the distance $s$ computed as

$$s\left( {\rm P}, {\rm M} \right) \hspace{-2mm} = \hspace{-2mm} \sqrt{ \left( {\phi_{}}_{\rm M} - {\phi_{}}_{\rm P} \right)^{2} +... ...bda_{}}_{\rm M} - {\lambda_{}}_{\rm P} \right)^{2} \cos^{2} {\phi_{}}_{\rm M} }$$ (12.3)

where $M$ corresponds to $A$, $B$, $C$ or $D$.

3.

Bilinear interpolation for a regular spaced grid. The interpolation is split into two 1D interpolations in the longitude and latitude directions, respectively.

4.

Bilinear remapping interpolation for a general grid. An iterative scheme that involves first mapping a quadrilateral cell into a cell with coordinates (0,0), (1,0), (0,1) and (1,1). This method is based on the SCRIP interpolation package [Jones, 1998].

Grid search

For many grids used by the NEMO model, such as the ORCA family, the horizontal grid coordinates $i$ and $j$ are not simple functions of latitude and longitude. Therefore, it is not always straightforward to determine the grid points surrounding any given observational position. Before the interpolation can be performed, a search algorithm is then required to determine the corner points of the quadrilateral cell in which the observation is located. This is the most difficult and time consuming part of the 2D interpolation procedure. A robust test for determining if an observation falls within a given quadrilateral cell is as follows. Let ${\rm P}({\lambda_{}}_{\rm P} ,{\phi_{}}_{\rm P} )$ denote the observation point, and let ${\rm A}({\lambda_{}}_{\rm A} ,{\phi_{}}_{\rm A} )$, ${\rm B}({\lambda_{}}_{\rm B} ,{\phi_{}}_{\rm B} )$, ${\rm C}({\lambda_{}}_{\rm C} ,{\phi_{}}_{\rm C} )$ and ${\rm D}({\lambda_{}}_{\rm D} ,{\phi_{}}_{\rm D} )$ denote the bottom left, bottom right, top left and top right corner points of the cell, respectively. To determine if P is inside the cell, we verify that the cross-products

$$\begin{array}{lllll} {{\bf r}_{}}_{\rm PA} \times {{\bf r}_{}... ...\; - \; {\phi_{}}_{\rm P} )] \; \widehat{\bf k} \\ \end{array}$$ (12.4)

point in the opposite direction to the unit normal $\widehat{\bf k}$ (i.e., that the coefficients of $\widehat{\bf k}$ are negative), where ${{\bf r}_{}}_{\rm PA}$, ${{\bf r}_{}}_{\rm PB}$, etc. correspond to the vectors between points P and A, P and B, etc.. The method used is similar to the method used in the SCRIP interpolation package [Jones, 1998].

In order to speed up the grid search, there is the possibility to construct a lookup table for a user specified resolution. This lookup table contains the lower and upper bounds on the $i$ and $j$ indices to be searched for on a regular grid. For each observation position, the closest point on the regular grid of this position is computed and the $i$ and $j$ ranges of this point searched to determine the precise four points surrounding the observation.

Parallel aspects of horizontal interpolation

For horizontal interpolation, there is the basic problem that the observations are unevenly distributed on the globe. In numerical models, it is common to divide the model grid into subgrids (or domains) where each subgrid is executed on a single processing element with explicit message passing for exchange of information along the domain boundaries when running on a massively parallel processor (MPP) system. This approach is used by NEMO.

For observations there is no natural distribution since the observations are not equally distributed on the globe. Two options have been made available: 1) geographical distribution; and 2) round-robin.

Geographical distribution of observations among processors

Figure 12.1: Example of the distribution of observations with the geographical distribution of observational data.

This is the simplest option in which the observations are distributed according to the domain of the grid-point parallelization. Figure 12.1 shows an example of the distribution of the in situ data on processors with a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2. The grid-point domain decomposition is clearly visible on the plot.

The advantage of this approach is that all information needed for horizontal interpolation is available without any MPP communication. Of course, this is under the assumption that we are only using a $2 \times 2$ grid-point stencil for the interpolation (e.g., bilinear interpolation). For higher order interpolation schemes this is no longer valid. A disadvantage with the above scheme is that the number of observations on each processor can be very different. If the cost of the actual interpolation is expensive relative to the communication of data needed for interpolation, this could lead to load imbalance.

Round-robin distribution of observations among processors

Figure 12.2: Example of the distribution of observations with the round-robin distribution of observational data.

An alternative approach is to distribute the observations equally among processors and use message passing in order to retrieve the stencil for interpolation. The simplest distribution of the observations is to distribute them using a round-robin scheme. Figure 12.2 shows the distribution of the in situ data on processors for the round-robin distribution of observations with a different colour for each observation on a given processor for a 4 $\times$ 2 decomposition with ORCA2 for the same input data as in Fig. 12.1. The observations are now clearly randomly distributed on the globe. In order to be able to perform horizontal interpolation in this case, a subroutine has been developed that retrieves any grid points in the global space.

Vertical interpolation operator

Vertical interpolation is achieved using either a cubic spline or linear interpolation. For the cubic spline, the top and bottom boundary conditions for the second derivative of the interpolating polynomial in the spline are set to zero. At the bottom boundary, this is done using the land-ocean mask.