Enforcing conservation of axial angular momentum in the atmospheric general circulation model CAM6

Numerical general circulation models of the atmosphere are generally required to conserve mass and energy for their application to climate studies. Here we draw attention to another conserved global integral, viz. the component of angular momentum (AM) along the Earth's axis of rotation, which tends to receive less consideration. We demonstrate the importance of global AM conservation in climate simulations with the example of the Community Atmosphere Model (CAM) with the finite-volume (FV) dynamical core, which produces a noticeable numerical sink of AM. We use a combination of mathematical analysis and numerical diagnostics to pinpoint the main source of AM non-conservation in CAM–FV. We then present a method to enforce global conservation of AM, and we discuss the results in a hierarchy of numerical simulations of the atmosphere of increasing complexity. In line with theoretical expectations, we show that even a crude, non-local enforcement of AM conservation in the simulations consistently results in the mitigation of certain persistent model biases.

1 Introduction

The atmosphere exchanges angular momentum (AM) with the material bodies at the surface, which are, to a good approximation, in a state of motion consisting of uniform rotation about the planetary axis connecting the poles. Per unit of mass, surface AM increases in quadratic proportion to its distance from the planetary axis of rotation, from zero at the poles to a maximum at the Equator. AM is a constant of motion of the dynamical (e.g. Newton's) equations so that as air travels meridionally, it carries a specific AM that increasingly differs from that of the Earth's surface. A variety of mechanisms redistribute atmospheric AM and eventually lead to an exchange of AM between the atmosphere and the surface, mainly as a result of low-level wind shear (“surface stress”) and of small-scale wave motions over steep surface topography (“form drag”).

The importance for the atmospheric circulation of the conservation of AM in the free troposphere and of AM exchange of air with the surface was recognised long ago. Already in 1735, George Hadley, Esq, F.R.S., noted that “without the Assistance of the diurnal Motion (i.e. rotation) of the Earth, Navigation (…) would be very tedious” (Hadley, 1735) due to the absence of the trade winds. This insight still lies at the core of modern conceptual models for the atmospheric circulation (Schneider, 1977; Held and Hou, 1980; Lindzen and Hou, 1988; Pauluis, 2004; Walker and Schneider, 2006). In the upper branch of the Hadley circulation (HC), the advection of planetary angular momentum determines a sharp acceleration of the zonal wind in the mid-latitudes linked to a front-like drop in air temperatures, marking the location of the subtropical jets (STJs). Partly by baroclinic instability, the mid-latitude circulation redistributes atmospheric AM vertically and produces intense surface westerlies, whereby the air loses AM to the surface. The equatorward return flow in the surface branch of the HC in turn results in easterly “trade” winds, whereby surface stresses replenish atmospheric AM until air is lifted in cumulus convection within the inter-tropical convergence zone (ITCZ).

This circulation is the object of numerical simulations with general circulation models (GCMs) used in meteorological forecasting and in climate modelling. They describe the atmosphere as a thin, density-stratified, rotating gaseous spherical shell. These properties allow the introduction of a convenient set of approximations in the equations of motion, which result in a system known as the hydrostatic primitive equations (HPEs). The reader is referred to White et al. (2005) for a detailed analysis and discussion. Given suitable boundary conditions, the HPEs guarantee the global conservation of three fundamental physical quantities: mass, energy, and AM along the Earth's rotation axis. Analytic expressions of these laws can be found e.g. in Laprise and Girard (1990). The three conservation laws determine the fundamental character of the large-scale circulation of the atmosphere, and virtually every climate application of GCMs is sensitive to their enforcement when the continuum equations are discretised in space and time. For example, the effects of changes in radiative forcing of ∼2 W m⁻² (e.g. IPCC, 2013, chap. 8, p. 697) can only be simulated if the model's energy conservation is significantly better than 1 %. Estimates based on ECMWF reanalysis data suggest that the conservation of AM of a similar precision is desirable for an accurate representation of the annual cycle and of interannual variations of the atmospheric circulation in model simulations (e.g. Egger and Hoinka, 2005).

CAM, the Community Atmosphere Model developed and maintained at the National Center for Atmospheric Research (NCAR) in Boulder, Colorado, is one of the atmospheric general circulations models (AGCMs) in most widespread use today. It also constitutes the core atmospheric component of NorESM, the Norwegian Earth System Model. Although it offers a choice of dynamical cores, the finite-volume (FV) dynamical core (Lin, 2004) has been, and in many instances still is, the default option. The FV dynamical core is exactly mass and vorticity conserving, and it has been employed in all model integrations submitted by NCAR and by the Norwegian Climate Centre (NCC) for the 5th phase of the Coupled Model Intercomparison Project (CMIP) contributing to the Assessment Report (AR) of the Intergovernmental Panel for Climate Change (IPCC, 2013); it is also expected to be used for phase 6 of CMIP by both institutions. Due to its high numerical efficiency, FV also continues to be the code of choice for all uses for which the overall availability of supercomputing resources is a limiting factor. This includes long historical or palaeoclimate simulations, studies with coupled chemistry and/or carbon cycle, seasonal-to-decadal coupled forecasts, academic research, and all model development efforts currently underway with NorESM.

In this paper, we employ CAM with the FV dynamical core at two standard CESM resolutions only, a coarser one of 1.9^∘ × 2.5^∘ in latitude and longitude, respectively (f19 for short), and a finer one of 0.9^∘ × 1.25^∘ (f09). In agreement with previous results (Lauritzen et al., 2014; Lebonnois et al., 2012), we find that all existing simulations with CAM FV, from CMIP5 to present development versions of CAM6, have a numerical sink of global AM with a magnitude of about 30 % of physical sources at f19 resolution and about 15 % at f09 resolution.

Figure 1 shows the spurious AM source in aquaplanet (AP; Neale and Hoskins, 2000; Blackburn et al., 2013) and Held–Suarez (HS; Held and Suarez, 1994) simulations with CAM FV and an otherwise identical simulation but with the global spectral dynamical core with T42 truncation. Although many other models also do not conserve AM, CAM FV is peculiar in producing a sink nearly everywhere, resulting in a particularly large global non-conservation.

Figure 1 Numerical torque in idealised CAM simulations. The vertically and zonally integrated apparent numerical torque is shown as a function of latitude for CAM simulations in aquaplanet (AP; a, b and c) and Held–Suarez (HS; d and e) configurations. The numerical torque here is obtained as a time-average residual of the tendency of angular momentum in each cylindrical shell of constant latitude of the model's domain, after subtracting the contributions from meridional convergence and from the surface stress torque. The details of the calculation are in Appendix A. Two simulations with the FV dynamical core are shown for each configuration, one at f19 resolution (i.e. on a regular latitude–longitude grid with spacing of 1.9^∘ × 2.5^∘; a and d) and one at f09 (i.e. with twice that resolution; b and e). For comparison, a CAM simulation in AP configuration with the global spectral dynamical core at quadratic triangular truncation T42 (roughly comparable to FV at f19 resolution) is also shown in (c). The dashed red line in each panel indicates the physical torque from surface stresses scaled by a factor 0.1. Positive values indicate an eastward torque acting on the atmosphere, and negative values indicate a westward torque acting on the atmosphere.

First principles (e.g. Held and Hou, 1980; Einstein, 1926) suggest that the dissipation of AM, equivalent to a body force acting on the fluid as a sink of zonal momentum, forces a secondary circulation with the same sign as the Hadley circulation. As a result, the simulated Hadley circulation may become too vigorous. Reduced meridional advection of zonal momentum may lead to mid-latitude westerlies that are too weak or displaced poleward. The zonal momentum lost to the non-physical sink must be balanced by a matching additional eastward torque, for example in an expanded or excessively intense area of tropical easterly surface winds. Model simulations with CAM FV consistently tend to reflect such phenomenology: for example, Feldl and Bordoni (2016) and Lipat et al. (2017) show that among CMIP5 models, those based on the FV dynamical core (GFDL-x, CCSM4, and NorESM-x) simulate both relatively large overturning mass flux in the HC and a high latitude of its edge.

It is useful to illustrate these effects of AM non-conservation by means of idealised AGCM experiments that do not include complicating factors such as orographic form drag or parameterised bulk stresses associated with gravity waves. Figure 2 shows the surface torques resulting from four solutions for the mean circulation with CAM in AP mode. One of these is obtained directly from integrations of CAM using the FV dynamical core at f19 resolution (black line). An otherwise identical integration with the global spectral-transform dynamical core at T42 spectral truncation (green line) is chosen for comparison as a bona fide example of an AM-conserving simulation (see Fig. 1).

Figure 2 Impact of AM sink in CAM–FV integrations. Meridional distribution of the surface stress torque (analogous to the dashed red lines in Fig. 1) in CAM simulations in AP configuration. Two integrations with the FV dynamical core (black and blue lines) and two simulations with the global spectral dynamical core (green and red lines) are shown. One of each pair of integrations is a control case (black and green lines); the other (blue and red lines) is an experiment in which an additional solid-body angular acceleration is applied to the entire atmosphere at each time step of the integration. The acceleration is diagnosed as the time mean of the ratio between the global total numerical torque in the FV control integration and the moment of inertia of the atmosphere. That acceleration is then applied with a negative sign in the FV experiment (blue curve), with the effect of compensating for the numerical torque and achieving approximate global AM conservation in that integration. For the experiment with the spectral dynamical core (red curve), the acceleration is applied with unchanged sign, causing a sink of AM approximately equal to that of the control FV integration. The numerical sink of the control spectral integration is nearly vanishing.

The other two integrations, represented by the blue and red lines, are perturbed in an identical but opposite manner. First, the global total numerical torque due to the FV dynamical core was diagnosed at every time step of the reference FV simulation and averaged in time afterwards. This was converted into a solid-body axial rotation tendency that was applied continuously everywhere as a constant sink of AM in a new integration with the spectral dynamical core, resulting in the simulation represented by the red curve. Vice versa, the opposite additional solid-body rotation tendency was applied to a new FV integration, thus compensating for its internal numerical sink. This integration produced the physical torque represented by the blue curve. Comparing the different curves, it may be seen that equatorward of about 23^∘ of latitude the simulated physical torque depends primarily on the global budget of atmospheric AM. In particular, notwithstanding the complications of interactive moist physics and the different spatial and temporal discretisations used in the two integrations, the stronger trade winds (in terms of surface stress) in the FV simulation compared with the T42 simulation can be explained entirely with the non-physical, numerical torque of the FV dynamical core. The result is insensitive to how that torque is in fact applied. Even at subtropical and middle latitudes, half of the difference between the two simulations, in terms of surface stresses, can be explained in this way. Similar results are found for the zonal-mean meridional circulation and for the surface pressure in the HC (Fig. S1 in the Supplement), confirming the strength and robustness of the Einstein (1926) “tea leaves” mechanism.

These results motivate us to address the issue of AM conservation in the CAM's FV dynamical core. One may speculate that systematic biases in surface stresses due to the numerical sink of AM must also impact coupled ocean–atmosphere climate simulations, with excessive Ekman and Sverdrup forcing of the subtropical gyres. The northward displacement of the mid-latitude westerlies may also result in excessive mechanical and thermal forcing of the subpolar gyres with possible implications for the Atlantic meridional overturning circulation.

In this paper, we propose ways to address the numerical dissipation of AM in CAM–FV simulations. Section 2 describes our main hypotheses as to the root cause of the error and our approaches towards rectification. Section 3 presents the result of our corrections in a set of idealised simulations. The impact on realistic simulations of the atmospheric circulation is discussed in Sect. 4. Conclusions are finally offered in Sect. 5.

2 Analysis of potential causes and approaches to correction

The FV dynamical core (Lin, 2004) solves the HPE by updating first the advective (C grid) and then the prognostic (D grid) winds in two steps. The first step represents pure advection, i.e. the increments associated with transport, including geometric and Coriolis terms. In this step, the scheme conserves absolute vorticity exactly for the D-grid winds (Lin and Rood, 1997; hereafter LR97). The second step calculates the wind increments associated with hydrostatic pressure forces. These are computed in a special way (Lin, 1997) that differs from most Arakawa and Lamb (1981) type schemes. Violations of AM conservation may occur in either sub-step.

2.1 Pressure-gradient force

We first analysed the Lin (1997) treatment of the pressure-gradient terms for conservation. A general discussion is given by Simmons and Burridge (1981), who introduce a set of hybrid-level dimensionless variables, a_k, defined as $a_k:=({\mathit{\varphi }}_k-{\mathit{\varphi }}_{k+\mathrm{1}/\mathrm{2}})/\mathrm{2}(\mathit{\alpha }p{)}_k$ (in Simmons and Burridge these variables are denoted by α_k; we change the notation here to avoid confusion), where ϕ is the geopotential, p the pressure, $\mathit{\alpha }:=-{\partial }_{\mathit{\eta }}\mathit{\varphi }/{\partial }_{\mathit{\eta }}p$ the specific volume, and η is the generalised or hybrid vertical coordinate. Here and in the following, the index k refers to the vertical level, or to half-levels as appropriate, and subscripts to the partial derivative symbol indicate differentiation with respect to the variable in subscript, ${\partial }_X\equiv \partial /\partial X$. The variables a_k need not be constants. Simmons and Burridge (1981) derive the discrete form that pressure and geopotential terms must take in general vertical coordinates in order to ensure conservation of axial angular momentum. Their Eq. (3.8) can be generalised to

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\left(\mathit{\alpha }{\partial }_{\mathit{\lambda }}p+{\partial }_{\mathit{\lambda }}\mathit{\varphi }\right)}_k& =-{\left({\displaystyle \frac{\mathrm{\Delta }\mathit{\varphi }}{\mathrm{\Delta }p}}\right)}_k{\partial }_{\mathit{\lambda }}{p}_{k-\mathrm{1}/\mathrm{2}}\\ & +{\partial }_{\mathit{\lambda }}{\mathit{\varphi }}_{k+\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{p}_k}}{\partial }_{\mathit{\lambda }}\left[a_k(\mathit{\alpha }p{)}_k\mathrm{\Delta }{p}_k\right],\end{array}\end{array}$$

