3.1 Eulerian Dynamical Core

The hybrid vertical coordinate that has been implemented in CAM 3.0 is described in this section. The hybrid coordinate was developed by Simmons and Strüfing [158] in order to provide a general framework for a vertical coordinate which is terrain following at the Earth's surface, but reduces to a pressure coordinate at some point above the surface. The hybrid coordinate is more general in concept than the modified $\sigma$ scheme of Sangster [155], which is used in the GFDL SKYHI model. However, the hybrid coordinate is normally specified in such a way that the two coordinates are identical.

The following description uses the same general development as Simmons and Strüfing [158], who based their development on the generalized vertical coordinate of Kasahara [84]. A specific form of the coordinate (the hybrid coordinate) is introduced at the latest possible point. The description here differs from Simmons and Strüfing [158] in allowing for an upper boundary at finite height (nonzero pressure), as in the original development by Kasahara. Such an upper boundary may be required when the equations are solved using vertical finite differences.

3.1.1 Generalized terrain-following vertical coordinates

Deriving the primitive equations in a generalized terrain-following vertical coordinate requires only that certain basic properties of the coordinate be specified. If the surface pressure is $\pi$, then we require the generalized coordinate $\eta(p,\pi)$ to satisfy:

1. $\eta(p,\pi)$ is a monotonic function of $p$.
2. $\eta(\pi,\pi)=1$
3. $\eta(0,\pi)=0$
4. $\eta(p_t,\pi)=\eta_t$ where $p_t$ is the top of the model.

The latter requirement provides that the top of the model will be a pressure surface, simplifying the specification of boundary conditions. In the case that $p_t=0$, the last two requirements are identical and the system reduces to that described in Simmons and Strüfing [158]. The boundary conditions that are required to close the system are:

$$\dot\eta(\pi,\pi) = 0,$$ (3.1)

$$\dot\eta(p_t,\pi) = \omega(p_t) = 0.$$ (3.2)

Given the above description of the coordinate, the continuous system of equations can be written following Kasahara [84] and Simmons and Strüfing [158]. The prognostic equations are:

$$\frac{\partial\zeta}{\partial t} = {\boldsymbol {k}}\cdot\nabla\times ({\boldsymbol {n}}/\cos\phi) + F_{\zeta_H} ,$$ (3.3)

$$\frac{\partial\delta}{\partial t} = \nabla\cdot ({\boldsymbol {n}}/\cos\phi) -\nabla^2\left(E+\Phi \right) + F_{\delta_H} ,$$ (3.4)

