New URL for NEMO forge!

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

Last change on this file since 11312 was 11312, checked in by acc, 5 years ago

Updated DOM chapter for v4.0 compatibility. Much reduced since many choices and options are now external to NEMO and moved to the DOMAINcfg tool (but they are not documented there). Should consider creating an appendix as a temporary home for some of the material

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