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_step.tex in branches/2017/dev_merge_2017/DOC/texfiles/chapters – NEMO

source: branches/2017/dev_merge_2017/DOC/texfiles/chapters/chap_step.tex @ 9373

Last change on this file since 9373 was 9373, checked in by nicolasmartin, 6 years ago

Global reorganisation of DOC directory: files/folders renaming

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