$$\frac{\partial T}{\partial t} = \frac{-1}{a\cos^2\phi} \left[\frac{\partial}{\partial\lambda} (UT... ...dot\eta \frac{\partial T}{\partial\eta} + \frac{R}{c_p^*}{T_v} \frac{\omega}{p}$$

$$\phantom{=} + Q + F_{T_H} + F_{F_H} ,$$ (3.5)

$$\frac{\partial q}{\partial t} = \frac{-1}{a\cos^2\phi} \left[ \frac{\partial}{\partial\lambda} (U... ...l\phi } (Vq) \right] + q\delta - \dot\eta \frac{\partial q}{\partial\eta} + S ,$$ (3.6)

$$\frac{\partial \pi}{\partial t} = \int_1^{\eta_t} {\boldsymbol{\nabla}\cdot} \left( \frac{\partial p}{\partial\eta} {\boldsymbol {V}} \right) d\eta .$$ (3.7)

The notation follows standard conventions, and the following terms have been introduced with ${\boldsymbol {n}} = (n_U,n_V)$:

$$n_U = + (\zeta+f)V - \dot\eta \frac{\partial U}{\partial\eta} R\frac{T_v}{p}\frac{1}{a} - \frac{\partial p}{\partial\lambda} + F_U ,$$ (3.8)

$$n_V = - (\zeta+f)U - \dot\eta \frac{\partial V}{\partial\eta} - R\frac{T_v}{p}\frac{\cos\phi}{a} \frac{\partial p}{\partial\phi} + F_V ,$$ (3.9)

$$E = \frac{U^2+V^2}{2\cos^2\phi} ,$$ (3.10)

$$(U,V) = (u,v)\cos\phi ,$$ (3.11)

$${T_v} = \left[ 1 + \left( \frac{R_v}{R} -1 \right) q \right] T ,$$ (3.12)

$$c_p^* = \left[ 1 + \left( \frac{c_{p_v}}{c_p} -1 \right) q \right] c_p .$$ (3.13)

The terms $F_U, F_V, Q$, and $S$ represent the sources and sinks from the parameterizations for momentum (in terms of $U$ and $V$), temperature, and moisture, respectively. The terms $F_{\zeta_H}$ and $F_{\delta_H}$ represent sources due to horizontal diffusion of momentum, while $F_{T_H}$ and $F_{F_H}$ represent sources attributable to horizontal diffusion of temperature and a contribution from frictional heating (see sections on horizontal diffusion and horizontal diffusion correction).

In addition to the prognostic equations, three diagnostic equations are required:

$$\Phi = \Phi_s + R\int_{p(\eta)}^{p(1)}{T_v} d\ln p ,$$ (3.14)

$$\dot\eta \frac{\partial p}{\partial\eta} = -\frac{\partial p}{\partial t} - \int^\eta_{\eta_t} {\boldsymbo... ...bla}\cdot}\left(\frac{\partial p}{\partial\eta}{\boldsymbol {V}}\right) d\eta ,$$ (3.15)

$$\omega = {\boldsymbol {V} \cdot\nabla}p - \int^\eta_{\eta_t} {\boldsymbo... ...\cdot} \left( \frac{\partial p}{\partial\eta} {\boldsymbol {V}} \right) d\eta .$$ (3.16)

Note that the bounds on the vertical integrals are specified as values of $\eta$ (e.g. $\eta_t$, 1) or as functions of $p$ (e.g. $p$ (1), which is the pressure at $\eta = 1$).

3.1.2 Conversion to final form

Equations (3.1)-(3.16) are the complete set which must be solved by a GCM. However, in order to solve them, the function $\eta(p,\pi)$ must be specified. In advance of actually specifying $\eta(p,\pi)$, the equations will be cast in a more convenient form. Most of the changes to the equations involve simple applications of the chain rule for derivatives, in order to obtain terms that will be easy to evaluate using the predicted variables in the model. For example, terms involving horizontal derivatives of $p$ must be converted to terms involving only $\partial p/\partial\pi$ and horizontal derivatives of $\pi$. The former can be evaluated once the function $\eta(p,\pi)$ is specified.

The vertical advection terms in (3.5), (3.6), (3.8), and (3.9) may be rewritten as:

$$\dot\eta \frac{\partial \psi}{\partial\eta} = \dot\eta \frac{\partial p}{\partial\eta} \frac{\partial\psi}{\partial p} ,$$ (3.17)

since $\dot\eta {\partial p/\partial\eta}$ is given by (3.15). Similarly, the first term on the right-hand side of (3.15) can be expanded as

$$\frac{\partial p}{\partial t} = \frac{\partial p}{\partial\pi} \frac{\partial\pi}{\partial t} ,$$ (3.18)

and (3.7) invoked to specify $\partial\pi/\partial t$.

The integrals which appear in (3.7), (3.15), and (3.16) can be written more conveniently by expanding the kernel as

$${\boldsymbol{\nabla}\cdot} \left( \frac{\partial p}{\partial\eta}... ...+ \frac{\partial p}{\partial\eta} {\boldsymbol{\nabla}\cdot\boldsymbol {V}} .$$ (3.19)

The second term in (3.19) is easily treated in vertical integrals, since it reduces to an integral in pressure. The first term is expanded to:

$${\boldsymbol {V}\cdot\nabla} \left(\frac{\partial p}{\partial\eta}\right) = {\boldsymbol {V}\cdot}\frac{\partial}{\partial\eta}\left(\nabla p\right)$$

$$= {\boldsymbol {V}\cdot}\frac{\partial}{\partial\eta} \left(\frac{\partial p}{\partial\pi}\nabla\pi\right)$$

$$= {\boldsymbol {V}\cdot}\frac{\partial}{\partial\eta} \left(\frac... ...partial p}{\partial\pi} \nabla\left(\frac{\partial\pi}{\partial\eta}\right) .$$ (3.20)

The second term in (3.20) vanishes because $\partial\pi/\partial\eta=0$, while the first term is easily treated once $\eta(p,\pi)$ is specified. Substituting (3.20) into (3.19), one obtains:

$${\boldsymbol{\nabla}\cdot} \left( \frac{\partial p}{\partial\eta}... ...t\nabla}\pi + \frac{\partial p}{\partial\eta} {\boldsymbol{\nabla}\cdot V} .$$ (3.21)

Using (3.21) as the kernel of the integral in (3.7), (3.15), and (3.16), one obtains integrals of the form

$$\int {\boldsymbol{\nabla}\cdot} \left( \frac{\partial p}{\partial\eta} {\boldsymbol {V}}\right) d\eta = \int \left[ \frac{\partial}{\partial\eta} \left(\frac{\partial ... ...pi + \frac{\partial p}{\partial\eta} {\boldsymbol{\nabla}\cdot V} \right] d\eta$$

$$= \int {\boldsymbol {V}\cdot\nabla}\pi d\left(\frac{\partial p}{\partial\pi}\right) + \int \delta dp.$$ (3.22)

The original primitive equations (3.3)-(3.7), together with (3.8), (3.9), and (3.14)-(3.16) can now be rewritten with the aid of (3.17), (3.18), and (3.22).

$$\frac{\partial\zeta}{\partial t} = {\boldsymbol {k}}\cdot\nabla\times ({\boldsymbol {n}}/\cos\phi) + F_{\zeta_H} ,$$ (3.23)

$$\frac{\partial\delta}{\partial t} = {\boldsymbol{\nabla}\cdot (\boldsymbol {n}/\cos\phi)} -\nabla^2\left(E+\Phi \right) + F_{\delta_H} ,$$ (3.24)

$$\frac{\partial T}{\partial t} = \frac{-1}{a\cos^2\phi} \left[ \frac{\partial}{\partial\lambda} (U... ...rtial\eta} \frac{\partial T}{\partial p} + \frac{R}{c_p^*}{T_v}\frac{\omega}{p}$$

$$\phantom{=} + Q + F_{T_H} + F_{F_H}$$ (3.25)

$$\frac{\partial q}{\partial t} = \frac{-1}{a\cos^2\phi} \left[ \frac{\partial}{\partial\lambda} (U... ... - \dot\eta \frac{\partial p}{\partial\eta} \frac{\partial q}{\partial p} + S ,$$ (3.26)

$$\frac{\partial \pi}{\partial t} = -\int_{(\eta_t)}^{(1)} {\boldsymbol {V}\cdot\nabla}\pi d\left(\frac{\partial p}{\partial\pi}\right) -\int_{p(\eta_t)}^{p(1)} \delta dp ,$$ (3.27)

$$n_U = + (\zeta+f)V - \dot\eta \frac{\partial p}{\partial\eta} \frac{\pa... ...}{p} \frac{\partial p}{\partial\pi} \frac{\partial\pi}{\partial\lambda} + F_U ,$$ (3.28)

$$n_V = - (\zeta+f)U - \dot\eta \frac{\partial p}{\partial\eta} \frac{\pa... ...c{1}{p} \frac{\partial p}{\partial\pi} \frac{\partial\pi}{\partial\phi} + F_V ,$$ (3.29)

$$\Phi = \Phi_s + R\int_{p(\eta)}^{p(1)}{T_v} d\ln p ,$$ (3.30)

$$\dot\eta \frac{\partial p}{\partial\eta} = \frac{\partial p}{\partial\pi} \left[ \int_{(\eta_t)}^{(1)} {\bol... ...frac{\partial p}{\partial\pi}\right) +\int_{p(\eta_t)}^{p(1)} \delta dp \right]$$ (3.31)

$$\phantom{=} -\int_{(\eta_t)}^{(\eta)} {\boldsymbol {V}}\cdot\nabla\pi d\left( \frac{\partial p}{\partial\pi}\right) -\int_{p(\eta_t)}^{p(\eta)} \delta dp ,$$

$$\omega = \frac{\partial p}{\partial\pi} {\boldsymbol {V} \cdot\nabla}\pi -... ...(\frac{\partial p}{\partial\pi}\right) - \int_{p(\eta_t)}^{p(\eta)} \delta dp .$$ (3.32)

Once $\eta(p,\pi)$ is specified, then $\partial p/\partial\pi$ can be determined and (3.23)-(3.32) can be solved in a GCM.

In the actual definition of the hybrid coordinate, it is not necessary to specify $\eta(p,\pi)$ explicitly, since (3.23)-(3.32) only requires that $p$ and $\partial p/\partial\pi$ be determined. It is sufficient to specify $p(\eta,\pi)$ and to let $\eta$ be defined implicitly. This will be done in section 3.1.7. In the case that $p(\eta,\pi)=\sigma\pi$ and $\eta_t=0$, (3.23)-(3.32) can be reduced to the set of equations solved by CCM1.

3.1.3 Continuous equations using $\partial\ln(\pi)/\partial t$

In practice, the solutions generated by solving the above equations are excessively noisy. This problem appears to arise from aliasing problems in the hydrostatic equation (3.30). The $\ln p$ integral introduces a high order nonlinearity which enters directly into the divergence equation (3.24). Large gravity waves are generated in the vicinity of steep orography, such as in the Pacific Ocean west of the Andes.

The noise problem is solved by converting the equations given above, which use $\pi$ as a prognostic variable, to equations using $\Pi=\ln(\pi)$. This results in the hydrostatic equation becoming only quadratically nonlinear except for moisture contributions to virtual temperature. Since the spectral transform method will be used to solve the equations, gradients will be obtained during the transform from wave to grid space. Outside of the prognostic equation for $\Pi$, all terms involving $\nabla\pi$ will then appear as $\pi\nabla\Pi$.

Equations (3.23)-(3.32) become:

$$\frac{\partial\zeta}{\partial t} = {\boldsymbol {k}\cdot\nabla\times (\boldsymbol {n}/\cos\phi)} + F_{\zeta_H} ,$$ (3.33)

$$\frac{\partial\delta}{\partial t} = {\boldsymbol{\nabla}\cdot (\boldsymbol {n}/\cos\phi)} -\nabla^2\left(E+\Phi \right) + F_{\delta_H} ,$$ (3.34)

$$\frac{\partial T}{\partial t} = \frac{-1}{a\cos^2\phi} \left[ \frac{\partial}{\partial\lambda} (U... ...rtial\eta} \frac{\partial T}{\partial p} + \frac{R}{c_p^*}{T_v}\frac{\omega}{p}$$ (3.35)

$$\phantom{=} + Q + F_{T_H} + F_{F_H} ,$$

$$\frac{\partial q}{\partial t} = \frac{-1}{a\cos^2\phi} \left[ \frac{\partial}{\partial\lambda} (U... ... - \dot\eta \frac{\partial p}{\partial\eta} \frac{\partial q}{\partial p} + S ,$$ (3.36)

$$\frac{\partial \Pi}{\partial t} = -\int_{(\eta_t)}^{(1)} {\boldsymbol {V}\cdot\nabla}\Pi d\left(\fr... ...artial p}{\partial\pi}\right) -\frac{1}{\pi}\int_{p(\eta_t)}^{p(1)} \delta dp ,$$ (3.37)

$$n_U = + (\zeta+f)V - \dot\eta \frac{\partial p}{\partial\eta} \frac{\pa... ...}{p} \frac{\partial p}{\partial\pi} \frac{\partial\Pi}{\partial\lambda} + F_U ,$$ (3.38)

$$n_V = - (\zeta+f)U - \dot\eta \frac{\partial p}{\partial\eta} \frac{\pa... ...\pi}{p} \frac{\partial p}{\partial\pi} \frac{\partial\Pi}{\partial\phi} + F_V ,$$ (3.39)

$$\Phi = \Phi_s + R\int_{p(\eta)}^{p(1)}{T_v} d\ln p ,$$ (3.40)

$$\dot\eta \frac{\partial p}{\partial\eta} = \frac{\partial p}{\partial\pi} \left[ \int_{(\eta_t)}^{(1)} \pi {... ...frac{\partial p}{\partial\pi}\right) +\int_{p(\eta_t)}^{p(1)} \delta dp \right]$$ (3.41)

$$\phantom{=} -\int_{(\eta_t)}^{(\eta)} \pi {\boldsymbol {V}}\cdot\nabla\Pi d\left(\frac{\partial p}{\partial\pi}\right) -\int_{p(\eta_t)}^{p(\eta)} \delta dp ,$$

$$\omega = \frac{\partial p}{\partial\pi} \pi {\boldsymbol {V} \cdot\nabla}\... ...(\frac{\partial p}{\partial\pi}\right) - \int_{p(\eta_t)}^{p(\eta)} \delta dp .$$ (3.42)

The above equations reduce to the standard $\sigma$ equations used in CCM1 if $\eta=\sigma$ and $\eta_t=0$. (Note that in this case $\partial p / \partial\pi = p/\pi = \sigma$.)

3.1.4 Semi-implicit formulation

The model described by (3.33)-(3.42), without the horizontal diffusion terms, together with boundary conditions (3.1) and (3.2), is integrated in time using the semi-implicit leapfrog scheme described below. The semi-implicit form of the time differencing will be applied to (3.34) and (3.36) without the horizontal diffusion sources, and to (3.37). In order to derive the semi-implicit form, one must linearize these equations about a reference state. Isolating the terms that will have their linear parts treated implicitly, the prognostic equations (3.33), (3.34), and (3.37) may be rewritten as:

$$\frac{\partial\delta}{\partial t} = - {R{T_v}} \nabla^2 \ln p -\nabla^2\Phi + X_1 ,$$ (3.43)

$$\frac{\partial T}{\partial t} = + \frac{R}{c_p^*}{T_v}\frac{\omega}{p} - \dot\eta \frac{\partial p}{\partial\eta} \frac{\partial T}{\partial p} + Y_1 ,$$ (3.44)

$$\frac{\partial \Pi}{\partial t} = - \frac{1}{\pi}\int_{p(\eta_t)}^{p(1)} \delta dp + Z_1 ,$$ (3.45)

where $X_1, Y_1, Z_1$ are the remaining nonlinear terms not explicitly written in (3.43)-(3.45). The terms involving $\Phi$ and $\omega$ may be expanded into vertical integrals using (3.40) and (3.42), while the $\nabla^2 \ln p$ term can be converted to $\nabla^2\Pi$, giving:

$$\frac{\partial\delta}{\partial t} = -{RT} \frac{\pi}{p}\frac{\partial p}{\partial \pi} \nabla^2 \Pi -R\nabla^2\int_{p(\eta)}^{p(1)}T d\ln p + X_2 ,$$ (3.46)

$$\frac{\partial T}{\partial t} = - \frac{R}{c_p}\frac{T}{p} \int_{p(\eta_t)}^{p(\eta)}\delta dp ... ...t_{p(\eta_t)}^{p(\eta)} \delta dp \right] \frac{\partial T}{\partial p} + Y_2 ,$$ (3.47)

$$\frac{\partial \Pi}{\partial t} = - \frac{1}{pi}\int_{p(\eta_t)}^{p(1)} \delta dp + Z_2 .$$ (3.48)

Once again, only terms that will be linearized have been explicitly represented in (3.46)-(3.48), and the remaining terms are included in $X_2$, $Y_2$, and $Z_2$. Anticipating the linearization, $T_v$ and $c_p^*$ have been replaced by $T$ and $c_p$ in (3.46) and (3.47). Furthermore, the virtual temperature corrections are included with the other nonlinear terms.

In order to linearize (3.46)-(3.48), one specifies a reference state for temperature and pressure, then expands the equations about the reference state:

$$T = T^r + T^\prime ,$$ (3.49)

$$\pi = \pi^r + \pi^\prime ,$$ (3.50)

$$p = p^r(\eta,\pi^r) + p^\prime .$$ (3.51)

In the special case that $p(\eta,\pi)=\sigma\pi$, (3.46)-(3.48) can be converted into equations involving only $\Pi=\ln\pi$ instead of $p$, and (3.50) and (3.51) are not required. This is a major difference between the hybrid coordinate scheme being developed here and the $\sigma$ coordinate scheme in CCM1.

Expanding (3.46)-(3.48) about the reference state (3.49)-(3.51) and retaining only the linear terms explicitly, one obtains:

$$\frac{\partial\delta}{\partial t} = -R\nabla^2 \left[ T^{r} \frac{\pi^r}{p^r} \left(\frac{\partial ... ... + \int_{p^\prime(\eta)}^{p^\prime(1)}\frac{T^r}{p^r} dp^\prime \right] + X_3 ,$$ (3.52)

$$\frac{\partial T}{\partial t} = - \frac{R}{c_p}\frac{T^r}{p^r} \int_{p^r(\eta_t)}^{p^r(\eta)}\d... ...ta_t)}^{p^r(\eta)} \delta dp^r \right] \frac{\partial T^r}{\partial p^r} + Y_3,$$ (3.53)

$$\frac{\partial \Pi}{\partial t} = - \frac{1}{\pi^r}\int_{p^r(\eta_t)}^{p^r(1)} \delta dp^r + Z_3 .$$ (3.54)

The semi-implicit time differencing scheme treats the linear terms in (3.52)-(3.54) by averaging in time. The last integral in (3.52) is reduced to purely linear form by the relation

$$dp^\prime = \pi^\prime d \left(\frac{\partial p}{\partial\pi}\right)^r + x .$$ (3.55)

In the hybrid coordinate described below, $p$ is a linear function of $\pi$, so $x$ above is zero.

We will assume that centered differences are to be used for the nonlinear terms, and the linear terms are to be treated implicitly by averaging the previous and next time steps. Finite differences are used in the vertical, and are described in the following sections. At this stage only some very general properties of the finite difference representation must be specified. A layering structure is assumed in which field values are predicted on $K$ layer midpoints denoted by an integer index, $\eta_k$ (see Figure 3.1). The interface between $\eta_k$ and $\eta_{k+1}$ is denoted by a half-integer index, $\eta_{k+1/2}$. The model top is at $\eta_{1/2}=\eta_t$, and the Earth's surface is at $\eta_{K+1/2}=1$. It is further assumed that vertical integrals may be written as a matrix (of order $K$) times a column vector representing the values of a field at the $\eta_k$ grid points in the vertical. The column vectors representing a vertical column of grid points will be denoted by underbars, the matrices will be denoted by bold-faced capital letters, and superscript $T$ will denote the vector transpose.

Figure 3.1: Vertical level structure of CAM 3.0

The finite difference forms of (3.52)-(3.54) may then be written down as:

$$\underline{\delta}^{n+1} = \underline{\delta}^{n-1} + 2\Delta t \underline{X}^n$$

$$\phantom{=} - 2\Delta t R \underline{b}^r \nabla^2 \left( \frac{\Pi^{n-1} + \Pi^{n+1}}{2} - \Pi^{n} \right)$$

$$\phantom{=} - 2\Delta t R{\boldsymbol {H}}^r \nabla^2 \left( \frac{(\underlin... ...)^{n-1} + (\underline{T}^\prime)^{n+1}}{2} - (\underline{T}^\prime)^{n} \right)$$

$$\phantom{=} - 2\Delta t R \underline{h}^r \nabla^2 \left( \frac{\Pi^{n-1} + \Pi^{n+1}}{2} - \Pi^{n} \right) ,$$ (3.56)

$$\underline{T}^{n+1} = \underline{T}^{n-1} + 2 \Delta t \underline{Y}^n - 2\Delta t {\bo... ...e{\delta}^{n-1} + \underline{\delta}^{n+1}}{2} - \underline{\delta}^n \right) ,$$ (3.57)

$$\Pi^{n+1} = \Pi^{n-1} + 2\Delta t Z^n - 2\Delta t \left( \frac{\underline{\de... ...}}{2} - \underline{\delta}^n \right)^T \frac{1}{\Pi^r} \underline{\Delta p}^r ,$$ (3.58)

where $()^n$ denotes a time varying value at time step $n$. The quantities $\underline{X}^n, \underline{Y}^n,$ and $Z^n$ are defined so as to complete the right-hand sides of (3.43)-(3.45). The components of $\underline{\Delta p}^r$ are given by $\Delta p^r_k = p^r_{k + \frac{1}{2}} - p^r_{k - \frac{1}{2}}$. This definition of the vertical difference operator $\Delta$ will be used in subsequent equations. The reference matrices ${\boldsymbol {H}}^r$ and ${\boldsymbol {D}}^r$, and the reference column vectors $$\underline{b}^r$$ and $$\underline{h}^r$$, depend on the precise specification of the vertical coordinate and will be defined later.

3.1.5 Energy conservation

We shall impose a requirement on the vertical finite differences of the model that they conserve the global integral of total energy in the absence of sources and sinks. We need to derive equations for kinetic and internal energy in order to impose this constraint. The momentum equations (more painfully, the vorticity and divergence equations) without the $F_U, F_V, F_{\zeta_H}$ and $F_{\delta_H}$ contributions, can be combined with the continuity equation

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta... ...ial}{\partial \eta} \left( \frac{\partial p}{\partial\eta} \dot\eta \right) = 0$$ (3.59)

to give an equation for the rate of change of kinetic energy:

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta} E \right) = -\nabla\cdot \left( \frac{\partial p}{\partial\eta} E {\boldsymbo... ...rtial}{\partial \eta} \left( \frac{\partial p}{\partial\eta} E \dot\eta \right)$$

$$\phantom{=} - \frac{R{T_v}}{p} \frac{\partial p}{\partial\eta} {\boldsymbol {... ...bla p - \frac{\partial p}{\partial\eta} {\boldsymbol {V}}\cdot\nabla\Phi - .$$ (3.60)

The first two terms on the right-hand side of (3.60) are transport terms. The horizontal integral of the first (horizontal) transport term should be zero, and it is relatively straightforward to construct horizontal finite difference schemes that ensure this. For spectral models, the integral of the horizontal transport term will not vanish in general, but we shall ignore this problem.

The vertical integral of the second (vertical) transport term on the right-hand side of (3.60) should vanish. Since this term is obtained from the vertical advection terms for momentum, which will be finite differenced, we can construct a finite difference operator that will ensure that the vertical integral vanishes.

The vertical advection terms are the product of a vertical velocity ( $\dot\eta \partial p/\partial\eta$) and the vertical derivative of a field ( $\partial\psi/\partial p$). The vertical velocity is defined in terms of vertical integrals of fields (3.42), which are naturally taken to interfaces. The vertical derivatives are also naturally taken to interfaces, so the product is formed there, and then adjacent interface values of the products are averaged to give a midpoint value. It is the definition of the average that must be correct in order to conserve kinetic energy under vertical advection in (3.60). The derivation will be omitted here, the resulting vertical advection terms are of the form:

$$\left( \dot\eta \frac{\partial p}{\partial\eta} \frac{\partial \psi}{\partial p} \right)_{k} = \frac{1}{2\Delta p_k} \left[ \left( \dot\eta \frac{\partial p}{\p... ...l p}{\partial\eta} \right)_{k-1/2} \left( \psi_k - \psi_{k-1} \right) \right] ,$$ (3.61)

$$\Delta p_k = p_{k+1/2} - p_{k-1/2} .$$ (3.62)

The choice of definitions for the vertical velocity at interfaces is not crucial to the energy conservation (although not completely arbitrary), and we shall defer its definition until later. The vertical advection of temperature is not required to use (3.61) in order to conserve mass or energy. Other constraints can be imposed that result in different forms for temperature advection, but we will simply use (3.61) in the system described below.

The last two terms in (3.60) contain the conversion between kinetic and internal (potential) energy and the form drag. Neglecting the transport terms, under assumption that global integrals will be taken, noting that $\nabla p/p = \frac{\pi}{p} \frac{\partial p}{\partial\pi} \nabla \Pi$, and substituting for the geopotential using (3.40), (3.60) can be written as:

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta} E \right) = - {R{T_v}} \frac{\partial p}{\partial\eta} {\boldsymbol {V}} \cdot \left( \frac{\pi}{p} \frac{\partial p}{\partial\pi} \nabla \Pi \right)$$ (3.63)

$$\phantom{=} - \frac{\partial p}{\partial\eta} {\boldsymbol {V}}\cdot\nabla\Ph... ...ta} {\boldsymbol {V}}\cdot\nabla \int_{p(\eta)}^{p(1)}R{T_v} d\ln p + \ldots$$

The second term on the right-hand side of (3.64) is a source (form drag) term that can be neglected as we are only interested in internal conservation properties. The last term on the right-hand side of (3.64) can be rewritten as

$$\frac{\partial p}{\partial\eta} {\boldsymbol {V}}\cdot\nabla \int... ...\partial\eta} {\boldsymbol {V}} \right) \int_{p(\eta)}^{p(1)}R{T_v} d\ln p .$$ (3.64)

The global integral of the first term on the right-hand side of (3.64) is obviously zero, so that (3.64) can now be written as:

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta... ...partial\eta} {\boldsymbol {V}} \right) \int_{p(\eta)}^{p(1)}R{T_v} d\ln p + ...$$ (3.65)

