%
%$Id: meanflowIntro.tex,v 20.0 2013/12/14 00:13:18 fms Exp $
%

\section{The mean flow model \label{sec:meanflowIntro}}

\subsection{Introduction}

This module contains the definitions of the most important
mean flow variables used in geophysical models. In GOTM, these
are
\begin{itemize}
  \item the mean horizontal velocity components, $U$ and $V$
  \item the mean potential temperature, $\Theta$, (or the mean buoyancy, $B$)
  \item the mean salinity, $S$
\end{itemize}
Note that in general a variable $\phi$ describing a turbulent
field can be decomposed into a mean and a fluctuating part. In GOTM, 
we use the notation
\begin{equation}
  \label{decomposition}
    \phi = \mean{\phi} + \phi'
   \comma
\end{equation}
where $\mean{\rule{3mm}{0mm}}$ denotes the ensemble mean and the prime 
the fluctuating part. In addition, for brevity, we use the following conventions:
\begin{equation}
  \label{decompositionConventions}
   \begin{array}{rcl}
     U      &=& \mean{u}      \quad \text{for the x-velocity}            \\	
     V      &=& \mean{v}      \quad \text{for the y-velocity}            \\
     P      &=& \mean{p}      \quad \text{for the pressure}              \\
     \Theta &=& \mean{\theta} \quad \text{for the potential temperature} \\
     B      &=& \mean{b}      \quad \text{for the buoyancy}              \\
     S      &=& \mean{s}      \quad \text{for the salinity}             
   \end{array}
\end{equation}
Note that, if not explicitly mentioned, GOTM uses the units kg, m, s,
K. Further conventions are introduced in the turbulence chapter
\sect{sec:turbulenceIntro}. All operations on these meanflow variables
are executed and coordinated in the {\tt meanflow} module.

\subsubsection{Physics}\label{sec:meanflowIntroPhysics}
Due to the one-dimensional character of GOTM, the state-variables
listed above are assumed to be horizontally homogeneous, depending
only on the vertical $z$-coordinate. As a consequence, all 
horizontal gradients have to be taken from observations, or they have
to be estimated, parameterised or neglected.

For example, the surface slopes $\partial_x\zeta$ and
$\partial_y\zeta$ representing the barotropic pressure-gradients may
be determined by means of local observations or results from
three-dimensional numerical models. It is also possible to prescribe a
time series of the near-bed velocity components for reconstructing the
barotropic pressure gradient, see \cite{Burchard99}.  The
implementation of these options for the external pressure gradient is
carried out in {\tt extpressure.F90}, described in
\sect{sec:extpressure}.  The internal pressure-gradient, which results
from horizontal density gradients, can be prescribed from observations
of horizontal gradients of $\Theta$ and $S$ or from three-dimensional
model results (see {\tt intpressure.F90} in \sect{sec:intpressure}).
These gradients may also be used for horizontally advecting $\Theta$
and $S$ (see \sect{sec:temperature} and \sect{sec:salinity}).

Another option in GOTM for parameterising the advection of $\Theta$
and $S$ is to relax the model results to observations. Evidently, this
raises questions about the physical consistency of the model, but it
might help to provide a more realistic density field for studies of
turbulence dynamics.  Nudging is also possible for the horizontal
velocity components.  This makes sense in order to initialise inertial
oscillations from observed velocity profiles, see \sect{sec:uequation}
and \sect{sec:vequation}.  In the momentum equations, advection and
horizontal diffusion terms are neglected.

In hydrostatic ocean models, the vertical velocity is calculated by
means of the continuity equation, where the horizontal gradients of
$U$ and $V$ are needed. Since these are not available or set to zero,
the assumption of zero vertical velocity would be consistent.  In many
applications however, a non-zero vertical velocity is needed in order
to reflect the vertical adiabatic motion of e.g.\ a thermocline.  In
GOTM, we have thus included the option of prescribing a vertical
velocity time series at one height level which might be vertically
moving.  Vertical velocities at the surface and at the bottom are
prescribed according to the kinematic boundary conditions ($w=0$ at
the bottom and $w=\partial_t\zeta$ at the surface), and between these
locations and the prescribed vertical velocity at a certain height,
linear interpolation is applied, see {\tt updategrid.F90} in
\sect{sec:updategrid}.  This vertical velocity is then used for the
vertical advection of all prognostic quantities.

Standard relations according to the law of the wall are used for
deriving bottom boundary conditions for the momentum equations (see
{\tt friction.F90} in \sect{sec:friction}). At the sea surface, they
have to be prescribed or calculated from meteorological observations
with the aid of bulk formulae using the simulated or observed sea
surface temperature (see \sect{sec:airsea}). In {\tt
stratification.F90} described in \sect{sec:stratification}, the
buoyancy $b$ as defined in equation \eq{DefBuoyancy} is calculated by
means of the UNESCO equation of state (\cite{FofonoffMillard83})
or its linearised version. In
special cases, the buoyancy may also be calculated from a simple
transport equation. {\tt stratification.F90} is also used for
calculating the Brunt-V\"ais\"al\"a frequency, $N$.

