4. Dynamics

4.1. Finite Volume Dynamical Core

4.1.1. Overview

This document describes the Finite-Volume (FV) dynamical core that was initially developed and used at the NASA Data Assimilation Office (DAO) for data assimilation, numerical weather predictions, and climate simulations. The finite-volume discretization is local and entirely in physical space. The horizontal discretization is based on a conservative “flux-form semi-Lagrangian” scheme described by [LR96] (hereafter LR96) and [LR97] (hereafter LR97). The vertical discretization can be best described as Lagrangian with a conservative re-mapping, which essentially makes it quasi-Lagrangian. The quasi-Lagrangian aspect of the vertical coordinate is transparent to model users or physical parameterization developers, and it functions exactly like the (a hybrid coordinate) used by other dynamical cores within CAM.

In the current implementation for use in CAM, the FV dynamics and physics are “time split” in the sense that all prognostic variables are updated sequentially by the “dynamics” and then the “physics”. The time integration within the FV dynamics is fully explicit, with sub-cycling within the 2D Lagrangian dynamics to stabilize the fastest wave (see section [FVvdisc]). The transport for tracers, however, can take a much larger time step (e.g., 30 minutes as for the physics).

4.1.2. The governing equations for the hydrostatic atmosphere

For reference purposes, we present the continuous differential equations for the hydrostatic 3D atmospheric flow on the sphere for a general vertical coordinate (e.g., [Kas74]). Using standard notations, the hydrostatic balance equation is given as follows:

(1)

where is the density of the air, p the pressure, and g the gravitational constant. Introducing the “pseudo-density” (i.e., the vertical pressure gradient in the general coordinate), from the hydrostatic balance equation the pseudo-density and the true density are related as follows:

(2)

where is the geopotential. Note that reduces to the “true density” if , and the “surface pressure” if (). The conservation of total air mass using as the prognostic variable can be written as

(3)

where . Similarly, the mass conservation law for tracer species (or water vapor) can be written as

(4)

where q is the mass mixing ratio (or specific humidity) of the tracers (or water vapor).

Choosing the (virtual) potential temperature as the thermodynamic variable, the first law of thermodynamics is written as

(5)

Letting denote the (longitude, latitude) coordinate, the momentum equations can be written in the “vector-invariant form” as follows:

(6)

(7)

where A is the radius of the earth, is the coefficient for the optional divergence damping, D is the horizontal divergence

and , the vertical component of the absolute vorticity, is defined as follows:

where is the angular velocity of the earth. Note that the last term in (6) and (7) vanishes if the vertical coordinate is a conservative quantity (e.g., entropy under adiabatic conditions [HA90] or an imaginary conservative tracer), and the 3D divergence operator becomes 2D along constant surfaces. The discretization of the 2D horizontal transport process is described in section [FVhdisc]. The complete dynamical system using the Lagrangian control-volume vertical discretization is described in section [FVvdisc] and section [sec:damp] describes the explicit diffusion operators available in CAM5. A mass, momentum, and total energy conservative mapping algorithm is described in section [FVmap] and in section [sec:geo] an alternative geopotential conserving vertical remapping method is described. Sections [FVqconserve] and [sec:neg] are on the adjusctment of pressure to include the change in mass of water vapor and on the negative tracer fixer in CAM, respectively. Last the global energy fixer is described (section [sec:Global-Energy-Fixer]).

4.1.3. Horizontal discretization of the transport process on the sphere

Since the vertical transport term would vanish after the introduction of the vertical Lagrangian control-volume discretization (see section [FVvdisc]), we shall present here only the 2D (horizontal) forms of the FFSL transport algorithm for the transport of density (3) and mixing ratio-like quantities (4) on the sphere. The governing equation for the pseudo-density (3) becomes

(8)

The finite-volume (integral) representation of the continuous field is defined as follows:

Given the exact 2D wind field the 2D integral representation of the conservation law for can be obtained by integrating (8) in time and in space

The above 2D transport equation is still exact for the finite-volume under consideration. To carry out the contour integral, certain approximations must be made. LR96 essentially decomposed the flux integral using two orthogonal 1D flux-form transport operators. Introducing the following difference operator

and assuming is the time-averaged (from time to time ) on the C-grid (e.g., Fig. 1 in LR96), the 1-D finite-volume flux-form transport operator F in the -direction is

(9)

where , the time-accumulated (from t to t+t) mass flux across the cell wall, is defined as follows,

(10)

and

(11)

can be interpreted as a time mean (from time to time :math:` t+Delta t`) pseudo-density value of all material that passed through the cell edge from the upwind direction.