We now turn to the internal energy equation, obtained by combining the thermodynamic equation (3.36), without the $Q$, $F_{T_H}$, and $F_{F_H}$ terms, and the continuity equation (3.59):

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta... ...ot\eta \right) + {R{T_v}} \frac{\partial p}{\partial\eta} \frac{\omega}{p} .$$ (3.66)

As in (3.60), the first two terms on the right-hand side are advection terms that can be neglected under global integrals. Using (3.16), (3.66) can be written as:

$$\frac{\partial}{\partial t} \left( \frac{\partial p}{\partial\eta... ...ot \left( \frac{\partial p}{\partial\eta} {\boldsymbol {V}} \right) d\eta + ...$$ (3.67)

The rate of change of total energy due to internal processes is obtained by adding (3.65) and (3.67) and must vanish. The first terms on the right-hand side of (3.65) and (3.67) obviously cancel in the continuous form. When the equations are discretized in the vertical, the terms will still cancel, providing that the same definition is used for $(1/p \partial p/\partial\pi)_k$ in the nonlinear terms of the vorticity and divergence equations (3.38) and (3.39), and in the $\omega$ term of (3.36) and (3.42).

The second terms on the right-hand side of (3.65) and (3.67) must also cancel in the global mean. This cancellation is enforced locally in the horizontal on the column integrals of (3.65) and (3.67), so that we require:

$$\int^1_{\eta_t} \left\{ \nabla\cdot \left( \frac{\partial p}{\par... ...p}{\partial\eta^\prime} {\boldsymbol {V}} \right) d\eta^\prime \right\} d\eta .$$ (3.68)

The inner integral on the left-hand side of (3.68) is derived from the hydrostatic equation (3.40), which we shall approximate as

$$\Phi_k = \Phi_s + R\sum_{\ell=k}^K H_{k\ell}{T_v}_\ell ,$$

$$= \Phi_s + R\sum_{\ell=1}^K H_{k\ell}{T_v}_\ell ,$$ (3.69)

$$\underline{\Phi} = \Phi_s \underline{1} + R {\boldsymbol {H}} \underline {T_v} ,$$ (3.70)

where $H_{k\ell}=0$ for $\ell<k$. The quantity $\underline{1}$ is defined to be the unit vector. The inner integral on the right-hand side of (3.68) is derived from the vertical velocity equation (3.42), which we shall approximate as

$$\left( \frac{\omega}{p} \right)_k = \left( \frac{\pi}{p} \frac{\p... ...\Pi \right) \Delta \left( \frac{\partial p}{\partial\pi} \right)_\ell \right] ,$$ (3.71)

where $C_{k\ell}=0$ for $\ell>k$, and $C_{k\ell}$ is included as an approximation to $1/p_k$ for $\ell \leq k$ and the symbol $\Delta$ is similarly defined as in (3.62). $C_{k\ell}$ will be determined so that $\omega$ is consistent with the discrete continuity equation following Williamson and Olson [185]. Using (3.69) and (3.71), the finite difference analog of (3.68) is

$${\sum_{k=1}^K \left\{ \frac{1}{\Delta\eta_k} \left[ \delta_k \Del... ... \right)_k \right] R\sum_{\ell=1}^K H_{k\ell}{T_v}_\ell \right\} \Delta\eta_k }$$

$$= \sum_{k=1}^K \left\{ R {T_v}_k \frac{\Delta p_k}{\Delta\eta_k} ... ...ft( \frac{\partial p}{\partial\pi} \right)_\ell \right] \right\} \Delta\eta_k ,$$ (3.72)

