PlaSim/Plasim_RM_16/atmos.tex

% \section{Primitive Equation Model of the Atmosphere \label{PUMA}}
The primitive equations, which represent the dynamical core of the 
atmospheric model,
consist of the conservation of momentum and mass, the first law of 
thermodynamics and the equation of state, simplified by the hydrostatic 
approximation. 


\section{A dimensionless set of differential equations \label{EQ}}
The prognostic equations for the horizontal velocities are 
transformed into equations of the vertical component of the vorticity $\zeta$ 
and the divergence $D$. 
A vertical coordinate system where the lower boundary
exactly coincides with a coordinate surface is defined
by $\sigma$ (the pressure normalized by 
the surface pressure).
Latitude $\varphi$ and longitude $\lambda$ represent the horizontal 
coordinates and the poleward convergence of the meridians is explicitly 
introduced re-writing the zonal ($u$) and meridional ($\nu$) velocities: 
$U= u \cos \varphi$ , $V= \nu \cos \varphi$ and $\mu = \sin \varphi$.
The implicitly treated gravity wave terms
are linearized about a reference profile $T_0$.
Therefore, prognostic equation for
temperature deviations $T'=T-T_0$ are derived; 
we use a constant reference temperature $T_0=250K$ for all
$\sigma$ levels.
The turbulent flux divergences due to prior Reynolds averaging
enter the dynamic and thermodynamic equations as parameterizations formally 
included in the terms: $P_\zeta , P_D , P_T$.\\


A dimensionless set of differential equations is derived
by scaling vorticity $\zeta$ and divergence $D$ 
by angular velocity of the earth 
$\Omega$, pressure $p$ by a 
constant surface
pressure $p_s$, temperatures $T$ and $T'$ by $a^2 \Omega^2 /R$ and 
the orography and geopotential $\psi$
by $a^2 \Omega^2 / g$ ($g$ is the acceleration of gravity
and $R$ the gas constant for dry air). 
The dimensionless primitive
equations in the $(\lambda , \mu, \sigma)$-coordinates
\cite{hoskins} are given by\\


Conservation of momentum (vorticity  and divergence equation)

\begin{equation}
{\displaystyle \frac{\partial \zeta + f}{\partial t} = \frac{1}{(1 - \mu^2)} \frac{\partial F_\nu}{\partial \lambda} - \frac{\partial F_u}{\partial \mu} + P_\zeta} \label{vort}
\end{equation}

\begin{equation}
{\displaystyle \frac{\partial D}{\partial t} = \frac{1}{(1 - \mu^2)} \frac{\partial F_u}{\partial \lambda} + \frac{\partial F_\nu}{\partial \mu} - \bigtriangledown^2 E - \bigtriangledown^2 (\phi + T_0 \ln p_s ) + P_D}\label{div}
\end{equation}


Hydrostatic approximation (using the equation of state)

\begin{equation}
{\displaystyle 0 = \frac{\partial \phi}{\partial \ln \sigma} + T} \label{hydroute}
\end{equation}


Conservation of mass (continuity equation)

\begin{equation}
{\displaystyle \frac{\partial \ln p_s}{\partial t} 
= - \int\limits_{0}^{1} A d \sigma}\label{konti}
\end{equation}


Thermodynamic equation