where the symbol Δ is employed to represent a difference between vertical levels, $\mathrm{\Delta }{p}_k:=p_{k+\mathrm{1}/\mathrm{2}}-p_{k-\mathrm{1}/\mathrm{2}}$ (and similarly for ϕ), and λ is the longitude.

Performing Lin's (1997) path integration around the finite-volume element on this expression yields the following form for the body force:

$$\begin{array}{}\text{(2)}& \oint \mathit{\varphi }\mathrm{d}p={\mathit{\delta }}_{\mathit{\lambda }}\left\{\left[{\mathit{\varphi }}_{k+\mathrm{1}/\mathrm{2}}+a_k(\mathit{\alpha }p{)}_k\right]\mathrm{\Delta }{p}_k\right\}-\mathrm{\Delta }{\left(\stackrel{\mathrm{‾}}{\mathit{\varphi }}{\mathit{\delta }}_{\mathit{\lambda }}p\right)}_k,\end{array}$$

where δ_λ is the finite-difference operator in the zonal direction, and $\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_{k\pm \mathrm{1}/\mathrm{2}}}$ is an average over λ. An expression identical in form to Lin's (1997) Eq. (11) is then recovered if the choices

$$\begin{array}{}\text{(3)}& a_k={\displaystyle \frac{\mathrm{\Delta }{\mathit{\varphi }}_k}{\mathrm{2}(\mathit{\alpha }p{)}_k}},\stackrel{\mathrm{‾}}{\mathit{\varphi }}={\displaystyle \frac{{\mathit{\varphi }}_{i+\mathrm{1}/\mathrm{2}}+{\mathit{\varphi }}_{i-\mathrm{1}/\mathrm{2}}}{\mathrm{2}}},\end{array}$$

are made, where i is the index corresponding to the longitude λ.

In other words, Lin's (1997) expression for the pressure-gradient term is consistent with the Simmons and Burridge (1981) prescription for AM conservation, provided that the physical pressure variable p is used in the integration in place of the general pressure function indicated by the symbol π in Lin (1997). This can be directly verified algebraically by summing all expressions of the form of the numerator in the right-hand side of Eq. (11) in Lin (1997) along all longitudes and levels. Provided ϕ is constant at one model boundary and p at the other, it always returns zero. This is the required result provided that the denominator on the right-hand side of Eq. (11) in Lin (1997) represents the inertial mass associated with the velocity points. They do so if π is the hydrostatic pressure.

Accordingly, we performed tests in which the integration variable in the relevant section of CAM–FV's dynamical core was replaced with true interface pressure. The effect was generally seen to be very small on the dynamical core's momentum conservation properties.

We note, however, that in the CAM implementation there may be an additional problem associated with the use of the D grid. The application of Lin's (1997) method would strictly require a C grid, with zonal velocity points interleaving pressure (scalar) points along the same latitude. Thus, in CAM pressure is interpolated to the grid cell corners before use. While the formal expressions for the pressure forces do not change, thus ensuring the Simmons and Burridge (1981) total torque constraints, the inertial mass associated with each D-grid zonal velocity point is in fact averaged over six scalar points surrounding it, with 1-2-1 weights along the zonal direction. This additional zonal smoothing effectively adds spurious terms to the zonal momentum equation of the form $-u{\partial }_x^{\mathrm{2}}\mathrm{\Delta }p$. This is a potential source of non-conservation. However, it is not expected to be systematic.

2.2 Geometry, polar filtering, and FFSL extension

AM conservation may be affected by the treatment of geometric terms in latitude–longitude coordinates, especially near the poles where such terms become large. Furthermore, convergence of the meridians forces filtering of the solution, and additional approximations need to be made. In particular, LR97 implement a flux-form semi-Lagrangian extension of Colella and Woodward's (1984) piece-wise parabolic method (PPM), which is used near the poles where Courant–Friedrichs–Lewy numbers (Courant et al., 1928) become large during the time integration. We performed several sensitivity tests on each of these aspects without being able to notice significant impacts on AM conservation.

