source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_STO.tex @ 10368

\documentclass[../tex_main/NEMO_manual]{subfiles}
\begin{document}
% ================================================================
% Chapter stochastic parametrization of EOS (STO)
% ================================================================
\chapter{Stochastic Parametrization of EOS (STO)}
\label{chap: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{sec: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
(\autoref{eq:autoreg}), and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$.
The parameters in (\autoref{eq:ord2}) are computed so that this recursive application of
(\autoref{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{sec:STO_thech_details}

%---------------------------------------namsbc--------------------------------------------------

\nlst{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 \autoref{eq:sto_pert} and specific piece of code in
  the equation of state implementing \autoref{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 (\autoref{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 \forcode{.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 \forcode{.true.},
$i.e.$ when the state of the random number generator is read in the restart file.


\end{document}