4.11 Vertical Diffusion and Boundary Layer Processes

The vertical diffusion parameterization in CAM 3.0 provides the interface to the turbulence parameterization, computes the molecular diffusivities (if necessary) and finally computes the tendencies of the input variables. The diffusion equations are actually solved implicitly, so the tendencies are computed from the difference between the final and initial profiles.

In the near future, the gravity wave parameterization will also be called from within the vertical diffusion. This will allow the turbulent and, especially, the molecular diffusivity to be passed to the gravity wave parameterization to damp vertically propagating waves. The gravity wave parameterization may return additional diffusivities and tendencies to be applied before the actual diffusion is applied.

As in CCM2 and CCM3, the turbulence parameterization in CAM 3.0 includes computation of diffusivities for the free atmosphere, based on the gradient Richardson number, and an explicit, non-local Atmospheric Boundary Layer (ABL) parameterization. The ABL parameterization includes a determination of the boundary layer depth. In practice, the free atmosphere diffusivities are calculated first at all levels. The ABL scheme then determines the ABL depth and diffusivities and replaces the free atmosphere values for all levels within the ABL, returning both the updated diffusivities and the non-local transport terms. The implementation of the ABL parameterization in CCM2 is discussed in [74], while the formalism only is discussed here. Following the ABL scheme, molecular diffusivities are computed if the model top extends above $\sim$90 km (0.1 Pa).

As described in Boville and Bretherton [25], a general vertical diffusion parameterization can be written in terms of the divergence of diffusive fluxes:

$$\frac{\partial}{\partial t} (u,v,q) = - \frac{1}{\rho} \frac{\partial}{\partial z} (F_u, F_v, F_q)$$ (4.448)

$$\frac{\partial}{\partial t} s = - \frac{1}{\rho} \frac{\partial}{\partial z} F_H + D$$ (4.449)

where $s = c_p T + g z$ is the dry static energy, $z$ is the geopotential height above the local surface (does not include the surface elevation) and $D$ is the heating rate due to the dissipation of resolved kinetic energy in the diffusion process. The diffusive fluxes are defined as:

$$F_{u,v} = -\rho K_m \frac{\partial}{\partial z}(u,v),$$ (4.450)

$$F_{q,H} = -\rho K_{q,H} \frac{\partial}{\partial z}(q,s) +\rho K_{q,H}^t\gamma_{q,H} .$$ (4.451)

The viscosity $K_m$ and diffusivities $K_{q,H}$ are the sums of: turbulent components $K_{m,q,H}^t$, which dominate below the mesopause; and molecular components $K_{m,q,H}$, which dominate above  120 km. The turbulent diffusivities are the sum of two components, free atmosphere and boundary layer diffusivities, defined below. In the future, these terms also may include effective diffusivities from the gravity wave parameterization. The non-local transport terms $\gamma_{q,H}$ are given by the ABL parameterization. Note that $F_q$, as defined in (4.451) and implemented in CAM 3.0, does not include the term which causes diffusive separation of constituents of differing molecular weights. The molecular diffusion in CAM 3.0 is currently incomplete and should be used with caution. The molecular viscosity and diffusivities are all currently defined as $3.55\times 10^{-7} T^{2/3}/\rho$. A more complete form, allowing separation of constituents, will be implemented later.

The kinetic energy dissipation term $D$ in (4.449) is determined by forming the equation for total energy from (4.448-4.449):

$$\frac{\partial E}{\partial t} = u\frac{\partial u}{\partial t} + v\frac{\partial v}{\partial t} + \frac{\partial s}{\partial t}$$ (4.452)

$$= -\frac{1}{\rho} \left(u\frac{\partial F_u}{\partial z} + v\frac{\partial F_v}{\partial z} + \frac{\partial F_H}{\partial z} \right) + D$$ (4.453)

$$= -\frac{1}{\rho} \left(\frac{\partial F_{KE}}{\partial z} + \frac{\partial F_H}{\partial z}\right).$$ (4.454)

The diffusive kinetic energy flux in (4.454) is

$$F_{KE} \equiv u F_u + v F_v$$ (4.455)

and the kinetic energy dissipation is

$$D \equiv -\frac{1}{\rho} \left( F_u\frac{\partial u}{\partial z} + F_v\frac{\partial v}{\partial z} \right).$$ (4.456)

To show that $D$ is positive definite, we use (4.450) to expand for $F_u$ and $F_v$:

$$D = (K_m^t + K_m^m) \left[ \left(\frac{\partial u}{\partial z}\right)^2 +\left(\frac{\partial v}{\partial z}\right)^2 \right] \ge 0.$$ (4.457)

