chap_DOM.tex

\documentclass[../main/NEMO_manual]{subfiles}

\begin{document}

\chapter{Space Domain (DOM)}
\label{chap:DOM}

% Missing things
% -    istate: description of the initial state   ==> this has to be put elsewhere..
%              perhaps in MISC ?  By the way the initialisation of T S and dynamics
%              should be put outside of DOM routine (better with TRC staff and off-line
%              tracers)
% - geo2ocean: how to switch from geographic to mesh coordinate
% -    domclo: closed sea and lakes....
%              management of closea sea area: specific to global cfg, both forced and coupled

\chaptertoc

\paragraph{Changes record} ~\\

{\footnotesize
  \begin{tabularx}{0.8\textwidth}{l||X|X}
    Release                                                                                 &
      Author(s)                                                                             &
        Modifications                                                                       \\
    \hline
    {\em 5.0                                                                              } &
      {\em Sibylle Techene, Daley Calvert, and Simon M\"{u}ller                           } &
        {\em Review, revision, and removal of a summary of model variables required to define %
               the domain configuration                                                       } \\
    {\em 4.0                                                                              } &
      {\em Simon M\"{u}ller and Andrew Coward \newline \newline \newline \newline
           Simona Flavoni and Tim Graham                                                  } &
        {\em Compatibility changes: many options moved to external domain configuration %
               tools (see \autoref{apdx:DOMCFG}) \newline
             Updates                                                                      } \\
    {\em 3.6                                                                              } &
      {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec       } &
        {\em Updates                                                                      } \\
    {\em $\leq$ 3.4                                                                       } &
      {\em Gurvan Madec and S\'{e}bastien Masson                                          } &
        {\em First version                                                                }
  \end{tabularx}
}

\clearpage

Having defined the continuous equations in \autoref{chap:MB} and
chosen a time discretisation \autoref{chap:TD},
we need to choose a grid for spatial discretisation and related numerical algorithms.
In the present chapter, we provide a general description of the staggered grid used in \NEMO\,
and other relevant information about the DOM (DOMain) source code modules.
Note that \mdl{istate} and \mdl{dtatsd} are located in the \path{./src/OCE/DOM} 
directory. 
However, they are described alongside the model's time domain chapter, together 
with the restart strategy.

%% =================================================================================================
\section{Fundamentals of the discretisation}
\label{sec:DOM_basics}

%% =================================================================================================
\subsection{Arrangement of variables}
\label{subsec:DOM_cell}

\begin{figure}
  \centering
  \includegraphics[width=0.33\textwidth]{DOM_cell}
  \caption[Arrangement of variables in the unit cell of space domain]{
    Arrangement of variables in the unit cell of space domain.
    $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.}
  \label{fig:DOM_cell}
\end{figure}

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 spatial directions.
The arrangement of variables is the same in all directions.
It consists of cells centred on scalar points ($T$, \eg\ teperature $T$, salinity $S$, pressure $p$, density $\rho$) with
vector points $(u, v, w)$ (velocity), whose components are defined in the centre of each cell face (\autoref{fig:DOM_cell}).
This is the generalisation to three dimensions of the well-known ``C'' grid in
Arakawa's classification \citep{mesinger.arakawa_bk76}.
The relative and planetary vorticity, $\zeta$ and $f$, are defined at 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 (\ie\ 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 values of $(i,j,k)$ as indicated on
\autoref{tab:DOM_cell}.
In all the following,
subscripts $T$, $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where
the scale factors $e_k$ are defined.
Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}.
As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
$\pd[]{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 \autoref{apdx:INVARIANTS}).
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 product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section).

\begin{table}
  \centering
  \begin{tabular}{|l|l|l|l|}
    \hline
    t   & $i      $ & $j      $ & $k      $ \\
    \hline
    u   & $i + 1/2$ & $j      $ & $k      $ \\
    \hline
    v   & $i      $ & $j + 1/2$ & $k      $ \\
    \hline
    w   & $i      $ & $j      $ & $k + 1/2$ \\
    \hline
    f   & $i + 1/2$ & $j + 1/2$ & $k      $ \\
    \hline
    uw  & $i + 1/2$ & $j      $ & $k + 1/2$ \\
    \hline
    vw  & $i      $ & $j + 1/2$ & $k + 1/2$ \\
    \hline
    fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\
    \hline
  \end{tabular}
  \caption[Location of grid-points]{
    Location of grid-points as a function of integer or
    integer and a half values of the column, row, or level.
    This indexing is only used for the writing of the semi-discrete equations.
    In the code, the indexing uses integer values only (the relative grid-cell location has to be inferred from the context) and
    is positive downwards in the vertical with $k=1$ at the surface.
    (see \autoref{subsec:DOM_Num_Index})}
  \label{tab:DOM_cell}
\end{table}

Note that the definition of the scale factors
(\ie\ as the analytical first derivative of the transformation that
results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
is specific to the \NEMO\ model \citep{marti.madec.ea_JGR92}.
As an example, a scale factor in the $i$ direction is defined locally at a $T$-point,
whereas many other models on a C grid choose to define such a scale factor as
the distance between the $u$-points on each side of the $T$-point.
Relying on an analytical transformation has two advantages:
firstly, there is no ambiguity in the scale factors appearing in the discrete equations,
since they are first introduced in the continuous equations;
secondly, analytical transformations encourage good practice by
the definition of smoothly varying grids
(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers)
\citep{treguier.dukowicz.ea_JGR96}.
An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}.
\begin{figure}
  \centering
  \usetikzlibrary{arrows.meta}
  \usetikzlibrary{calc}
  \newcommand{\fz}[1]{5*(#1)^3-45*(#1)^2+140*(#1)-150}
  \newcommand{\dfz}[1]{15*(#1)^2-90*(#1)+140}
  \newcommand{\zind}[1]{\pgfmathparse{array({"k-2","k-1","k","k+1"},#1)}\pgfmathresult}
  \newcommand{\wind}[1]{\pgfmathparse{array({"k-3/2","k-1/2","k+1/2","k+3/2"},#1)}\pgfmathresult}
  \begin{tikzpicture}[x=1mm,y=1mm,scale=0.7,font=\tiny,thick]
      \draw[-{Latex[length=6]}] ( 30,-155) -- ( 30,10) node[right] {$z$};
      \draw[-{Latex[length=6]}] (140,-155) -- (140,10) node[right] {$z$};
      \draw ( 30,-160) node {\textbf{(a)}};
      \draw (140,-160) node {\textbf{(b)}};
      \foreach \k in {1,...,3} {
        \draw let
          \n1 = {\fz{\k-0.5}},
          \n2 = {\dfz{\k-0.5}},
          \n3 = {0.5*(\fz{\k}+\fz{\k-1})},
          \n4 = {\fz{\k}-(\fz{\k-1})} in
        ( 15,\n3)              node[left]  {$\Delta_{\zind{\k}} = \n4$ m}
        ( 32,\n3) -- ( 28,\n3) node[left]  {$T_{\zind{\k}}$}
        (142,\n1) -- (138,\n1) node[left]  {$T_{\zind{\k}}$}
        (142,\n1)              node[right] {$e_{\zind{\k}} = \n2$ m};
        \draw[Latex-Latex] let
          \n1 = {\fz{\k}},
          \n2 = {\fz{\k-1}} in
        (15,\n2) -- (15,\n1);
      };
      \foreach \k in {0,...,3} {
        \draw let
          \n1 = {\fz{\k}},
          \n2 = {\dfz{\k}} in
        ( 10,\n1) -- ( 38,\n1) node[right] {$w_{\wind{\k}}$}
        ( 85,\n1)              node        {$z_{\wind{\k}} = \n1$ m}
        (130,\n1)              node[left]  {$w_{\wind{\k}}$}
        (130,\n1) -- (150,\n1) node[right] {$e_{\wind{\k}} = \n2$ m};
      };
  \end{tikzpicture}
  \caption[Comparison of grid-point position, vertical grid-size and scale factors]{
    Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
    and (b) analytically derived grid-point position and scale factors.
    For both grids here, the same $w$-point depth has been chosen but
    in (a) the $T$-points are set halfway between $w$-points while
    in (b) they are defined from an analytical function:
    $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
    Note the resulting difference between the value of the grid-size $\Delta_k$ and
    those of the scale factor $e_k$.}
  \label{fig:DOM_zgr_e3}
\end{figure}

%% =================================================================================================
\subsection{Discrete operators}
\label{subsec:DOM_operators}

Given the values of a variable $q$ at adjacent points,
the differencing and averaging operators at the midpoint between them are:
\begin{alignat*}{2}
  % \label{eq:DOM_di_mi}
  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\
  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2
\end{alignat*}

Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$.
Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap},
the gradient of a variable $q$ defined at a $T$-point has
its three components defined at $u$-, $v$- and $w$-points while
its Laplacian is defined at the $T$-point.
These operators have the following discrete forms in the curvilinear $s$-coordinates system:
\begin{gather*}
  % \label{eq:DOM_grad}
  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i
                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j
                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\
  % \label{eq:DOM_lap}
  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt)
                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt]
                  + \frac{1}{e_{3t}}
                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt]
\end{gather*}

Following \autoref{eq:MB_curl} and \autoref{eq:MB_div},
a vector $\vect A = (a_1,a_2,a_3)$ 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:
\begin{multline*}
% \label{eq:DOM_curl}
  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}}
                                 \Big[   \delta_{j + 1/2} (e_{3w} \, a_3)
                                       - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\
                               + \frac{1}{e_{2u} \, e_{3uw}}
                                 \Big[   \delta_{k + 1/2} (e_{1u} \, a_1)
                                       - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\
                               + \frac{1}{e_{1f} \, e_{2f}}
                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2)
                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k
\end{multline*}
\[
% \label{eq:DOM_div}
  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big]
                              + \frac{1}{e_{3t}} \delta_k (a_3)
\]

The vertical average over the whole water column is denoted by an overbar and
is for a land-masked field $q$ (\ie\ a quantity that is equal to zero at land points):
\begin{equation}
  \label{eq:DOM_bar}
  \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}
\end{equation}
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 $\sum \limits_k$ 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:
\begin{gather}
  \label{eq:DOM_curl_grad}
  \nabla \times \nabla q = \vect 0 \\
  \label{eq:DOM_div_curl}
  \nabla \cdot (\nabla \times \vect A) = 0
\end{gather}

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 $\vect A$ has its components defined at
vector points $(u,v,w)$.

Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas.
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{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and
$\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie
\begin{alignat}{5}
  \label{eq:DOM_di_adj}
  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\
  \label{eq:DOM_mi_adj}
  &\sum \limits_i a_i \; \overline b^{\, i} &\equiv &  &&\sum \limits_i \overline a ^{\, i + 1/2}     &b_{i + 1/2}
\end{alignat}

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

%% =================================================================================================
\subsection{Numerical indexing}
\label{subsec:DOM_Num_Index}

\begin{figure}
  \centering
  \includegraphics[width=0.33\textwidth]{DOM_index_hor}
  \caption[Horizontal integer indexing]{
    Horizontal integer indexing used in the \fortran\ code.
    The shaded area indicates the cell in which
    variables contained in arrays have the same $i$- and $j$-indices}
  \label{fig:DOM_index_hor}
\end{figure}

The array representation used in the \fortran\ code requires an integer indexing.
However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with
the use of integer values for $T$-points only while
all the other points involve integer and a half values.
Therefore, a specific integer indexing has been defined for points other than $T$-points
(\ie\ velocity and vorticity grid-points).
Furthermore, the direction of the vertical indexing has been reversed and
the surface level set at $k = 1$.

%% =================================================================================================
\subsubsection{Horizontal indexing}
\label{subsec:DOM_Num_Index_hor}

The indexing in the horizontal plane has been chosen such that the $i$ and $j$ indices increase towards
the east and north of the domain.
$u$-, $v$- and $f$-points are distributed such that a $T$-point has the same $i$ index ($j$ index) as its nearest
eastward $u$-point (northward $v$-point) and the same $i$-and $j$-indices as its nearest north-east $f$-point
(see the shaded area in \autoref{fig:DOM_index_hor}).

%% =================================================================================================
\subsubsection{Vertical indexing}
\label{subsec:DOM_Num_Index_vertical}

In the vertical, the chosen indexing requires special attention since
the direction of the $k$-axis in the \fortran\ code is the reverse of
that used in the semi-discrete equations given in \autoref{subsec:DOM_cell}.

The $w$-level at $k = 1$ corresponds to the sea surface, while the $w$-level at $k = jpk$ either corresponds to or is
below the ocean floor. $T$-levels are distributed such that the $t$-level at index $k$ is between the $w$-levels at
indices $k$ and $k + 1$. This is in contrast to the indexing on the horizontal plane, where for example the $T$-point
at index $i$ is between the $u$-points at indices $i$ and $i - 1$ (compare the shaded areas in \autoref{fig:DOM_index_hor}
and \autoref{fig:DOM_index_vert}).

Since the scale factors are chosen to be strictly positive,
a \textit{minus sign} is included in the \fortran\ implementations of
\textit{all the vertical derivatives} of the discrete equations given in this manual in order to
accommodate the opposing vertical index directions in the implementation and documentation.

\begin{figure}
  \centering
  \includegraphics[width=0.33\textwidth]{DOM_index_vert}
  \caption[Vertical integer indexing]{
    Vertical integer indexing used in the \fortran\ code.
    Note that the $k$-axis is oriented downward.
    The topmost shaded area indicates the cell in which
    variables contained in arrays have a common $k$-index.}
  \label{fig:DOM_index_vert}
\end{figure}

%% =================================================================================================
\section{Spatial domain configuration}
\label{subsec:DOM_config}
\begin{listing}
  \nlst{namcfg}
  \caption{\forcode{&namcfg}}
  \label{lst:namcfg}
\end{listing}
Two methods are available to specify the spatial domain configuration and are selected using the
namelist parameter \np{ln_read_cfg}{ln\_read\_cfg} in namelist \nam{cfg}{cfg}:

\begin{description}
\item [{\np[=.true.]{ln_read_cfg}{ln\_read\_cfg}}] \hfill \\
    The domain-specific parameters and fields are read from a NetCDF input file,
    whose name can be specified via the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}
\item [{\np[=.false.]{ln_read_cfg}{ln\_read\_cfg}}] \hfill \\
    The domain-specific parameters and fields are set as part of a user-defined configuration
    (modules \mdl{usrdef\_nam}, \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}-
    see the \href{https://sites.nemo-ocean.io/user-guide/}{user guide}).
\end{description}

From version 4.0 there are no longer any options for reading complex bathymetries and
performing a vertical discretisation at run-time.
Whilst it is occasionally convenient to have a common bathymetry file and, for example,
to run similar models with and without partial bottom boxes and/or sigma-coordinates,
supporting such choices leads to overly complex code.
Worse still is the difficulty of ensuring the model configurations intended to be identical are
indeed so when the model domain itself can be altered by runtime selections.
The code previously used to perform vertical discretisation has therefore been incorporated into
an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}.

\medskip
The following subsections cover several topics relating to the spatial domain configuration.

%% =================================================================================================
\subsection[Horizontal grid mesh (\textit{domhgr.F90})]{Horizontal grid mesh (\protect\mdl{domhgr})}
\label{subsec:DOM_hgr}

The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to
the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$,
evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position.
The calculation of the values of the horizontal scale factor arrays in general additionally involves
partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$,
evaluated for the same arguments as $\lambda$ and $\varphi$.
Longitudes, latitudes, and horizontal scale factors at $w$-points are exactly equal to
those at $T$-points, thus no specific arrays are defined at $w$-points.

\medskip
\NEMO\ can support the local reduction of key strait widths by
altering the values of horizontal scale factors at the appropriate locations.
This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points
(see \autoref{sec:MISC_strait} for illustrated examples).
The key is to reduce the surface area of $T$-cells without changing their volume, by altering the values of the
horizontal scale factors at $u$- and $v$-points.
Doing otherwise can lead to numerical stability issues.

Normally, the surface areas are computed from the product of the horizontal scale factors
($e_{1u} * e_{2u}$ and $e_{1v} * e_{2v}$). In cases where these need to be reduced,
the modified surface areas (variables \forcode{e1e2u} and \forcode{e1e2v} at $u$- and $v$-points respectively)
must either be read from the domain configuration file (see \autoref{subsec:DOM_config}) or calculated as
part of a user-defined configuration (module \mdl{usrdef\_hgr}- see the
\href{https://sites.nemo-ocean.io/user-guide/}{user guide} for more information).

Similar logic applies to the Coriolis parameter. This is normally calculated as $2 * \Omega * \sin(\varphi)$,
but when the horizontal grid mesh is not on a sphere it must instead be specified (through variables \forcode{ff_f} and \forcode{ff_t} at
$f$- and $T$-points respectively) using one of the above methods.

%% =================================================================================================
\subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})}
\label{subsec:DOM_zgr}

The \NEMO\ vertical mesh is defined using vertical scale factors, depths, and water heights,
with each of these variables structured differently depending on the chosen coordinate system.
It is important to distinguish between the temporal structure (how the coordinate system changes over time) 
and the spatial structure (how it varies across different locations) of this system.

A code substitution mechanism defined in \textit{domzgr\_substitute.h90}
is used to replace proxy arrays describing the grid point depths (\forcode{gdept}, \forcode{gdepw}), water column
heights (\forcode{ht}, \forcode{hu}, \forcode{hv}, \forcode{hf})
and vertical scale factors (\forcode{e3t}, \forcode{e3u}, \forcode{e3v}, \forcode{e3f}, \forcode{e3w}, \forcode{e3uw},
\forcode{e3vw}) with appropriate expressions.
These expressions always include a reference array that defines the spatial structure of the vertical grid. 
In case of a non-linear free-surface, \NEMO\ uses a quasi-eulerian coordinate and 
these proxy arrays take into account a dependency with the sea surface height; \key{qco} must be specified.
In case of a linear free-surface approximation, vertical proxy arrays do not include any sea surface height dependency;
\key{linssh} must be specified.

This formulation for the vertical coordinate arrays uses less memory than in previous versions of \NEMO\ .
It takes advantage of the use of specific keys (described in the following section) that
selectively enable certain processes and data structures, allowing for more efficient memory use without
compromising model accuracy where it is needed.

\subsubsection{Quasi-eulerian COordinate (QCO)}
\label{subsec:DOM_zgr_qco}

The quasi-eulerian coordinate (specified with \key{qco}) is the new name for the star framework ($z*$ or $s*$). 
The quasi-eulerian coordinate absorbs the divergence of horizontal barotropic velocities. 
Consequently, the vertical velocity resulting from the barotropic mode is converted into a variation in thickness. 
The vertical coordinate adapts to the time-varying free surface,
making the transformation time-dependent, represented as $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f).

The vertical mesh variables (grid point depths, water heights and vertical scale factors) are substituted by an expression
(see \autoref{fig:DOM_zgr_substitute}), which in the case of the vertical scale factor at $T$-points
(\forcode{e3t}) is:

\begin{align*}
  \mathrm{e3t}(i,j,k,t)   \leftarrow  \mathrm{e3t\_0}(i,j,k)
  \left(1+ r3t(i,j,t) \mathrm{tmask}(i,j,k) \right)
\end{align*}

Where $r3t(i,j,t) = \frac{\eta(i,j,t)}{\mathrm{ht\_0}(i,j)}$, 
the ratio of sea surface height to reference water column height,
is updated at every time step by \rou{domqco\_r3c} and where
\forcode{e3t_0}, \forcode{ht_0} and \forcode{tmask} are the $T$-point variants of the reference vertical scale
factor, the reference water column height, and the land-sea mask, respectively.
Similar expressions are applied to the scale factors at $u$/$v$/$f$-points using appropriate interpolation of $\eta$.
In the expressions for the scale factors at $w$-levels, grid point depths and water heights,
the ratio is instead unmasked.

With the non-linear free-surface, all the coordinates behave more like the $s$-coordinate in that
variations occur throughout the water column with displacements related to the sea surface height.
These variations are typically much smaller than those arising from terrain-following coordinates.
As such, the reference values of the grid point depths, water heights
and vertical scale factors can be considered as those arising from a flat sea
surface with zero elevation.

\subsubsection{Linear free surface}
\label{subsec:DOM_zgr_linssh}

In the case of the linear free-surface approximation (\key{linssh} is specified), 
the free-surface variation is neglected compared to ocean depths.
Then, vertical coordinates are independent from the sea surface height and they remain fixed over time.
This setup does not treat the ocean surface as a rigid lid as it enables vertical seawater movement across the top boundary.
In this case, the vertical mesh variables are simply substituted for their time-invariant reference
counterparts, \eg\:

\begin{align*}
  \mathrm{e3t}(i,j,k,t)   \leftarrow  \mathrm{e3t\_0}(i,j,k)
\end{align*}

The same substitution is applied to all other scale factors, grid point depths and water heights.

\subsubsection{Vertical coordinate system}
\label{subsec:DOM_zgr_space}

The model mesh is initially determined by four elements when setting up the configuration:

\begin{enumerate}
 \item the bathymetry specified in meters
 \item the number of levels of the model (\forcode{jpk})
 \item the analytical transformation $z(i,j,k)$ and the reference vertical scale factors $e_{3x}^0(i,j,k)$
  (derivatives of the transformation)
 \item the masking system, \ie\ the specification of the wet model levels at each $(i,j)$ location of the horizontal grid
\end{enumerate}

\begin{figure}
  \centering
  \includegraphics[width=0.75\textwidth]{DOM_z_zps_s_sps}
  \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{
    The ocean bottom as seen by the model: 
    \textit{(a)} $z$-coordinate with full step, 
    \textit{(b)} $z$-coordinate with partial step,
    \textit{(c)} $s$-coordinate: terrain following representation,
    \textit{(d)} hybrid $s-z$ coordinate,
    \textit{(e)} hybrid $s-z$ coordinate with partial step, and
    \textit{(f)} same as (e) but in the non-linear free surface (\key{qco}).
    Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).}
  \label{fig:DOM_z_zps_s_sps}
\end{figure}

\noindent The definition of the reference vertical scale factors ($e_{3x}^0$) is determined 
by the choice of vertical coordinate system.
The vertical scale factors ($e_{3x}^0$) form the basis for deriving all vertical variables. 
This choice is fixed for the duration of an experiment and cannot be modified midway.
As with other components of the spatial domain configuration (see \autoref{subsec:DOM_config}),
it must be specified either in the domain configuration file or as part of a user-defined 
configuration (module \mdl{usrdef\_zgr}). 
For further details, refer to the \href{https://sites.nemo-ocean.io/user-guide/}{user guide}.

In addition a set of optimization keys (\key{vco}*) determine whether 1-dimensional or 3-dimensional 
definitions (as illustrated in \autoref{fig:DOM_zgr_substitute})
substitute for reference vertical variables, the intention being to minimise memory use.
In the case of \key{vco\_1d} reference vertical variables are 1-dimensional vertical arrays, 
it can only be used with a uniform grid . 
In the case of \key{vco\_1d3d} reference vertical variables are 1-dimensional vertical arrays, 
except for reference scale factors at $T$-levels, it cannot be used with non uniform grid. 
In the case of \key{vco\_3d} reference vertical variables are 3-dimensional arrays, 
it can be used sub-optimally in all cases. 

Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c):

\begin{description}
  \item [\forcode{l_zco = .true.} recommended with \key{vco\_1d}] \hfill \\
    $z$-coordinate with full step bathymetry- $e_{3x}^0$ will be uniform across each horizontal level
  \item [\forcode{l_zps = .true.} recommended with \key{vco\_1d3d}] \hfill \\
    $z$-coordinate with partial step bathymetry ($zps$-coordinate)- $e_{3x}^0$ at $T$-levels at the bottom wet level
    (and, possibly, the top wet level if ice cavities are present) may vary from its horizontal neighbours
  \item [\forcode{l_sco = .true.} recommended with \key{vco\_3d}] \hfill \\
    Generalized $s$-coordinate- variations in $e_{3x}^0$ can occur throughout the water column
\end{description}

Hybrid combinations of the three main coordinates are also available such as
$s-z$ or $s-zps$ coordinates (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e).


\begin{figure}
  \centering
  \includegraphics[width=0.5\textwidth]{DOM_zgr_substitute}
  \caption[Overview of the substitution mechanism with respect to vertical keys]
  {Overview of the vertical arrays structure based on vertical space (\key{vco\_1d}, \key{vco\_1d3d} or \key{vco\_3d})
  and time keys (\key{qco} and \key{linssh})}
  \label{fig:DOM_zgr_substitute}
\end{figure}

\subsubsection{Partial cells description \forcode{l_zps = .true.}}
\label{subsec:DOM_zps}

In $zps$-coordinates reference levels are based on the same spatially uniform levels as in $z$-coordinates. 
At the bottom (and at the top) a partial cell volume varies in order to to take into account 
solid boundaries \ie\ the bathymetry (and the ice-shelf cavities) more accurately. 

In \NEMO\ v4.2 and previous versions, partial cells were vertically shrunk, causing the mass center and the $T$-point location to shift, 
as illustrated in \autoref{fig:DOM_partial_step_scheme}. 
All scale factors associated with these cells had to be adjusted accordingly.

In \NEMO\ v5.0, partial cells are modeled as porous, consisting of both solid and liquid fractions 
that are distributed homogeneously within the cell.
Cell properties now represent the average of the liquid fraction. 
Unlike in the previous approach, the height of the $T$-points remains unchanged.
This representation is described in \citep{kevlahan_GMD15}, it is
based on Brinkman penalization where a control parameter \ie\ the porosity modifies fluxes though penalized lateral surfaces. 
This parameter is encapsulated within vertical scale factors ($e3t^0$, $e3u^0$, $e3v^0$) as illustrated in \autoref{fig:DOM_partial_step_scheme}.
In case of ocean cavities, partial cells are also applied at the top interface using the same method as for the bottom interface but upside-down.

\begin{figure}
  \centering
  \includegraphics[width=0.6\textwidth]{DOM_partial_step_scheme}
  \caption[Ocean bottom regarding $zps$-coordinate systems]{Discretisation of the horizontal difference 
  and average of tracers in the $zps$-coordinate. In \NEMO\ v4.2 a linear interpolation is used to 
  to estimate $\tilde{T}$, the tracer value at the depth of the shallower tracer point of the two adjacent 
  bottom $T$-points. While in \NEMO\ v5.0 $zps$-coordinate takes advantage of partial cells, which are modeled 
  as porous.}
  \label{fig:DOM_partial_step_scheme}
\end{figure}

%% =================================================================================================
\subsubsection{Level bathymetry and mask}
\label{subsec:DOM_msk}

The \forcode{bottom_level} and \forcode{top_level} variables define
the bottom and top wet levels in each grid column.
The values of \forcode{top_level} depend on whether ice shelf cavities are used (\autoref{subsec:ISF_dom}):
without ice cavities, \forcode{top_level} is essentially a land mask (0 on land; 1 everywhere else);
with ice cavities, in locations below an overlying ice shelf \forcode{top_level} determines the topmost wet point instead.

Based on variables \forcode{top_level} and \forcode{bottom_level}, the mask variables are defined as follows:

\begin{align*}
  tmask(i,j,k) &=
  \begin{cases}
    0 &\text{if $                             k <    top\_level(i,j)$} \\
    1 &\text{if $     bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\
    0 &\text{if $k >  bottom\_level(i,j)                            $}
  \end{cases} \\
  umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\
  vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\
  fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\
  wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\
  \text{with~} wmask(i,j,1) &= tmask(i,j,1)  \\
  wumask(i,j,k) &= umask(i,j,k) * umask(i    ,j,k - 1) \\
  \text{with~} wumask(i,j,1) &= umask(i,j,1) \\
  wvmask(i,j,k) &= vmask(i,j,k) * vmask(i    ,j,k - 1) \\
  \text{with~} wvmask(i,j,1) &= vmask(i,j,1)
\end{align*}

Note that, without ice shelves cavities,
masks at $T$- and $w$-points are identical with the numerical indexing used
(see \autoref{subsec:DOM_Num_Index}).
Nevertheless, with ocean cavities,
$wmask$ are required to deal with the top boundary (ice shelf/ocean interface)
in exactly the same way as for the bottom boundary.

%% The specification of closed lateral boundaries requires that at least
%% the first and last rows and columns of the \textit{mbathy} array are set to zero.
%% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
%% the second one and its first column equal to the last but one (and so too the mask arrays)
%% (see \autoref{fig:LBC_jperio}).

%        Closed seas
%% =================================================================================================
\subsection{Closed seas}
\label{subsec:DOM_closea}

When a global ocean is coupled to an atmospheric model it is better to
represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if
the model resolution does not allow their communication with the rest of the ocean.
This is unnecessary when the ocean is forced by fixed atmospheric conditions,
so these seas can be removed from the ocean domain.

The available options to handle closed seas are explained in \autoref{sec:MISC_closea}, but
it should be noted here that their use requires the appropriate mask fields to be present
in the domain configuration file (see \autoref{subsec:DOM_config}).

Note that, the user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to
optionally decide on the fate of any freshwater imbalance over the area.

%% =================================================================================================
\subsection{Output grid files}
\label{subsec:DOM_meshmask}

\begin{listing}
  \nlst{namdom}
  \caption{\forcode{&namdom}}
  \label{lst:namdom}
\end{listing}

The model variables describing the spatial domain configuration (latitude/longitude, scale factors, etc)
can be written to a NetCDF file by using one of the following namelist parameters:

\begin{description}
\item [{\np[=.true.]{ln_write_cfg}{ln\_write\_cfg}}] - \,\, namelist \nam{cfg}{cfg} \hfill \\
    Produces a file whose name is set by the \np{cn_domcfg_out}{cn\_domcfg\_out} namelist parameter
\item [{\np[=.true.]{ln_meshmask}{ln\_meshmask}}] - \,\, namelist \nam{dom}{dom} \hfill \\
    Produces a file named \texttt{mesh\_mask}
\end{description}

Similar files are produced by both methods, but the \texttt{mesh\_mask} file will contain additional variables
describing the land-sea mask and depth coordinate that may be useful for post-processing applications.
Both files also contain a number of the same variables found in files generated by the \texttt{DOMAINcfg} tool
(see \autoref{apdx:DOMCFG}).



\subinc{\input{../../global/epilogue}}

\end{document}