3.3 Finite Volume Dynamical Core

3.3.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 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 $\eta -coordinate$ (a hybrid $\sigma -p$ 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).

3.3.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 $\zeta$ (e.g., Kasahara [84]). Using standard notations, the hydrostatic balance equation is given as follows:

$$\frac{1}{\rho }\frac{\partial p}{\partial z}+g=0 ,$$ (3.362)

where $\rho$ is the density of the air, p the pressure, and g the gravitational constant. Introducing the ``pseudo-density'' $\pi =\frac{\partial p}{\partial \zeta }$ (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:

$$\pi =-\frac{\partial \Phi }{\partial \zeta }\rho ,$$ (3.363)

where $\Phi =gz$ is the geopotential. Note that $\pi$ reduces to the ``true density'' if $\zeta =-gz$, and the ``surface pressure'' $P_{s}$ if $\zeta =\sigma$ ( $\sigma =\frac{p}{P_{s}}$). The conservation of total air mass using $\pi$ as the prognostic variable can be written as

$$\frac{\partial }{\partial t}\pi +\nabla \cdot \left(\overrightarrow{V}\pi \right) =0 ,$$ (3.364)

where $\overrightarrow{V}=(u,v,\frac{d\zeta }{dt})$. Similarly, the mass conservation law for tracer species (or water vapor) can be written as

$$\frac{\partial }{\partial t}(\pi q)+\nabla \cdot \left(\overrightarrow{V}\pi q\right) =0 ,$$ (3.365)

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

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

$$\frac{\partial }{\partial t}(\pi \Theta )+\nabla \cdot \left( \overrightarrow{V}\pi \Theta \right) =0 .$$ (3.366)

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

$$\frac{\partial }{\partial t}u=\Omega v-\frac{1}{Acos\theta }\left... ...rtial \lambda }p\right] -\frac{d\zeta }{dt}\frac{\partial u}{\partial \zeta } ,$$ (3.367)

$$\frac{\partial }{\partial t}v=-\Omega u-\frac{1}{A}\left[ \frac{\... ...artial \theta }p\right] -\frac{d\zeta }{dt}\frac{\partial v}{\partial \zeta } ,$$ (3.368)

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

$$D=\frac{1}{Acos\theta }\left[ \frac{\partial }{\partial \lambda }(u)+\frac{\partial }{\partial \theta }(v cos\theta )\right] ,$$

$$\kappa =\frac{1}{2}\left( u^{2}+v^{2}\right) ,$$

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

$$\Omega =2\omega sin\theta +\frac{1}{A cos\theta }\left[ \fra... ...\partial \lambda }v-\frac{\partial }{\partial \theta }(u cos\theta )\right] ,$$

where $\omega$ is the angular velocity of the earth. Note that the last term in (3.367) and (3.368) vanishes if the vertical coordinate $\zeta$ is a conservative quantity (e.g., entropy under adiabatic conditions [78] or an imaginary conservative tracer), and the 3D divergence operator becomes 2D along constant $\zeta$ surfaces. The discretization of the 2D horizontal transport process is described in section 3.3.3. The complete dynamical system using the Lagrangian control-volume vertical discretization is described in section 3.3.4. A mass, momentum, and total energy conservative mapping algorithm is described in section 3.3.5.

3.3.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 3.3.4), we shall present here only the 2D (horizontal) forms of the FFSL transport algorithm for the transport of density (3.364) and mixing ratio-like quantities (3.365) on the sphere. The governing equation for the pseudo-density (3.364) becomes

$$\frac{\partial }{\partial t}\pi +\frac{1}{Acos\theta }\left[ \fra... ...bda }(u\pi )+\frac{\partial }{\partial \theta }(v\pi cos\theta )\right] =0 .$$ (3.369)

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

$$\widetilde{\pi }(t)\equiv \frac{1}{A^{2}\Delta \theta \Delta \lambda cos\theta }\iint \pi (t;\lambda ,\theta )A^{2}cos\theta d\theta d\lambda .$$ (3.370)

Given the exact 2D wind field $\overrightarrow{V}(t;\lambda ,\theta )=(U,V)$ the 2D integral representation of the conservation law for $\widetilde{\pi }$ can be obtained by integrating (3.369) in time and in space

$$\widetilde{\pi }^{n+1}=\widetilde{\pi }^{n}-\frac{1}{A^{2}\Delta ... ...(t;\lambda ,\theta )\overrightarrow{V}\cdot \overrightarrow{n}dl\right] dt .$$ (3.371)

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