We show that energy is conserved in the column by integrating (4.454) in the vertical, from the surface ($z=0$), to the top of the model ( $z={z_{top}}$):

$$\int_0^{z_{top}} \rho\frac{\partial E}{\partial t} dz = (F_{KE} + F_H)\vert_{z_{top}}^0.$$ (4.458)

Therefore, the vertically integrated energy will only change because of the boundary fluxes of energy, of which only the surface heat flux, $F_H(z=0)$, is usually nonzero. It is typically assumed that the surface wind vanishes, even over oceans and sea ice, giving $F_{KE}(z=0)=0$. Then, the surface stress $F_{u,v}(z=0)$ does not change the total energy in the column, but does result in kinetic energy dissipation and heating near the surface (see below) For coupled models, nonzero surface velocities can be accommodated by including $F_{KE}$ on both sides of the surface interface.

4.11.1 Free atmosphere turbulent diffusivities

The free atmospheric turbulent diffusivities are typically taken as functions of length scales $\ell_c$ and local vertical gradients of wind and virtual potential temperature, e.g.

$$K_c = {{\ell_c}^2 S F_c(Ri)} .$$ (4.459)

Here $S$ is the local shear, defined by

$$S = \left\vert \frac{\partial {\boldsymbol {V}}}{\partial z} \right\vert,$$ (4.460)

and the mixing length $\ell_c$ is generally given by

$$\frac{1}{\ell_c} = \frac{1}{k z} + \frac{1}{\lambda_c},$$ (4.461)

where $k$ is the Von Karman constant, and $\lambda_c$ is the so-called asymptotic length scale, taken to be 30 m above the ABL. Since the lowest model level is always greater than 30 m in depth, $\ell_c$ is simply set to 30 m in CAM 3.0. Furthermore, $F_c(Ri)$ denotes a functional dependence of $K_c$ on the gradient Richardson number:

$$Ri = \frac{g}{\theta_v} \frac{\partial \theta_v / \partial z}{S^2},$$ (4.462)

where $\theta_v$ is the virtual potential temperature,

$$\theta_v = \theta \left[ 1 + \left( \frac{R_v}{R} -1 \right) q \right] .$$ (4.463)

For simplicity, in the free atmosphere, we specify the same stability functions $F_c$ for all $c$. For unstable conditions $(Ri < 0)$ we choose

$$F_c(Ri) = {\left( 1 - 18 Ri \right)^{1/2}} ,$$ (4.464)

$$(Ri > 0)$$

$$F_c(Ri) =\frac{1} {1 + 10 Ri\left(1 + 8 Ri \right)} ,$$ (4.465)

This means that no distinction is made between turbulent vertical diffusion of heat, scalars and momentum outside the boundary layer. However, separate coefficient arrays are maintained and other parameterizations (such as gravity wave drag) may provide distinct diffusivities. We also note the the turbulent diffusivity is the same for all constituents, even within the ABL. However, the molecular diffusivities differ for each constituent since they depend on it's molecular weight.

4.11.2 ``Non-local'' atmospheric boundary layer scheme

The free atmosphere turbulent diffusivities, described above, are an example of the local diffusion approach. In such an approach, the turbulent flux of a quantity is proportional to the local gradient of that quantity (e.g. (4.450)-(4.451)). In addition, the eddy diffusivity depends on local gradients of mean wind and mean virtual temperature (see (4.459)). These are reasonable assumptions when the length scale of the largest turbulent eddies is smaller than the size of the domain over which the turbulence extends. In the Atmospheric Boundary Layer (ABL) this is typically true for neutral and stable conditions only. For unstable and convective conditions, however, the largest transporting eddies may have a size similar to the boundary layer height itself, and the flux can be counter to the local gradient [76,48]. In such conditions a local diffusion approach is no longer appropriate, and the eddy diffusivity is better represented with turbulent properties characteristic of the ABL. We will refer to such an approach as non-local diffusion.

To account for ``non-local'' transport by convective turbulence in the ABL, the local diffusion term for constituent $c$ is modified as in (4.451):

$${\overline {w^\prime C^\prime}} = - {K_c \left( \frac{\partial C}{\partial z} - \gamma_ c\right)} ,$$ (4.466)

where $K_c$ is the non-local eddy diffusivity for the quantity of interest. The term $\gamma_c$ is a ``non-local'' transport term and reflects non-local transport due to dry convection. Eq. (4.466) applies to static energy, water vapor, and passive scalars. No countergradient term is applied to the wind components, so (4.450) does not contain these terms. For stable and neutral conditions the non-local term is not relevant for any of the quantities. The eddy diffusivity formalism is, however, modified for all conditions.