Particularly compelling is the comparison with the performance of a prototype implementation in CAM of the FV scheme on a cubed-sphere grid (FV3), which lacks any poles and does not require or use any of these special formulations (and is, in particular, run in pure Eulerian mode, i.e. without the flux-form semi-Lagrangian extension described in Lin and Rood, 1996). We ran an AP simulation on the C48 grid, viz. six pseudo-cubic faces with 48×48 grid cells each, for a total number of grid points identical to the standard 2^∘ FV configuration but a 25 % higher resolution at the Equator. The AM sink (Fig. S2 in the Supplement) is nevertheless comparable, i.e. about 25 % smaller, consistent with the scaling with the resolution of simulations with standard FV. We conclude that FV and FV3 suffer from the same problem independent of geometry or the FFSL extension of LR97.

In order to minimise the impact of other minor (and partly intentional) numerical sources and sinks of AM, in all idealised numerical tests presented in this paper we applied the following modifications: (1) the order of the advection scheme is kept the same (fourth) for all model layers, instead of reducing it to first in the top layer and to second up to the eighth layer; (2) an additional conservation check is applied in the vertical remapping of zonal wind , and column momentum is conserved in the moist-mass adjustment at the end of physics; (3) the surface stress residual resulting from closure of the diffusion operator (in physics) is applied in full rather than partially.

2.3 Discretisation of the kinetic energy term

The evidence from our theoretical and diagnostic analysis points at the advective, shallow-water part of the implementation of LR97 in CAM–FV as the root of the AM conservation error. Its “vector-invariant” formulation (Arakawa and Lamb, 1981) allows for different forms of the divergence to be used in the momentum and in the mass and tracer equations, resulting in inconsistent values for the divergence of the flux of planetary AM (associated with mass divergence) and of the flux of relative AM (associated with momentum divergence). In the momentum equations, the divergence is contained in a kinetic energy (KE) gradient term, which due to the presence of a numerical symmetric instability (Hollingsworth et al., 1983) is expressed as the local gradient of a Lagrangian-average KE. Its form violates the finite-volume approximations used for other quantities, e.g. vorticity. This feature is intrinsic to the LR97 numerical discretisation scheme and cannot be eliminated.

To address the resulting violation of AM conservation, we first note that even in AM-conserving schemes, conservation can only be guaranteed in the zonal average (Simmons and Burridge, 1981). We therefore do not attempt a local correction to the scheme, which is liable to numerical instabilities (Hollingsworth et al., 1983), and instead formulate a zonal-mean correction as follows. We enforce the AM conservation law,

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\int \mathrm{d}\mathit{\lambda }{\partial }_t\left(\mathrm{\Delta }pua{\cos }^{\mathrm{2}}\mathit{\phi }\right)& =-\int \mathrm{d}\mathit{\lambda }{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)\\ & +\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }pfva{\cos }^{\mathrm{2}}\mathit{\phi },\end{array}\end{array}$$

by adding a zonal-mean zonal-wind tendency term to the vector-invariant form:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\partial }_{t,c}\stackrel{\mathrm{‾}}{u}& ={\displaystyle \frac{\mathrm{1}}{\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }p}}\times \mathit{\{}\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }p\left({\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}K-\mathit{\zeta }v\right)\\ & -\int \mathrm{d}\mathit{\lambda }{\displaystyle \frac{\mathrm{1}}{a{\cos }^{\mathrm{2}}\mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)\\ & -\int \mathrm{d}\mathit{\lambda }u{\partial }_t\mathrm{\Delta }p\mathit{\}}.\end{array}\end{array}$$

Here, K is the KE plus the contribution from explicit divergence damping used in FV. In the continuum limit the expression on the right-hand side simply reduces to the mass-weighted zonal average of the zonal gradient of $K-(u^{\mathrm{2}}+v^{\mathrm{2}})/\mathrm{2}$.

In discrete form, the last two terms must be approximated. In the C–D grid formulation of the LR97 scheme the second one is especially problematic. Various possibilities were explored, which resulted in various degrees of accuracy and stability. The best compromise is to discretise it as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{a{\cos }^{\mathrm{2}}\mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)=\\ & {\displaystyle \frac{\mathrm{1}}{a{\cos }^{\mathrm{2}}\mathit{\phi }}}\left[\mathrm{\Delta }pv{\partial }_{\mathit{\phi }}\left(u\cos \mathit{\phi }\right)+u{\partial }_{\mathit{\phi }}\left(v\mathrm{\Delta }p\cos \mathit{\phi }\right)\right]\end{array},\end{array}$$

where u is the prognostic D-grid zonal wind and v is the advective (C-grid) meridional wind. The details of the derivation are given in Appendix B. Using the mass conservation equation, this approximation allows us to discretise the two last terms together and write the zonal-wind correction increment in a form consistent with LR97:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\delta }}_c\stackrel{\mathrm{‾}}{u}& ={\displaystyle \frac{\mathrm{1}}{\int \mathrm{d}\mathit{\lambda }\stackrel{\mathrm{‾}}{\mathrm{\Delta }{p}_{t+\mathit{\delta }t}}}}\int \mathrm{d}\mathit{\lambda }\mathit{\left\{}\stackrel{\mathrm{‾}}{\mathrm{\Delta }p}\right[{\displaystyle \frac{\mathit{\delta }t}{a\cos \mathit{\phi }\mathit{\delta }\mathit{\lambda }}}{\mathit{\delta }}_{\mathit{\lambda }}K\\ & -\stackrel{\mathrm{‾}}{\mathcal{Y}\left(v^{*},\mathit{\delta }t;{\mathit{\zeta }}_{\mathit{\lambda }}\right)}]+{\stackrel{\mathrm{‾}}{u}}^t\mathcal{F}\left(u^{*},\mathit{\delta }t;\stackrel{\mathrm{‾}}{\mathrm{\Delta }p}\right)+O\left(\mathit{\delta }{t}^{\mathrm{2}}\right)\mathit{\}}.\end{array}\end{array}$$

Here, ${\mathit{\zeta }}_{\mathit{\lambda }}:=\frac{\mathrm{1}}{a\cos \mathit{\phi }}{\partial }_{\mathit{\lambda }}v$, and the notation of LR97 is used for the discrete transport operators 𝒴 and ℱ for the meridional transport of ζ_λ and the zonal transport of mass, respectively. The first three terms in the integrand of Eq. (7) thus correspond to the first three terms (first and second lines) on the right-hand side of Eq. (B10) in Appendix B. The last symbol on the right-hand side of Eq. (7) represents higher-order terms (also detailed in Eq. B10). We will refer to this modification of the LR97 scheme as the “correction”.

2.4 Diagnostic tools and global conservation

Irrespective of whether the correction, as described above, is applied or not, for diagnostic purposes we calculate the apparent non-physical torque associated with the FV dynamical core advective tendencies only, i.e. excluding the increments associated with pressure gradients. These tendencies are diagnosed separately for each layer at every advective sub-step and integrated horizontally to yield the apparent numerical global total torque during the sub-step. At the same time, the layer effective moment of inertia over the sub-step is also computed.

The opposite of the ratio of these quantities gives an angular acceleration that, applied to the zonal wind in each layer at every advective sub-step, enforces conservation of the AM of that layer under advection. The application of this solid-body rotation increment at each dynamical time step and for each layer independently is what we call the “level” fixer. The details of the computation are given in Appendix C.

Irrespective of whether they are actually applied, the fixer's velocity increments (Eq. C2), are vertically interpolated and accumulated over the entire dynamic time step and written out diagnostically. In addition to the fixer, partial wind and pressure tendencies arising from the dynamical core are separately diagnosed and written to the standard output streams, providing additional diagnostic tools for cross-checking.

A variant of the fixer was tested in CAM simulations. This variant is a “global” fixer, which still acts by applying an increment to the zonal wind at each time step. In this fixer, the apparent torque and the moment of inertia are integrated over all levels within the domain over which strict overall angular momentum conservation is desired. The zonal-wind increments are then applied as a single solid-body rotational acceleration within this domain. Experimentation showed that such acceleration should not be applied in the stratosphere, where conservation errors are small and the impact of unphysical zonal accelerations large. The necessary limitation of the domain for the global fixer, however, introduces a certain degree of arbitrariness in its application. Although sometimes used for diagnostic purposes, we do not discuss this global fixer variant any further.

Lin's (2004) FV scheme conserves mass and absolute vorticity exactly. The AM modifications described above were explicitly designed not to alter the mass flux calculations and intervene only on the rotational component of the flow in the momentum equations. Other choices, involving alterations to the calculation for the divergent flow, would have been possible. However, we judged exact mass conservation more important for climate simulations than exact vorticity conservation. The AM modifications also change the kinetic energy of the flow and thus change the total energy budget of the model. However, the unmodified FV scheme does not conserve energy. CAM–FV therefore employs an energy “fixer” (analogous to out AM fixer), described e.g. in Williamson et al. (2015). The fixer diagnoses the energy non-conservation at each time step. This allowed us to monitor the impact of the AM mods on energy non-conservation in all our experiments. We found no systematic effect, either in sign or magnitude, of the AM modifications on the energy non-conservation of the model.

3 Numerical simulations and results

3.1 Dry baroclinic wave tests

