New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
chap_DOM.tex in NEMO/trunk/doc/latex/NEMO/subfiles – NEMO

source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex @ 11693

Last change on this file since 11693 was 11693, checked in by nicolasmartin, 4 years ago

Macros renaming

File size: 36.3 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5\chapter{Space Domain (DOM)}
6\label{chap:DOM}
7
8% Missing things
9% -    istate: description of the initial state   ==> this has to be put elsewhere..
10%              perhaps in MISC ?  By the way the initialisation of T S and dynamics
11%              should be put outside of DOM routine (better with TRC staff and off-line
12%              tracers)
13% - geo2ocean: how to switch from geographic to mesh coordinate
14% -    domclo: closed sea and lakes....
15%              management of closea sea area: specific to global cfg, both forced and coupled
16
17\thispagestyle{plain}
18
19\chaptertoc
20
21\paragraph{Changes record} ~\\
22
23{\footnotesize
24  \begin{tabularx}{0.8\textwidth}{l||X|X}
25    Release                                                                                 &
26    Author(s)                                                                               &
27    Modifications                                                                           \\
28    \hline
29    {\em 4.0                                                                              } &
30    {\em Simon M\"{u}ller \& Andrew Coward \newline \newline
31      Simona Flavoni and Tim Graham                                                       } &
32    {\em Compatibility changes: many options moved to external domain configuration tools
33      (see \autoref{apdx:DOMCFG}). \newline
34      Updates                                                                             } \\
35    {\em 3.6                                                                              } &
36    {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec         } &
37    {\em Updates                                                                          } \\
38    {\em $\leq$ 3.4                                                                       } &
39    {\em Gurvan Madec and S\'{e}bastien Masson                                            } &
40    {\em First version                                                                    }
41  \end{tabularx}
42}
43
44\clearpage
45
46Having defined the continuous equations in \autoref{chap:MB} and
47chosen a time discretisation \autoref{chap:TD},
48we need to choose a grid for spatial discretisation and related numerical algorithms.
49In the present chapter, we provide a general description of the staggered grid used in \NEMO,
50and other relevant information about the DOM (DOMain) source code modules.
51
52%% =================================================================================================
53\section{Fundamentals of the discretisation}
54\label{sec:DOM_basics}
55
56%% =================================================================================================
57\subsection{Arrangement of variables}
58\label{subsec:DOM_cell}
59
60\begin{figure}
61  \centering
62  \includegraphics[width=0.33\textwidth]{DOM_cell}
63  \caption[Arrangement of variables in the unit cell of space domain]{
64    Arrangement of variables in the unit cell of space domain.
65    $t$ indicates scalar points where
66    temperature, salinity, density, pressure and horizontal divergence are defined.
67    $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where
68    both relative and planetary vorticities are defined.}
69  \label{fig:DOM_cell}
70\end{figure}
71
72The numerical techniques used to solve the Primitive Equations in this model are based on
73the traditional, centred second-order finite difference approximation.
74Special attention has been given to the homogeneity of the solution in the three spatial directions.
75The arrangement of variables is the same in all directions.
76It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with
77vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}).
78This is the generalisation to three dimensions of the well-known ``C'' grid in
79Arakawa's classification \citep{mesinger.arakawa_bk76}.
80The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each
81vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying
82the $\zeta$ and $f$-points.
83
84The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by
85the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
86The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on
87\autoref{tab:DOM_cell}.
88In all the following,
89subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where
90the scale factors are defined.
91Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}.
92As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
93$\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity.
94Discrete partial derivatives are formulated by
95the traditional, centred second order finite difference approximation while
96the scale factors are chosen equal to their local analytical value.
97An important point here is that the partial derivative of the scale factors must be evaluated by
98centred finite difference approximation, not from their analytical expression.
99This preserves the symmetry of the discrete set of equations and
100therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}).
101A similar, related remark can be made about the domain size:
102when needed, an area, volume, or the total ocean depth must be evaluated as
103the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section).
104
105\begin{table}
106  \centering
107  \begin{tabular}{|l|l|l|l|}
108    \hline
109    t   & $i      $ & $j      $ & $k      $ \\
110    \hline
111    u   & $i + 1/2$ & $j      $ & $k      $ \\
112    \hline
113    v   & $i      $ & $j + 1/2$ & $k      $ \\
114    \hline
115    w   & $i      $ & $j      $ & $k + 1/2$ \\
116    \hline
117    f   & $i + 1/2$ & $j + 1/2$ & $k      $ \\
118    \hline
119    uw  & $i + 1/2$ & $j      $ & $k + 1/2$ \\
120    \hline
121    vw  & $i      $ & $j + 1/2$ & $k + 1/2$ \\
122    \hline
123    fw  & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\
124    \hline
125  \end{tabular}
126  \caption[Location of grid-points]{
127    Location of grid-points as a function of integer or
128    integer and a half value of the column, line or level.
129    This indexing is only used for the writing of the semi-discrete equations.
130    In the code, the indexing uses integer values only and
131    is positive downwards in the vertical with $k=1$ at the surface.
132    (see \autoref{subsec:DOM_Num_Index})}
133  \label{tab:DOM_cell}
134\end{table}
135
136Note that the definition of the scale factors
137(\ie\ as the analytical first derivative of the transformation that
138results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
139is specific to the \NEMO\ model \citep{marti.madec.ea_JGR92}.
140As an example, a scale factor in the $i$ direction is defined locally at a $t$-point,
141whereas many other models on a C grid choose to define such a scale factor as
142the distance between the $u$-points on each side of the $t$-point.
143Relying on an analytical transformation has two advantages:
144firstly, there is no ambiguity in the scale factors appearing in the discrete equations,
145since they are first introduced in the continuous equations;
146secondly, analytical transformations encourage good practice by
147the definition of smoothly varying grids
148(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers)
149\citep{treguier.dukowicz.ea_JGR96}.
150An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}.
151\begin{figure}
152  \centering
153  \includegraphics[width=0.5\textwidth]{DOM_zgr_e3}
154  \caption[Comparison of grid-point position, vertical grid-size and scale factors]{
155    Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
156    and (b) analytically derived grid-point position and scale factors.
157    For both grids here, the same $w$-point depth has been chosen but
158    in (a) the $t$-points are set half way between $w$-points while
159    in (b) they are defined from an analytical function:
160    $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
161    Note the resulting difference between the value of the grid-size $\Delta_k$ and
162    those of the scale factor $e_k$.}
163  \label{fig:DOM_zgr_e3}
164\end{figure}
165
166%% =================================================================================================
167\subsection{Discrete operators}
168\label{subsec:DOM_operators}
169
170Given the values of a variable $q$ at adjacent points,
171the differencing and averaging operators at the midpoint between them are:
172\begin{alignat*}{2}
173  % \label{eq:DOM_di_mi}
174  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\
175  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2
176\end{alignat*}
177
178Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$.
179Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap},
180the gradient of a variable $q$ defined at a $t$-point has
181its three components defined at $u$-, $v$- and $w$-points while
182its Laplacian is defined at the $t$-point.
183These operators have the following discrete forms in the curvilinear $s$-coordinates system:
184\begin{gather*}
185  % \label{eq:DOM_grad}
186  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i
187                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j
188                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k \\
189  % \label{eq:DOM_lap}
190  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
191                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt)
192                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt]
193                  + \frac{1}{e_{3t}}
194                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt]
195\end{gather*}
196
197Following \autoref{eq:MB_curl} and \autoref{eq:MB_div},
198a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has
199its three curl components defined at $vw$-, $uw$, and $f$-points, and
200its divergence defined at $t$-points:
201\begin{multline*}
202% \label{eq:DOM_curl}
203  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}}
204                                 \Big[   \delta_{j + 1/2} (e_{3w} \, a_3)
205                                       - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\
206                               + \frac{1}{e_{2u} \, e_{3uw}}
207                                 \Big[   \delta_{k + 1/2} (e_{1u} \, a_1)
208                                       - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\
209                               + \frac{1}{e_{1f} \, e_{2f}}
210                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2)
211                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k
212\end{multline*}
213\[
214% \label{eq:DOM_div}
215  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
216                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big]
217                              + \frac{1}{e_{3t}} \delta_k (a_3)
218\]
219
220The vertical average over the whole water column is denoted by an overbar and
221is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas):
222\begin{equation}
223  \label{eq:DOM_bar}
224  \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}
225\end{equation}
226where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points,
227$k^b$ and $k^o$ are the bottom and surface $k$-indices,
228and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in
229the direction indicated by the subscript (here $k$).
230
231In continuous form, the following properties are satisfied:
232\begin{gather}
233  \label{eq:DOM_curl_grad}
234  \nabla \times \nabla q = \vect 0 \\
235  \label{eq:DOM_div_curl}
236  \nabla \cdot (\nabla \times \vect A) = 0
237\end{gather}
238
239It is straightforward to demonstrate that these properties are verified locally in discrete form as
240soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at
241vector points $(u,v,w)$.
242
243Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas.
244It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and
245$\delta_k$) are skew-symmetric linear operators,
246and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and
247$\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie
248\begin{alignat}{5}
249  \label{eq:DOM_di_adj}
250  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\
251  \label{eq:DOM_mi_adj}
252  &\sum \limits_i a_i \; \overline b^{\, i} &\equiv &  &&\sum \limits_i \overline a ^{\, i + 1/2}     &b_{i + 1/2}
253\end{alignat}
254
255In other words,
256the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and
257$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively.
258These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to
259demonstrate integral conservative properties of the discrete formulation chosen.
260
261%% =================================================================================================
262\subsection{Numerical indexing}
263\label{subsec:DOM_Num_Index}
264
265\begin{figure}
266  \centering
267  \includegraphics[width=0.33\textwidth]{DOM_index_hor}
268  \caption[Horizontal integer indexing]{
269    Horizontal integer indexing used in the \fortran\ code.
270    The dashed area indicates the cell in which
271    variables contained in arrays have the same $i$- and $j$-indices}
272  \label{fig:DOM_index_hor}
273\end{figure}
274
275The array representation used in the \fortran\ code requires an integer indexing.
276However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with
277the use of integer values for $t$-points only while
278all the other points involve integer and a half values.
279Therefore, a specific integer indexing has been defined for points other than $t$-points
280(\ie\ velocity and vorticity grid-points).
281Furthermore, the direction of the vertical indexing has been reversed and
282the surface level set at $k = 1$.
283
284%% =================================================================================================
285\subsubsection{Horizontal indexing}
286\label{subsec:DOM_Num_Index_hor}
287
288The indexing in the horizontal plane has been chosen as shown in \autoref{fig:DOM_index_hor}.
289For an increasing $i$ index ($j$ index),
290the $t$-point and the eastward $u$-point (northward $v$-point) have the same index
291(see the dashed area in \autoref{fig:DOM_index_hor}).
292A $t$-point and its nearest north-east $f$-point have the same $i$-and $j$-indices.
293
294%% =================================================================================================
295\subsubsection{Vertical indexing}
296\label{subsec:DOM_Num_Index_vertical}
297
298In the vertical, the chosen indexing requires special attention since
299the direction of the $k$-axis in the \fortran\ code is the reverse of
300that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}.
301The sea surface corresponds to the $w$-level $k = 1$,
302which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}).
303The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while
304the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}).
305Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index
306(\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction),
307in contrast to the indexing on the horizontal plane where
308the $t$-point has the same index as the nearest velocity points in
309the positive direction of the respective horizontal axis index
310(compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}).
311Since the scale factors are chosen to be strictly positive,
312a \textit{minus sign} is included in the \fortran\ implementations of
313\textit{all the vertical derivatives} of the discrete equations given in this manual in order to
314accommodate the opposing vertical index directions in implementation and documentation.
315
316\begin{figure}
317  \centering
318  \includegraphics[width=0.33\textwidth]{DOM_index_vert}
319  \caption[Vertical integer indexing]{
320    Vertical integer indexing used in the \fortran\ code.
321    Note that the $k$-axis is oriented downward.
322    The dashed area indicates the cell in which
323    variables contained in arrays have a common $k$-index.}
324  \label{fig:DOM_index_vert}
325\end{figure}
326
327%% =================================================================================================
328\section{Spatial domain configuration}
329\label{subsec:DOM_config}
330
331Two typical methods are available to specify the spatial domain configuration;
332they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in
333namelist \nam{cfg}{cfg}.
334
335If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.},
336the domain-specific parameters and fields are read from a NetCDF input file,
337whose name (without its .nc suffix) can be specified as
338the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}.
339
340If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.},
341the domain-specific parameters and fields can be provided (\eg\ analytically computed) by
342subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}.
343These subroutines can be supplied in the \path{MY_SRC} directory of the configuration,
344and default versions that configure the spatial domain for the GYRE reference configuration are
345present in the \path{./src/OCE/USR} directory.
346
347In version 4.0 there are no longer any options for reading complex bathymetries and
348performing a vertical discretisation at run-time.
349Whilst it is occasionally convenient to have a common bathymetry file and, for example,
350to run similar models with and without partial bottom boxes and/or sigma-coordinates,
351supporting such choices leads to overly complex code.
352Worse still is the difficulty of ensuring the model configurations intended to be identical are
353indeed so when the model domain itself can be altered by runtime selections.
354The code previously used to perform vertical discretisation has been incorporated into
355an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}.
356
357The next subsections summarise the parameter and fields related to
358the configuration of the whole model domain.
359These represent the minimum information that must be provided either via
360the \np{cn_domcfg}{cn\_domcfg} file or
361set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines.
362The requirements are presented in three sections:
363the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}),
364and the vertical grid (\autoref{subsec:DOM_zgr}).
365
366%% =================================================================================================
367\subsection{Domain size}
368\label{subsec:DOM_size}
369
370The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and
371\jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively.
372Note, that the variables \texttt{jpi} and \texttt{jpj} refer to
373the size of each processor subdomain when the code is run in parallel using domain decomposition
374(\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}).
375
376The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg},
377and the nominal resolution through parameter \np{nn_cfg}{nn\_cfg}
378(unless in the input file both of variables \texttt{ORCA} and \texttt{ORCA\_index} are present,
379in which case \np{cn_cfg}{cn\_cfg} and \np{nn_cfg}{nn\_cfg} are set from these values accordingly).
380
381The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}.
382See \autoref{sec:LBC_jperio} for details on the available options and
383the corresponding values for \jp{jperio}.
384
385%% =================================================================================================
386\subsection[Horizontal grid mesh (\textit{domhgr.F90}]{Horizontal grid mesh (\protect\mdl{domhgr})}
387\label{subsec:DOM_hgr}
388
389%% =================================================================================================
390\subsubsection{Required fields}
391\label{sec:DOM_hgr_fields}
392
393The explicit specification of a range of mesh-related fields are required for
394the definition of a configuration.
395These include:
396
397\begin{clines}
398int    jpiglo, jpjglo, jpkglo     /* global domain sizes                                    */
399int    jperio                     /* lateral global domain b.c.                             */
400double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */
401double gphit, gphiu, gphiv, gphif /* geographic latitude                                    */
402double e1t, e1u, e1v, e1f         /* horizontal scale factors                               */
403double e2t, e2u, e2v, e2f         /* horizontal scale factors                               */
404\end{clines}
405
406The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to
407the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$,
408evaluated at the values as specified in \autoref{tab:DOM_cell} for the respective grid-point position.
409The calculation of the values of the horizontal scale factor arrays in general additionally involves
410partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$,
411evaluated for the same arguments as $\lambda$ and $\varphi$.
412
413%% =================================================================================================
414\subsubsection{Optional fields}
415
416\begin{clines}
417                        /* Optional:                                                 */
418int    ORCA, ORCA_index /* configuration name, configuration resolution              */
419double e1e2u, e1e2v     /* U and V surfaces (if grid size reduction in some straits) */
420double ff_f, ff_t       /* Coriolis parameter (if not on the sphere)                 */
421\end{clines}
422
423\NEMO\ can support the local reduction of key strait widths by
424altering individual values of e2u or e1v at the appropriate locations.
425This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points
426(see \autoref{sec:MISC_strait} for illustrated examples).
427The key is to reduce the faces of $T$-cell
428(\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but
429not the volume of the cells.
430Doing otherwise can lead to numerical instability issues.
431In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but
432in cases where a gridsize reduction is required,
433the unaltered surface areas at $u$ and $v$ grid points
434(\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}.
435If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and
436the internal computation is suppressed.
437Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should
438set the surface-area computation flag:
439\texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation.
440
441\smallskip
442Similar logic applies to the other optional fields:
443\texttt{ff\_f} and \texttt{ff\_t} which can be used to
444provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere.
445If present these fields will be read and used and
446the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed.
447Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should
448set the Coriolis computation flag:
449\texttt{iff} to a non-zero value to suppress their re-computation.
450
451Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to
452those of $t$ points, thus no specific arrays are defined at $w$ points.
453
454%% =================================================================================================
455\subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})}
456\label{subsec:DOM_zgr}
457
458\begin{listing}
459  \nlst{namdom}
460  \caption{\forcode{&namdom}}
461  \label{lst:namdom}
462\end{listing}
463
464In the vertical, the model mesh is determined by four things:
465\begin{enumerate}
466\item the bathymetry given in meters;
467\item the number of levels of the model (\jp{jpk});
468\item the analytical transformation $z(i,j,k)$ and the vertical scale factors
469  (derivatives of the transformation); and
470\item the masking system,
471  \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid.
472\end{enumerate}
473
474\begin{figure}
475  \centering
476  \includegraphics[width=0.5\textwidth]{DOM_z_zps_s_sps}
477  \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{
478    The ocean bottom as seen by the model:
479    \begin{enumerate*}[label=(\textit{\alph*})]
480    \item $z$-coordinate with full step,
481    \item $z$-coordinate with partial step,
482    \item $s$-coordinate: terrain following representation,
483    \item hybrid $s-z$ coordinate,
484    \item hybrid $s-z$ coordinate with partial step, and
485    \item same as (e) but in the non-linear free surface
486      (\protect\np[=.false.]{ln_linssh}{ln\_linssh}).
487  \end{enumerate*}
488  Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).}
489  \label{fig:DOM_z_zps_s_sps}
490\end{figure}
491
492The choice of a vertical coordinate is made when setting up the configuration;
493it is not intended to be an option which can be changed in the middle of an experiment.
494The one exception to this statement being the choice of linear or non-linear free surface.
495In v4.0 the linear free surface option is implemented as
496a special case of the non-linear free surface.
497This is computationally wasteful since it uses the structures for time-varying 3D metrics
498for fields that (in the linear free surface case) are fixed.
499However, the linear free-surface is rarely used and
500implementing it this way means a single configuration file can support both options.
501
502By default a non-linear free surface is used
503(\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}):
504the coordinate follow the time-variation of the free surface so that
505the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f).
506When a linear free surface is assumed
507(\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}),
508the vertical coordinates are fixed in time, but
509the seawater can move up and down across the $z_0$ surface
510(in other words, the top of the ocean in not a rigid lid).
511
512Note that settings:
513\np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav}
514mentioned in the following sections appear to be namelist options but
515they are no longer truly namelist options for \NEMO.
516Their value is written to and read from the domain configuration file and
517they should be treated as fixed parameters for a particular configuration.
518They are namelist options for the \texttt{DOMAINcfg} tool that can be used to
519build the configuration file and serve both to provide a record of the choices made whilst
520building the configuration and to trigger appropriate code blocks within \NEMO.
521These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file.
522
523\medskip
524The decision on these choices must be made when the \np{cn_domcfg}{cn\_domcfg} file is constructed.
525Three main choices are offered (\autoref{fig:DOM_z_zps_s_sps}a-c):
526
527\begin{itemize}
528\item $z$-coordinate with full step bathymetry (\np[=.true.]{ln_zco}{ln\_zco}),
529\item $z$-coordinate with partial step ($zps$) bathymetry (\np[=.true.]{ln_zps}{ln\_zps}),
530\item Generalized, $s$-coordinate (\np[=.true.]{ln_sco}{ln\_sco}).
531\end{itemize}
532
533Additionally, hybrid combinations of the three main coordinates are available:
534$s-z$ or $s-zps$ coordinate (\autoref{fig:DOM_z_zps_s_sps}d and \autoref{fig:DOM_z_zps_s_sps}e).
535
536A further choice related to vertical coordinate concerns
537the presence (or not) of ocean cavities beneath ice shelves within the model domain.
538A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that
539the domain contains ocean cavities,
540otherwise the top, wet layer of the ocean will always be at the ocean surface.
541This option is currently only available for $z$- or $zps$-coordinates.
542In the latter case, partial steps are also applied at the ocean/ice shelf interface.
543
544Within the model,
545the arrays describing the grid point depths and vertical scale factors are
546three set of three dimensional arrays $(i,j,k)$ defined at
547\textit{before}, \textit{now} and \textit{after} time step.
548The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.
549They are updated at each model time step.
550The initial fixed reference coordinate system is held in variable names with a $\_0$ suffix.
551When the linear free surface option is used (\np[=.true.]{ln_linssh}{ln\_linssh}),
552\textit{before}, \textit{now} and \textit{after} arrays are initially set to
553their reference counterpart and remain fixed.
554
555%% =================================================================================================
556\subsubsection{Required fields}
557\label{sec:DOM_zgr_fields}
558
559The explicit specification of a range of fields related to the vertical grid are required for
560the definition of a configuration.
561These include:
562
563\begin{clines}
564int    ln_zco, ln_zps, ln_sco            /* flags for z-coord, z-coord with partial steps and s-coord    */
565int    ln_isfcav                         /* flag  for ice shelf cavities                                 */
566double e3t_1d, e3w_1d                    /* reference vertical scale factors at T and W points           */
567double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */
568double e3uw_0, e3vw_0                    /* vertical scale factors 3D coordinate at UW and VW points     */
569int    bottom_level, top_level           /* last wet T-points, 1st wet T-points (for ice shelf cavities) */
570                                         /* For reference:                                               */
571float  bathy_metry                       /* bathymetry used in setting top and bottom levels             */
572\end{clines}
573
574This set of vertical metrics is sufficient to describe the initial depth and thickness of
575every gridcell in the model regardless of the choice of vertical coordinate.
576With constant z-levels, e3 metrics will be uniform across each horizontal level.
577In the partial step case each e3 at the \jp{bottom\_level}
578(and, possibly, \jp{top\_level} if ice cavities are present)
579may vary from its horizontal neighbours.
580And, in s-coordinates, variations can occur throughout the water column.
581With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that
582variations occur throughout the water column with displacements related to the sea surface height.
583These variations are typically much smaller than those arising from bottom fitted coordinates.
584The values for vertical metrics supplied in the domain configuration file can be considered as
585those arising from a flat sea surface with zero elevation.
586
587The \jp{bottom\_level} and \jp{top\_level} 2D arrays define
588the \jp{bottom\_level} and top wet levels in each grid column.
589Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else).
590With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf.
591
592%% =================================================================================================
593\subsubsection{Level bathymetry and mask}
594\label{subsec:DOM_msk}
595
596From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows:
597\begin{align*}
598  tmask(i,j,k) &=
599  \begin{cases}
600    0 &\text{if $                             k <    top\_level(i,j)$} \\
601    1 &\text{if $     bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\
602    0 &\text{if $k >  bottom\_level(i,j)                            $}
603  \end{cases} \\
604  umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) \\
605  vmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j + 1,k) \\
606  fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j,    k) * tmask(i,j,k) * tmask(i + 1,j,    k) \\
607  wmask(i,j,k) &= tmask(i,j,k) * tmask(i    ,j,k - 1) \\
608  \text{with~} wmask(i,j,1) &= tmask(i,j,1)
609\end{align*}
610
611Note that, without ice shelves cavities,
612masks at $t-$ and $w-$points are identical with the numerical indexing used
613(\autoref{subsec:DOM_Num_Index}).
614Nevertheless,
615$wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface)
616exactly in the same way as for the bottom boundary.
617
618%% The specification of closed lateral boundaries requires that at least
619%% the first and last rows and columns of the \textit{mbathy} array are set to zero.
620%% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
621%% the second one and its first column equal to the last but one (and so too the mask arrays)
622%% (see \autoref{fig:LBC_jperio}).
623
624%        Closed seas
625%% =================================================================================================
626\subsection{Closed seas}
627\label{subsec:DOM_closea}
628
629When a global ocean is coupled to an atmospheric model it is better to
630represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if
631the model resolution does not allow their communication with the rest of the ocean.
632This is unnecessary when the ocean is forced by fixed atmospheric conditions,
633so these seas can be removed from the ocean domain.
634The user has the option to
635set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to
636optionally decide on the fate of any freshwater imbalance over the area.
637The options are explained in \autoref{sec:MISC_closea} but
638it should be noted here that a successful use of these options requires
639appropriate mask fields to be present in the domain configuration file.
640Among the possibilities are:
641
642\begin{clines}
643int closea_mask     /* non-zero values in closed sea areas for optional masking                */
644int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */
645int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)   */
646\end{clines}
647
648%% =================================================================================================
649\subsection{Output grid files}
650\label{subsec:DOM_meshmask}
651
652Most of the arrays relating to a particular ocean model configuration discussed in this chapter
653(grid-point position, scale factors) can be saved in a file if
654namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to
655\forcode{.true.};
656the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}.
657This is only really useful if
658the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and
659checking or confirmation is required.
660
661Alternatively, all the arrays relating to a particular ocean model configuration
662(grid-point position, scale factors, depths and masks) can be saved in
663a file called \texttt{mesh\_mask} if
664namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to
665\forcode{.true.}.
666This file contains additional fields that can be useful for post-processing applications.
667
668%% =================================================================================================
669\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
670\label{sec:DOM_DTA_tsd}
671
672\begin{listing}
673  \nlst{namtsd}
674  \caption{\forcode{&namtsd}}
675  \label{lst:namtsd}
676\end{listing}
677
678Basic initial state options are defined in \nam{tsd}{tsd}.
679By default, the ocean starts from rest (the velocity field is set to zero) and
680the initialization of temperature and salinity fields is controlled through the \np{ln_tsd_init}{ln\_tsd\_init} namelist parameter.
681
682\begin{description}
683\item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on
684  the model grid itself or on their native input data grids.
685  In the latter case,
686  the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid
687  (see \autoref{subsec:SBC_iof}).
688  The information relating to the input files are specified in
689  the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures.
690  The computation is done in the \mdl{dtatsd} module.
691\item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via
692  a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}.
693  The default version sets horizontally uniform T and profiles as used in the GYRE configuration
694  (see \autoref{sec:CFGS_gyre}).
695\end{description}
696
697\subinc{\input{../../global/epilogue}}
698
699\end{document}
Note: See TracBrowser for help on using the repository browser.