Two flaws in the semi-Lagrangian algorithm originally implemented as an optional dynamical core in the NCAR Community Atmosphere Model (CAM3.1) are exposed by steady-state and baroclinic instability test cases. Remedies are demonstrated and have been incorporated in the dynamical core. One consequence of the first flaw is an erroneous damping of the speed of a zonally uniform zonal wind undergoing advection by a zonally uniform zonal flow field. It results from projecting the transported vector wind expressed in unit vectors at the arrival point to the surface of the sphere and is eliminated by rotating the vector to be parallel to the surface. The second flaw is the formulation of an a posteriori energy fixer that, although small, systematically affects the temperature field and leads to an incorrect evolution of the growing baroclinic wave. That fixer restores the total energy at each time step by changing the provisional forecast temperature proportionally to the magnitude of the temperature change at that time step. Two other fixers are introduced that do not exhibit the flaw. One changes the provisional temperature everywhere by an additive constant, and the other changes it proportionally by a multiplicative constant.
Testing the dynamical core of an atmospheric general circulation model (GCM) in isolation from the physical parameterization package is an essential stepping stone during the model development and evaluation phase. The dynamical core tests have the potential to assess not only the diffusive characteristics and convergence of the numerical schemes but also to reveal flaws in the dynamical core formulation. Such flaws might not be obvious in weather or climate simulations if physical parameterizations or model components, like the ocean or ice, compensate or mask the inaccuracies. The term dynamical core generally refers to the numerical approximation to the resolved fluid flow component of an atmospheric model. The hydrostatic primitive, nonhydrostatic shallow-atmosphere, and nonhydrostatic deep-atmosphere equation sets are most commonly used. The discussion here concentrates on the assessment of a hydrostatic dynamical core. In particular, we discuss two flaws in an early version of the semi-Lagrangian spectral transform dynamical core that is part of the Community Atmosphere Model (CAM3.1) developed at the National Center for Atmospheric Research (NCAR).
The flaws have been detected with the help of two recently developed steady-state and baroclinic wave test cases (Jablonowski and Williamson 2006a,b) and have since been eliminated in the standard release of CAM3.1 (Collins et al. 2004, 2006). The first flaw resulted in a damping of the speed of a zonally uniform zonal wind undergoing advection by a zonally uniform zonal flow field. As will be illustrated in section 3a, it was caused by the method used to obtain the vector wind at the arrival point, which systematically shortened the vector at every time step. The second flaw was the formulation of an a posteriori total energy fixer that, although small, systematically affected the temperature field and led to an incorrect evolution of the growing baroclinic wave. This paper serves as a precaution for other modeling groups using similar algorithms. It is organized as follows: section 2 briefly reviews the design and application of the steady-state and baroclinic wave test cases, section 3 discusses the flaw in the vector wind and a posteriori energy fixer, and section 4 provides a summary.
2. Review of the dynamical core test cases
The idealized steady-state and baroclinic instability test cases were suggested by Jablonowski (2004) and were applied to four very different dynamical cores by Jablonowski and Williamson (2006a,b). Among them were the three dynamical cores that are part of the NCAR CAM3.1 modeling framework. These are the spectral transform Eulerian (EUL) and semi-Lagrangian (SLD) models (Collins et al. 2004, 2006) as well as the finite-volume (FV) dynamical core developed by Lin (2004). In addition, the icosahedral finite-difference global model (GME) of the German Weather Service (Deutscher Wetterdienst) was tested (Majewski et al. 2002). These hydrostatic dynamical cores represent a broad range of numerical approaches and, at very high resolutions, provide independent reference solutions for the baroclinic instability test case. These are used to assess whether new model runs fall within the uncertainty range of the reference solutions. This is the primary measure for detecting flaws in the model formulation. In addition, the steady-state model runs are compared with the initial state, which is an analytic solution to the primitive equations.
The steady-state and baroclinic wave test cases have been developed for dry dynamical cores in spherical geometry with pressure-based vertical coordinates. The test cases are based on an analytically defined zonally symmetric zonal flow with two jets in the midlatitudes. In addition, a smooth surface geopotential, a constant surface pressure, and a quasi-realistic meridional and vertical temperature distribution are prescribed. The meridional wind is set to zero. All analytical expressions and further details can be found in Jablonowski and Williamson (2006a, hereinafter JW06a). These initial conditions are a balanced steady-state solution of the primitive equations and satisfy static, inertial, and symmetric stability properties. However, they are unstable with respect to baroclinic or barotropic instability mechanisms. A baroclinic wave can be triggered if the initial conditions are overlaid with a perturbation. In particular, we select a small-amplitude but wide Gaussian hill perturbation that overlays the zonal wind field at all vertical levels. This triggers the evolution of a baroclinic wave over the course of 10 days, with explosive cyclogenesis around day 8.
3. Exposed flaws in the semi-Lagrangian dynamical core
a. Semi-Lagrangian advection of the vector wind
The SLD dynamical core is initialized with the balanced initial conditions without a perturbation. This evaluates the model’s ability to maintain the steady-state solution and serves as a stringent debugging tool. The model is run for 30 days as suggested by Jablonowski and Williamson (2006a). However, as will be seen below, 10-day integrations are adequate to expose the errors discussed here, making it a very inexpensive test case. The resolution employed is T42L26 (triangular truncation) with varying time steps. This resolution corresponds to a horizontal grid spacing of approximately 2.8° × 2.8° in the underlying Gaussian grid and utilizes 26 hybrid levels in the vertical direction. The details of the vertical level placement can be found in Jablonowski and Williamson (2006b). No explicit diffusion is used. For these model runs error norms can be directly computed because the analytical solution is known.
We observed that the magnitude of the zonally uniform zonal wind in the SLD run decreased with time whereas both the EUL and FV dynamical cores maintained the zonal wind magnitude with very minor variations (JW06a). Figure 1 shows the time series of the global root-mean-square l2 error norms of the zonally averaged zonal wind for the SLD T42L26 simulations. The definition of the l2 error norm is given in JW06a’s Eq. (15). The figure depicts the global l2 error norms for the three time steps Δt = 900 s, Δt = 1800 s, and Δt = 3600 s. Both SLD versions with and without the flaw in the vector wind are shown. In the flawed simulations, the error is greater the longer the time step is. This contradicts the general perception that the damping arising from semi-Lagrangian interpolation is smaller with longer time steps because fewer interpolations are performed for a fixed elapsed time. The error in Fig. 1 is not caused by the damping from the interpolation. Because the flow is zonal and zonally uniform, the departure point is on the same latitude circle as the arrival point and the interpolation is between points in longitude that have the same value. The interpolation scheme should maintain a constant field in that case.
The damping is a result of the formulation of the vector momentum equation. In the original CAM3.1 formulation, Williamson and Olson (1994) developed the momentum equation approximations in vector form following Ritchie (1988, 1991) and Bates et al. (1990, 1993). In general, the vector wind that is parallel to the spherical surface at the departure point is no longer parallel when transported to the arrival point. This is schematically drawn in Fig. 2, which shows the wind vector at the departure (label 1) and arrival points (labels 2–4). The wind vectors at the arrival point depict the transported vector (label 2) that can either be projected onto the spherical surface (label 3) or rotated to be parallel to the surface (label 4). To calculate the vector components needed in a practical implementation, the original CAM3.1 formulation followed Bates et al. (1990) in which the equations for the individual vector components are obtained by relating the unit vectors at the departure points to those at the arrival points, ignoring the vertical components. The general forms of the equations are
with γ = 1. Here 𝒰 and 𝒱 denote combinations of terms and subscripts A and D denote the arrival and departure points, respectively. The α and β coefficients relate the unit vectors at the departure point to those at the arrival point. [See section 3a of Williamson and Olson (1994) for details.] Equation (1) with γ = 1 essentially projects the transported vector onto the spherical surface and thereby shortens it. Thus, the zonal wind decreases with time in this test case. With the longer time step, the wind at the arrival point is less parallel to the surface and is thus shortened more by the projection to the surface. In contrast, Ritchie (1988, 1991) rotated the vector to be parallel to the spherical surface, maintaining its length. In his case
which simply rescales the vector at the arrival point to have the same length as that at the departure point. Solutions with this rotated approach are also shown in Fig. 1 for the same three time steps. All errors from the rotated vector cases are clearly reduced and no longer exhibit the dependence on the time step length. They graphically overlay each other and now match those from the Eulerian dynamical core seen in Fig. 4a of JW06a. This revised formulation was implemented into CAM3.1 [Eq. (3.298) of Collins et al. (2004)] and used in the experiments reported in JW06a. We note that the original formulation always shortens the advected vector. The effect, however, is generally not seen in complete model integrations or in more complicated test cases such as Held and Suarez (1994). Other factors such as the parameterized forcing and the damping from the interpolations dominate the error. We emphasize that the error is dependent only on the time step and not on the vertical or horizontal grid size. The semi-Lagrangian curves presented here with the revised formulation are not the same as those presented in Fig. 4b of JW06a. Those errors were calculated with the model configured as it would be for an actual GCM climate simulation or weather forecast. Here we configured the SLD model to isolate the vector wind problem. In particular, we did not revert to a geodesic trajectory calculation near the poles (Williamson and Rasch 1989) because the flow is zonal. We also did not include a decentering mechanism to minimize orographic resonance (Collins et al. 2004). Both of those approximations contribute to the l2 error growth seen in Fig. 4a of JW06a, with the decentering mechanism being the dominant source.
b. Total energy fixer
The second test imposes a well-resolved Gaussian-hill perturbation on the balanced zonal flow. This perturbation grows with time and evolves into a baroclinic wave. We run the test for 30 days with both the spectral transform Eulerian and semi-Lagrangian dynamical cores at T170 triangular truncation and 26 hybrid levels. This resolution corresponds to a horizontal grid spacing of approximately 0.7° × 0.7° in the underlying Gaussian grid. The time steps are Δt = 300 s for EUL and Δt = 900 s for SLD. No explicit diffusion is used in the SLD simulations, whereas in EUL second- and fourth-order diffusion operators are applied to the prognostic variables (JW06a).
When we first ran this test with the SLD version we obtained unsatisfactory results. We isolated the source of the error in the semi-Lagrangian integration to an a posteriori total energy fixer that was included in that dynamical core. Boville (2000) and Williamson (2007) explain the desirability of including an energy fixer in the atmospheric component of a coupled climate system model. After identifying the source of the error, we considered two additional variants of an a posteriori total energy fixer. The details of these three energy fixers are explained below. In short, they artificially restore the global conservation of total energy in non-energy-conserving simulations. This is done by adjusting the temperature field after each time step. We label these different adjustment techniques simply as FIXER 1–3, FIXER 1 being the original, problematic fixer. Figure 3 shows the root-mean-square l2 difference of surface pressure between four semi-Lagrangian integrations and the Eulerian run. The four are with the three fixers and an additional run with no fixer. The definition of the l2 difference can be found in JW06a’s Eq. (16). The label “FIXER 1” (short dashed line) denotes the original implementation. Note that FIXER 1 lies outside the stippled region, which, at this resolution, indicates a flaw in the SLD model simulation. As argued in JW06a, the stippled region defines the uncertainty in the baroclinic wave reference solutions obtained from a set of four very different dynamical cores applied at high horizontal resolution. One of these reference solutions is that from an SLD model run with an improved fixer (FIXER 2). Differences between the reference solutions with different dynamical cores are introduced by different truncation errors. Although initially these errors might be very small at high resolutions, they do grow with time because the basic state is baroclinically unstable. This error growth limits how well the reference solutions are actually known. The stippled region in Fig. 3 is bounded by the envelope of differences between all combinations of the high-resolution solutions (JW06a). If the differences of one model against one of those reference solutions fall within the stippled region then the solution is captured to within the uncertainty. Differences outside the stippled region indicate that the solution is not as accurate as any of those reference solutions. This analysis technique is also explained in more detail in JW06a.
It is clear that the original semi-Lagrangian implementation with FIXER 1 does not produce an accurate evolution of the baroclinic wave. This can further be seen in Fig. 4, which plots the surface pressure at day 10 from three T170L26 SLD solutions and the T170L26 EUL run. Even the resolved scales, well away from the truncation limit, are incorrect in the semi-Lagrangian FIXER 1 solution. Plots for T340L26 SLD with FIXER 1 (not shown) look very similar to those of T170L26 SLD with FIXER 1, indicating that the flawed SLD simulation is converging to a different solution in comparison with models like EUL, FV, GME, and SLD (FIXER 2) at equivalent resolutions. The latter are shown at day 10 in Jablonowski and Williamson (2006b).
In the continuous primitive equations, total energy is conserved if
(Laprise and Girard 1990). This equation is valid in the absence of diabatic and frictional effects for hybrid vertical coordinate systems as used here. Variable pairs Φs and ps and Φtop and ptop are the geopotential and pressure at the surface and the model top, respectively; cp is the specific heat of dry air at constant pressure; and v stands for the horizontal velocity vector with zonal and meridional wind components u and υ. Furthermore, T symbolizes the temperature, p is the pressure, and t denotes the time. The integrals span the 3D and 2D domains where A symbolizes the horizontal area of the sphere and η stands for the hybrid vertical coordinate. The vertical integral is bounded by the value ηs at the surface and ηtop at the model top. Here, ηs is identical to unity and ηtop is set to ptop/p0 with p0 = 1000 hPa. Note that ηtop is nonzero for constant ptop > 0 hPa, which is the case here. A constant pressure at the model top ensures the global conservation of total energy in the continuous equations and simplifies the 2D integral. Equation (3) then becomes
where E symbolizes the global integral of the total energy as shown by the term in the angle bracket in Eq. (4). In the semidiscrete system with ∂p/∂η ≈ Δp/Δη and dη ≈ Δη, the domain-integrated total energy E is given by
The summation index k indicates the vertical index of a full model level with the maximum level number K near the surface. The pressure difference Δpk is defined as
where ΔAk = Ak+1/2 − Ak−1/2 and ΔBk = Bk+1/2 − Bk−1/2. The discrete positions of the hybrid coefficients Ak+1/2 and Bk+1/2 at the model interface levels are listed in Jablonowski and Williamson (2006b); Δηk is given by Δηk = ηk+1/2 − ηk−1/2 = ΔAk + ΔBk.
Let (T̂+, v̂+, p̂s+) denote the temperature, horizontal wind vector, and surface pressure at the end of a time step and (T −, v−, ps−) denote the values at the beginning of the time step. In energy-conserving model formulations, the residual
would be zero if there are no diabatic sources and sinks as is the case here. Collins et al. (2004) indicate how to include sources and sinks. Here, Ê+ is given by
The equation for E− is identical to Eq. (8) except that the plus-sign superscript is replaced by a minus-sign superscript and the circumflex or “hat” is removed from E, v, T, and ps.
In general, RES is not zero, because of implicit and explicit diffusion processes in the dynamical cores. Therefore, modifications can be made to the provisional forecast values (T̂+, v̂+, p̂s+). This adjustment yields updated values (T+, v+, ps+) that, if substituted for the provisional values in Eq. (7), yield a zero residual. This is the underlying concept of the energy fixer. The form of the energy fixer first used with the semi-Lagrangian model (FIXER 1) only modified the temperature. The future wind and surface pressure fields are set to v+ = v̂+ and ps+ = Mp̂s+ , where M symbolizes a mass fixer for non-mass-conserving model implementations. The details of the mass fixer are provided in Collins et al. (2004). They are unimportant for the discussion here. In short, the mass fixer ensures that the global integral of the surface pressure and thereby the total mass is conserved. The temperature modification was proportional to the magnitude of the local change in T at that time step and is given by (FIXER 1)
where λ and ϕ denote the longitudinal and latitudinal position in the spherical grid. This adjustment follows the spirit of the water vapor fixer developed by Rasch and Williamson (1991) and Williamson and Rasch (1994) for predecessors of the CAM3. The constant β1 is determined by replacing T̂+ with T+ in Eq. (8) and setting RES = 0 in Eq. (7). Equation (9) is then substituted for T+ in Eq. (8), and Eq. (8) and the equation for E− are plugged into Eq. (7), which is solved for β1.
After it was determined that FIXER 1 of Eq. (9) was the source of the error in the semi-Lagrangian integrations, two other fixers were tried, both of which also only modify the temperature field. FIXER 2 is given by
and FIXER 3 is formulated as
The energy fixer FIXER 2 changes the provisional temperature by a constant, whereas FIXER 3 changes it proportionally. Both constants β2 and β3 are determined as described above for FIXER 1. This yields the equation for β2:
FIXER 2 is adopted in the released version of CAM3.1 and is used operationally for climate applications. Figure 3 confirms that the SLD FIXER 2 and FIXER 3 simulations eliminate the problem with FIXER 1. The figure shows the difference curves for the T170L26 Eulerian simulation versus the T170L26 semi-Lagrangian runs with FIXER 2, FIXER 3, and with no fixer. These three curves all lie on top of each other and fall within the stippled region of the plot. Thus, all three SLD versions calculate the growth of the imposed perturbation to within the uncertainty of the reference solutions. The bottom two panels of Fig. 4 depict the surface pressure at day 10 from the solution with no fixer and with FIXER 2. The two are visually indistinguishable from one another and both are very similar to the T170L26 Eulerian reference solution. The solution with FIXER 3 (not shown) is also visually indistinguishable from the two other semi-Lagrangian solutions. In fact, the l2 difference between the solutions with FIXER 2 and FIXER 3 is 9.4 × 10−2 hPa at day 10, which is well below the differences of any of the three SLD simulations with the Eulerian solution.
As mentioned above, FIXER 1 was in the spirit of the water vapor fixer. The argument there was that if the change in water vapor was small then it was probably not responsible for a lack of conservation and therefore the fixer should be small. With total energy, most of the loss is associated with damping of the wind field and the reduction of kinetic energy. The damping is either explicit or is associated with the numerical approximations. With semi-Lagrangian approximations, it is due in large part to the interpolations fundamental to the method. In hindsight, there is no argument to justify making the fixer proportional to the change in temperature. Although the fixers are very small, FIXER 1 is systematic and apparently accumulates to a significant error. We note that this simple baroclinic instability test case was essential to identify this error. The problem was never noticed with the Held and Suarez (1994) test or with the complete model in which nonlinear feedbacks between dynamics and forcing can mask subtle problems.
The baroclinic instability test case developed by Jablonowski (2004) exposed two flaws in the semi-Lagrangian algorithms originally implemented as an optional dynamical core in NCAR’s CAM3.1 modeling framework. The first flaw resulted in a damping of the speed of a zonally uniform zonal wind undergoing advection by a zonally uniform zonal flow field. It was caused by projecting the transported wind expressed in unit vectors at the arrival point to the surface of the sphere following Bates et al. (1990, 1993), thus systematically shortening the vector at every time step. It is eliminated by rotating the vector to be parallel to the surface following Ritchie (1988, 1991), maintaining its length.
The second flaw was the formulation of an a posteriori energy fixer that, although small, systematically affected the temperature field and lead to an incorrect evolution of the baroclinically growing wave. That fixer changed the provisional temperature forecast at the end of the time step proportionally to the magnitude of the temperature change at that time step. Two other fixers were introduced that do not exhibit the flaw. The dynamical core without an energy fixer also does not exhibit the flaw. The two fixers change the provisional temperature by a constant, or proportionally to the temperature. The former has been adopted by CAM3.1 for all dynamical cores (Collins et al. 2004) and was used in the experiments reported in Jablonowski and Williamson (2006a).
This work was partially supported by the Office of Science (BER), U.S. Department of Energy, through Cooperative Agreement DE-FC02–97ER62402.
* The National Center for Atmospheric Research is sponsored by the National Science Foundation.
Corresponding author address: David L. Williamson, National Center for Atmospheric Research, P.O. Box 3000, Boulder, CO 80307-3000. Email: firstname.lastname@example.org