Initial tests were carried out for adiabatic dynamics and flat bottom topography from baroclinically unstable initial conditions, as defined in Jablonowski and Williamson (2006; JW06). Figure 3 shows the result in terms of the conservation of global AM for CAM–FV integrations at f19 resolution (1.9×2.5^∘ of latitude and longitude) and 30 hybrid levels.

Figure 3 AM correction and fixer in adiabatic, frictionless baroclinic wave tests. Three sets of curves are shown for each of four different simulations with CAM FV, indicating the time evolution of global AM (diamond shapes) and its two components of planetary AM (vertical crosses) and relative AM (× crosses) in each simulation. Total AM and each AM component are normalised to the initial total AM of the initial state, and differences with respect to initial values are shown, expressed in percentage. Standard CAM–FV is shown in black, CAM–FV with the AM correction only in blue, CAM–FV with the AM fixer only in yellow, and CAM–FV with both the AM correction and fixer in red. The inset panel on the lower right of the figure shows an enlargement for the initial evolution of total AM. Note that the four simulations are nearly indistinguishable before day 8, i.e. during the linear phase of the baroclinic wave. All simulations are run on the 2^∘ grid.

It may be seen that both the correction and the fixer are effective in reducing the systematic numerical sink of AM in these integrations. In particular, the fixer appears to remove it almost completely; in other words, the integration with the fixer conserves global AM in the time mean. This result is central to this paper, and it proves its two main conclusions. The first is that the systematic non-conservation of global AM in the FV dynamical core indeed resides in the advective wind increments of the shallow-water part of the dynamical core. The second is that, by virtue of its effectiveness and its formulation that is entirely independent of the model configuration or parameterisations (topography, physical momentum sources, etc.), the fixer is a useful and accurate general diagnostic tool that allows us to quantify the numerical torque in any CAM–FV integration. By virtue of this quality, the diagnosed time-averaged fixer tendencies were, for example, used for the perturbations in the experiments shown in Figs. 2 and S2.

The impact of the correction on conservation is generally smaller, and different dynamical regimes may be seen when the size and quality of that impact change. In the baroclinic instability tests in Fig. 3, the correction achieves good results in the linear and non-linear stages of baroclinic growth (up to day 30; see JW06) but is not able to correct the slow drift that sets in after zonalisation of the global flow; then wind speed decreases everywhere as a result of numerical dissipation (there are no external sources or sinks of either momentum or energy in these adiabatic simulations). This is a partly desirable behaviour, as the action of the correction should not change the dissipation properties of the scheme.

Aside from the conservation properties they are designed for, both the correction and the fixer represent a perturbation of the numerical solutions of the FV dynamical core. By arbitrarily modifying the relative vorticity associated with the zonal wind, both destroy one of the fundamental numerical properties of the LR97 formulation, viz. the conservation of absolute vorticity under advection. (In the case of the fixer, the vorticity input has a rigid dependency on latitude, ∼sin φ.) Figure 4a shows their impact on the accuracy of the JW06 baroclinic wave test in terms of the root mean square (RMS) of the differences in surface pressure from a nominal reference solution with original FV dynamical core. The latter is obtained at f19 resolution (0.9^∘ × 1.25^∘), which is sufficiently close to JW06's reference solution (see JW06, Section 5(e), points i and ii) for our purposes. It may be seen that on this measure the solutions with and without the AM corrections are virtually indistinguishable during the stages of both linear and non-linear baroclinic growth. A similar result holds for the phase (not shown).

Figure 4 AM correction and fixer in an adiabatic, frictionless baroclinic wave test. The simulations shown in Fig. 3 are compared with a standard CAM–FV simulation at 1^∘ resolution and against each other. Each panel shows seven curves, four of which nearly overlap and form the topmost set of lines (including the reference simulation with standard FV). These represent the time evolution of the RMS difference of surface pressure (a) and relative vorticity at 230 hPa (b) for each of the 2^∘ integrations and the control 1^∘ integration. Below that set of curves are two nearly overlapping curves, which show the RMS differences of the 2^∘ experiments with AM correction only and the control 2^∘ integration (blue lines), as well as of the experiment with both the AM correction and fixer and the control integration (red lines). Finally, the single yellow lines at the bottom in each panel show the RMS differences of the 2^∘ integration with the AM fixer only with the 2^∘ control integration.

It may be noted that the largest impact on the RMS of surface pressure arises from the correction. Within the first 30 d this impact is formally always well below significance (as defined in JW06; see their Fig. 10), but it increases in time and eventually becomes appreciable as a full global meridional circulation is established. Similar results hold for the vorticity field, as seen in Fig. 4b.

Other aspects of the solution besides RMS differences also show limited sensitivity to the application of the correction and the fixer. Figure 5 shows the evolution of the minimum pressure in the developing baroclinic wave. By this measure, the solutions only start to diverge with the filling of the primary cyclone and the deepening of the secondary wave after day 17. The solution with the fixer deepens the secondary cyclone more quickly so that the minimum pressure is seen to jump from the first to the second wave minimum between days 18 and 19; this occurs 1 d later with the unmodified dynamical core. A third transition after day 25 has higher central pressure in the solutions with the fixer; by this time, however, rapid cyclogenesis is occurring in the jet stream of the Southern Hemisphere, attaining a similar minimum pressure, which is slightly deeper in the solutions with the fixer. In any case the pressure differences of the minima remain of the order of a few hectopascals (hPa), and there is no systematic difference in their position.

Figure 5 AM correction and fixer in an adiabatic, frictionless baroclinic wave test. Evolution of minimum pressure (a) and its position (c, d) in the baroclinic wave evolution from the integrations shown in Fig. 3. The colour coding of the lines is the same as in Fig. 3. Panel (b) shows the differences in minimum pressure between the AM experiments and FV control, with the same colour coding as in the lower curves in Fig. 4.

3.2 Other idealised tests

Even if the impacts of the modifications of the FV dynamical core are relatively small on local circulations over subseasonal timescales, as shown above, the rationale for introducing them is the hope of achieving a better simulation of the state of the atmosphere in integrations under specified forcings. As explained in the Introduction, one particular expectation is that the subtropical easterlies should weaken without affecting the circulation elsewhere too heavily. In particular, the role of the correction, which alone does not ensure AM conservation, must be clarified and its eventual use justified. Here we document the results of two sets of idealised simulations that still have a simplified, equipotential lower boundary but include non-vanishing physical torques and heating tendencies.

The first set of such simulations adhere to the benchmark test of Held and Suarez (1994; HS henceforth), whereby the forcing has the form of a relaxation towards a specified three-dimensional atmospheric temperature field. Likewise, surface friction is represented by a damping of the winds within a set of levels near the bottom boundary. Apart from the small numerical diffusion, these stresses are communicated to the rest of the atmosphere by means of momentum advection in the mean circulation and of pressure fluctuation in resolved transient motions (including travelling waves). The second set of simulations follows the aquaplanet (AP) test first proposed by Neale and Hoskins (2000), wherein only a persistent field of bottom-boundary temperatures is prescribed (the QOBS profile of Neale and Hoskins, 2000), and the full set of moist atmospheric physical parameterisations of CAM6 are used to force the circulation (except for those specific to orographic processes). The bottom boundary is a notional static ocean with unlimited heat and water capacity. Surface stresses are computed by the coupler and passed to the moist atmospheric boundary-layer parameterisation, which then distributes those stresses vertically. Momentum is also transported in moist convection, where active, and further adjustments are made when the moist mass of the atmospheric column changes due to precipitation and surface evaporation processes. To simplify the analysis, the gravity-wave parameterisation of CAM6 was turned off in our AP tests. In both sets of tests, FV's advection scheme is used at PPM's standard fourth order at all levels; i.e. the numerical diffusion obtained in standard CAM–FV integrations by employing low-order calculations near the model top is avoided. For initial conditions, HS simulations are cold-started with uniform surface pressure and geopotential, as well as vanishing wind fields except for a westerly perturbation identical to that used in the dry baroclinic wave tests (necessary in order to break zonal symmetry and to allow for a non-vanishing correction). The AP simulations all take the same instantaneous atmospheric state from a previous spun-up run, even though this requires more adjustment for the corrected (fixed) simulations than for the control.

Figure 6 AM correction and fixer in Held–Suarez (HS) and aquaplanet (AP) integrations. Panel (a) shows the time evolution of total AM for each of the integrations, similar to Fig. 3 (diamond shapes) but normalised separately for each integration to the time-integrated physical (i.e. surface stress) torque at day 360. AP integrations are shown in solid, HS integrations in stippled lines. The colour coding is as in Fig. 3. Panel (b) shows the time-mean numerical torque, averaged over days 120–360, arising at each model level from advective increments, as diagnosed by the fixer and expressed as equatorial acceleration in a solid-body rotation required to compensate for the numerical sink. Line types and colours correspond to those shown in (a). The lists at the bottom of (b) indicate the time-mean equatorial accelerations of a global solid-body rotation, i.e. the increments shown by the lines but integrated vertically level by level, weighted with the appropriate moments of inertia.

