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_STP.tex in branches/NERC/dev_r5549_BDY_ZEROGRAD/DOC/TexFiles/Chapters – NEMO

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/DOC/TexFiles/Chapters/Chap_STP.tex @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

merge with trunk@6232 for consistency with SSB code

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