$$\delta _{x}q=q(x+\frac{\Delta x}{2})-q(x-\frac{\Delta x}{2}),$$

and assuming $(u^{*},v^{*})$ is the time-averaged (from time $t$ to time $t + \Delta t$) $\overrightarrow{V}$ on the C-grid (e.g., Fig. 1 in LR96), the 1-D finite-volume flux-form transport operator F in the $\lambda$-direction is

$$F(u^{*},\Delta t,\widetilde{\pi })=-\frac{1}{A\Delta \lambda cos\... ...mbda cos\theta } \delta _{\lambda }\left[ \chi (u^{*},\Delta t;\pi )\right] ,$$ (3.372)

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

$$\chi (u^{*},\Delta t;\pi )=\frac{1}{\Delta t}\int _{t}^{t+\Delta t}\pi U dt\equiv u^{*}\pi ^{*}(u^{*},\Delta t,\widetilde{\pi }) ,$$ (3.373)

and

$$\pi ^{*}(u^{*},\Delta t;\widetilde{\pi })\approx \frac{1}{\Delta t}\int _{t}^{t+\Delta t}\pi dt$$ (3.374)

can be interpreted as a time mean (from time $t$ to time $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 $t=t+\Delta t$ (the arrival point; (e.g., point B in Fig. 3 of LR96) back to time $t$ (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 $\widetilde{\pi }$, an approximated subgrid distribution of the true $\pi$ field, to enable an analytic integration of (3.374). Assuming there is no error in obtaining the time-mean wind $(u^{*})$, the only error produced by the 1D transport scheme would be solely due to the approximation to the continuous distribution of $\pi$ 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 $\pi ^{*}$. 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 $\lambda -$direction can be derived from its flux-form counterpart ($F)$ as follows:

$$f(u^{*},\Delta t,\widetilde{\pi })=F(u^{*},\Delta t,\widetilde{\p... ...v 1)=F(u^{*},\Delta t,\widetilde{\pi })+\widetilde{\pi } C_{def}^{\lambda } ,$$ (3.375)

$$C^{\lambda }_{def}=\frac{\Delta t \delta _{\lambda }u^{*}}{A\Delta \lambda cos\theta } ,$$ (3.376)

where $C_{def}^{\lambda }$ is a dimensionless number indicating the degree of the flow deformation in the $\lambda$-direction. The above derivation of $f$ 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 ($\theta$) direction is derived as follows:

$$G(v^{*},\Delta t,\widetilde{\pi })=-\frac{1}{A\Delta \theta cos\t... ...theta cos\theta } \delta _{\theta }\left[ v^{*}cos\theta \pi ^{*}\right] ,$$ (3.377)

and likewise the advective-form operator,

$$g(v^{*},\Delta t,\widetilde{\pi })=G(v^{*},\Delta t,\widetilde{\pi })+\widetilde{\pi } C_{def}^{\theta } ,$$ (3.378)

where

$$C^{\theta }_{def}=\frac{\Delta t \delta _{\theta }\left[ v^{*}cos\theta \right] }{A\Delta \theta cos\theta } .$$ (3.379)

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

$$( )^{\theta }=( )^{n}+\frac{1}{2}g\left[ v^{*},\Delta t, ( )^{n}\right] ,$$ (3.380)

$$( )^{\lambda }=( )^{n}+\frac{1}{2}f\left[ u^{*},\Delta t, ( )^{n}\right] .$$ (3.381)

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

$$\widetilde{\pi }^{n+1}=\widetilde{\pi }^{n}+F\left[ u^{*},\Delta ... ...}^{\theta }\right] +G\left[ v^{*},\Delta t,\widetilde{\pi }^{\lambda }\right] .$$ (3.382)

Using explicitly the mass fluxes $\left( \chi ,Y\right)$, (3.382) is rewritten as

$$\widetilde{\pi }^{n+1}=\widetilde{\pi }^{n}-\frac{\Delta t}{Acos\... ...t[ cos\theta Y(v^{*},\Delta t;\widetilde{\pi }^{\lambda })\right] \right\} ,$$ (3.383)

where $Y$, the mass flux in the meridional direction, is defined in a similar fashion as $\chi$ (3.373). It can be verified that in the special case of constant density flow ( $\widetilde{\pi }=constant)$ the above equation degenerates to the finite-difference representation of the incompressibility condition of the ``time mean'' wind field $(u^{*},v^{*})$, i.e.,

$$\frac{1}{\Delta \lambda }\delta _{\lambda }u^{*}+\frac{1}{\Delta \theta }\delta _{\theta }\left( v^{*}cos\theta \right) =0 .$$ (3.384)

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 $(\widetilde{q})$ the mass fluxes $(\chi ,Y)$ as defined previously should be used as follows

$$\widetilde{q}^{n+1}=\frac{1}{\widetilde{\pi }^{n+1}}\left[ \widet... ...elta t,\widetilde{q}^{\theta })+G(Y,\Delta t,\widetilde{q}^{\lambda })\right] .$$ (3.385)

Note that the above form of the tracer transport equation consistently degenerates to (3.382) if $\widetilde{q}\equiv 1$ (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 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. Starr [164]) 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 $\sigma$ or hybrid $\sigma$-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 $\zeta$ that is a monotonic function of height and constant on the initial reference coordinate surfaces (e.g., the value of ``$\eta$'' in the hybrid $\sigma -p$ 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 $\delta \zeta$, (3.364), the conservation law for the pseudo-density ( $\pi =\frac{\delta p}{\delta \zeta }$), becomes

$$\frac{\partial }{\partial t}\delta p+\frac{1}{Acos\theta }\left[ ... ...delta p)+\frac{\partial }{\partial \theta }(v\delta p cos\theta )\right] =0 ,$$ (3.386)

where the symbol $\delta$ represents the vertical difference between the two neighboring Lagrangian surfaces that bound the finite control-volume. From (3.362), the pressure thickness $\delta p$ of that control-volume is proportional to the total mass, i.e., $\delta p=-\rho g\delta z$. Therefore, it can be said that the Lagrangian control-volume vertical discretization has the hydrostatic balance built-in, and $\delta p$ can be regarded as the ``pseudo-density'' for the discretized Lagrangian vertical coordinate system.

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

$$\frac{\partial }{\partial t}(q\delta p)+\frac{1}{Acos\theta }\lef... ...delta p)+\frac{\partial }{\partial \theta }(vq\delta p cos\theta )\right] =0,$$ (3.387)

