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 @ 11597

Last change on this file since 11597 was 11597, checked in by nicolasmartin, 5 years ago

Continuation of coding rules application
Recovery of some sections deleted by the previous commit

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