\section{Pressure coordinate}

The primitive equations in the $(\lambda , \mu, p)$ -coordinates
without scaling. That means $D$ and $zeta$ in Appendix A and B have
the units: $s^-1$, $T$ is in $K$, $p$ in $Pa$,  $\phi$ in $m^2 s^{-2}$
and $\vec{\nu}$ in $m s^{-1}$.\\


Conservation of momentum (vorticity and divergence equation)

\begin{equation}
{\displaystyle \frac{\partial \zeta}{\partial t} = - \vec{\nu} \cdot \bigtriangledown (\zeta + f) - \omega \frac{\partial \zeta}{\partial p} - (\zeta + f) \bigtriangledown \cdot \vec{\nu} + \vec{k} \cdot (\frac{\partial \vec{\nu}}{\partial p} \times \bigtriangledown \omega) + P_\zeta} 
\end{equation}       

\begin{equation}
{\displaystyle \frac{\partial D}{\partial t} = \vec{k} \cdot \bigtriangledown \times (\zeta + f) \vec{\nu} - \bigtriangledown \cdot (\omega \frac{\partial \vec{\nu}}{\partial p}) - \bigtriangledown^2 (\phi + \frac{\vec{\nu}^2}{2}) + P_D} \end{equation}  

Hydrostatic approximation (using the equation of state)

\begin{equation}
{\displaystyle \frac{\partial \phi}{\partial p} = - \frac{1}{\rho} = - \frac{RT}{p}}     
\end{equation}     

Conservation of mass (continuity equation)

\begin{equation}
{\displaystyle \bigtriangledown \cdot \vec{\nu} + \frac{ \partial \omega}{\partial p} = 0}     
\end{equation} 

Thermodynamic equation ( J= diabatic heating per unit mass)

\begin{equation}
{\displaystyle \frac{d T}{d t} = \frac{\omega}{c_p \rho} + \frac{J}{c_p} + P_T}
\end{equation} 
    

\section{Sigma-system}

$\sigma=p/p_s$ ranges monotonically from zero at 
the top of the atmosphere to unity at the ground. For $\xi=x,y $ or $t$
\begin{equation}
{\displaystyle (\frac{\partial }{\partial \xi})_p
=\frac{\partial }{\partial \xi}
-\sigma \frac{\partial \ln p_s}{\partial \xi}
\frac{\partial }{\partial \sigma}}
\end{equation}

\begin{equation}
\frac{\partial }{\partial p}
=\frac{\partial \sigma }{\partial p}
\frac{\partial }{\partial \sigma}
=\frac{1 }{p_s}
\frac{\partial }{\partial \sigma}
\end{equation}

The vertical velocity in the p-coordinate system $\omega$
and in the new $\sigma$-coordinate 
system $\dot {\sigma}$
are given by \cite{phillips1957}

\begin{equation}
\omega=  \frac{p}{p_s}
[\vec{V} \cdot  \nabla  p_s
 -  \int\limits_{0}^{\sigma} \nabla \cdot  p_s \vec{V} d \sigma]
=  p [\vec{V} \cdot  \nabla \ln p_s]
 - p_s \int\limits_{0}^{\sigma} A d \sigma
\end{equation}

\begin{equation}
\dot {\sigma}=  \sigma \int\limits_{0}^{1} A d \sigma
-  \int\limits_{0}^{\sigma} A d \sigma
\end{equation}

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





The primitive equations in the $(\lambda , \mu, \sigma)$ -coordinates
without scaling

Conservation of momentum (vorticity and divergence equation)

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

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

Hydrostatic approximation (using the equation of state)

\begin{equation}
{\displaystyle \frac{\partial \phi}{\partial \ln \sigma} = - TR} 
\end{equation}

Conservation of mass (continuity equation)

