source: branches/DEV_r1826_DOC/DOC/TexFiles/Chapters/Chap_STP.tex @ 2164

% ================================================================
% Chapter 2 Ñ Time Domain (step.F90)
% ================================================================
\chapter{Time Domain (STP) }
\label{STP}
\minitoc

% Missing things:
%  - daymod: definition of the time domain (nit000, nitend andd the calendar)


\gmcomment{STEVEN :maybe a picture of the directory structure in the introduction 
which could be referred to here, would help  ==> to be added}
%%%%


\newpage
$\ $\newline    % force a new ligne


Having defined the continuous equations in Chap.~\ref{PE}, we need now to choose 
a time discretization. In the present chapter, we provide a general description of the \NEMO 
time stepping strategy and its consequences for the order in which the equations are
solved.

$\ $\newline    % force a new ligne
% ================================================================
% Time Discretisation
% ================================================================
\section{Time stepping environment}
\label{STP_environment}

The time stepping used in \NEMO is a three level scheme that can be 
represented as follows:
\begin{equation} \label{Eq_STP}
   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t-\rdt,\,t,\,t+\rdt}
\end{equation} 
where $x$ stands for $u$, $v$, $T$ or $S$; RHS is the Right-Hand-Side of the 
corresponding time evolution equation; $\rdt$ is the time step; and the 
superscripts indicate the time at which a quantity is evaluated. Each term of the 
RHS is evaluated at a specific time step(s) depending on the physics with which 
it is associated. 

The choice of the time step used for this evaluation is discussed below as 
well as the implications in terms of starting or restarting a model 
simulation. Note that the time stepping is generally performed in a one step 
operation. With such a complex and nonlinear system of equations it would be 
dangerous to let a prognostic variable evolve in time for each term separately.

The three level scheme requires three arrays for each prognostic variables. 
For each variable $x$ there is $x_b$ (before) and $x_n$ (now). The third array, 
although referred to as $x_a$ (after) in the code, is usually not the variable at 
the after time step; but rather it is used to store the time derivative (RHS in 
\eqref{Eq_STP}) prior to time-stepping the equation. Generally, the time 
stepping is performed once at each time step in \mdl{tranxt} and \mdl{dynnxt} 
modules, except for implicit vertical diffusion or sea surface height when 
time-splitting options are used.

% -------------------------------------------------------------------------------------------------------------
%        Non-Diffusive Part---Leapfrog Scheme
% -------------------------------------------------------------------------------------------------------------
\section{Non-Diffusive Part --- Leapfrog Scheme}
\label{STP_leap_frog}

The time stepping used for non-diffusive processes is the well-known leapfrog
scheme \citep{Mesinger_Arakawa_Bk76}. It is a time centred scheme, $i.e.$ 
the RHS in \eqref{Eq_STP} is evaluated at time step $t$, the now time step. 
It is only used for non-diffusive terms, that is momentum and tracer advection, 
pressure gradient, and Coriolis terms. This scheme is widely used for advective 
processes in low-viscosity fluids. It is an efficient method that achieves 
second-order accuracy with just one right hand side evaluation per time step. 
Moreover, it does not artificially damp linear oscillatory motion nor does it produce 
instability by amplifying the oscillations. These advantages are somewhat diminished 
by the large phase-speed error of the leapfrog scheme, and the unsuitability 
of leapfrog differencing for the representation of diffusive and Rayleigh 
damping processes. However, the scheme allows the coexistence of a numerical 
and a physical mode due to its leading third order dispersive error. In other words a 
divergence of odd and even time steps may occur. To prevent it, the leapfrog scheme 
is often used in association with a Robert-Asselin time filter (thereafter the LF-RA scheme). 
This filter, first designed by \citet{Robert_JMSJ66} and more comprehensively studied 
by \citet{Asselin_MWR72}, is a kind of laplacian diffusion in time that mixes odd and 
even time steps:
\begin{equation} \label{Eq_STP_asselin}
x_F^t  = x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right]
\end{equation} 
where the subscript $F$ denotes filtered values and $\gamma$ is the Asselin 
coefficient. $\gamma$ is initialized as \np{rn\_atfp} (namelist parameter). 
Its default value is \textit{now} (see \S~\ref{STP_mLF}) \np{rn\_atfp}=$10^{-3}$, 
causing only a weak dissipation of high frequency motions (\citep{Farge1987}). 
The addition of a time filter does, nevertheless, degrade the accuracy of the 
calculation from second to first order. However, the second order truncation 
error is proportional to $\gamma$, which is small compared to 1. Therefore, 
the LF-RA is a quasi second order accurate scheme. The LF-RA scheme 
has been preferred to other time differencing schemes such as 
predictor corrector or trapezoidal schemes, because the user has an explicit 
and simple control of the magnitude of the time diffusion of the scheme. 
In association with the 2nd order centred space discretisation of the 
advective terms in the momentum and tracer equations, it avoids implicit 
numerical diffusion in both the time and space discretisations of the 
advective term: they are both set explicitly by the user through the Robert-Asselin 
filter parameter and the viscous and diffusive coefficients.