In the non-local diffusion scheme the eddy diffusivity is given by

$$K_c = k w_t z\left(1 - \frac{z}{h}\right)^2 ,$$ (4.467)

where $w_t$ is a turbulent velocity scale and $h$ is the boundary layer height. Equation (4.467) applies for heat, water vapor and passive scalars. The eddy diffusivity of momentum $K_m$, is also defined as (4.467) but with $w_t$ replaced by another velocity scale $w_m$. With proper formulation of $w_t$ (or $w_m$) and $h$, it can be shown that equation (4.467) behaves well from very stable to very unstable conditions in horizontally homogeneous and quasi-stationary conditions. For unstable conditions $w_t$ and $w_m$ are proportional to the so-called convective velocity scale $w_\ast$, while for neutral and stable conditions $w_t$ and $w_m$ are proportional to the friction velocity $u_\ast$.

The major advantage of the present approach over the local eddy diffusivity approach is that large eddy transport in the ABL is accounted for and entrainment effects are treated implicitly. Above the ABL, $\gamma_c=0$ so (4.466) reduces to a local form with $K_c$ given by (4.459). Near the top of the ABL we use the maximum of the values by (4.459) and (4.467), although (4.467) almost always gives the larger value in practice.

The non-local transport term in (4.466), $\gamma_c$, represents non-local influences on the mixing by turbulence [48]. As such, this term is small in stable conditions, and is therefore neglected under these conditions. For unstable conditions, however, most heat and moisture transport is achieved by turbulent eddies with sizes on the order of the depth $h$ of the ABL. In such cases, a formulation for $\gamma_c$ consistent with the eddy formulation of (4.466) is given by

$$\gamma_c = a \frac{w_{\ast} (\overline{w^\prime C^\prime})_s} {{w_m}^2 h} ,$$ (4.468)

where $a$ is a constant and $(\overline{w^\prime C^\prime})_s$ is the surface flux (in kinematic units) of the transported scalar. The form of (4.468) is similar to the one proposed in Holtslag and Moeng [76]. The non-local correction vanishes under neutral conditions, for which $w_{\ast} = 0$.

The formulations of the eddy-diffusivity and the non-local terms are dependent on the boundary layer height $h$. The CCM2 configuration of this non-local scheme made use of a traditional approach to estimating the boundary layer depth by assuming a constant value for the bulk Richardson number across the boundary layer depth so that $h$ was iteratively determined using

$$h = \frac{Ri_{cr} \left\{u(h)^2 + v(h)^2\right\}}{(g / \theta_s) \left(\theta_v (h) - \theta_s\right)} ,$$ (4.469)

where $Ri_{cr}$ is a critical bulk Richardson number for the ABL; $u(h)$ and $v(h)$ are the horizontal velocity components at $h$; $g/\theta_s$ is the buoyancy parameter and $\theta_v (h)$ is the virtual temperature at $h$. The quantity $\theta_s$ is a measure of the surface air temperature, which under unstable conditions was given by

$$\theta_s = \theta_v \left(z_s\right) + b \frac{(\overline{w^\prime \theta^\prime_v})_s}{w_m} ,$$ (4.470)

where $b$ is a constant, $(\overline{w^\prime \theta^\prime_v})_s$ is the virtual heat flux at the surface, $\theta_v\left(z_s\right)$ is a virtual temperature in the atmospheric surface layer (nominally 10 m), $b {(\overline{w^\prime \theta^\prime_v})_s / {w_m}}$ represents a temperature excess (a measure of the strength of convective thermals in the lower part of the ABL) and unstable conditions are determined by $(\overline{w^\prime \theta^\prime_v})_s >0$. The quantity $\theta_v\left(z_s\right)$ was calculated from the temperature and moisture of the first model level and of the surface by applying the procedure in Geleyn [60]. The value of the critical bulk Richardson number $Ri_{cr}$ in (4.469), which generally depends on the vertical resolution of the model, was chosen as $Ri_{cr}$ = 0.5 for the CCM2.

Vogelezang and Holtslag [178] have recently studied the suitability of this formulation in the context of field observations, large-eddy simulations [128], and an $E-\epsilon$ turbulence closure model [53]. They propose a revised formulation which combines shear production in the outer region of the boundary layer with surface friction, where the Richardson number estimate is based on the differences in wind and virtual temperature between the top of the ABL and a lower height that is well outside the surface layer (i.e. 20 m - 80 m). In addition to providing more realistic estimates of boundary layer depth, the revised formulation provides a smoother transition between stable and neutral boundary layers. Consequently, CAM 3.0 employs the Vogelezang and Holtslag [178] formulation for estimating the atmospheric boundary layer height, which can be written as