\begin{equation}
{\displaystyle \frac{\partial \ln p_s}{\partial t} = - \frac{U}{a (1 - \mu^2)} \frac{\partial \ln p_s}{\partial \lambda} - \frac{V}{a} \frac{\partial \ln p_s}{\partial \mu} - D - \frac{\partial \dot{\sigma}}{\partial \sigma}
= - \int\limits_{0}^{1} (D+\vec{V} \cdot \nabla \ln p_s) d \sigma}
\end{equation}

Thermodynamic equation ( J= diabatic heating per unit mass)

\begin{equation}
{\displaystyle \frac{\partial T}{\partial t} = F_T  - \dot{\sigma} \frac{\partial T}{\partial \sigma} + \kappa T [\vec{V} \cdot  \nabla \ln p_s - \frac{1}{\sigma}\int\limits_{0}^{\sigma} A d \sigma] +\frac{J}{c_p} + P_T}
\end{equation}



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

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

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

${\displaystyle F_T =  - \frac{U}{a(1-\mu^2)} \frac{\partial T}{\partial \lambda} - \frac{V}{a} \frac{\partial T}{\partial \mu} } $

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


\section{Matrix {\em B}}


For the implicit scheme, fast (linear) gravity modes and
the slower non-linear terms are separated. \\

${\displaystyle\frac{ \partial  D }{\partial t}= 
{ N_D} - \bigtriangledown^2 (\phi + T_0 \ln p_s)} $\\

${\displaystyle \frac{\partial \ln p_s}{\partial t} 
= N_p - \int\limits_{0}^{1} D d \sigma}$\\


${\displaystyle \frac{\partial T'}{\partial t} = N_T- [ \sigma
  \int\limits_{0}^{1} D d \sigma -  \int\limits_{0}^{\sigma} D d \sigma ]
  \frac{\partial T_0}{\partial \sigma} + \kappa T_0 [-
  \int\limits_{0}^{\sigma} D d \ln  \sigma] }$\\

${\displaystyle \frac{\partial \phi}{\partial \ln \sigma} = - T} $\\


The set of differential equations are approximated 
by its finite difference analogues using the
notation (for each variable $D$, $T$, $\ln p_s$, and $\phi$)\\

${\displaystyle \overline{Q}^t = 0.5 (Q^{t + \Delta t} + Q^{t - \Delta t})
  =Q^{t - \Delta t} + \Delta t \delta_t Q}  $

and

${\displaystyle \delta_t Q = \frac{Q^{t + \Delta t} - Q^{t - \Delta t}}{2 \Delta t}}$\\



The hydrostatic approximation using an angular momentum
conserving finite-difference scheme
is solved at half levels\\

${\displaystyle  \phi_{r+0.5}-\phi_{r-0.5}=
 T_r \cdot \ln \frac{\sigma_{r+0.5}}{\sigma_{r-0.5}}}$\\

Full level values of geopotential are given by\\

${\displaystyle  \phi_{r}=\phi_{r+0.5}+\alpha_r
 T_r }$
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}$\\


Now, the implicit formulation for the divergence is derived 
using the conservation of mass, the hydrostatic approximation
and the thermodynamic equation at discrete time steps\\

${\displaystyle \delta_t { D} = { N_D} - \bigtriangledown^2 (\overline{\phi}^t + T_0 [\ln p_s^{t - \Delta t} + \Delta t \delta_t \ln p_s])} $\\


${\displaystyle \delta_t \ln p_s = N_p - L_p [D^{t - \Delta t} + \Delta t \delta_t D]}$\\


${\displaystyle \overline{ \phi - \phi_s}^t = L_{\phi} [T^{t - \Delta t} + \Delta t \delta_t T}]$\\


${\displaystyle \delta_t  T' =  N_T - L_T [D^{t - \Delta t} + \Delta t \delta_t D]} $\\

The set of differential equations
for each level $ k (k=1,..,n)$  written in
vector form leads to the matrix $ {\cal B}$ with n rows and
n columns. 
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}$
$\vec{N}_D$ and $\vec{N}_T$ are represented by column vectors with values
at each level. $L_p$, $L_T$ and $L_{\phi}$ contain the effect of the 
divergence (or the gravity waves) on 
the surface pressure tendency, the temperature tendency and the
geopotential.\\
 