% -------------------------------------------------------------------------------------------------------------
%        Diffusive Part---Forward or Backward Scheme
% -------------------------------------------------------------------------------------------------------------
\section{Diffusive Part --- Forward or Backward Scheme}
\label{STP_forward_imp}

The leapfrog differencing scheme is unsuitable for the representation of 
diffusive and damping processes. For a tendancy $D_x$, representing a 
diffusive term or a restoring term to a tracer climatology 
(when present, see \S~\ref{TRA_dmp}), a forward time differencing scheme
 is used :
\begin{equation} \label{Eq_STP_euler}
   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \ {D_x}^{t-\rdt}
\end{equation} 

This is diffusive in time and conditionally stable. For example, the 
conditions for stability of second and fourth order horizontal diffusion schemes are \citep{Griffies_Bk04}:
\begin{equation} \label{Eq_STP_euler_stability}
A^h < \left\{
\begin{aligned}
                    &\frac{e^2}{  8 \; \rdt }   &&\quad \text{laplacian diffusion}  \\
                    &\frac{e^4}{64 \; \rdt }    &&\quad \text{bilaplacian diffusion} 
            \end{aligned}
\right.
\end{equation}
where $e$ is the smallest grid size in the two horizontal directions and $A^h$ is 
the mixing coefficient. The linear constraint \eqref{Eq_STP_euler_stability} 
is a necessary condition, but not sufficient. If it is not satisfied, even mildly, 
then the model soon becomes wildly unstable. The instability can be removed 
by either reducing the length of the time steps or reducing the mixing coefficient.

For the vertical diffusion terms, a forward time differencing scheme can be 
used, but usually the numerical stability condition implies a strong 
constraint on the time step. Two solutions are available in \NEMO to overcome 
the stability constraint: $(a)$ a forward time differencing scheme using a 
time splitting technique (\np{ln\_zdfexp}=.true.) or $(b)$ a backward (or implicit) 
time differencing scheme by \np{ln\_zdfexp}=.false.). In $(a)$, the master 
time step $\Delta $t is cut into $N$ fractional time steps so that the 
stability criterion is reduced by a factor of $N$. The computation is done as 
follows:
\begin{equation} \label{Eq_STP_ts}
\begin{split}
& x_\ast ^{t-\rdt} = x^{t-\rdt}   \\
& x_\ast ^{t-\rdt+L\frac{2\rdt}{N}}=x_\ast ^{t-\rdt+\left( {L-1} 
\right)\frac{2\rdt}{N}}+\frac{2\rdt}{N}\;\text{DF}^{t-\rdt+\left( {L-1} \right)\frac{2\rdt}{N}}
        \quad \text{for $L=1$ to $N$}      \\
& x^{t+\rdt} = x_\ast^{t+\rdt}
\end{split}
\end{equation}
with DF a vertical diffusion term. The number of fractional time steps, $N$, is given 
by setting \np{n\_zdfexp}, (namelist parameter). The scheme $(b)$ is unconditionally 
stable but diffusive. It can be written as follows:
\begin{equation} \label{Eq_STP_imp}
   x^{t+\rdt} = x^{t-\rdt} + 2 \, \rdt \  \text{RHS}_x^{t+\rdt}
\end{equation} 

