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

Last change on this file since 7705 was 7705, checked in by flavoni, 4 years ago

update NEMO documentation in trunk

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