Changeset 11573 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_time_domain.tex
- Timestamp:
- 2019-09-19T11:18:03+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/doc/latex/NEMO/subfiles/chap_time_domain.tex
r11184 r11573 6 6 % Chapter 2 ——— Time Domain (step.F90) 7 7 % ================================================================ 8 \chapter{Time Domain (STP)}9 \label{chap: STP}10 \ minitoc8 \chapter{Time Domain} 9 \label{chap:TD} 10 \chaptertoc 11 11 12 12 % Missing things: … … 19 19 \newpage 20 20 21 Having defined the continuous equations in \autoref{chap: PE}, we need now to choose a time discretization,21 Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization, 22 22 a key feature of an ocean model as it exerts a strong influence on the structure of the computer code 23 (\ie on its flowchart).24 In the present chapter, we provide a general description of the \NEMO time stepping strategy and23 (\ie\ on its flowchart). 24 In the present chapter, we provide a general description of the \NEMO\ time stepping strategy and 25 25 the consequences for the order in which the equations are solved. 26 26 … … 29 29 % ================================================================ 30 30 \section{Time stepping environment} 31 \label{sec: STP_environment}32 33 The time stepping used in \NEMO is a three level scheme that can be represented as follows:34 \begin{equation} 35 \label{eq: STP}31 \label{sec:TD_environment} 32 33 The time stepping used in \NEMO\ is a three level scheme that can be represented as follows: 34 \begin{equation} 35 \label{eq:TD} 36 36 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t - \rdt, \, t, \, t + \rdt} 37 \end{equation} 37 \end{equation} 38 38 where $x$ stands for $u$, $v$, $T$ or $S$; 39 39 RHS is the Right-Hand-Side of the corresponding time evolution equation; … … 52 52 The third array, although referred to as $x_a$ (after) in the code, 53 53 is usually not the variable at the after time step; 54 but rather it is used to store the time derivative (RHS in \autoref{eq: STP}) prior to time-stepping the equation.55 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.54 but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation. 55 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. 56 56 57 57 % ------------------------------------------------------------------------------------------------------------- … … 59 59 % ------------------------------------------------------------------------------------------------------------- 60 60 \section{Non-diffusive part --- Leapfrog scheme} 61 \label{sec: STP_leap_frog}61 \label{sec:TD_leap_frog} 62 62 63 63 The time stepping used for processes other than diffusion is the well-known leapfrog scheme 64 64 \citep{mesinger.arakawa_bk76}. 65 65 This scheme is widely used for advection processes in low-viscosity fluids. 66 It is a time centred scheme, \ie the RHS in \autoref{eq:STP} is evaluated at time step $t$, the now time step.66 It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step. 67 67 It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms, 68 68 but not for diffusion terms. … … 81 81 is a kind of laplacian diffusion in time that mixes odd and even time steps: 82 82 \begin{equation} 83 \label{eq: STP_asselin}83 \label{eq:TD_asselin} 84 84 x_F^t = x^t + \gamma \, \lt[ x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt] 85 85 \end{equation} 86 86 where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient. 87 87 $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 88 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec: STP_mLF}),88 Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:TD_mLF}), 89 89 causing only a weak dissipation of high frequency motions (\citep{farge-coulombier_phd87}). 90 90 The addition of a time filter degrades the accuracy of the calculation from second to first order. … … 92 92 Therefore, the LF-RA is a quasi second order accurate scheme. 93 93 The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes, 94 because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. 94 because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme. 95 95 When used with the 2nd order space centred discretisation of the advection terms in 96 96 the momentum and tracer equations, LF-RA avoids implicit numerical diffusion: 97 diffusion is set explicitly by the user through the Robert-Asselin 97 diffusion is set explicitly by the user through the Robert-Asselin 98 98 filter parameter and the viscosity and diffusion coefficients. 99 99 … … 102 102 % ------------------------------------------------------------------------------------------------------------- 103 103 \section{Diffusive part --- Forward or backward scheme} 104 \label{sec: STP_forward_imp}104 \label{sec:TD_forward_imp} 105 105 106 106 The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes. … … 108 108 (when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used : 109 109 \[ 110 %\label{eq: STP_euler}110 %\label{eq:TD_euler} 111 111 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ D_x^{t - \rdt} 112 112 \] … … 115 115 The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}: 116 116 \begin{equation} 117 \label{eq: STP_euler_stability}117 \label{eq:TD_euler_stability} 118 118 A^h < 119 119 \begin{cases} … … 123 123 \end{equation} 124 124 where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient. 125 The linear constraint \autoref{eq: STP_euler_stability} is a necessary condition, but not sufficient.125 The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient. 126 126 If it is not satisfied, even mildly, then the model soon becomes wildly unstable. 127 127 The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient. 128 128 129 129 For the vertical diffusion terms, a forward time differencing scheme can be used, 130 but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a 130 but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a 131 131 backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows: 132 132 \begin{equation} 133 \label{eq: STP_imp}133 \label{eq:TD_imp} 134 134 x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t + \rdt} 135 135 \end{equation} … … 141 141 This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is: 142 142 \[ 143 % \label{eq: STP_imp_zdf}143 % \label{eq:TD_imp_zdf} 144 144 \frac{T(k)^{t + 1} - T(k)^{t - 1}}{2 \; \rdt} 145 145 \equiv … … 147 147 \] 148 148 where RHS is the right hand side of the equation except for the vertical diffusion term. 149 We rewrite \autoref{eq: STP_imp} as:150 \begin{equation} 151 \label{eq: STP_imp_mat}149 We rewrite \autoref{eq:TD_imp} as: 150 \begin{equation} 151 \label{eq:TD_imp_mat} 152 152 -c(k + 1) \; T^{t + 1}(k + 1) + d(k) \; T^{t + 1}(k) - \; c(k) \; T^{t + 1}(k - 1) \equiv b(k) 153 153 \end{equation} 154 where 155 \begin{align*} 154 where 155 \begin{align*} 156 156 c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k) \\ 157 157 d(k) &= e_{3t} (k) \, / \, (2 \rdt) + c_k + c_{k + 1} \\ … … 159 159 \end{align*} 160 160 161 \autoref{eq: STP_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal.161 \autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal. 162 162 Moreover, 163 163 $c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms, … … 169 169 % ------------------------------------------------------------------------------------------------------------- 170 170 \section{Surface pressure gradient} 171 \label{sec: STP_spg_ts}172 173 The leapfrog environment supports a centred in time computation of the surface pressure, \ie evaluated174 at \textit{now} time step. This refers to as the explicit free surface case in the code (\np{ln\_dynspg\_exp}\forcode{ = .true.}).175 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation 176 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface 177 (\np{ln\_dynspg\_ts}\forcode{ = .true.}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc178 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of 179 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:T imeStep_flowchart}).180 181 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 advantageous182 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication 183 time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. 184 185 %\gmcomment{ 171 \label{sec:TD_spg_ts} 172 173 The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated 174 at \textit{now} time step. This refers to as the explicit free surface case in the code (\np{ln\_dynspg\_exp}\forcode{=.true.}). 175 This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation 176 of external gravity waves. As a matter of fact, one rather use in a realistic setup, a split-explicit free surface 177 (\np{ln\_dynspg\_ts}\forcode{=.true.}) in which barotropic and baroclinic dynamical equations are solved separately with ad-hoc 178 time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of 179 the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}). 180 181 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 182 on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication 183 time. Fast barotropic motions (such as tides) are also simulated with a better accuracy. 184 185 %\gmcomment{ 186 186 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 187 187 \begin{figure}[!t] 188 \begin{center} 189 \includegraphics[width=\textwidth]{Fig_TimeStepping_flowchart_v4} 190 \caption{ 191 \protect\label{fig:TimeStep_flowchart} 192 Sketch of the leapfrog time stepping sequence in \NEMO with split-explicit free surface. The latter combined 193 with non-linear free surface requires the dynamical tendency being updated prior tracers tendency to ensure 194 conservation. Note the use of time integrated fluxes issued from the barotropic loop in subsequent calculations 195 of tracer advection and in the continuity equation. Details about the time-splitting scheme can be found 196 in \autoref{subsec:DYN_spg_ts}. 197 } 198 \end{center} 188 \centering 189 \includegraphics[width=0.66\textwidth]{Fig_TimeStepping_flowchart_v4} 190 \caption[Leapfrog time stepping sequence with split-explicit free surface]{ 191 Sketch of the leapfrog time stepping sequence in \NEMO\ with split-explicit free surface. 192 The latter combined with non-linear free surface requires the dynamical tendency being 193 updated prior tracers tendency to ensure conservation. 194 Note the use of time integrated fluxes issued from the barotropic loop in 195 subsequent calculations of tracer advection and in the continuity equation. 196 Details about the time-splitting scheme can be found in \autoref{subsec:DYN_spg_ts}.} 197 \label{fig:TD_TimeStep_flowchart} 199 198 \end{figure} 200 199 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 205 204 % ------------------------------------------------------------------------------------------------------------- 206 205 \section{Modified Leapfrog -- Asselin filter scheme} 207 \label{sec: STP_mLF}206 \label{sec:TD_mLF} 208 207 209 208 Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to … … 212 211 213 212 In a classical LF-RA environment, the forcing term is centred in time, 214 \ie it is time-stepped over a $2 \rdt$ period:213 \ie\ it is time-stepped over a $2 \rdt$ period: 215 214 $x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$, 216 and the time filter is given by \autoref{eq: STP_asselin} so that $Q$ is redistributed over several time step.215 and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step. 217 216 In the modified LF-RA environment, these two formulations have been replaced by: 218 217 \begin{gather} 219 \label{eq: STP_forcing}218 \label{eq:TD_forcing} 220 219 x^{t + \rdt} = x^{t - \rdt} + \rdt \lt( Q^{t - \rdt / 2} + Q^{t + \rdt / 2} \rt) \\ 221 \label{eq: STP_RA}220 \label{eq:TD_RA} 222 221 x_F^t = x^t + \gamma \, \lt( x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt) 223 222 - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt) 224 223 \end{gather} 225 The change in the forcing formulation given by \autoref{eq: STP_forcing} (see \autoref{fig:MLF_forcing})224 The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing}) 226 225 has a significant effect: 227 226 the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}. … … 231 230 Indeed, time filtering is no longer required on the forcing part. 232 231 The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter 233 (last term in \autoref{eq: STP_RA} compared to \autoref{eq:STP_asselin}).232 (last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}). 234 233 Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme, 235 234 the modified formulation becomes conservative \citep{leclair.madec_OM09}. 236 235 Second, the LF-RA becomes a truly quasi -second order scheme. 237 Indeed, \autoref{eq: STP_forcing} used in combination with a careful treatment of static instability238 (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) 236 Indeed, \autoref{eq:TD_forcing} used in combination with a careful treatment of static instability 237 (\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}) 239 238 (the two other main sources of time step divergence), 240 239 allows a reduction by two orders of magnitude of the Asselin filter parameter. … … 242 241 Note that the forcing is now provided at the middle of a time step: 243 242 $Q^{t + \rdt / 2}$ is the forcing applied over the $[t,t + \rdt]$ time interval. 244 This and the change in the time filter, \autoref{eq: STP_RA},243 This and the change in the time filter, \autoref{eq:TD_RA}, 245 244 allows for an exact evaluation of the contribution due to the forcing term between any two time steps, 246 245 even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term. … … 248 247 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 249 248 \begin{figure}[!t] 250 \ begin{center}251 \includegraphics[width=\textwidth]{Fig_MLF_forcing}252 \caption{253 \protect\label{fig:MLF_forcing}254 Illustration of forcing integration methods.255 (top) ''Traditional'' formulation:256 the forcing is defined at the same time as the variable to which it is applied257 (integer value of the time step index) and it is applied over a $2 \rdt$ period.258 (bottom) modified formulation:259 the forcing is defined in the middle of the time(integer and a half value of the time step index) and260 the mean of two successive forcing values ($n - 1 / 2$, $n + 1 / 2$) is applied over a $2 \rdt$ period.261 }262 \ end{center}249 \centering 250 \includegraphics[width=0.66\textwidth]{Fig_MLF_forcing} 251 \caption[Forcing integration methods for modified leapfrog (top and bottom)]{ 252 Illustration of forcing integration methods. 253 (top) ''Traditional'' formulation: 254 the forcing is defined at the same time as the variable to which it is applied 255 (integer value of the time step index) and it is applied over a $2 \rdt$ period. 256 (bottom) modified formulation: 257 the forcing is defined in the middle of the time 258 (integer and a half value of the time step index) and 259 the mean of two successive forcing values ($n - 1 / 2$, $n + 1 / 2$) is applied over 260 a $2 \rdt$ period.} 261 \label{fig:TD_MLF_forcing} 263 262 \end{figure} 264 263 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 268 267 % ------------------------------------------------------------------------------------------------------------- 269 268 \section{Start/Restart strategy} 270 \label{sec: STP_rst}269 \label{sec:TD_rst} 271 270 272 271 %--------------------------------------------namrun------------------------------------------- 273 \nlst{namrun} 272 \begin{listing} 273 \nlst{namrun} 274 \caption{\forcode{&namrun}} 275 \label{lst:namrun} 276 \end{listing} 274 277 %-------------------------------------------------------------------------------------------------------------- 275 278 … … 277 280 (Euler time integration): 278 281 \[ 279 % \label{eq: DOM_euler}282 % \label{eq:TD_DOM_euler} 280 283 x^1 = x^0 + \rdt \ \text{RHS}^0 281 284 \] 282 This is done simply by keeping the leapfrog environment (\ie the \autoref{eq:STP} three level time stepping) but285 This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but 283 286 setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and 284 using half the value of a leapfrog time step ($2 \rdt$). 287 using half the value of a leapfrog time step ($2 \rdt$). 285 288 286 289 It is also possible to restart from a previous computation, by using a restart file. … … 292 295 293 296 Note that the time step $\rdt$, is also saved in the restart file. 294 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step 295 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting 296 the namelist variable \np{nn\_euler}\forcode{=0}. Other options to control the time integration of the model 297 are defined through the \n gn{namrun} namelist variables.297 When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step 298 is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting 299 the namelist variable \np{nn\_euler}\forcode{=0}. Other options to control the time integration of the model 300 are defined through the \nam{run} namelist variables. 298 301 %%% 299 302 \gmcomment{ … … 302 305 add also the idea of writing several restart for seasonal forecast : how is it done ? 303 306 304 verify that all namelist pararmeters are truly described 307 verify that all namelist pararmeters are truly described 305 308 306 309 a word on the check of restart ..... … … 308 311 %%% 309 312 310 \gmcomment{ % add a subsection here 313 \gmcomment{ % add a subsection here 311 314 312 315 %------------------------------------------------------------------------------------------------------------- … … 314 317 % ------------------------------------------------------------------------------------------------------------- 315 318 \subsection{Time domain} 316 \label{subsec: STP_time}319 \label{subsec:TD_time} 317 320 %--------------------------------------------namrun------------------------------------------- 318 321 319 \nlst{namdom}320 322 %-------------------------------------------------------------------------------------------------------------- 321 323 322 Options are defined through the \n gn{namdom} namelist variables.324 Options are defined through the \nam{dom} namelist variables. 323 325 \colorbox{yellow}{add here a few word on nit000 and nitend} 324 326 … … 332 334 333 335 %% 334 \gmcomment{ % add implicit in vvl case and Crant-Nicholson scheme 336 \gmcomment{ % add implicit in vvl case and Crant-Nicholson scheme 335 337 336 338 Implicit time stepping in case of variable volume thickness.
Note: See TracChangeset
for help on using the changeset viewer.