where we have used the relation

$$\nabla\cdot {\boldsymbol {V}}(\partial p/\partial\eta )_k =[\delt... ...ot\nabla\Pi \right) \Delta \left(\partial p/\partial\pi\right)_k ]/\Delta\eta_k$$ (3.73)

(see 3.22). We can now combine the sums in (3.72) and simplify to give

$${\sum_{k=1}^K \sum_{\ell=1}^K \left\{ \left[ \delta_k \Delta p_k ... ... \frac{\partial p}{\partial\pi} \right)_k \right] H_{k\ell}{T_v}_\ell \right\}}$$

$$= \sum_{k=1}^K \sum_{\ell=1}^K \left\{ \left[ \delta_\ell \Delta ... ...ial p}{\partial\pi} \right)_\ell \right] \Delta p_k C_{k\ell}{T_v}_k \right\} .$$ (3.74)

Interchanging the indexes on the left-hand side of (3.74) will obviously result in identical expressions if we require that

$$H_{k\ell} = C_{\ell k} \Delta p_\ell .$$ (3.75)

Given the definitions of vertical integrals in (3.70) and (3.71) and of vertical advection in (3.61) and (3.62) the model will conserve energy as long as we require that $\boldsymbol {C}$ and $\boldsymbol {H}$ satisfy (3.75). We are, of course, still neglecting lack of conservation due to the truncation of the horizontal spherical harmonic expansions.

