- 3.1.1 Generalized terrain-following vertical coordinates
- 3.1.2 Conversion to final form
- 3.1.3 Continuous equations using
- 3.1.4 Semi-implicit formulation
- 3.1.5 Energy conservation
- 3.1.6 Horizontal diffusion
- 3.1.7 Finite difference equations
- 3.1.8 Time filter
- 3.1.9 Spectral transform
- 3.1.10 Spectral algorithm overview
- 3.1.11 Combination of terms
- 3.1.12 Transformation to spectral space
- 3.1.13 Solution of semi-implicit equations
- 3.1.14 Horizontal diffusion
- 3.1.15 Initial divergence damping
- 3.1.16 Transformation from spectral to physical space
- 3.1.17 Horizontal diffusion correction
- 3.1.18 Semi-Lagrangian Tracer Transport
- 3.1.19 Mass fixers
- 3.1.20 Energy Fixer
- 3.1.21 Statistics Calculations
- 3.1.22 Reduced grid

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 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.

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 , then we require the generalized coordinate to satisfy:

- is a monotonic function of .
- where is the top of the model.

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:

The notation follows standard conventions, and the following terms have been introduced with :

The terms , and represent the sources and sinks from the parameterizations for momentum (in terms of and ), temperature, and moisture, respectively. The terms and represent sources due to horizontal diffusion of momentum, while and 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:

