# source:branches/2013/dev_LOCEAN_2013/DOC/TexFiles/Chapters/Chap_STP.tex@4147 Last change on this file since 4147 was 4147, checked in by cetlod, 7 years ago

merge in dev_LOCEAN_2013, the 1st development branch dev_r3853_CNRS9_Confsetting, from its starting point ( r3853 ) on the trunk: see ticket #1169

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