This scheme is rather time consuming since it requires a matrix inversion, 
but it becomes attractive since a splitting factor of 3 or more is needed 
for the forward time differencing scheme. For example, the finite difference 
approximation of the temperature equation is:
\begin{equation} \label{Eq_STP_imp_zdf}
\frac{T(k)^{t+1}-T(k)^{t-1}}{2\;\rdt}\equiv \text{RHS}+\frac{1}{e_{3t} }\delta 
_k \left[ {\frac{A_w^{vT} }{e_{3w} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} 
\right]
\end{equation}
where RHS is the right hand side of the equation except for the vertical diffusion term. 
We rewrite \eqref{Eq_STP_imp} as:
\begin{equation} \label{Eq_STP_imp_mat}
-c(k+1)\;u^{t+1}(k+1) + d(k)\;u^{t+1}(k) - \;c(k)\;u^{t+1}(k-1) \equiv b(k)
\end{equation}
where 
\begin{align*} 
 c(k) &= A_w^{vm} (k) \, / \, e_{3uw} (k)     \\
 d(k) &= e_{3u} (k)       \, / \, (2\rdt) + c_k + c_{k+1}    \\
 b(k) &= e_{3u} (k) \; \left( u^{t-1}(k) \, / \, (2\rdt) + \text{RHS} \right)    
\end{align*}

\eqref{Eq_STP_imp_mat} is a linear system of equations which associated 
matrix is tridiagonal. Moreover, $c(k)$ and $d(k)$ are positive and the diagonal 
term is greater than the sum of the two extra-diagonal terms, therefore a special 
adaptation of the Gauss elimination procedure is used to find the solution 
(see for example \citet{Richtmyer1967}).



% -------------------------------------------------------------------------------------------------------------
%        Hydrostatic Pressure gradient
% -------------------------------------------------------------------------------------------------------------
\section{Hydrostatic Pressure Gradient --- semi-implicit scheme}
\label{STP_hpg_imp}

%\gmcomment{ 
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \label{Fig_TimeStep_flowchart}  \begin{center}
\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_TimeStepping_flowchart.pdf}
\caption{Sketch of the leapfrog time stepping sequence in \NEMO from \citet{Leclair_Madec_OM09}. 
The use of semi-implicit computation of the hydrostatic pressure gradient requires
the tracer equation to be step forward prior to the momentum equation. 
The need of the knowledge of the vertical scale factor (here denoted as $h$)
requires the sea surface height and the continuity equation to be step forward
prior to the computation of the tracer equation.
Note that how the surface pressure gradient $\nabla p_s$ is evaluated is not represented here 
(see \S~\ref{DYN_spg}). }
\end{center}   \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%}

The range of stability of the Leap-Frog scheme can be extended by a factor of two
by introducing a semi-implicit computation of the hydrostatic pressure gradient term
\citep{Brown_Campana_MWR78}. Instead of evaluating the pressure at $t$, a linear 
combination of $t-\rdt$, $t$ and $t+\rdt$ is used (see \S~\ref{DYN_hpg_imp}).  
This technique, controlled by  \np{nn\_dynhpg\_rst} namelist parameter, does not 
require a significant additional computational coast when tracers and thus density 
is time stepped before the dynamics, a choice made in \NEMO 
(Fig.\ref{Fig_TimeStep_flowchart}).


This technique, used in several GCMs (\NEMO, POP or MOM for instance), 
makes the Leap-Frog scheme as efficient 
\footnote{The efficiency is defined as the maximum allowed Courant number of the time 
stepping scheme divided by the number of computations of the right-hand side per time step.} 
as the Forward-Backward scheme used in MOM \citep{Griffies_al_OS05} and more 
efficient than the LF-AM3 scheme (leapfrog time stepping combined with a third order
Adams-Moulthon interpolation for the predictor phase) used in ROMS 
\citep{Shchepetkin_McWilliams_OM05}. 

In fact, this technique is efficient when the physical phenomenum that 
controls the time-step is internal gravity waves (IGWs). Indeed, it is 
equivalent to applying a time filter to the pressure gradient to eliminate high 
frequency IGWs. Obviously, the doubling of the time-step is achievable only 
if no other factors control the time-step, such as the stability limits associated 
with advection, diffusion or Coriolis terms. For example, it is useless in low resolution
global ocean configurations as inertial oscillations in vicinity of the North Pole 
are the limiting factor of the time step. It is also often useless in very high 
resolution configurations where the strong currents and small grid cells exert 
the strongest constraint on the time step.


