source: trunk/DOC/TexFiles/Chapters/Chap_DOM.tex @ 6497

Last change on this file since 6497 was 6497, checked in by gm, 5 years ago

#1720 - trunk: add Casimir tidal parameterization

File size: 54.0 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 lign
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{eqnarray} \label{Eq_DOM_div}
146\nabla \cdot \rm{\bf A} \equiv 
147    \frac{1}{e_{1t}\,e_{2t}\,e_{3t}} \left( \delta_i \left[e_{2u}\,e_{3u}\,a_1 \right]
148                                           +\delta_j \left[e_{1v}\,e_{3v}\,a_2 \right] \right)+\frac{1}{e_{3t} }\delta_k \left[a_3 \right]
149\end{eqnarray}
150
151The vertical average over the whole water column denoted by an overbar becomes
152for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area):
153\begin{equation} \label{DOM_bar}
154\bar q   =         \frac{1}{H}    \int_{k^b}^{k^o} {q\;e_{3q} \,dk} 
155      \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} }
156\end{equation}
157where $H_q$  is the ocean depth, which is the masked sum of the vertical scale
158factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices,
159and the symbol $k^o$ refers to a summation over all grid points of the same type
160in the direction indicated by the subscript (here $k$).
161
162In continuous form, the following properties are satisfied:
163\begin{equation} \label{Eq_DOM_curl_grad}
164\nabla \times \nabla q ={\rm {\bf {0}}}
165\end{equation}
166\begin{equation} \label{Eq_DOM_div_curl}
167\nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0
168\end{equation}
169
170It is straightforward to demonstrate that these properties are verified locally in
171discrete form as soon as the scalar $q$ is taken at $t$-points and the vector
172\textbf{A} has its components defined at vector points $(u,v,w)$.
173
174Let $a$ and $b$ be two fields defined on the mesh, with value zero inside
175continental area. Using integration by parts it can be shown that the differencing
176operators ($\delta_i$, $\delta_j$ and $\delta_k$) are skew-symmetric linear operators,
177and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$,
178$\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear
179operators, $i.e.$
180\begin{align} 
181\label{DOM_di_adj}
182\sum\limits_i { a_i \;\delta _i \left[ b \right]} 
183   &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} }      \\
184\label{DOM_mi_adj}
185\sum\limits_i { a_i \;\overline b^{\,i}} 
186   & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } 
187\end{align}
188
189In other words, the adjoint of the differencing and averaging operators are
190$\delta_i^*=\delta_{i+1/2}$ and
191${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively.
192These two properties will be used extensively in the Appendix~\ref{Apdx_C} to
193demonstrate integral conservative properties of the discrete formulation chosen.
194
195% -------------------------------------------------------------------------------------------------------------
196%        Numerical Indexing
197% -------------------------------------------------------------------------------------------------------------
198\subsection{Numerical Indexing}
199\label{DOM_Num_Index}
200
201%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
202\begin{figure}[!tb]  \begin{center}
203\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf}
204\caption{   \label{Fig_index_hor}   
205Horizontal integer indexing used in the \textsc{Fortran} code. The dashed area indicates
206the cell in which variables contained in arrays have the same $i$- and $j$-indices}
207\end{center}   \end{figure}
208%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
209
210The array representation used in the \textsc{Fortran} code requires an integer
211indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is
212associated with the use of integer values for $t$-points and both integer and
213integer and a half values for all the other points. Therefore a specific integer
214indexing must be defined for points other than $t$-points ($i.e.$ velocity and
215vorticity grid-points). Furthermore, the direction of the vertical indexing has
216been changed so that the surface level is at $k=1$.
217
218% -----------------------------------
219%        Horizontal Indexing
220% -----------------------------------
221\subsubsection{Horizontal Indexing}
222\label{DOM_Num_Index_hor}
223
224The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}.
225For an increasing $i$ index ($j$ index), the $t$-point and the eastward $u$-point
226(northward $v$-point) have the same index (see the dashed area in Fig.\ref{Fig_index_hor}).
227A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices.
228
229% -----------------------------------
230%        Vertical indexing
231% -----------------------------------
232\subsubsection{Vertical Indexing}
233\label{DOM_Num_Index_vertical}
234
235In the vertical, the chosen indexing requires special attention since the
236$k$-axis is re-orientated downward in the \textsc{Fortran} code compared
237to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}.
238The sea surface corresponds to the $w$-level $k=1$ which is the same index
239as $t$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$)
240either corresponds to the ocean floor or is inside the bathymetry while the last
241$t$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that
242for an increasing $k$ index, a $w$-point and the $t$-point just below have the
243same $k$ index, in opposition to what is done in the horizontal plane where
244it is the $t$-point and the nearest velocity points in the direction of the horizontal
245axis that have the same $i$ or $j$ index (compare the dashed area in
246Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are
247chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} 
248code \emph{before all the vertical derivatives} of the discrete equations given in
249this documentation.
250
251%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
252\begin{figure}[!pt]    \begin{center}
253\includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf}
254\caption{ \label{Fig_index_vert}     
255Vertical integer indexing used in the \textsc{Fortran } code. Note that
256the $k$-axis is orientated downward. The dashed area indicates the cell in
257which variables contained in arrays have the same $k$-index.}
258\end{center}   \end{figure}
259%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
260
261% -----------------------------------
262%        Domain Size
263% -----------------------------------
264\subsubsection{Domain Size}
265\label{DOM_size}
266
267The total size of the computational domain is set by the parameters \np{jpiglo},
268\np{jpjglo} and \np{jpkdta} in the $i$, $j$ and $k$ directions respectively. They are
269given as namelist variables in the \ngn{namcfg} namelist.
270
271Note that are other namelist variables in the \ngn{namcfg} namelist that refer to
272 the domain size.
273The two variables \np{jpidta} and \np{jpjdta} may be larger than \np{jpiglo}, \np{jpjglo}
274when the user wants to use only a sub-region of a given configuration. This is
275the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of
276the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters
277$jpi$ and $jpj$ refer to the size of each processor subdomain when the code is
278run in parallel using domain decomposition (\key{mpp\_mpi} defined, see
279\S\ref{LBC_mpp}).
280
281
282$\ $\newline    % force a new lign
283
284% ================================================================
285% Domain: Horizontal Grid (mesh)
286% ================================================================
287\section  [Domain: Horizontal Grid (mesh) (\textit{domhgr})]               
288      {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} }
289\label{DOM_hgr}
290
291% -------------------------------------------------------------------------------------------------------------
292%        Coordinates and scale factors
293% -------------------------------------------------------------------------------------------------------------
294\subsection{Coordinates and scale factors}
295\label{DOM_hgr_coord_e}
296
297The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined
298by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$.
299The grid-points are located at integer or integer and a half values of as indicated
300in Table~\ref{Tab_cell}. The associated scale factors are defined using the
301analytical first derivative of the transformation \eqref{Eq_scale_factors}. These
302definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which
303provide the horizontal and vertical meshes, respectively. This section deals with
304the horizontal mesh parameters.
305
306In a horizontal plane, the location of all the model grid points is defined from the
307analytical expressions of the longitude $\lambda$ and  latitude $\varphi$ as a
308function of  $(i,j)$. The horizontal scale factors are calculated using
309\eqref{Eq_scale_factors}. For example, when the longitude and latitude are
310function of a single value ($i$ and $j$, respectively) (geographical configuration
311of the mesh), the horizontal mesh definition reduces to define the wanted
312$\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the
313\mdl{domhgr} module. The model computes the grid-point positions and scale
314factors in the horizontal plane as follows:
315\begin{flalign*}
316\lambda_t &\equiv \text{glamt}= \lambda(i)     & \varphi_t &\equiv \text{gphit} = \varphi(j)\\
317\lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\
318\lambda_v &\equiv \text{glamv}= \lambda(i)       & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\
319\lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2)
320\end{flalign*}
321\begin{flalign*}
322e_{1t} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j)  |&
323e_{2t} &\equiv \text{e2t} = r_a |\varphi'(j)|  \\
324e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)   \; \cos\varphi(j)  |&
325e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\
326e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i)    \; \cos\varphi(j+1/2)  |&
327e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\
328e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2)  |&
329e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|
330\end{flalign*}
331where the last letter of each computational name indicates the grid point
332considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with
333all universal constants). Note that the horizontal position of and scale factors
334at $w$-points are exactly equal to those of $t$-points, thus no specific arrays
335are defined at $w$-points.
336
337Note that the definition of the scale factors ($i.e.$ as the analytical first derivative
338of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is
339specific to the \NEMO model \citep{Marti_al_JGR92}. As an example, $e_{1t}$ is defined
340locally at a $t$-point, whereas many other models on a C grid choose to define
341such a scale factor as the distance between the $U$-points on each side of the
342$t$-point. Relying on an analytical transformation has two advantages: firstly, there
343is no ambiguity in the scale factors appearing in the discrete equations, since they
344are first introduced in the continuous equations; secondly, analytical transformations
345encourage good practice by the definition of smoothly varying grids (rather than
346allowing the user to set arbitrary jumps in thickness between adjacent layers)
347\citep{Treguier1996}. An example of the effect of such a choice is shown in
348Fig.~\ref{Fig_zgr_e3}.
349%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
350\begin{figure}[!t]     \begin{center}
351\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf}
352\caption{ \label{Fig_zgr_e3}   
353Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical,
354and (b) analytically derived grid-point position and scale factors.
355For both grids here,  the same $w$-point depth has been chosen but in (a) the
356$t$-points are set half way between $w$-points while in (b) they are defined from
357an analytical function: $z(k)=5\,(k-1/2)^3 - 45\,(k-1/2)^2 + 140\,(k-1/2) - 150$.
358Note the resulting difference between the value of the grid-size $\Delta_k$ and
359those of the scale factor $e_k$. }
360\end{center}   \end{figure}
361%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
362
363% -------------------------------------------------------------------------------------------------------------
364%        Choice of horizontal grid
365% -------------------------------------------------------------------------------------------------------------
366\subsection{Choice of horizontal grid}
367\label{DOM_hgr_msh_choice}
368
369The user has three options available in defining a horizontal grid, which involve
370the namelist variable \np{jphgr\_mesh} of the \ngn{namcfg} namelist.
371\begin{description}
372\item[\np{jphgr\_mesh}=0]  The most general curvilinear orthogonal grids.
373The coordinates and their first derivatives with respect to $i$ and $j$ are provided
374in a input file (\ifile{coordinates}), read in \rou{hgr\_read} subroutine of the domhgr module.
375\item[\np{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below).
376For other analytical grids, the \mdl{domhgr} module must be modified by the user.
377\end{description}
378
379There are two simple cases of geographical grids on the sphere. With
380\np{jphgr\_mesh}=1, the grid (expressed in degrees) is regular in space,
381with grid sizes specified by parameters \np{ppe1\_deg} and \np{ppe2\_deg},
382respectively. Such a geographical grid can be very anisotropic at high latitudes
383because of the convergence of meridians (the zonal scale factors $e_1$ 
384become much smaller than the meridional scale factors $e_2$). The Mercator
385grid (\np{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale
386factors in the same way as the zonal ones. In this case, meridional scale factors
387and latitudes are calculated analytically using the formulae appropriate for
388a Mercator projection, based on \np{ppe1\_deg} which is a reference grid spacing
389at the equator (this applies even when the geographical equator is situated outside
390the model domain).
391%%%
392\gmcomment{ give here the analytical expression of the Mercator mesh}
393%%%
394In these two cases (\np{jphgr\_mesh}=1 or 4), the grid position is defined by the
395longitude and latitude of the south-westernmost point (\np{ppglamt0} 
396and \np{ppgphi0}). Note that for the Mercator grid the user need only provide
397an approximate starting latitude: the real latitude will be recalculated analytically,
398in order to ensure that the equator corresponds to line passing through $t$-
399and $u$-points. 
400
401Rectangular grids ignoring the spherical geometry are defined with
402\np{jphgr\_mesh} = 2, 3, 5. The domain is either an $f$-plane (\np{jphgr\_mesh} = 2,
403Coriolis factor is constant) or a beta-plane (\np{jphgr\_mesh} = 3, the Coriolis factor
404is linear in the $j$-direction). The grid size is uniform in meter in each direction,
405and given by the parameters \np{ppe1\_m} and \np{ppe2\_m} respectively.
406The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero
407with the first $t$-point. The meridional coordinate (gphi. arrays) is in kilometers,
408and the second $t$-point corresponds to coordinate $gphit=0$. The input
409variable \np{ppglam0} is ignored. \np{ppgphi0} is used to set the reference
410latitude for computation of the Coriolis parameter. In the case of the beta plane,
411\np{ppgphi0} corresponds to the center of the domain. Finally, the special case
412\np{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the
413GYRE configuration, representing a classical mid-latitude double gyre system.
414The rotation allows us to maximize the jet length relative to the gyre areas
415(and the number of grid points).
416
417The choice of the grid must be consistent with the boundary conditions specified
418by \np{jperio}, a parameter found in \ngn{namcfg} namelist (see {\S\ref{LBC}).
419
420% -------------------------------------------------------------------------------------------------------------
421%        Grid files
422% -------------------------------------------------------------------------------------------------------------
423\subsection{Output Grid files}
424\label{DOM_hgr_files}
425
426All the arrays relating to a particular ocean model configuration (grid-point
427position, scale factors, masks) can be saved in files if $\np{nn\_msh} \not= 0$ 
428(namelist variable in \ngn{namdom}). This can be particularly useful for plots and off-line
429diagnostics. In some cases, the user may choose to make a local modification
430of a scale factor in the code. This is the case in global configurations when
431restricting the width of a specific strait (usually a one-grid-point strait that
432happens to be too wide due to insufficient model resolution). An example
433is Gibraltar Strait in the ORCA2 configuration. When such modifications are done,
434the output grid written when $\np{nn\_msh} \not=0$ is no more equal to the input grid.
435
436$\ $\newline    % force a new lign
437
438% ================================================================
439% Domain: Vertical Grid (domzgr)
440% ================================================================
441\section  [Domain: Vertical Grid (\textit{domzgr})]
442      {Domain: Vertical Grid \small{(\mdl{domzgr} module)} }
443\label{DOM_zgr}
444%-----------------------------------------nam_zgr & namdom-------------------------------------------
445\namdisplay{namzgr} 
446\namdisplay{namdom} 
447%-------------------------------------------------------------------------------------------------------------
448
449Variables are defined through the \ngn{namzgr} and \ngn{namdom} namelists.
450In the vertical, the model mesh is determined by four things:
451(1) the bathymetry given in meters ;
452(2) the number of levels of the model (\jp{jpk}) ;
453(3) the analytical transformation $z(i,j,k)$ and the vertical scale factors
454(derivatives of the transformation) ;
455and (4) the masking system, $i.e.$ the number of wet model levels at each
456$(i,j)$ column of points.
457
458%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
459\begin{figure}[!tb]    \begin{center}
460\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf}
461\caption{  \label{Fig_z_zps_s_sps}   
462The ocean bottom as seen by the model:
463(a) $z$-coordinate with full step,
464(b) $z$-coordinate with partial step,
465(c) $s$-coordinate: terrain following representation,
466(d) hybrid $s-z$ coordinate,
467(e) hybrid $s-z$ coordinate with partial step, and
468(f) same as (e) but in the non-linear free surface (\np{ln\_linssh}=false).
469Note that the non-linear free surface can be used with any of the
4705 coordinates (a) to (e).}
471\end{center}   \end{figure}
472%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
473
474The choice of a vertical coordinate, even if it is made through \ngn{namzgr} namelist parameters,
475must be done once of all at the beginning of an experiment. It is not intended as an
476option which can be enabled or disabled in the middle of an experiment. Three main
477choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step
478bathymetry (\np{ln\_zco}~=~true), $z$-coordinate with partial step bathymetry
479(\np{ln\_zps}~=~true), or generalized, $s$-coordinate (\np{ln\_sco}~=~true).
480Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate
481(Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). By default a non-linear free surface is used:
482the coordinate follow the time-variation of the free surface so that the transformation is time dependent:
483$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). When a linear free surface is assumed (\np{ln\_linssh}=true),
484the vertical coordinate are fixed in time, but the seawater can move up and down across the z=0 surface
485(in other words, the top of the ocean in not a rigid-lid).
486The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain
487of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities,
488otherwise they are filled in. This option is currently only available in $z$- or $zps$-coordinate,
489and partial step are also applied at the ocean/ice shelf interface.
490
491Contrary to the horizontal grid, the vertical grid is computed in the code and no
492provision is made for reading it from a file. The only input file is the bathymetry
493(in meters) (\ifile{bathy\_meter}).
494\footnote{N.B. in full step $z$-coordinate, a \ifile{bathy\_level} file can replace the
495\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point
496in each water column is by-passed}.
497If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft
498(in meters) (\ifile{isf\_draft\_meter}) is needed.
499
500After reading the bathymetry, the algorithm for vertical grid definition differs
501between the different options:
502\begin{description}
503\item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$.
504\item[\textit{zps}] set a reference coordinate transformation $z_0 (k)$, and
505calculate the thickness of the deepest level at each $(i,j)$ point using the
506bathymetry, to obtain the final three-dimensional depth and scale factor arrays.
507\item[\textit{sco}] smooth the bathymetry to fulfil the hydrostatic consistency
508criteria and set the three-dimensional transformation.
509\item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfil the hydrostatic
510consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and
511possibly introduce masking of extra land points to better fit the original bathymetry file
512\end{description}
513%%%
514\gmcomment{   add the description of the smoothing:  envelop topography...}
515%%%
516
517Unless a linear free surface is used (\np{ln\_linssh}=false), the arrays describing
518the grid point depths and vertical scale factors are three set of three dimensional arrays $(i,j,k)$ 
519defined at \textit{before}, \textit{now} and \textit{after} time step. The time at which they are
520defined is indicated by a suffix:$\_b$, $\_n$, or $\_a$, respectively. They are updated at each model time step
521using a fixed reference coordinate system which computer names have a $\_0$ suffix.
522When the linear free surface option is used (\np{ln\_linssh}=true), \textit{before}, \textit{now} 
523and \textit{after} arrays are simply set one for all to their reference counterpart.
524
525
526% -------------------------------------------------------------------------------------------------------------
527%        Meter Bathymetry
528% -------------------------------------------------------------------------------------------------------------
529\subsection{Meter Bathymetry}
530\label{DOM_bathy}
531
532Three options are possible for defining the bathymetry, according to the
533namelist variable \np{nn\_bathy} (found in \ngn{namdom} namelist):
534\begin{description}
535\item[\np{nn\_bathy} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ 
536is given by the coordinate transformation. The domain can either be a closed
537basin or a periodic channel depending on the parameter \np{jperio}.
538\item[\np{nn\_bathy} = -1] a domain with a bump of topography one third of the
539domain width at the central latitude. This is meant for the "EEL-R5" configuration,
540a periodic or open boundary channel with a seamount.
541\item[\np{nn\_bathy} = 1] read a bathymetry and ice shelf draft (if needed).
542 The \ifile{bathy\_meter} file (Netcdf format) provides the ocean depth (positive, in meters)
543 at each grid point of the model grid. The bathymetry is usually built by interpolating a standard bathymetry product
544($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also
545defines the coastline: where the bathymetry is zero, no model levels are defined
546(all levels are masked).
547
548The \ifile{isfdraft\_meter} file (Netcdf format) provides the ice shelf draft (positive, in meters)
549 at each grid point of the model grid. This file is only needed if \np{ln\_isfcav}~=~true.
550Defining the ice shelf draft will also define the ice shelf edge and the grounding line position.
551\end{description}
552
553When a global ocean is coupled to an atmospheric model it is better to represent
554all large water bodies (e.g, great lakes, Caspian sea...) even if the model
555resolution does not allow their communication with the rest of the ocean.
556This is unnecessary when the ocean is forced by fixed atmospheric conditions,
557so these seas can be removed from the ocean domain. The user has the option
558to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the
559code has to be adapted to the user's configuration.
560
561% -------------------------------------------------------------------------------------------------------------
562%        z-coordinate  and reference coordinate transformation
563% -------------------------------------------------------------------------------------------------------------
564\subsection[$z$-coordinate (\np{ln\_zco}]
565        {$z$-coordinate (\np{ln\_zco}=true) and reference coordinate}
566\label{DOM_zco}
567
568%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
569\begin{figure}[!tb]    \begin{center}
570\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf}
571\caption{ \label{Fig_zgr}   
572Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level functions for
573(a) T-point depth and (b) the associated scale factor as computed
574from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.}
575\end{center}   \end{figure}
576%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
577
578The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ 
579and $gdepw_0$ for $t$- and $w$-points, respectively. As indicated on
580Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the
581ocean surface. There are at most \jp{jpk}-1 $t$-points inside the ocean, the
582additional $t$-point at $jk=jpk$ is below the sea floor and is not used.
583The vertical location of $w$- and $t$-levels is defined from the analytic expression
584of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the
585vertical scale factors. The user must provide the analytical expression of both
586$z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} 
587through statement functions, using parameters provided in the \ngn{namcfg} namelist.
588
589It is possible to define a simple regular vertical grid by giving zero stretching (\np{ppacr=0}).
590In that case, the parameters \jp{jpk} (number of $w$-levels) and \np{pphmax} 
591(total ocean depth in meters) fully define the grid.
592
593For climate-related studies it is often desirable to concentrate the vertical resolution
594near the ocean surface. The following function is proposed as a standard for a
595$z$-coordinate (with either full or partial steps):
596\begin{equation} \label{DOM_zgr_ana}
597\begin{split}
598 z_0 (k)    &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ 
599 e_3^0 (k)  &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right|
600\end{split}
601\end{equation}
602where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an
603expression allows us to define a nearly uniform vertical location of levels at the
604ocean top and bottom with a smooth hyperbolic tangent transition in between
605(Fig.~\ref{Fig_zgr}).
606
607If the ice shelf cavities are opened (\np{ln\_isfcav}=~true~}), the definition of $z_0$ is the same.
608However, definition of $e_3^0$ at $t$- and $w$-points is respectively changed to:
609\begin{equation} \label{DOM_zgr_ana}
610\begin{split}
611 e_3^T(k) &= z_W (k+1) - z_W (k)   \\
612 e_3^W(k) &= z_T (k)   - z_T (k-1) \\
613\end{split}
614\end{equation}
615This formulation decrease the self-generated circulation into the ice shelf cavity
616(which can, in extreme case, leads to blow up).\\
617
618 
619The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the
620surface (bottom) layers and a depth which varies from 0 at the sea surface to a
621minimum of $-5000~m$. This leads to the following conditions:
622\begin{equation} \label{DOM_zgr_coef}
623\begin{split}
624 e_3 (1+1/2)      &=10. \\ 
625 e_3 (jpk-1/2) &=500. \\ 
626 z(1)       &=0. \\ 
627 z(jpk)        &=-5000. \\ 
628\end{split}
629\end{equation}
630
631With the choice of the stretching $h_{cr} =3$ and the number of levels
632\jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in
633\eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is
634satisfied, through an optimisation procedure using a bisection method. For the first
635standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$,
636$h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and
637scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and
638given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters
639\np{ppsur}, \np{ppa0}, \np{ppa1}, \np{ppkth} in \ngn{namcfg} namelist.
640
641Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is
642possible to recalculate them. In that case the user sets
643\np{ppsur}=\np{ppa0}=\np{ppa1}=999999., in \ngn{namcfg} namelist,
644and specifies instead the four following parameters:
645\begin{itemize}
646\item    \np{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger
647\np{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual.
648\item    \np{ppkth}=$h_{th} $: is approximately the model level at which maximum
649stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk})
650\item    \np{ppdzmin}: minimum thickness for the top layer (in meters)
651\item    \np{pphmax}: total depth of the ocean (meters).
652\end{itemize}
653As an example, for the $45$ layers used in the DRAKKAR configuration those
654parameters are: \jp{jpk}=46, \np{ppacr}=9, \np{ppkth}=23.563, \np{ppdzmin}=6m,
655\np{pphmax}=5750m.
656
657%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
658\begin{table}     \begin{center} \begin{tabular}{c||r|r|r|r}
659\hline
660\textbf{LEVEL}& \textbf{gdept}& \textbf{gdepw}& \textbf{e3t }& \textbf{e3w  } \\ \hline
661&  \textbf{  5.00}   &       0.00 & \textbf{ 10.00} &  10.00 \\   \hline
662&  \textbf{15.00} &    10.00 &   \textbf{ 10.00} &  10.00 \\   \hline
663&  \textbf{25.00} &    20.00 &   \textbf{ 10.00} &     10.00 \\   \hline
664&  \textbf{35.01} &    30.00 &   \textbf{ 10.01} &     10.00 \\   \hline
665&  \textbf{45.01} &    40.01 &   \textbf{ 10.01} &  10.01 \\   \hline
666&  \textbf{55.03} &    50.02 &   \textbf{ 10.02} &     10.02 \\   \hline
667&  \textbf{65.06} &    60.04 &   \textbf{ 10.04} &  10.03 \\   \hline
668&  \textbf{75.13} &    70.09 &   \textbf{ 10.09} &  10.06 \\   \hline
669&  \textbf{85.25} &    80.18 &   \textbf{ 10.17} &  10.12 \\   \hline
67010 &  \textbf{95.49} &    90.35 &   \textbf{ 10.33} &  10.24 \\   \hline
67111 &  \textbf{105.97}   &   100.69 &   \textbf{ 10.65} &  10.47 \\   \hline
67212 &  \textbf{116.90}   &   111.36 &   \textbf{ 11.27} &  10.91 \\   \hline
67313 &  \textbf{128.70}   &   122.65 &   \textbf{ 12.47} &  11.77 \\   \hline
67414 &  \textbf{142.20}   &   135.16 &   \textbf{ 14.78} &  13.43 \\   \hline
67515 &  \textbf{158.96}   &   150.03 &   \textbf{ 19.23} &  16.65 \\   \hline
67616 &  \textbf{181.96}   &   169.42 &   \textbf{ 27.66} &  22.78 \\   \hline
67717 &  \textbf{216.65}   &   197.37 &   \textbf{ 43.26} &  34.30 \\ \hline
67818 &  \textbf{272.48}   &   241.13 &   \textbf{ 70.88} &  55.21 \\ \hline
67919 &  \textbf{364.30}   &   312.74 &   \textbf{116.11} &  90.99 \\ \hline
68020 &  \textbf{511.53}   &   429.72 &   \textbf{181.55} &    146.43 \\ \hline
68121 &  \textbf{732.20}   &   611.89 &   \textbf{261.03} &    220.35 \\ \hline
68222 &  \textbf{1033.22}&  872.87 &   \textbf{339.39} &    301.42 \\ \hline
68323 &  \textbf{1405.70}& 1211.59 & \textbf{402.26} &   373.31 \\ \hline
68424 &  \textbf{1830.89}& 1612.98 & \textbf{444.87} &   426.00 \\ \hline
68525 &  \textbf{2289.77}& 2057.13 & \textbf{470.55} &   459.47 \\ \hline
68626 &  \textbf{2768.24}& 2527.22 & \textbf{484.95} &   478.83 \\ \hline
68727 &  \textbf{3257.48}& 3011.90 & \textbf{492.70} &   489.44 \\ \hline
68828 &  \textbf{3752.44}& 3504.46 & \textbf{496.78} &   495.07 \\ \hline
68929 &  \textbf{4250.40}& 4001.16 & \textbf{498.90} &   498.02 \\ \hline
69030 &  \textbf{4749.91}& 4500.02 & \textbf{500.00} &   499.54 \\ \hline
69131 &  \textbf{5250.23}& 5000.00 &   \textbf{500.56} & 500.33 \\ \hline
692\end{tabular} \end{center} 
693\caption{ \label{Tab_orca_zgr}   
694Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration as computed
695from \eqref{DOM_zgr_ana} using the coefficients given in \eqref{DOM_zgr_coef}}
696\end{table}
697%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
698
699% -------------------------------------------------------------------------------------------------------------
700%        z-coordinate with partial step
701% -------------------------------------------------------------------------------------------------------------
702\subsection   [$z$-coordinate with partial step (\np{ln\_zps})]
703         {$z$-coordinate with partial step (\np{ln\_zps}=.true.)}
704\label{DOM_zps}
705%--------------------------------------------namdom-------------------------------------------------------
706\namdisplay{namdom} 
707%--------------------------------------------------------------------------------------------------------------
708
709In $z$-coordinate partial step, the depths of the model levels are defined by the
710reference analytical function $z_0 (k)$ as described in the previous
711section, \emph{except} in the bottom layer. The thickness of the bottom layer is
712allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a
713better representation of the bathymetry, especially in the case of small
714slopes (where the bathymetry varies by less than one level thickness from
715one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been
716defined in the absence of bathymetry. With partial steps, layers from 1 to
717\jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1)
718is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the
719maximum thickness allowed is $2*e_{3t}(jpk-1)$. This has to be kept in mind when
720specifying values in \ngn{namdom} namelist, as the maximum depth \np{pphmax} 
721in partial steps: for example, with
722\np{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth
723allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$).
724Two variables in the namdom namelist are used to define the partial step
725vertical grid. The mimimum water thickness (in meters) allowed for a cell
726partially filled with bathymetry at level jk is the minimum of \np{rn\_e3zps\_min} 
727(thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{rn\_e3zps\_rat}$ (a fraction,
728usually 10\%, of the default thickness $e_{3t}(jk)$).
729
730\gmcomment{ \colorbox{yellow}{Add a figure here of pstep especially at last ocean level }  }
731
732% -------------------------------------------------------------------------------------------------------------
733%        s-coordinate
734% -------------------------------------------------------------------------------------------------------------
735\subsection   [$s$-coordinate (\np{ln\_sco})]
736           {$s$-coordinate (\np{ln\_sco}=true)}
737\label{DOM_sco}
738%------------------------------------------nam_zgr_sco---------------------------------------------------
739\namdisplay{namzgr_sco} 
740%--------------------------------------------------------------------------------------------------------------
741Options are defined in \ngn{namzgr\_sco}.
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
759the sea bed and becomes a pseudo z-coordinate.
760The coordinate can also be hybridised by specifying \np{rn\_sbot\_min} and \np{rn\_sbot\_max} 
761as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
762
763Options for stretching the coordinate are provided as examples, but care must be taken to ensure
764that the vertical stretch used is appropriate for the application.
765
766The original default NEMO s-coordinate stretching is available if neither of the other options
767are specified as true (\np{ln\_s\_SH94}~=~false and \np{ln\_s\_SF12}~=~false).
768This uses a depth independent $\tanh$ function for the stretching \citep{Madec_al_JPO96}:
769
770\begin{equation}
771  z = s_{min}+C\left(s\right)\left(H-s_{min}\right)
772  \label{eq:SH94_1}
773\end{equation}
774
775where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and
776allows a $z$-coordinate to placed on top of the stretched coordinate,
777and $z$ is the depth (negative down from the asea surface).
778
779\begin{equation}
780  s = -\frac{k}{n-1} \quad \text{ and } \quad 0 \leq k \leq n-1
781  \label{eq:s}
782\end{equation}
783
784\begin{equation} \label{DOM_sco_function}
785\begin{split}
786C(s)  &\frac{ \left[   \tanh{ \left( \theta \, (s+b) \right)} 
787               - \tanh{ \left\theta \, b      \right)}  \right]}
788            {2\;\sinh \left( \theta \right)}
789\end{split}
790\end{equation}
791
792A stretching function, modified from the commonly used \citet{Song_Haidvogel_JCP94} 
793stretching (\np{ln\_s\_SH94}~=~true), is also available and is more commonly used for shelf seas modelling:
794
795\begin{equation}
796  C\left(s\right) =   \left(1 - b \right)\frac{ \sinh\left( \theta s\right)}{\sinh\left(\theta\right)} +      \\
797  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)}
798  \label{eq:SH94_2}
799\end{equation}
800
801%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
802\begin{figure}[!ht]    \begin{center}
803\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf}
804\caption{  \label{Fig_sco_function}   
805Examples of the stretching function applied to a seamount; from left to right:
806surface, surface and bottom, and bottom intensified resolutions}
807\end{center}   \end{figure}
808%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
809
810where $H_c$ is the critical depth (\np{rn\_hc}) at which the coordinate transitions from
811pure $\sigma$ to the stretched coordinate,  and $\theta$ (\np{rn\_theta}) and $b$ (\np{rn\_bb})
812are the surface and bottom control parameters such that $0\leqslant \theta \leqslant 20$, and
813$0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom
814increase of the vertical resolution (Fig.~\ref{Fig_sco_function}).
815
816Another example has been provided at version 3.5 (\np{ln\_s\_SF12}) that allows
817a fixed surface resolution in an analytical terrain-following stretching \citet{Siddorn_Furner_OM12}.
818In this case the a stretching function $\gamma$ is defined such that:
819
820\begin{equation}
821z = -\gamma h \quad \text{ with } \quad 0 \leq \gamma \leq 1
822\label{eq:z}
823\end{equation}
824
825The function is defined with respect to $\sigma$, the unstretched terrain-following coordinate:
826
827\begin{equation} \label{DOM_gamma_deriv}
828\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)
829\end{equation}
830
831Where:
832\begin{equation} \label{DOM_gamma}
833f\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} 
834\end{equation}
835
836This gives an analytical stretching of $\sigma$ that is solvable in $A$ and $B$ as a function of
837the user prescribed stretching parameter $\alpha$ (\np{rn\_alpha}) that stretches towards
838the surface ($\alpha > 1.0$) or the bottom ($\alpha < 1.0$) and user prescribed surface (\np{rn\_zs})
839and bottom depths. The bottom cell depth in this example is given as a function of water depth:
840
841\begin{equation} \label{DOM_zb}
842Z_b= h a + b
843\end{equation}
844
845where the namelist parameters \np{rn\_zb\_a} and \np{rn\_zb\_b} are $a$ and $b$ respectively.
846
847%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
848\begin{figure}[!ht]
849   \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/FIG_DOM_compare_coordinates_surface.pdf}
850        \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.}
851    \label{fig_compare_coordinates_surface}
852\end{figure}
853%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
854
855This 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.
856
857As 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$.
858
859Minimising 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.
860
861% -------------------------------------------------------------------------------------------------------------
862%        z*- or s*-coordinate
863% -------------------------------------------------------------------------------------------------------------
864\subsection{$z^*$- or $s^*$-coordinate (\np{ln\_linssh}=false) }
865\label{DOM_zgr_star}
866
867This option is described in the Report by Levier \textit{et al.} (2007), available on the \NEMO web site.
868
869%gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances
870
871% -------------------------------------------------------------------------------------------------------------
872%        level bathymetry and mask
873% -------------------------------------------------------------------------------------------------------------
874\subsection{level bathymetry and mask}
875\label{DOM_msk}
876
877Whatever the vertical coordinate used, the model offers the possibility of
878representing the bottom topography with steps that follow the face of the
879model cells (step like topography) \citep{Madec_al_JPO96}. The distribution of
880the steps in the horizontal is defined in a 2D integer array, mbathy, which
881gives the number of ocean levels ($i.e.$ those that are not masked) at each
882$t$-point. mbathy is computed from the meter bathymetry using the definiton of
883gdept as the number of $t$-points which gdept $\leq$ bathy.
884
885Modifications of the model bathymetry are performed in the \textit{bat\_ctl} 
886routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points
887that do not communicate with another ocean point at the same level are eliminated.
888
889As for the representation of bathymetry, a 2D integer array, misfdep, is created.
890misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked.
891By default, misfdep(:,:)=1 and no cells are masked.
892
893In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into
894the cavities are performed in the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked.
895All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked).
896If only one cell on the water column is opened at $t$-, $u$- or $v$-points, the bathymetry or the ice shelf draft is dug to fit this constrain.
897If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.\\ 
898
899From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows:
900\begin{align*}
901tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\
902                                \; 1&   \text{ if $misfdep(i,j) \leq k\leq mbathy(i,j)$  }    \\
903                                \; 0&   \text{ if $k > mbathy(i,j)$  }    \end{cases}     \\
904umask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
905vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\
906fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
907             &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
908wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1)
909\end{align*}
910
911Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with
912the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required
913with oceean cavities to deal with the top boundary (ice shelf/ocean interface)
914exactly in the same way as for the bottom boundary.
915
916The specification of closed lateral boundaries requires that at least the first and last
917rows and columns of the \textit{mbathy} array are set to zero. In the particular
918case of an east-west cyclical boundary condition, \textit{mbathy} has its last
919column equal to the second one and its first column equal to the last but one
920(and so too the mask arrays) (see \S~\ref{LBC_jperio}).
921
922
923% ================================================================
924% Domain: Initial State (dtatsd & istate)
925% ================================================================
926\section  [Domain: Initial State (\textit{istate and dtatsd})]
927      {Domain: Initial State \small{(\mdl{istate} and \mdl{dtatsd} modules)} }
928\label{DTA_tsd}
929%-----------------------------------------namtsd-------------------------------------------
930\namdisplay{namtsd} 
931%------------------------------------------------------------------------------------------
932
933Options are defined in \ngn{namtsd}.
934By default, the ocean start from rest (the velocity field is set to zero) and the initialization of
935temperature and salinity fields is controlled through the \np{ln\_tsd\_ini} namelist parameter.
936\begin{description}
937\item[ln\_tsd\_init = .true.]  use a T and S input files that can be given on the model grid itself or
938on their native input data grid. In the latter case, the data will be interpolated on-the-fly both in the
939horizontal and the vertical to the model grid (see \S~\ref{SBC_iof}). The information relative to the
940input files are given in the \np{sn\_tem} and \np{sn\_sal} structures.
941The computation is done in the \mdl{dtatsd} module.
942\item[ln\_tsd\_init = .false.] use constant salinity value of 35.5 psu and an analytical profile of temperature
943(typical of the tropical ocean), see \rou{istate\_t\_s} subroutine called from \mdl{istate} module.
944\end{description}
Note: See TracBrowser for help on using the repository browser.