% \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}.\\