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)$$

$$\phantom{=} - 2\Delta t \left( \frac{\underline{\delta}^{n-1} + \underline{\d... ...}}{2} - \underline{\delta}^n \right)^T \frac{1}{\pi^r} \underline{\Delta p}^r ,$$ (3.98)

$$\left( n_U \right)_k = \left( \zeta_k + f \right) V_k - R{T_v}_k \left( \frac{1}{p} \fra... ... p}{\partial\pi} \right)_k \pi \frac{1}{a} \frac{\partial \Pi}{\partial\lambda}$$

$$\phantom{=} - \frac{1}{2\Delta p_k} \left[ \left( \dot\eta \frac{\partial p}{... ...{\partial p}{\partial\eta} \right)_{k-1/2} \left( U_k - U_{k-1} \right) \right]$$

$$\phantom{=} + \left( F_U \right)_k ,$$ (3.99)

$$\left( n_V \right)_k = - \left( \zeta_k + f \right) U_k - R{T_v}_k \left( \frac{1}{p} \f... ...partial\pi} \right)_k \pi \frac{\cos \phi}{a} \frac{\partial \Pi}{\partial\phi}$$

$$\phantom{=} - \frac{1}{2\Delta p_k} \left[ \left( \dot\eta \frac{\partial p}{... ...{\partial p}{\partial\eta} \right)_{k-1/2} \left( V_k - V_{k-1} \right) \right]$$

$$\phantom{=} + \left( F_V \right)_k ,$$ (3.100)

$$\Gamma_k = T^{\prime}_k \delta_k + \frac{R{T_v}_k}{(c_p^*)_k} \left( \frac{\omega}{p} \right)_k - Q$$

$$\phantom{=} - \frac{1}{2\Delta p_k} \left[ \left( \dot\eta \frac{\partial p}{... ...partial p}{\partial\eta} \right)_{k-1/2} \left( T_k - T_{k-1} \right) \right] ,$$ (3.101)

$$E_k = \left(u_k \right)^2 + \left(v_k \right)^2 ,$$ (3.102)

$$\frac{R {T_v}_{k}}{(c^*_p)_k} = \frac{R}{c_p} \left( \frac{T^r_k + {T_v}_k^\prime} {1 + \left( \frac{c_{p_v}}{c_p}- 1 \right) q_k} \right) ,$$ (3.103)

$$\left( \dot\eta \frac{\partial p}{\partial\eta} \right)_{k+1/2} = B_{k+1/2} \sum^K_{\ell=1} \left[ \delta_\ell \Delta p_\ell + {\boldsymbol {V}}_\ell \cdot \pi \nabla \Pi \Delta B_\ell \right]$$

$$\phantom{=} - \sum^k_{\ell=1} \left[ \delta_\ell \Delta p_\ell + {\boldsymbol {V}}_\ell \cdot \pi \nabla \Pi \Delta B_\ell \right] ,$$ (3.104)

$$\left( \frac{\omega}{p} \right)_k = \left( \frac{1}{p} \frac{\partial p}{\partial\pi} \right)_k {\bol... ...ta p_\ell + {\boldsymbol {V}}_\ell \cdot \pi \nabla \Pi \Delta B_\ell \right] ,$$ (3.105)

$$C_{k\ell} = \left\{ \begin{array}{ll} \frac{1}{p_k}, & \ell < k \ [6pt] \frac{1}{2 p_k}, & \ell = k , \end{array} \right.$$ (3.106)

$$H_{k\ell} = C_{\ell k}\Delta p_\ell,$$ (3.107)

$$D^r_{k\ell} = {\Delta p^r_\ell} \frac{R}{c_p } T^r_k C^r_{\ell k} + \frac{\Delt... ...r_k} \left( T^r_k - T^r_{k-1} \right) \left(\epsilon_{k\ell+1}-B_{k-1/2}\right)$$

$$\phantom{=} + \frac{\Delta p^r_\ell}{2\Delta p^r_k} \left( T^r_{k+1} - T^r_k \right) \left(\epsilon_{k\ell}-B_{k+1/2}\right) ,$$ (3.108)

$$\frac{ \epsilon_{k\ell}}{R} = \left\{\begin{array}{ll} 1, & \ell \leq k \\ 0, & \ell > k, \end{array} \right.$$ (3.109)

where notation such as $\left( \underline{UT}^\prime \right)^n$ denotes a column vector with components $\left( U_k T_k^\prime \right)^n$. In order to complete the system, it remains to specify the reference vector $\underline{h}^r$, together with the term $(1/p \partial p/\partial\pi)$, which results from the pressure gradient terms and also appears in the semi-implicit reference vector $\underline{b}^r$:

$$\left( \frac{1}{p} \frac{\partial p}{\partial\pi} \right)_k = \left( \frac{1}{p} \right)_k \left( \frac{\partial p}{\partial\pi} \right)_k = \frac{B_k}{p_k} ,$$ (3.110)

$$\underline{b}^r = \underline{T}^r ,$$ (3.111)

$$\underline{h}^r = 0 .$$ (3.112)

The matrices ${\boldsymbol {C}}^n$ and ${\boldsymbol {H}}^n$ (i.e. with components $C_{k\ell}$ and $H_{k \ell}$) must be evaluated at each time step and each point in the horizontal. It is more efficient computationally to substitute the definitions of these matrices into (3.96) and (3.105) at the cost of some loss of generality in the code. The finite difference equations have been written in the form (3.95)-(3.112) because this form is quite general. For example, the equations solved by Simmons and Strüfing [158] at ECMWF can be obtained by changing only the vectors and hydrostatic matrix defined by (3.109)-(3.112).

3.1.8 Time filter

The time step is completed by applying a recursive time filter originally designed by [149] and later studied by [6].

$$\overline \psi^n = \psi^n + \alpha \left( \overline{\psi}^{n-1} - 2 \psi^n + \psi^{n+1} \right)$$ (3.113)

3.1.9 Spectral transform

The spectral transform method is used in the horizontal exactly as in CCM1. As shown earlier, the vertical and temporal aspects of the model are represented by finite-difference approximations. The horizontal aspects are treated by the spectral-transform method, which is described in this section. Thus, at certain points in the integration, the prognostic variables $\left(\zeta + f \right), \delta, T,$ and $\Pi$ are represented in terms of coefficients of a truncated series of spherical harmonic functions, while at other points they are given by grid-point values on a corresponding Gaussian grid. In general, physical parameterizations and nonlinear operations are carried out in grid-point space. Horizontal derivatives and linear operations are performed in spectral space. Externally, the model appears to the user to be a grid-point model, as far as data required and produced by it. Similarly, since all nonlinear parameterizations are developed and carried out in grid-point space, the model also appears as a grid-point model for the incorporation of physical parameterizations, and the user need not be too concerned with the spectral aspects. For users interested in diagnosing the balance of terms in the evolution equations, however, the details are important and care must be taken to understand which terms have been spectrally truncated and which have not. The algebra involved in the spectral transformations has been presented in several publications [24,46,118]. In this report, we present only the details relevant to the model code; for more details and general philosophy, the reader is referred to these earlier papers.

3.1.10 Spectral algorithm overview

The horizontal representation of an arbitrary variable $\psi$ consists of a truncated series of spherical harmonic functions,

$$\psi(\lambda, \mu) = \sum^M_{m=-M} \; \; \sum^{{\cal N} \left( m \right)}_{n=\vert m\vert} \psi^m_n P^m_n (\mu) e^{im \lambda} ,$$ (3.114)

where $\mu = \sin \phi, M$ is the highest Fourier wavenumber included in the east-west representation, and ${\cal N}\left( m \right)$ is the highest degree of the associated Legendre polynomials for longitudinal wavenumber $m$. The properties of the spherical harmonic functions used in the representation can be found in the review by Machenhauer [118]. The model is coded for a general pentagonal truncation, illustrated in Figure 3.2, defined by three parameters: $M, K$, and $N$, where $M$ is defined above, $K$ is the highest degree of the associated Legendre polynomials, and $N$ is the highest degree of the Legendre polynomials for $m= 0$. The common truncations are subsets of this pentagonal case:

$$\mathrm{Triangular:} \; M = N = K ,$$

$$\mathrm{Rhomboidal:} \; K = N + M ,$$ (3.115)

$$\mathrm{Trapezoidal:} \; N = K > M .$$

The quantity ${\cal N}\left( m \right)$ in (3.114) represents an arbitrary limit on the two-dimensional wavenumber $n$, and for the pentagonal truncation described above is simply given by
${\cal N} \left( m \right) = \min\left(N + \vert m \vert , K \right)$.

Figure 3.2: Pentagonal truncation parameters

The associated Legendre polynomials used in the model are normalized such that

$$\int^1_{-1} \left[P^m_n(\mu)\right]^2 d\mu = 1 .$$ (3.116)

With this normalization, the Coriolis parameter $f$ is

$$f = \frac{\Omega}{\sqrt {0.375}} P^o_1 ,$$ (3.117)

which is required for the absolute vorticity.

The coefficients of the spectral representation (3.114) are given by

$$\psi^m_n = \int^1_{-1} \frac{1}{2 \pi} \int^{2 \pi}_0 \psi (\lambda, \mu) e^{-im\lambda} d \lambda P^m_n (\mu) d \mu .$$ (3.118)

The inner integral represents a Fourier transform,

$$\psi^m (\mu) = \frac{1}{2 \pi} \int^{2 \pi}_0 \psi (\lambda, \mu) e^{-im\lambda} d \lambda ,$$ (3.119)

which is performed by a Fast Fourier Transform (FFT) subroutine. The outer integral is performed via Gaussian quadrature,

$$\psi^m_n = \sum^J_{j=1} \psi^m (\mu_j) P^m_n (\mu_j) w_j ,$$ (3.120)

where $\mu_j$ denotes the Gaussian grid points in the meridional direction, $w_j$ the Gaussian weight at point $\mu_j$, and $J$ the number of Gaussian grid points from pole to pole. The Gaussian grid points ($\mu_j$) are given by the roots of the Legendre polynomial $P_J(\mu)$, and the corresponding weights are given by

$$w_j = \frac{2(1 - \mu_j^2)}{\left[J\; P_{J-1} (\mu_j) \right]^2} .$$ (3.121)

The weights themselves satisfy

$$\sum^J_{j=1} w_j = 2.0 .$$ (3.122)

The Gaussian grid used for the north-south transformation is generally chosen to allow un-aliased computations of quadratic terms only. In this case, the number of Gaussian latitudes $J$ must satisfy

$$J \geq (2N + K + M + 1)/2 \quad \hbox{for}\> M \leq 2(K - N)\>,$$ (3.123)

$$J \geq (3K + 1)/2 \quad \hbox{for}\> M \geq 2(K - N)\> .$$ (3.124)

For the common truncations, these become

$$J \geq (3K + 1)/2 \quad \mathrm{for \; triangular \; and \; trapezoidal} ,$$ (3.125)

$$J \geq (3N + 2M + 1)/2 \quad \mathrm{for \; rhomboidal} .$$ (3.126)

In order to allow exact Fourier transform of quadratic terms, the number of points $P$ in the east-west direction must satisfy

$$P \geq 3M + 1 .$$ (3.127)

The actual values of $J$ and $P$ are often not set equal to the lower limit in order to allow use of more efficient transform programs.

Although in the next section of this model description, we continue to indicate the Gaussian quadrature as a sum from pole to pole, the code actually deals with the symmetric and antisymmetric components of variables and accumulates the sums from equator to pole only. The model requires an even number of latitudes to easily use the symmetry conditions. This may be slightly inefficient for some spectral resolutions. We define a new index, which goes from $-I$ at the point next to the south pole to $+I$ at the point next to the north pole and not including 0 (there are no points at the equator or pole in the Gaussian grid), i.e., let $I = J/2$ and $i = j - J/2$ for $j \geq J/2+1$ and $i = j - J/2 - 1$ for $j \leq J/2$; then the summation in (3.120) can be rewritten as

$$\psi^m_n = \sum \limits^{I}_{i = -I, \;i \neq 0} \psi^m (\mu_i) P^m_n (\mu_i) w_i .$$ (3.128)

The symmetric (even) and antisymmetric (odd) components of $\psi^m$ are defined by

$$\left(\psi_E\right)^m_i = \frac{1}{2} \left( \psi^m_i + \psi^m_{-i} \right),$$

$$\left(\psi_O\right)^m_i = \frac{1}{2} \left(\psi^m_i - \psi^m_{-i} \right) .$$ (3.129)

Since $w_i$ is symmetric about the equator, (3.128) can be rewritten to give formulas for the coefficients of even and odd spherical harmonics:

$$\psi^{m}_{n} = \begin{cases}\sum \limits ^I_{i=1} \left(\psi_E\ri... ...\psi_O\right)^m_i (\mu_i) P^m_n (\mu_i) 2w_i & \text{for $n-m$odd}. \end{cases}$$ (3.130)

The model uses the spectral transform method [118] for all nonlinear terms. However, the model can be thought of as starting from grid-point values at time $t$ (consistent with the spectral representation) and producing a forecast of the grid-point values at time $t + \Delta t$ (again, consistent with the spectral resolution). The forecast procedure involves computation of the nonlinear terms including physical parameterizations at grid points; transformation via Gaussian quadrature of the nonlinear terms from grid-point space to spectral space; computation of the spectral coefficients of the prognostic variables at time $t + \Delta t$ (with the implied spectral truncation to the model resolution); and transformation back to grid-point space. The details of the equations involved in the various transformations are given in the next section.

3.1.11 Combination of terms

In order to describe the transformation to spectral space, for each equation we first group together all undifferentiated explicit terms, all explicit terms with longitudinal derivatives, and all explicit terms with meridional derivatives appearing in the Dyn operator. Thus, the vorticity equation (3.95) is rewritten

$$\underline{\left(\zeta + f \right)}^{n+1} = \underline{\hbox{\sff... ...{\partial}{\partial \mu} (\underline{\hbox{\sffamily\slshape V}}_\mu) \right] ,$$ (3.131)

where the explicit forms of the vectors $\underline{\hbox{\sffamily\slshape V}}, \underline{\hbox{\sffamily\slshape V}}_\lambda,$ and $\underline{\hbox{\sffamily\slshape V}}_\mu$ are given as

$$\underline{\hbox{\sffamily\slshape V}} = \underline{(\zeta+f)}^{n-1} ,$$ (3.132)

$$\underline{\hbox{\sffamily\slshape V}}_\lambda = 2 \Delta t \underline{n}^{n}_{V} ,$$ (3.133)

$$\underline{\hbox{\sffamily\slshape V}}_{\mu} = 2 \Delta t \underline{n}^{n}_{U}.$$ (3.134)

The divergence equation (3.96) is

$$\underline{\delta}^{n+1} = \underline{\hbox{\sffamily\slshape D}} + \frac{1}{a(1 - \mu^2)} \... ...shape D}}_\mu) \right] - \nabla^2 \underline{\hbox{\sffamily\slshape D}}_\nabla$$

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

