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 | \vfill |
---|
21 | \begin{figure}[b] |
---|
22 | \subsubsection*{Changes record} |
---|
23 | \begin{tabular}{m{0.08\linewidth}||m{0.32\linewidth}|m{0.6\linewidth}} |
---|
24 | Release & Author(s) & Modifications \\ |
---|
25 | \hline |
---|
26 | {\em 4.0} & {\em Simon M{\"u}ller \& Andrew Coward} & {\em Compatibility changes for v4.0. Major simplication has moved many of the options to external domain configuration tools. For now this information has been retained in \autoref{apdx:DOMAINcfg} } \\ |
---|
27 | {\em 3.x} & {\em Sebastien Masson, Gurvan Madec \& Rashid Benshila } & {\em } \\ |
---|
28 | \end{tabular} |
---|
29 | \end{figure} |
---|
30 | |
---|
31 | \newpage |
---|
32 | |
---|
33 | Having defined the continuous equations in \autoref{chap:PE} and chosen a time discretization \autoref{chap:STP}, |
---|
34 | we need to choose a grid for spatial discretization and related numerical algorithms. |
---|
35 | In the present chapter, we provide a general description of the staggered grid used in \NEMO, |
---|
36 | and other relevant information about the DOM (DOMain) source-code modules . |
---|
37 | |
---|
38 | % ================================================================ |
---|
39 | % Fundamentals of the Discretisation |
---|
40 | % ================================================================ |
---|
41 | \section{Fundamentals of the discretisation} |
---|
42 | \label{sec:DOM_basics} |
---|
43 | |
---|
44 | % ------------------------------------------------------------------------------------------------------------- |
---|
45 | % Arrangement of Variables |
---|
46 | % ------------------------------------------------------------------------------------------------------------- |
---|
47 | \subsection{Arrangement of variables} |
---|
48 | \label{subsec:DOM_cell} |
---|
49 | |
---|
50 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
51 | \begin{figure}[!tb] |
---|
52 | \begin{center} |
---|
53 | \includegraphics[width=\textwidth]{Fig_cell} |
---|
54 | \caption{ |
---|
55 | \protect\label{fig:cell} |
---|
56 | Arrangement of variables. |
---|
57 | $t$ indicates scalar points where temperature, salinity, density, pressure and |
---|
58 | horizontal divergence are defined. |
---|
59 | $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where both relative and |
---|
60 | planetary vorticities are defined. |
---|
61 | } |
---|
62 | \end{center} |
---|
63 | \end{figure} |
---|
64 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
65 | |
---|
66 | The numerical techniques used to solve the Primitive Equations in this model are based on the traditional, |
---|
67 | centred second-order finite difference approximation. |
---|
68 | Special attention has been given to the homogeneity of the solution in the three spatial directions. |
---|
69 | The arrangement of variables is the same in all directions. |
---|
70 | It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in |
---|
71 | the centre of each face of the cells (\autoref{fig:cell}). |
---|
72 | This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification |
---|
73 | \citep{mesinger.arakawa_bk76}. |
---|
74 | The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and |
---|
75 | the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. |
---|
76 | |
---|
77 | The ocean mesh (\ie the position of all the scalar and vector points) is defined by the transformation that |
---|
78 | gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. |
---|
79 | The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:cell}. |
---|
80 | In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of |
---|
81 | the grid-point where the scale factors are defined. |
---|
82 | Each scale factor is defined as the local analytical value provided by \autoref{eq:scale_factors}. |
---|
83 | As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and |
---|
84 | $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. |
---|
85 | Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation |
---|
86 | while the scale factors are chosen equal to their local analytical value. |
---|
87 | An important point here is that the partial derivative of the scale factors must be evaluated by |
---|
88 | centred finite difference approximation, not from their analytical expression. |
---|
89 | This preserves the symmetry of the discrete set of equations and therefore satisfies many of |
---|
90 | the continuous properties (see \autoref{apdx:C}). |
---|
91 | A similar, related remark can be made about the domain size: |
---|
92 | when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors |
---|
93 | (see \autoref{eq:DOM_bar} in the next section). |
---|
94 | |
---|
95 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
96 | \begin{table}[!tb] |
---|
97 | \begin{center} |
---|
98 | \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} |
---|
99 | \hline |
---|
100 | t & $i $ & $j $ & $k $ \\ |
---|
101 | \hline |
---|
102 | u & $i + 1/2$ & $j $ & $k $ \\ |
---|
103 | \hline |
---|
104 | v & $i $ & $j + 1/2$ & $k $ \\ |
---|
105 | \hline |
---|
106 | w & $i $ & $j $ & $k + 1/2$ \\ |
---|
107 | \hline |
---|
108 | f & $i + 1/2$ & $j + 1/2$ & $k $ \\ |
---|
109 | \hline |
---|
110 | uw & $i + 1/2$ & $j $ & $k + 1/2$ \\ |
---|
111 | \hline |
---|
112 | vw & $i $ & $j + 1/2$ & $k + 1/2$ \\ |
---|
113 | \hline |
---|
114 | fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ |
---|
115 | \hline |
---|
116 | \end{tabular} |
---|
117 | \caption{ |
---|
118 | \protect\label{tab:cell} |
---|
119 | Location of grid-points as a function of integer or integer and a half value of the column, line or level. |
---|
120 | This indexing is only used for the writing of the semi -discrete equations. |
---|
121 | In the code, the indexing uses integer values only and is positive downwards in the vertical with $k=1$ at the surface. |
---|
122 | (see \autoref{subsec:DOM_Num_Index}) |
---|
123 | } |
---|
124 | \end{center} |
---|
125 | \end{table} |
---|
126 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
127 | |
---|
128 | Note that the definition of the scale factors |
---|
129 | (\ie as the analytical first derivative of the transformation that |
---|
130 | results in $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) |
---|
131 | is specific to the \NEMO model \citep{marti.madec.ea_JGR92}. |
---|
132 | As an example, a scale factor in the $i$ direction is defined locally at a $t$-point, |
---|
133 | whereas many other models on a C grid choose to define such a scale factor as |
---|
134 | the distance between the $u$-points on each side of the $t$-point. |
---|
135 | Relying on an analytical transformation has two advantages: |
---|
136 | firstly, there is no ambiguity in the scale factors appearing in the discrete equations, |
---|
137 | since they are first introduced in the continuous equations; |
---|
138 | secondly, analytical transformations encourage good practice by the definition of smoothly varying grids |
---|
139 | (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. |
---|
140 | An example of the effect of such a choice is shown in \autoref{fig:zgr_e3}. |
---|
141 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
142 | \begin{figure}[!t] |
---|
143 | \begin{center} |
---|
144 | \includegraphics[width=\textwidth]{Fig_zgr_e3} |
---|
145 | \caption{ |
---|
146 | \protect\label{fig:zgr_e3} |
---|
147 | Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, |
---|
148 | and (b) analytically derived grid-point position and scale factors. |
---|
149 | For both grids here, the same $w$-point depth has been chosen but |
---|
150 | in (a) the $t$-points are set half way between $w$-points while |
---|
151 | in (b) they are defined from an analytical function: |
---|
152 | $z(k) = 5 \, (k - 1/2)^3 - 45 \, (k - 1/2)^2 + 140 \, (k - 1/2) - 150$. |
---|
153 | Note the resulting difference between the value of the grid-size $\Delta_k$ and |
---|
154 | those of the scale factor $e_k$. |
---|
155 | } |
---|
156 | \end{center} |
---|
157 | \end{figure} |
---|
158 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
159 | |
---|
160 | % ------------------------------------------------------------------------------------------------------------- |
---|
161 | % Vector Invariant Formulation |
---|
162 | % ------------------------------------------------------------------------------------------------------------- |
---|
163 | \subsection{Discrete operators} |
---|
164 | \label{subsec:DOM_operators} |
---|
165 | |
---|
166 | Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at |
---|
167 | the midpoint between them are: |
---|
168 | \begin{alignat*}{2} |
---|
169 | % \label{eq:di_mi} |
---|
170 | \delta_i [q] &= & &q (i + 1/2) - q (i - 1/2) \\ |
---|
171 | \overline q^{\, i} &= &\big\{ &q (i + 1/2) + q (i - 1/2) \big\} / 2 |
---|
172 | \end{alignat*} |
---|
173 | |
---|
174 | Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. |
---|
175 | Following \autoref{eq:PE_grad} and \autoref{eq:PE_lap}, the gradient of a variable $q$ defined at |
---|
176 | a $t$-point has its three components defined at $u$-, $v$- and $w$-points while |
---|
177 | its Laplacian is defined at the $t$-point. |
---|
178 | These operators have the following discrete forms in the curvilinear $s$-coordinates system: |
---|
179 | \[ |
---|
180 | % \label{eq:DOM_grad} |
---|
181 | \nabla q \equiv \frac{1}{e_{1u}} \delta_{i + 1/2} [q] \; \, \vect i |
---|
182 | + \frac{1}{e_{2v}} \delta_{j + 1/2} [q] \; \, \vect j |
---|
183 | + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k |
---|
184 | \] |
---|
185 | \begin{multline*} |
---|
186 | % \label{eq:DOM_lap} |
---|
187 | \Delta q \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} |
---|
188 | \; \lt[ \delta_i \lt( \frac{e_{2u} \, e_{3u}}{e_{1u}} \; \delta_{i + 1/2} [q] \rt) |
---|
189 | + \delta_j \lt( \frac{e_{1v} \, e_{3v}}{e_{2v}} \; \delta_{j + 1/2} [q] \rt) \; \rt] \\ |
---|
190 | + \frac{1}{e_{3t}} |
---|
191 | \delta_k \lt[ \frac{1 }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] |
---|
192 | \end{multline*} |
---|
193 | |
---|
194 | Following \autoref{eq:PE_curl} and \autoref{eq:PE_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at |
---|
195 | vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and |
---|
196 | its divergence defined at $t$-points: |
---|
197 | \begin{multline} |
---|
198 | % \label{eq:DOM_curl} |
---|
199 | \nabla \times \vect A \equiv \frac{1}{e_{2v} \, e_{3vw}} |
---|
200 | \Big[ \delta_{j + 1/2} (e_{3w} \, a_3) |
---|
201 | - \delta_{k + 1/2} (e_{2v} \, a_2) \Big] \vect i \\ |
---|
202 | + \frac{1}{e_{2u} \, e_{3uw}} |
---|
203 | \Big[ \delta_{k + 1/2} (e_{1u} \, a_1) |
---|
204 | - \delta_{i + 1/2} (e_{3w} \, a_3) \Big] \vect j \\ |
---|
205 | + \frac{1}{e_{1f} \, e_{2f}} |
---|
206 | \Big[ \delta_{i + 1/2} (e_{2v} \, a_2) |
---|
207 | - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k |
---|
208 | \end{multline} |
---|
209 | \begin{equation} |
---|
210 | % \label{eq:DOM_div} |
---|
211 | \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} |
---|
212 | \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] |
---|
213 | + \frac{1}{e_{3t}} \delta_k (a_3) |
---|
214 | \end{equation} |
---|
215 | |
---|
216 | The vertical average over the whole water column is denoted by an overbar and is for |
---|
217 | a masked field $q$ (\ie a quantity that is equal to zero inside solid areas): |
---|
218 | \begin{equation} |
---|
219 | \label{eq:DOM_bar} |
---|
220 | \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} |
---|
221 | \end{equation} |
---|
222 | where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, |
---|
223 | $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over |
---|
224 | all grid points of the same type in the direction indicated by the subscript (here $k$). |
---|
225 | |
---|
226 | In continuous form, the following properties are satisfied: |
---|
227 | \begin{gather} |
---|
228 | \label{eq:DOM_curl_grad} |
---|
229 | \nabla \times \nabla q = \vect 0 \\ |
---|
230 | \label{eq:DOM_div_curl} |
---|
231 | \nabla \cdot (\nabla \times \vect A) = 0 |
---|
232 | \end{gather} |
---|
233 | |
---|
234 | It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as |
---|
235 | the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at |
---|
236 | vector points $(u,v,w)$. |
---|
237 | |
---|
238 | Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. |
---|
239 | It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) |
---|
240 | are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, |
---|
241 | $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie |
---|
242 | \begin{alignat}{4} |
---|
243 | \label{eq:DOM_di_adj} |
---|
244 | &\sum \limits_i a_i \; \delta_i [b] &\equiv &- &&\sum \limits_i \delta _{ i + 1/2} [a] &b_{i + 1/2} \\ |
---|
245 | \label{eq:DOM_mi_adj} |
---|
246 | &\sum \limits_i a_i \; \overline b^{\, i} &\equiv & &&\sum \limits_i \overline a ^{\, i + 1/2} &b_{i + 1/2} |
---|
247 | \end{alignat} |
---|
248 | |
---|
249 | In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and |
---|
250 | $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. |
---|
251 | These two properties will be used extensively in the \autoref{apdx:C} to |
---|
252 | demonstrate integral conservative properties of the discrete formulation chosen. |
---|
253 | |
---|
254 | % ------------------------------------------------------------------------------------------------------------- |
---|
255 | % Numerical Indexing |
---|
256 | % ------------------------------------------------------------------------------------------------------------- |
---|
257 | \subsection{Numerical indexing} |
---|
258 | \label{subsec:DOM_Num_Index} |
---|
259 | |
---|
260 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
261 | \begin{figure}[!tb] |
---|
262 | \begin{center} |
---|
263 | \includegraphics[width=\textwidth]{Fig_index_hor} |
---|
264 | \caption{ |
---|
265 | \protect\label{fig:index_hor} |
---|
266 | Horizontal integer indexing used in the \fortran code. |
---|
267 | The dashed area indicates the cell in which variables contained in arrays have the same $i$- and $j$-indices |
---|
268 | } |
---|
269 | \end{center} |
---|
270 | \end{figure} |
---|
271 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
272 | |
---|
273 | The array representation used in the \fortran code requires an integer indexing. |
---|
274 | However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of |
---|
275 | integer values for $t$-points only while all the other points involve integer and a half values. |
---|
276 | Therefore, a specific integer indexing has been defined for points other than $t$-points |
---|
277 | (\ie velocity and vorticity grid-points). |
---|
278 | Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. |
---|
279 | |
---|
280 | % ----------------------------------- |
---|
281 | % Horizontal Indexing |
---|
282 | % ----------------------------------- |
---|
283 | \subsubsection{Horizontal indexing} |
---|
284 | \label{subsec:DOM_Num_Index_hor} |
---|
285 | |
---|
286 | The indexing in the horizontal plane has been chosen as shown in \autoref{fig:index_hor}. |
---|
287 | For an increasing $i$ index ($j$ index), |
---|
288 | the $t$-point and the eastward $u$-point (northward $v$-point) have the same index |
---|
289 | (see the dashed area in \autoref{fig:index_hor}). |
---|
290 | A $t$-point and its nearest northeast $f$-point have the same $i$-and $j$-indices. |
---|
291 | |
---|
292 | % ----------------------------------- |
---|
293 | % Vertical indexing |
---|
294 | % ----------------------------------- |
---|
295 | \subsubsection{Vertical indexing} |
---|
296 | \label{subsec:DOM_Num_Index_vertical} |
---|
297 | |
---|
298 | In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in |
---|
299 | the \fortran code is the reverse of that used in the semi -discrete equations and |
---|
300 | given in \autoref{subsec:DOM_cell}. |
---|
301 | The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below |
---|
302 | (\autoref{fig:index_vert}). |
---|
303 | The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while |
---|
304 | the last $t$-level is always outside the ocean domain (\autoref{fig:index_vert}). |
---|
305 | Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index (\ie $t$-points and their |
---|
306 | nearest $w$-point neighbour in negative index direction), in contrast to the indexing on the horizontal plane where |
---|
307 | the $t$-point has the same index as the nearest velocity points in the positive direction of the respective horizontal axis index |
---|
308 | (compare the dashed area in \autoref{fig:index_hor} and \autoref{fig:index_vert}). |
---|
309 | Since the scale factors are chosen to be strictly positive, |
---|
310 | a \textit{minus sign} is included in the \fortran implementations of \textit{all the vertical derivatives} of |
---|
311 | the discrete equations given in this manual in order to accommodate the opposing vertical index directions in implementation and documentation. |
---|
312 | |
---|
313 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
314 | \begin{figure}[!pt] |
---|
315 | \begin{center} |
---|
316 | \includegraphics[width=\textwidth]{Fig_index_vert} |
---|
317 | \caption{ |
---|
318 | \protect\label{fig:index_vert} |
---|
319 | Vertical integer indexing used in the \fortran code. |
---|
320 | Note that the $k$-axis is oriented downward. |
---|
321 | The dashed area indicates the cell in which variables contained in arrays have a common $k$-index. |
---|
322 | } |
---|
323 | \end{center} |
---|
324 | \end{figure} |
---|
325 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
326 | |
---|
327 | % ------------------------------------------------------------------------------------------------------------- |
---|
328 | % Domain configuration |
---|
329 | % ------------------------------------------------------------------------------------------------------------- |
---|
330 | \section{Spatial domain configuration} |
---|
331 | \label{subsec:DOM_config} |
---|
332 | |
---|
333 | \nlst{namcfg} |
---|
334 | |
---|
335 | Two typical methods are available to specify the spatial domain |
---|
336 | configuration; they can be selected using parameter \np{ln\_read\_cfg} |
---|
337 | parameter in namelist \ngn{namcfg}. |
---|
338 | |
---|
339 | If \np{ln\_read\_cfg} is set to \forcode{.true.}, the domain-specific parameters |
---|
340 | and fields are read from a netCDF input file, whose name (without its .nc |
---|
341 | suffix) can be specified as the value of the \np{cn\_domcfg} parameter in |
---|
342 | namelist \ngn{namcfg}. |
---|
343 | |
---|
344 | If \np{ln\_read\_cfg} is set to \forcode{.false.}, the domain-specific |
---|
345 | parameters and fields can be provided (\eg analytically computed) by subroutines |
---|
346 | \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. These subroutines can be supplied in |
---|
347 | the \path{MY_SRC} directory of the configuration, and default versions that |
---|
348 | configure the spatial domain for the GYRE reference configuration are present in |
---|
349 | the \path{src/OCE/USR} directory. |
---|
350 | |
---|
351 | In version 4.0 there are no longer any options for reading complex bathmetries and |
---|
352 | performing a vertical discretization at run-time. Whilst it is occasionally convenient |
---|
353 | to have a common bathymetry file and, for example, to run similar models with and |
---|
354 | without partial bottom boxes and/or sigma-coordinates, supporting such choices leads to |
---|
355 | overly complex code. Worse still is the difficulty of ensuring the model configurations |
---|
356 | intended to be identical are indeed so when the model domain itself can be altered by runtime |
---|
357 | selections. The code previously used to perform vertical discretization has be incorporated |
---|
358 | into an external tool (\path{tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMAINcfg}. |
---|
359 | |
---|
360 | The next subsections summarise the parameter and fields related to the |
---|
361 | configuration of the whole model domain. These represent the minimum information |
---|
362 | that must be provided either via the \np{cn\_domcfg} file or set by code |
---|
363 | inserted into user-supplied versions of the \mdl{usrdef\_*} subroutines. The |
---|
364 | requirements are presented in three sections: the domain size |
---|
365 | (\autoref{subsec:DOM_size}), the horizontal mesh |
---|
366 | (\autoref{subsec:DOM_hgr}), and the vertical grid |
---|
367 | (\autoref{subsec:DOM_zgr}). |
---|
368 | |
---|
369 | % ----------------------------------- |
---|
370 | % Domain Size |
---|
371 | % ----------------------------------- |
---|
372 | \subsection{Domain size} |
---|
373 | \label{subsec:DOM_size} |
---|
374 | |
---|
375 | The total size of the computational domain is set by the parameters |
---|
376 | \np{jpiglo}, \np{jpjglo} and \np{jpkglo} for the $i$, $j$ and $k$ |
---|
377 | directions, respectively. Note, that the variables \forcode{jpi} and \forcode{jpj} |
---|
378 | refer to the size of each processor subdomain when the code is run in |
---|
379 | parallel using domain decomposition (\key{mpp\_mpi} defined, see |
---|
380 | \autoref{sec:LBC_mpp}). |
---|
381 | |
---|
382 | The name of the configuration is set through parameter \np{cn\_cfg}, |
---|
383 | and the nominal resolution through parameter \np{nn\_cfg} (unless in |
---|
384 | the input file both of variables \forcode{ORCA} and \forcode{ORCA_index} |
---|
385 | are present, in which case \np{cn\_cfg} and \np{nn\_cfg} are set from these |
---|
386 | values accordingly). |
---|
387 | |
---|
388 | The global lateral boundary condition type is selected from 8 options |
---|
389 | using parameter \np{jperio}. See \autoref{sec:LBC_jperio} for |
---|
390 | details on the available options and the corresponding values for |
---|
391 | \np{jperio}. |
---|
392 | |
---|
393 | % ================================================================ |
---|
394 | % Domain: Horizontal Grid (mesh) |
---|
395 | % ================================================================ |
---|
396 | \subsection{Horizontal grid mesh (\protect\mdl{domhgr})} |
---|
397 | \label{subsec:DOM_hgr} |
---|
398 | |
---|
399 | % ================================================================ |
---|
400 | % Domain: List of hgr-related fields needed |
---|
401 | % ================================================================ |
---|
402 | \subsubsection{Required fields} |
---|
403 | \label{sec:DOM_hgr_fields} |
---|
404 | The explicit specification of a range of mesh-related fields are required for the definition of a configuration. These include: |
---|
405 | |
---|
406 | \begin{Verbatim}[fontsize=\tiny] |
---|
407 | int jpiglo, jpjglo, jpkglo /* global domain sizes */ |
---|
408 | int jperio /* lateral global domain b.c. */ |
---|
409 | double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ |
---|
410 | double gphit, gphiu, gphiv, gphif /* geographic latitude */ |
---|
411 | double e1t, e1u, e1v, e1f /* horizontal scale factors */ |
---|
412 | double e2t, e2u, e2v, e2f /* horizontal scale factors */ |
---|
413 | \end{Verbatim} |
---|
414 | |
---|
415 | The values of the geographic longitude and latitude arrays at indices $i,j$ correspond to the analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a function of $(i,j)$, evaluated at the values as specified in Table \autoref{tab:cell} for the respective grid-point position. The calculation of the values of the horizontal scale factor arrays in general additionally involves partial derivatives of $\lambda$ and $\varphi$ with respect to $i$ and $j$, evaluated for the same arguments as $\lambda$ and $\varphi$. |
---|
416 | |
---|
417 | \subsubsection{Optional fields} |
---|
418 | \begin{Verbatim}[fontsize=\tiny] |
---|
419 | /* Optional: */ |
---|
420 | int ORCA, ORCA_index /* configuration name, configuration resolution */ |
---|
421 | double e1e2u, e1e2v /* U and V surfaces (if grid size reduction in some straits) */ |
---|
422 | double ff_f, ff_t /* Coriolis parameter (if not on the sphere) */ |
---|
423 | \end{Verbatim} |
---|
424 | |
---|
425 | NEMO can support the local reduction of key strait widths by altering individual values of |
---|
426 | e2u or e1v at the appropriate locations. This is particularly useful for locations such as |
---|
427 | Gibraltar or Indonesian Throughflow pinch-points (see \autoref{sec:MISC_strait} for |
---|
428 | illustrated examples). The key is to reduce the faces of $T$-cell (\ie change the value of |
---|
429 | the horizontal scale factors at $u$- or $v$-point) but not the volume of the cells. Doing |
---|
430 | otherwise can lead to numerical instability issues. In normal operation the surface areas |
---|
431 | are computed from $\texttt{e1u} * \texttt{e2u}$ and $\texttt{e1v} * \texttt{e2v}$ but in |
---|
432 | cases where a gridsize reduction is required, the unaltered surface areas at $u$ and $v$ |
---|
433 | grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed |
---|
434 | in \mdl{usrdef\_hgr}. If these arrays are present in the \np{cn\_domcfg} file they are |
---|
435 | read and the internal computation is suppressed. Versions of \mdl{usrdef\_hgr} which set |
---|
436 | their own values of \texttt{e1e2u} and \texttt{e1e2v} should set the surface-area |
---|
437 | computation flag: \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. |
---|
438 | |
---|
439 | \smallskip |
---|
440 | Similar logic applies to the other optional fields: \texttt{ff\_f} and \texttt{ff\_t} |
---|
441 | which can be used to provide the Coriolis parameter at F- and T-points respectively if the |
---|
442 | mesh is not on a sphere. If present these fields will be read and used and the normal |
---|
443 | calculation ($2*\Omega*\sin(\varphi)$) suppressed. Versions of \mdl{usrdef\_hgr} which set |
---|
444 | their own values of \texttt{ff\_f} and \texttt{ff\_t} should set the Coriolis computation |
---|
445 | flag: \texttt{iff} to a non-zero value to suppress their re-computation. |
---|
446 | |
---|
447 | Note that longitudes, latitudes, and scale factors at $w$ points are exactly |
---|
448 | equal to those of $t$ points, thus no specific arrays are defined at $w$ points. |
---|
449 | |
---|
450 | |
---|
451 | % ================================================================ |
---|
452 | % Domain: Vertical Grid (domzgr) |
---|
453 | % ================================================================ |
---|
454 | \subsection[Vertical grid (\textit{domzgr.F90})] |
---|
455 | {Vertical grid (\protect\mdl{domzgr})} |
---|
456 | \label{subsec:DOM_zgr} |
---|
457 | %-----------------------------------------namdom------------------------------------------- |
---|
458 | \nlst{namdom} |
---|
459 | %------------------------------------------------------------------------------------------------------------- |
---|
460 | |
---|
461 | In the vertical, the model mesh is determined by four things: |
---|
462 | \begin{enumerate} |
---|
463 | \item the bathymetry given in meters; |
---|
464 | \item the number of levels of the model (\jp{jpk}); |
---|
465 | \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and |
---|
466 | \item the masking system, \ie the number of wet model levels at each |
---|
467 | $(i,j)$ location of the horizontal grid. |
---|
468 | \end{enumerate} |
---|
469 | |
---|
470 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
471 | \begin{figure}[!tb] |
---|
472 | \begin{center} |
---|
473 | \includegraphics[width=\textwidth]{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 is made when setting up the configuration; |
---|
490 | it is not intended to be an option which can be changed in the middle of an |
---|
491 | experiment. The one exception to this statement being the choice of linear or |
---|
492 | non-linear free surface. In v4.0 the linear free surface option is implemented |
---|
493 | as a special case of the non-linear free surface. This is computationally |
---|
494 | wasteful since it uses the structures for time-varying 3D metrics for fields |
---|
495 | that (in the linear free surface case) are fixed. However, the linear |
---|
496 | free-surface is rarely used and implementing it this way means a single configuration |
---|
497 | file can support both options. |
---|
498 | |
---|
499 | By default a non-linear free surface is used (\np{ln\_linssh} set to \forcode{ = |
---|
500 | .false.} in \ngn{namdom}): the coordinate follow the time-variation of the free |
---|
501 | surface so that the transformation is time dependent: $z(i,j,k,t)$ |
---|
502 | (\eg \autoref{fig:z_zps_s_sps}f). When a linear free surface is assumed |
---|
503 | (\np{ln\_linssh} set to \forcode{ = .true.} in \ngn{namdom}), the vertical |
---|
504 | coordinates are fixed in time, but the seawater can move up and down across the |
---|
505 | $z_0$ surface (in other words, the top of the ocean in not a rigid lid). |
---|
506 | |
---|
507 | Note that settings: \np{ln\_zco}, \np{ln\_zps}, \np{ln\_sco} and \np{ln\_isfcav} mentioned |
---|
508 | in the following sections appear to be namelist options but they are no longer truly |
---|
509 | namelist options for NEMO. Their value is written to and read from the domain configuration file |
---|
510 | and they should be treated as fixed parameters for a particular configuration. They are |
---|
511 | namelist options for the \forcode{DOMAINcfg} tool that can be used to build the |
---|
512 | configuration file and serve both to provide a record of the choices made whilst building the |
---|
513 | configuration and to trigger appropriate code blocks within NEMO. |
---|
514 | These values should not be altered in the \np{cn\_domcfg} file. |
---|
515 | |
---|
516 | \medskip |
---|
517 | The decision on these choices must be made when the \np{cn\_domcfg} file is constructed. |
---|
518 | Three main choices are offered (\autoref{fig:z_zps_s_sps}a-c): |
---|
519 | |
---|
520 | \begin{itemize} |
---|
521 | \item $z$-coordinate with full step bathymetry (\np{ln\_zco}\forcode{ = .true.}), |
---|
522 | \item $z$-coordinate with partial step ($zps$) bathymetry (\np{ln\_zps}\forcode{ = .true.}), |
---|
523 | \item Generalized, $s$-coordinate (\np{ln\_sco}\forcode{ = .true.}). |
---|
524 | \end{itemize} |
---|
525 | |
---|
526 | Additionally, hybrid combinations of the three main coordinates are available: |
---|
527 | $s-z$ or $s-zps$ coordinate (\autoref{fig:z_zps_s_sps}d and \autoref{fig:z_zps_s_sps}e). |
---|
528 | |
---|
529 | A further choice related to vertical coordinate concerns the presence (or not) of ocean |
---|
530 | cavities beneath ice shelves within the model domain. A setting of \np{ln\_isfcav} as |
---|
531 | \forcode{.true.} indicates that the domain contains ocean cavities, otherwise the top, |
---|
532 | wet layer of the ocean will always be at the ocean surface. This option is currently only |
---|
533 | available for $z$- or $zps$-coordinates. In the latter case, partial steps are also applied |
---|
534 | at the ocean/ice shelf interface. |
---|
535 | |
---|
536 | Within the model, the arrays describing the grid point depths and vertical scale factors |
---|
537 | are three set of three dimensional arrays $(i,j,k)$ defined at \textit{before}, |
---|
538 | \textit{now} and \textit{after} time step. The time at which they are defined is |
---|
539 | indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. They are updated at each |
---|
540 | model time step. The initial fixed reference coordinate system is held in variable names |
---|
541 | with a $\_0$ suffix. When the linear free surface option is used |
---|
542 | (\np{ln\_linssh}\forcode{ = .true.}), \textit{before}, \textit{now} and \textit{after} |
---|
543 | arrays are initially set to their reference counterpart and remain fixed. |
---|
544 | |
---|
545 | \subsubsection{Required fields} |
---|
546 | \label{sec:DOM_zgr_fields} |
---|
547 | The explicit specification of a range of fields related to the vertical grid are required for the definition of a configuration. These include: |
---|
548 | |
---|
549 | \begin{Verbatim}[fontsize=\tiny] |
---|
550 | int ln_zco, ln_zps, ln_sco /* flags for z-coord, z-coord with partial steps and s-coord */ |
---|
551 | int ln_isfcav /* flag for ice shelf cavities */ |
---|
552 | double e3t_1d, e3w_1d /* reference vertical scale factors at T and W points */ |
---|
553 | double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 /* vertical scale factors 3D coordinate at T,U,V,F and W points */ |
---|
554 | double e3uw_0, e3vw_0 /* vertical scale factors 3D coordinate at UW and VW points */ |
---|
555 | int bottom_level, top_level /* last wet T-points, 1st wet T-points (for ice shelf cavities) */ |
---|
556 | /* For reference: */ |
---|
557 | float bathy_metry /* bathymetry used in setting top and bottom levels */ |
---|
558 | \end{Verbatim} |
---|
559 | |
---|
560 | This set of vertical metrics is sufficient to describe the initial depth and thickness of |
---|
561 | every gridcell in the model regardless of the choice of vertical coordinate. With constant |
---|
562 | z-levels, e3 metrics will be uniform across each horizontal level. In the partial step |
---|
563 | case each e3 at the \np{bottom\_level} (and, possibly, \np{top\_level} if ice cavities are |
---|
564 | present) may vary from its horizontal neighbours. And, in s-coordinates, variations can |
---|
565 | occur throughout the water column. With the non-linear free-surface, all the coordinates |
---|
566 | behave more like the s-coordinate in that variations occurr throughout the water column |
---|
567 | with displacements related to the sea surface height. These variations are typically much |
---|
568 | smaller than those arising from bottom fitted coordinates. The values for vertical metrics |
---|
569 | supplied in the domain configuration file can be considered as those arising from a flat |
---|
570 | sea surface with zero elevation. |
---|
571 | |
---|
572 | The \np{bottom\_level} and \np{top\_level} 2D arrays define the \np{bottom\_level} and top |
---|
573 | wet levels in each grid column. Without ice cavities, \np{top\_level} is essentially a land |
---|
574 | mask (0 on land; 1 everywhere else). With ice cavities, \np{top\_level} determines the |
---|
575 | first wet point below the overlying ice shelf. |
---|
576 | |
---|
577 | |
---|
578 | |
---|
579 | % ------------------------------------------------------------------------------------------------------------- |
---|
580 | % level bathymetry and mask |
---|
581 | % ------------------------------------------------------------------------------------------------------------- |
---|
582 | \subsubsection{Level bathymetry and mask} |
---|
583 | \label{subsec:DOM_msk} |
---|
584 | |
---|
585 | |
---|
586 | From \np{top\_level} and \np{bottom\_level} fields, the mask fields are defined as follows: |
---|
587 | \begin{alignat*}{2} |
---|
588 | tmask(i,j,k) &= & & |
---|
589 | \begin{cases} |
---|
590 | 0 &\text{if $ k < top\_level(i,j)$} \\ |
---|
591 | 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ |
---|
592 | 0 &\text{if $ k > bottom\_level(i,j)$} |
---|
593 | \end{cases} |
---|
594 | \\ |
---|
595 | umask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ |
---|
596 | vmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j + 1,k) \\ |
---|
597 | fmask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ |
---|
598 | & &* &tmask(i,j,k) * tmask(i + 1,j, k) \\ |
---|
599 | wmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j,k - 1) \\ |
---|
600 | \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) |
---|
601 | \end{alignat*} |
---|
602 | |
---|
603 | Note that, without ice shelves cavities, |
---|
604 | masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). |
---|
605 | Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) |
---|
606 | exactly in the same way as for the bottom boundary. |
---|
607 | |
---|
608 | %% The specification of closed lateral boundaries requires that at least |
---|
609 | %% the first and last rows and columns of the \textit{mbathy} array are set to zero. |
---|
610 | %% In the particular case of an east-west cyclical boundary condition, \textit{mbathy} has its last column equal to |
---|
611 | %% the second one and its first column equal to the last but one (and so too the mask arrays) |
---|
612 | %% (see \autoref{fig:LBC_jperio}). |
---|
613 | |
---|
614 | |
---|
615 | %------------------------------------------------------------------------------------------------- |
---|
616 | % Closed seas |
---|
617 | %------------------------------------------------------------------------------------------------- |
---|
618 | \subsection{Closed seas} \label{subsec:DOM_closea} |
---|
619 | |
---|
620 | When a global ocean is coupled to an atmospheric model it is better to represent all large |
---|
621 | water bodies (\eg great lakes, Caspian sea...) even if the model resolution does not allow |
---|
622 | their communication with the rest of the ocean. This is unnecessary when the ocean is |
---|
623 | forced by fixed atmospheric conditions, so these seas can be removed from the ocean |
---|
624 | domain. The user has the option to set the bathymetry in closed seas to zero (see |
---|
625 | \autoref{sec:MISC_closea}) and to optionally decide on the fate of any freshwater |
---|
626 | imbalance over the area. The options are explained in \autoref{sec:MISC_closea} but it |
---|
627 | should be noted here that a successful use of these options requires appropriate mask |
---|
628 | fields to be present in the domain configuration file. Among the possibilities are: |
---|
629 | |
---|
630 | \begin{Verbatim}[fontsize=\tiny] |
---|
631 | int closea_mask /* non-zero values in closed sea areas for optional masking */ |
---|
632 | int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ |
---|
633 | int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp) */ |
---|
634 | \end{Verbatim} |
---|
635 | |
---|
636 | % ------------------------------------------------------------------------------------------------------------- |
---|
637 | % Grid files |
---|
638 | % ------------------------------------------------------------------------------------------------------------- |
---|
639 | \subsection{Output grid files} |
---|
640 | \label{subsec:DOM_meshmask} |
---|
641 | |
---|
642 | \nlst{namcfg} |
---|
643 | |
---|
644 | Most of the arrays relating to a particular ocean model configuration dicussed in this |
---|
645 | chapter (grid-point position, scale factors) can be saved in a file if namelist parameter |
---|
646 | \np{ln\_write\_cfg} (namelist \ngn{namcfg}) is set to \forcode{.true.}; the output |
---|
647 | filename is set thorugh parameter \np{cn\_domcfg\_out}. This is only really useful |
---|
648 | if the fields are computed in subroutines \mdl{usrdef\_hgr} or \mdl{usrdef\_zgr} and |
---|
649 | checking or confirmation is required. |
---|
650 | |
---|
651 | \nlst{namdom} |
---|
652 | |
---|
653 | Alternatively, all the arrays relating to a particular ocean model configuration |
---|
654 | (grid-point position, scale factors, depths and masks) can be saved in a file called |
---|
655 | \texttt{mesh\_mask} if namelist parameter \np{ln\_meshmask} (namelist \ngn{namdom}) is set |
---|
656 | to \forcode{.true.}. This file contains additional fields that can be useful for |
---|
657 | post-processing applications |
---|
658 | |
---|
659 | % ================================================================ |
---|
660 | % Domain: Initial State (dtatsd & istate) |
---|
661 | % ================================================================ |
---|
662 | \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})] |
---|
663 | {Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} |
---|
664 | \label{sec:DTA_tsd} |
---|
665 | %-----------------------------------------namtsd------------------------------------------- |
---|
666 | \nlst{namtsd} |
---|
667 | %------------------------------------------------------------------------------------------ |
---|
668 | |
---|
669 | Basic initial state options are defined in \ngn{namtsd}. By default, the ocean starts |
---|
670 | from rest (the velocity field is set to zero) and the initialization of temperature and |
---|
671 | salinity fields is controlled through the \np{ln\_tsd\_init} namelist parameter. |
---|
672 | |
---|
673 | \begin{description} |
---|
674 | \item[\np{ln\_tsd\_init}\forcode{= .true.}] |
---|
675 | Use T and S input files that can be given on the model grid itself or on their native |
---|
676 | input data grids. In the latter case, the data will be interpolated on-the-fly both in |
---|
677 | the horizontal and the vertical to the model grid (see \autoref{subsec:SBC_iof}). The |
---|
678 | information relating to the input files are specified in the \np{sn\_tem} and |
---|
679 | \np{sn\_sal} structures. The computation is done in the \mdl{dtatsd} module. |
---|
680 | \item[\np{ln\_tsd\_init}\forcode{= .false.}] |
---|
681 | Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine |
---|
682 | contained in \mdl{userdef\_istate}. The default version sets horizontally uniform T and |
---|
683 | profiles as used in the GYRE configuration (see \autoref{sec:CFG_gyre}). |
---|
684 | \end{description} |
---|
685 | |
---|
686 | \biblio |
---|
687 | |
---|
688 | \pindex |
---|
689 | |
---|
690 | \end{document} |
---|