Fundamentals of the Discretisation

Arrangement of Variables

Figure 4.1: Arrangement of variables. $t$ indicates scalar points where temperature, salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) indicates vector points, and $f$ indicates vorticity points where both relative and planetary vorticities are defined

The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, centred second-order finite difference approximation. Special attention has been given to the homogeneity of the solution in the three space directions. The arrangement of variables is the same in all directions. It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u,v,w)$ defined in the centre of each face of the cells (Fig. 4.1). This is the generalisation to three dimensions of the well-known “C” grid in Arakawa's classification [Mesinger and Arakawa, 1976]. The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.

The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on Table 4.1. In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale factors are defined. Each scale factor is defined as the local analytical value provided by (2.6). As a result, the mesh on which partial derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and $\frac{\partial}{\partial z}$ are evaluated is a uniform mesh with a grid size of unity. Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation while the scale factors are chosen equal to their local analytical value. An important point here is that the partial derivative of the scale factors must be evaluated by centred finite difference approximation, not from their analytical expression. This preserves the symmetry of the discrete set of equations and therefore satisfies many of the continuous properties (see Appendix C). A similar, related remark can be made about the domain size: when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors (see (4.8)) in the next section).

Table: Location of grid-points as a function of integer or integer and a half value of the column, line or level. This indexing is only used for the writing of the semi-discrete equation. In the code, the indexing uses integer values only and has a reverse direction in the vertical (see §4.1.3)

T	$i$	$j$	$k$
u	$i+1/2$	$j$	$k$
v	$i$	$j+1/2$	$k$
w	$i$	$j$	$k+1/2$
f	$i+1/2$	$j+1/2$	$k$
uw	$i+1/2$	$j$	$k+1/2$
vw	$i$	$j+1/2$	$k+1/2$
fw	$i+1/2$	$j+1/2$	$k+1/2$

Discrete Operators

Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at the midpoint between them are:

Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and $k+1/2$. Following (2.7a) and (2.7d), the gradient of a variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$- and $w$-points while its Laplacien is defined at $t$-point. These operators have the following discrete forms in the curvilinear $s$-coordinate system:

$$\nabla q\equiv \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \; \mathbf... ...+1/2 } [q] \; \mathbf{j} + \frac{1}{e_{3w}} \delta _{k+1/2} [q] \; \mathbf{k}$$ (4.2)

$$\begin{multline} \Delta q\equiv \frac{1}{e_{1t} e_{2t} e_{3t} } \;\left( \del... ...}} \delta_k \left[ \frac{1}{e_{3w} } \;\delta_{k+1/2} [q] \right] \end{multline}$$

Following (2.7c) and (2.7b), a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and its divergence defined at $t$-points:

$$\nabla \times {\rm {\bf A}}\equiv \frac{1}{e_{2v} e_{3vw} } \left( \delta_{j +1/2} \left[e_{3w} a_3 \right] -\delta_{k+1/2} \left[e_{2v} a_2 \right] \right) \mathbf{i}$$ (4.3)

$$+ \frac{1}{e_{2u} e_{3uw}} \left( \delta_{k+1/2} \left[e_{1u} a_1 \right] -\delta_{i +1/2} \left[e_{3w} a_3 \right] \right) \mathbf{j}$$ (4.4)

$$+ \frac{1}{e_{1f} e_{2f} } \left( \delta_{i +1/2} \left[e_{2v} a_2 \right] -\delta_{j +1/2} \left[e_{1u} a_1 \right] \right) \mathbf{k}$$ (4.5)

$$\nabla \cdot \rm {\bf A}=\frac{1}{e_{1t} e_{2t} e_{3t}}\left( \... ...e_{1v} e_{3v} a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]$$ (4.6)

In the special case of a pure $z$-coordinate system, (4.3) and (4.7) can be simplified. In this case, the vertical scale factor becomes a function of the single variable $k$ and thus does not depend on the horizontal location of a grid point. For example (4.7) reduces to:

$$\nabla \cdot \rm {\bf A}=\frac{1}{e_{1t} e_{2t}} \left( \delta_i... ...left[e_{1v} a_2 \right] \right) +\frac{1}{e_{3t}} \delta_k \left[ a_3 \right]$$

The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):

$$\bar q = \frac{1}{H} \int_{k^b}^{k^o} {q\;e_{3q} dk} \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }$$ (4.7)

where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over all grid points of the same type in the direction indicated by the subscript (here $k$).

In continuous form, the following properties are satisfied:

$$\nabla \times \nabla q ={\rm {\bf {0}}}$$ (4.8)

$$\nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0$$ (4.9)

It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as the scalar $q$ is taken at $t$-points and the vector A has its components defined at vector points $(u,v,w)$.

Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area. Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators, and further that the averaging operators $\overline{ \cdot }^{ i}$, $\overline{ \cdot }^{ k}$ and $\overline{ \cdot }^{ k}$) are symmetric linear operators, $i.e.$

$$\sum\limits_i { a_i \;\delta _i \left[ b \right]} \equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }$$ (4.10)

$$\sum\limits_i { a_i \;\overline b^{ i}} \equiv \quad \sum\limits_i {\overline a ^{ i+1/2}\;b_{i+1/2} }$$ (4.11)

In other words, the adjoint of the differencing and averaging operators are $\delta_i^*=\delta_{i+1/2}$ and ${(\overline{ \cdot }^{ i})}^*= \overline{ \cdot }^{ i+1/2}$, respectively. These two properties will be used extensively in the Appendix C to demonstrate integral conservative properties of the discrete formulation chosen.

Numerical Indexing

Figure 4.2: Horizontal integer indexing used in the FORTRAN code. The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices

The array representation used in the FORTRAN code requires an integer indexing while the analytical definition of the mesh (see §4.1.1) is associated with the use of integer values for $t$-points and both integer and integer and a half values for all the other points. Therefore a specific integer indexing must be defined for points other than $t$-points ($i.e.$ velocity and vorticity grid-points). Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k=1$.

Horizontal Indexing

The indexing in the horizontal plane has been chosen as shown in Fig.4.2. For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point (northward $v$-point) have the same index (see the dashed area in Fig.4.2). A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices.

Vertical Indexing

In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in the FORTRAN code compared to the indexing used in the semi-discrete equations and given in §4.1.1. The sea surface corresponds to the $w$-level $k=1$ which is the same index as $t$-level just below (Fig.4.3). The last $w$-level ($k=jpk$) either corresponds to the ocean floor or is inside the bathymetry while the last $t$-level is always inside the bathymetry (Fig.4.3). Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index, in opposition to what is done in the horizontal plane where it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that have the same $i$ or $j$ index (compare the dashed area in Fig.4.2 and 4.3). Since the scale factors are chosen to be strictly positive, a minus sign appears in the FORTRAN code before all the vertical derivatives of the discrete equations given in this documentation.

Figure 4.3: Vertical integer indexing used in the FORTRAN code. Note that the $k$-axis is orientated downward. The dashed area indicates the cell in which variables contained in arrays have the same $k$-index.

Domain Size

The total size of the computational domain is set by the parameters jpiglo, jpjglo and jpkdta in the $i$, $j$ and $k$ directions respectively. They are given as namelist variables in the namcfg namelist.

Note that are other namelist variables in the namcfg namelist that refer to the domain size. The two variables jpidta and jpjdta may be larger than jpiglo, jpjglo when the user wants to use only a sub-region of a given configuration. This is the "zoom" capability described in §15.3. In most applications of the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is run in parallel using domain decomposition (key_ mpp_mpi defined, see §8.3).

$$