[707] | 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: |
---|
[817] | 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) |
---|
[707] | 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 | |
---|
[996] | 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} |
---|
[817] | 20 | %%%% |
---|
[707] | 21 | |
---|
| 22 | |
---|
[817] | 23 | \newpage |
---|
| 24 | $\ $\newline % force a new ligne |
---|
[707] | 25 | |
---|
| 26 | |
---|
[817] | 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. |
---|
[707] | 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} |
---|
[998] | 47 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_cell.pdf} |
---|
[817] | 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} |
---|
[707] | 52 | \end{center} \end{figure} |
---|
| 53 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 54 | |
---|
[817] | 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. |
---|
[707] | 66 | |
---|
[817] | 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 |
---|
[1224] | 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 |
---|
[817] | 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). |
---|
[707] | 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} |
---|
[817] | 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})} |
---|
[707] | 102 | \end{center} |
---|
| 103 | \end{table} |
---|
| 104 | |
---|
| 105 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 106 | % Vector Invariant Formulation |
---|
| 107 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 108 | \subsection{Discrete Operators} |
---|
| 109 | \label{DOM_operators} |
---|
| 110 | |
---|
[817] | 111 | Given the values of a variable $q$ at adjacent points, the differencing and |
---|
| 112 | averaging operators at the midpoint between them are: |
---|
[707] | 113 | \begin{subequations} \label{Eq_di_mi} |
---|
| 114 | \begin{align} |
---|
[817] | 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 |
---|
[707] | 117 | \end{align} |
---|
| 118 | \end{subequations} |
---|
| 119 | |
---|
[817] | 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: |
---|
[707] | 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 | |
---|
[817] | 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: |
---|
[707] | 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 | |
---|
[817] | 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*} |
---|
[707] | 170 | |
---|
[817] | 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): |
---|
[707] | 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} |
---|
[817] | 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$). |
---|
[707] | 181 | |
---|
[817] | 182 | In continuous form, the following properties are satisfied: |
---|
[707] | 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 | |
---|
[817] | 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)$. |
---|
[707] | 193 | |
---|
[817] | 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} |
---|
[707] | 208 | |
---|
[817] | 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 |
---|
[707] | 213 | demonstrate integral conservative properties of the discrete formulation chosen. |
---|
| 214 | |
---|
| 215 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 216 | % Numerical Indexing |
---|
[707] | 217 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 218 | \subsection{Numerical Indexing} |
---|
[707] | 219 | \label{DOM_Num_Index} |
---|
| 220 | |
---|
| 221 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 222 | \begin{figure}[!tb] \label{Fig_index_hor} \begin{center} |
---|
[998] | 223 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_index_hor.pdf} |
---|
[817] | 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} |
---|
[707] | 227 | \end{center} \end{figure} |
---|
| 228 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 229 | |
---|
[817] | 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$. |
---|
[707] | 237 | |
---|
| 238 | % ----------------------------------- |
---|
[817] | 239 | % Horizontal Indexing |
---|
[707] | 240 | % ----------------------------------- |
---|
[817] | 241 | \subsubsection{Horizontal Indexing} |
---|
[707] | 242 | \label{DOM_Num_Index_hor} |
---|
| 243 | |
---|
[817] | 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. |
---|
[707] | 248 | |
---|
| 249 | % ----------------------------------- |
---|
[817] | 250 | % Vertical indexing |
---|
[707] | 251 | % ----------------------------------- |
---|
[817] | 252 | \subsubsection{Vertical Indexing} |
---|
[707] | 253 | \label{DOM_Num_Index_vertical} |
---|
| 254 | |
---|
[817] | 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 |
---|
[1224] | 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. |
---|
[707] | 270 | |
---|
| 271 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 272 | \begin{figure}[!pt] \label{Fig_index_vert} \begin{center} |
---|
[998] | 273 | \includegraphics[width=.90\textwidth]{./TexFiles/Figures/Fig_index_vert.pdf} |
---|
[817] | 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.} |
---|
[707] | 277 | \end{center} \end{figure} |
---|
| 278 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 279 | |
---|
| 280 | % ----------------------------------- |
---|
[817] | 281 | % Domain Size |
---|
[707] | 282 | % ----------------------------------- |
---|
[817] | 283 | \subsubsection{Domain Size} |
---|
[707] | 284 | \label{DOM_size} |
---|
| 285 | |
---|
[817] | 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}). |
---|
[707] | 310 | |
---|
| 311 | % ================================================================ |
---|
| 312 | % Domain: Horizontal Grid (mesh) |
---|
| 313 | % ================================================================ |
---|
[817] | 314 | \section [Domain: Horizontal Grid (mesh) (\textit{domhgr})] |
---|
| 315 | {Domain: Horizontal Grid (mesh) \small{(\mdl{domhgr} module)} } |
---|
[707] | 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 | |
---|
[817] | 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. |
---|
[707] | 332 | |
---|
[817] | 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: |
---|
[707] | 342 | \begin{flalign*} |
---|
[817] | 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) |
---|
[707] | 347 | \end{flalign*} |
---|
| 348 | \begin{flalign*} |
---|
| 349 | e_{1T} &\equiv \text{e1t} = r_a |\lambda'(i) \; \cos\varphi(j) |& |
---|
[817] | 350 | e_{2T} &\equiv \text{e2t} = r_a |\varphi'(j)| \\ |
---|
[707] | 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*} |
---|
[817] | 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. |
---|
[707] | 363 | |
---|
[817] | 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}. |
---|
[707] | 376 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 377 | \begin{figure}[!t] \label{Fig_zgr_e3} \begin{center} |
---|
[998] | 378 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr_e3.pdf} |
---|
[817] | 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$. } |
---|
[707] | 386 | \end{center} \end{figure} |
---|
| 387 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 388 | |
---|
| 389 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 390 | % Choice of horizontal grid |
---|
[707] | 391 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 392 | \subsection{Choice of horizontal grid} |
---|
[707] | 393 | \label{DOM_hgr_msh_choice} |
---|
| 394 | |
---|
[817] | 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} |
---|
[707] | 404 | |
---|
[817] | 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. |
---|
[707] | 426 | |
---|
[817] | 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). |
---|
[707] | 442 | |
---|
[817] | 443 | The choice of the grid must be consistent with the boundary conditions specified |
---|
| 444 | by the parameter \jp{jperio} (see {\S\ref{LBC}). |
---|
[707] | 445 | |
---|
| 446 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 447 | % Grid files |
---|
| 448 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 449 | \subsection{Grid files} |
---|
| 450 | \label{DOM_hgr_files} |
---|
| 451 | |
---|
[817] | 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. |
---|
[707] | 461 | |
---|
| 462 | % ================================================================ |
---|
| 463 | % Domain: Vertical Grid (domzgr) |
---|
| 464 | % ================================================================ |
---|
[817] | 465 | \section [Domain: Vertical Grid (\textit{domzgr})] |
---|
| 466 | {Domain: Vertical Grid \small{(\mdl{domzgr} module)} } |
---|
[707] | 467 | \label{DOM_zgr} |
---|
| 468 | %-----------------------------------------nam_zgr & namdom------------------------------------------- |
---|
| 469 | \namdisplay{nam_zgr} |
---|
| 470 | \namdisplay{namdom} |
---|
| 471 | %------------------------------------------------------------------------------------------------------------- |
---|
| 472 | |
---|
[817] | 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. |
---|
[707] | 480 | |
---|
| 481 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 482 | \begin{figure}[!tb] \label{Fig_z_zps_s_sps} \begin{center} |
---|
[998] | 483 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_z_zps_s_sps.pdf} |
---|
[817] | 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).} |
---|
[707] | 493 | \end{center} \end{figure} |
---|
| 494 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 495 | |
---|
[817] | 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). |
---|
[707] | 509 | |
---|
[817] | 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: |
---|
[707] | 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 |
---|
[817] | 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 |
---|
[707] | 526 | \end{description} |
---|
[817] | 527 | %%% |
---|
| 528 | \gmcomment{ add the description of the smoothing: envelop topography...} |
---|
| 529 | %%% |
---|
[707] | 530 | |
---|
[817] | 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). |
---|
[707] | 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} |
---|
[817] | 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). |
---|
[707] | 564 | \end{description} |
---|
| 565 | |
---|
[817] | 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}). |
---|
[707] | 569 | |
---|
[817] | 570 | When a global ocean is coupled to an atmospheric model it is better to represent |
---|
[707] | 571 | all large water bodies (e.g, great lakes, Caspian sea...) even if the model |
---|
[817] | 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. |
---|
[707] | 577 | |
---|
| 578 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 579 | % z-coordinate and reference coordinate transformation |
---|
| 580 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 581 | \subsection[$z$-coordinate (\np{ln\_zco} or \key{zco})] |
---|
| 582 | {$z$-coordinate (\np{ln\_zco}=.true. or \key{zco}) and reference coordinate} |
---|
[707] | 583 | \label{DOM_zco} |
---|
| 584 | |
---|
| 585 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 586 | \begin{figure}[!tb] \label{Fig_zgr} \begin{center} |
---|
[998] | 587 | \includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_zgr.pdf} |
---|
[817] | 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.} |
---|
[707] | 591 | \end{center} \end{figure} |
---|
| 592 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 593 | |
---|
[817] | 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. |
---|
[707] | 604 | |
---|
[1224] | 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. |
---|
[707] | 608 | |
---|
[817] | 609 | For climate-related studies it is often desirable to concentrate the vertical resolution |
---|
[1224] | 610 | near the ocean surface. The following function is proposed as a standard for a |
---|
| 611 | $z$-coordinate (with either full or partial steps): |
---|
[707] | 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} |
---|
[817] | 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}). |
---|
[707] | 622 | |
---|
[817] | 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: |
---|
[707] | 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 | |
---|
[817] | 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}. |
---|
[707] | 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 |
---|
[817] | 647 | \pp{ppsur}=\pp{ppa0}=\pp{ppa1}=\pp{pp\_to\_be\_computed}, in \mdl{par\_oce}, |
---|
| 648 | and specifies instead the four following parameters: |
---|
[707] | 649 | \begin{itemize} |
---|
[817] | 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}) |
---|
[707] | 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} |
---|
[817] | 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. |
---|
[707] | 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 |
---|
[817] | 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 |
---|
[707] | 696 | 31 & \textbf{5250.23}& 5000.00 & \textbf{500.56} & 500.33 \\ \hline |
---|
| 697 | \end{tabular} \end{center} |
---|
[817] | 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}} |
---|
[707] | 701 | \end{table} |
---|
| 702 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 703 | |
---|
| 704 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 705 | % z-coordinate with partial step |
---|
| 706 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 707 | \subsection [$z$-coordinate with partial step (\np{ln\_zps})] |
---|
| 708 | {$z$-coordinate with partial step (\np{ln\_zps}=.true.)} |
---|
[707] | 709 | \label{DOM_zps} |
---|
[817] | 710 | %--------------------------------------------namdom------------------------------------------------------- |
---|
[707] | 711 | \namdisplay{namdom} |
---|
| 712 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 713 | |
---|
[817] | 714 | In $z$-coordinate partial step, the depths of the model levels are defined by the |
---|
[707] | 715 | reference analytical function $z_0 (k)$ as described in the previous |
---|
[817] | 716 | section, \emph{except} in the bottom layer. The thickness of the bottom layer is |
---|
[707] | 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 |
---|
[1224] | 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 |
---|
[707] | 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 |
---|
[1224] | 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 |
---|
[707] | 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 | % ------------------------------------------------------------------------------------------------------------- |
---|
[817] | 739 | \subsection [$s$-coordinate (\np{ln\_sco})] |
---|
| 740 | {$s$-coordinate (\np{ln\_sco}=true)} |
---|
[707] | 741 | \label{DOM_sco} |
---|
[817] | 742 | %------------------------------------------nam_zgr_sco--------------------------------------------------- |
---|
[707] | 743 | \namdisplay{nam_zgr_sco} |
---|
| 744 | %-------------------------------------------------------------------------------------------------------------- |
---|
[817] | 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: |
---|
[707] | 748 | \begin{equation} \label{DOM_sco_ana} |
---|
| 749 | \begin{split} |
---|
[817] | 750 | z(k) &= h(i,j) \; z_0(k) \\ |
---|
[707] | 751 | e_3(k) &= h(i,j) \; z_0'(k) |
---|
| 752 | \end{split} |
---|
| 753 | \end{equation} |
---|
[817] | 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 |
---|
[707] | 756 | surface to $1$ at the ocean bottom. The depth field $h$ is not necessary the ocean |
---|
[817] | 757 | depth, since a mixed step-like and bottom-following representation of the |
---|
[1224] | 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. |
---|
[707] | 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} |
---|
[817] | 771 | where $h_c$ is the thermocline depth and $\theta$ and $b$ are the surface and |
---|
[707] | 772 | bottom control parameters such that $0\leqslant \theta \leqslant 20$, and |
---|
[1224] | 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}). |
---|
[707] | 775 | |
---|
| 776 | %>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
| 777 | \begin{figure}[!tb] \label{Fig_sco_function} \begin{center} |
---|
[998] | 778 | \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_sco_function.pdf} |
---|
[1224] | 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} |
---|
[707] | 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 | |
---|
[817] | 790 | This option is described in the Report by Levier \textit{et al.} (2007), available on |
---|
| 791 | the \NEMO web site. |
---|
[707] | 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} |
---|
[817] | 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. |
---|
[707] | 815 | |
---|
[817] | 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. |
---|
[707] | 823 | |
---|
[817] | 824 | This is used to compute the island barotropic stream function used in the rigid lid |
---|
| 825 | computation (see \S\ref{MISC_solisl}). |
---|
| 826 | |
---|
[707] | 827 | From the \textit{mbathy} array, the mask fields are defined as follows: |
---|
| 828 | \begin{align*} |
---|
[817] | 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) \\ |
---|
[994] | 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) \\ |
---|
[817] | 834 | & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) |
---|
[707] | 835 | \end{align*} |
---|
| 836 | |
---|
[817] | 837 | \gmcomment{ STEVEN: are the dots multiplications?} |
---|
[707] | 838 | |
---|
[817] | 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}). |
---|
[707] | 846 | |
---|
[817] | 847 | %%% |
---|
| 848 | \gmcomment{ \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}. } |
---|
| 849 | %%% |
---|
| 850 | |
---|
[707] | 851 | % ================================================================ |
---|
[817] | 852 | % Time Discretisation |
---|
[707] | 853 | % ================================================================ |
---|
| 854 | \section{Time Discretisation} |
---|
| 855 | \label{DOM_nxt} |
---|
| 856 | |
---|
[817] | 857 | The time stepping used in \NEMO is a three level scheme that can be |
---|
| 858 | represented as follows: |
---|
[707] | 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} |
---|
[817] | 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 | |
---|
[707] | 868 | The choice of the time step used for this evaluation is discussed below as |
---|
[817] | 869 | well as the implications in terms of starting or restarting a model |
---|
[707] | 870 | simulation. Note that the time stepping is generally performed in a one step |
---|
[1224] | 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. |
---|
[707] | 873 | |
---|
[817] | 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. |
---|
[707] | 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 |
---|
[817] | 890 | leapfrog scheme. It is a time centred scheme, i.e. the RHS is evaluated at |
---|
[707] | 891 | time step $t$, the now time step. It is only used for non-diffusive terms, |
---|
[817] | 892 | that is momentum and tracer advection, pressure gradient, and Coriolis |
---|
[707] | 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 |
---|
[817] | 904 | to diverge into a physical and a computational mode. Time splitting can |
---|
[707] | 905 | be controlled through the use of an Asselin time filter (first designed by |
---|
[1224] | 906 | \citep{Robert1966} and more comprehensively studied by \citet{Asselin1972}), |
---|
| 907 | or by periodically reinitialising the leapfrog solution through a single |
---|
[817] | 908 | integration step with a two-level scheme. In \NEMO we follow the first |
---|
[707] | 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} |
---|
[817] | 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 |
---|
[707] | 920 | time filter has been preferred to other time differencing schemes such as |
---|
[817] | 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 |
---|
[707] | 924 | advective terms in the momentum and tracer equations, it avoids implicit |
---|
[817] | 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. |
---|
[707] | 928 | |
---|
[817] | 929 | \gmcomment{ |
---|
[707] | 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}: |
---|
[817] | 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: |
---|
[707] | 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} |
---|
[817] | 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: |
---|
[707] | 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} |
---|
[817] | 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} |
---|
[707] | 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 | |
---|
[817] | 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$: |
---|
[707] | 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 | |
---|
[817] | 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 : |
---|
[707] | 992 | \begin{equation} \label{Eq_DOM_nxt_euler} |
---|
[817] | 993 | x^{t+\Delta t} = x^{t-\Delta t} + 2 \, \Delta t \ {D_x}^{t-\Delta t} |
---|
[707] | 994 | \end{equation} |
---|
| 995 | |
---|
| 996 | This is diffusive in time and conditionally stable. For example, the |
---|
[817] | 997 | conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies2004}: |
---|
[707] | 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} |
---|
[1224] | 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. |
---|
[707] | 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 |
---|
[817] | 1014 | constraint on the time step. Two solutions are available in \NEMO to overcome |
---|
[707] | 1015 | the stability constraint: $(a)$ a forward time differencing scheme using a |
---|
[817] | 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 |
---|
[707] | 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} |
---|
[817] | 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: |
---|
[707] | 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} |
---|
[1224] | 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: |
---|
[707] | 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 | |
---|
[817] | 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}). |
---|
[707] | 1063 | |
---|
| 1064 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1065 | % Start/Restart strategy |
---|
| 1066 | % ------------------------------------------------------------------------------------------------------------- |
---|
| 1067 | \subsection{Start/Restart strategy} |
---|
| 1068 | \label{DOM_nxt_rst} |
---|
| 1069 | %--------------------------------------------namrun------------------------------------------- |
---|
| 1070 | \namdisplay{namrun} |
---|
| 1071 | %-------------------------------------------------------------------------------------------------------------- |
---|
| 1072 | |
---|
[817] | 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 | |
---|
[707] | 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, |
---|
[817] | 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. |
---|
[707] | 1086 | |
---|
[817] | 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 |
---|
[1224] | 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). |
---|
[817] | 1093 | %%% |
---|
| 1094 | \gmcomment{add here how to force the restart to contain only one time step for operational purposes} |
---|
| 1095 | %%% |
---|
[707] | 1096 | |
---|
[817] | 1097 | \gmcomment{ % add a subsection here |
---|
| 1098 | |
---|
[707] | 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 | |
---|
[817] | 1109 | } %% end add |
---|
[707] | 1110 | |
---|