Figure 7 Impact of the AM correction and fixer in Held–Suarez simulations. Time-mean latitude–pressure profiles of wind differences between HS simulations shown in the stippled lines in Fig. 6. Panel (a) shows the zonal-mean zonal-wind time-average (days 120–360) difference field of the integration with the AM fixer only and the control integration. Panel (b) shows the same field, but for the difference between the integration with the AM correction and control. Panel (c) shows the difference between the integration with both the AM correction and AM fixer and control; (d) the difference between the integration with both the AM correction and AM fixer and the integration with the AM fixer only. The contour interval is 0.6 m s⁻¹, with blue hues indicating negative values and red hues positive values. Values in the interval [−0.3, +0.3] m s⁻¹ are left in grey. The fields displayed have been symmetrised about the Equator, since departures from symmetry are very small in the time mean for these hemispherically symmetric simulations. Accordingly, only one hemisphere and the equatorial region are shown in each panel.

Figure 6a indicates that the global AM conservation properties of the simulations in these tests are broadly in line with the expectations from the previous discussion. Standard FV tests (black lines) show a steady loss of AM in the atmospheric circulation, with a magnitude of the order of 10 %–20 % of the physical flux of AM through the atmosphere. (We count eastward stress as positive, by which the atmosphere gains westerly momentum in the tropical surface easterlies and loses westerly momentum in the subtropical surface westerlies.) Use of the correction leads to an order-of-magnitude reduction of the numerical sink of AM in HS integrations, but it is of limited effectiveness in full-physics AP integrations (blue lines). Integrations with the fixer, with or without the correction (orange and red lines, respectively), maintain atmospheric AM in the time mean. In HS simulations, there appears to be a very small residual drift of AM notwithstanding the fixer. This is due to a small inconsistency in the application of the stress terms, which are calculated and diagnosed in the “physics” part of the model time stepping but applied later as velocity tendencies in the physics–dynamics interface on updated layer masses. This is an intrinsic feature of the time stepping of CAM–FV that we have not modified. More notably, AP simulations differ from HS simulations in that they show obvious fluctuations of total AM around the time mean or around the long-term drift when there is one. Such fluctuations are similar in all AP integrations, with a magnitude of a few percent of the physical sources, and depend on non-conservation in CAM's physics parameterisations. Fortunately, they are not systematic and do not produce a noticeable long-term drift.

Figure 8 Impact of the AM correction and fixer in aquaplanet simulations. Same as Fig. 7, but for the AP simulations shown in the solid lines in Fig. 6.

The effectiveness of the fixer in removing most of the AM drift confirms that the systematic sink of AM in CAM–FV integrations arises predominantly from the shallow-water advection calculations. The accuracy of the correction, by contrast, depends on the features of the circulation, with good accuracy for numerically well-resolved features, as in the HS tests, but a poorer one when grid-scale forcing associated with the water cycle occurs. Figure 6b gives more details on the effect of the correction. Here, the time-average AM sink due to the dynamical core is diagnosed using the fixer increments for the zonal velocity at the Equator at each model level. This diagnostic is produced irrespective of whether such increments are applied during the integration. Apart from the smaller increments in HS integrations than in AP integrations, which partly depend on the slower circulation (“surface” stresses are 1 order of magnitude larger in the HS set-up than in the AP set-up), the advective AM sink has a distinctive shape in pressure-level space, with a maximum in the upper troposphere and small values in the atmospheric boundary layer. This shape partly reflects the underlying global mean zonal-wind field, but the maximum sink lies below the maximum wind (at around 250 hPa rather than around 150 hPa). The profile of the impact of the correction, i.e. the reduction in fixer increments when the correction is applied, has again a similar shape but with an even lower position of the maximum, which better corresponds to the maximum in the vertical profile of the level-integral zonal momentum of the underlying flow. Combined with the offline diagnostic information for the apparent AM sink from Fig. 1, it can be deduced that the main loci of the time-mean AM sink in these simulations are found near the subtropical jet streams, where large zonal asymmetries occur in both the mass fields and the wind fields.

The effect on the mean circulation of applying the correction and/or the fixer is shown in Figs. 7 and 8 for HS and AP simulations, respectively. The zonal-mean zonal winds are shown, which is the quantity that both the correction and the fixer directly modify. Nevertheless, it should be remembered that the net effect is indirect, since the zonal winds remain in the time average close to geostrophic balance with the (equivalent) temperature field. In HS simulations, the local temperature differences between simulations are simply proportional to the difference in temperature advection by the meridional and vertical circulation, which is modified primarily through a “tea leaves” mechanism. As already seen in the Introduction, the leading-order effect of the fixer is a weakening of this circulation and thus of the associated advective temperature tendencies. These tend to cool the lower troposphere in the subtropical easterlies, cool the upper troposphere near the Equator, and warm the troposphere poleward of the jet streams. The effect of the fixer on the zonal-mean zonal wind shown in Fig. 7a is generally consistent with this expectation, with an equatorward retreat of the surface easterlies and weaker westerlies in the higher latitudes. There is, however, an additional large westerly difference near the equatorial tropopause; this is a direct consequence of the westerly forcing of the fixer, which is greatest at the Equator. This is clearly an undesirable effect of the fixer on the simulations. A more selective effect on the circulation is produced by the correction (Fig. 7b). As seen above, its main action is in where the greatest sink of AM is located, i.e. on the flanks of the subtropical jet stream. By correcting part of the AM non-conservation, it also acts to limit the action of the fixer (Fig. 7d). As a result, the combination of the correction and fixer together, as well as ensuring good global AM conservation, is less severe in terms of its upper-level equatorial westerly effect (Fig. 7d). This suggests that the fixer is best employed in combination with the correction.

In AP simulations, a slowdown of the meridional circulation is still expected and found, but the interaction between dynamical forcing by the fixer or the correction and the physics tendencies is much more complex and difficult to predict. The fixer now produces large westerly differences near the Equator at all levels and a marked weakening of the subtropical jet stream (Fig. 8a). The equatorial winds above 300 hPa become westerly. The correction is less effective overall than in HS simulations, and its impacts are mostly confined to levels close to the model lid or to the high latitudes (Fig. 8b). Nonetheless, its use is still beneficial in terms of limiting the action of the fixer, at least in the troposphere (Fig. 8d). The result of the combined correction and fixer can be seen in Fig. 8c. In terms of tropospheric impacts, it appears acceptable; equatorial winds remain easterly below 200 hPa and weak above. The weakening of the equatorial and tropical easterlies compared with the control simulation implies greater similarity with simulations with AM-conserving spectral models. Large changes, however, can be seen near the model lid, especially in the four model layers with pressures less than 25 hPa. This is a consequence of momentum accumulation within these layers. In CAM's default configuration, the order of FV's PPM advection scheme is reduced here, which results in large numerical dissipation. Effectively, these levels are used as sponge layers and are thus not part of the valid computational domain of the model. In full-model configurations it is therefore advised to keep the reduced order of advection and turn off both the correction and the fixer in these layers. The large mean-state changes seen near the top in Fig. 8d then vanish. Considering the troposphere only, the conclusion obtained from HS simulations can be seen to also hold for full-physics AP model simulations in that the combined application of the fixer and the correction results in smaller overall mean-state changes in the solution compared to default FV without modifications, while ensuring good conservation of AM.

Figure 9 Impact of the AM correction and fixer in F2000 simulations. Panels (a), (b), and (c) show maps of surface wind-stress vector differences (arrows) and wind-stress magnitude differences (colours) between F2000 simulations with CAM–FV at 1.9^∘ × 2.5^∘ resolution (f19) and a climatology obtained from satellite scatterometer observations (ERS; Quilfen et al., 1999). Panel (a) shows the annual-mean climatological bias in the f19 control integration, (b) for an f19 simulation with the AM fixer only, and (c) for an f19 simulation with both the AM correction and AM fixer. Panels (d) and (e) show the same fields, but for the differences between the simulation with the fixer only and the control and between the simulation with both the fixer and correction and that with the fixer only. The colour scale for all plots is on the right of (d) and (e). These plots were produced with the AMWG diagnostics package developed by the Atmospheric Model Working Group of the University Corporation for Atmospheric Research and the National Center for Atmospheric Research.

4 Simulations of the observed climatology

The relevance of the AM modifications to the FV dynamical core for CAM simulations in a realistic configuration is investigated here using F2000 cases, which are AMIP-type simulations (Gates, 1992) wherein sea surface temperatures (SSTs) and all compositional forcings are prescribed as a repeating annual cycle obtained from an observed climatology of the decade spanning the turn of the century. We test at two grid resolutions, one of 1.9^∘ × 2.5^∘ (f19) as in all integrations already discussed above and one of 0.9^∘ × 1.25^∘ (f09) to test the impacts of AM modifications in a case that is scientifically supported by NCAR at this time. The CESM model version used (here as above) is release 2.1.1¹.

