source: NEMO/trunk/doc/latex/NEMO/subfiles/chap_time_domain.tex @ 11567

Last change on this file since 11567 was 11567, checked in by nicolasmartin, 13 months ago

Style improvments

File size: 20.6 KB
Line 
1\documentclass[../main/NEMO_manual]{subfiles}
2
3\begin{document}
4
5% ================================================================
6% Chapter 2 ——— Time Domain (step.F90)
7% ================================================================
8\chapter{Time Domain}
9\label{chap:TD}
10\chaptertoc
11
12% Missing things:
13%  - daymod: definition of the time domain (nit000, nitend and the calendar)
14
15\gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here,
16  would help  ==> to be added}
17%%%%
18
19\newpage
20
21Having defined the continuous equations in \autoref{chap:MB}, we need now to choose a time discretization,
22a key feature of an ocean model as it exerts a strong influence on the structure of the computer code
23(\ie\ on its flowchart).
24In the present chapter, we provide a general description of the \NEMO\  time stepping strategy and
25the consequences for the order in which the equations are solved.
26
27% ================================================================
28% Time Discretisation
29% ================================================================
30\section{Time stepping environment}
31\label{sec:TD_environment}
32
33The time stepping used in \NEMO\ is a three level scheme that can be represented as follows:
34\begin{equation}
35  \label{eq:TD}
36  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t - \rdt, \, t, \, t + \rdt}
37\end{equation}
38where $x$ stands for $u$, $v$, $T$ or $S$;
39RHS is the Right-Hand-Side of the corresponding time evolution equation;
40$\rdt$ is the time step;
41and the superscripts indicate the time at which a quantity is evaluated.
42Each term of the RHS is evaluated at a specific time stepping depending on the physics with which it is associated.
43
44The choice of the time stepping used for this evaluation is discussed below as well as
45the implications for starting or restarting a model simulation.
46Note that the time stepping calculation is generally performed in a single operation.
47With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in
48time for each term separately.
49
50The three level scheme requires three arrays for each prognostic variable.
51For each variable $x$ there is $x_b$ (before), $x_n$ (now) and $x_a$.
52The third array, although referred to as $x_a$ (after) in the code,
53is usually not the variable at the after time step;
54but rather it is used to store the time derivative (RHS in \autoref{eq:TD}) prior to time-stepping the equation.
55The 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
57% -------------------------------------------------------------------------------------------------------------
58%        Non-Diffusive Part---Leapfrog Scheme
59% -------------------------------------------------------------------------------------------------------------
60\section{Non-diffusive part --- Leapfrog scheme}
61\label{sec:TD_leap_frog}
62
63The time stepping used for processes other than diffusion is the well-known leapfrog scheme
64\citep{mesinger.arakawa_bk76}.
65This scheme is widely used for advection processes in low-viscosity fluids.
66It is a time centred scheme, \ie\ the RHS in \autoref{eq:TD} is evaluated at time step $t$, the now time step.
67It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms,
68but not for diffusion terms.
69It is an efficient method that achieves second-order accuracy with
70just one right hand side evaluation per time step.
71Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by
72amplifying the oscillations.
73These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme,
74and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes.
75However, the scheme allows the coexistence of a numerical and a physical mode due to
76its leading third order dispersive error.
77In other words a divergence of odd and even time steps may occur.
78To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter
79(hereafter the LF-RA scheme).
80This filter, first designed by \citet{robert_JMSJ66} and more comprehensively studied by \citet{asselin_MWR72},
81is a kind of laplacian diffusion in time that mixes odd and even time steps:
82\begin{equation}
83  \label{eq:TD_asselin}
84  x_F^t = x^t + \gamma \, \lt[ x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt]
85\end{equation}
86where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient.
87$\gamma$ is initialized as \np{rn\_atfp} (namelist parameter).
88Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:TD_mLF}),
89causing only a weak dissipation of high frequency motions (\citep{farge-coulombier_phd87}).
90The addition of a time filter degrades the accuracy of the calculation from second to first order.
91However, the second order truncation error is proportional to $\gamma$, which is small compared to 1.
92Therefore, the LF-RA is a quasi second order accurate scheme.
93The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes,
94because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme.
95When used with the 2nd order space centred discretisation of the advection terms in
96the momentum and tracer equations, LF-RA avoids implicit numerical diffusion:
97diffusion is set explicitly by the user through the Robert-Asselin
98filter parameter and the viscosity and diffusion coefficients.
99
100% -------------------------------------------------------------------------------------------------------------
101%        Diffusive Part---Forward or Backward Scheme
102% -------------------------------------------------------------------------------------------------------------
103\section{Diffusive part --- Forward or backward scheme}
104\label{sec:TD_forward_imp}
105
106The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes.
107For a tendency $D_x$, representing a diffusion term or a restoring term to a tracer climatology
108(when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used :
109\[
110  %\label{eq:TD_euler}
111  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ D_x^{t - \rdt}
112\]
113
114This is diffusive in time and conditionally stable.
115The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{griffies_bk04}:
116\begin{equation}
117  \label{eq:TD_euler_stability}
118  A^h <
119  \begin{cases}
120    \frac{e^2}{ 8 \, \rdt} & \text{laplacian diffusion} \\
121    \frac{e^4}{64 \, \rdt} & \text{bilaplacian diffusion}
122  \end{cases}
123\end{equation}
124where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient.
125The linear constraint \autoref{eq:TD_euler_stability} is a necessary condition, but not sufficient.
126If it is not satisfied, even mildly, then the model soon becomes wildly unstable.
127The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient.
128
129For the vertical diffusion terms, a forward time differencing scheme can be used,
130but usually the numerical stability condition imposes a strong constraint on the time step. To overcome the stability constraint, a
131backward (or implicit) time differencing scheme is used. This scheme is unconditionally stable but diffusive and can be written as follows:
132\begin{equation}
133  \label{eq:TD_imp}
134  x^{t + \rdt} = x^{t - \rdt} + 2 \, \rdt \ \text{RHS}_x^{t + \rdt}
135\end{equation}
136
137%%gm
138%%gm   UPDATE the next paragraphs with time varying thickness ...
139%%gm
140
141This scheme is rather time consuming since it requires a matrix inversion. For example, the finite difference approximation of the temperature equation is:
142\[
143  % \label{eq:TD_imp_zdf}
144  \frac{T(k)^{t + 1} - T(k)^{t - 1}}{2 \; \rdt}
145  \equiv
146  \text{RHS} + \frac{1}{e_{3t}} \delta_k \lt[ \frac{A_w^{vT}}{e_{3w} } \delta_{k + 1/2} \lt[ T^{t + 1} \rt] \rt]
147\]
148where RHS is the right hand side of the equation except for the vertical diffusion term.
149We rewrite \autoref{eq:TD_imp} as:
150\begin{equation}
151  \label{eq:TD_imp_mat}
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\end{equation}
154where
155\begin{align*}
156  c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k)     \\
157  d(k) &= e_{3t}   (k)       \, / \, (2 \rdt) + c_k + c_{k + 1}    \\
158  b(k) &= e_{3t}   (k) \; \lt( T^{t - 1}(k) \, / \, (2 \rdt) + \text{RHS} \rt)
159\end{align*}
160
161\autoref{eq:TD_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal.
162Moreover,
163$c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms,
164therefore a special adaptation of the Gauss elimination procedure is used to find the solution
165(see for example \citet{richtmyer.morton_bk67}).
166
167% -------------------------------------------------------------------------------------------------------------
168%        Surface Pressure gradient
169% -------------------------------------------------------------------------------------------------------------
170\section{Surface pressure gradient}
171\label{sec:TD_spg_ts}
172
173The leapfrog environment supports a centred in time computation of the surface pressure, \ie\ evaluated
174at \textit{now} time step. This refers to as the explicit free surface case in the code (\np{ln\_dynspg\_exp}\forcode{=.true.}).
175This choice however imposes a strong constraint on the time step which should be small enough to resolve the propagation
176of 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
178time steps. The use of the time-splitting (in combination with non-linear free surface) imposes some constraints on the design of
179the overall flowchart, in particular to ensure exact tracer conservation (see \autoref{fig:TD_TimeStep_flowchart}).
180
181Compared 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
182on massively parallel computers. Indeed, no global computations are anymore required by the elliptic solver which saves a substantial amount of communication
183time. Fast barotropic motions (such as tides) are also simulated with a better accuracy.
184
185%\gmcomment{
186%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
187\begin{figure}[!t]
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}
198\end{figure}
199%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
200%}
201
202% -------------------------------------------------------------------------------------------------------------
203%        The Modified Leapfrog -- Asselin Filter scheme
204% -------------------------------------------------------------------------------------------------------------
205\section{Modified Leapfrog -- Asselin filter scheme}
206\label{sec:TD_mLF}
207
208Significant changes have been introduced by \cite{leclair.madec_OM09} in the LF-RA scheme in order to
209ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter.
210The modifications affect both the forcing and filtering treatments in the LF-RA scheme.
211
212In a classical LF-RA environment, the forcing term is centred in time,
213\ie\ it is time-stepped over a $2 \rdt$ period:
214$x^t = x^t + 2 \rdt Q^t$ where $Q$ is the forcing applied to $x$,
215and the time filter is given by \autoref{eq:TD_asselin} so that $Q$ is redistributed over several time step.
216In the modified LF-RA environment, these two formulations have been replaced by:
217\begin{gather}
218  \label{eq:TD_forcing}
219  x^{t + \rdt} = x^{t - \rdt} + \rdt \lt( Q^{t - \rdt / 2} + Q^{t + \rdt / 2} \rt\\
220  \label{eq:TD_RA}
221  x_F^t       = x^t + \gamma \, \lt( x_F^{t - \rdt} - 2 x^t + x^{t + \rdt} \rt)
222                    - \gamma \, \rdt \, \lt( Q^{t + \rdt / 2} - Q^{t - \rdt / 2} \rt)
223\end{gather}
224The change in the forcing formulation given by \autoref{eq:TD_forcing} (see \autoref{fig:TD_MLF_forcing})
225has a significant effect:
226the forcing term no longer excites the divergence of odd and even time steps \citep{leclair.madec_OM09}.
227% forcing seen by the model....
228This property improves the LF-RA scheme in two aspects.
229First, the LF-RA can now ensure the local and global conservation of tracers.
230Indeed, time filtering is no longer required on the forcing part.
231The influence of the Asselin filter on the forcing is explicitly removed by adding a new term in the filter
232(last term in \autoref{eq:TD_RA} compared to \autoref{eq:TD_asselin}).
233Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme,
234the modified formulation becomes conservative \citep{leclair.madec_OM09}.
235Second, the LF-RA becomes a truly quasi -second order scheme.
236Indeed, \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})
238(the two other main sources of time step divergence),
239allows a reduction by two orders of magnitude of the Asselin filter parameter.
240
241Note that the forcing is now provided at the middle of a time step:
242$Q^{t + \rdt / 2}$ is the forcing applied over the $[t,t + \rdt]$ time interval.
243This and the change in the time filter, \autoref{eq:TD_RA},
244allows for an exact evaluation of the contribution due to the forcing term between any two time steps,
245even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term.
246
247%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
248\begin{figure}[!t]
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}
262\end{figure}
263%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
264
265% -------------------------------------------------------------------------------------------------------------
266%        Start/Restart strategy
267% -------------------------------------------------------------------------------------------------------------
268\section{Start/Restart strategy}
269\label{sec:TD_rst}
270
271%--------------------------------------------namrun-------------------------------------------
272\begin{listing}
273  \nlst{namrun}
274  \caption{\forcode{&namrun}}
275  \label{lst:namrun}
276\end{listing}
277%--------------------------------------------------------------------------------------------------------------
278
279The first time step of this three level scheme when starting from initial conditions is a forward step
280(Euler time integration):
281\[
282  % \label{eq:TD_DOM_euler}
283  x^1 = x^0 + \rdt \ \text{RHS}^0
284\]
285This is done simply by keeping the leapfrog environment (\ie\ the \autoref{eq:TD} three level time stepping) but
286setting all $x^0$ (\textit{before}) and $x^1$ (\textit{now}) fields equal at the first time step and
287using half the value of a leapfrog time step ($2 \rdt$).
288
289It is also possible to restart from a previous computation, by using a restart file.
290The restart strategy is designed to ensure perfect restartability of the code:
291the user should obtain the same results to machine precision either by
292running the model for $2N$ time steps in one go,
293or by performing two consecutive experiments of $N$ steps with a restart.
294This requires saving two time levels and many auxiliary data in the restart files in machine precision.
295
296Note that the time step $\rdt$, is also saved in the restart file.
297When restarting, if the time step has been changed, or one of the prognostic variables at \textit{before} time step
298is missing, an Euler time stepping scheme is imposed. A forward initial step can still be enforced by the user by setting
299the namelist variable \np{nn\_euler}\forcode{=0}. Other options to control the time integration of the model
300are defined through the  \nam{run} namelist variables.
301%%%
302\gmcomment{
303add here how to force the restart to contain only one time step for operational purposes
304
305add also the idea of writing several restart for seasonal forecast : how is it done ?
306
307verify that all namelist pararmeters are truly described
308
309a word on the check of restart  .....
310}
311%%%
312
313\gmcomment{       % add a subsection here
314
315%-------------------------------------------------------------------------------------------------------------
316%        Time Domain
317% -------------------------------------------------------------------------------------------------------------
318\subsection{Time domain}
319\label{subsec:TD_time}
320%--------------------------------------------namrun-------------------------------------------
321
322%--------------------------------------------------------------------------------------------------------------
323
324Options are defined through the  \nam{dom} namelist variables.
325 \colorbox{yellow}{add here a few word on nit000 and nitend}
326
327 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj}
328
329add a description of daymod, and the model calandar (leap-year and co)
330
331}        %% end add
332
333
334
335%%
336\gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme
337
338Implicit time stepping in case of variable volume thickness.
339
340Tracer case (NB for momentum in vector invariant form take care!)
341
342\begin{flalign*}
343  &\frac{\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}}{2\rdt}
344  \equiv \text{RHS}+ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
345  \rt]      \\
346  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
347  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \lt[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k + 1/2} \lt[ {T^{t+1}} \rt]}
348  \rt]      \\
349  &\lt( e_{3t}\,T \rt)_k^{t+1}-\lt( e_{3t}\,T \rt)_k^{t-1}
350  \equiv 2\rdt \ \text{RHS}
351  + 2\rdt \ \lt\{ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2} [ T_{k +1}^{t+1} - T_k      ^{t+1} ]
352    - \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} [ T_k       ^{t+1} - T_{k -1}^{t+1} ]  \rt\}     \\
353  &\\
354  &\lt( e_{3t}\,T \rt)_k^{t+1}
355  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}                  T_{k +1}^{t+1}
356  + {2\rdt} \ \lt\{  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
357    +  \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}     \rt\}   T_{k    }^{t+1}
358  -  {2\rdt} \           \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                  T_{k -1}^{t+1}      \\
359  &\equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
360  %
361\end{flalign*}
362\begin{flalign*}
363  \allowdisplaybreaks
364  \intertext{ Tracer case }
365  %
366  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
367  \qquad \qquad \qquad  \qquad  T_{k +1}^{t+1}   \\
368  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k + 1/2}
369  +   \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
370  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ {2\rdt} \                          \lt[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \rt]_{k - 1/2}                          \quad \ \ T_{k -1}^{t+1}
371  \ \equiv \ \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
372  %
373\end{flalign*}
374\begin{flalign*}
375  \allowdisplaybreaks
376  \intertext{ Tracer content case }
377  %
378  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_{k +1}^{t+1}}  && \  \lt( e_{3t}\,T \rt)_{k +1}^{t+1}   &\\
379  & + {2\rdt} \ \lt[ 1  \rt.+ & \frac{(A_w^{vt})_{k + 1/2}} {(e_{3w})_{k + 1/2}^{t+1}\;(e_{3t})_k^{t+1}}
380  + & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \lt\rt& \lt( e_{3t}\,T \rt)_{k   }^{t+1}  &\\
381  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k - 1/2}} {(e_{3w})_{k - 1/2}^{t+1}\;(e_{3t})_{k -1}^{t+1}}     &\  \lt( e_{3t}\,T \rt)_{k -1}^{t+1}
382  \equiv \lt( e_{3t}\,T \rt)_k^{t-1} + {2\rdt} \ \text{RHS}  &
383\end{flalign*}
384
385%%
386}
387
388\biblio
389
390\pindex
391
392\end{document}
Note: See TracBrowser for help on using the repository browser.