1 | |
---|
2 | % ================================================================ |
---|
3 | % Chapter 2 Ñ Space and Time Domain (DOM) |
---|
4 | % ================================================================ |
---|
5 | \chapter{Space and Time Domain (DOM) } |
---|
6 | \label{DOM} |
---|
7 | \minitoc |
---|
8 | |
---|
9 | % Missing things: |
---|
10 | % - istate: description of the initial state ==> this has to be put elsewhere.. |
---|
11 | % perhaps in MISC ? By the way the initialisation of T S and dynamics |
---|
12 | % should be put outside of DOM routine (better with TRC staff and off-line |
---|
13 | % tracers) |
---|
14 | % - daymod: definition of the time domain (nit000, nitend andd the calendar) |
---|
15 | % -geo2ocean: how to switch from geographic to mesh coordinate |
---|
16 | % - domclo: closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled |
---|
17 | |
---|
18 | \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction |
---|
19 | which could be referred to here, would help ==> to be added} |
---|
20 | %%%% |
---|
21 | |
---|
22 | |
---|
23 | \newpage |
---|
24 | $\ $\newline % force a new ligne |
---|
25 | |
---|
26 | |
---|
27 | Having defined the continuous equations in Chap.~\ref{PE}, we need to choose a |
---|
28 | discretization on a grid, and numerical algorithms. In the present chapter, we |
---|
29 | provide a general description of the staggered grid used in \NEMO, and other |
---|
30 | information relevant to the main directory routines (time stepping, main program) |
---|
31 | as well as the DOM (DOMain) directory. |
---|
32 | |
---|
33 | % ================================================================ |
---|
34 | % Fundamentals of the Discretisation |
---|
35 | % ================================================================ |
---|
36 | \section{Fundamentals of the Discretisation} |
---|
37 | \label{DOM_basics} |
---|
38 | |
---|
39 | % ------------------------------------------------------------------------------------------------------------- |
---|
40 | % Arrangement of Variables |
---|
41 | % ------------------------------------------------------------------------------------------------------------- |
---|
42 | \subsection{Arrangement of Variables} |
---|
43 | \label{DOM_cell} |
---|
44 | |
---|
45 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
46 | \begin{figure}[!tb] \label{Fig_cell} \begin{center} |
---|
47 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} |
---|
48 | \caption{Arrangement of variables. $T$ indicates scalar points where temperature, |
---|
49 | salinity, density, pressure and horizontal divergence are defined. ($u$,$v$,$w$) |
---|
50 | indicates vector points, and $f$ indicates vorticity points where both relative and |
---|
51 | planetary vorticities are defined} |
---|
52 | \end{center} \end{figure} |
---|
53 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
54 | |
---|
55 | The numerical techniques used to solve the Primitive Equations in this model are |
---|
56 | based on the traditional, centred second-order finite difference approximation. |
---|
57 | Special attention has been given to the homogeneity of the solution in the three |
---|
58 | space directions. 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 |
---|
60 | points $(u, v, w)$ defined in the centre of each face of the cells (Fig. \ref{Fig_cell}). |
---|
61 | This is the generalisation to three dimensions of the well-known ``C'' grid in |
---|
62 | Arakawa's classification \citep{Mesinger_Arakawa_Bk76}. The relative and |
---|
63 | planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge |
---|
64 | and the barotropic stream function $\psi$ is defined at horizontal points overlying |
---|
65 | the $\zeta$ and $f$-points. |
---|
66 | |
---|
67 | The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined |
---|
68 | by the transformation that gives ($\lambda$ ,$\varphi$ ,$z$) as a function of $(i,j,k)$. |
---|
69 | The grid-points are located at integer or integer and a half value of $(i,j,k)$ as |
---|
70 | indicated on Table \ref{Tab_cell}. In all the following, subscripts $u$, $v$, $w$, |
---|
71 | $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where the scale |
---|
72 | factors are defined. Each scale factor is defined as the local analytical value |
---|
73 | provided by \eqref{Eq_scale_factors}. As a result, the mesh on which partial |
---|
74 | derivatives $\frac{\partial}{\partial \lambda}, \frac{\partial}{\partial \varphi}$, and |
---|
75 | $\frac{\partial}{\partial z} $ are evaluated is a uniform mesh with a grid size of unity. |
---|
76 | Discrete partial derivatives are formulated by the traditional, centred second order |
---|
77 | finite difference approximation while the scale factors are chosen equal to their |
---|
78 | local analytical value. An important point here is that the partial derivative of the |
---|
79 | scale factors must be evaluated by centred finite difference approximation, not |
---|
80 | from their analytical expression. This preserves the symmetry of the discrete set |
---|
81 | of equations and therefore satisfies many of the continuous properties (see |
---|
82 | Appendix~\ref{Apdx_C}). A similar, related remark can be made about the domain |
---|
83 | size: when needed, an area, volume, or the total ocean depth must be evaluated |
---|
84 | as the sum of the relevant scale factors (see \eqref{DOM_bar}) in the next section). |
---|
85 | |
---|
86 | \begin{table}[!tb] \label{Tab_cell} |
---|
87 | \begin{center} \begin{tabular}{|p{46pt}|p{56pt}|p{56pt}|p{56pt}|} |
---|
88 | \hline |
---|
89 | T &$i$ & $j$ & $k$ \\ \hline |
---|
90 | u & $i+1/2$ & $j$ & $k$ \\ \hline |
---|
91 | v & $i$ & $j+1/2$ & $k$ \\ \hline |
---|
92 | w & $i$ & $j$ & $k+1/2$ \\ \hline |
---|
93 | f & $i+1/2$ & $j+1/2$ & $k$ \\ \hline |
---|
94 | uw & $i+1/2$ & $j$ & $k+1/2$ \\ \hline |
---|
95 | vw & $i$ & $j+1/2$ & $k+1/2$ \\ \hline |
---|
96 | fw & $i+1/2$ & $j+1/2$ & $k+1/2$ \\ \hline |
---|
97 | \end{tabular} |
---|
98 | \caption{Location of grid-points as a function of integer or integer and a half value |
---|
99 | of the column, line or level. This indexing is only used for the writing of the semi- |
---|
100 | discrete equation. In the code, the indexing uses integer values only and has a |
---|
101 | reverse direction in the vertical (see \S\ref{DOM_Num_Index})} |
---|
102 | \end{center} |
---|
103 | \end{table} |
---|
104 | |
---|
105 | % ------------------------------------------------------------------------------------------------------------- |
---|
106 | % Vector Invariant Formulation |
---|
107 | % ------------------------------------------------------------------------------------------------------------- |
---|
108 | \subsection{Discrete Operators} |
---|
109 | \label{DOM_operators} |
---|
110 | |
---|
111 | Given the values of a variable $q$ at adjacent points, the differencing and |
---|
112 | averaging operators at the midpoint between them are: |
---|
113 | \begin{subequations} \label{Eq_di_mi} |
---|
114 | \begin{align} |
---|
115 | \delta _i [q] &= \ \ q(i+1/2) - q(i-1/2) \\ |
---|
116 | \overline q^{\,i} &= \left\{ q(i+1/2) + q(i-1/2) \right\} \; / \; 2 |
---|
117 | \end{align} |
---|
118 | \end{subequations} |
---|
119 | |
---|
120 | Similar operators are defined with respect to $i+1/2$, $j$, $j+1/2$, $k$, and |
---|
121 | $k+1/2$. Following \eqref{Eq_PE_grad} and \eqref{Eq_PE_lap}, the gradient of a |
---|
122 | variable $q$ defined at a $T$-point has its three components defined at $u$-, $v$- |
---|
123 | and $w$-points while its Laplacien is defined at $T$-point. These operators have |
---|
124 | the following discrete forms in the curvilinear $s$-coordinate system: |
---|
125 | \begin{equation} \label{Eq_DOM_grad} |
---|
126 | \nabla q\equiv \frac{1}{e_{1u} }\delta _{i+1/2} \left[ q \right]\;\,{\rm {\bf i}} |
---|
127 | + \frac{1}{e_{2v} }\delta _{j+1/2} \left[ q \right]\;\,{\rm {\bf j}} |
---|
128 | + \frac{1}{e_{3w} }\delta _{k+1/2} \left[ q \right]\;\,{\rm {\bf k}} |
---|
129 | \end{equation} |
---|
130 | \begin{multline} \label{Eq_DOM_lap} |
---|
131 | \Delta q\equiv \frac{1}{e_{1T} {\kern 1pt}e_{2T} {\kern 1pt}e_{3T} }\;\left( |
---|
132 | {\delta _i \left[ {\frac{e_{2u} e_{3u} }{e_{1u} }\;\delta _{i+1/2} |
---|
133 | \left[ q \right]} \right] |
---|
134 | +\delta _j \left[ {\frac{e_{1v} e_{3v} }{e_{2v} |
---|
135 | }\;\delta _{j+1/2} \left[ q \right]} \right]\;} \right) \\ |
---|
136 | +\frac{1}{e_{3T} }\delta _k \left[ {\frac{1}{e_{3w} }\;\delta _{k+1/2} |
---|
137 | \left[ q \right]} \right] |
---|
138 | \end{multline} |
---|
139 | |
---|
140 | Following \eqref{Eq_PE_curl} and \eqref{Eq_PE_div}, a vector ${\rm {\bf A}}=\left( a_1,a_2,a_3\right)$ defined at vector points $(u,v,w)$ has its three curl |
---|
141 | components defined at $vw$-, $uw$, and $f$-points, and its divergence defined |
---|
142 | at $T$-points: |
---|
143 | \begin{equation} \label{Eq_DOM_curl} |
---|
144 | \begin{split} |
---|
145 | \nabla \times {\rm {\bf A}}\equiv \frac{1}{e_{2v} {\kern 1pt}e_{3vw} |
---|
146 | }{\kern 1pt}\,\;\left( {\delta _{j+1/2} \left[ {e_{3w} a_3 } \right]-\delta |
---|
147 | _{k+1/2} \left[ {e_{2v} a_2 } \right]} \right) &\;\;{\rm {\bf i}} \\ |
---|
148 | +\frac{1}{e_{2u} {\kern 1pt}e_{3uw} }\;\left( {\delta _{k+1/2} \left[ {e_{1u} a_1 } |
---|
149 | \right]-\delta _{i+1/2} \left[ {e_{3w} a_3 } \right]} \right) &\;\;{\rm{\bf j}} \\ |
---|
150 | +\frac{e_{3f} }{e_{1f} {\kern 1pt}e_{2f} }\,{\kern 1pt}\;\left( {\delta |
---|
151 | _{i+1/2} \left[ {e_{2v} a_2 } \right]-\delta _{j+12} \left[ {e_{1u} a_1 } \right]} |
---|
152 | \right) &\;\;{\rm {\bf k}} |
---|
153 | \end{split} |
---|
154 | \end{equation} |
---|
155 | \begin{equation} \label{Eq_DOM_div} |
---|
156 | \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} e_{3T} }\left( {\delta |
---|
157 | _i \left[ {e_{2u} e_{3u} a_1 } \right]+\delta _j \left[ {e_{1v} e_{3v} a_2 } |
---|
158 | \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] |
---|
159 | \end{equation} |
---|
160 | |
---|
161 | In the special case of a pure $z$-coordinate system, \eqref{Eq_DOM_lap} and |
---|
162 | \eqref{Eq_DOM_div} can be simplified. In this case, the vertical scale factor |
---|
163 | becomes a function of the single variable $k$ and thus does not depend on the |
---|
164 | horizontal location of a grid point. For example \eqref{Eq_DOM_div} reduces to: |
---|
165 | \begin{equation*} |
---|
166 | \nabla \cdot {\rm {\bf A}}=\frac{1}{e_{1T} e_{2T} }\left( {\delta |
---|
167 | _i \left[ {e_{2u} a_1 } \right]+\delta _j \left[ {e_{1v} a_2 } |
---|
168 | \right]} \right)+\frac{1}{e_{3T} }\delta _k \left[ {a_3 } \right] |
---|
169 | \end{equation*} |
---|
170 | |
---|
171 | The vertical average over the whole water column denoted by an overbar becomes |
---|
172 | for a quantity $q$ which is a masked field (i.e. equal to zero inside solid area): |
---|
173 | \begin{equation} \label{DOM_bar} |
---|
174 | \bar q = \frac{1}{H}\int_{k^b}^{k^o} {q\;e_{3q} \,dk} |
---|
175 | \equiv \frac{1}{H_q }\sum\limits_k {q\;e_{3q} } |
---|
176 | \end{equation} |
---|
177 | where $H_q$ is the ocean depth, which is the masked sum of the vertical scale |
---|
178 | factors at $q$ points, $k^b$ and $k^o$ are the bottom and surface $k$-indices, |
---|
179 | and the symbol $k^o$ refers to a summation over all grid points of the same type |
---|
180 | in the direction indicated by the subscript (here $k$). |
---|
181 | |
---|
182 | In continuous form, the following properties are satisfied: |
---|
183 | \begin{equation} \label{Eq_DOM_curl_grad} |
---|
184 | \nabla \times \nabla q ={\rm {\bf {0}}} |
---|
185 | \end{equation} |
---|
186 | \begin{equation} \label{Eq_DOM_div_curl} |
---|
187 | \nabla \cdot \left( {\nabla \times {\rm {\bf A}}} \right)=0 |
---|
188 | \end{equation} |
---|
189 | |
---|
190 | It is straightforward to demonstrate that these properties are verified locally in |
---|
191 | discrete form as soon as the scalar $q$ is taken at $T$-points and the vector |
---|
192 | \textbf{A} has its components defined at vector points $(u,v,w)$. |
---|
193 | |
---|
194 | Let $a$ and $b$ be two fields defined on the mesh, with value zero inside |
---|
195 | continental area. Using integration by parts it can be shown that the differencing |
---|
196 | operators ($\delta_i$, $\delta_j$ and $\delta_k$) are anti-symmetric linear |
---|
197 | operators, and further that the averaging operators $\overline{\,\cdot\,}^{\,i}$, |
---|
198 | $\overline{\,\cdot\,}^{\,k}$ and $\overline{\,\cdot\,}^{\,k}$) are symmetric linear |
---|
199 | operators, $i.e.$ |
---|
200 | \begin{align} |
---|
201 | \label{DOM_di_adj} |
---|
202 | \sum\limits_i { a_i \;\delta _i \left[ b \right]} |
---|
203 | &\equiv -\sum\limits_i {\delta _{i+1/2} \left[ a \right]\;b_{i+1/2} } \\ |
---|
204 | \label{DOM_mi_adj} |
---|
205 | \sum\limits_i { a_i \;\overline b^{\,i}} |
---|
206 | & \equiv \quad \sum\limits_i {\overline a ^{\,i+1/2}\;b_{i+1/2} } |
---|
207 | \end{align} |
---|
208 | |
---|
209 | In other words, the adjoint of the differencing and averaging operators are |
---|
210 | $\delta_i^*=\delta_{i+1/2}$ and |
---|
211 | ${(\overline{\,\cdot \,}^{\,i})}^*= \overline{\,\cdot\,}^{\,i+1/2}$, respectively. |
---|
212 | These two properties will be used extensively in the Appendix~\ref{Apdx_C} to |
---|
213 | demonstrate integral conservative properties of the discrete formulation chosen. |
---|
214 | |
---|
215 | % ------------------------------------------------------------------------------------------------------------- |
---|
216 | % Numerical Indexing |
---|
217 | % ------------------------------------------------------------------------------------------------------------- |
---|
218 | \subsection{Numerical Indexing} |
---|
219 | \label{DOM_Num_Index} |
---|
220 | |
---|
221 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
222 | \begin{figure}[!tb] \label{Fig_index_hor} \begin{center} |
---|
223 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf} |
---|
224 | \caption{Horizontal integer indexing used in the \textsc{Fortran} code. The dashed |
---|
225 | area indicates the cell in which variables contained in arrays have the same |
---|
226 | $i$- and $j$-indices} |
---|
227 | \end{center} \end{figure} |
---|
228 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
229 | |
---|
230 | The array representation used in the \textsc{Fortran} code requires an integer |
---|
231 | indexing while the analytical definition of the mesh (see \S\ref{DOM_cell}) is |
---|
232 | associated with the use of integer values for $T$-points and both integer and |
---|
233 | integer and a half values for all the other points. Therefore a specific integer |
---|
234 | indexing must be defined for points other than $T$-points ($i.e.$ velocity and |
---|
235 | vorticity grid-points). Furthermore, the direction of the vertical indexing has |
---|
236 | been changed so that the surface level is at $k=1$. |
---|
237 | |
---|
238 | % ----------------------------------- |
---|
239 | % Horizontal Indexing |
---|
240 | % ----------------------------------- |
---|
241 | \subsubsection{Horizontal Indexing} |
---|
242 | \label{DOM_Num_Index_hor} |
---|
243 | |
---|
244 | The indexing in the horizontal plane has been chosen as shown in Fig.\ref{Fig_index_hor}. For an increasing $i$ index ($j$ index), the $T$-point |
---|
245 | and the eastward $u$-point (northward $v$-point) have the same index |
---|
246 | (see the dashed area in Fig.\ref{Fig_index_hor}). A $T$-point and its |
---|
247 | nearest northeast $f$-point have the same $i$-and $j$-indices. |
---|
248 | |
---|
249 | % ----------------------------------- |
---|
250 | % Vertical indexing |
---|
251 | % ----------------------------------- |
---|
252 | \subsubsection{Vertical Indexing} |
---|
253 | \label{DOM_Num_Index_vertical} |
---|
254 | |
---|
255 | In the vertical, the chosen indexing requires special attention since the |
---|
256 | $k$-axis is re-orientated downward in the \textsc{Fortran} code compared |
---|
257 | to the indexing used in the semi-discrete equations and given in \S\ref{DOM_cell}. |
---|
258 | The sea surface corresponds to the $w$-level $k=1$ which is the same index |
---|
259 | as $T$-level just below (Fig.\ref{Fig_index_vert}). The last $w$-level ($k=jpk$) |
---|
260 | either corresponds to the ocean floor or is inside the bathymetry while the last |
---|
261 | $T$-level is always inside the bathymetry (Fig.\ref{Fig_index_vert}). Note that |
---|
262 | for an increasing $k$ index, a $w$-point and the $T$-point just below have the |
---|
263 | same $k$ index, 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 |
---|
265 | axis that have the same $i$ or $j$ index (compare the dashed area in |
---|
266 | Fig.\ref{Fig_index_hor} and \ref{Fig_index_vert}). Since the scale factors are |
---|
267 | chosen to be strictly positive, a \emph{minus sign} appears in the \textsc{Fortran} |
---|
268 | code \emph{before all the vertical derivatives} of the discrete equations given in |
---|
269 | this documentation. |
---|
270 | |
---|
271 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
272 | \begin{figure}[!pt] \label{Fig_index_vert} \begin{center} |
---|
273 | \includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf} |
---|
274 | \caption{Vertical integer indexing used in the \textsc{Fortran } code. Note that |
---|
275 | the $k$-axis is orientated downward. The dashed area indicates the cell in |
---|
276 | which variables contained in arrays have the same $k$-index.} |
---|
277 | \end{center} \end{figure} |
---|
278 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
279 | |
---|
280 | % ----------------------------------- |
---|
281 | % Domain Size |
---|
282 | % ----------------------------------- |
---|
283 | \subsubsection{Domain Size} |
---|
284 | \label{DOM_size} |
---|
285 | |
---|
286 | The total size of the computational domain is set by the parameters \jp{jpiglo}, |
---|
287 | \jp{jpjglo} and \jp{jpk} in the $i$, $j$ and $k$ directions respectively. They are |
---|
288 | given as parameters in the \mdl{par\_oce} module\footnote{When a specific |
---|
289 | configuration is used (ORCA2 global ocean, etc...) the parameter are actually |
---|
290 | defined in additional files introduced by \mdl{par\_oce} module via CPP |
---|
291 | \textit{include} command. For example, ORCA2 parameters are set in |
---|
292 | \textit{par\_ORCA\_R2.h90} file}. The use of parameters rather than variables |
---|
293 | (together with dynamic allocation of arrays) was chosen because it ensured that |
---|
294 | the compiler would optimize the executable code efficiently, especially on vector |
---|
295 | machines (optimization may be less efficient when the problem size is unknown |
---|
296 | at the time of compilation). Nevertheless, it is possible to set up the code with full |
---|
297 | dynamical allocation by using the AGRIF packaged \citep{Debreu_al_CG2008}. |
---|
298 | % |
---|
299 | \gmcomment{ add the following ref |
---|
300 | \colorbox{yellow}{(ref part of the doc)} } |
---|
301 | % |
---|
302 | Note that are other parameters in \mdl{par\_oce} that refer to the domain size. |
---|
303 | The two parameters $jpidta$ and $jpjdta$ may be larger than $jpiglo$, $jpjglo$ |
---|
304 | when the user wants to use only a sub-region of a given configuration. This is |
---|
305 | the "zoom" capability described in \S\ref{MISC_zoom}. In most applications of |
---|
306 | the model, $jpidta=jpiglo$, $jpjdta=jpjglo$, and $jpizoom=jpjzoom=1$. Parameters |
---|
307 | $jpi$ and $jpj$ refer to the size of each processor subdomain when the code is |
---|
308 | run in parallel using domain decomposition (\key{mpp\_mpi} defined, see |
---|
309 | \S\ref{LBC_mpp}). |
---|
310 | |
---|
311 | % ================================================================ |
---|
312 | % Domain: Horizontal Grid (mesh) |
---|
313 | % ================================================================ |
---|
314 | \section [Domain: Horizontal Grid (mesh) (\textit{domhgr})] |
---|
315 | {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} } |
---|
316 | \label{DOM_hgr} |
---|
317 | |
---|
318 | % ------------------------------------------------------------------------------------------------------------- |
---|
319 | % Coordinates and scale factors |
---|
320 | % ------------------------------------------------------------------------------------------------------------- |
---|
321 | \subsection{Coordinates and scale factors} |
---|
322 | \label{DOM_hgr_coord_e} |
---|
323 | |
---|
324 | The ocean mesh ($i.e.$ the position of all the scalar and vector points) is defined |
---|
325 | by the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. |
---|
326 | The grid-points are located at integer or integer and a half values of as indicated |
---|
327 | in Table~\ref{Tab_cell}. The associated scale factors are defined using the |
---|
328 | analytical first derivative of the transformation \eqref{Eq_scale_factors}. These |
---|
329 | definitions are done in two modules, \mdl{domhgr} and \mdl{domzgr}, which |
---|
330 | provide the horizontal and vertical meshes, respectively. This section deals with |
---|
331 | the horizontal mesh parameters. |
---|
332 | |
---|
333 | In a horizontal plane, the location of all the model grid points is defined from the |
---|
334 | analytical expressions of the longitude $\lambda$ and latitude $\varphi$ as a |
---|
335 | function of $(i,j)$. The horizontal scale factors are calculated using |
---|
336 | \eqref{Eq_scale_factors}. For example, when the longitude and latitude are |
---|
337 | function of a single value ($i$ and $j$, respectively) (geographical configuration |
---|
338 | of the mesh), the horizontal mesh definition reduces to define the wanted |
---|
339 | $\lambda(i)$, $\varphi(j)$, and their derivatives $\lambda'(i)$ $\varphi'(j)$ in the |
---|
340 | \mdl{domhgr} module. The model computes the grid-point positions and scale |
---|
341 | factors in the horizontal plane as follows: |
---|
342 | \begin{flalign*} |
---|
343 | \lambda_T &\equiv \text{glamt}= \lambda(i) & \varphi_T &\equiv \text{gphit} = \varphi(j)\\ |
---|
344 | \lambda_u &\equiv \text{glamu}= \lambda(i+1/2)& \varphi_u &\equiv \text{gphiu}= \varphi(j)\\ |
---|
345 | \lambda_v &\equiv \text{glamv}= \lambda(i) & \varphi_v &\equiv \text{gphiv} = \varphi(j+1/2)\\ |
---|
346 | \lambda_f &\equiv \text{glamf }= \lambda(i+1/2)& \varphi_f &\equiv \text{gphif }= \varphi(j+1/2) |
---|
347 | \end{flalign*} |
---|
348 | \begin{flalign*} |
---|
349 | e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j) |& |
---|
350 | e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)| \\ |
---|
351 | e_{1u} &\equiv \text{e1t} = r_a |\lambda'(i+1/2) \; \cos\varphi(j) |& |
---|
352 | e_{2u} &\equiv \text{e2t} = r_a |\varphi'(j)|\\ |
---|
353 | e_{1v} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j+1/2) |& |
---|
354 | e_{2v} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)|\\ |
---|
355 | e_{1f} &\equiv \text{e1t} = r_a |\lambda'(i+1/2)\; \cos\varphi(j+1/2) |& |
---|
356 | e_{2f} &\equiv \text{e2t} = r_a |\varphi'(j+1/2)| |
---|
357 | \end{flalign*} |
---|
358 | where the last letter of each computational name indicates the grid point |
---|
359 | considered and $r_a$ is the earth radius (defined in \mdl{phycst} along with |
---|
360 | all universal constants). Note that the horizontal position of and scale factors |
---|
361 | at $w$-points are exactly equal to those of $T$-points, thus no specific arrays |
---|
362 | are defined at $w$-points. |
---|
363 | |
---|
364 | Note that the definition of the scale factors ($i.e.$ as the analytical first derivative |
---|
365 | of the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$) is |
---|
366 | specific to the \NEMO model \citep{Marti1992}. As an example, $e_{1T}$ is defined |
---|
367 | locally at a $T$-point, whereas many other models on a C grid choose to define |
---|
368 | such a scale factor as the distance between the $U$-points on each side of the |
---|
369 | $T$-point. Relying on an analytical transformation has two advantages: firstly, there |
---|
370 | is no ambiguity in the scale factors appearing in the discrete equations, since they |
---|
371 | are first introduced in the continuous equations; secondly, analytical transformations |
---|
372 | encourage good practice by the definition of smoothly varying grids (rather than |
---|
373 | allowing the user to set arbitrary jumps in thickness between adjacent layers) |
---|
374 | \citep{Treguier1996}. An example of the effect of such a choice is shown in |
---|
375 | Fig.~\ref{Fig_zgr_e3}. |
---|
376 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
377 | \begin{figure}[!t] \label{Fig_zgr_e3} \begin{center} |
---|
378 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf} |
---|
379 | \caption{Comparison of (a) traditional definitions of grid-point position and grid-size |
---|
380 | in the vertical, and (b) analytically derived grid-point position and scale factors. For |
---|
381 | both grids here, the same $w$-point depth has been chosen but in (a) the |
---|
382 | $T$-points are set half way between $w$-points while in (b) they are defined from |
---|
383 | an analytical function: $z(k)=5\,(i-1/2)^3 - 45\,(i-1/2)^2 + 140\,(i-1/2) - 150$. |
---|
384 | Note the resulting difference between the value of the grid-size $\Delta_k$ and |
---|
385 | those of the scale factor $e_k$. } |
---|
386 | \end{center} \end{figure} |
---|
387 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
388 | |
---|
389 | % ------------------------------------------------------------------------------------------------------------- |
---|
390 | % Choice of horizontal grid |
---|
391 | % ------------------------------------------------------------------------------------------------------------- |
---|
392 | \subsection{Choice of horizontal grid} |
---|
393 | \label{DOM_hgr_msh_choice} |
---|
394 | |
---|
395 | The user has three options available in defining a horizontal grid, which involve |
---|
396 | the parameter $jphgr\_mesh$ of the \mdl{par\_oce} module. |
---|
397 | \begin{description} |
---|
398 | \item[\jp{jphgr\_mesh}=0] The most general curvilinear orthogonal grids. |
---|
399 | The coordinates and their first derivatives with respect to $i$ and $j$ are |
---|
400 | provided in a file, read in \rou{hgr\_read} subroutine of the domhgr module. |
---|
401 | \item[\jp{jphgr\_mesh}=1 to 5] A few simple analytical grids are provided (see below). |
---|
402 | For other analytical grids, the \mdl{domhgr} module must be modified by the user. |
---|
403 | \end{description} |
---|
404 | |
---|
405 | There are two simple cases of geographical grids on the sphere. With |
---|
406 | \jp{jphgr\_mesh}=1, the grid (expressed in degrees) is regular in space, |
---|
407 | with grid sizes specified by parameters \pp{ppe1\_deg} and \pp{ppe2\_deg}, |
---|
408 | respectively. Such a geographical grid can be very anisotropic at high latitudes |
---|
409 | because of the convergence of meridians (the zonal scale factors $e_1$ |
---|
410 | become much smaller than the meridional scale factors $e_2$). The Mercator |
---|
411 | grid (\jp{jphgr\_mesh}=4) avoids this anisotropy by refining the meridional scale |
---|
412 | factors in the same way as the zonal ones. In this case, meridional scale factors |
---|
413 | and latitudes are calculated analytically using the formulae appropriate for |
---|
414 | a Mercator projection, based on \pp{ppe1\_deg} which is a reference grid spacing |
---|
415 | at the equator (this applies even when the geographical equator is situated outside |
---|
416 | the model domain). |
---|
417 | %%% |
---|
418 | \gmcomment{ give here the analytical expression of the Mercator mesh} |
---|
419 | %%% |
---|
420 | In these two cases (\jp{jphgr\_mesh}=1 or 4), the grid position is defined by the |
---|
421 | longitude and latitude of the south-westernmost point (\pp{ppglamt0} |
---|
422 | and \pp{ppgphi0}). Note that for the Mercator grid the user need only provide |
---|
423 | an approximate starting latitude: the real latitude will be recalculated analytically, |
---|
424 | in order to ensure that the equator corresponds to line passing through $T$- |
---|
425 | and $u$-points. |
---|
426 | |
---|
427 | Rectangular grids ignoring the spherical geometry are defined with |
---|
428 | \jp{jphgr\_mesh} = 2, 3, 5. The domain is either an $f$-plane (\jp{jphgr\_mesh} = 2, |
---|
429 | Coriolis factor is constant) or a beta-plane (\jp{jphgr\_mesh} = 3, the Coriolis factor |
---|
430 | is linear in the $j$-direction). The grid size is uniform in meter in each direction, |
---|
431 | and given by the parameters \pp{ppe1\_m} and \pp{ppe2\_m} respectively. |
---|
432 | The zonal grid coordinate (\textit{glam} arrays) is in kilometers, starting at zero |
---|
433 | with the first $T$-point. The meridional coordinate (gphi. arrays) is in kilometers, |
---|
434 | and the second $T$-point corresponds to coordinate $gphit=0$. The input |
---|
435 | parameter \pp{ppglam0} is ignored. \pp{ppgphi0} is used to set the reference |
---|
436 | latitude for computation of the Coriolis parameter. In the case of the beta plane, |
---|
437 | \pp{ppgphi0} corresponds to the center of the domain. Finally, the special case |
---|
438 | \jp{jphgr\_mesh}=5 corresponds to a beta plane in a rotated domain for the |
---|
439 | GYRE configuration, representing a classical mid-latitude double gyre system. |
---|
440 | The rotation allows us to maximize the jet length relative to the gyre areas |
---|
441 | (and the number of grid points). |
---|
442 | |
---|
443 | The choice of the grid must be consistent with the boundary conditions specified |
---|
444 | by the parameter \jp{jperio} (see {\S\ref{LBC}). |
---|
445 | |
---|
446 | % ------------------------------------------------------------------------------------------------------------- |
---|
447 | % Grid files |
---|
448 | % ------------------------------------------------------------------------------------------------------------- |
---|
449 | \subsection{Grid files} |
---|
450 | \label{DOM_hgr_files} |
---|
451 | |
---|
452 | All the arrays relating to a particular ocean model configuration (grid-point |
---|
453 | position, scale factors, masks) can be saved in files if $\np{nmsh} \not= 0$ |
---|
454 | (namelist parameter). This can be particularly useful for plots and off-line |
---|
455 | diagnostics. In some cases, the user may choose to make a local modification |
---|
456 | of a scale factor in the code. This is the case in global configurations when |
---|
457 | restricting the width of a specific strait (usually a one-grid-point strait that |
---|
458 | happens to be too wide due to insufficient model resolution). An example |
---|
459 | is Gibraltar Strait in the ORCA2 configuration. When such modifications are done, |
---|
460 | the output grid written when $\np{nmsh} \not=0$ is no more equal to the input grid. |
---|
461 | |
---|
462 | % ================================================================ |
---|
463 | % Domain: Vertical Grid (domzgr) |
---|
464 | % ================================================================ |
---|
465 | \section [Domain: Vertical Grid (\textit{domzgr})] |
---|
466 | {Domain: Vertical Grid \small{(\mdl{domzgr} module)} } |
---|
467 | \label{DOM_zgr} |
---|
468 | %-----------------------------------------nam_zgr & namdom------------------------------------------- |
---|
469 | \namdisplay{nam_zgr} |
---|
470 | \namdisplay{namdom} |
---|
471 | %------------------------------------------------------------------------------------------------------------- |
---|
472 | |
---|
473 | In the vertical, the model mesh is determined by four things: |
---|
474 | (1) the bathymetry given in meters ; |
---|
475 | (2) the number of levels of the model (\jp{jpk}) ; |
---|
476 | (3) the analytical transformation $z(i,j,k)$ and the vertical scale factors |
---|
477 | (derivatives of the transformation) ; |
---|
478 | and (4) the masking system, $i.e.$ the number of wet model levels at each |
---|
479 | $(i,j)$ column of points. |
---|
480 | |
---|
481 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
482 | \begin{figure}[!tb] \label{Fig_z_zps_s_sps} \begin{center} |
---|
483 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf} |
---|
484 | \caption{The ocean bottom as seen by the model: |
---|
485 | (a) $z$-coordinate with full step, |
---|
486 | (b) $z$-coordinate with partial step, |
---|
487 | (c) $s$-coordinate: terrain following representation, |
---|
488 | (d) hybrid $s-z$ coordinate, |
---|
489 | (e) hybrid $s-z$ coordinate with partial step, and |
---|
490 | (f) same as (e) but with variable volume associated with the non-linear free surface. |
---|
491 | Note that the variable volume option (\key{vvl}) can be used with any of the |
---|
492 | 5 coordinates (a) to (e).} |
---|
493 | \end{center} \end{figure} |
---|
494 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
495 | |
---|
496 | The choice of a vertical coordinate, even if it is made through a namelist parameter, |
---|
497 | must be done once of all at the beginning of an experiment. It is not intended as an |
---|
498 | option which can be enabled or disabled in the middle of an experiment. Three main |
---|
499 | choices are offered (Fig.~\ref{Fig_z_zps_s_sps}a to c): $z$-coordinate with full step |
---|
500 | bathymetry (\np{ln\_zco}=true), $z$-coordinate with partial step bathymetry |
---|
501 | (\np{ln\_zps}=true), or generalized, $s$-coordinate (\np{ln\_sco}=true). |
---|
502 | Hybridation of the three main coordinates are available: $s-z$ or $s-zps$ coordinate |
---|
503 | (Fig.~\ref{Fig_z_zps_s_sps}d and \ref{Fig_z_zps_s_sps}e). When using the variable |
---|
504 | volume option \key{vvl}) ($i.e.$ non-linear free surface), the coordinate follow the |
---|
505 | time-variation of the free surface so that the transformation is time dependent: |
---|
506 | $z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step |
---|
507 | bathymetry or $s$-coordinate (hybride and partial step coordinates have not |
---|
508 | yet been tested in NEMO v2.3). |
---|
509 | |
---|
510 | Contrary to the horizontal grid, the vertical grid is computed in the code and no |
---|
511 | provision is made for reading it from a file. The only input file is the bathymetry |
---|
512 | (in meters)\footnote{N.B. in full step $z$-coordinate, a \textit{bathy\_level} file can |
---|
513 | replace the \textit{bathy\_meter} file, so that the computation of the number of |
---|
514 | wet ocean point in each water column is by-passed}. After reading the bathymetry, |
---|
515 | the algorithm for vertical grid definition differs between the different options: |
---|
516 | \begin{description} |
---|
517 | \item[\textit{zco}] set a reference coordinate transformation $z_0 (k)$, and set $z(i,j,k,t)=z_0 (k)$. |
---|
518 | \item[\textit{zps}] set a reference coordinate transformation $z_0 (k)$, and |
---|
519 | calculate the thickness of the deepest level at each $(i,j)$ point using the |
---|
520 | bathymetry, to obtain the final three-dimensional depth and scale factor arrays. |
---|
521 | \item[\textit{sco}] smooth the bathymetry to fulfil the hydrostatic consistency |
---|
522 | criteria and set the three-dimensional transformation. |
---|
523 | \item[\textit{s-z} and \textit{s-zps}] smooth the bathymetry to fulfil the hydrostatic |
---|
524 | consistency criteria and set the three-dimensional transformation $z(i,j,k)$, and |
---|
525 | possibly introduce masking of extra land points to better fit the original bathymetry file |
---|
526 | \end{description} |
---|
527 | %%% |
---|
528 | \gmcomment{ add the description of the smoothing: envelop topography...} |
---|
529 | %%% |
---|
530 | |
---|
531 | Generally, the arrays describing the grid point depths and vertical scale factors |
---|
532 | are three dimensional arrays $(i,j,k)$. In the special case of $z$-coordinates with |
---|
533 | full step bottom topography, it is possible to define those arrays as one-dimensional, |
---|
534 | in order to save memory. This is performed by defining the \key{zco} |
---|
535 | C-Pre-Processor (CPP) key. To improve the code readability while providing this |
---|
536 | flexibility, the vertical coordinate and scale factors are defined as functions of |
---|
537 | $(i,j,k)$ with "fs" as prefix (examples: \textit{fsdeptht, fse3t,} etc) that can be either |
---|
538 | three-dimensional arrays, or a one dimensional array when \key{zco} is defined. |
---|
539 | These functions are defined in the file \hf{domzgr\_substitute} of the DOM directory. |
---|
540 | They are used throughout the code, and replaced by the corresponding arrays at |
---|
541 | the time of pre-processing (CPP capability). |
---|
542 | |
---|
543 | % ------------------------------------------------------------------------------------------------------------- |
---|
544 | % Meter Bathymetry |
---|
545 | % ------------------------------------------------------------------------------------------------------------- |
---|
546 | \subsection{Meter Bathymetry} |
---|
547 | \label{DOM_bathy} |
---|
548 | |
---|
549 | Three options are possible for defining the bathymetry, according to the |
---|
550 | namelist variable \np{ntopo}: |
---|
551 | \begin{description} |
---|
552 | \item[\np{ntopo} = 0] a flat-bottom domain is defined. The total depth $z_w (jpk)$ |
---|
553 | is given by the coordinate transformation. The domain can either be a closed |
---|
554 | basin or a periodic channel depending on the parameter \jp{jperio}. |
---|
555 | \item[\np{ntopo} = -1] a domain with a bump of topography one third of the |
---|
556 | domain width at the central latitude. This is meant for the "EEL-R5" configuration, |
---|
557 | a periodic or open boundary channel with a seamount. |
---|
558 | \item[\np{ntopo} = 1] read a bathymetry. The bathymetry file (Netcdf format) |
---|
559 | provides the ocean depth (positive, in meters) at each grid point of the model grid. |
---|
560 | The bathymetry is usually built by interpolating a standard bathymetry product |
---|
561 | ($e.g.$ ETOPO2) onto the horizontal ocean mesh. Defining the bathymetry also |
---|
562 | defines the coastline: where the bathymetry is zero, no model levels are defined |
---|
563 | (all levels are masked). |
---|
564 | \end{description} |
---|
565 | |
---|
566 | When using the rigid lid approximation (\key{dynspg\_rl} is defined) isolated land |
---|
567 | masses (islands) must be identified by negative integers in the input bathymetry file |
---|
568 | (see \S\ref{MISC_solisl}). |
---|
569 | |
---|
570 | When a global ocean is coupled to an atmospheric model it is better to represent |
---|
571 | all large water bodies (e.g, great lakes, Caspian sea...) even if the model |
---|
572 | resolution does not allow their communication with the rest of the ocean. |
---|
573 | This is unnecessary when the ocean is forced by fixed atmospheric conditions, |
---|
574 | so these seas can be removed from the ocean domain. The user has the option |
---|
575 | to set the bathymetry in closed seas to zero (see \S\ref{MISC_closea}), but the |
---|
576 | code has to be adapted to the user's configuration. |
---|
577 | |
---|
578 | % ------------------------------------------------------------------------------------------------------------- |
---|
579 | % z-coordinate and reference coordinate transformation |
---|
580 | % ------------------------------------------------------------------------------------------------------------- |
---|
581 | \subsection[$z$-coordinate (\np{ln\_zco} or \key{zco})] |
---|
582 | {$z$-coordinate (\np{ln\_zco}=.true. or \key{zco}) and reference coordinate} |
---|
583 | \label{DOM_zco} |
---|
584 | |
---|
585 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
586 | \begin{figure}[!tb] \label{Fig_zgr} \begin{center} |
---|
587 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf} |
---|
588 | \caption{Default vertical mesh for ORCA2: 30 ocean levels (L30). Vertical level |
---|
589 | functions for (a) T-point depth and (b) the associated scale factor as computed |
---|
590 | from \eqref{DOM_zgr_ana} using \eqref{DOM_zgr_coef} in $z$-coordinate.} |
---|
591 | \end{center} \end{figure} |
---|
592 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
593 | |
---|
594 | The reference coordinate transformation $z_0 (k)$ defines the arrays $gdept_0$ |
---|
595 | and $gdepw_0$ for $T$- and $w$-points, respectively. As indicated on |
---|
596 | Fig.\ref{Fig_index_vert} \jp{jpk} is the number of $w$-levels. $gdepw_0(1)$ is the |
---|
597 | ocean surface. There are at most \jp{jpk}-1 $T$-points inside the ocean, the |
---|
598 | additional $T$-point at $jk=jpk$ is below the sea floor and is not used. |
---|
599 | The vertical location of $w$- and $T$-levels is defined from the analytic expression |
---|
600 | of the depth $z_0(k)$ whose analytical derivative with respect to $k$ provides the |
---|
601 | vertical scale factors. The user must provide the analytical expression of both |
---|
602 | $z_0$ and its first derivative with respect to $k$. This is done in routine \mdl{domzgr} |
---|
603 | through statement functions, using parameters provided in the \textit{par\_oce.h90} file. |
---|
604 | |
---|
605 | It is possible to define a simple regular vertical grid by giving zero stretching (\pp{ppacr=0}). |
---|
606 | In that case, the parameters \jp{jpk} (number of $w$-levels) and \pp{pphmax} |
---|
607 | (total ocean depth in meters) fully define the grid. |
---|
608 | |
---|
609 | For climate-related studies it is often desirable to concentrate the vertical resolution |
---|
610 | near the ocean surface. The following function is proposed as a standard for a |
---|
611 | $z$-coordinate (with either full or partial steps): |
---|
612 | \begin{equation} \label{DOM_zgr_ana} |
---|
613 | \begin{split} |
---|
614 | z_0 (k) &= h_{sur} -h_0 \;k-\;h_1 \;\log \left[ {\,\cosh \left( {{(k-h_{th} )} / {h_{cr} }} \right)\,} \right] \\ |
---|
615 | e_3^0 (k) &= \left| -h_0 -h_1 \;\tanh \left( {{(k-h_{th} )} / {h_{cr} }} \right) \right| |
---|
616 | \end{split} |
---|
617 | \end{equation} |
---|
618 | where $k=1$ to \jp{jpk} for $w$-levels and $k=1$ to $k=1$ for $T-$levels. Such an |
---|
619 | expression allows us to define a nearly uniform vertical location of levels at the |
---|
620 | ocean top and bottom with a smooth hyperbolic tangent transition in between |
---|
621 | (Fig.~\ref{Fig_zgr}). |
---|
622 | |
---|
623 | The most used vertical grid for ORCA2 has $10~m$ ($500~m)$ resolution in the |
---|
624 | surface (bottom) layers and a depth which varies from 0 at the sea surface to a |
---|
625 | minimum of $-5000~m$. This leads to the following conditions: |
---|
626 | \begin{equation} \label{DOM_zgr_coef} |
---|
627 | \begin{split} |
---|
628 | e_3 (1+1/2) &=10. \\ |
---|
629 | e_3 (jpk-1/2) &=500. \\ |
---|
630 | z(1) &=0. \\ |
---|
631 | z(jpk) &=-5000. \\ |
---|
632 | \end{split} |
---|
633 | \end{equation} |
---|
634 | |
---|
635 | With the choice of the stretching $h_{cr} =3$ and the number of levels |
---|
636 | \jp{jpk}=$31$, the four coefficients $h_{sur}$, $h_{0}$, $h_{1}$, and $h_{th}$ in |
---|
637 | \eqref{DOM_zgr_ana} have been determined such that \eqref{DOM_zgr_coef} is |
---|
638 | satisfied, through an optimisation procedure using a bisection method. For the first |
---|
639 | standard ORCA2 vertical grid this led to the following values: $h_{sur} =4762.96$, |
---|
640 | $h_0 =255.58, h_1 =245.5813$, and $h_{th} =21.43336$. The resulting depths and |
---|
641 | scale factors as a function of the model levels are shown in Fig.~\ref{Fig_zgr} and |
---|
642 | given in Table \ref{Tab_orca_zgr}. Those values correspond to the parameters |
---|
643 | \pp{ppsur}, \pp{ppa0}, \pp{ppa1}, \pp{ppkth} in the parameter file \mdl{par\_oce}. |
---|
644 | |
---|
645 | Rather than entering parameters $h_{sur}$, $h_{0}$, and $h_{1}$ directly, it is |
---|
646 | possible to recalculate them. In that case the user sets |
---|
647 | \pp{ppsur}=\pp{ppa0}=\pp{ppa1}=\pp{pp\_to\_be\_computed}, in \mdl{par\_oce}, |
---|
648 | and specifies instead the four following parameters: |
---|
649 | \begin{itemize} |
---|
650 | \item \pp{ppacr}=$h_{cr} $: stretching factor (nondimensional). The larger |
---|
651 | \pp{ppacr}, the smaller the stretching. Values from $3$ to $10$ are usual. |
---|
652 | \item \pp{ppkth}=$h_{th} $: is approximately the model level at which maximum |
---|
653 | stretching occurs (nondimensional, usually of order 1/2 or 2/3 of \jp{jpk}) |
---|
654 | \item \pp{ppdzmin}: minimum thickness for the top layer (in meters) |
---|
655 | \item \pp{pphmax}: total depth of the ocean (meters). |
---|
656 | \end{itemize} |
---|
657 | As an example, for the $45$ layers used in the DRAKKAR configuration those |
---|
658 | parameters are: \jp{jpk}=46, \pp{ppacr}=9, \pp{ppkth}=23.563, \pp{ppdzmin}=6m, |
---|
659 | \pp{pphmax}=5750m. |
---|
660 | |
---|
661 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
662 | \begin{table} \label{Tab_orca_zgr} |
---|
663 | \begin{center} \begin{tabular}{c||r|r|r|r} |
---|
664 | \hline |
---|
665 | \textbf{LEVEL}& \textbf{GDEPT}& \textbf{GDEPW}& \textbf{E3T }& \textbf{E3W } \\ \hline |
---|
666 | 1 & \textbf{ 5.00} & 0.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
667 | 2 & \textbf{15.00} & 10.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
668 | 3 & \textbf{25.00} & 20.00 & \textbf{ 10.00} & 10.00 \\ \hline |
---|
669 | 4 & \textbf{35.01} & 30.00 & \textbf{ 10.01} & 10.00 \\ \hline |
---|
670 | 5 & \textbf{45.01} & 40.01 & \textbf{ 10.01} & 10.01 \\ \hline |
---|
671 | 6 & \textbf{55.03} & 50.02 & \textbf{ 10.02} & 10.02 \\ \hline |
---|
672 | 7 & \textbf{65.06} & 60.04 & \textbf{ 10.04} & 10.03 \\ \hline |
---|
673 | 8 & \textbf{75.13} & 70.09 & \textbf{ 10.09} & 10.06 \\ \hline |
---|
674 | 9 & \textbf{85.25} & 80.18 & \textbf{ 10.17} & 10.12 \\ \hline |
---|
675 | 10 & \textbf{95.49} & 90.35 & \textbf{ 10.33} & 10.24 \\ \hline |
---|
676 | 11 & \textbf{105.97} & 100.69 & \textbf{ 10.65} & 10.47 \\ \hline |
---|
677 | 12 & \textbf{116.90} & 111.36 & \textbf{ 11.27} & 10.91 \\ \hline |
---|
678 | 13 & \textbf{128.70} & 122.65 & \textbf{ 12.47} & 11.77 \\ \hline |
---|
679 | 14 & \textbf{142.20} & 135.16 & \textbf{ 14.78} & 13.43 \\ \hline |
---|
680 | 15 & \textbf{158.96} & 150.03 & \textbf{ 19.23} & 16.65 \\ \hline |
---|
681 | 16 & \textbf{181.96} & 169.42 & \textbf{ 27.66} & 22.78 \\ \hline |
---|
682 | 17 & \textbf{216.65} & 197.37 & \textbf{ 43.26} & 34.30 \\ \hline |
---|
683 | 18 & \textbf{272.48} & 241.13 & \textbf{ 70.88} & 55.21 \\ \hline |
---|
684 | 19 & \textbf{364.30} & 312.74 & \textbf{116.11} & 90.99 \\ \hline |
---|
685 | 20 & \textbf{511.53} & 429.72 & \textbf{181.55} & 146.43 \\ \hline |
---|
686 | 21 & \textbf{732.20} & 611.89 & \textbf{261.03} & 220.35 \\ \hline |
---|
687 | 22 & \textbf{1033.22}& 872.87 & \textbf{339.39} & 301.42 \\ \hline |
---|
688 | 23 & \textbf{1405.70}& 1211.59 & \textbf{402.26} & 373.31 \\ \hline |
---|
689 | 24 & \textbf{1830.89}& 1612.98 & \textbf{444.87} & 426.00 \\ \hline |
---|
690 | 25 & \textbf{2289.77}& 2057.13 & \textbf{470.55} & 459.47 \\ \hline |
---|
691 | 26 & \textbf{2768.24}& 2527.22 & \textbf{484.95} & 478.83 \\ \hline |
---|
692 | 27 & \textbf{3257.48}& 3011.90 & \textbf{492.70} & 489.44 \\ \hline |
---|
693 | 28 & \textbf{3752.44}& 3504.46 & \textbf{496.78} & 495.07 \\ \hline |
---|
694 | 29 & \textbf{4250.40}& 4001.16 & \textbf{498.90} & 498.02 \\ \hline |
---|
695 | 30 & \textbf{4749.91}& 4500.02 & \textbf{500.00} & 499.54 \\ \hline |
---|
696 | 31 & \textbf{5250.23}& 5000.00 & \textbf{500.56} & 500.33 \\ \hline |
---|
697 | \end{tabular} \end{center} |
---|
698 | \caption{Default vertical mesh in $z$-coordinate for 30 layers ORCA2 configuration |
---|
699 | as computed from \eqref{DOM_zgr_ana} using the coefficients given in |
---|
700 | \eqref{DOM_zgr_coef}} |
---|
701 | \end{table} |
---|
702 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
703 | |
---|
704 | % ------------------------------------------------------------------------------------------------------------- |
---|
705 | % z-coordinate with partial step |
---|
706 | % ------------------------------------------------------------------------------------------------------------- |
---|
707 | \subsection [$z$-coordinate with partial step (\np{ln\_zps})] |
---|
708 | {$z$-coordinate with partial step (\np{ln\_zps}=.true.)} |
---|
709 | \label{DOM_zps} |
---|
710 | %--------------------------------------------namdom------------------------------------------------------- |
---|
711 | \namdisplay{namdom} |
---|
712 | %-------------------------------------------------------------------------------------------------------------- |
---|
713 | |
---|
714 | In $z$-coordinate partial step, the depths of the model levels are defined by the |
---|
715 | reference analytical function $z_0 (k)$ as described in the previous |
---|
716 | section, \emph{except} in the bottom layer. The thickness of the bottom layer is |
---|
717 | allowed to vary as a function of geographical location $(\lambda,\varphi)$ to allow a |
---|
718 | better representation of the bathymetry, especially in the case of small |
---|
719 | slopes (where the bathymetry varies by less than one level thickness from |
---|
720 | one grid point to the next). The reference layer thicknesses $e_{3t}^0$ have been |
---|
721 | defined in the absence of bathymetry. With partial steps, layers from 1 to |
---|
722 | \jp{jpk}-2 can have a thickness smaller than $e_{3t}(jk)$. The model deepest layer (\jp{jpk}-1) |
---|
723 | is allowed to have either a smaller or larger thickness than $e_{3t}(jpk)$: the |
---|
724 | maximum thickness allowed is $2*e_{3t}(jpk-1)$. This has to be kept in mind when |
---|
725 | specifying the maximum depth \pp{pphmax} in partial steps: for example, with |
---|
726 | \pp{pphmax}$=5750~m$ for the DRAKKAR 45 layer grid, the maximum ocean depth |
---|
727 | allowed is actually $6000~m$ (the default thickness $e_{3t}(jpk-1)$ being $250~m$). |
---|
728 | Two variables in the namdom namelist are used to define the partial step |
---|
729 | vertical grid. The mimimum water thickness (in meters) allowed for a cell |
---|
730 | partially filled with bathymetry at level jk is the minimum of \np{e3zpsmin} |
---|
731 | (thickness in meters, usually $20~m$) or $e_{3t}(jk)*\np{e3zps\_rat}$ (a fraction, |
---|
732 | usually 10\%, of the default thickness $e_{3t}(jk)$). |
---|
733 | |
---|
734 | \colorbox{yellow}{Add a figure here of pstep especially at last ocean level } |
---|
735 | |
---|
736 | % ------------------------------------------------------------------------------------------------------------- |
---|
737 | % s-coordinate |
---|
738 | % ------------------------------------------------------------------------------------------------------------- |
---|
739 | \subsection [$s$-coordinate (\np{ln\_sco})] |
---|
740 | {$s$-coordinate (\np{ln\_sco}=true)} |
---|
741 | \label{DOM_sco} |
---|
742 | %------------------------------------------nam_zgr_sco--------------------------------------------------- |
---|
743 | \namdisplay{nam_zgr_sco} |
---|
744 | %-------------------------------------------------------------------------------------------------------------- |
---|
745 | In $s$-coordinate (\key{sco} is defined), the depth and thickness of the model |
---|
746 | levels are defined from the product of a depth field and either a stretching |
---|
747 | function or its derivative, respectively: |
---|
748 | \begin{equation} \label{DOM_sco_ana} |
---|
749 | \begin{split} |
---|
750 | z(k) &= h(i,j) \; z_0(k) \\ |
---|
751 | e_3(k) &= h(i,j) \; z_0'(k) |
---|
752 | \end{split} |
---|
753 | \end{equation} |
---|
754 | where $h$ is the depth of the last $w$-level ($z_0(k)$) defined at the $T$-point |
---|
755 | location in the horizontal and $z_0(k)$ is a function which varies from $0$ at the sea |
---|
756 | surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean |
---|
757 | depth, since a mixed step-like and bottom-following representation of the |
---|
758 | topography can be used (Fig.~\ref{Fig_z_zps_s_sps}d-e). In the example provided |
---|
759 | (\hf{zgr\_s} file) $h$ is a smooth envelope bathymetry and steps are used to represent |
---|
760 | sharp bathymetric gradients. |
---|
761 | |
---|
762 | A new flexible stretching function, modified from \citet{Song1994} is provided as an example: |
---|
763 | \begin{equation} \label{DOM_sco_function} |
---|
764 | \begin{split} |
---|
765 | z &= h_c +( h-h_c)\;c s) \\ |
---|
766 | c(s) &= \frac{ \left[ \tanh{ \left( \theta \, (s+b) \right)} |
---|
767 | - \tanh{ \left( \theta \, b \right)} \right]} |
---|
768 | {2\;\sinh \left( \theta \right)} |
---|
769 | \end{split} |
---|
770 | \end{equation} |
---|
771 | where $h_c$ is the thermocline depth and $\theta$ and $b$ are the surface and |
---|
772 | bottom control parameters such that $0\leqslant \theta \leqslant 20$, and |
---|
773 | $0\leqslant b\leqslant 1$. $b$ has been designed to allow surface and/or bottom |
---|
774 | increase of the vertical resolution (Fig.~\ref{Fig_sco_function}). |
---|
775 | |
---|
776 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
777 | \begin{figure}[!tb] \label{Fig_sco_function} \begin{center} |
---|
778 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf} |
---|
779 | \caption{Examples of the stretching function applied to a sea mont; from left to right: |
---|
780 | surface, surface and bottom, and bottom intensified resolutions} |
---|
781 | \end{center} \end{figure} |
---|
782 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
783 | |
---|
784 | % ------------------------------------------------------------------------------------------------------------- |
---|
785 | % z*- or s*-coordinate |
---|
786 | % ------------------------------------------------------------------------------------------------------------- |
---|
787 | \subsection{$z^*$- or $s^*$-coordinate (add \key{vvl}) } |
---|
788 | \label{DOM_zgr_vvl} |
---|
789 | |
---|
790 | This option is described in the Report by Levier \textit{et al.} (2007), available on |
---|
791 | the \NEMO web site. |
---|
792 | |
---|
793 | %gm% key advantage: minimise the diffusion/dispertion associated with advection in response to high frequency surface disturbances |
---|
794 | |
---|
795 | % ------------------------------------------------------------------------------------------------------------- |
---|
796 | % level bathymetry and mask |
---|
797 | % ------------------------------------------------------------------------------------------------------------- |
---|
798 | \subsection{level bathymetry and mask} |
---|
799 | \label{DOM_msk} |
---|
800 | |
---|
801 | Whatever the vertical coordinate used, the model offers the possibility of |
---|
802 | representing the bottom topography with steps that follow the face of the |
---|
803 | model cells (step like topography) \citep{Madec1996}. The distribution of |
---|
804 | the steps in the horizontal is defined in a 2D integer array, mbathy, which |
---|
805 | gives the number of ocean levels ($i.e.$ those that are not masked) at each |
---|
806 | $T$-point. mbathy is computed from the meter bathymetry using the definiton of |
---|
807 | gdept as the number of $T$-points which gdept $\leq$ bathy. Note that in version |
---|
808 | NEMO v2.3, the user still has to provide the "level" bathymetry in a NetCDF |
---|
809 | file when using the full step option (\np{ln\_zco}), rather than the bathymetry |
---|
810 | in meters: both will be allowed in future versions. |
---|
811 | |
---|
812 | Modifications of the model bathymetry are performed in the \textit{bat\_ctl} |
---|
813 | routine (see \mdl{domzgr} module) after mbathy is computed. Isolated grid points |
---|
814 | that do not communicate with another ocean point at the same level are eliminated. |
---|
815 | |
---|
816 | In the case of the rigid-lid approximation when islands occur in the computational |
---|
817 | domain (\np{ln\_dynspg\_rl}=.true. and \key{island} is defined), the \textit{mbathy} |
---|
818 | array must be provided and takes values from $-N$ to \jp{jpk}-1. It provides the |
---|
819 | following information: $mbathy(i,j) = -n, \ n \in \left] 0,N \right]$, $T$-points are |
---|
820 | land points on the $n^{th}$ island ; $mbathy(i,j) =0$, $T$-points are land points |
---|
821 | on the main land (continent) ; $mbathy(i,j) =k$, the first $k$ $T$- and $w$-points |
---|
822 | are ocean points, the others are points below the ocean floor. |
---|
823 | |
---|
824 | This is used to compute the island barotropic stream function used in the rigid lid |
---|
825 | computation (see \S\ref{MISC_solisl}). |
---|
826 | |
---|
827 | From the \textit{mbathy} array, the mask fields are defined as follows: |
---|
828 | \begin{align*} |
---|
829 | tmask(i,j,k) &= \begin{cases} \; 1& \text{ if $k\leq mbathy(i,j)$ } \\ |
---|
830 | \; 0& \text{ if $k\leq mbathy(i,j)$ } \end{cases} \\ |
---|
831 | umask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ |
---|
832 | vmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i,j+1,k) \\ |
---|
833 | fmask(i,j,k) &= \; tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ |
---|
834 | & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) |
---|
835 | \end{align*} |
---|
836 | |
---|
837 | \gmcomment{ STEVEN: are the dots multiplications?} |
---|
838 | |
---|
839 | Note that \textit{wmask} is not defined as it is exactly equal to \textit{tmask} with |
---|
840 | the numerical indexing used (\S~\ref{DOM_Num_Index}). Moreover, the |
---|
841 | specification of closed lateral boundaries requires that at least the first and last |
---|
842 | rows and columns of the \textit{mbathy} array are set to zero. In the particular |
---|
843 | case of an east-west cyclical boundary condition, \textit{mbathy} has its last |
---|
844 | column equal to the second one and its first column equal to the last but one |
---|
845 | (and so too the mask arrays) (see \S~\ref{LBC_jperio}). |
---|
846 | |
---|
847 | %%% |
---|
848 | \gmcomment{ \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}. } |
---|
849 | %%% |
---|
850 | |
---|
851 | % ================================================================ |
---|
852 | % Time Discretisation |
---|
853 | % ================================================================ |
---|
854 | \section{Time Discretisation} |
---|
855 | \label{DOM_nxt} |
---|
856 | |
---|
857 | The time stepping used in \NEMO is a three level scheme that can be |
---|
858 | represented as follows: |
---|
859 | \begin{equation} \label{Eq_DOM_nxt} |
---|
860 | x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ \text{RHS}_x^{t-\Delta t,t,t+\Delta t} |
---|
861 | \end{equation} |
---|
862 | where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the |
---|
863 | corresponding time evolution equation; $\Delta t$ is the time step; and the |
---|
864 | superscripts indicate the time at which a quantity is evaluated. Each term of the |
---|
865 | RHS is evaluated at a specific time step(s) depending on the physics with which |
---|
866 | it is associated. |
---|
867 | |
---|
868 | The choice of the time step used for this evaluation is discussed below as |
---|
869 | well as the implications in terms of starting or restarting a model |
---|
870 | simulation. Note that the time stepping is generally performed in a one step |
---|
871 | operation. With such a complex and nonlinear system of equations it would be |
---|
872 | dangerous to let a prognostic variable evolve in time for each term separately. |
---|
873 | |
---|
874 | The three level scheme requires three arrays for each prognostic variables. |
---|
875 | For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array, |
---|
876 | although referred to as $x_a$ (after) in the code, is usually not the variable at |
---|
877 | the next time step; but rather it is used to store the time derivative (RHS in |
---|
878 | \eqref{Eq_DOM_nxt}) prior to time-stepping the equation. Generally, the time |
---|
879 | stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt} |
---|
880 | modules, except for implicit vertical diffusion or sea surface height when |
---|
881 | time-splitting options are used. |
---|
882 | |
---|
883 | % ------------------------------------------------------------------------------------------------------------- |
---|
884 | % Non-Diffusive Part---Leapfrog Scheme |
---|
885 | % ------------------------------------------------------------------------------------------------------------- |
---|
886 | \subsection{Non-Diffusive Part --- Leapfrog Scheme} |
---|
887 | \label{DOM_nxt_leap_frog} |
---|
888 | |
---|
889 | The time stepping used for non-diffusive processes is the well-known |
---|
890 | leapfrog scheme. It is a time centred scheme, i.e. the RHS is evaluated at |
---|
891 | time step $t$, the now time step. It is only used for non-diffusive terms, |
---|
892 | that is momentum and tracer advection, pressure gradient, and Coriolis |
---|
893 | terms. This scheme is widely used for advective processes in low-viscosity |
---|
894 | fluids. It is an efficient method that achieves second-order accuracy with |
---|
895 | just one right hand side evaluation per time step. Moreover, it does not |
---|
896 | artificially damp linear oscillatory motion nor does it produce instability |
---|
897 | by amplifying the oscillations. These advantages are somewhat diminished by |
---|
898 | the large phase-speed error of the leapfrog scheme, and the unsuitability of |
---|
899 | leapfrog differencing for the representation of diffusive and Rayleigh |
---|
900 | damping processes. However, the most serious problem associated with the |
---|
901 | leapfrog scheme is a high-frequency computational noise called |
---|
902 | "time-splitting" \citep{Haltiner1980} that develops when the method |
---|
903 | is used to model non linear fluid dynamics: the even and odd time steps tend |
---|
904 | to diverge into a physical and a computational mode. Time splitting can |
---|
905 | be controlled through the use of an Asselin time filter (first designed by |
---|
906 | \citep{Robert1966} and more comprehensively studied by \citet{Asselin1972}), |
---|
907 | or by periodically reinitialising the leapfrog solution through a single |
---|
908 | integration step with a two-level scheme. In \NEMO we follow the first |
---|
909 | strategy: |
---|
910 | \begin{equation} \label{Eq_DOM_nxt_asselin} |
---|
911 | x_F^t = x^t + \gamma \, \left[ x_f^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] |
---|
912 | \end{equation} |
---|
913 | where the subscript $f$ denotes filtered values and $\gamma$ is the Asselin |
---|
914 | coefficient. $\gamma$ is initialized as \np{atfp} (namelist parameter). |
---|
915 | Its default value is \np{atfp}=0.1. This default value causes a significant dissipation |
---|
916 | of high frequency motions. Recommended values in idealized studies of shallow |
---|
917 | water turbulence are two orders of magnitude smaller (\citep{Farge1987}). |
---|
918 | Both strategies do, nevertheless, degrade the accuracy of the calculation from |
---|
919 | second to first order. The leapfrog scheme combined with a Robert-Asselin |
---|
920 | time filter has been preferred to other time differencing schemes such as |
---|
921 | predictor corrector or trapezoidal schemes, because the user has an explicit |
---|
922 | and simple control of the magnitude of the time diffusion of the scheme. |
---|
923 | In association with the 2nd order centred space discretisation of the |
---|
924 | advective terms in the momentum and tracer equations, it avoids implicit |
---|
925 | numerical diffusion in both the time and space discretisations of the |
---|
926 | advective term: they are both set explicitly by the user through the Robert-Asselin |
---|
927 | filter parameter and the viscous and diffusive coefficients. |
---|
928 | |
---|
929 | \gmcomment{ |
---|
930 | %gm - reflexion about leapfrog: ongoing work with Matthieu Leclair |
---|
931 | % to be updated latter with addition of new time stepping strategy |
---|
932 | \colorbox{yellow}{Note}: |
---|
933 | The Robert-Asselin time filter slightly departs from a simple second order time |
---|
934 | diffusive operator computed with a forward time stepping due to the presence of |
---|
935 | $x_f^{t-\Delta t}$ in the right hand side of \ref{Eq_DOM_nxt_asselin}. The original |
---|
936 | willing of Robert1966 and Asselin1972 was to design a time filter that allow much |
---|
937 | larger parameter than 0.5. is due to computer saving consideration. In the original |
---|
938 | asselin filter, $x^{t-\Delta t}$ is used instead: |
---|
939 | \begin{equation} \label{Eq_DOM_nxt_asselin_true} |
---|
940 | x_f^t = x^t + \gamma \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] |
---|
941 | \end{equation} |
---|
942 | Applying a "true" Asselin time filter is nothing more than adding a harmonic |
---|
943 | diffusive operator in time. Indeed, equations \ref{Eq_DOM_nxt} and |
---|
944 | \ref{Eq_DOM_nxt_asselin_true} can be rewritten together as: |
---|
945 | \begin{equation} \label{Eq_DOM_nxt2} |
---|
946 | \begin{split} |
---|
947 | \frac{ x^{t+\Delta t} - x^{t-\Delta t} } { 2 \,\Delta t } |
---|
948 | &= \text{RHS}_x^{t-\Delta t,t,t+\Delta t} + \frac{ x_f^t - x^t }{2 \,\Delta t} \\ |
---|
949 | &= \text{RHS}_x^{t-\Delta t,t,t+\Delta t} |
---|
950 | + \gamma\ \frac{ \, \left[ x^{t-\Delta t} - 2 x^t + x^{t+\Delta t} \right] }{2 \,\Delta t} \\ |
---|
951 | &= \text{RHS}_x^{t-\Delta t,t,t+\Delta t} |
---|
952 | + 2 \Delta t \ \gamma \ \frac{1}{{2 \Delta t}^2} |
---|
953 | \,\delta_{t-1}\,\left[ \delta_{t+1/2}\left[ x^t \right] \right] |
---|
954 | \end{split} |
---|
955 | \end{equation} |
---|
956 | expressing this again in a continuous form, one finds that the Asselin filter leads to : |
---|
957 | \begin{equation} \label{Eq_DOM_nxt3} |
---|
958 | \frac{ \partial x} { \partial t } = \text{RHS} + 2 \,\Delta t \ \gamma \ \frac{ {\partial}^2 x}{ \partial t ^2 } |
---|
959 | \end{equation} |
---|
960 | |
---|
961 | Equations \ref{Eq_DOM_nxt2} and \ref{Eq_DOM_nxt3} suggest several remarks. |
---|
962 | First the Asselin filter is definitively a second order time diffusive operator which is |
---|
963 | evaluated at centered time step. The magnitude of this diffusion is proportional to |
---|
964 | the time step (with $\gamma$ usually taken between $10^{-1}$ to $10^{-3}$). |
---|
965 | Second, this term has to be taken into account in all budgets of the equations |
---|
966 | (mass, heat content, salt content, kinetic energy...). Nevertheless,we stress here |
---|
967 | that it is small and does not create systematic biases. Indeed let us evaluate how |
---|
968 | it contributes to the time evolution of $x$ between $t_o$ and $t_1$: |
---|
969 | \begin{equation} \label{Eq_DOM_nxt4} |
---|
970 | \begin{split} |
---|
971 | t_1-t_o &= \int_{t_o}^{t_1} \frac{ \partial x} { \partial t } dt \\ |
---|
972 | &= \int_{t_o}^{t_1} \text{RHS} dt + 2 \,\Delta t \ \gamma \left( |
---|
973 | \left. \frac{ \partial x}{ \partial t } \right|_1 |
---|
974 | - \left. \frac{ \partial x}{ \partial t } \right|_o \right) |
---|
975 | \end{split} |
---|
976 | \end{equation} |
---|
977 | } |
---|
978 | |
---|
979 | Alternative time stepping schemes are currently under investigation. |
---|
980 | |
---|
981 | % ------------------------------------------------------------------------------------------------------------- |
---|
982 | % Diffusive Part---Forward or Backward Scheme |
---|
983 | % ------------------------------------------------------------------------------------------------------------- |
---|
984 | \subsection{Diffusive Part --- Forward or Backward Scheme} |
---|
985 | \label{DOM_nxt_forward_imp} |
---|
986 | |
---|
987 | The leapfrog differencing scheme is unsuitable for the representation of |
---|
988 | diffusive and damping processes. For a tendancy $D_x$, representing a |
---|
989 | diffusive term or a restoring term to a tracer climatology |
---|
990 | (when present, see \S~\ref{TRA_dmp}), a forward time differencing scheme |
---|
991 | is used : |
---|
992 | \begin{equation} \label{Eq_DOM_nxt_euler} |
---|
993 | x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ {D_x}^{t-\Delta t} |
---|
994 | \end{equation} |
---|
995 | |
---|
996 | This is diffusive in time and conditionally stable. For example, the |
---|
997 | conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}: |
---|
998 | \begin{equation} \label{Eq_DOM_nxt_euler_stability} |
---|
999 | A^h < \left\{ |
---|
1000 | \begin{aligned} |
---|
1001 | &\frac{e^2}{ 8 \; \Delta t } &&\quad \text{laplacian diffusion} \\ |
---|
1002 | &\frac{e^4}{64 \; \Delta t } &&\quad \text{bilaplacian diffusion} |
---|
1003 | \end{aligned} |
---|
1004 | \right. |
---|
1005 | \end{equation} |
---|
1006 | where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is |
---|
1007 | the mixing coefficient. The linear constraint \eqref{Eq_DOM_nxt_euler_stability} |
---|
1008 | is a necessary condition, but not sufficient. If it is not satisfied, even mildly, |
---|
1009 | then the model soon becomes wildly unstable. The instability can be removed |
---|
1010 | by either reducing the length of the time steps or reducing the mixing coefficient. |
---|
1011 | |
---|
1012 | For the vertical diffusion terms, a forward time differencing scheme can be |
---|
1013 | used, but usually the numerical stability condition implies a strong |
---|
1014 | constraint on the time step. Two solutions are available in \NEMO to overcome |
---|
1015 | the stability constraint: $(a)$ a forward time differencing scheme using a |
---|
1016 | time splitting technique (\np{ln\_zdfexp}=.true.) or $(b)$ a backward (or implicit) |
---|
1017 | time differencing scheme by \np{ln\_zdfexp}=.false.). In $(a)$, the master |
---|
1018 | time step $\Delta $t is cut into $N$ fractional time steps so that the |
---|
1019 | stability criterion is reduced by a factor of $N$. The computation is done as |
---|
1020 | follows: |
---|
1021 | \begin{equation} \label{Eq_DOM_nxt_ts} |
---|
1022 | \begin{split} |
---|
1023 | & u_\ast ^{t-\Delta t} = u^{t-\Delta t} \\ |
---|
1024 | & u_\ast ^{t-\Delta t+L\frac{2\Delta t}{N}}=u_\ast ^{t-\Delta t+\left( {L-1} |
---|
1025 | \right)\frac{2\Delta t}{N}}+\frac{2\Delta t}{N}\;\text{DF}^{t-\Delta t+\left( {L-1} \right)\frac{2\Delta t}{N}} |
---|
1026 | \quad \text{for $L=1$ to $N$} \\ |
---|
1027 | & u^{t+\Delta t} = u_\ast^{t+\Delta t} |
---|
1028 | \end{split} |
---|
1029 | \end{equation} |
---|
1030 | with DF a vertical diffusion term. The number of fractional time steps, $N$, is given |
---|
1031 | by setting \np{n\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally |
---|
1032 | stable but diffusive. It can be written as follows: |
---|
1033 | \begin{equation} \label{Eq_DOM_nxt_imp} |
---|
1034 | x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ \text{RHS}_x^{t+\Delta t} |
---|
1035 | \end{equation} |
---|
1036 | |
---|
1037 | This scheme is rather time consuming since it requires a matrix inversion, |
---|
1038 | but it becomes attractive since a splitting factor of 3 or more is needed |
---|
1039 | for the forward time differencing scheme. For example, the finite difference |
---|
1040 | approximation of the temperature equation is: |
---|
1041 | \begin{equation} \label{Eq_DOM_nxt_imp_zdf} |
---|
1042 | \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\Delta t}\equiv \text{RHS}+\frac{1}{e_{3T} }\delta |
---|
1043 | _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} |
---|
1044 | \right] |
---|
1045 | \end{equation} |
---|
1046 | where RHS is the right hand side of the equation except for the vertical diffusion term. |
---|
1047 | We rewrite \eqref{Eq_DOM_nxt_imp} as: |
---|
1048 | \begin{equation} \label{Eq_DOM_nxt_imp_mat} |
---|
1049 | -c(k+1)\;u^{t+1}(k+1)+d(k)\;u^{t+1}(k)-\;c(k)\;u^{t+1}(k-1) \equiv b(k) |
---|
1050 | \end{equation} |
---|
1051 | where |
---|
1052 | \begin{align*} |
---|
1053 | c(k) &= A_w^{vm} (k) \, / \, e_{3uw} (k) \\ |
---|
1054 | d(k) &= e_{3u} (k) \, / \, (2\Delta t) + c_k + c_{k+1} \\ |
---|
1055 | b(k) &= e_{3u} (k) \; \left( u^{t-1}(k) \, / \, (2\Delta t) + \text{RHS} \right) |
---|
1056 | \end{align*} |
---|
1057 | |
---|
1058 | \eqref{Eq_DOM_nxt_imp_mat} is a linear system of equations which associated |
---|
1059 | matrix is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal |
---|
1060 | term is greater than the sum of the two extra-diagonal terms, therefore a special |
---|
1061 | adaptation of the Gauss elimination procedure is used to find the solution |
---|
1062 | (see for example \citet{Richtmyer1967}). |
---|
1063 | |
---|
1064 | % ------------------------------------------------------------------------------------------------------------- |
---|
1065 | % Start/Restart strategy |
---|
1066 | % ------------------------------------------------------------------------------------------------------------- |
---|
1067 | \subsection{Start/Restart strategy} |
---|
1068 | \label{DOM_nxt_rst} |
---|
1069 | %--------------------------------------------namrun------------------------------------------- |
---|
1070 | \namdisplay{namrun} |
---|
1071 | %-------------------------------------------------------------------------------------------------------------- |
---|
1072 | |
---|
1073 | The first time step of this three level scheme when starting from initial conditions |
---|
1074 | is a forward step (Euler time integration): |
---|
1075 | \begin{equation} \label{Eq_DOM_euler} |
---|
1076 | x^1 = x^0 + \Delta t \ \text{RHS}^0 |
---|
1077 | \end{equation} |
---|
1078 | |
---|
1079 | It is also possible to restart from a previous computation, by using a |
---|
1080 | restart file. The restart strategy is designed to ensure perfect |
---|
1081 | restartability of the code: the user should obtain the same results to |
---|
1082 | machine precision either by running the model for $2N$ time steps in one go, |
---|
1083 | or by performing two consecutive experiments of $N$ steps with a restart. |
---|
1084 | This requires saving two time levels and many auxiliary data in the restart |
---|
1085 | files in machine precision. |
---|
1086 | |
---|
1087 | Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure |
---|
1088 | gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be |
---|
1089 | added in the restart file to ensure an exact restartability. This is done only optionally |
---|
1090 | via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of |
---|
1091 | restart file can be obtained when the restartability is not a key issue (operational |
---|
1092 | oceanography or ensemble simulation for seasonal forcast). |
---|
1093 | %%% |
---|
1094 | \gmcomment{add here how to force the restart to contain only one time step for operational purposes} |
---|
1095 | %%% |
---|
1096 | |
---|
1097 | \gmcomment{ % add a subsection here |
---|
1098 | |
---|
1099 | %------------------------------------------------------------------------------------------------------------- |
---|
1100 | % Time Domain |
---|
1101 | % ------------------------------------------------------------------------------------------------------------- |
---|
1102 | \subsection{Time domain} |
---|
1103 | \label{DOM_nxt_time} |
---|
1104 | |
---|
1105 | \colorbox{yellow}{add here a few word on nit000 and nitend} |
---|
1106 | |
---|
1107 | \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} |
---|
1108 | |
---|
1109 | } %% end add |
---|
1110 | |
---|