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

Last change on this file since 14257 was 14257, checked in by nicolasmartin, 3 years ago

Overall review of LaTeX sources (not tested completely as of now):

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