Changeset 11622 for NEMO/trunk/doc/latex/NEMO/subfiles
- Timestamp:
- 2019-10-01T13:10:55+02:00 (5 years ago)
- Location:
- NEMO/trunk/doc/latex/NEMO/subfiles
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_DOM.tex
r11598 r11622 6 6 \label{chap:DOM} 7 7 8 % Missing things: 9 % - istate: description of the initial state ==> this has to be put elsewhere.. 10 % perhaps in MISC ? By the way the initialisation of T S and dynamics 11 % should be put outside of DOM routine (better with TRC staff and off-line 12 % tracers) 13 % -geo2ocean: how to switch from geographic to mesh coordinate 14 % - domclo: closed sea and lakes.... management of closea sea area : specific to global configuration, both forced and coupled 15 16 % {\em 4.0} & {\em Simon M\"{u}ller \& Andrew Coward} & 17 % {\em 18 % Compatibility changes Major simplification has moved many of the options to external domain configuration tools. 19 % (see \autoref{apdx:DOMCFG}) 20 % } \\ 21 % {\em 3.x} & {\em Rachid Benshila, Gurvan Madec \& S\'{e}bastien Masson} & 22 % {\em First version} \\ 8 % Missing things 9 % - istate: description of the initial state ==> this has to be put elsewhere.. 10 % perhaps in MISC ? By the way the initialisation of T S and dynamics 11 % should be put outside of DOM routine (better with TRC staff and off-line 12 % tracers) 13 % - geo2ocean: how to switch from geographic to mesh coordinate 14 % - domclo: closed sea and lakes.... 15 % management of closea sea area: specific to global cfg, both forced and coupled 23 16 24 17 \thispagestyle{plain} … … 30 23 {\footnotesize 31 24 \begin{tabularx}{\textwidth}{l||X|X} 32 Release & Author(s) & Modifications \\ 33 \hline 34 {\em 4.0} & {\em ...} & {\em ...} \\ 35 {\em 3.6} & {\em ...} & {\em ...} \\ 36 {\em 3.4} & {\em ...} & {\em ...} \\ 37 {\em <=3.4} & {\em ...} & {\em ...} 25 Release & 26 Author(s) & 27 Modifications \\ 28 \hline 29 {\em 4.0 } & 30 {\em Simon M\"{u}ller \& Andrew Coward \newline \newline 31 Simona Flavoni and Tim Graham } & 32 {\em Compatibility changes: many options moved to external domain configuration tools 33 (see \autoref{apdx:DOMCFG}). \newline 34 Updates } \\ 35 {\em 3.6 } & 36 {\em Rachid Benshila, Christian \'{E}th\'{e}, Pierre Mathiot and Gurvan Madec } & 37 {\em Updates } \\ 38 {\em $\leq$ 3.4 } & 39 {\em Gurvan Madec and S\'{e}bastien Masson } & 40 {\em First version } 38 41 \end{tabularx} 39 42 } … … 41 44 \clearpage 42 45 43 Having defined the continuous equations in \autoref{chap:MB} and chosen a time discretisation \autoref{chap:TD}, 46 Having defined the continuous equations in \autoref{chap:MB} and 47 chosen a time discretisation \autoref{chap:TD}, 44 48 we need to choose a grid for spatial discretisation and related numerical algorithms. 45 49 In the present chapter, we provide a general description of the staggered grid used in \NEMO, … … 54 58 \label{subsec:DOM_cell} 55 59 56 \begin{figure} [!tb]60 \begin{figure} 57 61 \centering 58 \includegraphics[width=0. 66\textwidth]{Fig_cell}62 \includegraphics[width=0.33\textwidth]{Fig_cell} 59 63 \caption[Arrangement of variables in the unit cell of space domain]{ 60 64 Arrangement of variables in the unit cell of space domain. 61 65 $t$ indicates scalar points where 62 66 temperature, salinity, density, pressure and horizontal divergence are defined. 63 $(u,v,w)$ indicates vector points, 64 and $f$ indicates vorticity points where 67 $(u,v,w)$ indicates vector points, and $f$ indicates vorticity points where 65 68 both relative and planetary vorticities are defined.} 66 69 \label{fig:DOM_cell} 67 70 \end{figure} 68 71 69 The numerical techniques used to solve the Primitive Equations in this model are based on the traditional,70 centred second-order finite difference approximation.72 The numerical techniques used to solve the Primitive Equations in this model are based on 73 the traditional, centred second-order finite difference approximation. 71 74 Special attention has been given to the homogeneity of the solution in the three spatial directions. 72 75 The arrangement of variables is the same in all directions. 73 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with vector points $(u, v, w)$ defined in 74 the centre of each face of the cells (\autoref{fig:DOM_cell}). 75 This is the generalisation to three dimensions of the well-known ``C'' grid in Arakawa's classification 76 \citep{mesinger.arakawa_bk76}. 77 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each vertical edge and 78 the barotropic stream function $\psi$ is defined at horizontal points overlying the $\zeta$ and $f$-points. 79 80 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by the transformation that 81 gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 82 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on \autoref{tab:DOM_cell}. 83 In all the following, subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of 84 the grid-point where the scale factors are defined. 76 It consists of cells centred on scalar points ($t$, $S$, $p$, $\rho$) with 77 vector points $(u, v, w)$ defined in the centre of each face of the cells (\autoref{fig:DOM_cell}). 78 This is the generalisation to three dimensions of the well-known ``C'' grid in 79 Arakawa's classification \citep{mesinger.arakawa_bk76}. 80 The relative and planetary vorticity, $\zeta$ and $f$, are defined in the centre of each 81 vertical edge and the barotropic stream function $\psi$ is defined at horizontal points overlying 82 the $\zeta$ and $f$-points. 83 84 The ocean mesh (\ie\ the position of all the scalar and vector points) is defined by 85 the transformation that gives $(\lambda,\varphi,z)$ as a function of $(i,j,k)$. 86 The grid-points are located at integer or integer and a half value of $(i,j,k)$ as indicated on 87 \autoref{tab:DOM_cell}. 88 In all the following, 89 subscripts $u$, $v$, $w$, $f$, $uw$, $vw$ or $fw$ indicate the position of the grid-point where 90 the scale factors are defined. 85 91 Each scale factor is defined as the local analytical value provided by \autoref{eq:MB_scale_factors}. 86 92 As a result, the mesh on which partial derivatives $\pd[]{\lambda}$, $\pd[]{\varphi}$ and 87 93 $\pd[]{z}$ are evaluated is a uniform mesh with a grid size of unity. 88 Discrete partial derivatives are formulated by the traditional, centred second order finite difference approximation 89 while the scale factors are chosen equal to their local analytical value. 94 Discrete partial derivatives are formulated by 95 the traditional, centred second order finite difference approximation while 96 the scale factors are chosen equal to their local analytical value. 90 97 An important point here is that the partial derivative of the scale factors must be evaluated by 91 98 centred finite difference approximation, not from their analytical expression. 92 This preserves the symmetry of the discrete set of equations and therefore satisfies many of93 the continuous properties (see \autoref{apdx:INVARIANTS}).99 This preserves the symmetry of the discrete set of equations and 100 therefore satisfies many of the continuous properties (see \autoref{apdx:INVARIANTS}). 94 101 A similar, related remark can be made about the domain size: 95 when needed, an area, volume, or the total ocean depth must be evaluated as the product or sum of the relevant scale factors96 (see \autoref{eq:DOM_bar} in the next section).97 98 \begin{table} [!tb]102 when needed, an area, volume, or the total ocean depth must be evaluated as 103 the product or sum of the relevant scale factors (see \autoref{eq:DOM_bar} in the next section). 104 105 \begin{table} 99 106 \centering 100 \begin{tabular}{| p{46pt}|p{56pt}|p{56pt}|p{56pt}|}101 \hline 102 t & $i$ & $j $ & $k $ \\103 \hline 104 u 105 \hline 106 v & $i$ & $j + 1/2$ & $k $ \\107 \hline 108 w & $i$ & $j $ & $k + 1/2$ \\109 \hline 110 f 111 \hline 112 uw 113 \hline 114 vw & $i$ & $j + 1/2$ & $k + 1/2$ \\115 \hline 116 fw 107 \begin{tabular}{|l|l|l|l|} 108 \hline 109 t & $i $ & $j $ & $k $ \\ 110 \hline 111 u & $i + 1/2$ & $j $ & $k $ \\ 112 \hline 113 v & $i $ & $j + 1/2$ & $k $ \\ 114 \hline 115 w & $i $ & $j $ & $k + 1/2$ \\ 116 \hline 117 f & $i + 1/2$ & $j + 1/2$ & $k $ \\ 118 \hline 119 uw & $i + 1/2$ & $j $ & $k + 1/2$ \\ 120 \hline 121 vw & $i $ & $j + 1/2$ & $k + 1/2$ \\ 122 \hline 123 fw & $i + 1/2$ & $j + 1/2$ & $k + 1/2$ \\ 117 124 \hline 118 125 \end{tabular} … … 120 127 Location of grid-points as a function of integer or 121 128 integer and a half value of the column, line or level. 122 This indexing is only used for the writing of the semi 129 This indexing is only used for the writing of the semi-discrete equations. 123 130 In the code, the indexing uses integer values only and 124 131 is positive downwards in the vertical with $k=1$ at the surface. … … 137 144 firstly, there is no ambiguity in the scale factors appearing in the discrete equations, 138 145 since they are first introduced in the continuous equations; 139 secondly, analytical transformations encourage good practice by the definition of smoothly varying grids 140 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) \citep{treguier.dukowicz.ea_JGR96}. 146 secondly, analytical transformations encourage good practice by 147 the definition of smoothly varying grids 148 (rather than allowing the user to set arbitrary jumps in thickness between adjacent layers) 149 \citep{treguier.dukowicz.ea_JGR96}. 141 150 An example of the effect of such a choice is shown in \autoref{fig:DOM_zgr_e3}. 142 \begin{figure} [!t]151 \begin{figure} 143 152 \centering 144 \includegraphics[width=0. 66\textwidth]{Fig_zgr_e3}153 \includegraphics[width=0.5\textwidth]{Fig_zgr_e3} 145 154 \caption[Comparison of grid-point position, vertical grid-size and scale factors]{ 146 155 Comparison of (a) traditional definitions of grid-point position and grid-size in the vertical, … … 159 168 \label{subsec:DOM_operators} 160 169 161 Given the values of a variable $q$ at adjacent points, the differencing and averaging operators at162 the midpoint between them are:170 Given the values of a variable $q$ at adjacent points, 171 the differencing and averaging operators at the midpoint between them are: 163 172 \begin{alignat*}{2} 164 173 % \label{eq:DOM_di_mi} … … 168 177 169 178 Similar operators are defined with respect to $i + 1/2$, $j$, $j + 1/2$, $k$, and $k + 1/2$. 170 Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, the gradient of a variable $q$ defined at a $t$-point has 171 its three components defined at $u$-, $v$- and $w$-points while its Laplacian is defined at the $t$-point. 179 Following \autoref{eq:MB_grad} and \autoref{eq:MB_lap}, 180 the gradient of a variable $q$ defined at a $t$-point has 181 its three components defined at $u$-, $v$- and $w$-points while 182 its Laplacian is defined at the $t$-point. 172 183 These operators have the following discrete forms in the curvilinear $s$-coordinates system: 173 184 \[ … … 177 188 + \frac{1}{e_{3w}} \delta_{k + 1/2} [q] \; \, \vect k 178 189 \] 179 \ begin{multline*}190 \[ 180 191 % \label{eq:DOM_lap} 181 192 \Delta q \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} … … 184 195 + \frac{1}{e_{3t}} 185 196 \delta_k \lt[ \frac{1 }{e_{3w}} \; \delta_{k + 1/2} [q] \rt] 186 \end{multline*} 187 188 Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, a vector $\vect A = (a_1,a_2,a_3)$ defined at 189 vector points $(u,v,w)$ has its three curl components defined at $vw$-, $uw$, and $f$-points, and 197 \] 198 199 Following \autoref{eq:MB_curl} and \autoref{eq:MB_div}, 200 a vector $\vect A = (a_1,a_2,a_3)$ defined at vector points $(u,v,w)$ has 201 its three curl components defined at $vw$-, $uw$, and $f$-points, and 190 202 its divergence defined at $t$-points: 191 \begin{multline }203 \begin{multline*} 192 204 % \label{eq:DOM_curl} 193 205 \nabla \times \vect A \equiv \frac{1}{e_{2v} \, e_{3vw}} … … 200 212 \Big[ \delta_{i + 1/2} (e_{2v} \, a_2) 201 213 - \delta_{j + 1/2} (e_{1u} \, a_1) \Big] \vect k 202 \end{multline }203 \ begin{equation}214 \end{multline*} 215 \[ 204 216 % \label{eq:DOM_div} 205 217 \nabla \cdot \vect A \equiv \frac{1}{e_{1t} \, e_{2t} \, e_{3t}} 206 218 \Big[ \delta_i (e_{2u} \, e_{3u} \, a_1) + \delta_j (e_{1v} \, e_{3v} \, a_2) \Big] 207 219 + \frac{1}{e_{3t}} \delta_k (a_3) 208 \ end{equation}209 210 The vertical average over the whole water column is denoted by an overbar and is for211 a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas):220 \] 221 222 The vertical average over the whole water column is denoted by an overbar and 223 is for a masked field $q$ (\ie\ a quantity that is equal to zero inside solid areas): 212 224 \begin{equation} 213 225 \label{eq:DOM_bar} … … 215 227 \end{equation} 216 228 where $H_q$ is the ocean depth, which is the masked sum of the vertical scale factors at $q$ points, 217 $k^b$ and $k^o$ are the bottom and surface $k$-indices, and the symbol $\sum \limits_k$ refers to a summation over 218 all grid points of the same type in the direction indicated by the subscript (here $k$). 229 $k^b$ and $k^o$ are the bottom and surface $k$-indices, 230 and the symbol $\sum \limits_k$ refers to a summation over all grid points of the same type in 231 the direction indicated by the subscript (here $k$). 219 232 220 233 In continuous form, the following properties are satisfied: … … 226 239 \end{gather} 227 240 228 It is straightforward to demonstrate that these properties are verified locally in discrete form as soon as229 the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at241 It is straightforward to demonstrate that these properties are verified locally in discrete form as 242 soon as the scalar $q$ is taken at $t$-points and the vector $\vect A$ has its components defined at 230 243 vector points $(u,v,w)$. 231 244 232 245 Let $a$ and $b$ be two fields defined on the mesh, with a value of zero inside continental areas. 233 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and $\delta_k$) 234 are skew-symmetric linear operators, and further that the averaging operators $\overline{\cdots}^{\, i}$, 235 $\overline{\cdots}^{\, j}$ and $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 246 It can be shown that the differencing operators ($\delta_i$, $\delta_j$ and 247 $\delta_k$) are skew-symmetric linear operators, 248 and further that the averaging operators ($\overline{\cdots}^{\, i}$, $\overline{\cdots}^{\, j}$ and 249 $\overline{\cdots}^{\, k}$) are symmetric linear operators, \ie 236 250 \begin{alignat}{4} 237 251 \label{eq:DOM_di_adj} … … 241 255 \end{alignat} 242 256 243 In other words, the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 257 In other words, 258 the adjoint of the differencing and averaging operators are $\delta_i^* = \delta_{i + 1/2}$ and 244 259 $(\overline{\cdots}^{\, i})^* = \overline{\cdots}^{\, i + 1/2}$, respectively. 245 260 These two properties will be used extensively in the \autoref{apdx:INVARIANTS} to … … 250 265 \label{subsec:DOM_Num_Index} 251 266 252 \begin{figure} [!tb]267 \begin{figure} 253 268 \centering 254 \includegraphics[width=0. 66\textwidth]{Fig_index_hor}269 \includegraphics[width=0.33\textwidth]{Fig_index_hor} 255 270 \caption[Horizontal integer indexing]{ 256 271 Horizontal integer indexing used in the \fortran\ code. … … 261 276 262 277 The array representation used in the \fortran\ code requires an integer indexing. 263 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with the use of 264 integer values for $t$-points only while all the other points involve integer and a half values. 278 However, the analytical definition of the mesh (see \autoref{subsec:DOM_cell}) is associated with 279 the use of integer values for $t$-points only while 280 all the other points involve integer and a half values. 265 281 Therefore, a specific integer indexing has been defined for points other than $t$-points 266 282 (\ie\ velocity and vorticity grid-points). 267 Furthermore, the direction of the vertical indexing has been reversed and the surface level set at $k = 1$. 283 Furthermore, the direction of the vertical indexing has been reversed and 284 the surface level set at $k = 1$. 268 285 269 286 %% ================================================================================================= … … 281 298 \label{subsec:DOM_Num_Index_vertical} 282 299 283 In the vertical, the chosen indexing requires special attention since the direction of the $k$-axis in284 the \fortran\ code is the reverse of that used in the semi -discrete equations and285 given in \autoref{subsec:DOM_cell}.286 The sea surface corresponds to the $w$-level $k = 1$, which is the same index as the $t$-level just below287 (\autoref{fig:DOM_index_vert}).300 In the vertical, the chosen indexing requires special attention since 301 the direction of the $k$-axis in the \fortran\ code is the reverse of 302 that used in the semi-discrete equations and given in \autoref{subsec:DOM_cell}. 303 The sea surface corresponds to the $w$-level $k = 1$, 304 which is the same index as the $t$-level just below (\autoref{fig:DOM_index_vert}). 288 305 The last $w$-level ($k = jpk$) either corresponds to or is below the ocean floor while 289 306 the last $t$-level is always outside the ocean domain (\autoref{fig:DOM_index_vert}). 290 307 Note that a $w$-point and the directly underlaying $t$-point have a common $k$ index 291 308 (\ie\ $t$-points and their nearest $w$-point neighbour in negative index direction), 292 in contrast to the indexing on the horizontal plane where the $t$-point has the same index as 293 the nearest velocity points in the positive direction of the respective horizontal axis index 309 in contrast to the indexing on the horizontal plane where 310 the $t$-point has the same index as the nearest velocity points in 311 the positive direction of the respective horizontal axis index 294 312 (compare the dashed area in \autoref{fig:DOM_index_hor} and \autoref{fig:DOM_index_vert}). 295 313 Since the scale factors are chosen to be strictly positive, … … 298 316 accommodate the opposing vertical index directions in implementation and documentation. 299 317 300 \begin{figure} [!pt]318 \begin{figure} 301 319 \centering 302 \includegraphics[width=0. 66\textwidth]{Fig_index_vert}320 \includegraphics[width=0.33\textwidth]{Fig_index_vert} 303 321 \caption[Vertical integer indexing]{ 304 322 Vertical integer indexing used in the \fortran\ code. … … 314 332 315 333 Two typical methods are available to specify the spatial domain configuration; 316 they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in namelist \nam{cfg}{cfg}. 334 they can be selected using parameter \np{ln_read_cfg}{ln\_read\_cfg} parameter in 335 namelist \nam{cfg}{cfg}. 317 336 318 337 If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.true.}, 319 the domain-specific parameters and fields are read from a netCDF input file, 320 whose name (without its .nc suffix) can be specified as the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 338 the domain-specific parameters and fields are read from a NetCDF input file, 339 whose name (without its .nc suffix) can be specified as 340 the value of the \np{cn_domcfg}{cn\_domcfg} parameter in namelist \nam{cfg}{cfg}. 321 341 322 342 If \np{ln_read_cfg}{ln\_read\_cfg} is set to \forcode{.false.}, … … 324 344 subroutines \mdl{usrdef\_hgr} and \mdl{usrdef\_zgr}. 325 345 These subroutines can be supplied in the \path{MY_SRC} directory of the configuration, 326 and default versions that configure the spatial domain for the GYRE reference configuration are present in327 the \path{./src/OCE/USR} directory.346 and default versions that configure the spatial domain for the GYRE reference configuration are 347 present in the \path{./src/OCE/USR} directory. 328 348 329 349 In version 4.0 there are no longer any options for reading complex bathymetries and … … 332 352 to run similar models with and without partial bottom boxes and/or sigma-coordinates, 333 353 supporting such choices leads to overly complex code. 334 Worse still is the difficulty of ensuring the model configurations intended to be identical are indeed so when 335 the model domain itself can be altered by runtime selections. 336 The code previously used to perform vertical discretisation has been incorporated into an external tool 337 (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 338 339 The next subsections summarise the parameter and fields related to the configuration of the whole model domain. 340 These represent the minimum information that must be provided either via the \np{cn_domcfg}{cn\_domcfg} file or set by code 341 inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 354 Worse still is the difficulty of ensuring the model configurations intended to be identical are 355 indeed so when the model domain itself can be altered by runtime selections. 356 The code previously used to perform vertical discretisation has been incorporated into 357 an external tool (\path{./tools/DOMAINcfg}) which is briefly described in \autoref{apdx:DOMCFG}. 358 359 The next subsections summarise the parameter and fields related to 360 the configuration of the whole model domain. 361 These represent the minimum information that must be provided either via 362 the \np{cn_domcfg}{cn\_domcfg} file or 363 set by code inserted into user-supplied versions of the \texttt{usrdef\_*} subroutines. 342 364 The requirements are presented in three sections: 343 365 the domain size (\autoref{subsec:DOM_size}), the horizontal mesh (\autoref{subsec:DOM_hgr}), … … 348 370 \label{subsec:DOM_size} 349 371 350 The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and \jp{jpkglo} for351 the $i$, $j$ and $k$ directions, respectively.352 Note, that the variables \texttt{jpi} and \texttt{jpj} refer to the size of each processor subdomain when353 the code is run in parallel using domain decomposition (\key{mpp\_mpi} defined,354 see \autoref{sec:LBC_mpp}).372 The total size of the computational domain is set by the parameters \jp{jpiglo}, \jp{jpjglo} and 373 \jp{jpkglo} for the $i$, $j$ and $k$ directions, respectively. 374 Note, that the variables \texttt{jpi} and \texttt{jpj} refer to 375 the size of each processor subdomain when the code is run in parallel using domain decomposition 376 (\key{mpp\_mpi} defined, see \autoref{sec:LBC_mpp}). 355 377 356 378 The name of the configuration is set through parameter \np{cn_cfg}{cn\_cfg}, … … 360 382 361 383 The global lateral boundary condition type is selected from 8 options using parameter \jp{jperio}. 362 See \autoref{sec:LBC_jperio} for details on the available options and the corresponding values for \jp{jperio}. 384 See \autoref{sec:LBC_jperio} for details on the available options and 385 the corresponding values for \jp{jperio}. 363 386 364 387 %% ================================================================================================= … … 370 393 \label{sec:DOM_hgr_fields} 371 394 372 The explicit specification of a range of mesh-related fields are required for the definition of a configuration. 395 The explicit specification of a range of mesh-related fields are required for 396 the definition of a configuration. 373 397 These include: 374 398 375 399 \begin{clines} 376 int jpiglo, jpjglo, jpkglo /* global domain sizes*/377 int jperio /* lateral global domain b.c.*/378 double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively)*/379 double gphit, gphiu, gphiv, gphif /* geographic latitude*/380 double e1t, e1u, e1v, e1f /* horizontal scale factors*/381 double e2t, e2u, e2v, e2f /* horizontal scale factors*/400 int jpiglo, jpjglo, jpkglo /* global domain sizes */ 401 int jperio /* lateral global domain b.c. */ 402 double glamt, glamu, glamv, glamf /* geographic longitude (t,u,v and f points respectively) */ 403 double gphit, gphiu, gphiv, gphif /* geographic latitude */ 404 double e1t, e1u, e1v, e1f /* horizontal scale factors */ 405 double e2t, e2u, e2v, e2f /* horizontal scale factors */ 382 406 \end{clines} 383 407 … … 393 417 394 418 \begin{clines} 395 /* Optional:*/396 int ORCA, ORCA_index /* configuration name, configuration resolution*/397 double e1e2u, e1e2v /* U and V surfaces (if grid size reduction in some straits)*/398 double ff_f, ff_t /* Coriolis parameter (if not on the sphere)*/419 /* Optional: */ 420 int ORCA, ORCA_index /* configuration name, configuration resolution */ 421 double e1e2u, e1e2v /* U and V surfaces (if grid size reduction in some straits) */ 422 double ff_f, ff_t /* Coriolis parameter (if not on the sphere) */ 399 423 \end{clines} 400 424 … … 403 427 This is particularly useful for locations such as Gibraltar or Indonesian Throughflow pinch-points 404 428 (see \autoref{sec:MISC_strait} for illustrated examples). 405 The key is to reduce the faces of $T$-cell (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 429 The key is to reduce the faces of $T$-cell 430 (\ie\ change the value of the horizontal scale factors at $u$- or $v$-point) but 406 431 not the volume of the cells. 407 432 Doing otherwise can lead to numerical instability issues. 408 433 In normal operation the surface areas are computed from $e1u * e2u$ and $e1v * e2v$ but 409 434 in cases where a gridsize reduction is required, 410 the unaltered surface areas at $u$ and $v$ grid points (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or 411 pre-computed in \mdl{usrdef\_hgr}. 412 If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and the internal computation is suppressed. 413 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should set 414 the surface-area computation flag: 435 the unaltered surface areas at $u$ and $v$ grid points 436 (\texttt{e1e2u} and \texttt{e1e2v}, respectively) must be read or pre-computed in \mdl{usrdef\_hgr}. 437 If these arrays are present in the \np{cn_domcfg}{cn\_domcfg} file they are read and 438 the internal computation is suppressed. 439 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{e1e2u} and \texttt{e1e2v} should 440 set the surface-area computation flag: 415 441 \texttt{ie1e2u\_v} to a non-zero value to suppress their re-computation. 416 442 417 443 \smallskip 418 444 Similar logic applies to the other optional fields: 419 \texttt{ff\_f} and \texttt{ff\_t} which can be used to provide the Coriolis parameter at F- and T-points respectively if 420 the mesh is not on a sphere. 421 If present these fields will be read and used and the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 422 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should set 423 the Coriolis computation flag: 445 \texttt{ff\_f} and \texttt{ff\_t} which can be used to 446 provide the Coriolis parameter at F- and T-points respectively if the mesh is not on a sphere. 447 If present these fields will be read and used and 448 the normal calculation ($2 * \Omega * \sin(\varphi)$) suppressed. 449 Versions of \mdl{usrdef\_hgr} which set their own values of \texttt{ff\_f} and \texttt{ff\_t} should 450 set the Coriolis computation flag: 424 451 \texttt{iff} to a non-zero value to suppress their re-computation. 425 452 426 Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to those of $t$ points,427 th us no specific arrays are defined at $w$ points.453 Note that longitudes, latitudes, and scale factors at $w$ points are exactly equal to 454 those of $t$ points, thus no specific arrays are defined at $w$ points. 428 455 429 456 %% ================================================================================================= 430 457 \subsection[Vertical grid (\textit{domzgr.F90})]{Vertical grid (\protect\mdl{domzgr})} 431 458 \label{subsec:DOM_zgr} 459 432 460 \begin{listing} 433 461 \nlst{namdom} … … 438 466 In the vertical, the model mesh is determined by four things: 439 467 \begin{enumerate} 440 \item the bathymetry given in meters; 441 \item the number of levels of the model (\jp{jpk}); 442 \item the analytical transformation $z(i,j,k)$ and the vertical scale factors (derivatives of the transformation); and 443 \item the masking system, \ie\ the number of wet model levels at each 444 $(i,j)$ location of the horizontal grid. 468 \item the bathymetry given in meters; 469 \item the number of levels of the model (\jp{jpk}); 470 \item the analytical transformation $z(i,j,k)$ and the vertical scale factors 471 (derivatives of the transformation); and 472 \item the masking system, 473 \ie\ the number of wet model levels at each $(i,j)$ location of the horizontal grid. 445 474 \end{enumerate} 446 475 447 \begin{figure} [!tb]476 \begin{figure} 448 477 \centering 449 \includegraphics[width=0. 66\textwidth]{Fig_z_zps_s_sps}478 \includegraphics[width=0.5\textwidth]{Fig_z_zps_s_sps} 450 479 \caption[Ocean bottom regarding coordinate systems ($z$, $s$ and hybrid $s-z$)]{ 451 480 The ocean bottom as seen by the model: 452 (a) $z$-coordinate with full step, 453 (b) $z$-coordinate with partial step, 454 (c) $s$-coordinate: terrain following representation, 455 (d) hybrid $s-z$ coordinate, 456 (e) hybrid $s-z$ coordinate with partial step, and 457 (f) same as (e) but in the non-linear free surface (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 458 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 481 \begin{enumerate*}[label={(\alph*)}] 482 \item $z$-coordinate with full step, 483 \item $z$-coordinate with partial step, 484 \item $s$-coordinate: terrain following representation, 485 \item hybrid $s-z$ coordinate, 486 \item hybrid $s-z$ coordinate with partial step, and 487 \item same as (e) but in the non-linear free surface 488 (\protect\np[=.false.]{ln_linssh}{ln\_linssh}). 489 \end{enumerate*} 490 Note that the non-linear free surface can be used with any of the 5 coordinates (a) to (e).} 459 491 \label{fig:DOM_z_zps_s_sps} 460 492 \end{figure} … … 463 495 it is not intended to be an option which can be changed in the middle of an experiment. 464 496 The one exception to this statement being the choice of linear or non-linear free surface. 465 In v4.0 the linear free surface option is implemented as a special case of the non-linear free surface. 497 In v4.0 the linear free surface option is implemented as 498 a special case of the non-linear free surface. 466 499 This is computationally wasteful since it uses the structures for time-varying 3D metrics 467 500 for fields that (in the linear free surface case) are fixed. 468 However, the linear free-surface is rarely used and implementing it this way means 469 a single configuration file can support both options. 470 471 By default a non-linear free surface is used (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 472 the coordinate follow the time-variation of the free surface so that the transformation is time dependent: 473 $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 474 When a linear free surface is assumed (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 475 the vertical coordinates are fixed in time, but the seawater can move up and down across the $z_0$ surface 501 However, the linear free-surface is rarely used and 502 implementing it this way means a single configuration file can support both options. 503 504 By default a non-linear free surface is used 505 (\np{ln_linssh}{ln\_linssh} set to \forcode{=.false.} in \nam{dom}{dom}): 506 the coordinate follow the time-variation of the free surface so that 507 the transformation is time dependent: $z(i,j,k,t)$ (\eg\ \autoref{fig:DOM_z_zps_s_sps}f). 508 When a linear free surface is assumed 509 (\np{ln_linssh}{ln\_linssh} set to \forcode{=.true.} in \nam{dom}{dom}), 510 the vertical coordinates are fixed in time, but 511 the seawater can move up and down across the $z_0$ surface 476 512 (in other words, the top of the ocean in not a rigid lid). 477 513 478 514 Note that settings: 479 \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} mentioned in the following sections 480 appear to be namelist options but they are no longer truly namelist options for \NEMO. 515 \np{ln_zco}{ln\_zco}, \np{ln_zps}{ln\_zps}, \np{ln_sco}{ln\_sco} and \np{ln_isfcav}{ln\_isfcav} 516 mentioned in the following sections appear to be namelist options but 517 they are no longer truly namelist options for \NEMO. 481 518 Their value is written to and read from the domain configuration file and 482 519 they should be treated as fixed parameters for a particular configuration. 483 They are namelist options for the \texttt{DOMAINcfg} tool that can be used to build the configuration file and484 serve both to provide a record of the choices made whilst building the configuration and 485 to trigger appropriate code blocks within \NEMO.520 They are namelist options for the \texttt{DOMAINcfg} tool that can be used to 521 build the configuration file and serve both to provide a record of the choices made whilst 522 building the configuration and to trigger appropriate code blocks within \NEMO. 486 523 These values should not be altered in the \np{cn_domcfg}{cn\_domcfg} file. 487 524 … … 501 538 A further choice related to vertical coordinate concerns 502 539 the presence (or not) of ocean cavities beneath ice shelves within the model domain. 503 A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that the domain contains ocean cavities, 540 A setting of \np{ln_isfcav}{ln\_isfcav} as \forcode{.true.} indicates that 541 the domain contains ocean cavities, 504 542 otherwise the top, wet layer of the ocean will always be at the ocean surface. 505 543 This option is currently only available for $z$- or $zps$-coordinates. 506 544 In the latter case, partial steps are also applied at the ocean/ice shelf interface. 507 545 508 Within the model, the arrays describing the grid point depths and vertical scale factors are three set of 509 three dimensional arrays $(i,j,k)$ defined at \textit{before}, \textit{now} and \textit{after} time step. 546 Within the model, 547 the arrays describing the grid point depths and vertical scale factors are 548 three set of three dimensional arrays $(i,j,k)$ defined at 549 \textit{before}, \textit{now} and \textit{after} time step. 510 550 The time at which they are defined is indicated by a suffix: $\_b$, $\_n$, or $\_a$, respectively. 511 551 They are updated at each model time step. … … 534 574 \end{clines} 535 575 536 This set of vertical metrics is sufficient to describe the initial depth and thickness of every gridcell in537 the model regardless of the choice of vertical coordinate.576 This set of vertical metrics is sufficient to describe the initial depth and thickness of 577 every gridcell in the model regardless of the choice of vertical coordinate. 538 578 With constant z-levels, e3 metrics will be uniform across each horizontal level. 539 579 In the partial step case each e3 at the \jp{bottom\_level} … … 541 581 may vary from its horizontal neighbours. 542 582 And, in s-coordinates, variations can occur throughout the water column. 543 With the non-linear free-surface, all the coordinates behave more like the s-coordinate in 544 thatvariations occur throughout the water column with displacements related to the sea surface height.583 With the non-linear free-surface, all the coordinates behave more like the s-coordinate in that 584 variations occur throughout the water column with displacements related to the sea surface height. 545 585 These variations are typically much smaller than those arising from bottom fitted coordinates. 546 586 The values for vertical metrics supplied in the domain configuration file can be considered as 547 587 those arising from a flat sea surface with zero elevation. 548 588 549 The \jp{bottom\_level} and \jp{top\_level} 2D arrays define the \jp{bottom\_level} and top wet levels in each grid column. 589 The \jp{bottom\_level} and \jp{top\_level} 2D arrays define 590 the \jp{bottom\_level} and top wet levels in each grid column. 550 591 Without ice cavities, \jp{top\_level} is essentially a land mask (0 on land; 1 everywhere else). 551 592 With ice cavities, \jp{top\_level} determines the first wet point below the overlying ice shelf. … … 556 597 557 598 From \jp{top\_level} and \jp{bottom\_level} fields, the mask fields are defined as follows: 558 \begin{alignat*}{2} 559 tmask(i,j,k) &= & & 560 \begin{cases} 561 0 &\text{if $ k < top\_level(i,j)$} \\ 562 1 &\text{if $bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 563 0 &\text{if $ k > bottom\_level(i,j)$} 564 \end{cases} 565 \\ 566 umask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 567 vmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j + 1,k) \\ 568 fmask(i,j,k) &= & &tmask(i,j,k) * tmask(i + 1,j, k) \\ 569 & &* &tmask(i,j,k) * tmask(i + 1,j, k) \\ 570 wmask(i,j,k) &= & &tmask(i,j,k) * tmask(i ,j,k - 1) \\ 571 \text{with~} wmask(i,j,1) &= & &tmask(i,j,1) 572 \end{alignat*} 599 \begin{align*} 600 tmask(i,j,k) &= 601 \begin{cases} 602 0 &\text{if $ k < top\_level(i,j)$} \\ 603 1 &\text{if $ bottom\_level(i,j) \leq k \leq top\_level(i,j)$} \\ 604 0 &\text{if $k > bottom\_level(i,j) $} 605 \end{cases} \\ 606 umask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j, k) \\ 607 vmask(i,j,k) &= tmask(i,j,k) * tmask(i ,j + 1,k) \\ 608 fmask(i,j,k) &= tmask(i,j,k) * tmask(i + 1,j, k) * tmask(i,j,k) * tmask(i + 1,j, k) \\ 609 wmask(i,j,k) &= tmask(i,j,k) * tmask(i ,j,k - 1) \\ 610 \text{with~} wmask(i,j,1) &= tmask(i,j,1) 611 \end{align*} 573 612 574 613 Note that, without ice shelves cavities, 575 masks at $t-$ and $w-$points are identical with the numerical indexing used (\autoref{subsec:DOM_Num_Index}). 576 Nevertheless, $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 614 masks at $t-$ and $w-$points are identical with the numerical indexing used 615 (\autoref{subsec:DOM_Num_Index}). 616 Nevertheless, 617 $wmask$ are required with ocean cavities to deal with the top boundary (ice shelf/ocean interface) 577 618 exactly in the same way as for the bottom boundary. 578 619 … … 588 629 \label{subsec:DOM_closea} 589 630 590 When a global ocean is coupled to an atmospheric model it is better to represent all large water bodies591 (\eg\ Great Lakes, Caspian sea \dots) even if the model resolution does not allow their communication with 592 the rest of the ocean.631 When a global ocean is coupled to an atmospheric model it is better to 632 represent all large water bodies (\eg\ Great Lakes, Caspian sea, \dots) even if 633 the model resolution does not allow their communication with the rest of the ocean. 593 634 This is unnecessary when the ocean is forced by fixed atmospheric conditions, 594 635 so these seas can be removed from the ocean domain. 595 The user has the option to set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and 596 to optionally decide on the fate of any freshwater imbalance over the area. 597 The options are explained in \autoref{sec:MISC_closea} but it should be noted here that 598 a successful use of these options requires appropriate mask fields to be present in the domain configuration file. 636 The user has the option to 637 set the bathymetry in closed seas to zero (see \autoref{sec:MISC_closea}) and to 638 optionally decide on the fate of any freshwater imbalance over the area. 639 The options are explained in \autoref{sec:MISC_closea} but 640 it should be noted here that a successful use of these options requires 641 appropriate mask fields to be present in the domain configuration file. 599 642 Among the possibilities are: 600 643 601 644 \begin{clines} 602 int closea_mask /* non-zero values in closed sea areas for optional masking*/603 int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only)*/604 int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp)*/645 int closea_mask /* non-zero values in closed sea areas for optional masking */ 646 int closea_mask_rnf /* non-zero values in closed sea areas with runoff locations (precip only) */ 647 int closea_mask_emp /* non-zero values in closed sea areas with runoff locations (total emp) */ 605 648 \end{clines} 606 649 … … 610 653 611 654 Most of the arrays relating to a particular ocean model configuration discussed in this chapter 612 (grid-point position, scale factors) 613 can be saved in a file if 614 namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to\forcode{.true.};655 (grid-point position, scale factors) can be saved in a file if 656 namelist parameter \np{ln_write_cfg}{ln\_write\_cfg} (namelist \nam{cfg}{cfg}) is set to 657 \forcode{.true.}; 615 658 the output filename is set through parameter \np{cn_domcfg_out}{cn\_domcfg\_out}. 616 659 This is only really useful if … … 619 662 620 663 Alternatively, all the arrays relating to a particular ocean model configuration 621 (grid-point position, scale factors, depths and masks) 622 can be saved in a file called \texttt{mesh\_mask} if 623 namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to \forcode{.true.}. 664 (grid-point position, scale factors, depths and masks) can be saved in 665 a file called \texttt{mesh\_mask} if 666 namelist parameter \np{ln_meshmask}{ln\_meshmask} (namelist \nam{dom}{dom}) is set to 667 \forcode{.true.}. 624 668 This file contains additional fields that can be useful for post-processing applications. 625 669 … … 627 671 \section[Initial state (\textit{istate.F90} and \textit{dtatsd.F90})]{Initial state (\protect\mdl{istate} and \protect\mdl{dtatsd})} 628 672 \label{sec:DOM_DTA_tsd} 673 629 674 \begin{listing} 630 675 \nlst{namtsd} … … 638 683 639 684 \begin{description} 640 \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on the model grid itself or on their native input data grids. 641 In the latter case, the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 685 \item [{\np[=.true.]{ln_tsd_init}{ln\_tsd\_init}}] Use T and S input files that can be given on 686 the model grid itself or on their native input data grids. 687 In the latter case, 688 the data will be interpolated on-the-fly both in the horizontal and the vertical to the model grid 642 689 (see \autoref{subsec:SBC_iof}). 643 The information relating to the input files are specified in the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 690 The information relating to the input files are specified in 691 the \np{sn_tem}{sn\_tem} and \np{sn_sal}{sn\_sal} structures. 644 692 The computation is done in the \mdl{dtatsd} module. 645 \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 693 \item [{\np[=.false.]{ln_tsd_init}{ln\_tsd\_init}}] Initial values for T and S are set via 694 a user supplied \rou{usr\_def\_istate} routine contained in \mdl{userdef\_istate}. 646 695 The default version sets horizontally uniform T and profiles as used in the GYRE configuration 647 696 (see \autoref{sec:CFGS_gyre}). -
NEMO/trunk/doc/latex/NEMO/subfiles/chap_model_basics.tex
r11608 r11622 13 13 14 14 {\footnotesize 15 \begin{tabular x}{\textwidth}{l||X|X}16 Release & Author(s) & Modifications\\15 \begin{tabular}{l||l|l} 16 Release & Author(s) & Modifications \\ 17 17 \hline 18 {\em 4.0} & {\em Mike Bell } & {\em Update} \\19 {\em 3.6} & {\em Gurvan Madec } & {\em Minor changes} \\20 {\em <=3.4} & {\em Gurvan Madec and Steven Alderson} & {\em First version} \\21 \end{tabular x}18 {\em 4.0} & {\em Mike Bell } & {\em Review } \\ 19 {\em 3.6} & {\em Tim Graham and Gurvan Madec } & {\em Updates } \\ 20 {\em $\leq$ 3.4} & {\em Gurvan Madec and S\'{e}bastien Masson} & {\em First version} \\ 21 \end{tabular} 22 22 } 23 23 … … 1051 1051 \begin{equation} 1052 1052 \label{eq:MB_iso_slopes} 1053 r_1 = \frac{e_3}{e_1} 1054 r_2 = \frac{e_3}{e_2} 1053 r_1 = \frac{e_3}{e_1} \lt( \pd[\rho]{i} \rt) \lt( \pd[\rho]{k} \rt)^{-1} \quad 1054 r_2 = \frac{e_3}{e_2} \lt( \pd[\rho]{j} \rt) \lt( \pd[\rho]{k} \rt)^{-1} 1055 1055 \end{equation} 1056 1056 while in $s$-coordinates $\pd[]{k}$ is replaced by $\pd[]{s}$. -
NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex
r11599 r11622 13 13 14 14 {\footnotesize 15 \begin{tabular x}{\textwidth}{l||X|X}16 Release & Author(s) & Modifications\\15 \begin{tabular}{l||l|l} 16 Release & Author(s) & Modifications \\ 17 17 \hline 18 {\em 4.0} & {\em ...} & {\em ...} \\ 19 {\em 3.6} & {\em ...} & {\em ...} \\ 20 {\em 3.4} & {\em ...} & {\em ...} \\ 21 {\em <=3.4} & {\em ...} & {\em ...} 22 \end{tabularx} 18 {\em 4.0} & {\em J\'{e}r\^{o}me Chanut and Tim Graham} & {\em Review } \\ 19 {\em 3.6} & {\em Christian \'{E}th\'{e} } & {\em Update } \\ 20 {\em $\leq$ 3.4} & {\em Gurvan Madec } & {\em First version} \\ 21 \end{tabular} 23 22 } 24 23 … … 26 25 27 26 % Missing things: 28 % - daymod: definition of the time domain (nit000, nitend and the calendar) 29 30 \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here, 31 would help ==> to be added} 32 33 Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, 27 % - daymod: definition of the time domain (nit000, nitend and the calendar) 28 29 \gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which 30 could be referred to here, would help ==> to be added} 31 32 Having defined the continuous equations in \autoref{chap:MB}, 33 we need now to choose a time discretization, 34 34 a key feature of an ocean model as it exerts a strong influence on the structure of the computer code 35 35 (\ie\ on its flowchart). 36 In the present chapter, we provide a general description of the \NEMO\ 36 In the present chapter, we provide a general description of the \NEMO\ time stepping strategy and 37 37 the consequences for the order in which the equations are solved. 38 38 … … 47 47 \end{equation} 48 48 where $x$ stands for $u$, $v$, $T$ or $S$; 49 RHS is the Right-Hand-Side of the corresponding time evolution equation;49 RHS is the \textbf{R}ight-\textbf{H}and-\textbf{S}ide of the corresponding time evolution equation; 50 50 $\rdt$ is the time step; 51 51 and the superscripts indicate the time at which a quantity is evaluated. 52 Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated. 52 Each term of the RHS is evaluated at a specific time stepping depending on 53 the physics with which it is associated. 53 54 54 55 The choice of the time stepping used for this evaluation is discussed below as well as 55 56 the implications for starting or restarting a model simulation. 56 57 Note that the time stepping calculation is generally performed in a single operation. 57 With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in58 time for each term separately.58 With such a complex and nonlinear system of equations it would be dangerous to 59 let a prognostic variable evolve in time for each term separately. 59 60 60 61 The three level scheme requires three arrays for each prognostic variable. … … 62 63 The third array, although referred to as $x_a$ (after) in the code, 63 64 is usually not the variable at the after time step; 64 but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 65 The time stepping itself is performed once at each time step where implicit vertical diffusion is computed, \ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 65 but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) 66 prior to time-stepping the equation. 67 The time stepping itself is performed once at each time step where 68 implicit vertical diffusion is computed, 69 \ie\ in the \mdl{trazdf} and \mdl{dynzdf} modules. 66 70 67 71 %% ================================================================================================= … … 69 73 \label{sec:TD_leap_frog} 70 74 71 The time stepping used for processes other than diffusion is the well-known leapfrog scheme72 \citep{mesinger.arakawa_bk76}.75 The time stepping used for processes other than diffusion is 76 the well-known \textbf{L}eap\textbf{F}rog (LF) scheme \citep{mesinger.arakawa_bk76}. 73 77 This scheme is widely used for advection processes in low-viscosity fluids. 74 It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. 78 It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at 79 time step $t$, the now time step. 75 80 It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms, 76 81 but not for diffusion terms. 77 82 It is an efficient method that achieves second-order accuracy with 78 83 just one right hand side evaluation per time step. 79 Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by80 amplifying the oscillations.84 Moreover, it does not artificially damp linear oscillatory motion 85 nor does it produce instability by amplifying the oscillations. 81 86 These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme, 82 and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes. 87 and the unsuitability of leapfrog differencing for the representation of diffusion and 88 Rayleigh damping processes. 83 89 However, the scheme allows the coexistence of a numerical and a physical mode due to 84 90 its leading third order dispersive error. 85 91 In other words a divergence of odd and even time steps may occur. 86 To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter 87 (hereafter the LF-RA scheme). 88 This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 92 To prevent it, the leapfrog scheme is often used in association with 93 a \textbf{R}obert-\textbf{A}sselin time filter (hereafter the LF-RA scheme). 94 This filter, 95 first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72}, 89 96 is a kind of laplacian diffusion in time that mixes odd and even time steps: 90 97 \begin{equation} … … 99 106 However, the second order truncation error is proportional to $\gamma$, which is small compared to 1. 100 107 Therefore, the LF-RA is a quasi second order accurate scheme. 101 The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, 102 because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. 103 When used with the 2nd order space centred discretisation of the advection terms in 108 The LF-RA scheme is preferred to other time differencing schemes such as 109 predictor corrector or trapezoidal schemes, because the user has an explicit and simple control of 110 the magnitude of the time diffusion of the scheme. 111 When used with the 2$^nd$ order space centred discretisation of the advection terms in 104 112 the momentum and tracer equations, LF-RA avoids implicit numerical diffusion: 105 diffusion is set explicitly by the user through the Robert-Asselin 106 filter parameter andthe viscosity and diffusion coefficients.113 diffusion is set explicitly by the user through the Robert-Asselin filter parameter and 114 the viscosity and diffusion coefficients. 107 115 108 116 %% ================================================================================================= … … 110 118 \label{sec:TD_forward_imp} 111 119 112 The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. 120 The leapfrog differencing scheme is unsuitable for 121 the representation of diffusion and damping processes. 113 122 For a tendency $D_x$, representing a diffusion term or a restoring term to a tracer climatology 114 123 (when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : … … 119 128 120 129 This is diffusive in time and conditionally stable. 121 The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: 130 The conditions for stability of second and fourth order horizontal diffusion schemes are 131 \citep{griffies_bk04}: 122 132 \begin{equation} 123 133 \label{eq:TD_euler_stability} … … 128 138 \end{cases} 129 139 \end{equation} 130 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. 140 where $e$ is the smallest grid size in the two horizontal directions and 141 $A^h$ is the mixing coefficient. 131 142 The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient. 132 143 If it is not satisfied, even mildly, then the model soon becomes wildly unstable. 133 The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. 144 The instability can be removed by either reducing the length of the time steps or 145 reducing the mixing coefficient. 134 146 135 147 For the vertical diffusion terms, a forward time differencing scheme can be used, 136 but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a 137 backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: 148 but usually the numerical stability condition imposes a strong constraint on the time step. 149 To overcome the stability constraint, a backward (or implicit) time differencing scheme is used. 150 This scheme is unconditionally stable but diffusive and can be written as follows: 138 151 \begin{equation} 139 152 \label{eq:TD_imp} … … 145 158 %%gm 146 159 147 This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: 160 This scheme is rather time consuming since it requires a matrix inversion. 161 For example, the finite difference approximation of the temperature equation is: 148 162 \[ 149 163 % \label{eq:TD_imp_zdf} … … 165 179 \end{align*} 166 180 167 \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. 168 Moreover, 169 $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, 181 \autoref{eq:TD_imp_mat} is a linear system of equations with 182 an associated matrix which is tridiagonal. 183 Moreover, $c(k)$ and $d(k)$ are positive and 184 the diagonal term is greater than the sum of the two extra-diagonal terms, 170 185 therefore a special adaptation of the Gauss elimination procedure is used to find the solution 171 186 (see for example \citet{richtmyer.morton_bk67}). … … 175 190 \label{sec:TD_spg_ts} 176 191 177 The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 178 at \textit{now} time step. This refers to as the explicit free surface case in the code (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 179 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation 180 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface 181 (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc 182 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of 183 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 184 185 Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), the use of a split-explicit free surface is advantageous 186 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication 187 time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. 192 The leapfrog environment supports a centred in time computation of the surface pressure, 193 \ie\ evaluated at \textit{now} time step. 194 This refers to as the explicit free surface case in the code 195 (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}). 196 This choice however imposes a strong constraint on the time step which 197 should be small enough to resolve the propagation of external gravity waves. 198 As a matter of fact, one rather use in a realistic setup, 199 a split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}) in which 200 barotropic and baroclinic dynamical equations are solved separately with ad-hoc time steps. 201 The use of the time-splitting (in combination with non-linear free surface) imposes 202 some constraints on the design of the overall flowchart, 203 in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 204 205 Compared to the former use of the filtered free surface in \NEMO\ v3.6 (\citet{roullet.madec_JGR00}), 206 the use of a split-explicit free surface is advantageous on massively parallel computers. 207 Indeed, no global computations are anymore required by the elliptic solver which 208 saves a substantial amount of communication time. 209 Fast barotropic motions (such as tides) are also simulated with a better accuracy. 188 210 189 211 %\gmcomment{ 190 \begin{figure} [!t]212 \begin{figure} 191 213 \centering 192 214 \includegraphics[width=0.66\textwidth]{Fig_TimeStepping_flowchart_v4} 193 215 \caption[Leapfrog time stepping sequence with split-explicit free surface]{ 194 216 Sketch of the leapfrog time stepping sequence in \NEMO\ with split-explicit free surface. 195 The latter combined with non-linear free surface requires the dynamical tendency being196 updated prior tracers tendency to ensure conservation.217 The latter combined with non-linear free surface requires 218 the dynamical tendency being updated prior tracers tendency to ensure conservation. 197 219 Note the use of time integrated fluxes issued from the barotropic loop in 198 220 subsequent calculations of tracer advection and in the continuity equation. … … 203 225 204 226 %% ================================================================================================= 205 \section{Modified Leap frog -- Asselin filter scheme}227 \section{Modified LeapFrog -- Robert Asselin filter scheme (LF-RA)} 206 228 \label{sec:TD_mLF} 207 229 208 Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to 209 ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter. 230 Significant changes have been introduced by \cite{leclair.madec_OM09} in 231 the LF-RA scheme in order to ensure tracer conservation and to 232 allow the use of a much smaller value of the Asselin filter parameter. 210 233 The modifications affect both the forcing and filtering treatments in the LF-RA scheme. 211 234 212 In a classical LF-RA environment, the forcing term is centred in time,213 \ie\ it is time-stepped over a $2 \rdt$ period:235 In a classical LF-RA environment, 236 the forcing term is centred in time, \ie\ it is time-stepped over a $2 \rdt$ period: 214 237 $x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 215 and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. 238 and the time filter is given by \autoref{eq:TD_asselin} so that 239 $Q$ is redistributed over several time step. 216 240 In the modified LF-RA environment, these two formulations have been replaced by: 217 241 \begin{gather} … … 222 246 - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 223 247 \end{gather} 224 The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) 225 has a significant effect: 226 the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. 248 The change in the forcing formulation given by \autoref{eq:TD_forcing} 249 (see \autoref{fig:TD_MLF_forcing}) has a significant effect: 250 the forcing term no longer excites the divergence of odd and even time steps 251 \citep{leclair.madec_OM09}. 227 252 % forcing seen by the model.... 228 253 This property improves the LF-RA scheme in two aspects. 229 254 First, the LF-RA can now ensure the local and global conservation of tracers. 230 255 Indeed, time filtering is no longer required on the forcing part. 231 The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter232 (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}).256 The influence of the Asselin filter on the forcing is explicitly removed by 257 adding a new term in the filter (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 233 258 Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 234 259 the modified formulation becomes conservative \citep{leclair.madec_OM09}. 235 Second, the LF-RA becomes a truly quasi 260 Second, the LF-RA becomes a truly quasi-second order scheme. 236 261 Indeed, \autoref{eq:TD_forcing} used in combination with a careful treatment of static instability 237 262 (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) … … 245 270 even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term. 246 271 247 \begin{figure} [!t]272 \begin{figure} 248 273 \centering 249 274 \includegraphics[width=0.66\textwidth]{Fig_MLF_forcing} … … 271 296 \end{listing} 272 297 273 The first time step of this three level scheme when starting from initial conditions is a forward step274 (Euler time integration):298 The first time step of this three level scheme when starting from initial conditions is 299 a forward step (Euler time integration): 275 300 \[ 276 301 % \label{eq:TD_DOM_euler} 277 302 x^1 = x^0 + \rdt \ \text{RHS}^0 278 303 \] 279 This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but 304 This is done simply by keeping the leapfrog environment 305 (\ie\ the \autoref{eq:TD} three level time stepping) but 280 306 setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 281 307 using half the value of a leapfrog time step ($2 \rdt$). … … 286 312 running the model for $2N$ time steps in one go, 287 313 or by performing two consecutive experiments of $N$ steps with a restart. 288 This requires saving two time levels and many auxiliary data in the restart files in machine precision. 314 This requires saving two time levels and many auxiliary data in 315 the restart files in machine precision. 289 316 290 317 Note that the time step $\rdt$, is also saved in the restart file. 291 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step 292 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting 293 the namelist variable \np[=0]{nn_euler}{nn\_euler}. Other options to control the time integration of the model 294 are defined through the \nam{run}{run} namelist variables. 318 When restarting, if the time step has been changed, or 319 one of the prognostic variables at \textit{before} time step is missing, 320 an Euler time stepping scheme is imposed. 321 A forward initial step can still be enforced by the user by 322 setting the namelist variable \np[=0]{nn_euler}{nn\_euler}. 323 Other options to control the time integration of the model are defined through 324 the \nam{run}{run} namelist variables. 325 295 326 \gmcomment{ 296 327 add here how to force the restart to contain only one time step for operational purposes … … 298 329 add also the idea of writing several restart for seasonal forecast : how is it done ? 299 330 300 verify that all namelist para rmeters are truly described331 verify that all namelist parameters are truly described 301 332 302 333 a word on the check of restart ..... … … 309 340 \label{subsec:TD_time} 310 341 311 Options are defined through the 342 Options are defined through the\nam{dom}{dom} namelist variables. 312 343 \colorbox{yellow}{add here a few word on nit000 and nitend} 313 344 314 345 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj} 315 346 316 add a description of daymod, and the model cal andar (leap-year and co)317 318 } 347 add a description of daymod, and the model calendar (leap-year and co) 348 349 } %% end add 319 350 320 351 \gmcomment{ % add implicit in vvl case and Crant-Nicholson scheme
Note: See TracChangeset
for help on using the changeset viewer.