$$h = z_s + \frac{Ri_{cr} \left\{ (u(h) - u_{SL})^2 + (v(h) - v_{SL... ...}u^2_\ast \right\}}{(g / \theta_{SL}) \left(\theta_v(h) - \theta_{SL}\right)} .$$ (4.471)

The quantities $u_{SL}$, $v_{SL}$, and $\theta_{SL}$ represent the horizontal wind components and virtual potential temperature just above the surface layer (nominally 0.1$h$). In practice, the lowest model level values for these quantities are used to iteratively determine $h$ for all stability conditions, where the critical Richardson number, $Ri_{cr}$, is assumed to be 0.3. The disposable parameter ${\cal B}$ has been experimentally determined to be equal to 100 (see Vogelezang and Holtslag [178]). The computation starts by calculating the bulk Richardson number $Ri$ between the level of $\theta_{SL}$ and subsequent higher levels of the model. Once $Ri$ exceeds the critical value, the value of $h$ is derived by linear interpolation between the level with $Ri > Ri_{cr}$ and the level below.

Using the calculated value for $h$ and the surface fluxes, we calculate the velocity scales, the eddy diffusivities with (4.467), and the countergradient terms with (4.468), for each of the transported constituents. Subsequently, the new profiles for $\theta$, $q$, $u$, and $v$ are calculated using an implicit diffusion formulation.

The turbulent velocity scale of (4.467) depends primarily on the relative height $z/h$ ($h$ is boundary layer height), and the stability within the ABL. Here stability is defined with respect to the surface virtual heat flux $(\overline{w^\prime \theta^\prime_v})_s$. Secondly, the velocity scales are also generally dependent on the specific quantity of interest. We will assume that the velocity scales for mixing of passive scalars and specific humidity are equal to the one for heat, denoted by $w_t$. For the wind components, the velocity scale is different and denoted by $w_m$. The specification of $w_t$ and $w_m$ is given in detail by Troen and Mahrt [174]. Holtslag et al. [75] have rewritten the velocity scale, in terms of the more widely accepted profile functions of Dyer [54], and have given a new formulation for very stable conditions. Below we follow the latter approach.

For stable ( $(\overline{w^\prime \theta^\prime_v})_s < 0$) and neutral surface conditions ( $(\overline{w^\prime \theta^\prime_v})_s = 0$), the velocity scale for scalar transport is

$$w_t = \frac{u_\ast}{\phi_h} ,$$ (4.472)

where $u_\ast$ is the friction velocity defined by

$$u_\ast = \left[ (\overline{u^\prime w^\prime})_s^2 + (\overline{v^\prime w^\prime})_s^2 \right]^{1/4}.$$ (4.473)

Furthermore, $\phi_h$ is the dimensionless vertical temperature gradient given by Dyer [54],

$$\phi_h = 1 + 5 \frac{z}{L} ,$$ (4.474)

for $0 \leq z/L \leq 1$. Here $L$ is the Obukhov length, defined by

$$L = \frac{- u_{\ast}^3}{k (g/\theta_{v0}) (\overline{w^\prime\theta^\prime_v})_0}.$$ (4.475)

For $z/L > 1$,

$$\phi_h = 5 + \frac{z}{L} ,$$ (4.476)

which matches (4.474) for $z/L = 1$. Equation (4.476) is a simple means to prevent $\phi_h$ from becoming too large (and $K_c$ too small) in very stable conditions. In stable conditions, the exchange coefficients for heat and momentum are often found to be similar. Therefore we may use $w_m$=$w_t$.

For unstable conditions $(\overline{w^\prime \theta^\prime_v})_s >0$, we have that $w_t$ and $w_m$ differ in the surface layer $(z/h \leq 0.1)$ and in the outer layer of the ABL $(z/h > 0.1)$. For the surface layer, $w_t$ is given by (4.472) with

$$\phi_h = \left(1 - 15 \frac{z}{L}\right)^{-1/2} .$$ (4.477)

Similarly, $w_m$ is written as

$$w_m = \frac{u_\ast}{\phi_m} ,$$ (4.478)

where $\phi_m$ is the dimensionless wind gradient given by

$$\phi_m = \left(1 - 15 \frac{z}{L}\right)^{-1/3} .$$ (4.479)

In the surface layer, the scalar flux is normally given by

$$(\overline{w^\prime c^\prime})_0 = - \frac{k u_{\ast} z}{\phi_h} \left( \frac{\partial C}{\partial z} \right) .$$ (4.480)

Comparison with (4.466) and (4.467) shows that, in the surface layer, we should have $a = 0$ in (4.468) for consistency.

For the outer layer, $w_t$ and $w_m$ are given by

$$w_t = w_m/Pr ,$$ (4.481)
where

$$w_m = \left(u_\ast^3 + c_1 w_\ast^3\right)^{1/3} ,$$ (4.482)
and

$$w_\ast = \left(\left(g/\theta_{v0}\right) (\overline{w^\prime\theta^\prime_v})_0 h\right)^{1/3}$$ (4.483)

is the convective velocity scale. Furthermore, $Pr$ is the turbulent Prandtl number and $c_1$ is a constant. The latter is obtained by evaluating the dimensionless vertical wind gradient $\phi_m$ by (4.479) at the top of the surface layer, as discussed by Troen and Mahrt [174]. This results in $c_1 = 0.6$. For very unstable conditions ($h \gg -L$ or $w_{\ast}/u_{\ast} \gg 0)$, it can be shown with (4.481) that $w_m$ is proportional to 0.85 $w_\ast$, while for the neutral case $w_m = u_\ast$. The turbulent Prandtl number $Pr$ $(= K_m/K_h = w_m/w_t)$ of (4.481) is evaluated from

$$Pr = \frac{\phi_h}{\phi_m} \left( \frac{z}{L}\right) + a k \frac{z}{h} \frac{w_\ast}{w_m}$$ (4.484)

for $z = 0.1 h$. Equation (4.484) arises from matching (4.466), (4.467), (4.468), and (4.480) at the top of the surface layer. As in Troen and Mahrt we assume that $Pr$ is independent of height in the unstable outer layer. Its value decreases from $Pr = 1$ for the neutral case ($z/L =0$ and $w_\ast = 0$), to $Pr = 0.6$ for $w_\ast/u_\ast \simeq 10$ in very unstable conditions.

In very unstable conditions, the countergradient term of (4.468) approaches

$$\gamma_c = d \frac{\overline{wC}_0}{w_{\ast} h} ,$$ (4.485)

where $d \simeq a/0.85^2$, because for very unstable conditions we obtain $w_m \simeq 0.85 w_{\ast}$. Since typically $d \simeq 10$ Troen and Mahrt [174], we have $a = 7.2$. Similarly, the temperature excess of (4.470) reads in this limit as $d (\overline{w^\prime\theta^\prime_v})_0/w_{\ast}$. This leads to $b$ (= 0.85 $d$) = 8.5 in (4.470).

Finally, using the velocity scales described above, the flux equation (4.466) is continuous in relative height ($z/h$) and in the boundary layer stability parameter ($h/L$ or $w_{\ast}/u_{\ast}$).

4.11.3 Discretization of the vertical diffusion equations

In CAM 3.0, as in previous version of the CCM, (4.448-4.451) are cast in pressure coordinates, using

$$dp = -\rho g dz,$$ (4.486)

and discretized in a time-split form using an Euler backward time step. Before describing the numerical solution of the diffusion equations, we define a compact notation for the discrete equations. For an arbitrary variable $\psi$, let a subscript denote a discrete time level, with current step $\psi_n$ and next step $\psi_{n+1}$. The model has $L$ layers in the vertical, with indexes running from top to bottom. Let $\psi^k$ denote a layer midpoint quantity and let $\psi^{k-}$ denote the value on the upper interface of layer $k$ while $\psi^{k+}$ denotes the value on the lower interface. The relevant quantities, used below, are then:

$$\psi^{k+} = (\psi^{k}+\psi^{k+1})/2,\quad k\in (1,2,3,...,L-1)$$

$$\psi^{k-} = (\psi^{k-1}+\psi^{k})/2,\quad k\in (2,3,4...,L)$$

$$\delta^{k}\psi = \psi^{k+}-\psi^{k-},$$

$$\delta^{k+}\psi = \psi^{k+1}-\psi^{k},$$

$$\delta^{k-}\psi = \psi^{k}-\psi^{k-1},$$ (4.487)

$$\psi_{n+} = (\psi_{n}+\psi_{n+1})/2,$$

$$\delta_n\psi = \psi_{n+1}-\psi_{n},$$

$$\delta t = t_{n+1}-t_{n},$$

$$\Delta^{k,l} = 1, k=l,$$

$$= 0, k\neq l.$$

Like the continuous equations, the discrete equations are required to conserve momentum, total energy and constituents. The discrete forms of (4.448-4.449) are:

$$\frac{\delta_n (u,v,q)^k}{\delta t} = g \frac{\delta^k F_{u,v,q}}{\delta^k p}$$ (4.488)

$$\frac{\delta_n s^k}{\delta t} = g \frac{\delta^k F_{H}}{\delta^k p} + D^k.$$ (4.489)

For interior interfaces, $1\le k \le L-1$,

$$F_{u,v}^{k+} = \left(g\rho^2 K_m\right)_n^{k+} \frac{\delta^{k+} (u,v)_{n+1}} {\delta^{k+} p}$$ (4.490)

$$F_{q,H}^{k+} = \left(g\rho^2 K_{q,H}\right)_n^{k+} \frac{\delta^{k+} (u,v)_{n+1}} {\delta^{k+} p}$$

$$+ \left(\rho K_{q,H}^t \gamma_{q,H}\right)_n^{k+}.$$ (4.491)

Surface fluxes $F_{u,v,q,H}^{L+}$ are provided explicitly at time $n$ by separate surface models for land, ocean, and sea ice while the top boundary fluxes are usually $F_{u,v,q,H}^{1-}=0$. The turbulent diffusion coefficients $K_{m,q,H}^{t}$ and non-local transport terms $\gamma_{q,H}$ are calculated for time $n$ by the turbulence model described above, which is identical to CCM3. The molecular diffusion coefficients, described earlier, are only included if the model top is above $\sim 90$ km, in which case nonzero top boundary fluxes may be included for heat and some constituents.

The free atmosphere turbulent diffusivities ${K}_{n}^{k+}$, given by (4.459-4.465), are discretized as

$${K}_{n}^{k+} = {K}_N^{k+} \cdot F_c\! \left(R_{I}^{k+}\right) \ge 0.01.$$ (4.492)

The stability function is:

$$F_c(R_I) = \begin{cases}1/\left(1 + 10 R_I[1+8 R_I]\right) & \tex... ...1ex] \sqrt{1 - 18 R_I} & \text{for} R_I < 0 \text{(unstable)} , \end{cases}$$ (4.493)