3.1.6 Horizontal diffusion

CAM 3.0 contains a horizontal diffusion term for $T, \zeta$, and $\delta$ to prevent spectral blocking and to provide reasonable kinetic energy spectra. The horizontal diffusion operator in CAM 3.0 is also used to ensure that the CFL condition is not violated in the upper layers of the model. The horizontal diffusion is a linear $\nabla^2$ form on $\eta$ surfaces in the top three levels of the model and a linear $\nabla^4$ form with a partial correction to pressure surfaces for temperature elsewhere. The $\nabla^2$ diffusion near the model top is used as a simple sponge to absorb vertically propagating planetary wave energy and also to control the strength of the stratospheric winter jets. The $\nabla^2$ diffusion coefficient has a vertical variation which has been tuned to give reasonable Northern and Southern Hemisphere polar night jets.

In the top three model levels, the $\nabla^2$ form of the horizontal diffusion is given by

$$F_{\zeta_H} = K^{(2)} \left[ \nabla^2 \left(\zeta + f \right) + 2\left(\zeta + f \right)/a^2 \right] ,$$ (3.76)

$$F_{\delta_H} = K^{(2)} \left[ \nabla^2 \delta + 2(\delta/a^2)\right] ,$$ (3.77)

$$F_{T_H} = K^{(2)} \nabla^2T .$$ (3.78)

