- 3.3.1 Overview
- 3.3.2 The governing equations for the hydrostatic atmosphere
- 3.3.3 Horizontal discretization of the transport process on the sphere
- 3.3.4 A
*vertically Lagrangian*and*horizontally Eulerian*control-volume discretization of the hydrodynamics - 3.3.5 A mass, momentum, and total energy conserving mapping algorithm
- 3.3.6 Adjustment of specific humidity to conserve water
- 3.3.7 Further discussion

3.3 Finite Volume Dynamical Core

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 Lin and Rood [109] (hereafter LR96) and Lin and Rood [110]
(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 3.3.4). The transport
for tracers, however, can take a much larger time step (*e.g.*,
30 minutes as for the physics).

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*.,
Kasahara [84]). Using standard notations, the hydrostatic balance
equation is given as follows:

where is the density of the air,

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

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

where

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

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

where

3.3.3 Horizontal discretization of the transport process on the sphere

The finite-volume (

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

where , the time-accumulated (from

and

can be interpreted as a time mean (from time to time ) 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
, an approximated subgrid distribution of the true
field, to enable an analytic integration of (3.374). 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. 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 (3.374) 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 (3.374) 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*., Lin et al. [108]) 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 (3.374). The Piecewise Parabolic Method (PPM,
Colella and Woodward [36]) 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. (An extension of the standard
PPM by S.-J. Lin has also been documented in
Machenhauer [119]). 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. A consistent advective-form operator in the direction can be derived from its flux-form counterpart ( as follows:

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 (3.375) 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 ()
direction is derived as follows:*

and likewise the advective-form operator,

where

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

The 2D transport algorithm (

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

where , the mass flux in the meridional direction, is defined in a similar fashion as (3.373). It can be verified that in the special case of constant density flow ( the above equation degenerates to the finite-difference representation of the

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

Note that the above form of the tracer transport equation consistently
degenerates to (3.382) 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.

3.3.4 A

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.364), the conservation law for the pseudo-density (
), becomes

where the symbol represents the vertical difference between the two neighboring Lagrangian surfaces that bound the finite control-volume. From (3.362), the pressure thickness of that control-volume is proportional to the total mass,

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

the thermodynamic equation, (3.366), becomes

and (3.367) and (3.368), the momentum equations, are reduced to

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:

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 (3.391). 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 3.3.5).

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 (3.386) for the pressure thickness . The ideal gas law, the mass conservation law for air mass, the conservation law for the potential temperature (3.388), together with the modified momentum equations (3.389) and (3.390) close the 2D Lagrangian dynamical system, which are vertically coupled only by the hydrostatic relation (see (3.406), section 3.3.5).

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*, (3.380) and
(3.381)):

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:

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 (3.394) and (3.395) are computed as follows:

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 (3.392-3.395), together with their counterparts on the C-grid are cycled times using the fractional time step , which are followed by the tracer transport using (3.387) 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

where the time-accumulated mass fluxes are computed as

The time-averaged winds , defined as follows, are to be used as input for the computations of and

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.

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
Suarez and Takacs [168] for the C-grid Arakawa type dynamical core
(*e.g*., Arakawa and Lamb [5]). 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. However, this would require a significant
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*.,
Ringler et al. [147]). 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. Furthermore, it may be possible to further increase the
size of the small time step via vertical mode decomposition. This
approach is one of the algorithm design issues we plan to revisit.

3.3.5 A mass, momentum, and total energy conserving mapping algorithm

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 (3.391)) as the mapping coordinate. We outline the mapping procedure as follows.

Step 1: Define a suitable Eulerian reference 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 [90] 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 (uandv), and total energy () in the Lagrangian control-volume coordinate based on the Piece-wise Parabolic Method (PPM, Colella and Woodward [36]). The total energy is computed as the sum of the finite-volume integrated geopotential , internal energy , and the kinetic energy () as follows:

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

where is the layer mean temperature, is the 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. 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.

Step 3: Compute kinetic energy in the Eulerian coordinate system for each layer. Substituting kinetic energy and the hydrostatic relationship into (3.404), the layer mean temperature for layer in the Eulerian coordinate is then retrieved from the reconstructed total energy (done in Step 2) by a fully explicit integration procedure starting from the surface up to the model top as follows:

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

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

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 3.3.4, 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 [68], 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.

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. This implies a change of dry air mass at the end of the physics updates. We impose a restriction that dry air mass and water mass be conserved as follows:

The total pressure is

(3.409) |

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

(3.410) |

We define a layer thickness as , so

(3.411) |

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

(3.412) |

or

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

Substituting (3.414) into (3.413),

(3.415) |

(3.416) |

which yields a modified specific humidity for the dynamics:

(3.417) |

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.

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. However, it should be noted that even without the ``energy fixer'' the loss in total energy (in flux unit) is found to be less than 2 ) with the 2 degrees resolution, and much smaller with higher resolution. In the future, we may consider using the total energy as a transported prognostic variable so that the total energy could be automatically conserved.