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