Since these terms are linear, they are easily calculated in spectral space. The undifferentiated correction term is added to the vorticity and divergence diffusion operators to prevent damping of uniform $(n=1)$ rotations [24,133]. The $\nabla^2$ form of the horizontal diffusion is applied only to pressure surfaces in the standard model configuration.

The horizontal diffusion operator is better applied to pressure surfaces than to terrain-following surfaces (applying the operator on isentropic surfaces would be still better). Although the governing system of equations derived above is designed to reduce to pressure surfaces above some level, problems can still occur from diffusion along the lower surfaces. Partial correction to pressure surfaces of harmonic horizontal diffusion ( $\partial\xi/\partial t = K\nabla^2\xi$) can be included using the relations:

$$\nabla_p\xi = \nabla_\eta\xi - p \frac{\partial\xi}{\partial p} \nabla_\eta \ln p$$

$$\nabla^2_p\xi = \nabla^2_\eta\xi - p \frac{\partial\xi}{\partial p} \nabla^2_\e... ...ht)\cdot\nabla_\eta p + \frac{\partial^2\xi}{\partial^2 p} \nabla^2_\eta p .$$ (3.79)

Retaining only the first two terms above gives a correction to the $\eta$ surface diffusion which involves only a vertical derivative and the Laplacian of log surface pressure,