Figure 10 Impact of the AM correction and fixer in F2000 simulations. Latitude–pressure maps of zonal-mean zonal-wind climatologies for boreal winter (DJF). Panels (a), (b), and (c) show total fields for the CAM–FV f19 control simulation: (a) for the f19 simulation with both the AM fixer and AM correction (b), as well as for the ERA40 reanalysis (Uppala et al., 2005). The colour scale is at the bottom of (a). Panels (d) and (e) show the differences of each of the two f19 integrations and ERA40, and (f) shows the differences between the two f19 simulations. The colour scale is on the right of (f). Panels (g), (h), and (i) are analogous to (d), (e), and (f), respectively, but for CAM–FV simulations at 0.9^∘ × 1.25^∘ resolution. These plots were produced with the AMWG diagnostics package developed by the Atmospheric Model Working Group of the University Corporation for Atmospheric Research and the National Center for Atmospheric Research.

Figure 11 Impact of the AM correction and fixer in F2000 simulations. Panels (a) and (b) show Taylor (2001) diagrams for the validation of the CAM–FV F2000 simulations at f19 (a) and at f09 (b) resolution against observations for a standard set of diagnostic fields, which are listed in the panels. Black symbols represent RMS differences to observations for the control simulations without modifications, and red symbols are for the simulations using both the AM correction and the AM fixer. For the overall RMSE and bias scores, those from the control simulations are used as normalisation. Panel (c) summarises the correlation values between simulated and observed diagnostic fields, as listed in the central table. Green fields mark all instances in which one of the AM-modified simulations represents an improvement over the respective control simulation. These plots were produced with the AMWG diagnostics package developed by the Atmospheric Model Working Group of the University Corporation for Atmospheric Research and the National Center for Atmospheric Research.

Figure 9 illustrates the effects of the fixer and the correction on f19 simulations. The control simulation shows a characteristic easterly surface wind-stress bias throughout the tropics (Fig. 9a). In addition, there are excessive westerlies at southern high latitudes. The effect of the fixer is to reduce the tropical biases (Fig. 9b), with an evident westerly effect on the simulations nearly symmetrically about the Equator (Fig. 9d). By that same token, however, the high-latitude westerly errors are enhanced (Fig. 9b). The application of the correction in addition to the fixer not only brings further improvements in the tropics, but also corrects the westerly effect of the fixer in high latitudes (Fig. 9e). The result is a significant improvement in the simulation of the surface wind-stress field over the entire ocean domain.

In general, we obtain a similar conclusion as for the AP simulations. The impact of the correction on the global conservation of AM is modest, removing only about 15 % of the sink at f19 resolution. However, its action is stronger on upper-level winds (see Fig. 6b), which leads to proportionally reduced fixer increments at those levels and thus to smaller impacts by the fixer on areas affected by baroclinic instability.

Figures 10 and S3 in the Supplement show the seasonally resolved impacts on the zonal-mean zonal winds from applying the combination of the fixer and correction in F2000 simulations at both f19 and f09 resolutions (see also Fig. S3 in the Supplement for JJA). In all cases, the reduction of biases in both easterly and westerly wind regimes is noticeable, the latter especially at the subpolar latitudes of the winter hemisphere.

More in detail, it may be noted that the benefits of the AM modifications appear more clearly for the winds in the simulation at the lower resolution, for which the numerical sink of AM is indeed larger. These benefits, however, are not limited to the zonal-mean zonal winds, and they are also appreciable at the f09 resolution. Most notable is the reduction in the strength of the Hadley circulations (see Fig. S4 in the Supplement), which is expected from the arguments set out in the Introduction. This has consequences for many aspects of the global circulation. Figure 11 shows a summary of the impacts on the quality of the simulations in relation to the observed climatology. The improvements at f09 seem particularly remarkable considering that the unmodified simulation is a scientifically supported case that has been fully tuned for a best match to observations. It may be noted that no additional tuning whatsoever is involved in the simulation with AM modifications shown here and that the AM modifications themselves have no free parameters as they follow directly from an effort to reduce the numerical sink stemming from the FV dynamical core. The better quality of this simulation thus follows entirely from better adherence of the solution to a fundamental property of the equations of motion.

5 Summary and conclusions

AM conservation in CAM–FV has been substantially improved by means of a correction that reduces the zonal-mean numerical sink of Lin and Rood's (1997) shallow-water scheme and a fixer that ensures the conservation of global angular momentum under advection. The effectiveness of these modifications in terms of AM conservation in the simulations presented here is summarised in Table 1. We show that aside from global AM conservation, they have other significant impacts on the simulations, consistent with the “tea leaves” mechanism (Einstein, 1926) that rapidly redistributes pressure forces in a rotating fluid in response to zonal accelerations. The most notable effect is a reduction of the excessive easterlies of the model, with a concomitant slowdown of the Hadley circulation. As a result of such changes, the simulations of the observed climatology show marked improvements.

The zonal-mean correction of the shallow-water scheme is not necessary for enforcing global conservation, as this can be achieved by the fixer alone. Indeed, the correction is quite ineffective in realistic simulations of the atmosphere in terms of global conservation. Nevertheless, we find that its concomitant application with the fixer has positive impacts on the simulations. In particular, it reduces the effects of the fixer in the mid-latitudes. This can be explained with the greater effectiveness of the correction in the baroclinically unstable regions around the subtropical jet streams, where the zonal-mean numerical sink appears to be largest. Even so, because of its potentially large local effects, the utilisation of the correction under different set-ups should be tested on a case-by-case basis according to its impacts on the results.

Improving the quality of the simulation of the global distribution of surface wind stress should be expected to bring particular benefits to coupled atmosphere–ocean simulations. An adequate discussion of such a coupled simulation would exceed the scope of the present paper, which is aimed primarily at presenting the method. In particular, due to their computational expense, at the present time it is not possible to produce well spun-up coupled simulations that can provide an assessment of the impact of the AM modifications.

Table 1 Simulation set-ups and the effect of AM modifications. The percentage figures represent the numerical source (negative for sink) of global total atmospheric AM relative to the global total physical eastward torque acting on the atmosphere (terms T_x and C_λ in Eq. (A1), when only the positive parts of the integrands are summed). The column “Experiments” indicates which modifications to CAM–FV are used (the relevant sections of this paper are indicated in the footnotes). The three columns under “Simulations” are for results obtained with model integrations in Held–Suarez mode (Held and Suarez, 1994), in aquaplanet mode (Neale and Hoskins, 2000), and in F2000 mode, i.e. an AMIP-type (Gates, 1992) simulation with annually repeating present-day climatological SSTs.

^a Sections 2.2 and 2.1. ^b Section 2.3. ^c Sections 2.3 and 2.4.

The modification to the FV dynamical core that we describe and utilise is relatively crude and causes local loss of accuracy due to violation of vorticity conservation under advection. Nevertheless, the associated detrimental impacts appear to be fairly limited, with insignificant differences under standard tests such as the Jablonowski and Williamson (2006) baroclinic wave test, which should be sensitive to local conservation. Even so, it is clear from the very same tests that simulations over weather timescales are not sensitive to AM conservation, so for such an application it is not advisable to trade enforcing such conservation for a loss of accuracy. On the longer timescales of climate simulations, by contrast, our results demonstrate the importance of the global conservation of atmospheric AM in order to obtain a realistic global circulation.

Appendix A: Offline diagnostics of numerical torque in model simulations

The diagnosis of the residual torque that violates AM conservation in CAM simulations follows from the hydrostatic primitive equations (see White et al., 2005). In our zonally and vertically integrated diagnostics such as in Fig. 1 the AM source is calculated as

$$\begin{array}{}\text{(A1)}& S_M={\partial }_t{L}_r+D_L-T_x-C_{\mathit{\lambda }},\end{array}$$

where the first term on the right-hand side represents the tendency of relative atmospheric AM, the second term represents the divergence of the flux of relative AM, the third the external torque (which in all simulations presented in Sects. 1, 2, and 3, when non-vanishing, is exclusively due to surface stresses or linear friction in the planetary boundary layer – PBL), and the last term is the tendency of planetary atmospheric AM due to the vertically integrated divergence of atmospheric mass. The formulas are as follows.

$$\begin{array}{c}{\displaystyle }{L}_r=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{p_{*}}{\overset{p_{\mathrm{top}}}{\int }}\left(ua\cos \mathit{\phi }\right){\displaystyle \frac{\mathrm{d}p}{g}}a\cos \mathit{\phi }\mathrm{d}\mathit{\lambda }\\ {\displaystyle }{D}_L={\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{p_{*}}{\overset{p_{\mathrm{top}}}{\int }}\left(uva\cos \mathit{\phi }\right){\displaystyle \frac{\mathrm{d}p}{g}}a\cos \mathit{\phi }\mathrm{d}\mathit{\lambda }\\ {\displaystyle }{T}_x=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left({\mathit{\tau }}_{x}a\cos \mathit{\phi }\right)a\cos \mathit{\phi }\mathrm{d}\mathit{\lambda }\\ {\displaystyle }{C}_{\mathit{\lambda }}=-{\displaystyle \frac{a\mathrm{\Omega }\sin \mathrm{2}\mathit{\phi }}{g}}{\partial }_t\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\phi }}{\int }}{p}_{*}{a}^{\mathrm{2}}\cos {\mathit{\phi }}^{\prime }\mathrm{d}{\mathit{\phi }}^{\prime }\mathrm{d}\mathit{\lambda }\end{array}$$

