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

Last change on this file since 11571 was 11571, checked in by nicolasmartin, 13 months ago

Cosmetic modifications on section titles

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