The hydrostatic primitive equations (HPE’s) in p-coordinates are:
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}}_{h}+\mathbf{\nabla }_{p}\phi = \vec{\mathbf{\mathcal{F}}}
\frac{\partial \phi }{\partial p}+\alpha = 0
\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \omega }{\partial p} = 0
p\alpha = RT
c_{v}\frac{DT}{Dt}+p\frac{D\alpha }{Dt} = \mathcal{Q}
where \vec{\mathbf{v}}_{h}=(u,v,0) is the ‘horizontal’ (on pressure surfaces) component of velocity, \frac{D}{Dt}=\frac{\partial}{\partial t}+\vec{\mathbf{v}}_{h}\cdot \mathbf{\nabla }_{p}+\omega \frac{\partial }{\partial p} is the total derivative, f=2\Omega \sin \varphi is the Coriolis parameter, \phi =gz is the geopotential, \alpha =1/\rho is the specific volume, \omega =\frac{Dp }{Dt} is the vertical velocity in the p-coordinate. Equation :eq:`atmos-heat` is the first law of thermodynamics where internal energy e=c_{v}T, T is temperature, Q is the rate of heating per unit mass and p\frac{D\alpha }{Dt} is the work done by the fluid in compressing.
It is convenient to cast the heat equation in terms of potential temperature \theta so that it looks more like a generic conservation law. Differentiating :eq:`atmos-eos` we get:
p\frac{D\alpha }{Dt}+\alpha \frac{Dp}{Dt}=R\frac{DT}{Dt}
which, when added to the heat equation :eq:`atmos-heat` and using c_{p}=c_{v}+R, gives:
c_{p}\frac{DT}{Dt}-\alpha \frac{Dp}{Dt}=\mathcal{Q}
Potential temperature is defined:
\theta =T(\frac{p_{c}}{p})^{\kappa }
where p_{c} is a reference pressure and \kappa =R/c_{p}. For convenience we will make use of the Exner function \Pi (p) which is defined by:
\Pi (p)=c_{p}(\frac{p}{p_{c}})^{\kappa }
The following relations will be useful and are easily expressed in terms of the Exner function:
c_{p}T=\Pi \theta \;\;;\;\;\frac{\partial \Pi }{\partial p}=\frac{\kappa \Pi
}{p}\;\;;\;\;\alpha =\frac{\kappa \Pi \theta }{p}=\frac{\partial \ \Pi }{
\partial p}\theta \;\;;\;\;\frac{D\Pi }{Dt}=\frac{\partial \Pi }{\partial p}
\frac{Dp}{Dt}
where b=\frac{\partial \ \Pi }{\partial p}\theta is the buoyancy.
The heat equation is obtained by noting that
c_{p}\frac{DT}{Dt}=\frac{D(\Pi \theta )}{Dt}=\Pi \frac{D\theta }{Dt}+\theta
\frac{D\Pi }{Dt}=\Pi \frac{D\theta }{Dt}+\alpha \frac{Dp}{Dt}
and on substituting into :eq:`eq-p-heat-interim` gives:
\Pi \frac{D\theta }{Dt}=\mathcal{Q}
which is in conservative form.
For convenience in the model we prefer to step forward :eq:`potential-temperature-equation` rather than :eq:`atmos-heat`.
The upper and lower boundary conditions are:
\begin{aligned}\mbox{at the top:}\;\;p=0 &\text{, }\omega =\frac{Dp}{Dt}=0\end{aligned}
\begin{aligned}\mbox{at the surface:}\;\;p=p_{s} &\text{, }\phi =\phi _{topo}=g~Z_{topo}\end{aligned}
In p-coordinates, the upper boundary acts like a solid boundary (\omega=0 ); in z-coordinates the lower boundary is analogous to a free surface (\phi is imposed and \omega \neq 0).
For the purposes of initialization and reducing round-off errors, the model deals with perturbations from reference (or ‘standard’) profiles. For example, the hydrostatic geopotential associated with the resting atmosphere is not dynamically relevant and can therefore be subtracted from the equations. The equations written in terms of perturbations are obtained by substituting the following definitions into the previous model equations:
\theta = \theta _{o}+\theta ^{\prime }
\alpha = \alpha _{o}+\alpha ^{\prime }
\phi = \phi _{o}+\phi ^{\prime }
The reference state (indicated by subscript ‘o’) corresponds to horizontally homogeneous atmosphere at rest (\theta _{o},\alpha _{o},\phi_{o}) with surface pressure p_{o}(x,y) that satisfies \phi_{o}(p_{o})=g~Z_{topo}, defined:
\theta _{o}(p) = f^{n}(p) \\
\alpha _{o}(p) = \Pi _{p}\theta _{o} \\
\phi _{o}(p) = \phi _{topo}-\int_{p_{0}}^{p}\alpha _{o}dp
The final form of the HPE’s in p-coordinates is then:
\frac{D\vec{\mathbf{v}}_{h}}{Dt}+f\hat{\mathbf{k}}\times \vec{\mathbf{v}}
_{h}+\mathbf{\nabla }_{p}\phi ^{\prime } = \vec{\mathbf{\mathcal{F}}}
\frac{\partial \phi ^{\prime }}{\partial p}+\alpha ^{\prime } = 0
\mathbf{\nabla }_{p}\cdot \vec{\mathbf{v}}_{h}+\frac{\partial \omega }{
\partial p} = 0
\frac{\partial \Pi }{\partial p}\theta ^{\prime } = \alpha ^{\prime }
\frac{D\theta }{Dt} = \frac{\mathcal{Q}}{\Pi }