\begin{equation}
{\displaystyle \frac{\partial T'}{\partial t} = F_T - \dot{\sigma}
  \frac{\partial T}{\partial \sigma} + \kappa W T   + \frac{J}{c_p} + P_{T}}
\label{temp} 
\end{equation}


with the notations\\


${\displaystyle F_u = ( \zeta + f ) V - \dot{\sigma} \frac{\partial U}{\partial \sigma} - T' \frac{\partial \ln p_s}{\partial \lambda}}$


${\displaystyle F_\nu = - (\zeta + f)U - \dot{\sigma} \frac{\partial V}{\partial\sigma} - (1 - \mu^2) T' \frac{\partial \ln p_s}{\partial \mu}}$



${\displaystyle F_T =  - \frac{1}{(1 - \mu^2)} \frac{\partial (U T')}{\partial \lambda} - \frac{\partial (VT')}{\partial \mu} + D T'} $

${\displaystyle E = \frac{U^2 + V^2}{2(1 - \mu^2)}}$


${\displaystyle
\dot {\sigma}=  \sigma \int\limits_{0}^{1} A d \sigma
-  \int\limits_{0}^{\sigma} A d \sigma  }$\\


${\displaystyle
W= \frac {\omega} {p} =
\vec{V} \cdot  \nabla \ln p_s - \frac{1}{\sigma}\int\limits_{0}^{\sigma} A d
\sigma  }$\\

$A=D+\vec{V} \cdot \nabla \ln p_s 
= \frac{1}{p_s} \nabla \cdot  p_s  \vec{V}$.\\


Here is $\dot{\sigma}$ the 
vertical velocity in the $\sigma$ system, $J$ the diabatic heating 
per unit mass and $E$  
the kinetic energy per unit mass.
The streamfunction $\psi$ and the velocity potential $\chi$ represent the 
nondivergent and the irrotational part of the velocity field\\
${\displaystyle U = -(1-\mu^2) \frac{\partial \psi}{\partial \mu} + \frac{\partial \chi}{\partial \lambda}}$
and 
${\displaystyle V = \frac{\partial \psi}{\partial \lambda} + (1-\mu^2) \frac{\partial \chi}{\partial \mu}}$
with $ \zeta = \nabla^2 \psi $ 
and 
$ D =  \nabla^2 \chi $.

% \section{Moisture \label{Moist}}
%  ** Frank **

\section{Mode splitting \label{Mode}}


The fast gravity wave modes are linearized around
a reference temperature profile $\vec{T_0}$.
Now, the differential equations (\ref{vort}-\ref{temp}) can be separated into 
fast (linear) gravity modes and
the slower non-linear terms ($N_D, N_p, N_T $).  
The linear terms of the equations contain the effect of the 
divergence (or the gravity waves) on 
the surface pressure tendency, the temperature tendency and the
geopotential.
A discussion of the impact of the reference profile
on the stability of the semi-implicit numerical scheme is
presented by \cite{simmons1978}.

\begin{equation}
{\displaystyle\frac{ \partial  D }{\partial t}= 
{ N_D} - \bigtriangledown^2 (\phi + T_0 \ln p_s)}\label{Mdiv}
\end{equation}


\begin{equation}
{\displaystyle \frac{\partial \phi}{\partial \ln \sigma} = - T}\label{Mhydro}
\end{equation}

\begin{equation}
{\displaystyle \frac{\partial \ln p_s}{\partial t} 
= N_p - \int\limits_{0}^{1} D d \sigma}\label{Mkonti}
\end{equation}

\begin{equation}
{\displaystyle \frac{\partial T'}{\partial t} = N_T- 
\dot {\sigma}_L
  \frac{\partial T_0}{\partial \sigma} + \kappa W_L T_0 } \label{Mtemp}
\end{equation}\\




with the non-linear terms\\




${\displaystyle 
{\displaystyle N_D = \frac{1}{(1 - \mu^2)} \frac{\partial F_u}{\partial \lambda} + \frac{F_\nu}{\partial \mu} - \bigtriangledown^2 E  + P_D}
}$

${\displaystyle 
{\displaystyle N_p 
= - \int\limits_{0}^{1} [A-D] d \sigma}
}$

${\displaystyle 
{\displaystyle N_T = F_T 
- \dot{\sigma}_N \frac{\partial T_0}{\partial \sigma} 
- \dot{\sigma} \frac{\partial T'}{\partial \sigma} 
+ \kappa W_N T_0   + \kappa W T'   
+ \frac{J}{c_p} + P_{T}}
}$\\


and the notations\\

${\displaystyle
\dot {\sigma}_L =  \sigma \int\limits_{0}^{1} D d \sigma
-  \int\limits_{0}^{\sigma} D d \sigma  }$\\

${\displaystyle
\dot {\sigma}_N =  \sigma \int\limits_{0}^{1} [A-D] d \sigma
-  \int\limits_{0}^{\sigma}  [A-D] d \sigma  }$\\


${\displaystyle
W_L  =  - \frac{1}{\sigma} \int\limits_{0}^{\sigma} D d \sigma  }$\\

${\displaystyle
W_N  =
\vec{V} \cdot  \nabla \ln p_s - \frac{1}{\sigma}\int\limits_{0}^{\sigma} 
[A  - D] d\sigma  }$\\

$A-D=\vec{V} \cdot \nabla \ln p_s$ \\


${\displaystyle
\dot {\sigma}=\dot {\sigma_L} + \dot {\sigma_N}
=  \sigma \int\limits_{0}^{1} A d \sigma
-  \int\limits_{0}^{\sigma} A d \sigma  }$\\


${\displaystyle
W= W_L+W_N =
\vec{V} \cdot  \nabla \ln p_s - \frac{1}{\sigma}\int\limits_{0}^{\sigma} A d
\sigma  }$\\



%A combination of the linear terms, the matrix $ {\cal B} = {\cal L}_{\phi} {\cal L}_T + \vec{T}_0 \vec{L}_p = {\cal B}(\sigma , \kappa , \vec{T}_0)$ is constant in time.

The index $L$ denote the linear and $N$ the non-linear part in  
the vertical advection ($\dot {\sigma} \frac{\partial T}{\partial \sigma}$) 
and the adiabatic
heating or cooling ($\kappa W T$ with $W= \frac {\omega} {p}$).
The non-linear terms are solve explicitly in the physical space (on the
Gaussian grid; section \ref{Trans}) and the linear terms
are calculated implicitly in the spectral space (for the spherical harmonics; 
see section \ref{Trans}).



\section{Numerics \label{Num}}

Solving the equations requires a suitable numerical representation of the 
spatial fields and their time change. A conventional approach is 
spectral representation in the horizontal using the
transform method, finite differences in the vertical, 
and a semi-implicit time stepping. 
 

\subsection{Spectral Transform method \label{Trans}}
The spectral method used in the computation of the nonlinear terms 
involves storing of a large number of so-called
interaction coefficients, the number of which
increases very fast with increasing resolution.
The computing time and storing
space requirements exceed all practical limits for high
resolution models. Furthermore, there are problems to incorporate 
locally dependent physical
processes, such as release of precipitation or a convective 
adjustment.
Therefore, the equations are solved using the spectral 
transform method \cite{orszag1970, eliassen1970}. 
This method uses an auxiliary grid in the physical space
where point values of the dependent variables are computed.\\


The prognostic variables are represented in the horizontal by truncated series of spherical harmonics
($Q$ stands for $\zeta, D,T$ and $ \ln p_s$)
\begin{equation}
{\displaystyle Q(\lambda , \mu , \sigma , t) = \sum_{m = - M}^{M} \sum_{n =
    |m|}^{M} Q_n^m (\sigma , t) P_n^m (\mu) e ^{im\lambda} }
\end{equation}

\begin{center}
$=  Q_n^0 (\sigma , t) P_n^0 (\mu) 
+ 2   \sum_{m = 1}^{M} \sum_{n =
    m}^{M} Q_n^m (\sigma , t) P_n^m (\mu) e ^{im\lambda}$ 
\end{center}


For each variable the spectral coefficient is defined by
\begin{equation}
{\displaystyle Q_n^m (\sigma , t) = \frac{1}{4 \pi} \int\limits_{-1}^{1}
  \int\limits_{0}^{2 \pi} Q(\lambda , \mu , \sigma , t) P_n^m (\mu) e^{ -i m
    \lambda} d \lambda d \mu}
\end{equation}


The spectral coefficients $Q_n^m (\sigma , t)$
are obtained by Gaussian quadrature 
of the Fourier coefficients $F^m$ at each latitude $\varphi$ which 
are calculated by Fast Fourier 
Transformation with \\
${\displaystyle F^m (\mu, \sigma , t) = \frac{1}{4 \pi}  \int\limits_{0}^{2 \pi} Q(\lambda , \mu , \sigma , t) e^{ -i m \lambda} d \lambda}$\\
The auxiliary grid in the physical space (Gaussian grid)
is defined by $M_g$ equally spaced longitudes 
and $J_g$ Gaussian latitudes with
$M_g \ge 3 M + 1$ and $J_g  \ge 0.5 (3 M + 1)$.



\subsection{Vertical discretization} 
The prognostic variables vorticity, temperature and 
divergence are calculated at full levels and the vertical velocity at half 
levels. Therefore, the vertical advection for the level r  is calculated
($Q$ stands for $\zeta, D,T$ and $ \ln p_s$)

\begin{equation}
{\displaystyle ( \dot{\sigma} \frac{\partial Q}{\partial \sigma} ) 
\hat{=} 
\frac{1}{2\Delta \sigma_r } [\dot{\sigma}_{r + 0.5} (Q_{r+1} - Q_r) + \dot{\sigma}_{r - 0.5} (Q_r - Q_{r-1})]} 
\end{equation}


For the hydrostatic approximation (3)
an angular momentum
conserving finite-difference scheme 
\cite{simmons1981}
is used which solves the equation at half levels 
($r+0.5 ;r=1,...,n; n=$ number of levels)
\begin{equation}
\displaystyle  
\frac{\partial \phi}{\partial \ln \sigma}+T
\hat{=} 
\phi_{r+0.5}-\phi_{r-0.5}+
 T_r \cdot \ln \frac{\sigma_{r+0.5}}{\sigma_{r-0.5}}
\end{equation}

Full level values ($r$) 
of geopotential are given by
\begin{equation}
 \phi_{r}=\phi_{r+0.5}+\alpha_r T_r 
\end{equation}
with
${\displaystyle  \alpha_r=1-\frac{\sigma_{r-0.5}}{\Delta \sigma_r}
 \ln \frac{\sigma_{r+0.5}}{\sigma_{r-0.5}}}$ 
and
$ \Delta \sigma_r=\sigma_{r+0.5} - \sigma_{r-0.5}$\\


%The analogues for the integral formulation in (4) and (5)are given by \\

%$ \int\limits_{0}^{1} D d \sigma $  and 
%$ \int\limits_{0}^{\sigma} D d \sigma $ 
%$\rightarrow $
%${\displaystyle  \sum_{k=1}^{n} D_k \Delta \sigma_k  }$ and
%${\displaystyle  \sum_{k=1}^{r} D_k \Delta \sigma_k  }$
%with $ \Delta \sigma_k=\sigma_{k+0.5} - \sigma_{k-0.5}$\\



\subsection{Semi-implicit time stepping}


Sound waves are filtered by the hydrostatic 
approximation (filter for vertical sound waves) and the lower boundary 
condition in pressure or sigma-coordinates 
(vanishing vertical velocity at the surface, i.e. 
the total derivative of the surface pressure is zero; 
filter for horizontal sound waves).
But the fast propagation of the gravity 
waves strongly reduce the time step of explicit
numerical schemes, therefore mode splitting is used
(section \ref{Mode}) and an implicit scheme for the
divergence is applied (see below). 
The vorticity equation is computed by an explicit scheme (leap frog) and the 
common Robert/Asselin time filter is used \cite{haltiner1982}.\\


The implicit formulation for the divergence is derived 
using the conservation of mass, the hydrostatic approximation
and the thermodynamic equation (eq. \ref{Mdiv}-\ref{Mtemp}) approximated 
by its finite difference analogues in time ($t$) using the
notation (for each variable $D$, $T$, $\ln p_s$, and $\phi$)\\

${\displaystyle \delta_t Q = \frac{Q^{t + \Delta t} - Q^{t - \Delta t}}{2
\Delta t}}$ \hspace{0.8cm}
and \hspace{0.8cm}
${\displaystyle \overline{Q}^t = 0.5 (Q^{t + \Delta t} + Q^{t - \Delta t})
  =Q^{t - \Delta t} + \Delta t \delta_t Q}  $\\



The divergence is calculated by the 
non-linear term at time step $t$ and the 
linearized term which is a function of the geopotential
(or the temperature tendency) and  
the surface pressure tendency.

\begin{equation}
\delta_t { D} = { N_D}^t - \bigtriangledown^2 (\overline{\phi}^t + T_0 [\ln
p_s^{t - \Delta t} + \Delta t \; \delta_t \ln p_s])
\end{equation}

\begin{equation}
\overline{\phi-\phi_s}^t = L_{\phi} [T^{t - \Delta t} + \Delta t \delta_t T]
= L_{\phi} [T^{t - \Delta t} + \Delta t \; \delta_t T']
\end{equation}

\begin{equation}
\delta_t \ln p_s = {N_p}^t - L_p [D^{t - \Delta t} + \Delta t \; \delta_t D]
\end{equation}

\begin{equation}
\delta_t  T' =  {N_T}^t - L_T [D^{t - \Delta t} + \Delta t \; \delta_t D]
\end{equation}


The implicit formulation of the divergence equation is derived from the
finite difference analogues of the new time step 
$t+\Delta t$ applied for each level $ r$ ($r=1,...n$) which can also
formulated as a vector $\vec{D}$ with the n components.\\


${\left(\begin{array}{*{4}{c}}
            1-b_{11} & b_{21}   & \cdots & b_{n1}   \\
            b_{12}   & 1-b_{22} & \ddots & \vdots   \\
            \vdots   & \vdots   & \ddots & \vdots   \\
            b_{1n}   & b_{2n}   & \cdots & 1-b_{nn} \\
            \end{array}
            \right)
\left(\begin{array}{*{1}{c}}
            D_{1}^{t+\Delta t}   \\
            D_{2}^{t+\Delta t}     \\
            \vdots   \\
            D_{n}^{t+\Delta t}    \\
            \end{array}
            \right)
=
\left(\begin{array}{*{1}{c}}
            D_{1}^{t-\Delta t}    \\
            D_{2}^{t-\Delta t}     \\
            \vdots   \\
            D_{n}^{t-\Delta t}      \\
            \end{array}
            \right)
+2 \Delta t
\left(\begin{array}{*{1}{c}}
            R_{1}    \\
            R_{2}     \\
            \vdots     \\
            R_{n}     \\
            \end{array}
            \right)}$\\


In matrix formulation\\

$({\cal I} - {\cal B} \Delta t^2 \bigtriangledown^2) \vec{D}^{t + \Delta t} =
\vec{D}^{t-\Delta t}
 + 2 \Delta t [ \vec{N}_D - \bigtriangledown^2 
(\vec{\phi}^{t - \Delta t} + \vec{T}_0 \ln p_s^{t- \Delta t})]$ \\
\begin{equation}
 - 2 \Delta t^2 \bigtriangledown^2 ({\cal L}_{\phi} \vec{N}_T + \vec{T}_0 N_p )
\end{equation}\\

The matrix $ {\cal B} = {\cal L}_{\phi} {\cal L}_T + \vec{T}_0 \vec{L}_p = 
{\cal B}(\sigma , \kappa , \vec{T}_0)$ 
is constant in time. 
The variables ${\vec{D},\vec{T},\vec{T}',\vec{\phi}-\vec{\phi}_s}$
are represented by column vectors with values
at each layer, as are also $\vec{N}_D$ and $\vec{N}_T$.
${\cal L}_{\phi}$ and ${\cal L}_T$ are constant matrices, $\vec{L}_p$ 
is a row vector (see Appendix C). 
The matrix  $ {\cal B}$ can be calculated seperately for
each spectral coefficient because in the linearized part
the spectral modes are independent of each other.
\begin{equation}
(\vec{D}_n^{m})^{t + \Delta t} =
({\cal I} + {\cal B} \Delta t^2 c_n)^{-1} 
[(\vec{D}_n^{m})^{t - \Delta t}
+ 2 \Delta t \vec{R}]
\end{equation}

\begin{equation}
(\vec{D}_n^{m})^{t + \Delta t} =
(\frac{{\cal I}}{c_n} + {\cal B} \Delta t^2)^{-1} 
[\frac{1}{c_n}(\vec{D}_n^{m})^{t - \Delta t}
+ \frac{2 \Delta t}{c_n} \vec{R}]
\end{equation}\\


with $\bigtriangledown^2 (P_n^m (\mu) e^{ -i m \lambda})=- n(n+1) P_n^m (\mu) e^{
  -i m \lambda} =- c_n P_n^m (\mu) e^{ -i m \lambda} $.
% Numerical Recipes routines
% (lubksb and ludcmp) are used for matrix manipulation \cite [] {Press:1992}.\\