The mean component of the temperature is not included in the next-to-last term since the Laplacian of it is zero. The thermodynamic equation (3.98) is

$$\underline{T}^{\prime n+1} = \underline{\hbox{\sffamily\slshap... ...T}}_\mu) - \right] - \Delta t {\boldsymbol {D}}^r \; \underline{\delta}^{n+1} .$$ (3.136)

The surface-pressure tendency (3.98) is

$$\Pi^{n+1} = {\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}... ...ta t}{\pi^r} \left( \underline{\Delta p}^r \right)^T \underline{\delta}^{n+1} .$$ (3.137)

The grouped explicit terms in (3.135)-(3.137) are given as follows. The terms of (3.135) are

$$\underline{\hbox{\sffamily\slshape D}} = \underline{\delta}^{n-1} ,$$ (3.138)

$$\underline{\hbox{\sffamily\slshape D}}_\lambda = 2 \Delta t \underline{n}^{n}_{U} ,$$ (3.139)

$$\underline{\hbox{\sffamily\slshape D}}_\mu = 2 \Delta t \underline{n}^{n}_{V} ,$$ (3.140)

$${\underline{\hbox{\sffamily\slshape D}}_\nabla = 2 \Delta t \left... ...hi_s \underline{1} + R {\boldsymbol {H}}^{r} \underline{{\cal T}}^{'n} \right]}$$

$$+ \Delta t \left[ R {\boldsymbol {H}}^{r} \left( {( \underline{T}... ... \underline{h}^{r} \right) \left( {\Pi}^{n-1} - 2 {\Pi}^{n} \right) \right] .$$ (3.141)

The terms of (3.136) are

$$\underline{\hbox{\sffamily\slshape T}} = \left(\underline{T}'\right)^{n-1} + 2 \Delta t \underline{\Gam... ...l {D}}^{r} \left[\underline{\delta}^{n-1} - 2\underline{\delta}^{n} \right] ,$$ (3.142)

$$\underline{\hbox{\sffamily\slshape T}}_{\lambda} = 2 \Delta t \underline{\left(UT'\right)}^n ,$$ (3.143)

$$\underline{\hbox{\sffamily\slshape T}}_\mu = 2 \Delta t \underline{\left(VT'\right)}^n .$$ (3.144)

The nonlinear term in (3.137) is

$${\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}} = \Pi^{n-1... ...ne{\boldsymbol {V}}^n \right)^T \nabla \Pi^n \pi^n \underline{\Delta B} \right]$$

$$- \Delta t \left[ \left(\underline{\Delta p}^r \right)^T \frac{1}... ...r} \right] \left[ \underline{\delta}^{n-1} - 2 \underline{\delta}^n \right] .$$ (3.145)

3.1.12 Transformation to spectral space

Formally, Equations (3.131)-(3.137) are transformed to spectral space by performing the operations indicated in (3.146) to each term. We see that the equations basically contain three types of terms, for example, in the vorticity equation the undifferentiated term $\underline{\hbox{\sffamily\slshape V}}$, the longitudinally differentiated term $\underline{\hbox{\sffamily\slshape V}}_\lambda$, and the meridionally differentiated term $\underline{\hbox{\sffamily\slshape V}}_\mu$. All terms in the original equations were grouped into one of these terms on the Gaussian grid so that they could be transformed at once.

Transformation of the undifferentiated term is obtained by straightforward application of (3.118)-(3.120),

$$\left\{ \underline{\hbox{\sffamily\slshape V}} \right\}^m_n = \sum^J_{j=1} \underline{\hbox{\sffamily\slshape V}}^m(\mu_j) P^m_n(\mu_j) w_j ,$$ (3.146)

where $\underline{\hbox{\sffamily\slshape V}}^m(\mu_j)$ is the Fourier coefficient of $\underline{\hbox{\sffamily\slshape V}}$ with wavenumber $m$ at the Gaussian grid line $\mu_j$. The longitudinally differentiated term is handled by integration by parts, using the cyclic boundary conditions,

$$\left\{ \frac{\partial}{\partial \lambda} (\underline{\hbox{\sffamily\slshape V}}_\lambda) \right\}^m = \frac{1}{2 \pi} \int ^{2\pi}_o \frac{\partial \underline{\hbox{\sffamily\slshape V}}_\lambda} {\partial \lambda} e^{-im\lambda} d \lambda ,$$ (3.147)

$$= im \frac{1}{2 \pi} \int^{2 \pi}_o \underline{\hbox{\sffamily\slshape V}}_\lambda e^{-im\lambda} d\lambda ,$$ (3.148) (3.149)

so that the Fourier transform is performed first, then the differentiation is carried out in spectral space. The transformation to spherical harmonic space then follows (3.152):

$$\left \{ \frac{1}{a(1 - \mu^2)} \frac{\partial}{\partial \lambda}... ...ffamily\slshape V}}_\lambda^m (\mu_j) \frac{P^m_n(\mu_j)}{a(1 - \mu^2_j)} w_j ,$$ (3.150)

where $\underline{\hbox{\sffamily\slshape V}}_\lambda^m (\mu_j)$ is the Fourier coefficient of $\underline{\hbox{\sffamily\slshape V}}_\lambda$ with wavenumber $m$ at the Gaussian grid line $\mu_j$.

The latitudinally differentiated term is handled by integration by parts using zero boundary conditions at the poles:

$$\left \{ \frac{1}{a(1 - \mu^2)} (1 - \mu^2) \frac{\partial}{\partial \mu} (\underline{\hbox{\sffamily\slshape V}}_\mu) \right \}^m_n = \int^1_{-1} \frac{1}{a(1 - \mu^2)} (1 - \mu^2) \frac{\partial} {\partial \mu} (\underline{\hbox{\sffamily\slshape V}}_\mu)^m P^m_n d\mu ,$$ (3.151)

$$= - \int^1_{-1} \frac{1}{a(1 - \mu^2)} (\underline{\hbox{\sffamily\slshape V}}_\mu)^m (1 - \mu^2) \frac{dP^m_n}{d\mu} d\mu .$$ (3.152)

Defining the derivative of the associated Legendre polynomial by

$$H^m_n = (1 - \mu^2) \frac{dP^m_n}{d\mu} ,$$ (3.153)

(3.155) can be written

$$\left \{ \frac{1}{a(1 - \mu^2)} (1 - \mu^2) \frac{\partial}{\part... ...ne{\hbox{\sffamily\slshape V}}_\mu)^m \frac{H^m_n(\mu_j)}{a(1 - \mu^2_j)} w_j .$$ (3.154)

Similarly, the $\nabla^2$ operator in the divergence equation can be converted to spectral space by sequential integration by parts and then application of the relationship

$$\nabla^2 P^m_n( \mu) e^{im \lambda} = \frac{-n (n+1)}{a^2} P^m_n(\mu) e^{im \lambda},$$ (3.155)

to each spherical harmonic function individually so that

$$\left \{ \nabla^2 \underline{\hbox{\sffamily\slshape D}}_\nabla \... ...j=1} \underline{\hbox{\sffamily\slshape D}}_\nabla^m (\mu_j) P^m_n(\mu_j) w_j ,$$ (3.156)

where $\underline{\hbox{\sffamily\slshape D}}_\nabla^m (\mu)$ is the Fourier coefficient of the original grid variable $\underline{\hbox{\sffamily\slshape D}}_\nabla$.

3.1.13 Solution of semi-implicit equations

The prognostic equations can be converted to spectral form by summation over the Gaussian grid using (3.146), (3.150), and (3.154). The resulting equation for absolute vorticity is

$$\underline{\left(\zeta + f \right)}^m_n = \underline {\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}_n^m ,$$ (3.157)

where $\underline{\left(\zeta + f \right)}_n^m$ denotes a spherical harmonic coefficient of $\underline{\left(\zeta + f \right)}^{n+1}$, and the form of $\underline{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}_n^m$, as a summation over the Gaussian grid, is given as

$$\underline{\hbox{\sffamily\slshape V}\hbox{\sffamily\slshape S}}^... ...ly\slshape V}}^m_\mu (\mu_j) \frac{H^m_n (\mu_j)}{a(1 - \mu^2_j)} \right] w_j .$$ (3.158)

The spectral form of the divergence equation (3.135) becomes

$$\underline{\delta}^m_n = \underline{\hbox{\sffamily\slshape D}\hb... ... m}_n + R \left( \underline{b}^r + \underline{h}^r \right) \Pi^m_n \right] ,$$ (3.159)

where $\underline{\delta}^m_n, \;\underline{T}'^m_n, \;$ and $\Pi^m_n$ are spectral coefficients of $\underline{\delta}^{n+1}, \; \underline{T}^{\prime n+1}$, and $\Pi^{n+1}$. The Laplacian of the total temperature in (3.135) is replaced by the equivalent Laplacian of the perturbation temperature in (3.159). $\underline{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}^m_n$ is given by

$$\underline{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}}^... ... \underline {\hbox{\sffamily\slshape D}}^m_\nabla (\mu_j) \right] P^m_n (\mu_j)$$