The neutral $K_N$ is calculated by

$${K}_N^{k+} = \ell^2 \frac{\left[\left( \delta^{k+}{u}_{n} \right)^2 + \left( \delta^{k+}{v}_{n} \right)^2 \right]^{1/2}}{\delta^{k+} z_{n}},$$ (4.494)

with $\ell=30$ m. The Richardson number in the free atmosphere is calculated from

$$R_{I}^{k+} = \frac{g} {\theta_{v}^{k+}} \times \frac{ \delta^{k+} z_n \delta... ...ta_{v}} {\left(\delta^{k+} u_{n} \right)^2 + \left(\delta^{k+} v_{n} \right)^2}$$ (4.495)
where

$$\theta^{k}_{v} = \theta_{n}^{k} \left( 1.0 + \left( \frac{R_v}{R} -1 \right) q_{n}^{k} \right) .$$ (4.496)

Similarly to the continuous form (4.456), $D^k$ is determined by separating the kinetic energy change over a time step into the kinetic energy flux divergence and the kinetic energy dissipation. The discrete system is required to conserve energy exactly:

$${\sum_{k=1}^L \left[(u_{n+1}^k)^2 + (v_{n+1}^k)^2 + s_{n+1}^k\right]\delta^k p =}$$

$$\sum_{k=1}^L \left[(u_{n}^k)^2 + (v_{n}^k)^2 + s_{n}^k\right]\delta^k p + \delta t (F_H^{L+} + F_H^{1-}),$$

