source: NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/doc/latex/NEMO/subfiles/chap_DOM.tex @ 11263

Last change on this file since 11263 was 11263, checked in by smasson, 15 months ago

dev_r10984_HPC-13 : merge with trunk@11242, see #2285

File size: 55.0 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\minitoc
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\newpage
21
22Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP},
23we need to choose a discretization on a grid, and numerical algorithms.
24In the present chapter, we provide a general description of the staggered grid used in \NEMO,
25and other information relevant to the main directory routines as well as the DOM (DOMain) directory.
26
27% ================================================================
28% Fundamentals of the Discretisation
29% ================================================================
30\section{Fundamentals of the discretisation}
31\label{sec:DOM_basics}
32
33% -------------------------------------------------------------------------------------------------------------
34%        Arrangement of Variables
35% -------------------------------------------------------------------------------------------------------------
36\subsection{Arrangement of variables}
37\label{subsec:DOM_cell}
38
39%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
40\begin{figure}[!tb]
41  \begin{center}
42    \includegraphics[width=\textwidth]{Fig_cell}
43    \caption{
44      \protect\label{fig:cell}
45      Arrangement of variables.
46      $t$ indicates scalar points where temperature, salinity, density, pressure and
47      horizontal divergence are defined.
48      $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and
49      planetary vorticities are defined.
50    }
51  \end{center}
52\end{figure}
53%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
54
55The numerical techniques used to solve the Primitive Equations in this model are based on the traditional,
56centred second-order finite difference approximation.
57Special attention has been given to the homogeneity of the solution in the three space directions.
58The arrangement of variables is the same in all directions.
59It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in
60the centre of each face of the cells (\autoref{fig:cell}).
61This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification
62\citep{mesinger.arakawa_bk76}.
63The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and
64the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points.
65
66The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that
67gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
68The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}.
69In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of
70the grid-point where the scale factors are defined.
71Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}.
72As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and
73$\pd[]{z}$ are evaluated in a uniform mesh with a grid size of unity.
74Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation
75while the scale factors are chosen equal to their local analytical value.
76An important point here is that the partial derivative of the scale factors must be evaluated by
77centred finite difference approximation, not from their analytical expression.
78This preserves the symmetry of the discrete set of equations and therefore satisfies many of
79the continuous properties (see \autoref{apdx:C}).
80A similar, related remark can be made about the domain size:
81when needed, an area, volume, or the total ocean depth must be evaluated as the sum of the relevant scale factors
82(see \autoref{eq:DOM_bar} in the next section).
83
84%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
85\begin{table}[!tb]
86  \begin{center}
87    \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
88      \hline
89      T  & $i      $ & $j      $ & $k      $ \\
90      \hline
91      u  & $i + 1/2$ & $j      $ & $k      $ \\
92      \hline
93      v  & $i      $ & $j + 1/2$ & $k      $ \\
94      \hline
95      w  & $i      $ & $j      $ & $k + 1/2$ \\
96      \hline
97      f  & $i + 1/2$ & $j + 1/2$ & $k      $ \\
98      \hline
99      uw & $i + 1/2$ & $j      $ & $k + 1/2$ \\
100      \hline
101      vw & $i      $ & $j + 1/2$ & $k + 1/2$ \\
102      \hline
103      fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\
104      \hline
105    \end{tabular}
106    \caption{
107      \protect\label{tab:cell}
108      Location of grid-points as a function of integer or integer and a half value of the column, line or level.
109      This indexing is only used for the writing of the semi -discrete equation.
110      In the code, the indexing uses integer values only and has a reverse direction in the vertical
111      (see \autoref{subsec:DOM_Num_Index})
112    }
113  \end{center}
114\end{table}
115%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
116
117% -------------------------------------------------------------------------------------------------------------
118%        Vector Invariant Formulation
119% -------------------------------------------------------------------------------------------------------------
120\subsection{Discrete operators}
121\label{subsec:DOM_operators}
122
123Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at
124the midpoint between them are:
125\begin{alignat*}{2}
126  % \label{eq:di_mi}
127  \delta_i [q]      &= &       &q (i + 1/2) - q (i - 1/2) \\
128  \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2
129\end{alignat*}
130
131Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$.
132Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at
133a $t$-point has its three components defined at $u$-, $v$- and $w$-points while
134its Laplacian is defined at $t$-point.
135These operators have the following discrete forms in the curvilinear $s$-coordinates system:
136\[
137  % \label{eq:DOM_grad}
138  \nabla q \equiv   \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i
139                  + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j
140                  + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k
141\]
142\begin{multline*}
143  % \label{eq:DOM_lap}
144  \Delta q \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
145                    \; \lt[   \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt)
146                            + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\
147                  + \frac{1}{e_{3t}}
148                              \delta_k \lt[ \frac{1              }{e_{3w}} \; \delta_{k + 1/2} [q] \rt]
149\end{multline*}
150
151Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at
152vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and
153its divergence defined at $t$-points:
154\begin{multline}
155% \label{eq:DOM_curl}
156  \nabla \times \vect A \equiv   \frac{1}{e_{2v} \, e_{3vw}}
157                                 \Big[   \delta_{j + 1/2} (e_{3w} \, a_3)
158                                       - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\
159                               + \frac{1}{e_{2u} \, e_{3uw}}
160                                 \Big[   \delta_{k + 1/2} (e_{1u} \, a_1)
161                                       - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\
162                               + \frac{1}{e_{1f} \, e_{2f}}
163                                 \Big[   \delta_{i + 1/2} (e_{2v} \, a_2)
164                                       - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k
165\end{multline}
166\begin{equation}
167% \label{eq:DOM_div}
168  \nabla \cdot \vect A \equiv   \frac{1}{e_{1t} \, e_{2t} \, e_{3t}}
169                                \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big]
170                              + \frac{1}{e_{3t}} \delta_k (a_3)
171\end{equation}
172
173The vertical average over the whole water column denoted by an overbar becomes for a quantity $q$ which
174is a masked field (i.e. equal to zero inside solid area):
175\begin{equation}
176  \label{eq:DOM_bar}
177  \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}
178\end{equation}
179where $H_q$  is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points,
180$k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $k^o$ refers to a summation over
181all grid points of the same type in the direction indicated by the subscript (here $k$).
182
183In continuous form, the following properties are satisfied:
184\begin{gather}
185  \label{eq:DOM_curl_grad}
186  \nabla \times \nabla q = \vect 0 \\
187  \label{eq:DOM_div_curl}
188  \nabla \cdot (\nabla \times \vect A) = 0
189\end{gather}
190
191It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as
192the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at
193vector points $(u,v,w)$.
194
195Let $a$ and $b$ be two fields defined on the mesh, with value zero inside continental area.
196Using integration by parts it can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$)
197are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$,
198$\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie
199\begin{alignat}{4}
200  \label{eq:DOM_di_adj}
201  &\sum \limits_i a_i \; \delta_i [b]      &\equiv &- &&\sum \limits_i \delta      _{   i + 1/2} [a] &b_{i + 1/2} \\
202  \label{eq:DOM_mi_adj}
203  &\sum \limits_i a_i \; \overline b^{\, i} &\equiv &  &&\sum \limits_i \overline a ^{\, i + 1/2}     &b_{i + 1/2}
204\end{alignat}
205
206In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and
207$(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively.
208These two properties will be used extensively in the \autoref{apdx:C} to
209demonstrate integral conservative properties of the discrete formulation chosen.
210
211% -------------------------------------------------------------------------------------------------------------
212%        Numerical Indexing
213% -------------------------------------------------------------------------------------------------------------
214\subsection{Numerical indexing}
215\label{subsec:DOM_Num_Index}
216
217%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
218\begin{figure}[!tb]
219  \begin{center}
220    \includegraphics[width=\textwidth]{Fig_index_hor}
221    \caption{
222      \protect\label{fig:index_hor}
223      Horizontal integer indexing used in the \fortran code.
224      The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices
225    }
226  \end{center}
227\end{figure}
228%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
229
230The array representation used in the \fortran code requires an integer indexing while
231the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of
232integer values for $t$-points and both integer and integer and a half values for all the other points.
233Therefore a specific integer indexing must be defined for points other than $t$-points
234(\ie velocity and vorticity grid-points).
235Furthermore, the direction of the vertical indexing has been changed so that the surface level is at $k = 1$.
236
237% -----------------------------------
238%        Horizontal Indexing
239% -----------------------------------
240\subsubsection{Horizontal indexing}
241\label{subsec:DOM_Num_Index_hor}
242
243The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}.
244For an increasing $i$ index ($j$ index),
245the $t$-point and the eastward $u$-point (northward $v$-point) have the same index
246(see the dashed area in \autoref{fig:index_hor}).
247A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices.
248
249% -----------------------------------
250%        Vertical indexing
251% -----------------------------------
252\subsubsection{Vertical indexing}
253\label{subsec:DOM_Num_Index_vertical}
254
255In the vertical, the chosen indexing requires special attention since the $k$-axis is re-orientated downward in
256the \fortran code compared to the indexing used in the semi -discrete equations and
257given in \autoref{subsec:DOM_cell}.
258The sea surface corresponds to the $w$-level $k = 1$ which is the same index as $t$-level just below
259(\autoref{fig:index_vert}).
260The last $w$-level ($k = jpk$) either corresponds to the ocean floor or is inside the bathymetry while
261the last $t$-level is always inside the bathymetry (\autoref{fig:index_vert}).
262Note that for an increasing $k$ index, a $w$-point and the $t$-point just below have the same $k$ index,
263in opposition to what is done in the horizontal plane where
264it is the $t$-point and the nearest velocity points in the direction of the horizontal axis that
265have the same $i$ or $j$ index
266(compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}).
267Since the scale factors are chosen to be strictly positive,
268a \textit{minus sign} appears in the \fortran code \textit{before all the vertical derivatives} of
269the discrete equations given in this documentation.
270
271%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
272\begin{figure}[!pt]
273  \begin{center}
274    \includegraphics[width=\textwidth]{Fig_index_vert}
275    \caption{
276      \protect\label{fig:index_vert}
277      Vertical integer indexing used in the \fortran code.
278      Note that the $k$-axis is orientated downward.
279      The dashed area indicates the cell in which variables contained in arrays have the same $k$-index.
280    }
281  \end{center}
282\end{figure}
283%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
284
285% -----------------------------------
286%        Domain Size
287% -----------------------------------
288\subsubsection{Domain size}
289\label{subsec:DOM_size}
290
291The total size of the computational domain is set by the parameters \np{jpiglo},
292\np{jpjglo} and \np{jpkglo} in the $i$, $j$ and $k$ directions respectively.
293Parameters $jpi$ and $jpj$ refer to the size of each processor subdomain when
294the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined,
295see \autoref{sec:LBC_mpp}).
296
297% ================================================================
298% Domain: List of fields needed
299% ================================================================
300\section{Needed fields}
301\label{sec:DOM_fields}
302The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that
303gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
304The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}.
305The associated scale factors are defined using the analytical first derivative of the transformation
306\autoref{eq:scale_factors}.
307Necessary fields for configuration definition are:
308
309\begin{itemize}
310\item
311  Geographic position:
312  longitude with \texttt{glamt}, \texttt{glamu}, \texttt{glamv}, \texttt{glamf} and
313  latitude  with \texttt{gphit}, \texttt{gphiu}, \texttt{gphiv}, \texttt{gphif}
314  (all respectively at T, U, V and F point)
315\item
316  Coriolis parameter (if domain not on the sphere): \texttt{ff\_f} and \texttt{ff\_t}
317  (at T and F point)
318\item
319  Scale factors:
320  \texttt{e1t}, \texttt{e1u}, \texttt{e1v} and \texttt{e1f} (on i direction),
321  \texttt{e2t}, \texttt{e2u}, \texttt{e2v} and \texttt{e2f} (on j direction) and
322  \texttt{ie1e2u\_v}, \texttt{e1e2u}, \texttt{e1e2v}. \\
323  \texttt{e1e2u}, \texttt{e1e2v} are u and v surfaces (if gridsize reduction in some straits),
324  \texttt{ie1e2u\_v} is to flag set u and v surfaces are neither read nor computed.
325\end{itemize}
326 
327These fields can be read in an domain input file which name is setted in \np{cn\_domcfg} parameter specified in
328\ngn{namcfg}.
329
330\nlst{namcfg}
331
332Or they can be defined in an analytical way in \path{MY_SRC} directory of the configuration.
333For Reference Configurations of NEMO input domain files are supplied by NEMO System Team.
334For analytical definition of input fields two routines are supplied: \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}.
335They are an example of GYRE configuration parameters, and they are available in \path{src/OCE/USR} directory,
336they provide the horizontal and vertical mesh.
337% -------------------------------------------------------------------------------------------------------------
338%        Needed fields
339% -------------------------------------------------------------------------------------------------------------
340%\subsection{List of needed fields to build DOMAIN}
341%\label{subsec:DOM_fields_list}
342
343
344% ================================================================
345% Domain: Horizontal Grid (mesh)
346% ================================================================
347\section[Horizontal grid mesh (\textit{domhgr.F90})]
348{Horizontal grid mesh (\protect\mdl{domhgr})}
349\label{sec:DOM_hgr}
350
351% -------------------------------------------------------------------------------------------------------------
352%        Coordinates and scale factors
353% -------------------------------------------------------------------------------------------------------------
354\subsection{Coordinates and scale factors}
355\label{subsec:DOM_hgr_coord_e}
356
357The ocean mesh (\ie the position of all the scalar and vector points) is defined by
358the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
359The grid-points are located at integer or integer and a half values of as indicated in \autoref{tab:cell}.
360The associated scale factors are defined using the analytical first derivative of the transformation
361\autoref{eq:scale_factors}.
362These definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr},
363which provide the horizontal and vertical meshes, respectively.
364This section deals with the horizontal mesh parameters.
365
366In a horizontal plane, the location of all the model grid points is defined from
367the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$.
368The horizontal scale factors are calculated using \autoref{eq:scale_factors}.
369For example, when the longitude and latitude are function of a single value
370($i$ and $j$, respectively) (geographical configuration of the mesh),
371the horizontal mesh definition reduces to define the wanted $\lambda(i)$, $\varphi(j)$,
372and their derivatives $\lambda'(i) \ \varphi'(j)$ in the \mdl{domhgr} module.
373The model computes the grid-point positions and scale factors in the horizontal plane as follows:
374\begin{align*}
375   \lambda_t &\equiv \text{glamt} =      \lambda (i      )
376  &\varphi_t &\equiv \text{gphit} =      \varphi (j      ) \\
377   \lambda_u &\equiv \text{glamu} =      \lambda (i + 1/2)
378  &\varphi_u &\equiv \text{gphiu} =      \varphi (j      ) \\
379   \lambda_v &\equiv \text{glamv} =      \lambda (i      )
380  &\varphi_v &\equiv \text{gphiv} =      \varphi (j + 1/2) \\
381   \lambda_f &\equiv \text{glamf} =      \lambda (i + 1/2)
382  &\varphi_f &\equiv \text{gphif} =      \varphi (j + 1/2) \\
383   e_{1t}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j      ) |
384  &e_{2t}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
385   e_{1u}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j      ) |
386  &e_{2u}    &\equiv \text{e2t}   = r_a |\varphi'(j      )                         | \\
387   e_{1v}    &\equiv \text{e1t}   = r_a |\lambda'(i      ) \; \cos\varphi(j + 1/2) |
388  &e_{2v}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         | \\
389   e_{1f}    &\equiv \text{e1t}   = r_a |\lambda'(i + 1/2) \; \cos\varphi(j + 1/2) |
390  &e_{2f}    &\equiv \text{e2t}   = r_a |\varphi'(j + 1/2)                         |
391\end{align*}
392where the last letter of each computational name indicates the grid point considered and
393$r_a$ is the earth radius (defined in \mdl{phycst} along with all universal constants).
394Note that the horizontal position of and scale factors at $w$-points are exactly equal to those of $t$-points,
395thus no specific arrays are defined at $w$-points.
396
397Note that the definition of the scale factors
398(\ie as the analytical first derivative of the transformation that
399gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$)
400is specific to the \NEMO model \citep{marti.madec.ea_JGR92}.
401As an example, $e_{1t}$ is defined locally at a $t$-point,
402whereas many other models on a C grid choose to define such a scale factor as
403the distance between the $U$-points on each side of the $t$-point.
404Relying on an analytical transformation has two advantages:
405firstly, there is no ambiguity in the scale factors appearing in the discrete equations,
406since they are first introduced in the continuous equations;
407secondly, analytical transformations encourage good practice by the definition of smoothly varying grids
408(rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}.
409An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}.
410%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
411\begin{figure}[!t]
412  \begin{center}
413    \includegraphics[width=\textwidth]{Fig_zgr_e3}
414    \caption{
415      \protect\label{fig:zgr_e3}
416      Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
417      and (b) analytically derived grid-point position and scale factors.
418      For both grids here, the same $w$-point depth has been chosen but
419      in (a) the $t$-points are set half way between $w$-points while
420      in (b) they are defined from an analytical function:
421      $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$.
422      Note the resulting difference between the value of the grid-size $\Delta_k$ and
423      those of the scale factor $e_k$.
424    }
425  \end{center}
426\end{figure}
427%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
428
429% -------------------------------------------------------------------------------------------------------------
430%        Choice of horizontal grid
431% -------------------------------------------------------------------------------------------------------------
432\subsection{Choice of horizontal grid}
433\label{subsec:DOM_hgr_msh_choice}
434
435% -------------------------------------------------------------------------------------------------------------
436%        Grid files
437% -------------------------------------------------------------------------------------------------------------
438\subsection{Output grid files}
439\label{subsec:DOM_hgr_files}
440
441All the arrays relating to a particular ocean model configuration (grid-point position, scale factors, masks)
442can be saved in files if \np{nn\_msh} $\not = 0$ (namelist variable in \ngn{namdom}).
443This can be particularly useful for plots and off-line diagnostics.
444In some cases, the user may choose to make a local modification of a scale factor in the code.
445This is the case in global configurations when restricting the width of a specific strait
446(usually a one-grid-point strait that happens to be too wide due to insufficient model resolution).
447An example is Gibraltar Strait in the ORCA2 configuration.
448When such modifications are done,
449the output grid written when \np{nn\_msh} $\not = 0$ is no more equal to the input grid.
450
451% ================================================================
452% Domain: Vertical Grid (domzgr)
453% ================================================================
454\section[Vertical grid (\textit{domzgr.F90})]
455{Vertical grid (\protect\mdl{domzgr})}
456\label{sec:DOM_zgr}
457%-----------------------------------------nam_zgr & namdom-------------------------------------------
458%
459%\nlst{namzgr}
460
461\nlst{namdom} 
462%-------------------------------------------------------------------------------------------------------------
463
464Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists.
465In the vertical, the model mesh is determined by four things:
466(1) the bathymetry given in meters;
467(2) the number of levels of the model (\jp{jpk});
468(3) the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and
469(4) the masking system, \ie the number of wet model levels at each
470$(i,j)$ column of points.
471
472%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
473\begin{figure}[!tb]
474  \begin{center}
475    \includegraphics[width=\textwidth]{Fig_z_zps_s_sps}
476    \caption{
477      \protect\label{fig:z_zps_s_sps}
478      The ocean bottom as seen by the model:
479      (a) $z$-coordinate with full step,
480      (b) $z$-coordinate with partial step,
481      (c) $s$-coordinate: terrain following representation,
482      (d) hybrid $s-z$ coordinate,
483      (e) hybrid $s-z$ coordinate with partial step, and
484      (f) same as (e) but in the non-linear free surface (\protect\np{ln\_linssh}\forcode{ = .false.}).
485      Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).
486    }
487  \end{center}
488\end{figure}
489%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
490
491The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,
492must be done once of all at the beginning of an experiment.
493It is not intended as an option which can be enabled or disabled in the middle of an experiment.
494Three main choices are offered (\autoref{fig:z_zps_s_sps}):
495$z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}),
496$z$-coordinate with partial step bathymetry (\np{ln\_zps}\forcode{ = .true.}),
497or generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}).
498Hybridation of the three main coordinates are available:
499$s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps} and \autoref{fig:z_zps_s_sps}).
500By default a non-linear free surface is used: the coordinate follow the time-variation of the free surface so that
501the transformation is time dependent: $z(i,j,k,t)$ (\autoref{fig:z_zps_s_sps}).
502When a linear free surface is assumed (\np{ln\_linssh}\forcode{ = .true.}),
503the vertical coordinate are fixed in time, but the seawater can move up and down across the $z_0$ surface
504(in other words, the top of the ocean in not a rigid-lid).
505The last choice in terms of vertical coordinate concerns the presence (or not) in
506the model domain of ocean cavities beneath ice shelves.
507Setting \np{ln\_isfcav} to true allows to manage ocean cavities, otherwise they are filled in.
508This option is currently only available in $z$- or $zps$-coordinate,
509and partial step are also applied at the ocean/ice shelf interface.
510
511Contrary to the horizontal grid, the vertical grid is computed in the code and no provision is made for
512reading it from a file.
513The only input file is the bathymetry (in meters) (\ifile{bathy\_meter})
514\footnote{
515  N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the \ifile{bathy\_meter} file,
516  so that the computation of the number of wet ocean point in each water column is by-passed}.
517If \np{ln\_isfcav}\forcode{ = .true.}, an extra file input file (\ifile{isf\_draft\_meter}) describing
518the ice shelf draft (in meters) is needed.
519
520After reading the bathymetry, the algorithm for vertical grid definition differs between the different options:
521\begin{description}
522\item[\textit{zco}]
523  set a reference coordinate transformation $z_0(k)$, and set $z(i,j,k,t) = z_0(k)$.
524\item[\textit{zps}]
525  set a reference coordinate transformation $z_0(k)$, and calculate the thickness of the deepest level at
526  each $(i,j)$ point using the bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
527\item[\textit{sco}]
528  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
529  set the three-dimensional transformation.
530\item[\textit{s-z} and \textit{s-zps}]
531  smooth the bathymetry to fulfill the hydrostatic consistency criteria and
532  set the three-dimensional transformation $z(i,j,k)$,
533  and possibly introduce masking of extra land points to better fit the original bathymetry file.
534\end{description}
535%%%
536\gmcomment{   add the description of the smoothing:  envelop topography...}
537%%%
538
539Unless a linear free surface is used (\np{ln\_linssh}\forcode{ = .false.}),
540the arrays describing the grid point depths and vertical scale factors are three set of
541three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step.
542The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively.
543They are updated at each model time step using a fixed reference coordinate system which
544computer names have a $\_0$ suffix.
545When the linear free surface option is used (\np{ln\_linssh}\forcode{ = .true.}), \textit{before},
546\textit{now} and \textit{after} arrays are simply set one for all to their reference counterpart.
547
548% -------------------------------------------------------------------------------------------------------------
549%        Meter Bathymetry
550% -------------------------------------------------------------------------------------------------------------
551\subsection{Meter bathymetry}
552\label{subsec:DOM_bathy}
553
554Three options are possible for defining the bathymetry, according to the namelist variable \np{nn\_bathy}
555(found in \ngn{namdom} namelist):
556\begin{description}
557\item[\np{nn\_bathy}\forcode{ = 0}]:
558  a flat-bottom domain is defined.
559  The total depth $z_w (jpk)$ is given by the coordinate transformation.
560  The domain can either be a closed basin or a periodic channel depending on the parameter \np{jperio}.
561\item[\np{nn\_bathy}\forcode{ = -1}]:
562  a domain with a bump of topography one third of the domain width at the central latitude.
563  This is meant for the "EEL-R5" configuration, a periodic or open boundary channel with a seamount.
564\item[\np{nn\_bathy}\forcode{ = 1}]:
565  read a bathymetry and ice shelf draft (if needed).
566  The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters) at
567  each grid point of the model grid.
568  The bathymetry is usually built by interpolating a standard bathymetry product (\eg ETOPO2) onto
569  the horizontal ocean mesh.
570  Defining the bathymetry also defines the coastline: where the bathymetry is zero,
571  no model levels are defined (all levels are masked).
572
573  The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters) at
574  each grid point of the model grid.
575  This file is only needed if \np{ln\_isfcav}\forcode{ = .true.}.
576  Defining the ice shelf draft will also define the ice shelf edge and the grounding line position.
577\end{description}
578
579When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies
580(\eg great lakes, Caspian sea...) even if the model resolution does not allow their communication with
581the rest of the ocean.
582This is unnecessary when the ocean is forced by fixed atmospheric conditions,
583so these seas can be removed from the ocean domain.
584The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}),
585but the code has to be adapted to the user's configuration.
586
587% -------------------------------------------------------------------------------------------------------------
588%        z-coordinate  and reference coordinate transformation
589% -------------------------------------------------------------------------------------------------------------
590\subsection[$Z$-coordinate (\forcode{ln_zco = .true.}) and ref. coordinate]
591{$Z$-coordinate (\protect\np{ln\_zco}\forcode{ = .true.}) and reference coordinate}
592\label{subsec:DOM_zco}
593
594%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
595\begin{figure}[!tb]
596  \begin{center}
597    \includegraphics[width=\textwidth]{Fig_zgr}
598    \caption{
599      \protect\label{fig:zgr}
600      Default vertical mesh for ORCA2: 30 ocean levels (L30).
601      Vertical level functions for (a) T-point depth and (b) the associated scale factor as computed from
602      \autoref{eq:DOM_zgr_ana_1} using \autoref{eq:DOM_zgr_coef} in $z$-coordinate.
603    }
604  \end{center}
605\end{figure}
606%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
607
608The reference coordinate transformation $z_0(k)$ defines the arrays $gdept_0$ and $gdepw_0$ for $t$- and $w$-points,
609respectively.
610As indicated on \autoref{fig:index_vert} \jp{jpk} is the number of $w$-levels.
611$gdepw_0(1)$ is the ocean surface.
612There are at most \jp{jpk}-1 $t$-points inside the ocean,
613the additional $t$-point at $jk = jpk$ is below the sea floor and is not used.
614The vertical location of $w$- and $t$-levels is defined from the analytic expression of the depth $z_0(k)$ whose
615analytical derivative with respect to $k$ provides the vertical scale factors.
616The user must provide the analytical expression of both $z_0$ and its first derivative with respect to $k$.
617This is done in routine \mdl{domzgr} through statement functions,
618using parameters provided in the \ngn{namcfg} namelist.
619
620It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr}\forcode{ = 0}).
621In that case, the parameters \jp{jpk} (number of $w$-levels) and
622\np{pphmax} (total ocean depth in meters) fully define the grid.
623
624For climate-related studies it is often desirable to concentrate the vertical resolution near the ocean surface.
625The following function is proposed as a standard for a $z$-coordinate (with either full or partial steps):
626\begin{gather}
627  \label{eq:DOM_zgr_ana_1}
628    z_0  (k) = h_{sur} - h_0 \; k - \; h_1 \; \log  \big[ \cosh ((k - h_{th}) / h_{cr}) \big] \\
629    e_3^0(k) = \lt|    - h_0      -    h_1 \; \tanh \big[        (k - h_{th}) / h_{cr}  \big] \rt|
630\end{gather}
631where $k = 1$ to \jp{jpk} for $w$-levels and $k = 1$ to $k = 1$ for $T-$levels.
632Such an expression allows us to define a nearly uniform vertical location of levels at the ocean top and bottom with
633a smooth hyperbolic tangent transition in between (\autoref{fig:zgr}).
634
635If the ice shelf cavities are opened (\np{ln\_isfcav}\forcode{ = .true.}), the definition of $z_0$ is the same.
636However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to:
637\begin{equation}
638  \label{eq:DOM_zgr_ana_2}
639  \begin{split}
640    e_3^T(k) &= z_W (k + 1) - z_W (k    ) \\
641    e_3^W(k) &= z_T (k    ) - z_T (k - 1)
642  \end{split}
643\end{equation}
644This formulation decrease the self-generated circulation into the ice shelf cavity
645(which can, in extreme case, leads to blow up).\\
646 
647The most used vertical grid for ORCA2 has $10~m$ ($500~m$) resolution in the surface (bottom) layers and
648a depth which varies from 0 at the sea surface to a minimum of $-5000~m$.
649This leads to the following conditions:
650\begin{equation}
651  \label{eq:DOM_zgr_coef}
652  \begin{array}{ll}
653    e_3 (1   + 1/2) =  10. & z(1  ) =     0. \\
654    e_3 (jpk - 1/2) = 500. & z(jpk) = -5000.
655  \end{array}
656\end{equation}
657
658With the choice of the stretching $h_{cr} = 3$ and the number of levels \jp{jpk}~$= 31$,
659the four coefficients $h_{sur}$, $h_0$, $h_1$, and $h_{th}$ in
660\autoref{eq:DOM_zgr_ana_2} have been determined such that
661\autoref{eq:DOM_zgr_coef} is satisfied, through an optimisation procedure using a bisection method.
662For the first standard ORCA2 vertical grid this led to the following values:
663$h_{sur} = 4762.96$, $h_0 = 255.58, h_1 = 245.5813$, and $h_{th} = 21.43336$.
664The resulting depths and scale factors as a function of the model levels are shown in
665\autoref{fig:zgr} and given in \autoref{tab:orca_zgr}.
666Those values correspond to the parameters \np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.
667
668Rather than entering parameters $h_{sur}$, $h_0$, and $h_1$ directly, it is possible to recalculate them.
669In that case the user sets \np{ppsur}~$=$~\np{ppa0}~$=$~\np{ppa1}~$= 999999$.,
670in \ngn{namcfg} namelist, and specifies instead the four following parameters:
671\begin{itemize}
672\item
673  \np{ppacr}~$= h_{cr}$: stretching factor (nondimensional).
674  The larger \np{ppacr}, the smaller the stretching.
675  Values from $3$ to $10$ are usual.
676\item
677  \np{ppkth}~$= h_{th}$: is approximately the model level at which maximum stretching occurs
678  (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
679\item
680  \np{ppdzmin}: minimum thickness for the top layer (in meters).
681\item
682  \np{pphmax}: total depth of the ocean (meters).
683\end{itemize}
684As an example, for the $45$ layers used in the DRAKKAR configuration those parameters are:
685\jp{jpk}~$= 46$, \np{ppacr}~$= 9$, \np{ppkth}~$= 23.563$, \np{ppdzmin}~$= 6~m$, \np{pphmax}~$= 5750~m$.
686
687%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
688\begin{table}
689  \begin{center}
690    \begin{tabular}{c||r|r|r|r}
691      \hline
692      \textbf{LEVEL} & \textbf{gdept\_1d} & \textbf{gdepw\_1d} & \textbf{e3t\_1d } & \textbf{e3w\_1d} \\
693      \hline
694      1              & \textbf{     5.00} &               0.00 & \textbf{   10.00} &            10.00 \\
695      \hline
696      2              & \textbf{    15.00} &              10.00 & \textbf{   10.00} &            10.00 \\
697      \hline
698      3              & \textbf{    25.00} &              20.00 & \textbf{   10.00} &            10.00 \\
699      \hline
700      4              & \textbf{    35.01} &              30.00 & \textbf{   10.01} &            10.00 \\
701      \hline
702      5              & \textbf{    45.01} &              40.01 & \textbf{   10.01} &            10.01 \\
703      \hline
704      6              & \textbf{    55.03} &              50.02 & \textbf{   10.02} &            10.02 \\
705      \hline
706      7              & \textbf{    65.06} &              60.04 & \textbf{   10.04} &            10.03 \\
707      \hline
708      8              & \textbf{    75.13} &              70.09 & \textbf{   10.09} &            10.06 \\
709      \hline
710      9              & \textbf{    85.25} &              80.18 & \textbf{   10.17} &            10.12 \\
711      \hline
712      10             & \textbf{    95.49} &              90.35 & \textbf{   10.33} &            10.24 \\
713      \hline
714      11             & \textbf{   105.97} &             100.69 & \textbf{   10.65} &            10.47 \\
715      \hline
716      12             & \textbf{   116.90} &             111.36 & \textbf{   11.27} &            10.91 \\
717      \hline
718      13             & \textbf{   128.70} &             122.65 & \textbf{   12.47} &            11.77 \\
719      \hline
720      14             & \textbf{   142.20} &             135.16 & \textbf{   14.78} &            13.43 \\
721      \hline
722      15             & \textbf{   158.96} &             150.03 & \textbf{   19.23} &            16.65 \\
723      \hline
724      16             & \textbf{   181.96} &             169.42 & \textbf{   27.66} &            22.78 \\
725      \hline
726      17             & \textbf{   216.65} &             197.37 & \textbf{   43.26} &            34.30 \\
727      \hline
728      18             & \textbf{   272.48} &             241.13 & \textbf{   70.88} &            55.21 \\
729      \hline
730      19             & \textbf{   364.30} &             312.74 & \textbf{  116.11} &            90.99 \\
731      \hline
732      20             & \textbf{   511.53} &             429.72 & \textbf{  181.55} &           146.43 \\
733      \hline
734      21             & \textbf{   732.20} &             611.89 & \textbf{  261.03} &           220.35 \\
735      \hline
736      22             & \textbf{  1033.22} &             872.87 & \textbf{  339.39} &           301.42 \\
737      \hline
738      23             & \textbf{  1405.70} &            1211.59 & \textbf{  402.26} &           373.31 \\
739      \hline
740      24             & \textbf{  1830.89} &            1612.98 & \textbf{  444.87} &           426.00 \\
741      \hline
742      25             & \textbf{  2289.77} &            2057.13 & \textbf{  470.55} &           459.47 \\
743      \hline
744      26             & \textbf{  2768.24} &            2527.22 & \textbf{  484.95} &           478.83 \\
745      \hline
746      27             & \textbf{  3257.48} &            3011.90 & \textbf{  492.70} &           489.44 \\
747      \hline
748      28             & \textbf{  3752.44} &            3504.46 & \textbf{  496.78} &           495.07 \\
749      \hline
750      29             & \textbf{  4250.40} &            4001.16 & \textbf{  498.90} &           498.02 \\
751      \hline
752      30             & \textbf{  4749.91} &            4500.02 & \textbf{  500.00} &           499.54 \\
753      \hline
754      31             & \textbf{  5250.23} &            5000.00 & \textbf{  500.56} &           500.33 \\
755      \hline
756    \end{tabular}
757  \end{center}
758  \caption{
759    \protect\label{tab:orca_zgr}
760    Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed from
761    \autoref{eq:DOM_zgr_ana_2} using the coefficients given in \autoref{eq:DOM_zgr_coef}
762  }
763\end{table}
764%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
765
766% -------------------------------------------------------------------------------------------------------------
767%        z-coordinate with partial step
768% -------------------------------------------------------------------------------------------------------------
769\subsection[$Z$-coordinate with partial step (\forcode{ln_zps = .true.})]
770{$Z$-coordinate with partial step (\protect\np{ln\_zps}\forcode{ = .true.})}
771\label{subsec:DOM_zps}
772%--------------------------------------------namdom-------------------------------------------------------
773
774\nlst{namdom} 
775%--------------------------------------------------------------------------------------------------------------
776
777In $z$-coordinate partial step,
778the depths of the model levels are defined by the reference analytical function $z_0(k)$ as described in
779the previous section, \textit{except} in the bottom layer.
780The thickness of the bottom layer is allowed to vary as a function of geographical location $(\lambda,\varphi)$ to
781allow a better representation of the bathymetry, especially in the case of small slopes
782(where the bathymetry varies by less than one level thickness from one grid point to the next).
783The reference layer thicknesses $e_{3t}^0$ have been defined in the absence of bathymetry.
784With partial steps, layers from 1 to \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$.
785The model deepest layer (\jp{jpk}-1) is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$:
786the maximum thickness allowed is $2*e_{3t}(jpk - 1)$.
787This has to be kept in mind when specifying values in \ngn{namdom} namelist,
788as the maximum depth \np{pphmax} in partial steps:
789for example, with \np{pphmax}~$= 5750~m$ for the DRAKKAR 45 layer grid,
790the maximum ocean depth allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk - 1)$ being $250~m$).
791Two variables in the namdom namelist are used to define the partial step vertical grid.
792The mimimum water thickness (in meters) allowed for a cell partially filled with bathymetry at level jk is
793the minimum of \np{rn\_e3zps\_min} (thickness in meters, usually $20~m$) or $e_{3t}(jk)*$\np{rn\_e3zps\_rat}
794(a fraction, usually 10\%, of the default thickness $e_{3t}(jk)$).
795
796\gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  }
797
798% -------------------------------------------------------------------------------------------------------------
799%        s-coordinate
800% -------------------------------------------------------------------------------------------------------------
801\subsection[$S$-coordinate (\forcode{ln_sco = .true.})]
802{$S$-coordinate (\protect\np{ln\_sco}\forcode{ = .true.})}
803\label{subsec:DOM_sco}
804%------------------------------------------nam_zgr_sco---------------------------------------------------
805%
806%\nlst{namzgr_sco}
807%--------------------------------------------------------------------------------------------------------------
808Options are defined in \ngn{namzgr\_sco}.
809In $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}), the depth and thickness of the model levels are defined from
810the product of a depth field and either a stretching function or its derivative, respectively:
811
812\begin{align*}
813  % \label{eq:DOM_sco_ana}
814  z(k)   &= h(i,j) \; z_0 (k) \\
815  e_3(k) &= h(i,j) \; z_0'(k)
816\end{align*}
817
818where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point location in the horizontal and
819$z_0(k)$ is a function which varies from $0$ at the sea surface to $1$ at the ocean bottom.
820The depth field $h$ is not necessary the ocean depth,
821since a mixed step-like and bottom-following representation of the topography can be used
822(\autoref{fig:z_zps_s_sps}) or an envelop bathymetry can be defined (\autoref{fig:z_zps_s_sps}).
823The namelist parameter \np{rn\_rmax} determines the slope at which
824the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate.
825The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as
826the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
827
828Options for stretching the coordinate are provided as examples,
829but care must be taken to ensure that the vertical stretch used is appropriate for the application.
830
831The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true
832(\np{ln\_s\_SH94}\forcode{ = .false.} and \np{ln\_s\_SF12}\forcode{ = .false.}).
833This uses a depth independent $\tanh$ function for the stretching \citep{madec.delecluse.ea_JPO96}:
834
835\[
836  z = s_{min} + C (s) (H - s_{min})
837  % \label{eq:SH94_1}
838\]
839
840where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and
841allows a $z$-coordinate to placed on top of the stretched coordinate,
842and $z$ is the depth (negative down from the asea surface).
843\begin{gather*}
844  s = - \frac{k}{n - 1} \quad \text{and} \quad 0 \leq k \leq n - 1
845  % \label{eq:DOM_s}
846 \\
847  % \label{eq:DOM_sco_function}
848  C(s) = \frac{[\tanh(\theta \, (s + b)) - \tanh(\theta \, b)]}{2 \; \sinh(\theta)}
849\end{gather*}
850
851A stretching function,
852modified from the commonly used \citet{song.haidvogel_JCP94} stretching (\np{ln\_s\_SH94}\forcode{ = .true.}),
853is also available and is more commonly used for shelf seas modelling:
854
855\[
856  C(s) =   (1 - b) \frac{\sinh(\theta s)}{\sinh(\theta)}
857         + b       \frac{\tanh \lt[ \theta \lt(s + \frac{1}{2} \rt) \rt] -   \tanh \lt( \frac{\theta}{2} \rt)}
858                        {                                                  2 \tanh \lt( \frac{\theta}{2} \rt)}
859  % \label{eq:SH94_2}
860\]
861
862%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
863\begin{figure}[!ht]
864  \begin{center}
865    \includegraphics[width=\textwidth]{Fig_sco_function}
866    \caption{
867      \protect\label{fig:sco_function}
868      Examples of the stretching function applied to a seamount;
869      from left to right: surface, surface and bottom, and bottom intensified resolutions
870    }
871  \end{center}
872\end{figure}
873%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
874
875where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to
876the stretched coordinate, and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and
877bottom control parameters such that $0 \leqslant \theta \leqslant 20$, and $0 \leqslant b \leqslant 1$.
878$b$ has been designed to allow surface and/or bottom increase of the vertical resolution
879(\autoref{fig:sco_function}).
880
881Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows a fixed surface resolution in
882an analytical terrain-following stretching \citet{siddorn.furner_OM13}.
883In this case the a stretching function $\gamma$ is defined such that:
884
885\begin{equation}
886  z = - \gamma h \quad \text{with} \quad 0 \leq \gamma \leq 1
887  % \label{eq:z}
888\end{equation}
889
890The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
891
892\begin{gather*}
893  % \label{eq:DOM_gamma_deriv}
894  \gamma =   A \lt( \sigma   - \frac{1}{2} (\sigma^2     + f (\sigma)) \rt)
895           + B \lt( \sigma^3 - f           (\sigma) \rt) + f (\sigma)       \\
896  \intertext{Where:}
897  % \label{eq:DOM_gamma}
898  f(\sigma) = (\alpha + 2) \sigma^{\alpha + 1} - (\alpha + 1) \sigma^{\alpha + 2}
899  \quad \text{and} \quad \sigma = \frac{k}{n - 1}
900\end{gather*}
901
902This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
903the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards
904the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and
905user prescribed surface (\np{rn\_zs}) and bottom depths.
906The bottom cell depth in this example is given as a function of water depth:
907
908\[
909  % \label{eq:DOM_zb}
910  Z_b = h a + b
911\]
912
913where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively.
914
915%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
916\begin{figure}[!ht]
917  \includegraphics[width=\textwidth]{Fig_DOM_compare_coordinates_surface}
918  \caption{
919    A comparison of the \citet{song.haidvogel_JCP94} $S$-coordinate (solid lines),
920    a 50 level $Z$-coordinate (contoured surfaces) and
921    the \citet{siddorn.furner_OM13} $S$-coordinate (dashed lines) in the surface $100~m$ for
922    a idealised bathymetry that goes from $50~m$ to $5500~m$ depth.
923    For clarity every third coordinate surface is shown.
924  }
925  \label{fig:fig_compare_coordinates_surface}
926\end{figure}
927 % >>>>>>>>>>>>>>>>>>>>>>>>>>>>
928
929This gives a smooth analytical stretching in computational space that is constrained to
930given specified surface and bottom grid cell thicknesses in real space.
931This is not to be confused with the hybrid schemes that
932superimpose geopotential coordinates on terrain following coordinates thus
933creating a non-analytical vertical coordinate that
934therefore may suffer from large gradients in the vertical resolutions.
935This stretching is less straightforward to implement than the \citet{song.haidvogel_JCP94} stretching,
936but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.
937
938As with the \citet{song.haidvogel_JCP94} stretching the stretch is only applied at depths greater than
939the critical depth $h_c$.
940In this example two options are available in depths shallower than $h_c$,
941with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false
942(the z-coordinate being equal to the depths of the stretched coordinate at $h_c$).
943
944Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as
945large slopes lead to hydrostatic consistency.
946A hydrostatic consistency parameter diagnostic following \citet{haney_JPO91} has been implemented,
947and is output as part of the model mesh file at the start of the run.
948
949% -------------------------------------------------------------------------------------------------------------
950%        z*- or s*-coordinate
951% -------------------------------------------------------------------------------------------------------------
952\subsection[\zstar- or \sstar-coordinate (\forcode{ln_linssh = .false.})]
953{\zstar- or \sstar-coordinate (\protect\np{ln\_linssh}\forcode{ = .false.})}
954\label{subsec:DOM_zgr_star}
955
956This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.
957
958%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
959
960% -------------------------------------------------------------------------------------------------------------
961%        level bathymetry and mask
962% -------------------------------------------------------------------------------------------------------------
963\subsection{Level bathymetry and mask}
964\label{subsec:DOM_msk}
965
966Whatever the vertical coordinate used, the model offers the possibility of representing the bottom topography with
967steps that follow the face of the model cells (step like topography) \citep{madec.delecluse.ea_JPO96}.
968The distribution of the steps in the horizontal is defined in a 2D integer array, mbathy, which
969gives the number of ocean levels (\ie those that are not masked) at each $t$-point.
970mbathy is computed from the meter bathymetry using the definiton of gdept as the number of $t$-points which
971gdept $\leq$ bathy.
972
973Modifications of the model bathymetry are performed in the \textit{bat\_ctl} routine (see \mdl{domzgr} module) after
974mbathy is computed.
975Isolated grid points that do not communicate with another ocean point at the same level are eliminated.
976
977As for the representation of bathymetry, a 2D integer array, misfdep, is created.
978misfdep defines the level of the first wet $t$-point.
979All the cells between $k = 1$ and $misfdep(i,j) - 1$ are masked.
980By default, $misfdep(:,:) = 1$ and no cells are masked.
981
982In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into
983the cavities are performed in the \textit{zgr\_isf} routine.
984The compatibility between ice shelf draft and bathymetry is checked.
985All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded (\ie masked).
986If only one cell on the water column is opened at $t$-, $u$- or $v$-points,
987the bathymetry or the ice shelf draft is dug to fit this constrain.
988If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.
989
990From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows:
991\begin{alignat*}{2}
992  tmask(i,j,k) &= &  &
993    \begin{cases}
994                  0 &\text{if $                  k  <    misfdep(i,j)$} \\
995                  1 &\text{if $misfdep(i,j) \leq k \leq   mbathy(i,j)$} \\
996                  0 &\text{if $                  k  >     mbathy(i,j)$}
997    \end{cases}
998  \\
999  umask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\
1000  vmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j + 1,k) \\
1001  fmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i + 1,j,    k) \\
1002               &  &* &tmask(i,j,k) * tmask(i + 1,j,    k) \\
1003  wmask(i,j,k) &= &  &tmask(i,j,k) * tmask(i    ,j,k - 1) \\
1004  \text{with~} wmask(i,j,1) &= & &tmask(i,j,1)
1005\end{alignat*}
1006
1007Note that, without ice shelves cavities,
1008masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}).
1009Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface)
1010exactly in the same way as for the bottom boundary.
1011
1012The specification of closed lateral boundaries requires that at least
1013the first and last rows and columns of the \textit{mbathy} array are set to zero.
1014In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to
1015the second one and its first column equal to the last but one (and so too the mask arrays)
1016(see \autoref{fig:LBC_jperio}).
1017
1018% ================================================================
1019% Domain: Initial State (dtatsd & istate)
1020% ================================================================
1021\section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]
1022{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})}
1023\label{sec:DTA_tsd}
1024%-----------------------------------------namtsd-------------------------------------------
1025
1026\nlst{namtsd} 
1027%------------------------------------------------------------------------------------------
1028
1029Options are defined in \ngn{namtsd}.
1030By default, the ocean start from rest (the velocity field is set to zero) and the initialization of temperature and
1031salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
1032\begin{description}
1033\item[\np{ln\_tsd\_init}\forcode{ = .true.}]
1034  use a T and S input files that can be given on the model grid itself or on their native input data grid.
1035  In the latter case,
1036  the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid
1037  (see \autoref{subsec:SBC_iof}).
1038  The information relative to the input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.
1039  The computation is done in the \mdl{dtatsd} module.
1040\item[\np{ln\_tsd\_init}\forcode{ = .false.}]
1041  use constant salinity value of $35.5~psu$ and an analytical profile of temperature
1042  (typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
1043\end{description}
1044
1045\biblio
1046
1047\pindex
1048
1049\end{document}
Note: See TracBrowser for help on using the repository browser.