$$+ im \underline {\hbox{\sffamily\slshape D}}^m_\lambda (\mu_j) \f... ...slshape D}}^m_\mu (\mu_j) \frac{H^m_n (\mu_j)} {a(1 - \mu^2_j)} \Biggr \} w_j .$$ (3.160)

The spectral thermodynamic equation is

$$\underline{T}'^m_n = \underline{\hbox{\sffamily\slshape T}\hbox{\sffamily\slshape S}}^m_n - \Delta t {\boldsymbol {D}}^r \underline{\delta}^m_n ,$$ (3.161)

with $\underline{\hbox{\sffamily\slshape T}\hbox{\sffamily\slshape S}}^m_n$ defined as

$$\underline{\hbox{\sffamily\slshape T}\hbox{\sffamily\slshape S}}^... ...\slshape T}}^m_\mu (\mu_j) \frac{H^m_n (\mu_j)}{a (1 - \mu^2_j)} \right ] w_j ,$$ (3.162)

while the surface pressure equation is

$$\Pi^m_n = {\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}}^... ...ine{\delta}^m_n \left(\underline{\Delta p}^r \right)^T \frac{\Delta t}{\pi^r} ,$$ (3.163)

where ${\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}}^m_n$ is given by

$${\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}}^m_n = \sum... ...x{\sffamily\slshape P}\hbox{\sffamily\slshape S}}^m (\mu_j) P^m_n (\mu_j) w_j .$$ (3.164)