where we have assumed zero boundary fluxes for kinetic energy. This leads to

$$D^k = \frac{g}{2\delta^k p}(d_u^{k+} + d_u^{k-} + d_v^{k+} + d_v^{k-})$$ (4.498)

$$d_{u,v}^{k+} = \delta^{k+}(u,v)_{n+} F_{u,v}^{k+},\quad 1\le k\le L-1$$ (4.499)

$$d_{u,v}^{L+} = -2 (u,v)_{n+}^L F_{u,v}^{L+}$$ (4.500)

According to (4.498), the internal dissipation of kinetic energy in each layer $D^k$ is the average of the dissipation on the bounding interfaces $d_{u,v}^{k\pm}$, given by (4.499) and (4.500). Expanding (4.499) using (4.490) and recalling that $u_{n+} = (u_{n+1} + u_n)/2$,

$$d_u^{k+} = \frac{\left(g\rho^2 K_m\right)^{k+}}{2\delta^{k+}p} \left[\left(\delta^{k+}u_{n+1}\right)^2 + \delta^{k+}u_{n+1}\delta^{k+}u_{n}\right],$$ (4.501)

for $1\le k \le L-1$ and similarly for $d_v^{k+}$. The discrete kinetic energy dissipation is not positive definite, because the last term in (4.501) is the product of the vertical difference of momentum at two time levels. Although $d_{u,v}^{k+}$ will almost always be $>0$, values $\le 0$ may occur occasionally. The kinetic energy dissipation at the surface is