the thermodynamic equation, (3.366), becomes

$$\frac{\partial }{\partial t}(\Theta \delta p)+\frac{1}{Acos\theta... ...p)+\frac{\partial }{\partial \theta }(v\Theta \delta p cos\theta )\right] =0,$$ (3.388)

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

$$\frac{\partial }{\partial t}u=\Omega v-\frac{1}{Acos\theta }\left... ...Phi -\nu D\right) +\frac{1}{\rho }\frac{\partial }{\partial \lambda }p\right] ,$$ (3.389)

$$\frac{\partial }{\partial t}v=-\Omega u-\frac{1}{A}\left[ \frac{\... ...\Phi -\nu D\right) +\frac{1}{\rho }\frac{\partial }{\partial \theta }p\right] .$$ (3.390)

Given the prescribed pressure at the model top $P_{\infty }$, the position of each Lagrangian surface $P_{l}$ (horizontal subscripts omitted) is determined in terms of the hydrostatic pressure as follows:

$$P_{l}=P_{\infty }+\sum ^{l}_{k=1}\delta P_{k}, (for l=1, 2, 3, ..., N) ,$$ (3.391)

where the subscript $l$ is the vertical index ranging from 1 at the lower bounding Lagrangian surface of the first (the highest) layer to $N$ at the Earth's surface. There are $N$+1 Lagrangian surfaces to define a total number of $N$ Lagrangian layers. The surface pressure, which is the pressure at the lowest Lagrangian surface, is easily computed as $P_{N}$ 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 $h$ in the shallow water system of LR97 is replaced by (3.386) for the pressure thickness $\delta p$. 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 ($\Delta t)$ 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, $\Delta \tau =\Delta t/m$, where $m$ 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 $\left[ \delta p,\Theta ,u,v;q\right]$ on the D-grid. Discretization on the C-grid for obtaining the diagnostic variables, the time-averaged winds $(u^{*},v^{*})$, is analogous to that of the D-grid (see also LR97).

Introducing the following short hand notations (cf, (3.380) and (3.381)):

$$( )_{i}^{\theta }=( )^{n+\frac{i-1}{m}}+\frac{1}{2}g[v_{i}^{*},\Delta \tau ,( )^{n+\frac{i-1}{m}}],$$

