New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
chap_time_domain.tex in NEMO/branches/2019/dev_r11351_fldread_with_XIOS/doc/latex/NEMO/subfiles – NEMO

source: NEMO/branches/2019/dev_r11351_fldread_with_XIOS/doc/latex/NEMO/subfiles/chap_time_domain.tex @ 13463

Last change on this file since 13463 was 13463, checked in by andmirek, 4 years ago

Ticket #2195:update to trunk 13461

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