source: branches/2012/dev_MERGE_2012/DOC/TexFiles/Chapters/Chap_DOM.tex @ 3764

Last change on this file since 3764 was 3764, checked in by smasson, 9 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012

File size: 53.1 KB
Line 
1% ================================================================
2% Chapter 2 Ñ Space and Time Domain (DOM)
3% ================================================================
4\chapter{Space Domain (DOM) }
5\label{DOM}
6\minitoc
7
8% Missing things:
9%  - istate: description of the initial state   ==> this has to be put elsewhere..
10%                  perhaps in MISC ?  By the way the initialisation of T S and dynamics
11%                  should be put outside of DOM routine (better with TRC staff and off-line
12%                  tracers)
13%  -geo2ocean:  how to switch from geographic to mesh coordinate
14%     - domclo:  closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled
15
16
17\newpage
18$\ $\newline    % force a new ligne
19
20Having defined the continuous equations in Chap.~\ref{PE} and chosen a time
21discretization Chap.~\ref{STP}, we need to choose a discretization on a grid,
22and numerical algorithms. In the present chapter, we provide a general description
23of the staggered grid used in \NEMO, and other information relevant to the main
24directory routines as well as the DOM (DOMain) directory.
25
26$\ $\newline    % force a new ligne
27
28% ================================================================
29% Fundamentals of the Discretisation
30% ================================================================
31\section{Fundamentals of the Discretisation}
32\label{DOM_basics}
33
34% -------------------------------------------------------------------------------------------------------------
35%        Arrangement of Variables
36% -------------------------------------------------------------------------------------------------------------
37\subsection{Arrangement of Variables}
38\label{DOM_cell}
39
40%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
41\begin{figure}[!tb]    \begin{center}
42\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf}
43\caption{ \label{Fig_cell}   
44Arrangement of variables. $t$ indicates scalar points where temperature,
45salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$)
46indicates vector points, and $f$ indicates vorticity points where both relative and
47planetary vorticities are defined}
48\end{center}   \end{figure}
49%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
50
51The numerical techniques used to solve the Primitive Equations in this model are
52based on the traditional, centred second-order finite difference approximation.
53Special attention has been given to the homogeneity of the solution in the three
54space directions. The arrangement of variables is the same in all directions.
55It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector
56points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}).
57This is the generalisation to three dimensions of the well-known ``C'' grid in
58Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and
59planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge
60and the barotropic stream function $\psi$ is defined at horizontal points overlying
61the $\zeta$ and $f$-points.
62
63The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined
64by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$.
65The grid-points are located at integer or integer and a half value of $(i,j,k)$ as
66indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$,
67$f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale
68factors are defined. Each scale factor is defined as the local analytical value
69provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial
70derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and
71$\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity.
72Discrete partial derivatives are formulated by the traditional, centred second order
73finite difference approximation while the scale factors are chosen equal to their
74local analytical value. An important point here is that the partial derivative of the
75scale factors must be evaluated by centred finite difference approximation, not
76from their analytical expression. This preserves the symmetry of the discrete set
77of equations and therefore satisfies many of the continuous properties (see
78Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain
79size: when needed, an area, volume, or the total ocean depth must be evaluated
80as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section).
81
82%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
83\begin{table}[!tb]
84\begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|}
85\hline
86&$i$     & $j$    & $k$     \\ \hline
87& $i+1/2$   & $j$    & $k$    \\ \hline
88& $i$    & $j+1/2$   & $k$    \\ \hline
89& $i$    & $j$    & $k+1/2$   \\ \hline
90& $i+1/2$   & $j+1/2$   & $k$    \\ \hline
91uw & $i+1/2$   & $j$    & $k+1/2$   \\ \hline
92vw & $i$    & $j+1/2$   & $k+1/2$   \\ \hline
93fw & $i+1/2$   & $j+1/2$   & $k+1/2$   \\ \hline
94\end{tabular}
95\caption{ \label{Tab_cell}
96Location of grid-points as a function of integer or integer and a half value of the column,
97line or level. This indexing is only used for the writing of the semi-discrete equation.
98In the code, the indexing uses integer values only and has a reverse direction
99in the vertical (see \S\ref{DOM_Num_Index})}
100\end{center}
101\end{table}
102%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
103
104% -------------------------------------------------------------------------------------------------------------
105%        Vector Invariant Formulation
106% -------------------------------------------------------------------------------------------------------------
107\subsection{Discrete Operators}
108\label{DOM_operators}
109
110Given the values of a variable $q$ at adjacent points, the differencing and
111averaging operators at the midpoint between them are:
112\begin{subequations} \label{Eq_di_mi}
113\begin{align}
114 \delta _i [q]       &\  \    q(i+1/2)  - q(i-1/2)    \\
115 \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2
116\end{align}
117\end{subequations}
118
119Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and
120$k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a
121variable $q$ defined at a $t$-point has its three components defined at $u$-, $v$-
122and $w$-points while its Laplacien is defined at $t$-point. These operators have
123the following discrete forms in the curvilinear $s$-coordinate system:
124\begin{equation} \label{Eq_DOM_grad}
125\nabla q\equiv    \frac{1}{e_{1u} } \delta _{i+1/2 } [q] \;\,\mathbf{i}
126      +  \frac{1}{e_{2v} } \delta _{j+1/2 } [q] \;\,\mathbf{j}
127      +  \frac{1}{e_{3w}} \delta _{k+1/2} [q] \;\,\mathbf{k}
128\end{equation}
129\begin{multline} \label{Eq_DOM_lap}
130\Delta q\equiv \frac{1}{e_{1t}\,e_{2t}\,e_{3t} }
131       \;\left(          \delta_\left[ \frac{e_{2u}\,e_{3u}} {e_{1u}} \;\delta_{i+1/2} [q] \right]
132+                        \delta_\left[ \frac{e_{1v}\,e_{3v}}  {e_{2v}} \;\delta_{j+1/2} [q] \right] \;  \right)      \\
133+\frac{1}{e_{3t}} \delta_k \left[ \frac{1}{e_{3w} }                     \;\delta_{k+1/2} [q] \right]
134\end{multline}
135
136Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ 
137defined at vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$,
138and $f$-points, and its divergence defined at $t$-points:
139\begin{eqnarray}  \label{Eq_DOM_curl}
140 \nabla \times {\rm {\bf A}}\equiv &
141      \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} \\ 
142 +& \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} \\
143 +& \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}
144 \end{eqnarray}
145\begin{equation} \label{Eq_DOM_div}
146\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}\,e_{3t}}\left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
147                                                                                         +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]
148\end{equation}
149
150In the special case of a pure $z$-coordinate system, \eqref{Eq_DOM_lap} and
151\eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor
152becomes a function of the single variable $k$ and thus does not depend on the
153horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to:
154\begin{equation*}
155\nabla \cdot \rm{\bf A}=\frac{1}{e_{1t}\,e_{2t}} \left( \delta_i \left[e_{2u}\,a_1 \right] 
156                                                                              +\delta_j \left[e_{1v}\, a_2 \right]  \right)
157                                                     +\frac{1}{e_{3t}} \delta_k \left[             a_3 \right]
158\end{equation*}
159
160The vertical average over the whole water column denoted by an overbar becomes
161for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):
162\begin{equation} \label{DOM_bar}
163\bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk} 
164      \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }
165\end{equation}
166where $H_q$  is the ocean depth, which is the masked sum of the vertical scale
167factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices,
168and the symbol $k^o$ refers to a summation over all grid points of the same type
169in the direction indicated by the subscript (here $k$).
170
171In continuous form, the following properties are satisfied:
172\begin{equation} \label{Eq_DOM_curl_grad}
173\nabla \times \nabla q ={\rm {\bf {0}}}
174\end{equation}
175\begin{equation} \label{Eq_DOM_div_curl}
176\nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0
177\end{equation}
178
179It is straightforward to demonstrate that these properties are verified locally in
180discrete form as soon as the scalar $q$ is taken at $t$-points and the vector
181\textbf{A} has its components defined at vector points $(u,v,w)$.
182
183Let $a$ and $b$ be two fields defined on the mesh, with value zero inside
184continental area. Using integration by parts it can be shown that the differencing
185operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear
186operators, and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,
187$\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear
188operators, $i.e.$
189\begin{align} 
190\label{DOM_di_adj}
191\sum\limits_i { a_i \;\delta _i \left[ b \right]} 
192   &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\
193\label{DOM_mi_adj}
194\sum\limits_i { a_i \;\overline b^{\,i}} 
195   & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 
196\end{align}
197
198In other words, the adjoint of the differencing and averaging operators are
199$\delta_i^*=\delta_{i+1/2}$ and
200${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively.
201These two properties will be used extensively in the Appendix~\ref{Apdx_C} to
202demonstrate integral conservative properties of the discrete formulation chosen.
203
204% -------------------------------------------------------------------------------------------------------------
205%        Numerical Indexing
206% -------------------------------------------------------------------------------------------------------------
207\subsection{Numerical Indexing}
208\label{DOM_Num_Index}
209
210%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
211\begin{figure}[!tb]  \begin{center}
212\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf}
213\caption{   \label{Fig_index_hor}   
214Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates
215the cell in which variables contained in arrays have the same $i$- and $j$-indices}
216\end{center}   \end{figure}
217%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
218
219The array representation used in the \textsc{Fortran} code requires an integer
220indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is
221associated with the use of integer values for $t$-points and both integer and
222integer and a half values for all the other points. Therefore a specific integer
223indexing must be defined for points other than $t$-points ($i.e.$ velocity and
224vorticity grid-points). Furthermore, the direction of the vertical indexing has
225been changed so that the surface level is at $k=1$.
226
227% -----------------------------------
228%        Horizontal Indexing
229% -----------------------------------
230\subsubsection{Horizontal Indexing}
231\label{DOM_Num_Index_hor}
232
233The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}.
234For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point
235(northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}).
236A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices.
237
238% -----------------------------------
239%        Vertical indexing
240% -----------------------------------
241\subsubsection{Vertical Indexing}
242\label{DOM_Num_Index_vertical}
243
244In the vertical, the chosen indexing requires special attention since the
245$k$-axis is re-orientated downward in the \textsc{Fortran} code compared
246to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}.
247The sea surface corresponds to the $w$-level $k=1$ which is the same index
248as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)
249either corresponds to the ocean floor or is inside the bathymetry while the last
250$t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that
251for an increasing $k$ index, a $w$-point and the $t$-point just below have the
252same $k$ index, in opposition to what is done in the horizontal plane where
253it is the $t$-point and the nearest velocity points in the direction of the horizontal
254axis that have the same $i$ or $j$ index (compare the dashed area in
255Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are
256chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} 
257code \emph{before all the vertical derivatives} of the discrete equations given in
258this documentation.
259
260%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
261\begin{figure}[!pt]    \begin{center}
262\includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf}
263\caption{ \label{Fig_index_vert}     
264Vertical integer indexing used in the \textsc{Fortran } code. Note that
265the $k$-axis is orientated downward. The dashed area indicates the cell in
266which variables contained in arrays have the same $k$-index.}
267\end{center}   \end{figure}
268%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
269
270% -----------------------------------
271%        Domain Size
272% -----------------------------------
273\subsubsection{Domain Size}
274\label{DOM_size}
275
276The total size of the computational domain is set by the parameters \jp{jpiglo},
277\jp{jpjglo} and \jp{jpk} in the $i$, $j$ and $k$ directions respectively. They are
278given as parameters in the \mdl{par\_oce} module\footnote{When a specific
279configuration is used (ORCA2 global ocean, etc...) the parameter are actually
280defined in additional files introduced by \mdl{par\_oce} module via CPP
281\textit{include} command. For example, ORCA2 parameters are set in
282\textit{par\_ORCA\_R2.h90} file}. The use of parameters rather than variables
283(together with dynamic allocation of arrays) was chosen because it ensured that
284the compiler would optimize the executable code efficiently, especially on vector
285machines (optimization may be less efficient when the problem size is unknown
286at the time of compilation). Nevertheless, it is possible to set up the code with full
287dynamical allocation by using the AGRIF packaged \citep{Debreu_al_CG2008}.
288%
289\gmcomment{  add the following ref
290\colorbox{yellow}{(ref part of the doc)} } 
291%
292Note that are other parameters in \mdl{par\_oce} that refer to the domain size.
293The two parameters $jpidta$ and $jpjdta$ may be larger than $jpiglo$, $jpjglo$ 
294when the user wants to use only a sub-region of a given configuration. This is
295the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of
296the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters
297$jpi$ and $jpj$ refer to the size of each processor subdomain when the code is
298run in parallel using domain decomposition (\key{mpp\_mpi} defined, see
299\S\ref{LBC_mpp}).
300
301
302$\ $\newline    % force a new ligne
303
304% ================================================================
305% Domain: Horizontal Grid (mesh)
306% ================================================================
307\section  [Domain: Horizontal Grid (mesh) (\textit{domhgr})]               
308      {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} }
309\label{DOM_hgr}
310
311% -------------------------------------------------------------------------------------------------------------
312%        Coordinates and scale factors
313% -------------------------------------------------------------------------------------------------------------
314\subsection{Coordinates and scale factors}
315\label{DOM_hgr_coord_e}
316
317The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined
318by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
319The grid-points are located at integer or integer and a half values of as indicated
320in Table~\ref{Tab_cell}. The associated scale factors are defined using the
321analytical first derivative of the transformation \eqref{Eq_scale_factors}. These
322definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which
323provide the horizontal and vertical meshes, respectively. This section deals with
324the horizontal mesh parameters.
325
326In a horizontal plane, the location of all the model grid points is defined from the
327analytical expressions of the longitude $\lambda$ and  latitude $\varphi$ as a
328function of  $(i,j)$. The horizontal scale factors are calculated using
329\eqref{Eq_scale_factors}. For example, when the longitude and latitude are
330function of a single value ($i$ and $j$, respectively) (geographical configuration
331of the mesh), the horizontal mesh definition reduces to define the wanted
332$\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the
333\mdl{domhgr} module. The model computes the grid-point positions and scale
334factors in the horizontal plane as follows:
335\begin{flalign*}
336\lambda_t &\equiv \text{glamt}= \lambda(i)     & \varphi_t &\equiv \text{gphit} = \varphi(j)\\
337\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\
338\lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\
339\lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2)
340\end{flalign*}
341\begin{flalign*}
342e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |&
343e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\
344e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |&
345e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\
346e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j+1/2)  |&
347e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\
348e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2)  |&
349e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|
350\end{flalign*}
351where the last letter of each computational name indicates the grid point
352considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with
353all universal constants). Note that the horizontal position of and scale factors
354at $w$-points are exactly equal to those of $t$-points, thus no specific arrays
355are defined at $w$-points.
356
357Note that the definition of the scale factors ($i.e.$ as the analytical first derivative
358of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is
359specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined
360locally at a $t$-point, whereas many other models on a C grid choose to define
361such a scale factor as the distance between the $U$-points on each side of the
362$t$-point. Relying on an analytical transformation has two advantages: firstly, there
363is no ambiguity in the scale factors appearing in the discrete equations, since they
364are first introduced in the continuous equations; secondly, analytical transformations
365encourage good practice by the definition of smoothly varying grids (rather than
366allowing the user to set arbitrary jumps in thickness between adjacent layers)
367\citep{Treguier1996}. An example of the effect of such a choice is shown in
368Fig.~\ref{Fig_zgr_e3}.
369%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
370\begin{figure}[!t]     \begin{center}
371\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf}
372\caption{ \label{Fig_zgr_e3}   
373Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
374and (b) analytically derived grid-point position and scale factors.
375For both grids here,  the same $w$-point depth has been chosen but in (a) the
376$t$-points are set half way between $w$-points while in (b) they are defined from
377an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$.
378Note the resulting difference between the value of the grid-size $\Delta_k$ and
379those of the scale factor $e_k$. }
380\end{center}   \end{figure}
381%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
382
383% -------------------------------------------------------------------------------------------------------------
384%        Choice of horizontal grid
385% -------------------------------------------------------------------------------------------------------------
386\subsection{Choice of horizontal grid}
387\label{DOM_hgr_msh_choice}
388
389The user has three options available in defining a horizontal grid, which involve
390the parameter $jphgr\_mesh$ of the \mdl{par\_oce} module.
391\begin{description}
392\item[\jp{jphgr\_mesh}=0]  The most general curvilinear orthogonal grids.
393The coordinates and their first derivatives with respect to $i$ and $j$ are provided
394in a input file (\ifile{coordinates}), read in \rou{hgr\_read} subroutine of the domhgr module.
395\item[\jp{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below).
396For other analytical grids, the \mdl{domhgr} module must be modified by the user.
397\end{description}
398
399There are two simple cases of geographical grids on the sphere. With
400\jp{jphgr\_mesh}=1, the grid (expressed in degrees) is regular in space,
401with grid sizes specified by parameters \pp{ppe1\_deg} and \pp{ppe2\_deg},
402respectively. Such a geographical grid can be very anisotropic at high latitudes
403because of the convergence of meridians (the zonal scale factors $e_1$ 
404become much smaller than the meridional scale factors $e_2$). The Mercator
405grid (\jp{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale
406factors in the same way as the zonal ones. In this case, meridional scale factors
407and latitudes are calculated analytically using the formulae appropriate for
408a Mercator projection, based on \pp{ppe1\_deg} which is a reference grid spacing
409at the equator (this applies even when the geographical equator is situated outside
410the model domain).
411%%%
412\gmcomment{ give here the analytical expression of the Mercator mesh}
413%%%
414In these two cases (\jp{jphgr\_mesh}=1 or 4), the grid position is defined by the
415longitude and latitude of the south-westernmost point (\pp{ppglamt0} 
416and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide
417an approximate starting latitude: the real latitude will be recalculated analytically,
418in order to ensure that the equator corresponds to line passing through $t$-
419and $u$-points. 
420
421Rectangular grids ignoring the spherical geometry are defined with
422\jp{jphgr\_mesh} = 2, 3, 5. The domain is either an $f$-plane (\jp{jphgr\_mesh} = 2,
423Coriolis factor is constant) or a beta-plane (\jp{jphgr\_mesh} = 3, the Coriolis factor
424is linear in the $j$-direction). The grid size is uniform in meter in each direction,
425and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively.
426The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero
427with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers,
428and the second $t$-point corresponds to coordinate $gphit=0$. The input
429parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference
430latitude for computation of the Coriolis parameter. In the case of the beta plane,
431\pp{ppgphi0} corresponds to the center of the domain. Finally, the special case
432\jp{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the
433GYRE configuration, representing a classical mid-latitude double gyre system.
434The rotation allows us to maximize the jet length relative to the gyre areas
435(and the number of grid points).
436
437The choice of the grid must be consistent with the boundary conditions specified
438by the parameter \jp{jperio} (see {\S\ref{LBC}).
439
440% -------------------------------------------------------------------------------------------------------------
441%        Grid files
442% -------------------------------------------------------------------------------------------------------------
443\subsection{Output Grid files}
444\label{DOM_hgr_files}
445
446All the arrays relating to a particular ocean model configuration (grid-point
447position, scale factors, masks) can be saved in files if $\np{nn\_msh} \not= 0$ 
448(namelist parameter). This can be particularly useful for plots and off-line
449diagnostics. In some cases, the user may choose to make a local modification
450of a scale factor in the code. This is the case in global configurations when
451restricting the width of a specific strait (usually a one-grid-point strait that
452happens to be too wide due to insufficient model resolution). An example
453is Gibraltar Strait in the ORCA2 configuration. When such modifications are done,
454the output grid written when $\np{nn\_msh} \not=0$ is no more equal to the input grid.
455
456$\ $\newline    % force a new ligne
457
458% ================================================================
459% Domain: Vertical Grid (domzgr)
460% ================================================================
461\section  [Domain: Vertical Grid (\textit{domzgr})]
462      {Domain: Vertical Grid \small{(\mdl{domzgr} module)} }
463\label{DOM_zgr}
464%-----------------------------------------nam_zgr & namdom-------------------------------------------
465\namdisplay{namzgr} 
466\namdisplay{namdom} 
467%-------------------------------------------------------------------------------------------------------------
468
469In the vertical, the model mesh is determined by four things:
470(1) the bathymetry given in meters ;
471(2) the number of levels of the model (\jp{jpk}) ;
472(3) the analytical transformation $z(i,j,k)$ and the vertical scale factors
473(derivatives of the transformation) ;
474and (4) the masking system, $i.e.$ the number of wet model levels at each
475$(i,j)$ column of points.
476
477%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
478\begin{figure}[!tb]    \begin{center}
479\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf}
480\caption{  \label{Fig_z_zps_s_sps}   
481The ocean bottom as seen by the model:
482(a) $z$-coordinate with full step,
483(b) $z$-coordinate with partial step,
484(c) $s$-coordinate: terrain following representation,
485(d) hybrid $s-z$ coordinate,
486(e) hybrid $s-z$ coordinate with partial step, and
487(f) same as (e) but with variable volume associated with the non-linear free surface.
488Note that the variable volume option (\key{vvl}) can be used with any of the
4895 coordinates (a) to (e).}
490\end{center}   \end{figure}
491%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
492
493The choice of a vertical coordinate, even if it is made through a namelist parameter,
494must be done once of all at the beginning of an experiment. It is not intended as an
495option which can be enabled or disabled in the middle of an experiment. Three main
496choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step
497bathymetry (\np{ln\_zco}~=~true), $z$-coordinate with partial step bathymetry
498(\np{ln\_zps}~=~true), or generalized, $s$-coordinate (\np{ln\_sco}~=~true).
499Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate
500(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable
501volume option \key{vvl} ($i.e.$ non-linear free surface), the coordinate follow the
502time-variation of the free surface so that the transformation is time dependent:
503$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step
504bathymetry or $s$-coordinate (hybrid and partial step coordinates have not
505yet been tested in NEMO v2.3).
506
507Contrary to the horizontal grid, the vertical grid is computed in the code and no
508provision is made for reading it from a file. The only input file is the bathymetry
509(in meters) (\ifile{bathy\_meter})
510\footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the
511\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point
512in each water column is by-passed}.
513After reading the bathymetry, the algorithm for vertical grid definition differs
514between the different options:
515\begin{description}
516\item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$.
517\item[\textit{zps}] set a reference coordinate transformation $z_0 (k)$, and
518calculate the thickness of the deepest level at each $(i,j)$ point using the
519bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
520\item[\textit{sco}] smooth the bathymetry to fulfil the hydrostatic consistency
521criteria and set the three-dimensional transformation.
522\item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfil the hydrostatic
523consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and
524possibly introduce masking of extra land points to better fit the original bathymetry file
525\end{description}
526%%%
527\gmcomment{   add the description of the smoothing:  envelop topography...}
528%%%
529
530The arrays describing the grid point depths and vertical scale factors
531are three dimensional arrays $(i,j,k)$ even in the case of $z$-coordinate with
532full step bottom topography. In non-linear free surface (\key{vvl}), their knowledge
533is required at \textit{before}, \textit{now} and \textit{after} time step, while they
534do not vary in time in linear free surface case.
535To improve the code readability while providing this flexibility, the vertical coordinate
536and scale factors are defined as functions of
537$(i,j,k)$ with "fs" as prefix (examples: \textit{fse3t\_b, fse3t\_n, fse3t\_a,} 
538for the  \textit{before}, \textit{now} and \textit{after} scale factors at $t$-point)
539that can be either three different arrays when \key{vvl} is defined, or a single fixed arrays.
540These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory.
541They are used throughout the code, and replaced by the corresponding arrays at
542the time of pre-processing (CPP capability).
543
544% -------------------------------------------------------------------------------------------------------------
545%        Meter Bathymetry
546% -------------------------------------------------------------------------------------------------------------
547\subsection{Meter Bathymetry}
548\label{DOM_bathy}
549
550Three options are possible for defining the bathymetry, according to the
551namelist variable \np{nn\_bathy}:
552\begin{description}
553\item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ 
554is given by the coordinate transformation. The domain can either be a closed
555basin or a periodic channel depending on the parameter \jp{jperio}.
556\item[\np{nn\_bathy} = -1] a domain with a bump of topography one third of the
557domain width at the central latitude. This is meant for the "EEL-R5" configuration,
558a periodic or open boundary channel with a seamount.
559\item[\np{nn\_bathy} = 1] read a bathymetry. The \ifile{bathy\_meter} file (Netcdf format)
560provides the ocean depth (positive, in meters) at each grid point of the model grid.
561The bathymetry is usually built by interpolating a standard bathymetry product
562($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also
563defines the coastline: where the bathymetry is zero, no model levels are defined
564(all levels are masked).
565\end{description}
566
567When a global ocean is coupled to an atmospheric model it is better to represent
568all large water bodies (e.g, great lakes, Caspian sea...) even if the model
569resolution does not allow their communication with the rest of the ocean.
570This is unnecessary when the ocean is forced by fixed atmospheric conditions,
571so these seas can be removed from the ocean domain. The user has the option
572to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the
573code has to be adapted to the user's configuration.
574
575% -------------------------------------------------------------------------------------------------------------
576%        z-coordinate  and reference coordinate transformation
577% -------------------------------------------------------------------------------------------------------------
578\subsection[$z$-coordinate (\np{ln\_zco}]
579        {$z$-coordinate (\np{ln\_zco}=true) and reference coordinate}
580\label{DOM_zco}
581
582%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
583\begin{figure}[!tb]    \begin{center}
584\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf}
585\caption{ \label{Fig_zgr}   
586Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for
587(a) T-point depth and (b) the associated scale factor as computed
588from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.}
589\end{center}   \end{figure}
590%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
591
592The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ 
593and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on
594Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the
595ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the
596additional $t$-point at $jk=jpk$ is below the sea floor and is not used.
597The vertical location of $w$- and $t$-levels is defined from the analytic expression
598of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the
599vertical scale factors. The user must provide the analytical expression of both
600$z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} 
601through statement functions, using parameters provided in the \textit{par\_oce.h90} file.
602
603It is possible to define a simple regular vertical grid by giving zero stretching (\pp{ppacr=0}).
604In that case, the parameters \jp{jpk} (number of $w$-levels) and \pp{pphmax} 
605(total ocean depth in meters) fully define the grid.
606
607For climate-related studies it is often desirable to concentrate the vertical resolution
608near the ocean surface. The following function is proposed as a standard for a
609$z$-coordinate (with either full or partial steps):
610\begin{equation} \label{DOM_zgr_ana}
611\begin{split}
612 z_0 (k)    &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 
613 e_3^0 (k)  &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right|
614\end{split}
615\end{equation}
616where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an
617expression allows us to define a nearly uniform vertical location of levels at the
618ocean top and bottom with a smooth hyperbolic tangent transition in between
619(Fig.~\ref{Fig_zgr}).
620
621The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the
622surface (bottom) layers and a depth which varies from 0 at the sea surface to a
623minimum of $-5000~m$. This leads to the following conditions:
624\begin{equation} \label{DOM_zgr_coef}
625\begin{split}
626 e_3 (1+1/2)      &=10. \\ 
627 e_3 (jpk-1/2) &=500. \\ 
628 z(1)       &=0. \\ 
629 z(jpk)        &=-5000. \\ 
630\end{split}
631\end{equation}
632
633With the choice of the stretching $h_{cr} =3$ and the number of levels
634\jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in
635\eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is
636satisfied, through an optimisation procedure using a bisection method. For the first
637standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$,
638$h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and
639scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and
640given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters
641\pp{ppsur}, \pp{ppa0}, \pp{ppa1}, \pp{ppkth} in the parameter file \mdl{par\_oce}.
642
643Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is
644possible to recalculate them. In that case the user sets
645\pp{ppsur}=\pp{ppa0}=\pp{ppa1}=\pp{pp\_to\_be\_computed}, in \mdl{par\_oce},
646and specifies instead the four following parameters:
647\begin{itemize}
648\item    \pp{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger
649\pp{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual.
650\item    \pp{ppkth}=$h_{th} $: is approximately the model level at which maximum
651stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
652\item    \pp{ppdzmin}: minimum thickness for the top layer (in meters)
653\item    \pp{pphmax}: total depth of the ocean (meters).
654\end{itemize}
655As an example, for the $45$ layers used in the DRAKKAR configuration those
656parameters are: \jp{jpk}=46, \pp{ppacr}=9, \pp{ppkth}=23.563, \pp{ppdzmin}=6m,
657\pp{pphmax}=5750m.
658
659%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
660\begin{table}     \begin{center} \begin{tabular}{c||r|r|r|r}
661\hline
662\textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline
663&  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline
664&  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline
665&  \textbf{25.00} &    20.00 &   \textbf{ 10.00} &     10.00 \\   \hline
666&  \textbf{35.01} &    30.00 &   \textbf{ 10.01} &     10.00 \\   \hline
667&  \textbf{45.01} &    40.01 &   \textbf{ 10.01} &  10.01 \\   \hline
668&  \textbf{55.03} &    50.02 &   \textbf{ 10.02} &     10.02 \\   \hline
669&  \textbf{65.06} &    60.04 &   \textbf{ 10.04} &  10.03 \\   \hline
670&  \textbf{75.13} &    70.09 &   \textbf{ 10.09} &  10.06 \\   \hline
671&  \textbf{85.25} &    80.18 &   \textbf{ 10.17} &  10.12 \\   \hline
67210 &  \textbf{95.49} &    90.35 &   \textbf{ 10.33} &  10.24 \\   \hline
67311 &  \textbf{105.97}   &   100.69 &   \textbf{ 10.65} &  10.47 \\   \hline
67412 &  \textbf{116.90}   &   111.36 &   \textbf{ 11.27} &  10.91 \\   \hline
67513 &  \textbf{128.70}   &   122.65 &   \textbf{ 12.47} &  11.77 \\   \hline
67614 &  \textbf{142.20}   &   135.16 &   \textbf{ 14.78} &  13.43 \\   \hline
67715 &  \textbf{158.96}   &   150.03 &   \textbf{ 19.23} &  16.65 \\   \hline
67816 &  \textbf{181.96}   &   169.42 &   \textbf{ 27.66} &  22.78 \\   \hline
67917 &  \textbf{216.65}   &   197.37 &   \textbf{ 43.26} &  34.30 \\ \hline
68018 &  \textbf{272.48}   &   241.13 &   \textbf{ 70.88} &  55.21 \\ \hline
68119 &  \textbf{364.30}   &   312.74 &   \textbf{116.11} &  90.99 \\ \hline
68220 &  \textbf{511.53}   &   429.72 &   \textbf{181.55} &    146.43 \\ \hline
68321 &  \textbf{732.20}   &   611.89 &   \textbf{261.03} &    220.35 \\ \hline
68422 &  \textbf{1033.22}&  872.87 &   \textbf{339.39} &    301.42 \\ \hline
68523 &  \textbf{1405.70}& 1211.59 & \textbf{402.26} &   373.31 \\ \hline
68624 &  \textbf{1830.89}& 1612.98 & \textbf{444.87} &   426.00 \\ \hline
68725 &  \textbf{2289.77}& 2057.13 & \textbf{470.55} &   459.47 \\ \hline
68826 &  \textbf{2768.24}& 2527.22 & \textbf{484.95} &   478.83 \\ \hline
68927 &  \textbf{3257.48}& 3011.90 & \textbf{492.70} &   489.44 \\ \hline
69028 &  \textbf{3752.44}& 3504.46 & \textbf{496.78} &   495.07 \\ \hline
69129 &  \textbf{4250.40}& 4001.16 & \textbf{498.90} &   498.02 \\ \hline
69230 &  \textbf{4749.91}& 4500.02 & \textbf{500.00} &   499.54 \\ \hline
69331 &  \textbf{5250.23}& 5000.00 &   \textbf{500.56} & 500.33 \\ \hline
694\end{tabular} \end{center} 
695\caption{ \label{Tab_orca_zgr}   
696Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed
697from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}}
698\end{table}
699%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
700
701% -------------------------------------------------------------------------------------------------------------
702%        z-coordinate with partial step
703% -------------------------------------------------------------------------------------------------------------
704\subsection   [$z$-coordinate with partial step (\np{ln\_zps})]
705         {$z$-coordinate with partial step (\np{ln\_zps}=.true.)}
706\label{DOM_zps}
707%--------------------------------------------namdom-------------------------------------------------------
708\namdisplay{namdom} 
709%--------------------------------------------------------------------------------------------------------------
710
711In $z$-coordinate partial step, the depths of the model levels are defined by the
712reference analytical function $z_0 (k)$ as described in the previous
713section, \emph{except} in the bottom layer. The thickness of the bottom layer is
714allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a
715better representation of the bathymetry, especially in the case of small
716slopes (where the bathymetry varies by less than one level thickness from
717one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been
718defined in the absence of bathymetry. With partial steps, layers from 1 to
719\jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1)
720is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the
721maximum thickness allowed is $2*e_{3t}(jpk-1)$. This has to be kept in mind when
722specifying the maximum depth \pp{pphmax} in partial steps: for example, with
723\pp{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth
724allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$).
725Two variables in the namdom namelist are used to define the partial step
726vertical grid. The mimimum water thickness (in meters) allowed for a cell
727partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min} 
728(thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction,
729usually 10\%, of the default thickness $e_{3t}(jk)$).
730
731 \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }
732
733% -------------------------------------------------------------------------------------------------------------
734%        s-coordinate
735% -------------------------------------------------------------------------------------------------------------
736\subsection   [$s$-coordinate (\np{ln\_sco})]
737           {$s$-coordinate (\np{ln\_sco}=true)}
738\label{DOM_sco}
739%------------------------------------------nam_zgr_sco---------------------------------------------------
740\namdisplay{namzgr_sco} 
741%--------------------------------------------------------------------------------------------------------------
742In $s$-coordinate (\np{ln\_sco}~=~true), the depth and thickness of the model
743levels are defined from the product of a depth field and either a stretching
744function or its derivative, respectively:
745
746\begin{equation} \label{DOM_sco_ana}
747\begin{split}
748 z(k)       &= h(i,j) \; z_0(k)  \\
749 e_3(k)  &= h(i,j) \; z_0'(k)
750\end{split}
751\end{equation}
752
753where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $t$-point
754location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea
755surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean
756depth, since a mixed step-like and bottom-following representation of the
757topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e) or an envelop bathymetry can be defined (Fig.~\ref{Fig_z_zps_s_sps}f).
758The namelist parameter \np{rn\_rmax} determines the slope at which the terrain-following coordinate intersects the sea bed and becomes a pseudo z-coordinate. The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
759
760Options for stretching the coordinate are provided as examples, but care must be taken to ensure that the vertical stretch used is appropriate for the application.
761
762The original default NEMO s-coordinate stretching is available if neither of the other options are specified as true (\np{ln\_sco\_SH94}~=~false and \np{ln\_sco\_SF12}~=~false.) This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
763
764\begin{equation}
765  z = s_{min}+C\left(s\right)\left(H-s_{min}\right)
766  \label{eq:SH94_1}
767\end{equation}
768
769where $s_{min}$ is the depth at which the s-coordinate stretching starts and allows a z-coordinate to placed on top of the stretched coordinate, and z is the depth (negative down from the asea surface).
770
771\begin{equation}
772  s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1
773  \label{eq:s}
774\end{equation}
775
776\begin{equation} \label{DOM_sco_function}
777\begin{split}
778C(s)  &\frac{ \left[   \tanh{ \left( \theta \, (s+b) \right)} 
779               - \tanh{ \left\theta \, b      \right)}  \right]}
780            {2\;\sinh \left( \theta \right)}
781\end{split}
782\end{equation}
783
784A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} stretching (\np{ln\_sco\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling:
785
786\begin{equation}
787  C\left(s\right) =   \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} +      \\
788  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)}
789  \label{eq:SH94_2}
790\end{equation}
791
792%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
793\begin{figure}[!ht]    \begin{center}
794\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
795\caption{  \label{Fig_sco_function}   
796Examples of the stretching function applied to a seamount; from left to right:
797surface, surface and bottom, and bottom intensified resolutions}
798\end{center}   \end{figure}
799%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
800
801where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from pure $\sigma$ to the stretched coordinate,  and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb}) are the surface and
802bottom control parameters such that $0\leqslant \theta \leqslant 20$, and
803$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom
804increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
805
806Another example has been provided at version 3.5 (\np{ln\_sco\_SF12}) that allows a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}. In this case the a stretching function $\gamma$ is defined such that:
807
808\begin{equation}
809z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1
810\label{eq:z}
811\end{equation}
812
813The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
814
815\begin{equation} \label{DOM_gamma_deriv}
816\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)
817\end{equation}
818
819Where:
820\begin{equation} \label{DOM_gamma}
821f\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} 
822\end{equation}
823
824This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs}) and bottom depths. The bottom cell depth in this example is given as a function of water depth:
825
826\begin{equation} \label{DOM_zb}
827Z_b= h a + b
828\end{equation}
829
830where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively.
831
832%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
833\begin{figure}[!ht]
834   \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/FIG_DOM_compare_coordinates_surface.pdf}
835        \caption{A comparison of the \citet{Song_Haidvogel_JCP94} $S$-coordinate (solid lines), a 50 level $Z$-coordinate (contoured surfaces) and the \citet{Siddorn_Furner_OM12} $S$-coordinate (dashed lines) in the surface 100m for a idealised bathymetry that goes from 50m to 5500m depth. For clarity every third coordinate surface is shown.}
836    \label{fig_compare_coordinates_surface}
837\end{figure}
838%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
839
840This gives a smooth analytical stretching in computational space that is constrained to given specified surface and bottom grid cell thicknesses in real space. This is not to be confused with the hybrid schemes that superimpose geopotential coordinates on terrain following coordinates thus creating a non-analytical vertical coordinate that therefore may suffer from large gradients in the vertical resolutions. This stretching is less straightforward to implement than the \citet{Song_Haidvogel_JCP94} stretching, but has the advantage of resolving diurnal processes in deep water and has generally flatter slopes.
841
842As with the \citet{Song_Haidvogel_JCP94} stretching the stretch is only applied at depths greater than the critical depth $h_c$. In this example two options are available in depths shallower than $h_c$, with pure sigma being applied if the \np{ln\_sigcrit} is true and pure z-coordinates if it is false (the z-coordinate being equal to the depths of the stretched coordinate at $h_c$.
843
844Minimising the horizontal slope of the vertical coordinate is important in terrain-following systems as large slopes lead to hydrostatic consistency. A hydrostatic consistency parameter diagnostic following \citet{Haney1991} has been implemented, and is output as part of the model mesh file at the start of the run.
845
846% -------------------------------------------------------------------------------------------------------------
847%        z*- or s*-coordinate
848% -------------------------------------------------------------------------------------------------------------
849\subsection{$z^*$- or $s^*$-coordinate (add \key{vvl}) }
850\label{DOM_zgr_vvl}
851
852This option is described in the Report by Levier \textit{et al.} (2007), available on
853the \NEMO web site.
854
855%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
856
857% -------------------------------------------------------------------------------------------------------------
858%        level bathymetry and mask
859% -------------------------------------------------------------------------------------------------------------
860\subsection{level bathymetry and mask}
861\label{DOM_msk}
862
863Whatever the vertical coordinate used, the model offers the possibility of
864representing the bottom topography with steps that follow the face of the
865model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of
866the steps in the horizontal is defined in a 2D integer array, mbathy, which
867gives the number of ocean levels ($i.e.$ those that are not masked) at each
868$t$-point. mbathy is computed from the meter bathymetry using the definiton of
869gdept as the number of $t$-points which gdept $\leq$ bathy.
870
871Modifications of the model bathymetry are performed in the \textit{bat\_ctl} 
872routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points
873that do not communicate with another ocean point at the same level are eliminated.
874
875From the \textit{mbathy} array, the mask fields are defined as follows:
876\begin{align*}
877tmask(i,j,k) &= \begin{cases}   \; 1&   \text{ if $k\leq mbathy(i,j)$  }    \\
878                                                \; 0&   \text{ if $k\leq mbathy(i,j)$  }    \end{cases}     \\
879umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
880vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\
881fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
882                   & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k)
883\end{align*}
884
885Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with
886the numerical indexing used (\S~\ref{DOM_Num_Index}). Moreover, the
887specification of closed lateral boundaries requires that at least the first and last
888rows and columns of the \textit{mbathy} array are set to zero. In the particular
889case of an east-west cyclical boundary condition, \textit{mbathy} has its last
890column equal to the second one and its first column equal to the last but one
891(and so too the mask arrays) (see \S~\ref{LBC_jperio}).
892
893%%%
894\gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}}
895%%%
896
897% ================================================================
898% Domain: Initial State (dtatsd & istate)
899% ================================================================
900\section  [Domain: Initial State (\textit{istate and dtatsd})]
901      {Domain: Initial State \small{(\mdl{istate} and \mdl{dtatsd} modules)} }
902\label{DTA_tsd}
903%-----------------------------------------namtsd-------------------------------------------
904\namdisplay{namtsd} 
905%------------------------------------------------------------------------------------------
906
907By default, the ocean start from rest (the velocity field is set to zero) and the initialization of
908temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
909\begin{description}
910\item[ln\_tsd\_init = .true.]  use a T and S input files that can be given on the model grid itself or
911on their native input data grid. In the latter case, the data will be interpolated on-the-fly both in the
912horizontal and the vertical to the model grid (see \S~\ref{SBC_iof}). The information relative to the
913input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.
914The computation is done in the \mdl{dtatsd} module.
915\item[ln\_tsd\_init = .false.] use constant salinity value of 35.5 psu and an analytical profile of temperature
916(typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
917\end{description}
Note: See TracBrowser for help on using the repository browser.