$$( )_{i}^{\lambda }=( )^{n+\frac{i-1}{m}}+\frac{1}{2}f[u_{i}^{*},\Delta \tau ,( )^{n+\frac{i-1}{m}}],$$

and applying directly (3.383), the update of ``pressure thickness'' $\delta p$, using the fractional time step $\Delta \tau =\Delta t/m$, can be written as

$$\delta p^{n+\frac{i}{m}}=\delta p^{n+\frac{i-1}{m}}-\frac{\Delta ... ...ta y_{i}^{*}(v_{i}^{*},\Delta \tau ;\delta p_{i}^{\lambda })\right] \right\}$$ (3.392)

$$(for i=1,...,m),$$

where $\left[ x_{i}^{*},y_{i}^{*}\right]$ are the background air mass fluxes, which are then used as input to Eq. 24 for transport of the potential temperature $\Theta$:

$$\Theta ^{n+\frac{i}{m}}=\frac{1}{\delta p^{n+\frac{i}{m}}}\left[ ... ...Theta _{i}^{\theta })+G(y_{i}^{*},\Delta \tau ,\Theta _{i}^{\lambda })\right] .$$ (3.393)

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:

$$u^{n+\frac{i}{m}}=u^{n+\frac{i-1}{m}}+\Delta \tau \left[ y_{i}... ...heta }\delta _{\lambda }(\kappa ^{*}-\nu D^{*})+\widehat{P_{\lambda }}\right] ,$$ (3.394)

$$v^{n+\frac{i}{m}}=v^{n+\frac{i-1}{m}}-\Delta \tau \left[ x_{i}... ...\theta }\delta _{\theta }(\kappa ^{*}-\nu D^{*})-\widehat{P_{\theta }}\right] ,$$ (3.395)

where $\kappa ^{*}$ is the upwind-biased ``kinetic energy'' (as defined by Eq. 18 in LR97), and $D^{*}$, the horizontal divergence on the D-grid, is discretized as follows:

$$D^{*}=\frac{1}{Acos\theta }\left[ \frac{1}{\Delta \lambda }\delta... ...\theta }\delta _{\theta }\left( v^{n+\frac{i-1}{m}}cos\theta \right) \right] .$$

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

$$\widehat{P_{\lambda }}=\frac{\oint _{\Pi \rightleftharpoons \lamb... ...i d\Pi }{Acos\theta \oint _{\Pi \rightleftharpoons \lambda }\Pi d\lambda } ,$$ (3.396)

$$\widehat{P_{\theta }}=\frac{\oint _{\Pi \rightleftharpoons \theta }\phi d\Pi }{A \oint _{\Pi \rightleftharpoons \theta }\Pi d\theta } ,$$ (3.397)

where $\Pi =p^{\kappa } (\kappa =R/C_{p})$, and the symbols `` $\Pi \rightleftharpoons \lambda$'' and ``$\Pi$ $\rightleftharpoons \theta$'' indicate that the contour integrations are to be carried out, using the finite-volume algorithm described in L97, in the $(\Pi ,\lambda )$ and $(\Pi ,\theta )$ space, respectively.

To complete one time step, equations (3.392-3.395), together with their counterparts on the C-grid are cycled $m$ times using the fractional time step $\Delta \tau$, which are followed by the tracer transport using (3.387) with the large-time-step $\Delta t$.

Mass fluxes $(x^{*},y^{*})$ and the winds $(u^{*},v^{*})$ on the C-grid are accumulated for the large-time-step transport of tracer species (including water vapor) $q$ as

$$q_{}^{n+1}=\frac{1}{\delta p^{n+1}}\left[ q_{}^{n}\delta p^{n}+F(X^{*},\Delta t,q_{}^{\theta })+G(Y^{*},\Delta t,q_{}^{\lambda })\right] ,$$ (3.398)

where the time-accumulated mass fluxes $(X^{*},Y^{*})$ are computed as

$$X^{*}=\sum ^{m}_{i=1}x_{i}^{*}(u_{i}^{*}, \Delta \tau , \delta p_{i}^{\theta }) ,$$ (3.399)

$$Y^{*}=\sum _{i=1}^{m}y_{i}^{*}(v_{i}^{*}, \Delta \tau , \delta p_{i}^{\lambda }) .$$ (3.400)

The time-averaged winds $(U^{*},V^{*})$, defined as follows, are to be used as input for the computations of $q^{\lambda }$ and $q^{\theta }:$

$$U^{*}=\frac{1}{m}\sum ^{m}_{i=1}u_{i}^{*} ,$$ (3.401)

$$V^{*}=\frac{1}{m}\sum ^{m}_{i=1}v_{i}^{*} .$$ (3.402)

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 $\Delta \tau$. 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

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 $\delta p$), 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 ($\delta p$) 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 $\eta -coordinate$ 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 ($q$), zonal and meridional winds (u and v), and total energy ($\Gamma$) in the Lagrangian control-volume coordinate based on the Piece-wise Parabolic Method (PPM, Colella and Woodward [36]). The total energy $\Gamma$ is computed as the sum of the finite-volume integrated geopotential $\phi$, internal energy $(C_{v}T)$, and the kinetic energy ($K$) as follows:

