# source:NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_DOM.tex@10419

Last change on this file since 10419 was 10419, checked in by smasson, 22 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

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