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