Note that the bounds on the vertical integrals are specified as values of (

Equations (3.1)-(3.16) are the complete set which must be solved by a GCM. However, in order to solve them, the function must be specified. In advance of actually specifying , 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 must be converted to terms involving only and horizontal derivatives of . The former can be evaluated once the function is specified.

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

since is given by (3.15). Similarly, the first term on the right-hand side of (3.15) can be expanded as

and (3.7) invoked to specify .

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

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:

The second term in (3.20) vanishes because , while the first term is easily treated once is specified. Substituting (3.20) into (3.19), one obtains:

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

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

Once is specified, then 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 explicitly, since (3.23)-(3.32) only requires that and be determined. It is sufficient to specify and to let be defined implicitly. This will be done in section 3.1.7. In the case that and , (3.23)-(3.32) can be reduced to the set of equations solved by CCM1.

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 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 as a prognostic variable, to equations using . 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 , all terms involving will then appear as .

Equations (3.23)-(3.32) become:

The above equations reduce to the standard equations used in CCM1 if and . (Note that in this case .)

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:

where are the remaining nonlinear terms not explicitly written in (3.43)-(3.45). The terms involving and may be expanded into vertical integrals using (3.40) and (3.42), while the term can be converted to , giving:

Once again, only terms that will be linearized have been explicitly represented in (3.46)-(3.48), and the remaining terms are included in , , and . Anticipating the linearization, and have been replaced by and 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:

In the special case that , (3.46)-(3.48) can be converted into equations involving only instead of , and (3.50) and (3.51) are not required. This is a major difference between the hybrid coordinate scheme being developed here and the 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:

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

In the hybrid coordinate described below, is a linear function of , so 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 layer midpoints denoted by an integer index, (see Figure 3.1). The interface between and is denoted by a half-integer index, . The model top is at , and the Earth's surface is at . It is further assumed that vertical integrals may be written as a matrix (of order ) times a column vector representing the values of a field at the 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 will denote the vector transpose.

The finite difference forms of (3.52)-(3.54) may then be written down as:where denotes a time varying value at time step . The quantities and are defined so as to complete the right-hand sides of (3.43)-(3.45). The components of are given by . This definition of the vertical difference operator will be used in subsequent equations. The reference matrices and , and the reference column vectors and , depend on the precise specification of the vertical coordinate and will be defined later.

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
and
contributions, can be combined with the continuity
equation

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

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
(
) and the vertical derivative of a
field (
). 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:

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
, and substituting for the geopotential
using (3.40), (3.60) can be written as:

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

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:

We now turn to the internal energy equation, obtained by combining the thermodynamic equation (3.36), without the , , and terms, and the continuity equation (3.59):

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:

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 in the nonlinear terms of the vorticity and divergence equations (3.38) and (3.39), and in the 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:

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

where for . The quantity 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

where for , and is included as an approximation to for and the symbol is similarly defined as in (3.62). will be determined so that 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

where we have used the relation

(3.73) |

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

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

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 and satisfy (3.75). We are, of course, still neglecting lack of conservation due to the truncation of the horizontal spherical harmonic expansions.

CAM 3.0 contains a horizontal diffusion term for , and 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 form on surfaces in the top three levels of the model and a linear form with a partial correction to pressure surfaces for temperature elsewhere. The 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 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 form of the horizontal diffusion is given by

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 rotations [24,133]. The form of the horizontal diffusion is applied

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 ( ) can be included using the relations:

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

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

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

The second term in consists of the leading term in the transformation of the operator to pressure surfaces. It is included to offset partially a spurious diffusion of over mountains. As with the form, the operator can be conveniently calculated in spectral space. The correction term is then completed after transformation of and back to grid-point space. As with the 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

for and , and

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

for and , and

for , where in the standard model only takes the value 2 in (3.92). The first step from to includes the transformation to spectral coefficients. The second step from to for and , or to for , is done on the spectral coefficients, and the final step from to for 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 anywhere. The quantities that must be known are and at the grid points. Therefore the coordinate is defined implicitly through the relation:

which gives

A set of levels may be specified by specifying and , such that , 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:

where notation such as denotes a column vector with components . In order to complete the system, it remains to specify the reference vector , together with the term , which results from the pressure gradient terms and also appears in the semi-implicit reference vector :

The matrices and (

(3.113) |

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 and 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.

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

where is the highest Fourier wavenumber included in the east-west representation, and is the highest degree of the associated Legendre polynomials for longitudinal wavenumber . 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: , and , where is defined above, is the highest degree of the associated Legendre polynomials, and is the highest degree of the Legendre polynomials for . The common truncations are subsets of this pentagonal case:

The quantity in (3.114) represents an arbitrary limit on the two-dimensional wavenumber , and for the pentagonal truncation described above is simply given by

.

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

With this normalization, the Coriolis parameter is

which is required for the absolute vorticity.

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

The inner integral represents a Fourier transform,

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

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

The weights themselves satisfy

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 must satisfy

For the common truncations, these become

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

The actual values of and 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 at the point
next to the south pole to 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 and
for
and
for
; then the summation
in (3.120) can be rewritten as

The symmetric (even) and antisymmetric (odd) components of are defined by

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

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 (consistent with the spectral representation) and producing a forecast of the grid-point values at time (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 (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.

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

where the explicit forms of the vectors and are given as

The divergence equation (3.96) is

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

The surface-pressure tendency (3.98) is

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

The terms of (3.136) are

The nonlinear term in (3.137) is

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 , the longitudinally differentiated term , and the meridionally differentiated term . 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),

where is the Fourier coefficient of with wavenumber at the Gaussian grid line . The longitudinally differentiated term is handled by integration by parts, using the cyclic boundary conditions,

(3.147) | ||

(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):

where is the Fourier coefficient of with wavenumber at the Gaussian grid line .

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

(3.151) | ||

(3.152) |

Defining the derivative of the associated Legendre polynomial by

(3.155) can be written

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

to each spherical harmonic function individually so that

where is the Fourier coefficient of the original grid variable .

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

where denotes a spherical harmonic coefficient of , and the form of , as a summation over the Gaussian grid, is given as

The spectral form of the divergence equation (3.135) becomes

where and are spectral coefficients of , and . The Laplacian of the total temperature in (3.135) is replaced by the equivalent Laplacian of the perturbation temperature in (3.159). is given by

The spectral thermodynamic equation is

with defined as

while the surface pressure equation is

where is given by

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 :

which is simply a set of simultaneous equations for the coefficients with given wavenumbers () at each level and is solved by inverting . In order to prevent the accumulation of round-off error in the global mean divergence (which if exactly zero initially, should remain exactly zero) is set to the null matrix rather than the identity, and the formal application of (3.165) then always guarantees . Once is known, and can be computed from (3.161) and (3.163), respectively, and all prognostic variables are known at time as spherical harmonic coefficients. Note that the mean component is not necessarily zero since the perturbations are taken with respect to a specified .

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 and equations have a similar form, so we write only the equation:

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

and are both set to 1 for = 0. The quantity represents the ``Courant number limiter'', normally set to 1. However, 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 in that level for all , where . This condition is applied whenever the wind speed is large enough that , 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 is not complete at this stage. In order to make the partial correction from to in (3.82) local, it is not included until grid-point values are available. This requires that also be transformed from spectral to grid-point space. The values of the coefficients and for the standard T42 resolution are msec and msec, respectively.

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 and linearly decreases to 0 over a specified number of days, , usually set to be 2. The damping is computed implicitly via time splitting after the horizontal diffusion.

After the prognostic variables are completed at time in spectral space , , , they are transformed to grid space. For a variable , the transformation is given by

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

In addition, the derivatives of are needed on the grid for the terms involving and ,

These required derivatives are given by

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

in which the only nonzero is and

Thus, the direct transformation is

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

using or 2 as appropriate for the or forms. These coefficients are transformed to grid space following (3.114) for the term and (3.186) and (3.187) for vorticity and divergence. Thus, the vorticity and divergence diffusion tendencies are converted to equivalent and diffusion tendencies.

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 diffusion from to surfaces is applied to . The frictional heating rate is calculated from the kinetic energy tendency produced by the momentum diffusion

where , and are the momentum equivalent diffusion tendencies, determined from and just as and are determined from and , and

These heating rates are then combined with the correction,

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

(3.194) | ||

(3.195) | ||

(3.196) |

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

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

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

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

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

Equation (3.202) represents the horizontal interpolation of at the departure point calculated assuming . Equation (3.203) represents the vertical interpolation of at the departure point, assuming .

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

where subscript denotes the arrival (Gaussian grid) point and subscript the midpoint of the trajectory. The velocity components at 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

where the subscript 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 latitude. Poleward of 70 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
and
, by an
angle so the equator goes through
. The longitude of the transformed system is chosen
to be zero at the arrival point. If the local geodesic system is
denoted by
, with
velocities
, the two systems are
related by

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 by design and and . The interpolations are always done in global spherical coordinates.

The interpolants are most easily defined on the interval 0 . Define

where is either or and the departure point falls within the interval . Following (23) of [146] with the Hermite cubic interpolant is given by

where is the value at the grid point , is the derivative estimate given below, and .

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

where

where can represent either or , or their counterparts in the geodesic coordinate system.

The derivative approximations used in (3.214) for are obtained by differentiating (3.215) with respect to , replacing by and evaluating the result at equal and . With these derivative estimates, the Hermite cubic interpolant (3.214) is equivalent to the Lagrangian (3.215). If we denote the four point stencil by the cubic derivative estimates are

The two dimensional interpolant is obtained as a tensor product application of the one-dimensional interpolants, with interpolations done first. Assume the departure point falls in the grid box and . Four interpolations are performed to find values at , , , and . This is followed by one interpolation in using these four values to obtain the value at . 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 and to wavenumber 1 components for and . The row across the pole is filled with the values from the first row below the pole shifted in longitude for and minus the value shifted by in longitude for and .

Once the departure point is known, the constituent value of 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 continuity (SCMO) described below. Define by

First, if then

Then, if either

or

is violated, or is brought to the appropriate bound of the relationship. These conditions ensure that the Hermite cubic interpolant is monotonic in the interval .

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

with . Note, this is the only place that the model actually requires an explicit specification of . The mid-point of the vertical trajectory is found by iteration

Note, the arrival point is a mid-level point where is carried, while the used for the interpolation to mid-points is at interfaces. We restrict by

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

with the restriction

The appropriate values of and are determined by interpolation (3.214), with the derivative estimates given by (3.215) and (3.216) for to . At the top and bottom we assume a zero derivative (which is consistent with (3.231) and (3.233)), for the interval , and for the interval . The estimate at the interior end of the first and last grid intervals is determined from an uncentered cubic approximation; that is at the interval is equal to from the interval, and at the interval is equal to at the interval. The monotonic conditions (3.227) to (3.228) are applied to the derivative estimates.

3.1.19 Mass fixers

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

Let , and 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.)

, and are the values after fixers are applied at the end of the time step.

, and are the values after the parameterizations have updated the moisture field and tracers.

Since the physics parameterizations do not change the surface pressure, and 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

where is the dry mass of the atmosphere. From the definition of the vertical coordinate,

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

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

In (3.237) and (3.238) the 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 factor is included to make the changes approximately proportional to mass per unit volume [141]. Satisfying (3.234) and (3.235) gives

and

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 denote the mixing ratio of constituents. Historically we have used the following relationship for conservation:

The term 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)

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

where the following shorthand notation is adopted:

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
, but the surface pressure is not
allowed to change. Since
, 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 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

There is a corresponding change in the first term of the numerator of
(3.243) in which is replace by . 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

(3.246) | ||

(3.247) |

(3.248) |

where is the net source of energy from the parameterizations. is the net downward solar flux at the model top, is the net upward longwave flux at the model top, is the net downward solar flux at the surface, is the net upward longwave flux at the surface, is the surface sensible heat flux, and is the total precipitation during the time step. From equation (3.237)

(3.249) |

and from (3.236)

(3.250) |

The energy fixer is chosen to have the form

(3.251) | |||

(3.252) | |||

(3.253) |

Then

(3.254) |

3.1.21 Statistics Calculations

where recall that

The quantities monitored are:

global rms | (3.258) | |

global rms | (3.259) | |

global rms | (3.260) | |

global average mass times | (3.261) | |

global average mass of moisture | (3.262) |