source: branches/UKMO/v3_6_extra_CMIP6_diagnostics/DOC/TexFiles/Chapters/Chap_DOM.tex @ 7492

Last change on this file since 7492 was 7492, checked in by timgraham, 5 years ago

Final merge of head of v3.6 STABLE branch into this branch

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