Equation (3.157) for vorticity is explicit and complete at this point. However, the remaining equations (3.159)-(3.163) are coupled. They are solved by eliminating all variables except $\underline{\delta}^m_n$:

$${\boldsymbol {A}}_n \underline{\delta}^m_n = \underline{\hbox{\sffamily\slshape D}\hbox{\sffamily\slshape S}... ...\right) ({\hbox{\sffamily\slshape P}\hbox{\sffamily\slshape S}})^m_n \right ] ,$$ (3.165)
where

$${\boldsymbol {A}}_n = {\boldsymbol {I}} + \Delta t^2 \frac{n(n+1)}{a^2} \left [ R {\b... ...( \left( \underline{\Delta p^r} \right)^T \frac{1}{\pi^r} \right) \right ] ,$$ (3.166)

which is simply a set of $K$ simultaneous equations for the coefficients with given wavenumbers ($m,n$) at each level and is solved by inverting ${\boldsymbol {A}}_n$. In order to prevent the accumulation of round-off error in the global mean divergence (which if exactly zero initially, should remain exactly zero) $\left({\boldsymbol {A}}_o\right)^{-1}$ is set to the null matrix rather than the identity, and the formal application of (3.165) then always guarantees $\underline{\delta}^o_o = 0$. Once $\delta^m_n$ is known, $\underline{T}'^m_n$ and $\Pi^m_n$ can be computed from (3.161) and (3.163), respectively, and all prognostic variables are known at time $n\!+\!1$ as spherical harmonic coefficients. Note that the mean component $\underline{T}'^o_o$ is not necessarily zero since the perturbations are taken with respect to a specified $\underline{T}^r$.

3.1.14 Horizontal diffusion

As mentioned earlier, the horizontal diffusion in (3.88) and (3.91) is computed implicitly via time splitting after the transformations into spectral space and solution of the semi-implicit equations. In the following, the $\zeta$ and $\delta$ equations have a similar form, so we write only the $\delta$ equation:

$$\left(\delta^*\right)^m_n = \left(\delta^{n+1}\right)^m_n - \left( -1 \right)^i 2 \Delta t ... ... \left( - 1 \right)^i \left(\delta^*\right)^m_n \left(2/a^2 \right)^i \right] ,$$ (3.167)

$$\left(T^*\right)^m_n = \left(T^{n+1}\right)^m_n - \left( -1\right)^i 2 \Delta t K^{(2i)} \left[ \nabla^{2i} \left(T^*\right)^m_n \right] .$$ (3.168)

The extra term is present in (3.167), (3.171) and (3.173) to prevent damping of uniform rotations. The solutions are just

$$\left(\delta^*\right)^m_n = K^{(2i)}_n \left(\delta\right) \left(\delta^{n+1}\right)^m_n ,$$ (3.169)

$$\left(T^*\right)^m_n = K^{(2i)}_n \left(T \right) \left(T^{n+1}\right)^m_n ,$$ (3.170)

$$K^{(2)}_n \left(\delta \right) = \left\{ 1 + 2 \Delta t D_n K^{(2)} \left[ \left( \frac{n(n + 1)}{a^2} \right) - \frac{2}{a^2} \right] \right\}^{-1} ,$$ (3.171)

$$K^{(2)}_n \left( T \right) = \left\{ 1 + 2\Delta t D_n K^{(2)} \left( \frac{n(n + 1)}{a^2} \right) \right\}^{-1} ,$$ (3.172)

$$K^{(4)}_n \left(\delta\right) = \left\{ 1 + 2 \Delta t D_n K^{(4)} \left[ \left( \frac{n(n+1)}{a^2} \right)^2 - \frac{4}{a^4} \right] \right\}^{-1} ,$$ (3.173)

$$K^{(4)}_n \left(T \right) = \left\{ 1 + 2 \Delta t D_n K^{(4)} \left( \frac{n(n+1)}{a^2} \right)^2 \right\}^{-1} .$$ (3.174)

$K^{(2)}_n \left(\delta \right)$ and $K^{(4)}_n \left(\delta \right)$ are both set to 1 for $n$ = 0. The quantity $D_n$ represents the ``Courant number limiter'', normally set to 1. However, $D_n$ is modified to ensure that the CFL criterion is not violated in selected upper levels of the model. If the maximum wind speed in any of these upper levels is sufficiently large, then $D_n = 1000$ in that level for all $n > n_c$, where $n_c = a \Delta t \big/ \max \vert {\boldsymbol {V}} \vert$. This condition is applied whenever the wind speed is large enough that $n_c < K$, the truncation parameter in (3.115), and temporarily reduces the effective resolution of the model in the affected levels. The number of levels at which this ``Courant number limiter'' may be applied is user-selectable, but it is only used in the top level of the 26 level CAM 3.0 control runs.

The diffusion of $T$ is not complete at this stage. In order to make the partial correction from $\eta$ to $p$ in (3.82) local, it is not included until grid-point values are available. This requires that $\nabla^4 \Pi$ also be transformed from spectral to grid-point space. The values of the coefficients $K^{(2)}$ and $K^{(4)}$ for the standard T42 resolution are $2.5 \times 10^5$m$^2$sec$^{-1}$ and $1.0 \times 10^{16}$m$^4$sec$^{-1}$, respectively.

3.1.15 Initial divergence damping

Occasionally, with poorly balanced initial conditions, the model exhibits numerical instability during the beginning of an integration because of excessive noise in the solution. Therefore, an optional divergence damping is included in the model to be applied over the first few days. The damping has an initial e-folding time of $\Delta t$ and linearly decreases to 0 over a specified number of days, $t_D$, usually set to be 2. The damping is computed implicitly via time splitting after the horizontal diffusion.

$$r = \max \left[ \frac{1}{\Delta t} (t_D - t) / t_D ,  0 \right]$$ (3.175)