Here, a is the Earth's radius, φ the latitude, λ the longitude, g the gravitational acceleration in Earth's surface, Ω the angular speed of Earth's rotation, and u, v, p_*, and τ_x are the zonal-wind component, the meridional wind component, the surface pressure, and the zonal component of the surface or frictional stress acting on the air in the model simulations. Note that to obtain C_λ the continuity equation was used. Note that for the time-average values of S_M, the time differentials become increments between the initial and the final state; terms T_x and C_λ are linear in the wind stress and the surface pressure, respectively. Terms L_r and D_L are bilinear and trilinear in the model prognostic quantities u, v, and p_*, so an online computation of the time averages of the integrands is required for these terms. CAM provides time-mean diagnostics of the zonal wind u and of the product of the wind components uv conservatively interpolated onto standard pressure levels, and the integrals in Eq. (A1) are computed with their help.

Appendix B: Formulation and approximations for the AM correction in CAM–FV

The local conservation equation for the shallow-water equations is

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& {\partial }_t\left[\mathrm{\Delta }p\left(ua\cos \mathit{\phi }+\mathrm{\Omega }{a}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\phi }\right)\right]=\\ & -{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left[\mathrm{\Delta }p\left(ua\cos \mathit{\phi }+\mathrm{\Omega }{a}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\phi }\right)v\cos \mathit{\phi }\right]\\ & -{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}\left[\mathrm{\Delta }p\left(ua\cos \mathit{\phi }+\mathrm{\Omega }{a}^{\mathrm{2}}{\cos }^{\mathrm{2}}\mathit{\phi }\right)u\right],\end{array}\end{array}$$

where φ and λ are latitude and longitude, respectively, Δp is the layer thickness in terms of hydrostatic pressure, u and v are the zonal and meridional wind components, a is the Earth's radius, and Ω the Earth's angular velocity. Note that we are ignoring pressure and geopotential terms here, as we focus exclusively on the process of advection. Accordingly, Δp, i.e. the layer under consideration, may be arbitrary, except that it satisfies the shallow-water mass conservation equation; i.e. we follow Lin's (2004) “vertically Lagrangian” approach by following the vertical motion of the layer. Integrating Eq. (B1) over longitude, we obtain

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}\int \mathrm{d}\mathit{\lambda }{\partial }_t\left(\mathrm{\Delta }pua{\cos }^{\mathrm{2}}\mathit{\phi }\right)& =-\int \mathrm{d}\mathit{\lambda }{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)\\ & +\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }pfva{\cos }^{\mathrm{2}}\mathit{\phi },\end{array}\end{array}$$

where f is the Coriolis parameter. To address the FV scheme's violation of this conservation, we apply an additional, zonally uniform increment of the zonal wind, $\stackrel{\mathrm{‾}}{\mathit{\delta }u}$, such that over each shallow-water sub-step δt (we shall refer to this simply as the “time step” in this section) of the dynamical core,

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathit{\delta }t}}\int \mathrm{d}\mathit{\lambda }\cos \mathit{\phi }\left[\mathrm{\Delta }{p}_{\mathrm{n}}\left(u_{\mathrm{n}}+\stackrel{\mathrm{‾}}{\mathit{\delta }u}\right)-\mathrm{\Delta }{p}_{\mathrm{o}}{u}_{\mathrm{o}}\right]\cos \mathit{\phi }=\\ & -\int \mathrm{d}\mathit{\lambda }\cos \mathit{\phi }{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)\\ & +\int \mathrm{d}\mathit{\lambda }{\cos }^{\mathrm{2}}\mathit{\phi }\mathrm{\Delta }pfv.\end{array}\end{array}$$

Here, “old” prognostic quantities (i.e. valid at the beginning of the time step) and “new” prognostic quantities (i.e. valid at the end of the time step, before any correction) are indicated by the subscripts “o” and “n”, respectively; quantities without subscripts are intended to be time-centred, representing advective fluxes over the time step. To obtain the correction, we solve this equation for the required increment $\stackrel{\mathrm{‾}}{\mathit{\delta }u}$ and substitute for u_n the actual FV zonal-wind increment over the time step:

$$\begin{array}{}\text{(B4)}& u_{\mathrm{n}}=u_{\mathrm{o}}+\left({\mathit{\xi }}_{\mathrm{o}}v-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}K\right)\mathit{\delta }t,\end{array}$$

where ξ is the absolute vorticity, and K is the kinetic energy term as discretised in LR97's scheme. The result is

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& \left(\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }{p}_{\mathrm{n}}\right)\stackrel{\mathrm{‾}}{\mathit{\delta }u}=-\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }{p}_{\mathrm{n}}\left({\mathit{\zeta }}_{\mathrm{o}}v-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}K\right)\mathit{\delta }t\\ & -\int \mathrm{d}\mathit{\lambda }\left(\mathrm{\Delta }{p}_{\mathrm{n}}-\mathrm{\Delta }{p}_{\mathrm{o}}\right)\left[u_{\mathrm{o}}+\left({\mathit{\xi }}_{\mathrm{o}}v-{\mathit{\zeta }}_{\mathrm{o}}v\right)\mathit{\delta }t\right]\\ & -\int \mathrm{d}\mathit{\lambda }{\displaystyle \frac{\mathrm{1}}{a{\cos }^{\mathrm{2}}\mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)\mathit{\delta }t.\end{array}\end{array}$$

The term in the second line on the right-hand side representing the advection of planetary vorticity is written in a roundabout way for later convenience.

We note two aspects of this expression. First, there is a significant numerical cancellation between the second and the third lines on the right-hand side. Second, all advective terms in the first two lines on the right-hand side can be easily discretised according to the standard LR97 prescription and are thus automatically defined on D-grid zonal velocity points, i.e. where required for $\stackrel{\mathrm{‾}}{\mathit{\delta }u}$. However, all mass factors are defined on scalar points, i.e. on the A grid. Furthermore, the integrand in the third line on the right-hand side has no natural expression in LR97's discretisation, and both zonal and meridional winds in that expression need to be interpolated onto the A grid. Hence, additional interpolation is required for these terms. Notwithstanding these issues, we found that this correction, when implemented, gave accurate conservation of AM. However, it also proved to cause numerical instability such that the integration crashed within seven or eight time steps. Analysis suggested that the last term on the right-hand side had to be recast in a different form.

We therefore chose to approximate the last term as follows:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{a{\cos }^{\mathrm{2}}\mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }puv{\cos }^{\mathrm{2}}\mathit{\phi }\right)& \approx \left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }pv\cos \mathit{\phi }\right)\right]u\\ & +\left[{\displaystyle \frac{v}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(u\cos \mathit{\phi }\right)\right]\mathrm{\Delta }p.\end{array}\end{array}$$

The approximation here consists of using C-grid (advective) fluxes in the partial differentials on the right-hand side. Considering this as a calculation for the advective fluxes of zonal momentum, which is its physical meaning, this appears to be a valid interpretation for v. For the values of Δp and u outside the operators, we adopt the substitutions

$$\begin{array}{ll}{\displaystyle }u& {\displaystyle }=:u_{\mathrm{o}}+{\mathit{\delta }}_{h}u+{\mathit{\delta }}^{\prime \prime }u\\ {\displaystyle }\mathrm{\Delta }p& {\displaystyle }=:\mathrm{\Delta }{p}_{\mathrm{n}}-{\mathit{\delta }}_h\mathrm{\Delta }p+{\mathit{\delta }}^{\prime \prime }\mathrm{\Delta }p,\end{array}$$

where

$$\begin{array}{}\text{(B7)}& {\mathit{\delta }}_h\mathrm{\Delta }p:={\displaystyle \frac{\mathrm{\Delta }{p}_{\mathrm{n}}-\mathrm{\Delta }{p}_{\mathrm{o}}}{\mathrm{2}}},{\mathit{\delta }}_{h}u:={\displaystyle \frac{u_{\mathrm{n}}-u_{\mathrm{o}}}{\mathrm{2}}},\end{array}$$

and δ^′′u and δ^′′Δp are formally o(δt). The increments are still understood as advective only, i.e. they exclude pressure force terms. By further using the identities

$$\begin{array}{}\text{(B8)}& {\displaystyle }-{\displaystyle \frac{\mathit{\delta }t}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(\mathrm{\Delta }pv\cos \mathit{\phi }\right)=\mathrm{\Delta }{p}_{\mathrm{n}}-\mathrm{\Delta }{p}_{\mathrm{o}}+{\displaystyle \frac{\mathit{\delta }t}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}\left(\mathrm{\Delta }pu\right),\text{(B9)}& {\displaystyle }-\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\phi }}\left(u_{\mathrm{o}}\cos \mathit{\phi }\right)\right]v\mathit{\delta }t=\left({\mathit{\zeta }}_{\mathrm{o}}-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}{v}_{\mathrm{o}}\right)v\mathit{\delta }t,\end{array}$$