$$\Gamma =\frac{1}{\delta p}\int \left[ C_{v}T+\phi +\frac{1}{2}\left( u^{2}+v^{2}\right) \right] dp .$$ (3.403)

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

$$\Gamma =C_{p}\overline{T}+\frac{1}{\delta p}\delta \left( p\phi \right) +K ,$$ (3.404)

where $\overline{T}$ is the layer mean temperature, $K$ is the kinetic energy, $p$ is the pressure at layer edges, and $C_{v}$ and $C_{p}$ are the specific heat of the air at constant volume and at constant pressure, respectively. Layer mean values of $q$, (u, v), and $\Gamma$ 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 $\overline{T}_{k}$ for layer $k$ 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:

$$\overline{T}_{k}=\frac{\Gamma _{k}-K_{k}-\phi _{k+\frac{1}{2}}}{C... ...{1}{2}}-ln p_{k-\frac{1}{2}}}{p_{k+\frac{1}{2}}-p_{k-\frac{1}{2}}}\right] } .$$ (3.405)

To convert the potential temperature $\Theta$ to the layer mean temperature the conversion factor is obtained by equating the following two equivalent forms of the hydrostatic relation for $\Theta$ and $\overline{T}:$

$$\delta \phi =-C_{p}\Theta \delta \Pi ,$$ (3.406)

$$\delta \phi =-R\overline{T} \delta ln \Pi ,$$ (3.407)

where $\Pi =p^{\kappa }$. The conversion formula between layer mean temperature and layer mean potential temperature is obtained as follows:

$$\Theta =\kappa \frac{\delta ln\Pi }{\delta \Pi }\overline{T} .$$ (3.408)

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

3.3.6 Adjustment of specific humidity to conserve water

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 $p$ is

$$p = d + e .$$ (3.409)

with dry pressure $d$, water vapor pressure $e$. The specific humidity is

$$q = \frac{e}{p} = \frac{e}{d+e}, \qquad d = (1-q) p .$$ (3.410)

We define a layer thickness as $\delta^k p \equiv p^{k+1/2} - p^{k-1/2}$, so

$$\delta^k d = (1-q^k)\delta^k p .$$ (3.411)

We are concerned about 3 time levels: $q_n$ is input to physics, $q_{n*}$ is output from physics, $q_{n+1}$ is the adjusted value for dynamics.

Dry mass is the same at $n$ and $n+1$ but not at $n*$. To conserve dry mass, we require that

$$\delta^k d_n = \delta^k d_{n+1}$$ (3.412)

or

$$(1-q^k_n)\delta^k p_n = (1-q^k_{n+1})\delta^k p_{n+1} .$$ (3.413)

Water mass is the same at $n*$ and $n+1$, but not at $n$. To conserve water mass, we require that

$$q^k_{n*} \delta^k p_n = q^k_{n+1}\delta^k p_{n+1}.$$ (3.414)

Substituting (3.414) into (3.413),

$$(1-q^k_n)\delta^k p_n = \delta^k p_{n+1} - q^k_{n*} \delta^k p_n$$ (3.415)

$$\delta^k p_{n+1} = (1 - q^k_n + q^k_{n*})\delta^k p_n$$ (3.416)

which yields a modified specific humidity for the dynamics:

$$q^k_{n+1} = q^k_n \frac{\delta^k p_n}{\delta^k p_{n+1}} = \frac{q^k_{n*}}{1 - q^k_n + q^k_{n*}} .$$ (3.417)

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

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 $(W/m^{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.