$$\nabla^2_p\xi = \nabla^2_\eta\xi - \pi \frac{\partial\xi}{\partial p} \frac{\partial p}{\partial\pi} \nabla^2 \Pi + \ldots$$ (3.80)

Similarly, biharmonic diffusion can be partially corrected to pressure surfaces as:

$$\nabla^4_p\xi = \nabla^4_\eta\xi - \pi \frac{\partial\xi}{\partial p} \frac{\partial p}{\partial\pi} \nabla^4 \Pi + \ldots$$ (3.81)

The bi-harmonic $\nabla^4$ form of the diffusion operator is applied at all other levels (generally throughout the troposphere) as

$$F_{\zeta_H} = -K^{(4)} \left[\nabla^4 \left(\zeta + f \right) - \left(\zeta + f \right) \left( 2/a^2\right)^2 \right] ,$$ (3.82)

$$F_{\delta_H} = -K^{(4)} \left[\nabla^4 \delta - \delta (2/a^2)^2 \right],$$ (3.83)

$$F_{T_H} = -K^{(4)} \left[ \nabla^4 T - \pi \frac{\partial T}{\partial p} \frac{\partial p}{\partial \pi} \nabla^4 \Pi \right] .$$ (3.84)

The second term in $F_{T_H}$ consists of the leading term in the transformation of the $\nabla^4$ operator to pressure surfaces. It is included to offset partially a spurious diffusion of $T$ over mountains. As with the $\nabla^2$ form, the $\nabla^4$ operator can be conveniently calculated in spectral space. The correction term is then completed after transformation of $T$ and $\nabla^4 \Pi$ back to grid-point space. As with the $\nabla^2$ form, an undifferentiated term is added to the vorticity and divergence diffusion operators to prevent damping of uniform rotations.

3.1.7 Finite difference equations

The governing equations are solved using the spectral method in the horizontal, so that only the vertical and time differences are presented here. The dynamics includes horizontal diffusion of $T, (\zeta + f)$, and $\delta$. Only $T$ has the leading term correction to pressure surfaces. Thus, equations that include the terms in this time split sub-step are of the form

