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