$$d_{u,v}^{L+} = - \left[(u,v)_{n+1}^L + (u,v)_n^L\right]F_{u,v}^{L+}.$$ (4.502)

Since the surface stress is opposed to the bottom level wind, the surface layer is heated by the frictional dissipation. However, $d_{u,v}^{L+}$ is not guaranteed to be positive, since it involves the bottom level wind at two time levels.

Note that it has been assumed that the pressure does not change within the vertical diffusion, even though there are boundary fluxes of constituents, including water. This assumption has been made in all versions of the CCM and is still made in CAM 3.0. This assumption will be removed in a future version of CAM 3.0, since the implied horizontal fluxes of dry air, to compensate for the boundary flux of water, cause implied fluxes of other constituents.

4.11.4 Solution of the vertical diffusion equations

A series of time-split operators is actually defined by (4.488-4.491) and (4.498-4.500). Once the diffusivities ($K_{m,q,H}$) and the non-local transport terms ( $\gamma_{q,H}$) have been determined, the solution of (4.488-4.491), proceeds in several steps.

1. update the $q$ and $s$ profiles using $\gamma_{q,H};$
2. update the bottom level values of $u$, $v$, $q$ and $s$ using the surface fluxes;
3. invert (4.488) and (4.490) for $u,v_{n+1}$;
4. compute $D$ and use to update the $s$ profile;
5. invert (4.488,4.489) and (4.491) for $s_{n+1}$ and $q_{n+1}$;

Note that since all parameterizations in CAM 3.0 return tendencies rather than modified profiles, the actual quantities returned by the vertical diffusion are $\delta_n (u,v,s,q) / \delta t$.

The non-local transport terms, $\gamma_{q,H}$, given by (4.468), cannot be treated implicitly because they depend on the surface flux, the boundary layer depth and the velocity scale, but not explicitly on the profile of the transported quantity. Therefore, application of $\gamma_q$ is not guaranteed to give a positive value for $q$ and negative values may not be removed by the subsequent implicit diffusion step. This problem is not strictly numerical; it arises under highly non-stationary conditions for which the ABL formulation is not strictly applicable. In practice, we evaluate

$$q_{n*} = q_n + \frac{g\delta t}{\delta^k p} \delta^k\left[\rho K_{q}^t \gamma_{q}\right]_n$$ (4.503)

and check the $q_{n*}$ profile for negative values (actually for $q_{n*}^k < q_{min}$, where $q_{min}$ may be $>0$). If any negative values are found, we set $q_{n*} = q_n$ for that constituent profile (but not for other constituents at the same point).

Equations (4.488-4.491) constitute a set of four tridiagonal systems of the form

$$-A^k \psi^{k+1}_{n+1} + B^k\psi^k_{n+1} -C^k\psi^{k-1}_{n+1} = \psi^k_{n\prime},$$ (4.504)

where $\psi_{n\prime}$ indicates $u$, $v$, $q$, or $s$ after updating from time $n$ values with the nonlocal and boundary fluxes. The super-diagonal ($A^k$), diagonal ($B^k$) and sub-diagonal ($C^k$) elements of (4.504) are:

$$A^k = \frac{1}{\delta^{k} p} \frac{\delta t}{\delta^{k+}p}\left(g^2\rho^2 K\right)_n^{k+},$$ (4.505)

$$B^k = 1 + A^k + C^k,$$ (4.506)

$$C^k = \frac{1}{\delta^{k} p} \frac{\delta t}{\delta^{k-}p}\left(g^2\rho^2 K\right)_n^{k-}.$$ (4.507)

The solution of (4.504) has the form

$$\psi^k_{n+1} = E^k \psi^{k-1}_{n+1} + F^k,$$ (4.508)

or,

$$\psi^{k+1}_{n+1} = E^{k+1} \psi^{k}_{n+1} + F^{k+1}.$$ (4.509)

Substituting (4.509) into (4.504),

$$\psi^{k}_{n+1} = \frac{C^k}{B^k - A^k E^{k+1}} \psi^{k-1}_{n+1} + \frac{\psi^k_{n\prime} + A^k F^{k+1}}{B^k - A^k E^{k+1}}.$$ (4.510)