we finally arrive at the expression for our approximate angular-momentum-conserving zonal-mean zonal-wind correction:

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}& \left(\int \mathrm{d}\mathit{\lambda }\mathrm{\Delta }{p}_{\mathrm{n}}\right)\stackrel{\mathrm{‾}}{\mathit{\delta }u}=\int \mathrm{d}\mathit{\lambda }\left(\mathrm{\Delta }{p}_{\mathrm{n}}-{\mathit{\delta }}_h\mathrm{\Delta }p\right)\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}K-{\mathit{\zeta }}_{\mathit{\lambda }o}v\right]\mathit{\delta }t\\ & +\int \mathrm{d}\mathit{\lambda }\left[{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}\left(\mathrm{\Delta }pu\right)\mathit{\delta }t\right]\left(u_{\mathrm{o}}+{\mathit{\delta }}_{h}u\right)\\ & +\int \mathrm{d}\mathit{\lambda }\left[\mathrm{2}{\mathit{\delta }}_h\mathrm{\Delta }p+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\phi }}}{\partial }_{\mathit{\lambda }}\left(\mathrm{\Delta }pu\right)\mathit{\delta }t\right]{\mathit{\delta }}^{\prime \prime }u\\ & +\int \mathrm{d}\mathit{\lambda }{\mathit{\delta }}^{\prime \prime }\mathrm{\Delta }p\left[{\mathit{\xi }}_{\mathrm{o}}v-{\mathit{\zeta }}_{\mathit{\lambda }o}v\right]\mathit{\delta }t,\end{array}\end{array}$$

where we have used the shorthand ${\mathit{\zeta }}_{\mathit{\lambda }o}:=\frac{\mathrm{1}}{a\cos \mathit{\phi }}{\partial }_{\mathit{\lambda }}{v}_{\mathrm{o}}$.

We note that setting the higher-order terms to zero implies that the correction has no effect on a zonally symmetric flow. If, in addition, the flow is in an exact steady state, then the correction always vanishes identically, regardless of these terms. It can further be shown that, if the term in K is the true gradient of the kinetic energy in the original scheme, for any values of δ^′′u and δ^′′Δp that are first order in δt or higher, the correction Eq. (B10) is formally third order in δt or higher. In other words, the correction will not affect solutions that are already locally angular momentum conserving.

In Eq. (B10), all mass terms must be averaged over φ; by contrast, all advective terms (in square brackets) represent fluxes as discretised according to the standard LR97 algorithm. The discretised expression of Eq. (B10) thus corresponds to Eq. (7). The only additional PPM calculation required to calculate this correction is the meridional advection of the partial relative vorticity, ζ_λ, with a minimal additional computational cost that is hardly detectable in CAM simulations.

Appendix C: Formulation and implementation of the AM fixer in CAM–FV

As we explain in Sect. 2.4, the fixer is based on diagnosing the global change in atmospheric AM due to advective increments only, which should vanish identically according to the continuous equations. When applied, the fixer counteracts that change at every advective sub-step; irrespectively, its time-mean increments can always be used to diagnose AM non-conservation in the simulations in a manner that is completely independent of the physics parameterisations or boundary conditions used and hence independent of the particular configuration of the simulations itself. All the calculations related to the fixer and the quantification of the numerical (advective) AM source are internal to the dynamical core only, indeed of its shallow-water part.

So, for each time step and at each level k, we require the advective shallow-water equation increments to satisfy

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}& \mathit{\delta }\mathit{\left\{}\sum _{i,j}\right[u_{i,j}\cos e_j+u_{i,j+\mathrm{1}}\cos e_{j+\mathrm{1}}\\ & +a\mathrm{\Omega }\left({\cos }^{\mathrm{2}}{e}_j+{\cos }^{\mathrm{2}}{e}_{j+\mathrm{1}}\right)]\cos c_j\mathrm{\Delta }{p}_{i,j}{\mathit{\}}}_k=\mathrm{0}\end{array},\end{array}$$

where the indices i and j refer to longitude and latitude, respectively, e_j represents the latitudes of the zonal velocity points of the D grid, and c_j represents the latitudes of the scalar points (A grid). The other symbols have the same meaning as in the previous section, and δ represents the purely advective increment obtained in the dynamical core, which may include the correction discussed above. The action of the fixer in this context is represented by an additional increment δϖ_k so that the total increment of the zonal wind becomes $\mathit{\delta }{u}_{i,j,k}+a\mathit{\delta }{\mathit{\varpi }}_k\cos e_j$. We obtain the following:

$$\begin{array}{}\text{(C2)}& \mathit{\delta }{\mathit{\varpi }}_k=-{\displaystyle \frac{T_k}{I_k}},\end{array}$$

where the numerical torque is

$$\begin{array}{}\text{(C3)}& \begin{array}{rl}{T}_k& =a\sum _{i,j}\cos e_j\left(\cos c_j+\cos c_{j-\mathrm{1}}\right)\mathit{\left\{}\mathit{\delta }{u}_{i,j}{\stackrel{\mathrm{‾}}{\mathrm{\Delta }{p}_{i,j}}}^{\mathit{\phi }}\right(t+\mathrm{\Delta }t)\\ & +\left[u_{i,j}\left(t\right)+a\mathrm{\Omega }\cos e_j\right]\mathit{\delta }{\stackrel{\mathrm{‾}}{\mathrm{\Delta }p}}_{i,j}^{\mathit{\phi }}{\mathit{\}}}_k\end{array},\end{array}$$

and the moment of inertia is

$$\begin{array}{}\text{(C4)}& I_k=a^{\mathrm{2}}\sum _{i,j}{\cos }^{\mathrm{2}}{e}_j\left(\cos c_j+\cos c_{j-\mathrm{1}}\right){\stackrel{\mathrm{‾}}{\mathrm{\Delta }p}}_{i,j,k}^{\mathit{\phi }}(t+\mathrm{\Delta }t).\end{array}$$

In these expressions,

$$\begin{array}{}\text{(C5)}& {\stackrel{\mathrm{‾}}{\mathrm{\Delta }p}}_{i,j,k}^{\mathit{\phi }}:={\displaystyle \frac{\mathrm{\Delta }{p}_{i,j,k}\cos c_j+\mathrm{\Delta }{p}_{i,j-\mathrm{1},k}\cos c_{j-\mathrm{1}}}{\cos c_j+\cos c_{j-\mathrm{1}}}}.\end{array}$$

Equation (C2) gives the required angular acceleration of the entire atmospheric shell at model level k. The action of the level fixer is therefore to add an increment to the zonal wind:

$$\begin{array}{}\text{(C6)}& {\mathit{\delta }}^f{u}_{i,j,k}=a\mathit{\delta }{\mathit{\varpi }}_k\cos e_j.\end{array}$$

In some regions of the model domain, it is not desirable to apply a fixer, since dissipation is explicitly built into in the dynamical core formulation. This is the case near the upper boundary of CAM's domain (the lower boundary in pressure space), where the fixer is accordingly switched off. In general, a weight w_k≤1 can be applied at each level so that Eq. (C2) becomes

$$\begin{array}{}\text{(C7)}& \mathit{\delta }{\mathit{\varpi }}_k=-w_k{\displaystyle \frac{T_k}{I_k}},\end{array}$$

where only a fraction w_k of the numerical torque at level k is compensated for by the fixer at that level.

The global fixer applies the same solid-body rotation increment to all levels within the domain where it is required. When all weights are unity, this is simply

$$\begin{array}{}\text{(C8)}& \mathit{\delta }{\mathit{\varpi }}_g=-{\displaystyle \frac{{\sum }_i{T}_i}{{\sum }_j{I}_j}};\end{array}$$

when $\exists k:w_k<\mathrm{1}$, the vertical integrals must be weighted accordingly and the weights applied to the correction at each level so that

$$\begin{array}{}\text{(C9)}& \mathit{\delta }{\mathit{\varpi }}_{g,k}=-w_k{\displaystyle \frac{{\sum }_i{w}_i{T}_i}{{\sum }_j{w}_j{I}_j}}.\end{array}$$

It can be seen that ${\sum }_k{I}_k\mathit{\delta }{\mathit{\varpi }}_{g,k}=-{\sum }_k{w}_k{T}_k$ so that the numerical torque associated with the domain of interest is also fully compensated for by this fixer. Experimentation has shown that tapering the global fixer so as to exclude its action from levels in the stratosphere was necessary in order to avoid distortions of the dynamics in layers where it is sensitive to small amounts of zonal acceleration and where, moreover, thanks to the predominance of solenoidal dynamics (before gravity-wave drag, which is applied in the physics parameterisations), the dynamical core performs well in terms of AM conservation. For the latter reason, no tapering (i.e. any weights other than 1 in the valid domain and 0 in the filtered layers near the model lid) is in fact required for the level fixer.

For diagnostic purposes, fixer increments are always calculated as in Eq. (C2) and provided in output. Use of the increments in Eq. (C2) leads to the conservation of total AM in idealised spin-up or spin-down experiments with no physical sources or sinks of momentum (see Fig. 3), as well as an accurate balance of the surface torques in Held–Suarez or aquaplanet simulations wherein only surface stresses are present (and accurately diagnosed). Hence, we obtain two important conclusions. First, all numerical sources of AM indeed reside in the advective wind increments of the shallow-water part of the dynamical core; second, the fixer diagnostics return an accurate estimate of the apparent numerical AM source for any CAM–FV integration, irrespective of physics parameterisations or boundary fluxes (including orographic form drag).