Spatial operators

To calculate the various terms and perform the time-stepping, discrete difference and interpolation operators must be designed from which all the terms, such as momentum advection and Laplacian diffusion, may be constructed. Much of the material in this section is derived from Marshall et al. (1997).

Differences

Difference operators act as the discrete form of the derivative operator. Care must be taken when calculating differences on a staggered grid. For example, the the difference of a cell-centered variable such as temperature lies on the faces in the direction of the difference, and vice versa. In principle, there are three difference operators, one for each direction

where the and subscripts indicate that the value is evaluated the eastern or western wall of the cell, and indicate the northern and southern walls, and and indicate the top and bottom walls.

Additionally, two operators must be defined for each direction to account for the staggered nature of the grid. One for taking the difference of a cell-centered variable and projecting it onto the cell faces

and another for taking the difference of a face-centered variable and projecting it onto the cell centers

The three superscript characters indicate the cell location of the output in the , , and dimensions. "f" stands for "face", "c" for "center", and "a" for "any". For example, "faa" in indicates that the output lies on the cell face in the dimension but remains at their original positions ("any") in the and dimensions.

Interpolation

In order to add or multiply variables that are defined at different points they are interpolated. In our case, linear interpolation or averaging is employed. Once again, there are two averaging operators, one for each direction,

Additionally, three averaging operators must be defined for each direction. One for taking the average of a cell-centered variable and projecting it onto the cell faces

and another for taking the average of a face-centered variable and projecting it onto the cell centers

Divergence and flux divergence

The divergence of the flux of a cell-centered quantity over the cell can be calculated as

where   is the flux with components defined normal to the faces, and is the volume of the cell. The presence of a solid boundary is indicated by setting the appropriate flux normal to the boundary to zero.

A similar divergence operator can be defined for a face-centered quantity. The divergence of, e.g., the flux of over a cell,  , is then

where is interpolated onto the cell faces where it can be multiplied by the velocities, which are then differenced and projected onto the cell centers where they added together.

Momentum advection

The advection terms that appear in model equations can be rewritten using the incompressibility (  ) as, e.g,

which can then be discretized similarly to the flux divergence operator, however, they must be discretized differently for each direction.

For example, the -momentum advection operator is discretized as

where is the average of the volumes of the cells on either side of the face in question. Calculating can be performed by interpolating and onto the cell centers then multiplying them and differencing them back onto the faces. However, in the case of the the two other terms, and , the two variables must be interpolated onto the cell edges to be multiplied then differenced back onto the cell faces.

Discretization of isotropic diffusion operators

An isotropic viscosity operator acting on vertical momentum is discretized via

where is the kinematic viscosity.

An isotropic diffusion operator acting on a tracer , on the other hand, is discretized via

Vertical integrals

Vertical integrals are converted into sums along each column. For example, the hydrostatic pressure anomaly is

where is the buoyancy perturbation. Converting it into a sum that we compute from the top downwards we get

where we converted the sum into a recursive definition for in terms of   so that the integral may be computed with operations by a single thread.

The vertical velocity may be computed from and via the continuity equation

to satisfy the incompressibility condition    to numerical precision. This also involves computing a vertical integral, in this case evaluated from the bottom up