The turbulent fluxes are calculated by means of various different
turbulence closure models described in great detail in the {\tt
turbulence} module, see \sect{sec:turbulence}.  As a simplifying
alternative, mixing can be computed according to the so-called
`convective adjustment' algorithm, see \sect{sec:convective}.

Furthermore, the vertical grid is also defined in the meanflow module
(see {\tt updategrid.F90} in \sect{sec:updategrid}). Choices for the
numerical grid are so-called $\sigma$-coordinates with layers heights
having a fixed portion of the water depth throughout the
simulation. Equidistant and non-equidistant grids are possible. 

\subsubsection{Numerics}\label{SectionNumericsMean}

For the spatial discretisation, the water column is divided into $N_i$
layers of not necessarily equal thickness $h_i$,
\begin{equation}
  \label{grid}
   h_i=(\gamma_i-\gamma_{i-1})D, \qquad i=1,\dots,N_i
   \comma
\end{equation}
with nondimensional interfaces $\gamma_i$ with $\gamma_0=-1$,
$\gamma_{i-1}< \gamma_i$ and $\gamma_{N_i}=0$,
see \cite{BurchardPetersen97}.

The discrete values for the mean flow quantities $U$, $V$, $\Theta$,
and $S$ represent interval means and are therefore located at the
centres of the intervals, and the turbulent quantities like $k$, $L$,
$\epsilon$, $\nu_t$, $\nu'_t$, $N$, $P$, $G$, $c_{\mu}$, and
$c_{\mu}'$ are positioned at the interfaces of the intervals (see
\sect{sec:turbulence}).  The indexing is such, that the interface
above an interval has the same index as the interval itself.  This
means that mean flow quantities range from $i=1,..,N_i$ while
turbulent quantities range from $i=0,..,N_i$ (see \fig{FigGrid}).
\begin{figure}[!h]	
  \begin{center}	
    \scalebox{0.5}{\includegraphics{figures/gridvert.eps}}
    \caption{Spatial organisation and indexing of the numerical grid.\label{FigGrid}}
  \end{center}
\end{figure}
The staggering of the grid allows for a
straight-forward discretisation of the vertical fluxes of momentum
and tracers without averaging.  However, for the vertical fluxes of
e.g.\ $k$ and $\epsilon$, averaging of the eddy diffusivities is
necessary. This is only problematic for the fluxes near the surface
and the bottom, where viscosities at the boundaries have to be
considered for the averaging. These can however be derived from the
law of the wall.
\begin{figure}
  \begin{center}
   \scalebox{0.5}{\includegraphics{figures/gridtime.eps}}
   \caption{Temporal organisation and indexing of the numerical grid.
            Here, a time stepping slightly more implicit than the 
            \index{Crank-Nicolson scheme} \cite{CrankNicolson47} 
            scheme with $\sigma=0.6$ is shown.\label{FigGridTime}}
\end{center}
\end{figure}

The time stepping is equidistant, based on two time levels and not
limited by Courant numbers, because of the absence of advection and an
implicit treatment of vertical diffusion, see \fig{FigGridTime}.  In
the following, the discretisation of a simple diffusion equation,
\begin{equation}
  \label{simpleDiffusion}
  \partder{X}{t} - \partder{}{z} \left( \nu \partder{X}{z} \right) = 0
  \comma
\end{equation}
will be illustrated for Neumann-type
boundary conditions
\begin{equation}
  \nu \partder{X}{z} = F_s
  \qquad \mbox{for } z=\zeta,\qquad
\end{equation}
and
\begin{equation}
  \nu \partder{X}{z} = F_b
  \qquad \mbox{for } z=-H.\qquad
\end{equation}

The semi-implicit discretisation of \eq{simpleDiffusion}
can then be written as
 \begin{equation}\label{sigmafirst}
 \displaystyle
 \frac{X^{n+1}_{N_i}-X^n_{N_i}}{\Delta t}
 -\frac{F_s
 -\nu^n_{N_i-1}\frac{X^{n+\sigma}_{N_i}-X^{n+\sigma}_{N_i-1}}{0.5(h^{n+1}_{N_i}+h^{n+1}_{N_i-1})}}{h^{n
 +1}_{N_i}}
 =
  \comma
 \end{equation}

\begin{equation}\label{Xdiscrete}
\displaystyle
\frac{X^{n+1}_i-X^n_i}{\Delta t}
-\frac{\nu^n_i\frac{X^{n+\sigma}_{i+1}-X^{n+\sigma}_{i}}{0.5(h^{n+1}_{i+1}+h^{n+1}_i)}
-\nu^n_{i-1}\frac{X^{n+\sigma}_{i}-X^{n+\sigma}_{i-1}}{0.5(h^{n+1}_i+h^{n+1}_{i-1})}}{h^{n
+1}_i}
=0
\comma
\end{equation}

\begin{equation}\label{sigmalast}
\displaystyle
\frac{X^{n+1}_1-X^n_1}{\Delta t}
-\frac{\nu^n_1\frac{X^{n+\sigma}_{2}-X^{n+\sigma}_{1}}{0.5(h^{n+1}_{2}+h^{n+1}_1)}
-F_b}{h^{n+1}_1}
=0
\comma
\end{equation}
for $1<i<N_i$. Here, the semi-implicit time level is defined by
\begin{equation}
  X^{n+\sigma}=\sigma X^{n+1}+(1-\sigma)X^n.
\end{equation}
Thus, for $\sigma=0$, a fully explicit, for $\sigma=1$ a fully
implicit, and for $\sigma=0.5$ the \cite{CrankNicolson47}
second-order scheme are obtained. \Fig{FigGridTime} shows an
example for $\sigma=0.6$. It should be noted that often a time
stepping is preferable which is slightly more implicit than the
\cite{CrankNicolson47} scheme in order to obtain
asymptotic stability.  The resulting linear system of equations
(\ref{sigmafirst}) -- (\ref{sigmalast}) with tri-diagonal matrix
structure is solved by means of the simplified Gaussian elimination. 

With the same strategy, a very similar system of equations can be
derived for variables located at the interfaces of the grid cells,
i.e. variables describing turbulence.

