source: branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_STO.tex @ 6347

% ================================================================
% Chapter stochastic parametrization of EOS (STO)
% ================================================================
\chapter{Stochastic parametrization of EOS (STO)}
\label{STO}

Authors: P.-A. Bouttier

\minitoc

\newpage


The stochastic parametrization module aims to explicitly simulate uncertainties in the model. 
More particularly, \cite{Brankart_OM2013} has shown that, 
because of the nonlinearity of the seawater equation of state, unresolved scales represent 
a major source of uncertainties in the computation of the large scale horizontal density gradient 
(from T/S large scale fields), and that the impact of these uncertainties can be simulated 
by random processes representing unresolved T/S fluctuations.

The stochastic formulation of the equation of state can be written as:
\begin{equation}
 \label{eq:eos_sto}
  \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \}
\end{equation}
where $p_o(z)$ is the reference pressure depending on the depth and, 
$\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product 
of the respective local T/S gradients with random walks $\mathbf{\xi}$:
\begin{equation}
 \label{eq:sto_pert}
 \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S
\end{equation}
$\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with 
a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$. 
$\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical.


\section{Stochastic processes}
\label{STO_the_details}

The starting point of our implementation of stochastic parameterizations
in NEMO is to observe that many existing parameterizations are based
on autoregressive processes, which are used as a basic source of randomness
to transform a deterministic model into a probabilistic model.
A generic approach is thus to add one single new module in NEMO,
generating processes with appropriate statistics
to simulate each kind of uncertainty in the model
(see \cite{Brankart_al_GMD2015} for more details).

In practice, at every model grid point, independent Gaussian autoregressive
processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated
using the same basic equation:

\begin{equation}
\label{eq:autoreg}
\xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)}
\end{equation}

\noindent
where $k$ is the index of the model timestep; and
$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining
the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$)
and correlation timescale ($\tau^{(i)}$) of each process:

\begin{itemize}
\item for order~1 processes, $w^{(i)}$ is a Gaussian white noise,
with zero mean and standard deviation equal to~1, and the parameters
$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:

\begin{equation}
\label{eq:ord1}
\left\{
\begin{array}{l}
a^{(i)} = \varphi \\
b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 } 
 \qquad\qquad\mbox{with}\qquad\qquad
\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
\end{array}
\right.
\end{equation}

\item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process,
with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale
equal to~$\tau^{(i)}$; and the parameters $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:

\begin{equation}
\label{eq:ord2}
\left\{
\begin{array}{l}
a^{(i)} = \varphi \\
b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 } 
 \qquad\qquad\mbox{with}\qquad\qquad
\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
\end{array}
\right.
\end{equation}

\end{itemize}

\noindent
In this way, higher order processes can be easily generated recursively using 
the same piece of code implementing Eq.~(\ref{eq:autoreg}), 
and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$.
The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application
of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation
and correlation timescale, with the additional condition that
the $n-1$ first derivatives of the autocorrelation function
are equal to zero at~$t=0$, so that the resulting processes
become smoother and smoother as $n$ is increased.

Overall, this method provides quite a simple and generic way of generating 
a wide class of stochastic processes. 
However, this also means that new model parameters are needed to specify each of 
these stochastic processes. As in any parameterization of lacking physics, 
a very important issues then to tune these new parameters using either first principles, 
model simulations, or real-world observations.

\section{Implementation details}
\label{STO_thech_details}

%---------------------------------------namsbc--------------------------------------------------
\namdisplay{namsto}
%--------------------------------------------------------------------------------------------------------------

The computer code implementing stochastic parametrisations can be found in the STO directory.
It involves three modules : 
\begin{description}
\item[\mdl{stopar}] : define the Stochastic parameters and their time evolution.
\item[\mdl{storng}] : a random number generator based on (and includes) the 64-bit KISS 
                      (Keep It Simple Stupid) random number generator distributed by George Marsaglia 
                      (see \href{https://groups.google.com/forum/#!searchin/comp.lang.fortran/64-bit$20KISS$20RNGs}{here})
\item[\mdl{stopts}] : stochastic parametrisation associated with the non-linearity of the equation of seawater, 
 implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}.
\end{description}

The \mdl{stopar} module has 3 public routines to be called by the model (in our case, NEMO):

The first routine (\rou{sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}),
applied at each model grid point (in 2D or 3D), 
and called at each model time step ($k$) to update
every autoregressive process ($i=1,\ldots,m$).
This routine also includes a filtering operator, applied to $w^{(i)}$,
to introduce a spatial correlation between the stochastic processes.

The second routine (\rou{sto\_par\_init}) is an initialization routine mainly dedicated
to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$
for each autoregressive process, as a function of the statistical properties
required by the model user (mean, standard deviation, time correlation,
order of the process,\ldots). 

Parameters for the processes can be specified through the following \ngn{namsto} namelist parameters:
\begin{description}
   \item[\np{nn\_sto\_eos}]   : number of independent random walks 
   \item[\np{rn\_eos\_stdxy}] : random walk horz. standard deviation (in grid points)
   \item[\np{rn\_eos\_stdz}]  : random walk vert. standard deviation (in grid points)
   \item[\np{rn\_eos\_tcor}]  : random walk time correlation (in timesteps)
   \item[\np{nn\_eos\_ord}]   : order of autoregressive processes
   \item[\np{nn\_eos\_flt}]   : passes of Laplacian filter
   \item[\np{rn\_eos\_lim}]   : limitation factor (default = 3.0)
\end{description}
This routine also includes the initialization (seeding) of the random number generator.

The third routine (\rou{sto\_rst\_write}) writes a restart file (which suffix name is 
given by \np{cn\_storst\_out} namelist parameter) containing the current value of 
all autoregressive processes to allow restarting a simulation from where it has been interrupted.
This file also contains the current state of the random number generator.
When \np{ln\_rststo} is set to \textit{true}), the restart file (which suffix name is 
given by \np{cn\_storst\_in} namelist parameter) is read by the initialization routine 
(\rou{sto\_par\_init}). The simulation will continue exactly as if it was not interrupted
only  when \np{ln\_rstseed} is set to \textit{true}, $i.e.$ when the state of 
the random number generator is read in the restart file.