Note that the above time integration is to be carried out along the backward-in-time trajectory of the cell edge position from (the arrival point; (e.g., point B in Fig. 3 of LR96) back to time (the departure point; e.g., point B’ in Fig. 3 of LR96). The very essence of the 1D finite-volume algorithm is to construct, based on the given initial cell-mean values of :math:` widetilde{pi }`, an approximated subgrid distribution of the true field, to enable an analytic integration of (11). Assuming there is no error in obtaining the time-mean wind , the only error produced by the 1D transport scheme would be solely due to the approximation to the continuous distribution of within the subgrid under consideration (this is not the case in 2D; Lauritzen, Ullrich, and Nair (2010)). From this perspective, it can be said that the 1D finite-volume transport algorithm combines the time-space discretization in the approximation of the time-mean cell-edge values . The physically correct way of approximating the integral (11) must be “upwind”, in the sense that it is integrated along the backward trajectory of the cell edges. For example, a center difference approximation to (11) would be physically incorrect, and consequently numerically unstable unless artificial numerical diffusion is added.

Central to the accuracy and computational efficiency of the finite-volume algorithms is the degrees of freedom that describe the subgrid distribution. The first order upwind scheme, for example, has zero degrees of freedom within the volume as it is assumed that the subgrid distribution is piecewise constant having the same value as the given volume-mean. The second order finite-volume scheme (e.g., [LCSW94]) assumes a piece-wise linear subgrid distribution, which allows one degree of freedom for the specification of the “slope” of the linear distribution to improve the accuracy of integrating (11). The Piecewise Parabolic Method (PPM, [CW84]) has two degrees of freedom in the construction of the second order polynomial within the volume, and as a result, the accuracy is significantly enhanced. The PPM appears to strike a good balance between computational efficiency and accuracy. Therefore, the PPM is the basic 1D scheme we chose (see, e.g., [Mac98]). Note that the subgrid PPM distributions are compact, and do not extend beyond the volume under consideration. The accuracy is therefore significantly better than the order of the chosen polynomials implies. While the PPM scheme possesses all the desirable attributes (mass conserving, monotonicity preserving, and high-order accuracy) in 1D, it is important that a solution be found to avoid the directional splitting in the multi-dimensional problem of modeling the dynamics and transport processes of the Earth’s atmosphere.

The first step for reducing the splitting error is to apply the two orthogonal 1D flux-form operators in a directionally symmetric way. After symmetry is achieved, the “inner operators” are then replaced with corresponding advective-form operators (in CAM5 the “inner operators” are based on constant cell-average values and not the PPM). A stability analysis of the consequences of using different inner and outer operators in the LR96 scheme is given in Lauritzen (2007). A consistent advective-form operator in the direction can be derived from its flux-form counterpart ( as follows:

(12)

(13)

where is a dimensionless number indicating the degree of the flow deformation in the -direction. The above derivation of is slightly different from LR96’s approach, which adopted the traditional 1D advective-form semi-Lagrangian scheme. The advantage of using (12) is that computation of winds at cell centers (Eq. 2.25 in LR96) are avoided.

Analogously, the 1D flux-form transport operator G in the latitudinal (:math:`theta`) direction is derived as follows:

(14)

and likewise the advective-form operator,

(15)

where

(16)

To complete the construction of the 2D algorithm on the sphere, we introduce the following short hand notations:

(17)

(18)

The 2D transport algorithm (cf, Eq. 2.24 in LR96) can then be written as

(19)

Using explicitly the mass fluxes , (19) is rewritten as

(20)

where , the mass flux in the meridional direction, is defined in a similar fashion as (10). The ability of the LR96 scheme to approximate the exact geometry of the fluxes for deformational flows is discussed in Machenhauer, Kaas, and Lauritzen (2009) and Lauritzen, Ullrich, and Nair (2010).

It can be verified that in the special case of constant density flow ( the above equation degenerates to the finite-difference representation of the incompressibility condition of the “time mean” wind field , i.e.,

(21)

The fulfillment of the above incompressibility condition for constant density flows is crucial to the accuracy of the 2D flux-form formulation. For transport of volume mean mixing ratio-like quantities the mass fluxes as defined previously should be used as follows

(22)

Note that the above form of the tracer transport equation consistently degenerates to (19) if (i.e., the tracer density equals to the background air density), which is another important condition for a flux-form transport algorithm to be able to avoid generation of noise (e.g., creation of artificial gradients) and to maintain mass conservation.

4.1.4. A vertically Lagrangian and horizontally Eulerian control-volume discretization of the hydrodynamics

The very idea of using Lagrangian vertical coordinate for formulating governing equations for the atmosphere is not entirely new. [Sta45]) is likely the first to have formulated, in the continuous differential form, the governing equations using a Lagrangian coordinate. Starr did not make use of the discrete Lagrangian control-volume concept for discretization nor did he present a solution to the problem of computing the pressure gradient forces. In the finite-volume discretization to be described here, the Lagrangian surfaces are treated as the bounding material surfaces of the Lagrangian control-volumes within which the finite-volume algorithms developed in LR96, LR97, and L97 will be directly applied.

To use a vertical Lagrangian coordinate system to reduce the 3D governing equations to the 2D forms, one must first address the issue of whether it is an inertial coordinate or not. For hydrostatic flows, it is. This is because both the right-hand-side and the left-hand-side of the vertical momentum equation vanish for purely hydrostatic flows.

Realizing that the earth’s surface, for all practical modeling purposes, can be regarded as a non-penetrable material surface, it becomes straightforward to construct a terrain-following Lagrangian control-volume coordinate system. In fact, any commonly used terrain-following coordinates can be used as the starting reference (i.e., fixed, Eulerian coordinate) of the floating Lagrangian coordinate system. To close the coordinate system, the model top (at a prescribed constant pressure) is also assumed to be a Lagrangian surface, which is the same assumption being used by practically all global hydrostatic models.

The basic idea is to start the time marching from the chosen terrain-following Eulerian coordinate (e.g., pure or hybrid -p), treating the initial coordinate surfaces as material surfaces, the finite-volumes bounded by two coordinate surfaces, i.e., the Lagrangian control-volumes, are free vertically, to float, compress, or expand with the flow as dictated by the hydrostatic dynamics.

By choosing an imaginary conservative tracer that is a monotonic function of height and constant on the initial reference coordinate surfaces (e.g., the value of “” in the hybrid coordinate used in CAM), the 3D governing equations written for the general vertical coordinate in section 1.2 can be reduced to 2D forms. After factoring out the constant , (3), the conservation law for the pseudo-density (), becomes

(23)

where the symbol represents the vertical difference between the two neighboring Lagrangian surfaces that bound the finite control-volume. From (1), the pressure thickness of that control-volume is proportional to the total mass, i.e., . Therefore, it can be said that the Lagrangian control-volume vertical discretization has the hydrostatic balance built-in, and can be regarded as the “pseudo-density” for the discretized Lagrangian vertical coordinate system.

Similarly, (4), the mass conservation law for all tracer species, is

(24)

the thermodynamic equation, (5), becomes

(25)

and (6) and (7), the momentum equations, are reduced to

(26)

(27)

Given the prescribed pressure at the model top , the position of each Lagrangian surface (horizontal subscripts omitted) is determined in terms of the hydrostatic pressure as follows:

(28)

where the subscript is the vertical index ranging from 1 at the lower bounding Lagrangian surface of the first (the highest) layer to at the Earth’s surface. There are +1 Lagrangian surfaces to define a total number of Lagrangian layers. The surface pressure, which is the pressure at the lowest Lagrangian surface, is easily computed as using (28). The surface pressure is needed for the physical parameterizations and to define the reference Eulerian coordinate for the mapping procedure (to be described in section A mass, momentum, and total energy conserving mapping algorithm).

With the exception of the pressure-gradient terms and the addition of a thermodynamic equation, the above 2D Lagrangian dynamical system is the same as the shallow water system described in LR97. The conservation law for the depth of fluid in the shallow water system of LR97 is replaced by (23) for the pressure thickness . The ideal gas law, the mass conservation law for air mass, the conservation law for the potential temperature (25), together with the modified momentum equations (26) and (27) close the 2D Lagrangian dynamical system, which are vertically coupled only by the hydrostatic relation (see (51)), section [FVmap].

The time marching procedure for the 2D Lagrangian dynamics follows closely that of the shallow water dynamics fully described in LR97. For computational efficiency, we shall take advantage of the stability of the FFSL transport algorithm by using a much larger time step ( for the transport of all tracer species (including water vapor). As in the shallow water system, the Lagrangian dynamics uses a relatively small time step, , where is the number of the sub-cycling needed to stabilize the fastest wave in the system. We shall describe here this time-split procedure for the prognostic variables on the D-grid. Discretization on the C-grid for obtaining the diagnostic variables, the time-averaged winds , is analogous to that of the D-grid (see also LR97).

Introducing the following short hand notations (cf, (17) and (18)):

and applying directly (20), the update of “pressure thickness” , using the fractional time step , can be written as

(29)

where are the background air mass fluxes, which are then used as input to Eq. 24 for transport of the potential temperature :

(30)

The discretized momentum equations for the shallow water system (cf, Eq. 16 and Eq. 17 in LR97) are modified for the pressure gradient terms as follows:

(31)

(32)

where is the upwind-biased “kinetic energy” (as defined by Eq. 18 in LR97), and , the horizontal divergence on the D-grid, is discretized as follows:

The finite-volume mean pressure-gradient terms in (31) and (32) are computed as follows:

(33)

(34)

where , and the symbols “” and “” indicate that the contour integrations are to be carried out, using the finite-volume algorithm described in L97, in the and space, respectively.

To complete one time step, equations (29)-[v], together with their counterparts on the C-grid are cycled times using the fractional time step , which are followed by the tracer transport using (24) with the large-time-step .

Mass fluxes and the winds on the C-grid are accumulated for the large-time-step transport of tracer species (including water vapor) as

(35)

where the time-accumulated mass fluxes are computed as

(36)

(37)

The time-averaged winds , defined as follows, are to be used as input for the computations of and :math:` q^{theta }:`

(38)

(39)

The use of the time accumulated mass fluxes and the time-averaged winds for the large-time-step tracer transport in the manner described above ensures the conservation of the tracer mass and maintains the highest degree of consistency possible given the time split integration procedure. A graphical illustration of the different levels of sub-cycling in CAM5 is given on Figure [fig:subc].

Figure 3.1: A graphical illustration of the different levels of sub-cycling in CAM5.

The algorithm described here can be readily applied to a regional model if appropriate boundary conditions are supplied. There is formally no Courant number related time step restriction associated with the transport processes. There is, however, a stability condition imposed by the gravity-wave processes. For application on the whole sphere, it is computationally advantageous to apply a polar filter to allow a dramatic increase of the size of the small time step . The effect of the polar filter is to stabilize the short-in-wavelength (and high-in-frequency) gravity waves that are being unnecessarily and unidirectionally resolved at very high latitudes in the zonal direction. To minimize the impact to meteorologically significant larger scale waves, the polar filter is highly scale selective and is applied only to the diagnostic variables on the auxiliary C-grid and the tendency terms in the D-grid momentum equations. No polar filter is applied directly to any of the prognostic variables.

The design of the polar filter follows closely that of [ST95] for the C-grid Arakawa type dynamical core (e.g., [AL81]). For the CAM5.0 the fast-fourier transform component of the polar filtering has replaced the algebraic form at all filtering latitudes. Because our prognostic variables are computed on the D-grid and the fact that the FFSL transport scheme is stable for Courant number greater than one, in realistic test cases the maximum size of the time step is about two to three times larger than a model based on Arakawa and Lamb’s C-grid differencing scheme. It is possible to avoid the use of the polar filter if, for example, the “Cubed grid” is chosen, instead of the current latitude-longitude grid. rewrite of the rest of the model codes including physics parameterizations, the land model, and most of the post processing packages.

The size of the small time step for the Lagrangian dynamics is only a function of the horizontal resolution. Applying the polar filter, for the 2-degree horizontal resolution, a small-time-step size of 450 seconds can be used for the Lagrangian dynamics. From the large-time-step transport perspective, the small-time-step integration of the 2D Lagrangian dynamics can be regarded as a very accurate iterative solver, with m iterations, for computing the time mean winds and the mass fluxes, analogous in functionality to a semi-implicit algorithm’s elliptic solver (e.g., [RHaDAR00]). Besides accuracy, the merit of an “explicit” versus “semi-implicit” algorithm ultimately depends on the computational efficiency of each approach. In light of the advantage of the explicit algorithm in parallelization, we do not regard the explicit algorithm for the Lagrangian dynamics as an impedance to computational efficiency, particularly on modern parallel computing platforms.

4.1.5. Optional diffusion operators in CAM5

The ‘CD’-grid discretization method used in the CAM finite-volume dynamical core provides explicit control over the rotational modes at the grid scale, due to monotonicity constraint in the PPM-based advection, but there is no explicit control over the divergent modes at the grid scale Skamarock (see, e.g., 2010). Therefore divergence damping terms appear on the right-hand side of the momentum equations (26) and (27):

(40)

and

(41)

respectively, where the strength of the divergence damping is controlled by the coefficient given by

(42)

where throughout the atmosphere except in the top model levels where it monotonically increases to approximately at the top of the atmosphere. The divergence damping described above is referred to as ‘second-order’ divergence damping as it effectively damps divergence with a operator.

In CAM5 optional ‘fourth-order’ divergence damping has been implemented where the divergence is effectively damped with a -operator which is usually more scale selective than ‘second-order’ damping operators. For ‘fourth-order’ divergence damping the terms

(43)

and

(44)

are added to the right-hand side of ([u-lcv]) and ([v-lcv]), respectively. The horizontal Laplacian -operator in spherical coordinates for a scalar is given by

The fourth-order divergence damping coefficient is given by

(45)

Since divergence damping is added explicitly to the equations of motion it is unstable if the time-step is too large or the damping coefficients ( or ) are too large. To stabilize the fourth-order divergence damping the winds used to compute the divergence are filtered using the same FFT filtering which is applied to stabilize the gravity waves.

To control potentially excessive polar night jets in high-resolution configurations of CAM, Laplacian damping of the wind components has been added as an option in CAM5. That is, the terms

(46)

and

(47)

are added to the right-hand side of the momentum equations ([u-lcv]) and ([v-lcv]), respectively. The damping coefficient is zero throughout the atmosphere except in the top layers where it increases monotonically and smoothly from zero to approximately four times a user-specified damping coefficient at the top of the atmosphere (the user-specified damping coefficient is typically on the order of msec).

4.1.6. A mass, momentum, and total energy conserving mapping algorithm

The Lagrangian surfaces that bound the finite-volume will eventually deform, particularly in the presence of persistent diabatic heating/cooling, in a time scale of a few hours to a day depending on the strength of the heating and cooling, to a degree that it will negatively impact the accuracy of the horizontal-to-Lagrangian-coordinate transport and the computation of the pressure gradient forces. Therefore, a key to the success of the Lagrangian control-volume discretization is an accurate and conservative algorithm for mapping the deformed Lagrangian coordinate back to a fixed reference Eulerian coordinate.

There are some degrees of freedom in the design of the vertical mapping algorithm. To ensure conservation, our current (and recommended) mapping algorithm is based on the reconstruction of the “mass” (pressure thickness ), zonal and meridional “winds”, “tracer mixing ratios”, and “total energy” (volume integrated sum of the internal, potential, and kinetic energy), using the monotonic Piecewise Parabolic sub-grid distributions with the hydrostatic pressure (as defined by ([L-coord])) as the mapping coordinate. We outline the mapping procedure as follows.

Step 1: Define a suitable Eulerian reference coordinate as a target coordinate. The mass in each layer () is then distributed vertically according to the chosen Eulerian coordinate. The surface pressure typically plays an “anchoring” role in defining the terrain following Eulerian vertical coordinate. The hybrid used in the NCAR CCM3 ([KHB+96]) is adopted in the current model setup.

Step 2: Construct the piece-wise continuous vertical subgrid profiles of tracer mixing ratios (), zonal and meridional winds (u and v), and total energy () in the Lagrangian control-volume coordinate, or the source coordinate. The total energy is computed as the sum of the finite-volume integrated geopotential , internal energy , and the kinetic energy () as follows:

(48)

Applying integration by parts and the ideal gas law, the above integral can be rewritten as

(49)

where is the layer mean virtual temperature, is the layer mean kinetic energy, is the pressure at layer edges, and and are the specific heat of the air at constant volume and at constant pressure, respectively. The total energy in each grid cell is calculated as

The method employed to create subgrid profiles is set by the flag . For = 0 (default), the Piece-wise Parabolic Method (PPM, [CW84]) over a pressure coordinate is used and for a cublic spline over a logarithmic pressure coordinate is used.

Step 3: Layer mean values of , (u, v), and in the Eulerian coordinate system are obtained by integrating analytically the sub-grid distributions, in the vertical direction, from model top to the surface, layer by layer. Since the hydrostatic pressure is chosen as the mapping coordinate, tracer mass, momentum, and total energy are locally and globally conserved. In mapping a variable from the source coordinate to the target coordinate, different limiter constraints may be used and they are controlled by two flags, and . For winds on D-grid, should be set to -1. For tracers, should be set to 0. For all others, . directly controls which limiter constraint is used. For , Huynh’s 2nd constraint is used. If , the original quasi-monotonic constraint is used. If , a full monotonic constraint is used. If is less than 7, the variable, , is determined by the following:

If , a standard PPM constraint is used. If , an improved full monotonicity constraint is used. If , a positive definite constraint is used. If , the algorithm will do nothing.

Step 4: Retrieve virtual temperature in the Eulerian (target) coordinate. Start by computing kinetic energy in the Eulerian coordinate system for each layer. Then substitute kinetic energy and the hydrostatic relationship into ([TE:sub:fv]). The layer mean temperature for layer in the Eulerian coordinate is then retrieved from the reconstructed total energy (done in Step 3) by a fully explicit integration procedure starting from the surface up to the model top as follows:

(50)

where and is the gas constant for dry air.

To convert the potential virtual temperature to the layer mean temperature the conversion factor is obtained by equating the following two equivalent forms of the hydrostatic relation for and

(51)

(52)

where . The conversion formula between layer mean temperature and layer mean potential temperature is obtained as follows:

(53)

The physical implication of retrieving the layer mean temperature from the total energy as described in Step 3 is that the dissipated kinetic energy, if any, is locally converted into internal energy via the vertically sub-grid mixing (dissipation) processes. Due to the monotonicity preserving nature of the sub-grid reconstruction the column-integrated kinetic energy inevitably decreases (dissipates), which leads to local frictional heating. The frictional heating is a physical process that maintains the conservation of the total energy in a closed system.

As viewed by an observer riding on the Lagrangian surfaces, the mapping procedure essentially performs the physical function of the relative-to-the-Eulerian-coordinate vertical transport, by vertically redistributing (air and tracer) mass, momentum, and total energy from the Lagrangian control-volume back to the Eulerian framework.

As described in section [FVvdisc], the model time integration cycle consists of small time steps for the 2D Lagrangian dynamics and one large time step for tracer transport. The mapping time step can be much larger than that used for the large-time-step tracer transport. In tests using the Held-Suarez forcing ([HS94]), a three-hour mapping time interval is found to be adequate. In the full model integration, one may choose the same time step used for the physical parameterizations so as to ensure the input state variables to physical parameterizations are in the usual “Eulerian” vertical coordinate. In CAM5, vertical remapping takes place at each physics time step.

4.1.7. A geopotential conserving mapping algorithm

An alternative vertical mapping approach is available in CAM5. Instead of retrieving temperature by remapped total energy in the Eulerian coordinate, the alternative approach maps temperature directly from the Lagrangian coordinate to the Eulerian coordinate. Since geopotential is defined as

mapping over or over preserves the geopotential at the model lid. This approach prevents the mapping procedure from generating spurious pressure gradient forces at the model lid. Unlike the energy-conserving algorithm which could produce substantial temperature fluctuations at the model lid, the geopotential conserving approach guarantees a smooth (potential) temperature profile. However, the geopotential conserving does not conserve total energy in the remapping procedure. This may be resolved by a global energy fixer already implemented in the model (see section [sec:Global-Energy-Fixer]).

4.1.8. Adjustment of pressure to include change in mass of water vapor

The physics parameterizations operate on a model state provided by the dynamics, and are allowed to update specific humidity. However, the surface pressure remains fixed throughout the physics updates, and since there is an explicit relationship between the surface pressure and the air mass within each layer, the total air mass must remain fixed as well throughout the physics updates. If no further correction were made, this would imply that the dry air mass changed if the water vapor mass changed in the physics updates. Therefore the pressure field is changed to include the change in water vapor mass due to the physics updates. We impose the restrictions that dry air mass and water mass are conserved as follows:

The total pressure is

with dry pressure , water vapor pressure . The specific humidity is

We define a layer thickness as , so

We are concerned about 3 time levels: is input to physics, is output from physics, is the adjusted value for dynamics.

Dry mass is the same at and but not at . To conserve dry mass, we require that

or

(54)

Water mass is the same at and , but not at . To conserve water mass, we require that

(55)

Substituting ([eq:wetdif]) into ([eq:drydif]),

which yields a modified specific humidity for the dynamics:

We note that this correction as implemented makes a small change to the water vapor as well. The pressure correction could be formulated to leave the water vapor unchanged.

4.1.9. Negative Tracer Fixer

In the Finite Volume dynamical core, neither the monotonic transport nor the conservative vertical remapping guarantee that tracers will remain positive definite. Thus the Finite Volume dynamical core includes a negative tracer fixer applied before the parameterizations are calculated. For negative mixing ratios produced by horizontal transport, the model will attempt to borrow mass from the east and west neighboring cells. In practice, most negative values are introduced by the vertical remapping which does not guarantee positive definiteness in the first and last layer of the vertical column.

A minimum value is defined for each tracer. If the tracer falls below that minimum value, it is set to that minimum value. If there is enough mass of the tracer in the layer immediately above, tracer mass is removed from that layer to conserve the total mass in the column. If there is not enough mass in the layer immediately above, no compensation is applied, violating conservation. Usually such computational sources are very small.

The amount of tracer needed from the layer above to bring up to is

where is the vertical index, increasing downward. After the filling

Currently for water vapor, for CLDLIQ, CLDICE, NUMLIQ and NUMICE, and for the remaining constituents.

4.1.10. Global Energy Fixer

The finite-volume dynamical core as implemented in CAM and described here conserves the dry air and all other tracer mass exactly without a “mass fixer”. The vertical Lagrangian discretization and the associated remapping conserves the total energy exactly. The only remaining issue regarding conservation of the total energy is the horizontal discretization and the use of the “diffusive” transport scheme with monotonicity constraint. To compensate for the loss of total energy due to horizontal discretization, we apply a global fixer to add the loss in kinetic energy due to “diffusion” back to the thermodynamic equation so that the total energy is conserved. The loss in total energy (in flux unit) is found to be around 2 ) with the 2 degrees resolution.

The energy fixer is applied following the negative tracer fixer. The fixer is applied on the unstaggered physics grid rather than on the staggered dynamics grid. The energies on these two grids are difficult to relate because of the nonlinear terms in the energy definition and the interpolation of the state variables between the grids. The energy is calculated in the parameterization suite before the state is passed to the finite volume core as described in the beginning of Chapter [chap:model:sub:physics]. The fixer is applied just before the parameterizations are calculated. The fixer is a simplification of the fixer in the Eulerian dynamical core described in section [energyfixer].

Let minus sign superscript denote the values at the beginning of the dynamics time step, i.e. after the parameterizations are applied, let a plus sign superscript denote the values after fixer is applied, and let a hat denote the provisional value before adjustment. The total energy over the entire computational domain after the fixer is

where is the latent heat of vaporation, is the latent heat of fusion, is water vapor mixing ratio, and is cloud water mixing ratio. should equal the energy at the beginning of the dynamics time step

Let denote the energy of the provisional state provided by the dynamical core before the adjustment.

Thus, the total energy added into the system by the dynamical core is . The energy fixer then changes dry static energy () by a constant amount over each grid cell to conserve total energy in the entire computational domain. The dry static energy added to each grid cell may be expressed as

Therefore,

or

(56)

This will ensure .

By hydrostatic approximation, the geopotential equation is

and for any arbitrary point between and the geopotential may be written as

The geopotential at the mid point of a model layer between and , or the layer mean, is

For layer , the energy fixer will solve the following equation based on ([dry:sub:static_eqn]),

Since the energy fixer will not alter the water vapor mixing ratio and the pressure field,

Therefore,

The energy fixer starts from the Earth’s surface and works its way up to the model top in adjusting the temperature field. At the surface layer, . After the temperature is adjusted in a grid cell, the geopotential at the upper interface of the cell is updated which is needed for the temperature adjustment in the grid cell above.

4.1.11. Further discussion

There are still aspects of the numerical formulation in the finite volume dynamical core that can be further improved. For example, the choice of the horizontal grid, the computational efficiency of the split-explicit time marching scheme, the choice of the various monotonicity constraints, and how the conservation of total energy is achieved.

The impact of the non-linear diffusion associated with the monotonicity constraint is difficult to assess. All discrete schemes must address the problem of subgrid-scale mixing. The finite-volume algorithm contains a non-linear diffusion that mixes strongly when monotonicity principles are locally violated. However, the effect of nonlinear diffusion due to the imposed monotonicity constraint diminishes quickly as the resolution matches better to the spatial structure of the flow. In other numerical schemes, however, an explicit (and tunable) linear diffusion is often added to the equations to provide the subgrid-scale mixing as well as to smooth and/or stabilize the time marching.

4.1.12. Specified Dynamics Option

In CAM4 the capability included to perform simulations using specified dynamics, where offline meteorological fields are nudged to the online calculated meteorology. This procedure was originally used in the Model of Atmospheric Transport and Chemistry (MATCH) (Rasch et al., 1997). In this procedure the horizontal wind components, air temperature, surface temperature, surface pressure, sensible and latent heat flux, and wind stress are read into the model simulation from the input meteorological dataset. The nudging coefficient can be chosen to be 1 (for 100% nudging) or smaller. The desired percentage of the offline meteorology and the remaining percent from the internally calcuated meteorology is used every timestep to prescribe the meteorological parameters. In addition, the model solves the model internal advection equations for the mass flux every sub-step. In this way, some inconsistencies between the inserted and model-computed velocity and mass fields subsequently used for tracer transport are dampended. The mass flux at each sub-step is accumulated to produce the net mass flux over the entire time step. A graphical explanation of the sub-cycling is given in Lauritzen et al. (2011).

A nudging coefficent of 100 can be used to allow for more precise comparisons between measurements of atmospheric composition and model output for example using CAM-Chem (Lamarque et al., 2012). A reduced nudging coefficent is used for instant for WACCM simulations, if more of the internal transport parameters needs to be contained, while the meteorology is still close to the analysied fields (e.g., Brakebusch et al., 2012).

Currently, we recommend for input offline meteorology interpolated from 0.5x0.6 degree fields of the NASA Goddard Global Modeling and Assimilation Office (GMAO) GEOS-5 and Modern Era Retrospective-Analysis For Research And Applications (MERRA) generated meteorology. These fields are available on the Earth System Grid (http://www.earthsystemgrid.org/home.htm) for the CAM resolution of 1.9x2.5. These files were generated from the original resolution by using a conservative regridding procedure based on the same 1-D operators as used in the transport scheme of the finite-volume dynamical core used in GEOS-5 and CAM (S.-J. Lin, personal communication, 2009). Note that because of a difference in the sign convention of the surface wind stress (TAUX and TAUY) between CESM and GEOS5/MERRA, these fields in the interpolated datasets have been reversed from the original files supplied by GMAO. In addition, it is important for users to recognize the importance of specifying the correct surface geopotential height (PHIS) to ensure consistency with the input dynamical fields, which is important to prevent unrealistic vertical mixing.

4.1.13. Further discussion

4.2. Spectral Element Dynamical Core

The CAM includes an optional dynamical core from HOMME, NCAR’s High-Order Method Modeling Environment [DFS+05]. The stand-alone HOMME is used for research in several different types of dynamical cores. The dynamical core incorporated into CAM4 uses HOMME’s continuous Galerkin spectral finite element method [TTI97],:cite:fournier04,:cite:thomas05,:cite:wang07,:cite:taylor10b, here abbreviated to the spectral element method (SEM). This method is designed for fully unstructured quadrilateral meshes. The current configurations in the CAM are based on the cubed-sphere grid. The main motivation for the inclusion of HOMME is to improve the scalability of the CAM by introducing quasi-uniform grids which require no polar filters [TEStCyr08]. HOMME is also the first dynamical core in the CAM which locally conserves energy in addition to mass and two-dimensional potential vorticity [Tay10].

HOMME represents a large change in the horizontal grid as compared to the other dynamical cores in CAM. Almost all other aspects of HOMME are based on a combination of well-tested approaches from the Eulerian and FV dynamical cores. For tracer advection, HOMME is modeled as closely as possible on the FV core. It uses the same conservation form of the transport equation and the same vertically Lagrangian discretization [Lin04]. The HOMME dynamics are modeled as closely as possible on Eulerian core. They share the same vertical coordinate, vertical discretization, hyper-viscosity based horizontal diffusion, top-of-model dissipation, and solve the same moist hydrostatic equations. The main differences are that HOMME advects the surface pressure instead of its logarithm (in order to conserve mass and energy), and HOMME uses the vector-invariant form of the momentum equation instead of the vorticity-divergence formulation. Several dry dynamical cores including HOMME are evaluated in [LJTN09] using a grid-rotated version of the baroclinic instability test case [JW06].

The timestepping in HOMME is a form of dynamics/tracer/physics subcycling, achieved through the use of multi-stage 2nd order accurate Runge-Kutta methods. The tracers and dynamics use the same timestep which is controlled by the maximum anticipated wind speed, but the dynamics uses more stages than the tracers in order to maintain stability in the presence of gravity waves. The forcing is applied using a time-split approach. The optimal forcing strategy in HOMME has not yet been determined, so HOMME supports several options. The first option is modeled after the FV dynamical core and the forcing is applied as an adjustment at each physics timestep. The second option is to convert all forcings into tendencies which are applied at the end of each dynamics/tracer timestep. If the physics timestep is larger than the tracer timestep, then the tendencies are held fixed and only updated at each physics timestep. Finally, a hybrid approach can be used where the tracer tendencies are applied as in the first option and the dynamics tendencies are applied as in the second option.

4.2.1. Continuum Formulation of the Equations

HOMME uses a conventional vector-invariant form of the moist primitive equations. For the vertical discretization it uses the hybrid pressure vertical coordinate system modeled after [eul:terrain] The formulation here differs only in that surface pressure is used as a prognostic variable as opposed to its logarithm.

In the -coordinate system, the pressure is given by

The hydrostatic approximation is used to replace the mass density by an -coordinate pseudo-density . The material derivative in -coordinates can be written (e.g. (e.g. [Sat04], Sec.3.3),

where the operator (as well as and below) is the two-dimensional gradient on constant -surfaces, is the vertical derivative, is a vertical flow velocity and is the horizontal velocity component (tangent to constant -surfaces, not -surfaces).

The -coordinate atmospheric primitive equations, neglecting dissipation and forcing terms can then be written as

(57)

These are prognostic equations for , the temperature , density , and where is the specific humidity. The prognostic variables are functions of time , vertical coordinate and two coordinates describing the surface of the sphere. The unit vector normal to the surface of the sphere is denoted by . This formulation has already incorporated the hydrostatic equation and the ideal gas law, . There is a no-flux () boundary condition at and . The vorticity is denoted by , is a Coriolis term and is the pressure vertical velocity. The virtual temperature and variable-of-convenience are defined as in [eul:terrain].

The diagnostic equations for the geopotential height field is

(58)

where is the prescribed surface geopotential height (given at ). To complete the system, we need diagnostic equations for and , which come from integrating with respect to . In fact, can be replaced by a diagnostic equation for and a prognostic equation for surface pressure

(59)

where is evaluated at the model bottom () after using that and . Using Eq [E:PEcont2c], we can derive a diagnostic equation for the pressure vertical velocity ,

Finally, we rewrite as

(60)

4.2.2. Conserved Quantities

The equations have infinitely many conserved quantities, including mass, tracer mass, potential temperature defined by

with ( or ) and the total moist energy defined by

(61)

where is the spherical area measure. To compute these quantities in their traditional units they should be divided by the constant of gravity . We have omitted this scaling since has also been scaled out from –. We note that in this formulation of the primitive equations, the pressure is a moist pressure, representing the effects of both dry air and water vapor. The unforced equations conserve both the moist air mass ( above) and the dry air mass ( ). However, in the presence of a forcing term in (representing sources and sinks of water vapor as would be present in a full model) a corresponding forcing term must be added to to ensure that dry air mass is conserved.

The energy is specific to the hydrostatic equations. We have omitted terms from the physical total energy which are constant under the evolution of the unforced hydrostatic equations ([SWG03]). It can be converted into a more universal form involving , with defined similarly to , so that where and are the specific heats of dry air and water vapor defined at constant volume. We note that and so that . Expanding with this expression, integrating by parts with respect to and making use of the fact that the model top is at a constant pressure

and thus

(62)

The model top boundary term in vanishes if . Otherwise it must be included to be consistent with the hydrostatic equations. It is present due to the fact that the hydrostatic momentum equation neglects the vertical pressure gradient.

4.2.3. Horizontal Discretization: Functional Spaces

In the finite element method, instead of constructing discrete approximations to derivative operators, one constructs a discrete functional space, and then finds the function in this space which solves the equations of interest in a minimum residual sense. As compared to finite volume methods, there is less choice in how one constructs the discrete derivative operators in this setting, since functions in the discrete space are represented in terms of known basis functions whose derivatives are known, often analytically.

Let and be the Cartesian coordinates and position vector of a point in the reference square and let and be the coordinates and position vector of a point on the surface of the sphere, denoted by . We mesh using the cubed-sphere grid (Fig. Figure 3.2: Tiling the surface of the sphere with quadrilaterals. An inscribed cube is projected to the surface of the sphere. The faces of the cubed sphere are further subdivided to form a quadrilateral grid of the desired resolution. Coordinate lines from the gnomonic equal-angle projection are shown.) first used in [Sad72]. Each cube face is mapped to the surface of the sphere with the equal-angle gnomonic projection ([RanvcicPM96]). The map from the reference element to the cube face is a translation and scaling. The composition of these two maps defines a map from the spherical elements to the reference element . We denote this map and its inverse by

(63)

Figure 3.2: Tiling the surface of the sphere with quadrilaterals. An inscribed cube is projected to the surface of the sphere. The faces of the cubed sphere are further subdivided to form a quadrilateral grid of the desired resolution. Coordinate lines from the gnomonic equal-angle projection are shown.

We now define the discrete space used by the SEM. First we denote the space of polynomials up to degree in by

Todo

The above equation is not correct - the following part did not translate correctly and needs to be fixed {mathop{mathrm{span}}}limits_{{smash[t]{vec{imath}}}inmathbb{I}}phi_{{smash[t]{vec{imath}}}}({{smash[t]{vec{x}}}),

where contains all the degrees and , , are the cardinal functions, namely polynomials that interpolate the tensor-product of degree- Gauss-Lobatto-Legendre (GLL) nodes . The GLL nodes used within an element for are shown in Fig. [f:GLLnodes]. The cardinal-function expansion coefficients of a function are its GLL nodal values, so we have

(64)

We can now define the piecewise-polynomial SEM spaces and as

(65)

Functions in are polynomial within each element but may be discontinuous at element boundaries and is the subspace of continuous function in . We take , and . We then construct a set of unique points by

(66)

For every point , there exists at least one element and at least one GLL node . In 2D, if belongs to exactly one it is an element-interior node. If it belongs to exactly two s, it is an element-edge interior node. Otherwise it is a vertex node.

[l]2.0in

2.5in

We also define similar spaces for 2D vectors. We introduce two families of spaces, with a subscript of either con or cov, denoting if the contravariant or covariant components of the vectors are piecewise polynomial, respectively.

where are the contravariant components of defined below. Vectors in are globally continuous and their contravariant components are polynomials in each element. Similarly,

The SEM is a Galerkin method with respect to the subspace and it can be formulated solely in terms of functions in . In CAM-HOMME, the typical configuration is to run with which achieves a 4th order accurate horizontal discretization ([TF10]). All variables in the CAM-HOMME initial condition and history files as well as variables passed to the physics routines are represented by their grid point values at the points . However, for some intermediate quantities and internally in the dynamical core it is useful to consider the larger space, where variables are represented by their grid point values at the mapped GLL nodes. This later representation can also be considered as the cardinal-function expansion of a function local to each element,

(67)

since the expansion coefficients are the function values at the mapped GLL nodes. Functions in can be multiple-valued at GLL nodes that are redundant (i.e., shared by more than one element), while for , the values at any redundant points must all be the same.

4.2.4. Horizontal Discretization: Differential Operators

We use the standard curvilinear coordinate formulas for vector operators following [Hei01]. Given the Jacobian of the the mapping from to , we denote its determinant-magnitude by

(68)

A vector may be written in terms of physical or covariant or contravariant components, or or ,

(69)

that are related by and , where is a contravariant basis vector and is a covariant basis vector.

The dot product and contravariant components of the cross product are [Hei01] Table 1

(70)

where is the Levi-Civita symbol. The divergence, covariant coordinates of the gradient and contravariant coordinates of the curl are [Hei01] [eqs.2.1.1, 2.1.4 and 2.1.6]

(71)

In the SEM, these operators are all computed in terms of the derivatives with respect to in the reference element, computed exactly (to machine precision) by differentiating the local element expansion . For the gradient, the covariant coordinates of are thus computed exactly within each element. Note that , but may not be in even for due to the fact that its components will be multi-valued at element boundaries because computed in adjacent elements will not necessarily agree along their shared boundary. In the case where is constant within each element, the SEM curl of and the divergence of will also be exact, but as with the gradient, multiple-valued at element boundaries.

For non-constant , these operators may not be computed exactly by the SEM due to the Jacobian factors in the operators and the Jacobian factors that appear when converting between covariant and contravariant coordinates. We follow [TL00] and evaluate these operators in the form shown in . The quadratic terms that appear are first projected into via interpolation at the GLL nodes and then this interpolant is differentiated exactly using . For example, to compute the divergence of , we first compute the interpolant of , where the GLL interpolant of a product derives simply from the product of the GLL nodal values of and . This operation is just a reinterpretation of the nodal values and is essentially free in the SEM. The derivatives of this interpolant are then computed exactly from . The sum of partial derivatives are then divided by at the GLL nodal values and thus the SEM divergence operator is given by

(72)

Similarly, the gradient and curl are approximated by

(73)

with and . The SEM is well known for being quite efficient in computing these types of operations. The SEM divergence, gradient and curl can all be evaluated at the GLL nodes within each element in operations per node using the tensor-product property of these points [DFM02],:cite:karniadakis05.

4.2.5. Horizontal Discretization: Discrete Inner-Product

Instead of using exact integration of the basis functions as in a traditional finite-element method, the SEM uses a GLL quadrature approximation for the integral over , that we denote by . We can write this integral as a sum of area-weighted integrals over the set of elements used to decompose the domain,

The integral over a single element is written as an integral over by

(74)

where we approximate the integral over by GLL quadrature,

(75)

The SEM approximation to the global integral is then naturally defined as

(76)

When applied to the product of functions , the quadrature approximation defines a discrete inner-product in the usual manner.

4.2.6. Horizontal Discretization: The Projection Operators

Let be the unique orthogonal (self-adjoint) projection operator from onto w.r.t. the SEM discrete inner product . The operation is essentially the same as the common procedure in the SEM described as assembly [KS05] [p.7] or direct stiffness summation [DFM02] [eq.4.5.8]. Thus the SEM assembly procedure is not an ad-hoc way to remove the redundant degrees of freedom in , but is in fact the natural projection operator . Applying the projection operator in a finite element method requires inverting the finite element mass matrix. A remarkable fact about the SEM is that with the GLL based discrete inner product and the careful choice of global basis functions, the mass matrix is diagonal [MP87]. The resulting projection operator then has a very simple form: at element interior points, it leaves the nodal values unchanged, while at element boundary points shared by multiple elements it is a Jacobian-weighted average over all redundant values ([TF10]).

To apply the projection to vectors , one cannot project the covariant components since the corresponding basis vectors and do not necessarily agree along element faces. Instead we must define the projection as acting on the components using a globally continuous basis such as the latitude-longitude unit vectors and ,

4.2.7. Horizontal Discretization: Galerkin Formulation

The SEM solves a Galerkin formulation of the equations of interest. Given the discrete differential operators described above, the primitive equations can be written as an ODE for a generic prognostic variable and right-hand-side (RHS) terms

The SEM solves this equation in integral form with respect to the SEM inner product. That is, for a , the SEM finds the unique such that

As the prognostic variable is assumed to belong to , the RHS will in general belong to since it contains derivatives of the prognostic variables, resulting in the loss of continuity at the element boundaries. If one picks a suitable basis for , this discrete integral equation results in a system of equations for the expansion coefficients of . The SEM solves these equations exactly, and the solution can be written in terms of the SEM projection operator as

The projection operator commutes with any time-stepping scheme, so the equations can be solved in a two step process, illustrated here for simplicity with the forward Euler method

• Step 1:

  

• Step 2:

  

For compactness of notation, we will denote this two step procedure in what follows by

Note that maps a dimensional space into a dimensional space , so here denotes the left inverse of . This inverse will never be computed, it is only applied as in step 2 above.

This two step Galerkin solution process represents a natural separation between computation and communication for the implementation of the SEM on a parallel computer. The computations in step 1 are all local to the data contained in a single element. Assuming an element-based decomposition so that each processor contains at least one element, no inter-processor communication is required in step 1. All inter-processor communication in HOMME is isolated to the projection operator step, in which element boundary data must be exchanged between adjacent elements.

4.2.8. Vertical Discretization

The vertical coordinate system uses a Lorenz staggering of the variables as shown in Figure 3.1: A graphical illustration of the different levels of sub-cycling in CAM5.. Let be the total number of layers, with variables at layer mid points denoted by . We denote layer interfaces by , so that and . The -integrals will be replaced by sums. We will use to denote the discrete operator. The operator uses centered differences to compute derivatives with respect to at layer mid point from layer interface values, . We will use the over-bar notation for vertical averaging, . We also introduce the symbol to denote the discrete pseudo-density given by

.

We will use to denote the discrete form of the operator. We use the discretization given in [eul:econserve]. This operator acts on quantities defined at layer mid-points and returns a result also at layer mid-points,

(77)

where . We use the over-bar notation since the formula can be seen as a -weighted average of a layer interface centered difference approximation to . This formulation was constructed in [SB81] in order to ensure mass and energy conservation. Here we will use an equivalent expression that can be written in terms of ,

(78)

4.2.9. Discrete formulation: Dynamics

We discretize the equations exactly in the form shown in , , and , obtaining

(79)

We consider a single quantity given at layer interfaces and defined by . The no-flux boundary condition is . In , we used a midpoint quadrature rule to evaluate the indefinite integral from . In practice can be eliminated from the discrete equations by scaling , but here we retain them so as to have a direct correspondence with the continuum form of the equations written in terms of .

Finally we give the approximations for the diagnostic equations. We first integrate to layer interface using the same mid-point rule as used to derive , and then add an additional term representing the integral from to :

(80)

where

and similar for ,

(81)

where

Similar to [eul:econserve], we note that

(82)

which ensures energy conservation ([Tay10]).

4.2.10. Consistency

It is important that the discrete equations be as consistent as possible. In particular, we need a discrete version of , the non-vertically averaged continuity equation. Equation implicitly implies such an equation. To see this, apply to and using that then we can derive, at layer mid-points,

(83)

A second type of consistency that has been identified as important is that , the discrete equation for , be consistent with , the discrete continuity equation ([WO94b]). The two discrete equations should imply a reasonable discretization of . To show this, we take the average of at layers and and combine this with (at layer mid-points ) and assuming that we obtain

which, since is given at layer mid-points and at layer interfaces, is the SEM discretization of .

4.2.11. Time Stepping

Applying the SEM discretization to - results in a system of ODEs. These are solved with an -stage Runge-Kutta method. This method allows for a gravity-wave based CFL number close to , (normalized so that the largest stable timestep of the Robert filtered Leapfrog method has a CFL number of 1.0). The value of is chosen large enough so that the dynamics will be stable at the same timestep used by the tracer advection scheme. To determine , we first note that the tracer advection scheme uses a less efficient (in terms of maximum CFL) strong stability preserving Runge-Kutta method described below. It is stable at an advective CFL number of 1.4. Let be a maximum wind speed and be the maximum gravity wave speed. The gravity wave and advective CFL conditions are

In the case where is chosen as the largest stable timestep for advection, then we require for a stable dynamics timestep. Using a typical values m/s and m/s gives . CAM places additional restrictions on the timestep (such as that the physics timestep must be an integer multiple of ) which also influence the choice of and .

4.2.12. Dissipation

A horizontal hyper-viscosity operator, modeled after [eul:hdiff] is applied to the momentum and temperature equations. It is applied in a time-split manor after each dynamics timestep. The hyper-viscosity step for vectors can be written as

An integral form of this equation suitable for the SEM is obtained using a mixed finite element formulation ([Gir99]) which writes the equation as a system of equations involving only first derivatives. We start by introduced an auxiliary vector and using the identity ,

(84)

Integrating the gradient and curl operators by parts gives

(85)

The SEM Galerkin solution of this integral equation is most naturally written in terms of an inverse mass matrix instead of the projection operator. It can be written in terms of the SEM projection operator by first testing with the product of the element cardinal functions and the contravariant basis vector . With this type of test function, the RHS of can be defined as a weak Laplacian operator . The covariant components of given by are then

Then the SEM solution to and is given by

Because of the SEM tensor product decomposition, the expression for can be evaluated in only operations per grid point, and in CAM-HOMME typically .

Following [eul:hdiff], a correction term is added so the hyper-viscosity does not damp rigid rotation. The hyper-viscosity formulation used for scalars such as is much simpler, since instead of the vector Laplacian identity we use . Otherwise the approach is identical to that used above so we omit the details. The correction for terrain following coordinates given in [eul:hdiff] is not yet implemented in CAM-HOMME.

4.2.13. Discrete formulation: Tracer Advection

All tracers, including specific humidity, are advected with a discretized version of . HOMME uses the vertically Lagrangian approach (see [FVvdisc]) from [Lin04]. At the beginning of each timestep, the tracers are assumed to be given on the -coordinate layer mid points. The tracers are advanced in time on a moving vertical coordinate system defined so that . At the end of the timestep, the tracers are remapped back to the -coordinate layer mid points using the monotone remap algorithm from [ZWS05].

The horizontal advection step consists of using the SEM to solve

(86)

on the surfaces defined by the layer mid points. The quantity is the mean flux computed during the dynamics update. The mean flux used in , combined with a suitable mean vertical flux used in the remap stage allows HOMME to preserve mass/tracer-mass consistency: The tracer advection of with will be identical to the advection of implied from . The mass/tracer-mass consistency capability is not in the version of HOMME included in CAM 4.0, but should be in all later versions.

The equation is discretized in time using the optimal 3 stage strong stability preserving (SSP) second order Runge-Kutta method from [SR02]. The RK-SSP method is chosen because it will preserve the monotonicity properties of the horizontal discretization. RK-SSP methods are convex combinations of forward-Euler timesteps, so each stage of the RK-SSP timestep looks like

(87)

Simply discretizing this equation with the SEM will result in locally conservative, high-order accurate but oscillatory transport scheme. A limiter is added to reduce or eliminate these oscillations ([TStCyrF09]). HOMME supports both monotone and sign-preserving limiters, but the most effective limiter for HOMME has not yet been determined. The default configuration in CAM4 is to use the sign-preserving limiter to prevent negative values of coupled with a sign-preserving hyper-viscosity operator which dissipates .

4.2.14. Conservation and Compatibility

The SEM is compatible, meaning it has a discrete version of the divergence theorem, Stokes theorem and curl/gradient annihilator properties [TF10]. The divergence theorem is the key property of the horizontal discretization that is needed to show conservation. For an arbitrary scalar and vector at layer mid-points, the divergence theorem (or the divergence/gradient adjoint relation) can be written

The discrete version obeyed by the SEM discretization, using , is given by

(88)

The discrete divergence and Stokes theorem apply locally at the element with the addition of an element boundary integral. The local form is used to show local conservation of mass and that the horizontal advection operator locally conserves the two-dimensional potential vorticity (citep{taylor10b}).

In the vertical, [SB81] showed that the and operators needed to satisfy two integral identities to ensure conservation. For any layer interface velocity which satisfies and arbitrary functions of layer mid points. The first identity is the adjoint property (compatibility) for and ,

(89)

which follows directly from the definition of the difference operator given in . The second identity we write in terms of ,

(90)

which is a discrete integrated-by-parts analog of :math:` partial (fg) = f partial g + g partial f. ` Construction of methods with both properties on a staggered unequally spaced grid is the reason behind the complex definition for in .

The energy conservation properties of CAM-HOMME were studied in [Tay10]) using the aqua planet test case ([NH01b][NH01a]). CAM-HOMME uses

as the discretization of the total moist energy . The conservation of is semi-discrete, meaning that the only error in conservation is the time truncation error. In the adiabatic case (with no hyper-viscosity and no limiters), running from a fully spun up initial condition, the error in conservation decreases to machine precision at a second-order rate with decreasing timestep. In the full non-adiabatic case with a realistic timestep, .

The CAM physics conserve a dry energy from [BB03] which is not conserved by the moist primitive equations. Although is small, adiabatic processes in the primitive equations result in a net heating ([Tay10]). If it is desired that the dynamical core conserve instead of , HOMME uses the energy fixer from Energy Fixer.

4.3. Eulerian Dynamical Core

The hybrid vertical coordinate that has been implemented in {cam} is described in this section. The hybrid coordinate was developed by [SStrufing81] 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 [San60], 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 [SStrufing81], who based their development on the generalized vertical coordinate of [Kas74]. A specific form of the coordinate (the hybrid coordinate) is introduced at the latest possible point. The description here differs from [SStrufing81] 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.

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

1. is a monotonic function of .
2. 
3. 
4. where 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 , the last two requirements are identical and the system reduces to that described in [SStrufing81]. The boundary conditions that are required to close the system are:

(91)

(92)

Given the above description of the coordinate, the continuous system of equations can be written following citet{kasahara74} and [SStrufing81]. The prognostic equations are:

(93)

(94)

(95)

(96)

(97)

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

(98)

(99)

(100)

(101)

(102)

(103)

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:

(104)

(105)

(106)

Note that the bounds on the vertical integrals are specified as values of (, 1) or as functions of ( (1), which is the pressure at ).

4.3.2. Conversion to final form

Equations (91) - (106) 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 (95), (96), (98), and (99) may be rewritten as:

(107)

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

(108)

and (97) invoked to specify .

The integrals which appear in (97) , (105) , and (106) can be written more conveniently by expanding the kernel as

(109)

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

(110)

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

(111)

Using (111) as the kernel of the integral in (97), (105), and (106), one obtains integrals of the form

(112)

The original primitive equations (93) -(97), together with (98), (99), and (104) -(106) can now be rewritten with the aid of (107), (108), and (112).

(113)

(114)

(115)

(116)

(117)

(118)

(119)

(120)

(121)

(122)

Once is specified, then can be determined and (113) - (122) can be solved in a GCM.

In the actual definition of the hybrid coordinate, it is not necessary to specify explicitly, since (113) -(122) only requires that and be determined. It is sufficient to specify and to let be defined implicitly. This will be done in section [ssec:finitediffeqs]. In the case that and , (113) - (122) can be reduced to the set of equations solved by CCM1.

4.3.3. Continuous equations using

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 (120). The integral introduces a high order nonlinearity which enters directly into the divergence equation (114). 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 (113) -(122) become:

(123)

(124)

(125)

(126)

(127)

(128)

(129)

(130)

(131)

(132)

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

4.3.4. Semi-implicit formulation

The model described by (123) -(132) , without the horizontal diffusion terms, together with boundary conditions (91) and (92) , is integrated in time using the semi-implicit leapfrog scheme described below. The semi-implicit form of the time differencing will be applied to (124) and (125) without the horizontal diffusion sources, and to (127) . 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 (123) , (124) , and (127) may be rewritten as:

(133)

(134)

(135)

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

(136)

(137)

(138)

Once again, only terms that will be linearized have been explicitly represented in (136) -(138) , and the remaining terms are included in , , and . Anticipating the linearization, and have been replaced by and in (136) and (137). Furthermore, the virtual temperature corrections are included with the other nonlinear terms.

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

(139)

(140)

(141)

In the special case that , (136) - (138) can be converted into equations involving only instead of , and (140) and (141) are not required. This is a major difference between the hybrid coordinate scheme being developed here and the coordinate scheme in CCM1.

Expanding (136) -(138) about the reference state (139) -(141) and retaining only the linear terms explicitly, one obtains:

(142)

(143)

(144)

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

(145)

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 Figure 3.1: A graphical illustration of the different levels of sub-cycling in CAM5.). 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.

Vertical structure of CAM5.0

The finite difference forms of (142) -(144) may then be written down as:

(146)

where denotes a time varying value at time step . The quantities and are defined so as to complete the right-hand sides of (133) -(135). 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.

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

(147)

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

(148)

The first two terms on the right-hand side of (148) 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 (148) 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 (131) , 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 (148) . The derivation will be omitted here, the resulting vertical advection terms are of the form:

(149)

(150)

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 (149) 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 (149) in the system described below.

The last two terms in (148) 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 (130) , (148) can be written as:

(151)

The second term on the right-hand side of (151) 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 (151) can be rewritten as

(152)

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

(153)

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

(154)

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

(155)

The rate of change of total energy due to internal processes is obtained by adding (153) and (155) and must vanish. The first terms on the right-hand side of (153) and (155) 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 (128) and (129) , and in the term of (125) and (132) .

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

(156)

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

(157)

(158)

where for . The quantity is defined to be the unit vector. The inner integral on the right-hand side of (156) is derived from the vertical velocity equation (132) , which we shall approximate as

(159)

where for , and is included as an approximation to for and the symbol is similarly defined as in (150) . will be determined so that is consistent with the discrete continuity equation following [WO94a]. Using (157) and (159) , the finite difference analog of (156) is

(160)

where we have used the relation

(see [3.a.22]). We can now combine the sums in (160) and simplify to give

(161)

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

(162)

Given the definitions of vertical integrals in (158) and (159) and of vertical advection in (149) and (150) the model will conserve energy as long as we require that and satisfy (162) . We are, of course, still neglecting lack of conservation due to the truncation of the horizontal spherical harmonic expansions.

4.3.6. Horizontal diffusion

CAM5.0 contains a horizontal diffusion term for , and to prevent spectral blocking and to provide reasonable kinetic energy spectra. The horizontal diffusion operator in CAM5.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

(163)

(164)

(165)

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 ([Ors74][BMPT77]). The 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 () can be included using the relations:

(166)

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,

(167)

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

(168)

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

(169)

(170)

(171)

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.

4.3.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 , and . Only has the leading term correction to pressure surfaces. Thus, equations that include the terms in this time split sub-step are of the form

(172)

for and , and

(173)

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

(174)

(175)

(176)

for and , and

(177)

(178)

(179)

for , where in the standard model only takes the value 2 in (179) . 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 (174) and (177) . 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 (123) -(132) . 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:

(180)

which gives

(181)

A set of levels may be specified by specifying and , such that , and difference forms of (123) -(132) may be derived.

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

(182)

(183)

(184)

(185)

(186)

(187)

(188)

(189)

(190)

(191)

(192)

(193)

(194)

(195)

(196)

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 :

(197)

(198)

(199)

The matrices and ( with components and ) 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 (183) and (192) at the cost of some loss of generality in the code. The finite difference equations have been written in the form (182) -(199) because this form is quite general. For example, the equations solved by [SStrufing81] at ECMWF can be obtained by changing only the vectors and hydrostatic matrix defined by (196) -(199) .

4.3.8. Time filter

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

4.3.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 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 [DGHS76][BMPT77][Mac79]. 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.

4.3.10. Spectral algorithm overview

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

(200)

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 [Mac79]. The model is coded for a general pentagonal truncation, illustrated in Figure [figure: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:

(201)

The quantity in (200) 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

(202)

With this normalization, the Coriolis parameter is

(203)

which is required for the absolute vorticity.

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

(204)

The inner integral represents a Fourier transform,

(205)

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

(206)

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

(207)

The weights themselves satisfy

(208)

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

(209)

(210)

For the common truncations, these become

(211)

(212)

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

(213)

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 (206) can be rewritten as

(214)

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

(215)

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

(216)

The model uses the spectral transform method (citep{machenhauer79}) 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.

4.3.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 (182) is rewritten

(217)

where the explicit forms of the vectors and are given as

(218)

(219)

(220)

The divergence equation (183) is

(221)

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 (184) is

(222)

The surface–pressure tendency (185) is

(223)

The grouped explicit terms in (221) –(223) are given as follows. The terms of (221) are

(224)

(225)

(226)

(227)

The terms of (222) are

(228)

(229)

(230)

The nonlinear term in (223) is

(231)

4.3.12. Transformation to spectral space

Formally, Equations (217) -(223) are transformed to spectral space by performing the operations indicated in (232) 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 (204) -(206) ,

(232)

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,

(233)

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

(234)

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:

(235)

Defining the derivative of the associated Legendre polynomial by

(236)

(238) can be written

(237)

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

(238)

to each spherical harmonic function individually so that

(239)

where is the Fourier coefficient of the original grid variable .

4.3.13. Solution of semi-implicit equations

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

(240)

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

(241)

The spectral form of the divergence equation (221) becomes

(242)

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

(243)

The spectral thermodynamic equation is

(244)

with defined as

(245)

while the surface pressure equation is

(246)

where is given by

(247)

Equation (240) for vorticity is explicit and complete at this point. However, the remaining equations (242) –(246) are coupled. They are solved by eliminating all variables except :

(248)

(249)

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 (248) then always guarantees . Once is known, and can be computed from (244) and (246) , 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 .

4.3.14. Horizontal diffusion

As mentioned earlier, the horizontal diffusion in (175) and (178) 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:

(250)

(251)

The extra term is present in (250) , (254) and (256) to prevent damping of uniform rotations. The solutions are just

(252)

(253)

(254)

(255)

(256)

(257)

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 (201) , 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 CAM5.0 control runs.

The diffusion of is not complete at this stage. In order to make the partial correction from to in (169) 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.

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

(258)

(259)

4.3.16. Transformation from spectral to physical space

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

(260)

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

(261)

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

(262)

These required derivatives are given by

(263)

(264)

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

(265)

(266)

in which the only nonzero is and

(267)

Thus, the direct transformation is

(268)

(269)

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

(270)

(271)

(272)

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

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

(273)

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

(274)

These heating rates are then combined with the correction,

(275)

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

(276)

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

(277)

4.3.18. Semi-Lagrangian Tracer Transport

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

(278)

(279)

Equation (279) 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

(280)

(281)

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

(282)

(283)

Equation (282) represents the horizontal interpolation of at the departure point calculated assuming . Equation (283) 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

(284)

(285)

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 (284) and (285) is complete, the departure point is given by

(286)

(287)

where the subscript denotes the departure point.

The form given by (284) -(287) 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 [WR89]. 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

(288)

(289)

(290)

(291)

(292)

The calculation of the departure point in the local geodesic system is identical to (284) -(287) 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

(293)

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

(294)

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

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

(295)

where

(296)

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

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

(297)

(298)

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 (282) by Hermite cubic interpolation (294) , with cubic derivative estimates (295) and (296) modified to satisfy the Sufficient Condition for Monotonicity with C continuity (SCMO) described below. Define by

(299)

First, if then

(300)

Then, if either

(301)

or

(302)

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 (282) is followed by the vertical step (283) . The vertical velocity is obtained from that diagnosed in the dynamical calculations (181) by

(303)

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

(304)

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

(305)

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

(306)

with the restriction

(307)

The appropriate values of and are determined by interpolation (294) , with the derivative estimates given by (295) and (296) for to . At the top and bottom we assume a zero derivative (which is consistent with (305) and (307) ), 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 (301) to (302) are applied to the derivative estimates.

4.3.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 [WO94a] are

(308)

(309)

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

(310)

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

(311)

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

(312)

In (311) and (312) the denotes the provisional value before adjustment. The form (312) 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 ([RBB95]). Satisfying (308) and (309) gives

(313)

and

(314)

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:

(315)

The term defines the dry air mass in a layer. Following [RBB95] the change made by the fixer has the same form as (312)

(316)

Substituting (316) into (315) and using (311) through (314) gives

(317)

where the following shorthand notation is adopted:

(318)

We note that there is a small error in (315) . 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 (315) with

(319)

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

4.3.20. Energy Fixer

Following notation in section [massfixers], the total energy integrals are

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

and from (310)

The energy fixer is chosen to have the form

Then

4.3.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 denote a global and vertical average and a horizontal global average. For an arbitrary variable , these are defined by

(320)

(321)

where recall that

(322)

The quantities monitored are:

(323)

4.3.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 ([WR00]). This option provides a saving of computer time of up to 25%.

4.4. Semi-Lagrangian Dynamical Core

4.4.1. Introduction

The two-time-level semi-implicit semi-Lagrangian spectral transform dynamical core in CAM5.0 evolved from the three-time-level CCM2 semi-Lagrangian version detailed in [WO94a] hereafter referred to as W&O94. As a first approximation, to convert from a three-time-level scheme to a two-time-level scheme, the time level index n-1 becomes n, the time level index n becomes n+, and becomes . Terms needed at n+ are extrapolated in time using time n and n-1 terms, except the Coriolis term which is implicit as the average of time n and n+1. This leads to a more complex semi-implicit equation to solve. Additional changes have been made in the scheme to incorporate advances in semi-Lagrangian methods developed since W&O94. In the following, reference is made to changes from the scheme developed in W&O94. The reader is referred to that paper for additional details of the derivation of basic aspects of the semi-Lagrangian approximations. Only the details of the two-time-level approximations are provided here.

4.4.2. Vertical coordinate and hydrostatic equation

The semi-Lagrangian dynamical core adopts the same hybrid vertical coordinate () as the Eulerian core defined by

(324)

where is pressure, is surface pressure, and is a specified constant reference pressure. The coefficients and specify the actual coordinate used. As mentioned by [SB81] and implemented by [SStrufing81] and [SStrufing83], the coefficients and are defined only at the discrete model levels. This has implications in the continuity equation development which follows.

In the system the hydrostatic equation is approximated in a general way by

(325)

where k is the vertical grid index running from 1 at the top of the model to at the first model level above the surface, is the geopotential at level , is the surface geopotential, is the virtual temperature, and R is the gas constant. The matrix , referred to as the hydrostatic matrix, represents the discrete approximation to the hydrostatic integral and is left unspecified for now. It depends on pressure, which varies from horizontal point to point.

4.4.3. Semi-implicit reference state

The semi-implicit equations are linearized about a reference state with constant and . We choose

4.4.4. Perturbation surface pressure prognostic variable

To ameliorate the mountain resonance problem, [RT96] introduce a perturbation surface pressure prognostic variable

The perturbation surface pressure, , is never actually used as a grid point variable in the CAM5.0 code. It is only used for the semi-implicit development and solution. The total is reclaimed in spectral space from the spectral coefficients of immediately after the semi-implicit equations are solved, and transformed back to spectral space along with its derivatives. This is in part because is needed for the horizontal diffusion correction to pressure surfaces. However the semi-Lagrangian CAM5.0 default is to run with no horizontal diffusion.

4.4.5. Extrapolated variables

Variables needed at time () are obtained by extrapolation

4.4.6. Interpolants

Lagrangian polynomial quasi-cubic interpolation is used in the prognostic equations for the dynamical core. Monotonic Hermite quasi-cubic interpolation is used for tracers. Details are provided in the Eulerian Dynamical Core description. The trajectory calculation uses tri-linear interpolation of the wind field.

4.4.7. Continuity Equation

The discrete semi-Lagrangian, semi-implicit continuity equation is obtained from (16) of W&O94 modified to be spatially uncentered by a fraction , and to predict

where

(326)

denotes a vertical difference, denotes the vertical level, denotes the arrival point, the departure point from horizontal (two-dimensional) advection, and the midpoint of that trajectory.

The surface pressure forecast equation is obtained by summing over all levels and is related to (18) of W&O94 but is spatially uncentered and uses

The corresponding equation for the semi-implicit development follows and is related to (19) of W&O94, again spatially uncentered and using .

This is not the actual equation used to determine in the code. The equation actually used in the code to calculate involves only the divergence at time () with eliminated.

The combination

is treated as a unit, and follows from (326) .

4.4.8. Thermodynamic Equation

The thermodynamic equation is obtained from (25) of W&O94 modified to be spatially uncentered and to use . In addition Hortal’s modification (citep{tempertonetal01}) is included, in which

is subtracted from both sides of the temperature equation. This is akin to horizontal diffusion which includes the first order term converting horizontal derivatives from eta to pressure coordinates, with replaced by , and taken as a global average so it is invariant with time and can commute with the differential operators.

Note that represents the heating calculated to advance from time to time and is valid over the interval.

The calculation of follows that of the ECMWF (Research Manual 3, ECMWF Forecast Model, Adiabatic Part, ECMWF Research Department, 2nd edition, 1/88, pp 2.25-2.26) Consider a constant lapse rate atmosphere

4.4.9. Momentum equations

The momentum equations follow from (3) of W&O94 modified to be spatially uncentered, to use , and with the Coriolis term implicit following [CoteS88] and [Tem97]. The semi-implicit, semi-Lagrangian momentum equation at level (but with the level subscript suppressed) is

The gradient of the geopotential is more complex than in the system because the hydrostatic matrix depends on the local pressure:

where is and is the gas constant for water vapor. The gradient of is calculated from the spectral representation and that of from a discrete cubic approximation that is consistent with the interpolation used in the semi-Lagrangian water vapor advection. In general, the elements of are functions of pressure at adjacent discrete model levels

The gradient is then a function of pressure and the pressure gradient

The pressure gradient is available from (324) and the surface pressure gradient calculated from the spectral representation

4.4.10. Development of semi-implicit system equations

The momentum equation can be written as

(327)

where contains known terms at times () and ().

By combining terms, [vmeq] can be written in general as

(328)

where and denote the spherical unit vectors in the longitudinal and latitudinal directions, respectively, at the points indicated by the subscripts, and and denote the appropriate combinations of terms in [vmeq]. Note that is distinct from the . Following [BSHB90], equations for the individual components are obtained by relating the unit vectors at the departure points (,) to those at the arrival points (,):

in which the vertical components () are ignored. The dependence of ’s and ’s on the latitudes and longitudes of the arrival and departure points is given in the Appendix of [BSHB90].

W&O94 followed [BSHB90] which ignored rotating the vector to remain parallel to the earth’s surface during translation. We include that factor by keeping the length of the vector written in terms of the same as the length of the vector written in terms of . Thus, (10) of W&O94 becomes

where

After the momentum equation is written in a common set of unit vectors

Drop the from the notation, define

and transform to vorticity and divergence

Note that

Then the vorticity and divergence equations become

Transform to spectral space as described in the description of the Eulerian spectral transform dynamical core. Note, from (4.5b) and (4.6) on page 177 of [Mac79]

and from (4.5a) on page 177 of [Mac79]

Then the equations for the spectral coefficients at time at each vertical level are

The underbar denotes a vector over vertical levels. Rewrite the vorticity and divergence equations in terms of vectors over vertical levels.

Define by

Then the vorticity and divergence equations are

Note that these equations are uncoupled in the vertical, i.e. each vertical level involves variables at that level only. The equation for however couples all levels.

Define and so that

Let denote the eigenvalues of with corresponding eigenvectors and is the matrix with columns

and the diagonal matrix of corresponding eigenvalues

Then transform

Since is diagonal, all equations are now uncoupled in the vertical.

For each vertical mode, i.e. element of , and for each Fourier wavenumber we have a system of equations in to solve. In following we drop the Fourier index and the modal element index from the notation.

The modal index was included in the above equation on only as a reminder, but will also be dropped in the following.

Substitute and into the equation.

which is just two tri-diagonal systems of equations, one for the even and one for the odd ’s, and

At the end of the system, the boundary conditions are

the term is not present, and from the underlying truncation

For each and we have the general systems of equations

Assume solutions of the form

then

Divergence in physical space is obtained from the vertical mode coefficients by

The remaining variables are obtained in physical space by

4.4.11. Trajectory Calculation

The trajectory calculation follows [Hor99] Let denote the position vector of the parcel,

which can be approximated in general by

Hortal’s method is based on a Taylor’s series expansion

or substituting for

Approximate

for the trajectory equation.

4.4.12. Mass and energy fixers and statistics calculations

The semi-Lagrangian dynamical core applies the same mass and energy fixers and statistical calculations as the Eulerian dynamical core. These are described in sections [massfixers], [energyfixer], and [statcalc].