$$\left(\delta^*\right)^m_n = \frac{1}{1 + 2\Delta t r} \left(\delta^*\right)^m_n$$ (3.176)

3.1.16 Transformation from spectral to physical space

After the prognostic variables are completed at time $n+1$ in spectral space $\left(\underline{(\zeta + f)^*}\right)^m_n$, $\left(\underline{\delta^*}\right)^m_n$, $\left(\underline{T^*}\right)^m_n$, $\left(\Pi^{n+1}\right)^m_n$ they are transformed to grid space. For a variable $\psi$, the transformation is given by

$$\psi( \lambda, \mu) = \sum^M_{m=-M} \left [ \sum^{{\cal N}(m)}_{n=\vert m\vert} \psi^m_n P^m_n (\mu) \right ] e^{im\lambda} .$$ (3.177)

The inner sum is done essentially as a vector product over $n$, and the outer is again performed by an FFT subroutine. The term needed for the remainder of the diffusion terms, $\nabla^4 \Pi$, is calculated from

$$\nabla^4 \Pi^{n+1} = \sum^M_{m=-M} \left [\sum^{{\cal N}(m)}_{n=\... ...^2} \right )^2 \left(\Pi^{n+1} \right)^m_n P^m_n (\mu) \right ] e^{im\lambda} .$$ (3.178)

In addition, the derivatives of $\Pi$ are needed on the grid for the terms involving $\nabla \Pi$ and ${\boldsymbol {V}} \cdot \nabla \Pi$,

$${\boldsymbol {V}} \cdot \nabla \Pi = \frac{U}{a(1 - \mu^2)} \frac... ...lambda} + \frac{V}{a(1 - \mu^2)} (1 - \mu^2) \frac{\partial \Pi}{\partial \mu}.$$ (3.179)

These required derivatives are given by

$$\frac{\partial \Pi}{\partial \lambda} = \sum^M_{m=-M} im \left[ \sum^{{\cal N}(m)}_{n=\vert m \vert} \Pi^m_n P^m_n (\mu) \right] e^{im\lambda} ,$$ (3.180)
and using (3.153),

$$(1 - \mu^2) \frac{\partial \Pi}{\partial \mu} = \sum^M_{m=-M} \left [ \sum^{{\cal N}(m)}_{n=\vert m \vert} \Pi^m_n H^m_n (\mu) \right ] e^{im \lambda} ,$$ (3.181)

which involve basically the same operations as (3.178). The other variables needed on the grid are $U$ and $V$. These can be computed directly from the absolute vorticity and divergence coefficients using the relations

$$\left(\zeta + f \right)^m_n = - \frac{n(n+1)}{a^2} \psi^m_n + f^m_n ,$$ (3.182)

$$\delta^m_n = - \frac{n(n+1)}{a^2} \chi^m_n ,$$ (3.183)

in which the only nonzero $f^m_n$ is $f^o_1 = \Omega/ \sqrt{.375},$ and

$$U = \frac{1}{a} \frac{\partial \chi}{\partial \lambda} - \frac{(1 - \mu^2)}{a} \frac{\partial \psi}{\partial \mu} ,$$ (3.184)

$$V = \frac{1}{a} \frac{\partial \psi}{\partial \lambda} + \frac{(1 - \mu^2)} {a} \frac{\partial \chi}{\partial \mu} .$$ (3.185)

Thus, the direct transformation is

$$U = - \sum^M_{m=-M} a \sum^{{\cal N}(m)}_{n=\vert m\vert} \left [ \fr... ...m_n (\mu) - \frac{1}{n(n+1)} (\zeta + f)^m_n H^m_n (\mu) \right] e^{im \lambda}$$

$$\phantom{=} - \>\frac{a}{2} \frac{\Omega}{\sqrt{0.375}} H^o_1 ,$$ (3.186)

$$V = - \sum^M_{m=-M} a \sum^{{\cal N}(m)}_{n=\vert m\vert} \bigg [\fra... ..._n P^m_n (\mu) + \frac{1}{n(n+1)} \delta^m_n H^m_n (\mu) \bigg ] e^{im\lambda}.$$ (3.187)

The horizontal diffusion tendencies are also transformed back to grid space. The spectral coefficients for the horizontal diffusion tendencies follow from (3.167) and (3.168):

$$F_{T_H}\left(T^*\right)^m_n = \left(-1\right)^{i+1} K^{2i} \left[ \nabla^{2i} \left(T^*\right) \right]^m_n ,$$ (3.188)