$\vec{L}_p =(\Delta \sigma_1, ..., \Delta \sigma_n)$
is a row vector with
$ \Delta \sigma_n=\sigma_{n+0.5} - \sigma_{n-0.5}$.\\


${\cal L}_{\phi}=
{\left(\begin{array}{*{5}{c}}
            1 & \alpha_{21} &\alpha_{31}& \cdots & \alpha_{n1}    \\
            0 &\alpha_{22}  &\alpha_{32}& \ddots & \vdots   \\
            \vdots & \vdots & \vdots    & \ddots & \vdots   \\
            0      & 0      &  \cdots   & 0      & \alpha_{nn} \\
            \end{array}
            \right)} $\\


For $i=j: {\displaystyle \alpha_{jj}=1- [
 \frac{\sigma_{j-0.5}}{\sigma_{j+0.5}-\sigma_{j-0.5}}
(\ln \sigma_{j+0.5} - \ln \sigma_{j-0.5})]}$\\

$i>j: \alpha_{ij}=\ln \sigma_{j+0.5} - \ln \sigma_{j-0.5}$\\

$i<j: \alpha_{ij}=0$.\\



${\cal L}_{T}
=  {\left(\begin{array}{*{4}{c}}
\kappa (T_0)_1 \alpha_{11}&\kappa (T_0)_1 \alpha_{21}&\cdots&\kappa (T_0)_1 \alpha_{n1} \\
\kappa (T_0)_2 \alpha_{12}&\kappa (T_0)_2 \alpha_{22}&\ddots&   \vdots \\
            \vdots   & \vdots   & \ddots & \vdots   \\
\kappa (T_0)_n \alpha_{1n}&\kappa (T_0)_n \alpha_{2n}&\cdots&\kappa (T_0)_n \alpha_{nn}    \\
            \end{array}
            \right)}
+ {\left(\begin{array}{*{5}{c}}
            \gamma_{11}   &\gamma_{21}& \gamma_{31} & \cdots & \gamma_{n1}  \\
            \gamma_{12}   &\gamma_{22}& \gamma_{32} & \ddots & \vdots   \\
            \vdots      & \vdots  & \vdots    & \ddots & \vdots   \\
            \gamma_{1n} & \gamma_{2n} & \cdots    & \cdots & \gamma_{nn} \\
            \end{array}
            \right)} $\\

$ \tau_{ij}=\kappa (T_0)_j \alpha_{ij}+\gamma_{ij}$
with
$ \Delta T_{n+0.5}=(T_0)_{n+1} - (T_0)_{n}$\\

for $j=1$ and\\
$i=j$:  $ \gamma_{jj}=  \frac{1}{2} [\Delta T_{0.5} (\sigma_{1}-1)  ] $ \\

$i>j$: $ \gamma_{ij}= \frac{1}{2} \Delta \sigma_{i} [ \Delta T_{0.5}  \sigma_{1}  ]  $ \\


for $j>1$ and\\
$i=j$:  $ \gamma_{jj}= \frac{1}{2} [\Delta T_{j-0.5} \sigma_{j-0.5}
+ \Delta T_{j+0.5} (\sigma_{j+0.5}-1)  ]  $ \\

$i<j$:  $ \gamma_{ij}= \frac{\Delta \sigma_{i} }{2 \Delta \sigma_{j} } 
[\Delta T_{j-0.5} (\sigma_{j-0.5}-1)
+ \Delta T_{j+0.5}  (\sigma_{j+0.5}-1)  ]   $  \\

$i>j$: $ \gamma_{ij}= \frac{\Delta \sigma_{i}}{2 \Delta \sigma_{j} } 
[\Delta T_{j-0.5}  \sigma_{j-0.5}
+ \Delta T_{j+0.5}  \sigma_{j+0.5}  ] $ \\






%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: