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