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