Comparing (4.508) and (4.510), we find

$$E^k = \frac{C^k} {B^k - A^k E^{k+1}}, \quad L>k>1,$$ (4.511)

$$F^k = \frac{\psi^k_{n\prime} + A^k F^{k+1}}{B^k - A^k E^{k+1}}, \quad L>k>1.$$ (4.512)

The terms $E^k$ and $F^k$ can be determined upward from $k=L$, using the boundary conditions

$$E^{L+1} = F^{L+1} = A^L = 0.$$ (4.513)

Finally, (4.510) can be solved downward for $\psi^{k}_{n+1}$, using the boundary condition

$$C^1 = 0 \Rightarrow E^1 = 0.$$ (4.514)

CCM1-3 used the same solution method, but with the order of the solution reversed, which merely requires writing (4.509) for $\psi^{k-1}_{n+1}$ instead of $\psi^{k+1}_{n+1}$. The order used here is particularly convenient because the turbulent diffusivities for heat and all constituents are the same but their molecular diffusivities are not. Since the terms in (4.511-4.512) are determined from the bottom upward, it is only necessary to recalculate $A^k$, $C^k$, $E^k$ and $1/({B^k - A^k E^{k+1}})$ for each constituent within the region where molecular diffusion is important. Note that including the diffusive separation term for constituents (which will be in the next version of CAM 3.0) adds additional terms to the definitions of $A^k$, $B^k$, and $C_k$, but does not otherwise change the solution method.

4.11.5 Discrete equations for $s$, $T$, and $z$

The dry static energy at step $n$ and level $k$ is

$$s_{n}^k = c_p^d T_{n}^k + g z^k,$$ (4.515)

which can be calculated from $T_{n}$ by integrating the hydrostatic equation using the perfect gas law.

$$g z \equiv \Phi = \Phi_s + \int_{p_s}^p R T d\ln p^\prime,$$ (4.516)

where $\Phi$ is the geopotential, $\Phi_s$ is the geopotential at the Earth's surface and $p_s$ is the surface pressure. A fairly arbitrary discretization of (4.516) can be represented using a triangular hydrostatic matrix $H^{kl}$,

$$\Phi^k = \Phi_s + \sum_{l=L}^{k} R^l H^{kl} T^l.$$ (4.517)

Note that (4.517) is often written in terms of the virtual temperature $T_v = T R/R^d$. The apparent gas constant $R$ includes the effect of water vapor and is defined as

$$R = R^d + (R^w - R^d) q,$$ (4.518)

where $R^d$ is the apparent gas constant for dry air and $R^w$ is the gas constant for water vapor.

Using (4.517) in (4.515), we have

$$s_{n}^k = c_p^d T_{n}^k + \sum_{l=L}^{k} R^l H^{kl} T_{n}^l,$$ (4.519)

$$= \left(c_p^d + R^k H^{kk} \right) T_{n}^k + \Phi_n^{k+}.$$ (4.520)

The interface geopotential in (4.520) is defined as

$$\Phi^{k+} = \sum_{l=L}^{k+1} R^k H^{kl} T^l,$$ (4.521)

and $R^k$ is evaluated from (4.518), using $q^k_{n}$. Although the correct boundary condition on (4.521) is $\Phi_{L+}=\Phi_s$, within the parameterization suite it is usually sufficient to take $\Phi_{L+}=0$.

The definition of the hydrostatic matrix $H$ depends on the numerical method used in the dynamics and is subject to constraints from energy and mass conservation. The definitions of $H$ for the three dynamical methods used in CAM 3.0 are given in the dynamics descriptions.

After $s_n$ is modified by diabatic heating in a time split process, the new $s_{n+1}=s_n+Q_n\delta t$ can be converted into $T_{n+1}$ and $\Phi_{n+1}$ using (4.520):

$$s_{n+1}^k = \left(c_p^d + R^k H^{kk} \right) T_{n+1}^k + \Phi_{n+1}^{k+}$$ (4.522)

$$T_{n+1}^k = \left(s_{n+1}^k - \Phi_{n+1}^{k+}\right) \left(c_p + R^k H^{kk}\right)^{-1}$$ (4.523)

with $R^k$ evaluated from using $q_{n+1}^k$. Once $H$ is defined, (4.521) and (4.523) can be solved for $T_{n+1}$ and $\Phi_{n+1}$ from the bottom up. Since the latter must normally recalculated if $T$ is modified, calculating $T$ and $\Phi$ from $s$ involves the same amount of computation as calculating $\Phi$ and $s$ from $T$.