% -------------------------------------------------------------------------------------------------------------
%        The Modified Leapfrog -- Asselin Filter scheme
% -------------------------------------------------------------------------------------------------------------
\section{The Modified Leapfrog -- Asselin Filter scheme}
\label{STP_mLF}

Significant changes have been introduced by \cite{Leclair_Madec_OM09} in the 
LF-RA scheme in order to ensure the tracer conservation and to allow the use of 
a much smaller value of the Asselin filter parameter. The modifications concern 
both the forcing and filtering treatment in the LF-RA scheme.

In a classical LF-RA environment, the forcing term is centred in time, $i.e.$ 
it is time-stepped over a $2\rdt$ period:  $x^t  = x^t + 2\rdt Q^t $ where $Q$ 
is the forcing applied on $x$, and the filter is given by \eqref{Eq_STP_asselin}. 
In the modified LF-RA environment, these two formulations have been replaced by:
\begin{align} 
x^{t+\rdt}  &= x^{t-\rdt} + \rdt \left( Q^{t-\rdt/2} + Q^{t+\rdt/2} \right)                   \label{Eq_STP_forcing} \\
%
x_F^t  &= x^t + \gamma \, \left[ x_F^{t-\rdt} - 2 x^t + x^{t+\rdt} \right] 
           - \gamma\,\rdt \, \left[ Q^{t+\rdt/2} -  Q^{t-\rdt/2} \right]                          \label{Eq_STP_RA}
\end{align}

The change in the forcing formulation given by \eqref{Eq_STP_forcing} 
(see Fig.\ref{Fig_MLF_forcing}) has a paramount effect: the forcing term no 
longer excites the divergence of odd and even time steps \citep{Leclair_Madec_OM09}. 
This property allows to improve the LF-RA scheme in two aspects. 
First, the LF-RA becomes a truly quasi-second order scheme. Indeed, 
\eqref{Eq_STP_forcing} used in combination with a careful treatment of the static 
instabilities (\S\ref{ZDF_evd}) and of the TKE physics (\S\ref{ZDF_tke_ene}),
the two main other sources of time step divergence, allows a reduction by two order 
of magnitude of the Asselin filter parameter. 
Second, the LF-RA can now ensure the local and global conservation of tracers.
Indeed, the time filtering is no more required on the forcing part. The influence of 
the forcing in the Asselin filter can be removed by adding a new term in the filter
(last term in \eqref{Eq_STP_RA} compared to \eqref{Eq_STP_asselin}). Since 
the filtering of the forcing was the source of non-conservation in the LF-RA
scheme, it becomes conservative  \citep{Leclair_Madec_OM09}.

Note that the forcing is now provided at the middle of a time step: $Q^{t+\rdt/2}$ 
is the forcing applied over the $[t,t+\rdt]$ time interval. This and the change 
in the time filter, \eqref{Eq_STP_RA}, allows an exact evaluation of the 
contribution due to the forcing term between any two time steps, 
even if separated by only $\rdt$ as the time filter is no more applied on the
forcing term.

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
\begin{figure}[!t] \label{Fig_MLF_forcing}  \begin{center}
\includegraphics[width=0.90\textwidth]{./TexFiles/Figures/Fig_MLF_forcing.pdf}
\caption{Illustration of forcing integration methods. ''Traditional'' formulation (top) 
where a centred forcing is applied over a $2\rdt$ period and modified formulation 
(bottom) where a mean forcing over the two successive time step is applied over a $2\rdt$ period.}
\end{center}   \end{figure}
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>


% -------------------------------------------------------------------------------------------------------------
%        Start/Restart strategy
% -------------------------------------------------------------------------------------------------------------
\section{Start/Restart strategy}
\label{STP_rst}
%--------------------------------------------namrun-------------------------------------------
\namdisplay{namrun}         
%--------------------------------------------------------------------------------------------------------------

The first time step of this three level scheme when starting from initial conditions 
is a forward step (Euler time integration): 
\begin{equation} \label{Eq_DOM_euler}
   x^1 = x^0 + \rdt \ \text{RHS}^0
\end{equation}
This is done simply by keeping the leapfrog environment but setting at the first time step
only both a half $\rdt$ and the equality of all \textit{before} and \textit{now} fields.

