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

Last change on this file since 10442 was 10442, checked in by nicolasmartin, 2 years ago

Front page edition, cleaning in custom LaTeX commands and add index for single subfile compilation

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