$$F_{\zeta_H} \left(\left(\zeta + f \right)^* \right)^m_n = \left(-1\right)^{i+1} K^{2i} \left \{\nabla^{2i} \left(\zeta + ... ...- \left(-1\right)^i \left(\zeta + f \right)^* \left(2/a^2 \right)^i \right \} ,$$ (3.189)

$$F_{\delta_H} \left(\delta^*\right)^m_n = \left(-1\right) K^{2i} \left\{ \nabla^{2i} \left(\delta ^*\right) - \left(-1\right)^i \delta^* \left(2/a^2 \right)^i \right\},$$ (3.190)

using $i = 1$ or 2 as appropriate for the $\nabla^2$ or $\nabla^4$ forms. These coefficients are transformed to grid space following (3.114) for the $T$ term and (3.186) and (3.187) for vorticity and divergence. Thus, the vorticity and divergence diffusion tendencies are converted to equivalent $U$ and $V$ diffusion tendencies.

3.1.17 Horizontal diffusion correction

After grid-point values are calculated, frictional heating rates are determined from the momentum diffusion tendencies and are added to the temperature, and the partial correction of the $\nabla^4$ diffusion from $\eta$ to $p$ surfaces is applied to $T$. The frictional heating rate is calculated from the kinetic energy tendency produced by the momentum diffusion

$$F_{F_H} = -u^{n-1} F_{u_H} (u^*)/c_p^* -v^{n-1} F_{v_H} (v^*)/c^*_p ,$$ (3.191)

where $F_{u_H}$, and $F_{v_H}$ are the momentum equivalent diffusion tendencies, determined from $F_{\zeta_H}$ and $F_{\delta_H}$ just as $U$ and $V$ are determined from $\zeta$ and $\delta$, and

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

These heating rates are then combined with the correction,

$$\hat {T}^{n+1}_k = T^*_k + \left( 2 \Delta t F_{F_H}\right)_k + 2... ...t( \pi B \frac{\partial T^*}{\partial p} \right)_k K^{(4)} \nabla^4 \Pi^{n+1} .$$ (3.193)

The vertical derivatives of $T^*$ (where the $^*$ notation is dropped for convenience) are defined by

$$\left( \pi B \frac{\partial T}{\partial p} \right)_1 = \frac{\pi} {2\Delta p_1} \left[ B_{1+\frac{1}{2}} \left( T_2 - T_1 \right) \right] ,$$ (3.194)

$$\left(\pi B \frac{\partial T}{\partial p} \right)_k = \frac{\pi} {2\Delta p_k} \left[ B_{k + \frac{1}{2}} \left( T_{k+1} - T_k \right) + B_{k - \frac{1}{2}} \left( T_k - T_{k-1} \right) \right] ,$$ (3.195)

$$\left(\pi B \frac{\partial T}{\partial p} \right)_K = \frac{\pi} {2\Delta p_K} \left[ B_{K - \frac{1}{2}} \left( T_K - T_{K-1} \right) \right] .$$ (3.196)

The corrections are added to the diffusion tendencies calculated earlier (3.188) to give the total temperature tendency for diagnostic purposes:

$$\hat{F}_{T_H}(T^*)_k = F_{T_H}(T^*)_k + \left( 2 \Delta t F_{F_H}... ...left(\pi \frac{\partial T^*}{\partial p} \right)_k K^{(4)} \nabla^4 \Pi^{n+1} .$$ (3.197)

3.1.18 Semi-Lagrangian Tracer Transport

The forecast equation for water vapor specific humidity and constituent mixing ratio in the $\eta$ system is from (3.36) excluding sources and sinks.

$$\frac{dq}{dt} = \frac{\partial q}{\partial t} + {\boldsymbol {V}} \cdot \nabla ... ...t\eta \frac{\partial p}{\partial \eta} \frac{\partial q} {\partial p} = 0$$ (3.198)
or

$$\frac{dq}{dt} = \frac{\partial q}{\partial t} + {\boldsymbol {V}} \cdot \nabla q + \dot\eta \frac{\partial q}{\partial \eta} = 0 .$$ (3.199)

Equation (3.199) is more economical for the semi-Lagrangian vertical advection, as $\Delta \eta$ does not vary in the horizontal, while $\Delta p$ does. Written in this form, the $\eta$ advection equations look exactly like the $\sigma$ equations.

The parameterizations are time-split in the moisture equation. The tendency sources have already been added to the time level $\left( n - 1 \right)$. The semi-Lagrangian advection step is subdivided into horizontal and vertical advection sub-steps, which, in an Eulerian form, would be written

$$q^* = q^{n-1} + 2\Delta t \left( {\boldsymbol {V}} \cdot \nabla q \right)^n$$ (3.200)
and

$$q^{n+1} = q^* + 2\Delta t \left( \dot\eta \frac{\partial q}{\partial n} \right)^n .$$ (3.201)

In the semi-Lagrangian form used here, the general form is

$$q^* = {\rm L}_{\lambda \varphi} \left( q^{n-1} \right) ,$$ (3.202)

$$q^{n+1} = {\rm L}_\eta \left( q^* \right) .$$ (3.203)

Equation (3.202) represents the horizontal interpolation of $q^{n-1}$ at the departure point calculated assuming $\dot\eta = 0$. Equation (3.203) represents the vertical interpolation of $q^*$ at the departure point, assuming ${\boldsymbol {V}} = 0$.

The horizontal departure points are found by first iterating for the mid-point of the trajectory, using winds at time $n$, and a first guess as the location of the mid-point of the previous time step

$$\lambda^{k+1}_M = \lambda_A - \Delta t u^n \left( \lambda^k_M, \varphi^k_M \right) \big/a \cos \varphi^k_M ,$$ (3.204)

$$\varphi^{k+1}_M = \varphi_A - \Delta t v^n \left( \lambda^k_M, \varphi^k_M \right)/a ,$$ (3.205)

where subscript $A$ denotes the arrival (Gaussian grid) point and subscript $M$ the midpoint of the trajectory. The velocity components at $\left( \lambda_M^k, \varphi^k_M \right)$ are determined by Lagrange cubic interpolation. For economic reasons, the equivalent Hermite cubic interpolant with cubic derivative estimates is used at some places in this code. The equations will be presented later.

Once the iteration of (3.204) and (3.205) is complete, the departure point is given by

$$\lambda_D = \lambda_A - 2 \Delta t u^n \left( \lambda_M, \varphi_M \right) \big/a \cos \varphi_M ,$$ (3.206)

$$\varphi_D = \lambda_A - 2 \Delta t v^n \left( \lambda_M, \varphi_M \right)/a ,$$ (3.207)

where the subscript $D$ denotes the departure point.

The form given by (3.204)-(3.207) is inaccurate near the poles and thus is only used for arrival points equatorward of 70$^\circ$ latitude. Poleward of 70$^\circ$ we transform to a local geodesic coordinate for the calculation at each arrival point. The local geodesic coordinate is essentially a rotated spherical coordinate system whose equator goes through the arrival point. Details are provided in Williamson and Rasch [186]. The transformed system is rotated about the axis through $\left( \lambda_A - \frac{\pi}{2}, 0 \right)$ and $\left( \lambda_A + \frac{\pi}{2}, 0 \right)$, by an angle $\varphi_A$ so the equator goes through $\left( \lambda_A, \varphi_A \right)$. The longitude of the transformed system is chosen to be zero at the arrival point. If the local geodesic system is denoted by $\left( \lambda^\prime, \varphi^\prime \right)$, with velocities $\left( u^\prime, v^\prime \right)$, the two systems are related by

$$\sin \phi^\prime = \sin \phi \cos \phi_A - \cos \phi \sin \phi_A \cos \left( \lambda_A - \lambda \right) ,$$ (3.208)

$$\sin \phi = \sin \phi^\prime \cos \phi_A + \cos \phi^\prime \sin \prime_A \cos \lambda^\prime ,$$ (3.209)

$$\sin \lambda^\prime \cos \phi^\prime = -\sin \left( \lambda_A -\lambda \right) \cos \phi ,$$ (3.210)

$$v^\prime \cos \phi^\prime = v \left[ \cos \phi \cos \phi_A + \sin \phi \sin \phi_A \cos \left( \lambda_A - \lambda \right) \right]$$

$$\phantom{=} - u \sin \phi_A \sin \left( \lambda_A - \lambda \right),$$ (3.211)

$$u^\prime \cos \lambda^\prime - v^\prime \sin \lambda^\prime \sin \phi^\prime = u \cos \left( \lambda_A - \lambda \right) + v \sin \phi \sin \left( \lambda_A - \lambda \right) .$$ (3.212)

The calculation of the departure point in the local geodesic system is identical to (3.204)-(3.207) with all variables carrying a prime. The equations can be simplified by noting that $\left( \lambda^\prime_A, \varphi_A^\prime \right) = \left( 0, 0 \right)$ by design and $u^\prime \left( \lambda^\prime_A, \varphi_A^\prime \right) = u\left( \lambda_A, \varphi_A \right)$ and $v^\prime \left( \lambda^\prime_A, \varphi_A^\prime \right) = v\left( \lambda_A, \varphi_A \right)$. The interpolations are always done in global spherical coordinates.

The interpolants are most easily defined on the interval 0 $\leq \theta \leq 1$. Define

$$\theta = \left( x_D - x_i \right) \big/ \left( x_{i+1} - x_i \right),$$ (3.213)

where $x$ is either $\lambda$ or $\varphi$ and the departure point $x_D$ falls within the interval $\left( x_i, x_{i+1} \right)$. Following (23) of [146] with $r_i=3$ the Hermite cubic interpolant is given by

$$q_D = q_{i+1} \left[ 3-2\theta \right]\theta^2 - d_{i+1} \left[ h_i \theta^2 \left( 1-\theta \right) \right]$$

$$\phantom{=} + q_i \left[ 3-2\left(1-\theta\right) \right] \left(1-\theta \right)^2 + d_i \left[ h_i \theta \left( 1-\theta \right)^2 \right]$$ (3.214)

where $q_i$ is the value at the grid point $x_i$, $d_i$ is the derivative estimate given below, and $h_i = x_{i+1} - x_i$.

Following (3.2.12) and (3.2.13) of Hildebrand [71], the Lagrangian cubic polynomial interpolant used for the velocity interpolation, is given by

$$f_{D\phantom{\dot D}} = \sum^2_{j=-1} \ell_j \left( x_D \right) f_{i+j}$$ (3.215)

where

$$\ell_j \left( x_D \right) = \frac{ \left( x_D - x_{i-1} \right) \... ...t) \left( x_{i+j} - x_{i+j+1} \right) \ldots \left( x_{i+j} - x_{i+2} \right) }$$ (3.216)

where $f$ can represent either $u$ or $v$, or their counterparts in the geodesic coordinate system.

The derivative approximations used in (3.214) for $q$ are obtained by differentiating (3.215) with respect to $x_D$, replacing $f$ by $q$ and evaluating the result at $x_D$ equal $x_{i}$ and $x_{i+1}$. With these derivative estimates, the Hermite cubic interpolant (3.214) is equivalent to the Lagrangian (3.215). If we denote the four point stencil $\left( x_{i-1},x_{i},x_{i+1},x_{i+2} \right)$ by $\left( x_1,x_2,x_3,x_4, \right)$ the cubic derivative estimates are

$$d_2 = \left[ \frac{(x_2 - x_3)( x_2 - x_4 )} {(x_1 - x_2)(x_1 - x_3)(x_1 - x_4)} \right] q_1$$ (3.217)

$$- \left[ \frac{1}{(x_1 - x_2)} - \frac{1}{(x_2 - x_3)} - \frac{1}{(x_2 - x_4)} \right] q_2$$ (3.218)

$$+ \left[ \frac{(x_2 - x_1)(x_2 - x_4)} {(x_1 - x_3)(x_2 - x_3)(x_3 - x_4)}\right] q_3$$ (3.219)

$$- \left[ \frac{(x_2 - x_1)(x_2 - x_3)} {(x_1 - x_4)(x_2 - x_4)(x_3 - x_4)} \right] q_4$$ (3.220)
and

$$d_3 = \left[ \frac{(x_3 - x_2)(x_3 - x_4)} {(x_1 - x_2)(x_1 - x_3)(x_1 - x_4)} \right] q_1$$ (3.221)

$$- \left[ \frac{(x_3 - x_1)(x_3 - x_4)} {(x_1 - x_2)(x_2 - x_3)(x_2 - x_4)} \right] q_2$$ (3.222)

$$- \left[ \frac{1}{(x_1 - x_3)} + \frac{1}{(x_2 - x_3)} - \frac{1}{(x_3 - x_4)} \right] q_3$$ (3.223)

$$- \left[ \frac{(x_3 - x_1)(x_3 - x_2)} {(x_1 - x_4)(x_2 - x_4)(x_3 - x_4)} \right] q_4$$ (3.224)

The two dimensional $\left( \lambda, \varphi \right)$ interpolant is obtained as a tensor product application of the one-dimensional interpolants, with $\lambda$ interpolations done first. Assume the departure point falls in the grid box $(\lambda_i,\lambda_{i+1})$ and $(\varphi_i,\varphi_{i+1})$. Four $\lambda$ interpolations are performed to find $q$ values at $(\lambda_D,\varphi_{j-1})$, $(\lambda_D,\varphi_j)$, $(\lambda_D,\varphi_{j+1})$, and $(\lambda_D, \varphi_{j+2})$. This is followed by one interpolation in $\varphi$ using these four values to obtain the value at $(\lambda_D, \varphi_D)$. Cyclic continuity is used in longitude. In latitude, the grid is extended to include a pole point (row) and one row across the pole. The pole row is set equal to the average of the row next to the pole for $q$ and to wavenumber 1 components for $u$ and $v$. The row across the pole is filled with the values from the first row below the pole shifted $\pi$ in longitude for $q$ and minus the value shifted by $\pi$ in longitude for $u$ and $v$.

Once the departure point is known, the constituent value of $q^* = q^{n-1}_D$ is obtained as indicated in (3.202) by Hermite cubic interpolation (3.214), with cubic derivative estimates (3.215) and (3.216) modified to satisfy the Sufficient Condition for Monotonicity with C$^\circ$ continuity (SCMO) described below. Define $\Delta_i q$ by

$$\Delta_i q = \frac{q_{i+1} - q_i}{x_{i+1} - x_i} .$$ (3.225)

First, if $\Delta_i q= 0$ then

$$d_i = d_{i+1} = 0 .$$ (3.226)

Then, if either

$$0 \leq \frac{d_i}{\Delta_i q} \leq 3$$ (3.227)

or

$$0 \leq \frac{d_{i+1}}{\Delta_i q} \leq 3$$ (3.228)

is violated, $d_i$ or $d_{i+1}$ is brought to the appropriate bound of the relationship. These conditions ensure that the Hermite cubic interpolant is monotonic in the interval $\left[ x_i, x_{i+1} \right]$.

The horizontal semi-Lagrangian sub-step (3.202) is followed by the vertical step (3.203). The vertical velocity $\dot \eta$ is obtained from that diagnosed in the dynamical calculations (3.94) by

$$\left( \dot\eta \right)_{k+\frac{1}{2}} = \left( \dot\eta \fra... ...{k+\frac{1}{2}} \Bigg/ \left(\frac{p_{k+1} -p_k}{\eta_{k+1} - \eta_k} \right) ,$$ (3.229)

with $\eta_k= A_k + B_k$. Note, this is the only place that the model actually requires an explicit specification of $\eta$. The mid-point of the vertical trajectory is found by iteration

$$\eta^{k+1}_M = \eta_A - \Delta t \dot\eta^n \left( \eta^k_M \right) .$$ (3.230)

Note, the arrival point $\eta_A$ is a mid-level point where $q$ is carried, while the $\dot \eta$ used for the interpolation to mid-points is at interfaces. We restrict $\eta_M$ by

$$\eta_1 \leq \eta_M \leq \eta_K ,$$ (3.231)

which is equivalent to assuming that $q$ is constant from the surface to the first model level and above the top $q$ level. Once the mid-point is determined, the departure point is calculated from

$$\eta_D = \eta_A - 2 \Delta t \dot\eta^n \left( \eta_M \right) ,$$ (3.232)

with the restriction

$$\eta_1 \leq \eta_D \leq \eta_K .$$ (3.233)

The appropriate values of $\dot \eta$ and $q$ are determined by interpolation (3.214), with the derivative estimates given by (3.215) and (3.216) for $i = 2$ to $K - 1$. At the top and bottom we assume a zero derivative (which is consistent with (3.231) and (3.233)), $d_i = 0$ for the interval $k=1$, and $\delta_{i+1} = 0$ for the interval $k = K - 1$. The estimate at the interior end of the first and last grid intervals is determined from an uncentered cubic approximation; that is $d_{i+1}$ at the $k=1$ interval is equal to $d_i$ from the $k=2$ interval, and $d_i$ at the $k = K - 1$ interval is equal to $d_{i+1}$ at the $k = K -2$ interval. The monotonic conditions (3.227) to (3.228) are applied to the $q$ derivative estimates.

3.1.19 Mass fixers

This section describes original and modified fixers used for the Eulerian and semi-Lagrangian dynamical cores.

Let $\pi^0$, $\Delta p^0$ and $q^0$ denote the values of air mass, pressure intervals, and water vapor specific humidity at the beginning of the time step (which are the same as the values at the end of the previous time step.)

$\pi^+$, $\Delta p^+$ and $q^+$ are the values after fixers are applied at the end of the time step.

$\pi^-$, $\Delta p^-$ and $q^-$ are the values after the parameterizations have updated the moisture field and tracers.

Since the physics parameterizations do not change the surface pressure, $\pi^-$ and $\Delta p^-$ are also the values at the beginning of the time step.

The fixers which ensure conservation are applied to the dry atmospheric mass, water vapor specific humidity and constituent mixing ratios. For water vapor and atmospheric mass the desired discrete relations, following Williamson and Olson [185] are

$$\int\limits_2 \pi^+ - \int\limits_3 q^+ \Delta p^+ = {\boldsymbol {P}} ,$$ (3.234)

$$\int\limits_3 q^+ \Delta p^+ = \int\limits_3 q^- \Delta p^- ,$$ (3.235)

where $\boldsymbol {P}$ is the dry mass of the atmosphere. From the definition of the vertical coordinate,

$$\Delta p = p_0 \Delta A + \pi \Delta B,$$ (3.236)

and the integral $\int\limits_2$ denotes the normal Gaussian quadrature while $\int\limits_3$ includes a vertical sum followed by Gaussian quadrature. The actual fixers are chosen to have the form

$$\pi^+ \left( \lambda, \varphi \right) = {\boldsymbol {M}} \hat \pi^+ \left( \lambda, \varphi \right) ,$$ (3.237)

preserving the horizontal gradient of $\Pi$, which was calculated earlier during the inverse spectral transform, and

$$q^+ \left( \lambda, \varphi, \eta \right) = \hat q^+ + \alpha\eta \hat q^+ \vert \hat q^+ - q^- \vert .$$ (3.238)

In (3.237) and (3.238) the $\hat{\left(\hskip 10pt \right)}$ denotes the provisional value before adjustment. The form (3.238) forces the arbitrary corrections to be small when the mixing ratio is small and when the change made to the mixing ratio by the advection is small. In addition, the $\eta$ factor is included to make the changes approximately proportional to mass per unit volume [141]. Satisfying (3.234) and (3.235) gives

$$\alpha = \frac{\int\limits_3 q^- \Delta p^- - \int\limits_3 ... ... A + M \int\limits_3\eta\hat q^+ \vert \hat q^+ - q^- \vert \hat\pi^+ \Delta B}$$ (3.239)

and

$${\boldsymbol {M}} = \left({\boldsymbol {P}} + \int\limits_3 q^- \Delta p^- \right) \Bigg/ \int\limits_2 \hat \pi^+ .$$ (3.240)

Note that water vapor and dry mass are corrected simultaneously. Additional advected constituents are treated as mixing ratios normalized by the mass of dry air. This choice was made so that as the water vapor of a parcel changed, the constituent mixing ratios would not change. Thus the fixers which ensure conservation involve the dry mass of the atmosphere rather than the moist mass as in the case of the specific humidity above. Let $\chi$ denote the mixing ratio of constituents. Historically we have used the following relationship for conservation:

$$\int\limits_3\chi^+ (1-q^+) \Delta p^+ = \int\limits_3 \chi^- (1-q^-)\Delta p^- .$$ (3.241)

The term $(1-q)\Delta p$ defines the dry air mass in a layer. Following Rasch et al. [141] the change made by the fixer has the same form as (3.238)

$$\chi^+ \left( \lambda, \varphi, \eta \right) = \hat \chi^+ + \alpha_\chi \eta \hat \chi^+ \vert \hat \chi^+ - \chi^- \vert .$$ (3.242)

Substituting (3.242) into (3.241) and using (3.237) through (3.240) gives

$$\alpha_\chi = \frac{ \int\limits_3 \chi^- (1-q^-) \Delta p^- - \i... ...vert\hat\chi^+ - \chi^-\vert \eta\hat q^+ \vert \hat q^+ - q^-\vert\Delta p} ,$$ (3.243)

where the following shorthand notation is adopted:

$$\int\limits_{A,B} (  )\Delta p = \int\limits_3 (  ) p_0 \Delta A + M\int\limits_3 (  ) p_s \Delta B .$$ (3.244)

We note that there is a small error in (3.241). Consider a situation in which moisture is transported by a physical parameterization, but there is no source or sink of moisture. Under this circumstance $q^- \ne q^0$, but the surface pressure is not allowed to change. Since $(1- q^-)\Delta p^- \ne (1-q^0)\Delta p^0$, there is an implied change of dry mass of dry air in the layer, and even in circumstances where there is no change of dry mixing ratio $\chi$ there would be an implied change in mass of the tracer. The solution to this inconsistency is to define a dry air mass only once within the model time step, and use it consistently throughout the model. In this revision, we have chosen to fix the dry air mass in the model time step where the surface pressure is updated, e.g. at the end of the model time step. Therefore, we now replace (3.241) with

$$\int\limits_3\chi^+ (1-q^+) \Delta p^+ = \int\limits_3 \chi^- (1-q^0)\Delta p^0 .$$ (3.245)

There is a corresponding change in the first term of the numerator of (3.243) in which $q^-$ is replace by $q^0$. CAM 3.0uses (3.243) for water substances and constituents affecting the temperature field to prevent changes to the IPCC simulations. In the future, constituent fields may use a corrected version of (3.243).

3.1.20 Energy Fixer

Following notation in section 3.1.19, the total energy integrals are

$$\int\limits_3 {1 \over g} \left[ c_p T^+ + \Phi_s + {1 \over 2} \left( {u^+}^2 + {v^+}^2 \right) \right] \Delta p^+ = {\boldsymbol {E}}$$ (3.246)

$${\boldsymbol {E}} = \int\limits_3 {1 \over g} \left[ c_p T^- + \Phi_s + {1 \over 2} \left( {u^-}^2 + {v^-}^2 \right) \right] \Delta p^- + {\boldsymbol {S}}$$ (3.247)

$${\boldsymbol {S}} = \int\limits_2 \left[ \left( FSNT - FLNT \right) - \left( FSNS - FLNS - SHFLX - \rho_{H_2O} L_v PRECT \right) - \right] \Delta t$$ (3.248)

where ${\boldsymbol {S}}$ is the net source of energy from the parameterizations. $FSNT$ is the net downward solar flux at the model top, $FLNT$ is the net upward longwave flux at the model top, $FSNS$ is the net downward solar flux at the surface, $FLNS$ is the net upward longwave flux at the surface, $SHFLX$ is the surface sensible heat flux, and $PRECT$ is the total precipitation during the time step. From equation (3.237)

$$\pi^+ \left( \lambda, \varphi \right) = {\boldsymbol {M}} \hat \pi^+ \left( \lambda, \varphi \right)$$ (3.249)

and from (3.236)

$$\Delta p = p_0 \Delta A + \pi \Delta B$$ (3.250)

The energy fixer is chosen to have the form

$$T^+ \left( \lambda, \varphi, \eta \right) = \hat T^+ + \beta$$ (3.251)

$$u^+ \left( \lambda, \varphi, \eta \right) = \hat u^+$$ (3.252)

$$v^+ \left( \lambda, \varphi, \eta \right) = \hat v^+$$ (3.253)

Then

$$\beta = { g{\boldsymbol {E}} - \int\limits_3 \left[ c_p \hat T... ...3 c_p p_0\Delta A + {\boldsymbol {M}} \int\limits_3 c_p \hat\pi^+ \Delta B }$$ (3.254)

3.1.21 Statistics Calculations

At each time step, selected global average statistics are computed for diagnostic purposes when the model is integrated with the Eulerian and semi-Lagrangian dynamical cores. Let $\int_3$ denote a global and vertical average and $\int_2$ a horizontal global average. For an arbitrary variable $\psi$, these are defined by

$$\int_3 \psi d V = \sum^K_{k=1} \sum^J_{j=1} \sum^I_{i=1} \psi_{ijk} w_j \left(\frac{\Delta p_k}{\pi} \right) \bigg/ 2I ,$$ (3.255)
and

$$[-1.0em]$$

$$[-2.0em] \int_2 \psi dA = \sum^J_{j=1} \sum^I_{i=1} \psi_{ijk} w_j/2I ,$$ (3.256)

where recall that

$$\sum^J_{j=1} w_j = 2 .$$ (3.257)

The quantities monitored are:

$$\; (\zeta+f) (\mathrm{s}^{-1}) = \left[ \int_3 (\zeta^n + f)^2 dV \right]^{1/2} ,$$ (3.258)

$$\; \delta (\mathrm{s}^{-1}) = \left [\int_3 (\delta^n)^2 dV \right]^{1/2} ,$$ (3.259)

$$\; T \; (\mathrm{K}) \; = \left[ \int_3 (T^r + T'^n)^2 dV \right]^{1/2} ,$$ (3.260)

$$\; g \; (\mathrm{Pa}) = \int_2 \pi^{n} dA ,$$ (3.261)

$$\; (\mathrm{kg \; m}^{-2}) = \int_3 \pi^{n} q^{n}/g dV .$$ (3.262)

3.1.22 Reduced grid

The Eulerian core and semi-Lagrangian tracer transport can be run on reduced grids. The term reduced grid generally refers to a grid based on latitude and longitude circles in which the longitudinal grid increment increases at latitudes approaching the poles so that the longitudinal distance between grid points is reasonably constant. Details are provided in [187]. This option provides a saving of computer time of up to 25%.