$$\frac{\partial \psi}{\partial t} = {\rm Dyn} \left( \psi \right) - (-1)^i K^{(2i)} \nabla^{2i}_\eta \psi ,$$ (3.85)

for $(\zeta + f)$ and $\delta$, and

$$\frac{\partial T}{\partial t} = {\rm Dyn} \left( T \right) - (-1)... ...T} {\partial p} \frac{\partial p}{\partial \pi} \nabla^{2i} \Pi \right\} ,$$ (3.86)

where $i = 1$ in the top few model levels and $i = 2$ elsewhere (generally within the troposphere). These equations are further subdivided into time split components:

$$\psi^{n+1} = \psi^{n-1} + 2\Delta t {\rm Dyn} \left( \psi^{n+1}, \psi^n, \psi^{n-1} \right) ,$$ (3.87)

$$\psi^* = \psi^{n+1} - 2\Delta t (-1)^i K^{(2i)} \nabla^{2i}_\eta \left( {\psi}^{*n+1} \right) ,$$ (3.88)

$$\hat \psi^{n+1} = \psi^* ,$$ (3.89)

for $\left( \zeta + f \right)$ and $\delta$, and

$$T^{n+1} = T^{n-1} + 2\Delta t {\rm Dyn} \left( T^{n+1}, T^n, T^{n-1} \right)$$ (3.90)

$$T^* = T^{n+1} - 2\Delta t \left( -1 \right)^i K^{(2i)} \nabla^{2i}\eta \left( T^* \right) ,$$ (3.91)

$$\hat T^{n+1} = T^* + 2\Delta t \left( -1 \right)^i K^{(2i)} \pi \frac{\pa... ...l T^*}{\partial p} \frac{\partial p}{\partial \pi} \nabla^{2i} \Pi ,$$ (3.92)

for $T$, where in the standard model $i$ only takes the value 2 in (3.92). The first step from $\left( \hskip 5pt \right)^{n-1}$ to $\left( \hskip 5pt \right)^{n+1}$ includes the transformation to spectral coefficients. The second step from $\left( \hskip 5pt \right)^{n+1}$ to $\left( \hat{\hskip 5pt} \right)^{n+1}$ for $\delta$ and $\zeta$, or $\left( \hskip 5pt \right)^{n+1}$ to $\left( \hskip 5pt \right)^*$ for $T$, is done on the spectral coefficients, and the final step from $\left( \hskip 5pt \right)^*$ to $\left( \hat{\hskip 5pt} \right)^{n+1}$ for $T$ is done after the inverse transform to the grid point representation.

The following finite-difference description details only the forecast given by (3.87) and (3.90). The finite-difference form of the forecast equation for water vapor will be presented later in Section 3c. The general structure of the complete finite difference equations is determined by the semi-implicit time differencing and the energy conservation properties described above. In order to complete the specification of the finite differencing, we require a definition of the vertical coordinate. The actual specification of the generalized vertical coordinate takes advantage of the structure of the equations (3.33)-(3.42). The equations can be finite-differenced in the vertical and, in time, without having to know the value of $\eta$ anywhere. The quantities that must be known are $p$ and $\partial p/\partial\pi$ at the grid points. Therefore the coordinate is defined implicitly through the relation:

$$p(\eta,\pi) = A(\eta)p_0 + B(\eta)\pi ,$$ (3.93)

which gives

$$\frac{\partial p}{\partial\pi} = B(\eta) .$$ (3.94)

A set of levels $\eta_k$ may be specified by specifying $A_k$ and $B_k$, such that $\eta_k \equiv A_k + B_k$, and difference forms of (3.33)-(3.42) may be derived.

The finite difference forms of the Dyn operator (3.33)-(3.42), including semi-implicit time integration are:

$$\underline{\zeta}^{n+1} = \underline{\zeta}^{n-1} + 2\Delta t {\boldsymbol {k}\cdot\nabla\times}\left(\underline{\boldsymbol {n}}^n/\cos\phi\right) ,$$ (3.95)

$$\underline{\delta}^{n+1} = \underline{\delta}^{n-1} + 2\Delta t \left[ {\nabla\cdot \left(\u... ..._s \underline{1} + R{\boldsymbol {H}}^n (\underline{T_v}^{'})^n \right) \right]$$

$$\phantom{=} - 2\Delta t R{\boldsymbol {H}}^r \nabla^2 \left( \frac{(\underlin... ...)^{n-1} + (\underline{T}^\prime)^{n+1}}{2} - (\underline{T}^\prime)^{n} \right)$$

$$\phantom{=} - 2\Delta t R\left( \underline{b}^r + \underline{h}^r \right) \nabla^2 \left( \frac{\Pi^{n-1} + \Pi^{n+1}}{2} - \Pi^{n} \right) ,$$ (3.96)

$$(\underline{T}^{'})^{n+1} = (\underline{T}^{'})^{n-1} - 2 \Delta t \left[ \frac{1}{a\cos^2\ph... ...tial\phi} \left( \underline{VT}^\prime \right)^n - \underline{\Gamma}^n \right]$$ (3.97)

$$\phantom{=} - 2\Delta t {\boldsymbol {D}}^r \left( \frac{\underline{\delta}^{n-1} + \underline{\delta}^{n+1}}{2} - \underline{\delta}^n \right)$$

$$\Pi^{n+1} = \Pi^{n-1} - 2\Delta t \frac{1}{\pi^n} \left( \left(\underline{\de... ...ldsymbol {V}}^n \right)^T \cdot \nabla \Pi^n \pi^n \underline{\Delta B} \right)$$