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

Last change on this file since 10442 was 10442, checked in by nicolasmartin, 2 years ago

Front page edition, cleaning in custom LaTeX commands and add index for single subfile compilation

  • Use \thanks storing cmd to refer to the ST members list for 2018 in an footnote on the cover page
  • NEMO and Fortran in small capitals
  • Removing of unused or underused custom cmds, move local cmds to their respective .tex file
  • Addition of new ones (\zstar, \ztilde, \sstar, \stilde, \ie, \eg, \fortran, \fninety)
  • Fonts for indexed items: italic font for files (modules and .nc files), preformat for code (CPP keys, routines names and namelists content)
File size: 20.8 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 (STP) }
9\label{chap:STP}
10
11\minitoc
12
13% Missing things:
14%  - daymod: definition of the time domain (nit000, nitend andd the calendar)
15
16\gmcomment{STEVEN :maybe a picture of the directory structure in the introduction which could be referred to here,
17  would help  ==> to be added}
18%%%%
19
20\newpage
21
22Having defined the continuous equations in \autoref{chap:PE}, we need now to choose a time discretization,
23a key feature of an ocean model as it exerts a strong influence on the structure of the computer code
24(\ie on its flowchart).
25In the present chapter, we provide a general description of the \NEMO time stepping strategy and
26the consequences for the order in which the equations are solved.
27
28% ================================================================
29% Time Discretisation
30% ================================================================
31\section{Time stepping environment}
32\label{sec:STP_environment}
33
34The time stepping used in \NEMO is a three level scheme that can be represented as follows:
35\begin{equation}
36  \label{eq:STP}
37  x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t-\rdt,\,t,\,t+\rdt}
38\end{equation} 
39where $x$ stands for $u$, $v$, $T$ or $S$;
40RHS is the Right-Hand-Side of the corresponding time evolution equation;
41$\rdt$ is the time step;
42and the superscripts indicate the time at which a quantity is evaluated.
43Each term of the RHS is evaluated at a specific time step depending on the physics with which it is associated.
44
45The choice of the time step used for this evaluation is discussed below as well as
46the implications for starting or restarting a model simulation.
47Note that the time stepping calculation is generally performed in a single operation.
48With such a complex and nonlinear system of equations it would be dangerous to let a prognostic variable evolve in
49time for each term separately.
50
51The three level scheme requires three arrays for each prognostic variable.
52For each variable $x$ there is $x_b$ (before), $x_n$ (now) and $x_a$.
53The third array, although referred to as $x_a$ (after) in the code,
54is usually not the variable at the after time step;
55but rather it is used to store the time derivative (RHS in \autoref{eq:STP}) prior to time-stepping the equation.
56Generally, the time stepping is performed once at each time step in the \mdl{tranxt} and \mdl{dynnxt} modules,
57except when using implicit vertical diffusion or calculating sea surface height in which
58case time-splitting options are used.
59
60% -------------------------------------------------------------------------------------------------------------
61%        Non-Diffusive Part---Leapfrog Scheme
62% -------------------------------------------------------------------------------------------------------------
63\section{Non-diffusive part --- Leapfrog scheme}
64\label{sec:STP_leap_frog}
65
66The time stepping used for processes other than diffusion is the well-known leapfrog scheme
67\citep{Mesinger_Arakawa_Bk76}.
68This scheme is widely used for advection processes in low-viscosity fluids.
69It is a time centred scheme, \ie the RHS in \autoref{eq:STP} is evaluated at time step $t$, the now time step.
70It may be used for momentum and tracer advection, pressure gradient, and Coriolis terms,
71but not for diffusion terms.
72It is an efficient method that achieves second-order accuracy with
73just one right hand side evaluation per time step.
74Moreover, it does not artificially damp linear oscillatory motion nor does it produce instability by
75amplifying the oscillations.
76These advantages are somewhat diminished by the large phase-speed error of the leapfrog scheme,
77and the unsuitability of leapfrog differencing for the representation of diffusion and Rayleigh damping processes.
78However, the scheme allows the coexistence of a numerical and a physical mode due to
79its leading third order dispersive error.
80In other words a divergence of odd and even time steps may occur.
81To prevent it, the leapfrog scheme is often used in association with a Robert-Asselin time filter
82(hereafter the LF-RA scheme).
83This filter, first designed by \citet{Robert_JMSJ66} and more comprehensively studied by \citet{Asselin_MWR72},
84is a kind of laplacian diffusion in time that mixes odd and even time steps:
85\begin{equation}
86  \label{eq:STP_asselin}
87  x_F^t  = x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]
88\end{equation} 
89where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin coefficient.
90$\gamma$ is initialized as \np{rn\_atfp} (namelist parameter).
91Its default value is \np{rn\_atfp}\forcode{ = 10.e-3} (see \autoref{sec:STP_mLF}),
92causing only a weak dissipation of high frequency motions (\citep{Farge1987}).
93The addition of a time filter degrades the accuracy of the calculation from second to first order.
94However, the second order truncation error is proportional to $\gamma$, which is small compared to 1.
95Therefore, the LF-RA is a quasi second order accurate scheme.
96The LF-RA scheme is preferred to other time differencing schemes such as predictor corrector or trapezoidal schemes,
97because the user has an explicit and simple control of the magnitude of the time diffusion of the scheme.
98When used with the 2nd order space centred discretisation of the advection terms in
99the momentum and tracer equations, LF-RA avoids implicit numerical diffusion:
100diffusion is set explicitly by the user through the Robert-Asselin
101filter parameter and the viscosity and diffusion coefficients.
102
103% -------------------------------------------------------------------------------------------------------------
104%        Diffusive Part---Forward or Backward Scheme
105% -------------------------------------------------------------------------------------------------------------
106\section{Diffusive part --- Forward or backward scheme}
107\label{sec:STP_forward_imp}
108
109The leapfrog differencing scheme is unsuitable for the representation of diffusion and damping processes.
110For a tendancy $D_x$, representing a diffusion term or a restoring term to a tracer climatology
111(when present, see \autoref{sec:TRA_dmp}), a forward time differencing scheme is used :
112\[
113  % \label{eq:STP_euler}
114   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \ {D_x}^{t-\rdt}
115\] 
116
117This is diffusive in time and conditionally stable.
118The conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}:
119\begin{equation}
120  \label{eq:STP_euler_stability}
121  A^h < \left\{
122    \begin{aligned}
123      &\frac{e^2}{  8 \; \rdt }  &&\quad \text{laplacian diffusion}  \\
124      &\frac{e^4}{64 \; \rdt }   &&\quad \text{bilaplacian diffusion}
125    \end{aligned}
126  \right.
127\end{equation}
128where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is the mixing coefficient.
129The linear constraint \autoref{eq:STP_euler_stability} is a necessary condition, but not sufficient.
130If it is not satisfied, even mildly, then the model soon becomes wildly unstable.
131The instability can be removed by either reducing the length of the time steps or reducing the mixing coefficient.
132
133For the vertical diffusion terms, a forward time differencing scheme can be used,
134but usually the numerical stability condition imposes a strong constraint on the time step.
135Two solutions are available in \NEMO to overcome the stability constraint:
136$(a)$ a forward time differencing scheme using a time splitting technique (\np{ln\_zdfexp}\forcode{ = .true.}) or
137$(b)$ a backward (or implicit) time differencing scheme                   (\np{ln\_zdfexp}\forcode{ = .false.}).
138In $(a)$, the master time step $\Delta $t is cut into $N$ fractional time steps so that
139the stability criterion is reduced by a factor of $N$.
140The computation is performed as follows:
141\[
142  % \label{eq:STP_ts}
143  \begin{split}
144    & x_\ast ^{t-\rdt} = x^{t-\rdt}   \\
145    & x_\ast ^{t-\rdt+L\frac{2\rdt}{N}}=x_\ast ^{t-\rdt+\left( {L-1}
146      \right)\frac{2\rdt}{N}}+\frac{2\rdt}{N}\;\text{DF}^{t-\rdt+\left( {L-1} \right)\frac{2\rdt}{N}}
147    \quad \text{for $L=1$ to $N$}      \\
148    & x^{t+\rdt} = x_\ast^{t+\rdt}
149  \end{split}
150\]
151with DF a vertical diffusion term.
152The number of fractional time steps, $N$, is given by setting \np{nn\_zdfexp}, (namelist parameter).
153The scheme $(b)$ is unconditionally stable but diffusive. It can be written as follows:
154\begin{equation}
155  \label{eq:STP_imp}
156  x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t+\rdt}
157\end{equation} 
158
159%%gm
160%%gm   UPDATE the next paragraphs with time varying thickness ...
161%%gm
162
163This scheme is rather time consuming since it requires a matrix inversion,
164but it becomes attractive since a value of 3 or more is needed for N in the forward time differencing scheme.
165For example, the finite difference approximation of the temperature equation is:
166\[
167  % \label{eq:STP_imp_zdf}
168  \frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\rdt}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta
169  _k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
170  \right]
171\]
172where RHS is the right hand side of the equation except for the vertical diffusion term.
173We rewrite \autoref{eq:STP_imp} as:
174\begin{equation}
175  \label{eq:STP_imp_mat}
176  -c(k+1)\;T^{t+1}(k+1) + d(k)\;T^{t+1}(k) - \;c(k)\;T^{t+1}(k-1) \equiv b(k)
177\end{equation}
178where
179\begin{align*}
180  c(k) &= A_w^{vT} (k) \, / \, e_{3w} (k)     \\
181  d(k) &= e_{3t} (k)       \, / \, (2\rdt) + c_k + c_{k+1}    \\
182  b(k) &= e_{3t} (k) \; \left( T^{t-1}(k) \, / \, (2\rdt) + \text{RHS} \right)   
183\end{align*}
184
185\autoref{eq:STP_imp_mat} is a linear system of equations with an associated matrix which is tridiagonal.
186Moreover,
187$c(k)$ and $d(k)$ are positive and the diagonal term is greater than the sum of the two extra-diagonal terms,
188therefore a special adaptation of the Gauss elimination procedure is used to find the solution
189(see for example \citet{Richtmyer1967}).
190
191
192
193% -------------------------------------------------------------------------------------------------------------
194%        Surface Pressure gradient
195% -------------------------------------------------------------------------------------------------------------
196\section{Surface pressure gradient}
197\label{sec:STP_spg_ts}
198
199===>>>>  TO BE written....  :-)
200
201%\gmcomment{
202%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
203\begin{figure}[!t]
204  \begin{center}
205    \includegraphics[width=0.7\textwidth]{Fig_TimeStepping_flowchart}
206    \caption{
207      \protect\label{fig:TimeStep_flowchart}
208      Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}.
209      The use of a semi-implicit computation of the hydrostatic pressure gradient requires the tracer equation to
210      be stepped forward prior to the momentum equation.
211      The need for knowledge of the vertical scale factor (here denoted as $h$) requires the sea surface height and
212      the continuity equation to be stepped forward prior to the computation of the tracer equation.
213      Note that the method for the evaluation of the surface pressure gradient $\nabla p_s$ is not presented here
214      (see \autoref{sec:DYN_spg}).
215    }
216  \end{center}
217\end{figure}
218%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
219%}
220
221% -------------------------------------------------------------------------------------------------------------
222%        The Modified Leapfrog -- Asselin Filter scheme
223% -------------------------------------------------------------------------------------------------------------
224\section{Modified Leapfrog -- Asselin filter scheme}
225\label{sec:STP_mLF}
226
227Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the LF-RA scheme in order to ensure tracer conservation and to allow the use of a much smaller value of the Asselin filter parameter.
228The modifications affect both the forcing and filtering treatments in the LF-RA scheme.
229
230In a classical LF-RA environment, the forcing term is centred in time,
231\ie it is time-stepped over a $2\rdt$ period:
232$x^t  = x^t + 2\rdt Q^t $ where $Q$ is the forcing applied to $x$,
233and the time filter is given by \autoref{eq:STP_asselin} so that $Q$ is redistributed over several time step.
234In the modified LF-RA environment, these two formulations have been replaced by:
235\begin{align}
236  x^{t+\rdt}  &= x^{t-\rdt} + \rdt \left( Q^{t-\rdt/2} + Q^{t+\rdt/2} \right)                   \label{eq:STP_forcing} \\
237  %
238  x_F^&= x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]
239           - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]                          \label{eq:STP_RA}
240\end{align}
241The change in the forcing formulation given by \autoref{eq:STP_forcing} (see \autoref{fig:MLF_forcing})
242has a significant effect:
243the forcing term no longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}.
244% forcing seen by the model....
245This property improves the LF-RA scheme in two respects.
246First, the LF-RA can now ensure the local and global conservation of tracers.
247Indeed, time filtering is no longer required on the forcing part.
248The influence of the Asselin filter on the forcing is be removed by adding a new term in the filter
249(last term in \autoref{eq:STP_RA} compared to \autoref{eq:STP_asselin}).
250Since the filtering of the forcing was the source of non-conservation in the classical LF-RA scheme,
251the modified formulation becomes conservative \citep{Leclair_Madec_OM09}.
252Second, the LF-RA becomes a truly quasi-second order scheme.
253Indeed, \autoref{eq:STP_forcing} used in combination with a careful treatment of static instability
254(\autoref{subsec:ZDF_evd}) and of the TKE physics (\autoref{subsec:ZDF_tke_ene}),
255the two other main sources of time step divergence,
256allows a reduction by two orders of magnitude of the Asselin filter parameter.
257
258Note that the forcing is now provided at the middle of a time step:
259$Q^{t+\rdt/2}$ is the forcing applied over the $[t,t+\rdt]$ time interval.
260This and the change in the time filter, \autoref{eq:STP_RA},
261allows an exact evaluation of the contribution due to the forcing term between any two time steps,
262even if separated by only $\rdt$ since the time filter is no longer applied to the forcing term.
263
264%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
265\begin{figure}[!t]
266  \begin{center}
267    \includegraphics[width=0.90\textwidth]{Fig_MLF_forcing}
268    \caption{
269      \protect\label{fig:MLF_forcing}
270      Illustration of forcing integration methods.
271      (top) ''Traditional'' formulation:
272      the forcing is defined at the same time as the variable to which it is applied
273      (integer value of the time step index) and it is applied over a $2\rdt$ period.
274      (bottom)  modified formulation:
275      the forcing is defined in the middle of the time (integer and a half value of the time step index) and
276      the mean of two successive forcing values ($n-1/2$, $n+1/2$) is applied over a $2\rdt$ period.
277    }
278  \end{center}
279\end{figure}
280%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
281
282% -------------------------------------------------------------------------------------------------------------
283%        Start/Restart strategy
284% -------------------------------------------------------------------------------------------------------------
285\section{Start/Restart strategy}
286\label{sec:STP_rst}
287
288%--------------------------------------------namrun-------------------------------------------
289\nlst{namrun}
290%--------------------------------------------------------------------------------------------------------------
291
292The first time step of this three level scheme when starting from initial conditions is a forward step
293(Euler time integration):
294\[
295  % \label{eq:DOM_euler}
296  x^1 = x^0 + \rdt \ \text{RHS}^0
297\]
298This is done simply by keeping the leapfrog environment (\ie the \autoref{eq:STP} three level time stepping) but
299setting all $x^0$ (\textit{before}) and $x^{1}$ (\textit{now}) fields equal at the first time step and
300using half the value of $\rdt$.
301
302It is also possible to restart from a previous computation, by using a restart file.
303The restart strategy is designed to ensure perfect restartability of the code:
304the user should obtain the same results to machine precision either by
305running the model for $2N$ time steps in one go,
306or by performing two consecutive experiments of $N$ steps with a restart.
307This requires saving two time levels and many auxiliary data in the restart files in machine precision.
308
309Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure gradient
310(see \autoref{subsec:DYN_hpg_imp}), an extra three-dimensional field has to
311be added to the restart file to ensure an exact restartability.
312This is done optionally via the \np{nn\_dynhpg\_rst} namelist parameter,
313so that the size of the restart file can be reduced when restartability is not a key issue
314(operational oceanography or in ensemble simulations for seasonal forecasting).
315
316Note the size of the time step used, $\rdt$, is also saved in the restart file.
317When restarting, if the the time step has been changed, a restart using an Euler time stepping scheme is imposed.
318Options are defined through the \ngn{namrun} namelist variables.
319%%%
320\gmcomment{
321add here how to force the restart to contain only one time step for operational purposes
322
323add also the idea of writing several restart for seasonal forecast : how is it done ?
324
325verify that all namelist pararmeters are truly described
326
327a word on the check of restart  .....
328}
329%%%
330
331\gmcomment{       % add a subsection here 
332
333%-------------------------------------------------------------------------------------------------------------
334%        Time Domain
335% -------------------------------------------------------------------------------------------------------------
336\subsection{Time domain}
337\label{subsec:STP_time}
338%--------------------------------------------namrun-------------------------------------------
339
340\nlst{namdom}         
341%--------------------------------------------------------------------------------------------------------------
342
343Options are defined through the  \ngn{namdom} namelist variables.
344 \colorbox{yellow}{add here a few word on nit000 and nitend}
345
346 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj}
347
348add a description of daymod, and the model calandar (leap-year and co)
349
350}        %% end add
351
352
353
354%%
355\gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme   
356
357Implicit time stepping in case of variable volume thickness.
358
359Tracer case (NB for momentum in vector invariant form take care!)
360
361\begin{flalign*}
362  &\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\rdt}
363  \equiv \text{RHS}+ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
364  \right]      \\
365  &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
366  \equiv {2\rdt} \ \text{RHS}+ {2\rdt} \ \delta_k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta_{k+1/2} \left[ {T^{t+1}} \right]}
367  \right]      \\
368  &\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
369  \equiv 2\rdt \ \text{RHS}
370  + 2\rdt \ \left\{ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} [ T_{k+1}^{t+1} - T_k      ^{t+1} ]
371    - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\
372  &\\
373  &\left( e_{3t}\,T \right)_k^{t+1}
374  -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1}
375  + {2\rdt} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
376    +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1}
377  -  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\
378  &\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
379  %
380\end{flalign*}
381
382\begin{flalign*}
383  \allowdisplaybreaks
384  \intertext{ Tracer case }
385  %
386  &  \qquad \qquad  \quad   -  {2\rdt}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
387  \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\
388  &+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}
389  +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
390  & \qquad \qquad  \qquad \qquad \qquad \quad \ \ {2\rdt} \                          \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                          \quad \ \ T_{k-1}^{t+1}
391  \ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
392  %
393\end{flalign*}
394\begin{flalign*}
395  \allowdisplaybreaks
396  \intertext{ Tracer content case }
397  %
398  & -  {2\rdt} \              & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_{k+1}^{t+1}}  && \  \left( e_{3t}\,T \right)_{k+1}^{t+1}   &\\
399  & + {2\rdt} \ \left[ 1  \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}}
400  + & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_k^{t+1}}  \left\right& \left( e_{3t}\,T \right)_{k   }^{t+1}  &\\
401  & -  {2\rdt} \               & \frac{(A_w^{vt})_{k -1/2}} {(e_{3w})_{k -1/2}^{t+1}\;(e_{3t})_{k-1}^{t+1}}     &\  \left( e_{3t}\,T \right)_{k-1}^{t+1}
402  \equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  &
403\end{flalign*}
404
405%%
406}
407%%
408\biblio
409
410\pindex
411
412\end{document}
Note: See TracBrowser for help on using the repository browser.