It is also possible to restart from a previous computation, by using a 
restart file. The restart strategy is designed to ensure perfect 
restartability of the code: the user should obtain the same results to 
machine precision either by running the model for $2N$ time steps in one go, 
or by performing two consecutive experiments of $N$ steps with a restart. 
This requires saving two time levels and many auxiliary data in the restart 
files in machine precision. 

Note that when a semi-implicit scheme is used to evaluate the hydrostatic pressure 
gradient (see \S\ref{DYN_hpg_imp}), an extra three-dimensional field has to be 
added in the restart file to ensure an exact restartability. This is done only optionally 
via the namelist parameter \np{nn\_dynhpg\_rst}, so that a reduction of the size of 
restart file can be obtained when the restartability is not a key issue (operational 
oceanography or ensemble simulation for seasonal forcast).

Note the size of the time step used, $\rdt$, is also saved in the restart file. 
When restarting, if the the time step has been changed, a restart using an Euler time 
stepping is imposed. 
%%%
\gmcomment{
add here how to force the restart to contain only one time step for operational purposes

add also the idea of writing several restart for seasonal forecast : how is it done ?

verify that all namelist pararmeters are truly described 

a word on the check of restart  .....
}
%%%

\gmcomment{       % add a subsection here  

%-------------------------------------------------------------------------------------------------------------
%        Time Domain
% -------------------------------------------------------------------------------------------------------------
\subsection{Time domain}
\label{STP_time}
%--------------------------------------------namrun-------------------------------------------
\namdisplay{namdom}         
%--------------------------------------------------------------------------------------------------------------


 \colorbox{yellow}{add here a few word on nit000 and nitend}

 \colorbox{yellow}{Write documentation on the calendar and the key variable adatrj}

add a description of daymod, and the model calandar (leap-year and co)

}        %% end add



%%
\gmcomment{       % add implicit in vvl case  and Crant-Nicholson scheme   

Implicit time stepping in case of variable volume thickness.

Tracer case (NB for momentum in vector invariant form take care!)

\begin{flalign*}
&\frac{\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}}{2\rdt}
\equiv \text{RHS}+ \delta _k \left[ {\frac{A_w^{vt} }{e_{3w}^{t+1} }\delta _{k+1/2} \left[ {T^{t+1}} \right]} 
\right]      \\
&\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
\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]} 
\right]      \\
&\left( e_{3t}\,T \right)_k^{t+1}-\left( e_{3t}\,T \right)_k^{t-1}
\equiv 2\rdt \ \text{RHS}
+ 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} ]
                          - \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} [ T_k       ^{t+1} - T_{k-1}^{t+1} ]  \right\}     \\
&\\
&\left( e_{3t}\,T \right)_k^{t+1}
-  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}                  T_{k+1}^{t+1} 
+ {2\rdt} \ \left\{  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2} 
                            +  \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}     \right\}   T_{k    }^{t+1}
-  {2\rdt} \           \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2}                  T_{k-1}^{t+1}      \\
&\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}    \\
%
\end{flalign*}

\begin{flalign*}
\allowdisplaybreaks
\intertext{ Tracer case }
%
&  \qquad \qquad  \quad   -  {2\rdt}                  \ \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}   
                                                                                                      \qquad \qquad \qquad  \qquad  T_{k+1}^{t+1}   \\
&+ {2\rdt} \ \biggl\{  (e_{3t})_{k   }^{t+1}  \bigg. +    \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k+1/2}  
                                                                               +   \left[ \frac{A_w^{vt}}{e_{3w}^{t+1}} \right]_{k-1/2} \bigg. \biggr\}  \ \ \ T_{k   }^{t+1}  &&\\
& \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}     
\ \equiv \ \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  \\
%
\end{flalign*}
\begin{flalign*}
\allowdisplaybreaks
\intertext{ Tracer content case }
%
& -  {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}   &\\
& + {2\rdt} \ \left[ 1  \right.+ & \frac{(A_w^{vt})_{k+1/2}} {(e_{3w})_{k+1/2}^{t+1}\;(e_{3t})_k^{t+1}} 
                                    + & \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}  &\\
& -  {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}   
\equiv \left( e_{3t}\,T \right)_k^{t-1} + {2\rdt} \ \text{RHS}  &
\end{flalign*}

%%
}
%%