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  (hereafter LR96) and Lin and Rood  (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 ). Using standard notations, the hydrostatic balance equation is given as follows:
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:
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
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. ) 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 ) 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 ). 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:
Analogously, the 1D flux-form transport operator G in the latitudinal () direction is derived as follows:
To complete the construction of the 2D algorithm on the sphere, we introduce the following short hand notations:
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.
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
Similarly, (3.365), the mass conservation law for all tracer species, is
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:
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:
The finite-volume mean pressure-gradient terms in (3.394) and (3.395) are computed as follows:
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
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  for the C-grid Arakawa type dynamical core (e.g., Arakawa and Lamb ). 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. ). 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.
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  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 based on the Piece-wise Parabolic Method (PPM, Colella and Woodward ). 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
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 , 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
Dry mass is the same at and but not at . To conserve dry mass, we require that
Water mass is the same at and , but not at . To conserve water mass, we require that
Substituting (3.414) into (3.413),
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.