- Articles & preprints
- Submission
- Policies
- Peer-review process
- Editorial board
- About
- EGU publications
- Manuscript tracking

Journal cover
Journal topic
**Geoscientific Model Development**
An interactive open-access journal of the European Geosciences Union

Journal topic

- Articles & preprints
- Submission
- Policies
- Peer-review process
- Editorial board
- About
- EGU publications
- Manuscript tracking

- Articles & preprints
- Submission
- Policies
- Peer-review process
- Editorial board
- About
- EGU publications
- Manuscript tracking

- Abstract
- Introduction
- Analysis of potential causes and approaches to correction
- Numerical simulations and results
- Simulations of the observed climatology
- Summary and conclusions
- Appendix A: Offline diagnostics of numerical torque in model simulations
- Appendix B: Formulation and approximations for the AM correction in CAM–FV
- Appendix C: Formulation and implementation of the AM fixer in CAM–FV
- Code and data availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement

GMD | Articles | Volume 13, issue 2

Geosci. Model Dev., 13, 685–705, 2020

https://doi.org/10.5194/gmd-13-685-2020

© Author(s) 2020. This work is distributed under

the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/gmd-13-685-2020

© Author(s) 2020. This work is distributed under

the Creative Commons Attribution 4.0 License.

Special issue: The Norwegian Earth System Model: NorESM; basic development,...

**Development and technical paper**
21 Feb 2020

**Development and technical paper** | 21 Feb 2020

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

^{1}NORCE Klima and Bjerknes Centre for Climate Research, Bergen, Norway^{2}Department of Meteorology (MISU), Stockholm University, Stockholm, Sweden^{3}National Center for Atmospheric Research, Boulder, Colorado, USA^{4}Department of Climate and Space Sciences and Engineering, University of Michigan, Ann Arbor, Michigan, USA

^{1}NORCE Klima and Bjerknes Centre for Climate Research, Bergen, Norway^{2}Department of Meteorology (MISU), Stockholm University, Stockholm, Sweden^{3}National Center for Atmospheric Research, Boulder, Colorado, USA^{4}Department of Climate and Space Sciences and Engineering, University of Michigan, Ann Arbor, Michigan, USA

**Correspondence**: Thomas Toniazzo (thomas.toniazzo@uni.no)

**Correspondence**: Thomas Toniazzo (thomas.toniazzo@uni.no)

Abstract

Back to toptop
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.

How to cite

Back to top
top
How to cite.

Toniazzo, T., Bentsen, M., Craig, C., Eaton, B. E., Edwards, J., Goldhaber, S., Jablonowski, C., and Lauritzen, P. H.: Enforcing conservation of axial angular momentum in the atmospheric general circulation model CAM6, Geosci. Model Dev., 13, 685–705, https://doi.org/10.5194/gmd-13-685-2020, 2020.

1 Introduction

Back to toptop
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^{−2} (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.

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

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

Back to toptop
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.

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(\phantom{\rule{0.125em}{0ex}}\mathit{\alpha}\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\lambda}}p\phantom{\rule{0.125em}{0ex}}+\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\lambda}}\mathit{\varphi}\phantom{\rule{0.125em}{0ex}}\right)}_{k}& =-{\left({\displaystyle \frac{\mathrm{\Delta}\phantom{\rule{-0.125em}{0ex}}\mathit{\varphi}}{\mathrm{\Delta}\phantom{\rule{-0.125em}{0ex}}p}}\right)}_{k}\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\lambda}}{p}_{k-\mathrm{1}/\mathrm{2}}\\ & +\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\lambda}}{\mathit{\varphi}}_{k+\mathrm{1}/\mathrm{2}}\phantom{\rule{0.125em}{0ex}}+\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta}\phantom{\rule{-0.125em}{0ex}}{p}_{k}}}{\partial}_{\mathit{\lambda}}\left[{a}_{k}(\mathit{\alpha}p{)}_{k}\mathrm{\Delta}\phantom{\rule{-0.125em}{0ex}}{p}_{k}\right],\end{array}\end{array}$$

where the symbol Δ is employed to represent a difference
between vertical levels, $\mathrm{\Delta}\phantom{\rule{-0.125em}{0ex}}{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}\phantom{\rule{0.125em}{0ex}}\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}\phantom{\rule{-0.125em}{0ex}}{p}_{k}\right\}-\mathrm{\Delta}{\left(\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}{\mathit{\delta}}_{\mathit{\lambda}}p\right)}_{k},\end{array}$$

where *δ*_{λ} is the finite-difference operator in the zonal
direction, and $\stackrel{\mathrm{\u203e}}{{\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}}}\phantom{\rule{0.25em}{0ex}},\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}={\displaystyle \frac{{\mathit{\varphi}}_{i+\mathrm{1}/\mathrm{2}}+{\mathit{\varphi}}_{i-\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}},\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.

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.

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 \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\partial}_{t}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}ua{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)& =-\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\\ & +\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}fva{\mathrm{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{\u203e}}{u}& ={\displaystyle \frac{\mathrm{1}}{\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}p}}\times \mathit{\{}\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}\left({\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}K-\mathit{\zeta}v\right)\\ & -\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{\mathrm{1}}{a{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\\ & -\int \phantom{\rule{0.125em}{0ex}}\mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}u{\partial}_{t}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}\mathit{\}}\phantom{\rule{0.125em}{0ex}}.\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{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)=\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{\mathrm{1}}{a{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}}}\left[\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}v{\partial}_{\mathit{\phi}}\left(u\mathrm{cos}\mathit{\phi}\right)+u{\partial}_{\mathit{\phi}}\left(v\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}\mathrm{cos}\mathit{\phi}\right)\right]\end{array}\phantom{\rule{0.125em}{0ex}},\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{\u203e}}{u}& ={\displaystyle \frac{\mathrm{1}}{\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}{p}_{t+\mathit{\delta}t}}}}\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathit{\left\{}\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}p}\right[{\displaystyle \frac{\mathit{\delta}t}{a\mathrm{cos}\mathit{\phi}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}\mathit{\lambda}}}{\mathit{\delta}}_{\mathit{\lambda}}\phantom{\rule{-0.125em}{0ex}}K\\ & -\stackrel{\mathrm{\u203e}}{\mathcal{Y}\left({v}^{*},\mathit{\delta}t;{\mathit{\zeta}}_{\mathit{\lambda}}\right)}]+{\stackrel{\mathrm{\u203e}}{u}}^{t}\mathcal{F}\left({u}^{*},\mathit{\delta}t;\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}p}\right)+O\left(\mathit{\delta}{t}^{\mathrm{2}}\right)\mathit{\}}\phantom{\rule{0.125em}{0ex}}.\end{array}\end{array}$$

Here, ${\mathit{\zeta}}_{\mathit{\lambda}}:=\frac{\mathrm{1}}{a\mathrm{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”.

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

Back to toptop
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.

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

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.

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

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.

4 Simulations of the observed climatology

Back to toptop
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^{1}.

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

Back to toptop
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.

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

Back to toptop
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\mathrm{cos}\mathit{\phi}\right){\displaystyle \frac{\mathrm{d}p}{g}}a\mathrm{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\mathrm{cos}\mathit{\phi}\right){\displaystyle \frac{\mathrm{d}p}{g}}a\mathrm{cos}\mathit{\phi}\mathrm{d}\mathit{\lambda}\\ {\displaystyle}{T}_{x}=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi}}{\int}}\left({\mathit{\tau}}_{x}a\mathrm{cos}\mathit{\phi}\right)a\mathrm{cos}\mathit{\phi}\mathrm{d}\mathit{\lambda}\\ {\displaystyle}{C}_{\mathit{\lambda}}=-{\displaystyle \frac{a\mathrm{\Omega}\mathrm{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}}\mathrm{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 *u**v*
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

Back to toptop
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\mathrm{cos}\mathit{\phi}+\mathrm{\Omega}{a}^{\mathrm{2}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\right]=\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}-{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left[\mathrm{\Delta}p\left(ua\mathrm{cos}\mathit{\phi}+\mathrm{\Omega}{a}^{\mathrm{2}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)v\mathrm{cos}\mathit{\phi}\right]\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}-{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}\left[\mathrm{\Delta}p\left(ua\mathrm{cos}\mathit{\phi}+\mathrm{\Omega}{a}^{\mathrm{2}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)u\right]\phantom{\rule{0.125em}{0ex}},\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}\phantom{\rule{0.125em}{0ex}}{\partial}_{t}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}ua{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)& =-\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\\ & +\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}fva{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\phantom{\rule{0.125em}{0ex}},\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{\u203e}}{\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}\phantom{\rule{0.125em}{0ex}}\mathrm{cos}\mathit{\phi}\left[\mathrm{\Delta}{p}_{\mathrm{n}}\left({u}_{\mathrm{n}}+\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}\right)-\mathrm{\Delta}{p}_{\mathrm{o}}{u}_{\mathrm{o}}\right]\mathrm{cos}\mathit{\phi}=\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}-\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{cos}\mathit{\phi}{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}puv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}fv\phantom{\rule{0.125em}{0ex}}.\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{\u203e}}{\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\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}K\right)\mathit{\delta}t\phantom{\rule{0.125em}{0ex}},\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}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{p}_{\mathrm{n}}\right)\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}=-\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{p}_{\mathrm{n}}\phantom{\rule{0.125em}{0ex}}\left({\mathit{\zeta}}_{\mathrm{o}}v-{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}K\right)\mathit{\delta}t\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}-\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{\Delta}{p}_{\mathrm{n}}\phantom{\rule{0.125em}{0ex}}-\mathrm{\Delta}{p}_{\mathrm{o}}\phantom{\rule{0.125em}{0ex}}\right)\left[{u}_{\mathrm{o}}+\left({\mathit{\xi}}_{\mathrm{o}}v-{\mathit{\zeta}}_{\mathrm{o}}v\right)\mathit{\delta}t\right]\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}-\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{\mathrm{1}}{a{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)\mathit{\delta}t\phantom{\rule{0.125em}{0ex}}.\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{\u203e}}{\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{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}uv{\mathrm{cos}}^{\mathrm{2}}\mathit{\phi}\right)& \approx \left[{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}v\mathrm{cos}\mathit{\phi}\right)\right]u\\ & +\left[{\displaystyle \frac{v}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(u\mathrm{cos}\mathit{\phi}\right)\right]\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}.\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}\phantom{\rule{-0.125em}{0ex}}u\\ {\displaystyle}\mathrm{\Delta}p& {\displaystyle}=:\mathrm{\Delta}{p}_{\mathrm{n}}-{\mathit{\delta}}_{h}\mathrm{\Delta}p+{\mathit{\delta}}^{\prime \prime}\phantom{\rule{-0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}},\end{array}$$

where

$$\begin{array}{}\text{(B7)}& {\mathit{\delta}}_{h}\mathrm{\Delta}p:={\displaystyle \frac{\mathrm{\Delta}{p}_{\mathrm{n}}\phantom{\rule{0.125em}{0ex}}-\mathrm{\Delta}{p}_{\mathrm{o}}}{\mathrm{2}}}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}{\mathit{\delta}}_{h}u:={\displaystyle \frac{{u}_{\mathrm{n}}-{u}_{\mathrm{o}}}{\mathrm{2}}}\phantom{\rule{0.125em}{0ex}},\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\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}v\mathrm{cos}\mathit{\phi}\right)=\mathrm{\Delta}{p}_{\mathrm{n}}\phantom{\rule{0.125em}{0ex}}-\mathrm{\Delta}{p}_{\mathrm{o}}+{\displaystyle \frac{\mathit{\delta}t}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}u\right),\text{(B9)}& {\displaystyle}-\left[{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\phi}}\left({u}_{\mathrm{o}}\mathrm{cos}\mathit{\phi}\right)\right]v\mathit{\delta}t=\left({\mathit{\zeta}}_{\mathrm{o}}-{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}{v}_{\mathrm{o}}\right)v\mathit{\delta}t\phantom{\rule{0.125em}{0ex}},\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}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{p}_{\mathrm{n}}\right)\stackrel{\mathrm{\u203e}}{\mathit{\delta}u}=\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\left(\mathrm{\Delta}{p}_{\mathrm{n}}-{\mathit{\delta}}_{h}\mathrm{\Delta}p\right)\left[{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}K-{\mathit{\zeta}}_{\mathit{\lambda}o}v\right]\mathit{\delta}t\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\left[{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}u\right)\mathit{\delta}t\right]\left({u}_{\mathrm{o}}+{\mathit{\delta}}_{h}u\right)\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\left[\mathrm{2}{\mathit{\delta}}_{h}\mathrm{\Delta}p+{\displaystyle \frac{\mathrm{1}}{a\mathrm{cos}\mathit{\phi}}}{\partial}_{\mathit{\lambda}}\left(\mathrm{\Delta}p\phantom{\rule{0.125em}{0ex}}u\right)\mathit{\delta}t\right]\phantom{\rule{0.125em}{0ex}}{\mathit{\delta}}^{\prime \prime}\phantom{\rule{-0.125em}{0ex}}u\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+\int \mathrm{d}\mathit{\lambda}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}{\mathit{\delta}}^{\prime \prime}\phantom{\rule{-0.125em}{0ex}}\mathrm{\Delta}p\phantom{\rule{0.25em}{0ex}}\left[{\mathit{\xi}}_{\mathrm{o}}v-{\mathit{\zeta}}_{\mathit{\lambda}o}v\right]\mathit{\delta}t\phantom{\rule{0.125em}{0ex}},\end{array}\end{array}$$

where we have used the shorthand ${\mathit{\zeta}}_{\mathit{\lambda}o}:=\frac{\mathrm{1}}{a\mathrm{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

Back to toptop
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}\mathrm{cos}{e}_{j}+{u}_{i,j+\mathrm{1}}\mathrm{cos}{e}_{j+\mathrm{1}}\\ & \phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}+a\mathrm{\Omega}\left({\mathrm{cos}}^{\mathrm{2}}{e}_{j}+{\mathrm{cos}}^{\mathrm{2}}{e}_{j+\mathrm{1}}\right)]\mathrm{cos}{c}_{j}\phantom{\rule{0.125em}{0ex}}\mathrm{\Delta}{p}_{i,j}\phantom{\rule{0.125em}{0ex}}{\mathit{\}}}_{k}=\mathrm{0}\end{array}\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}},\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}\mathrm{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}\mathrm{cos}{e}_{j}\left(\mathrm{cos}{c}_{j}+\mathrm{cos}{c}_{j-\mathrm{1}}\right)\mathit{\left\{}\mathit{\delta}{u}_{i,j}{\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}{p}_{i,j}}}^{\mathit{\phi}}\right(t+\mathrm{\Delta}t)\\ & +\left[{u}_{i,j}\left(t\right)+a\mathrm{\Omega}\mathrm{cos}{e}_{j}\right]\mathit{\delta}{\stackrel{\mathrm{\u203e}}{\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}{\mathrm{cos}}^{\mathrm{2}}{e}_{j}\left(\mathrm{cos}{c}_{j}+\mathrm{cos}{c}_{j-\mathrm{1}}\right){\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}p}}_{i,j,k}^{\mathit{\phi}}(t+\mathrm{\Delta}t)\phantom{\rule{0.125em}{0ex}}.\end{array}$$

In these expressions,

$$\begin{array}{}\text{(C5)}& {\stackrel{\mathrm{\u203e}}{\mathrm{\Delta}p}}_{i,j,k}^{\mathit{\phi}}:={\displaystyle \frac{\mathrm{\Delta}{p}_{i,j,k}\mathrm{cos}{c}_{j}+\mathrm{\Delta}{p}_{i,j-\mathrm{1},k}\mathrm{cos}{c}_{j-\mathrm{1}}}{\mathrm{cos}{c}_{j}+\mathrm{cos}{c}_{j-\mathrm{1}}}}\phantom{\rule{0.125em}{0ex}}.\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\phantom{\rule{0.125em}{0ex}}\mathit{\delta}{\mathit{\varpi}}_{k}\phantom{\rule{0.125em}{0ex}}\mathrm{cos}{e}_{j}\phantom{\rule{0.125em}{0ex}}.\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}}}\phantom{\rule{0.125em}{0ex}},\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}}}\phantom{\rule{0.125em}{0ex}};\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}}}\phantom{\rule{0.125em}{0ex}}.\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).

Code and data availability

Back to toptop
Code and data availability.

The code used in the numerical simulations of this paper is available
under: https://doi.org/10.5281/zenodo.3529202 (Goldhaber et al., 2019).
CAM6 is published
in the open-access CESM ESCOMP git repository, which is freely available under:
https://github.com/ESCOMP (Goldhaber and Craig, 2020). The AM options can be switched on by
setting standard CAM namelist parameters to non-default values (i.e.
*T* instead of *F*; there are no free numerical parameters). Apart from
these switches, all atmosphere model configurations presented in this
paper are standard CESM cases that can be set up and run using the
scripts provided in the repository. Users can obtain technical support
if requested.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-685-2020-supplement.

Author contributions

Back to toptop
Author contributions.

TT conceived the idea, proposed the work, made the calculations, implemented the code, ran the simulations, evaluated them, produced all figures, and wrote the paper. MB supported this activity through national infrastructure projects of the Norwegian Research Council. CC, BEE, and JE revised the code and included it in the official ESCOMP CESM repository. SG gave technical advice on CAM code and simulations. PHL, MB, and CJ were at hand for critical discussions of the scientific ideas and helped provide the initial impetus for this work. MB, JE, SG, and PHL also provided useful comments and suggestions on the draft paper.

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements

Back to toptop
Acknowledgements.

Warm thanks go to Christoph Heinze for his unbending dedication to model development in the NorESM consortium and in particular for allowing this work to go forward despite the lack of dedicated funding. We are grateful to Alok Gupta at NORCE and Cecile Hannay at NCAR for their assistance with NorESM and CESM development simulations.

Financial support

Back to toptop
Financial support.

This work was partially supported by the Norwegian Research Council (EVA (grant no. 229771) and INES (grant no. 270061)) as well as by the US National Science Foundation Cooperative Agreement (grant no. 1852977).

Review statement

Back to toptop
Review statement.

This paper was edited by Paul Ullrich and reviewed by two anonymous referees.

References

Back to toptop
Arakawa, A. and Lamb, V. R.: A potential enstrophy and energy conserving scheme for the shallow water equations, Mon. Weather Rev., 109, 18–36, 1981.

Blackburn, M., Williamson, D. L., Nakajima, K., Ohfuchi, W., Takahashi, Y. O., Hayashi, Y.-Y., Nakamura, H., Ishiwatari, M., McGregor, J. L., Borth, H., Wirth, V., Frank, H., Bechtold, P., Wedi, N. P., Tomita, H., Satoh, M., Zhao, M., Held, I. M., Suarez, M. J., Lee, M.-I., Watanabe, M., Kimoto, M., Liu, Y., Wang, Z., Molod, A., Rajendran, K., Kitoh, A., and Stratton, R.: The Aqua-Planet Experiment (APE): CONTROL SST simulation, J. Meteorol. Soc. Jpn., 91A, 17–56, https://doi.org/10.2151/jmsj.2013-A02, 2013.

Colella, P. and Woodward, P. R.: The piecewise parabolic method (PPM) for gas-dynamical simulations, J. Comput. Phys., 54, 174–201, 1984.

Courant, R., Friedrichs, K., and Lewy, H.: Über die partiellen Differenzengleichungen der mathematischen Physik, Math. Ann., 100, 3274, https://doi.org/10.1007/BF01448839, 1928.

Egger, J. and Hoinka, K.-P.: The annual cycle of the axial angular momentum of the atmosphere, J. Climate, 18, 757–771, 2005.

Einstein, A.: Die Ursache der Mäanderbildung der Flußläufe und des sogenannten Baerschen Gesetzes, Die Naturwissenschaften, 11, 223–224, 1926.

Feldl, N. and Bordoni, S.: Characterizing the Hadley circulation response through regional climate feedbacks, J. Climate, 29, 613–622, https://doi.org/10.1175/JCLI-D-15-0424.1, 2016.

Gates, W. L.: AMIP: The Atmospheric Model Intercomparison Project, B. Am. Meteorol. Soc., 73, 1962–1970, 1992.

Goldhaber, S. and Craig, C.: CAM: The Community Atmosphere Model, available at: https://github.com/ESCOMP/CAM, last access: 17 February 2020.

Goldhaber, S., Craig, C., and Toniazzo, T.: tto061/CAM: v6.1.029.GMD_AM_paper (Version cam6_1_029), Zenodo, https://doi.org/10.5281/zenodo.3529202, 2019.

Hadley, G.: Concerning the cause of the general trade-winds, Philos. T. R. Soc., 39, 58–62, 1735.

Held, I. M. and Hou, A. Y.: Nonlinear axially symmetric circulations in a nearly inviscid atmosphere, J. Atmos. Sci., 37, 515–533, 1980.

Held, I. M. and Suarez, M. J.: A proposal for the intercomparison of the dynamical cores of atmospheric general circulation models, B. Am. Meteorol. Soc., 75, 1825–1830, 1994.

Hollingsworth, A,, Kållberg, P., Renner, V., and Burridge, D. M.: An internal symmetric computational instability, Q. J. Roy. Meteor. Soc., 109, 417–428, 1983.

IPCC: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, UK and New York, NY, USA, 1535 pp., 2013.

Jablonowski, C. and Williamson, D. L.: A Baroclinic Instability Test Case for Atmospheric Model Dynamical Cores, Q. J. Roy. Meteor. Soc., 132, 2943–2975, 2006.

Laprise, R. and Girard, C.: A spectral general circulation model using a piecewise-constant finite-element representation on a hybrid vertical coordinate system, J. Climate, 3, 32–52, 1990.

Lauritzen, P. H., Bacmeister, J. T., Dubos, T., Lebonnois, S., and Taylor, M. A.: Held-Suarez simulations with the community atmosphere model spectral element (CAM-SE) dynamical core: a global axial angular momentum analysis using Eulerian and floating Lagrangean vertical coordinates, J. Adv. Model. Earth Sy., 6, 129–140, https://doi.org/10.1002/2013MS000268, 2014.

Lebonnois, S., Covey, C., Grossman, A., Parish, H., Schubert, G., Walterscheid, R., Lauritzen, P. H., and Jablonowski, C.: Angular momentum budget in general circulation models of superrotating atmospheres: A critical diagnostic, J. Geophys. Res., 117, E12004, https://doi.org/10.1029/2012JE004223, 2012.

Lin, S. J.: A finite-volume integration method for computing pressure gradient force in general vertical coordinates, Q. J. Roy. Meteor. Soc., 123, 1749–1762, 1997.

Lin, S.-J.: A “vertically lagrangian” finite-volume dynamical core for global models, Mon. Weather Rev., 132, 2293–2307, 2004.

Lin, S. J. and Rood, R. B.: Multidimensional ﬂux-form semi-Lagrangian transport schemes, Mon. Weather Rev., 124, 2046–2070, 1996.

Lin, S. J. and Rood, R. B.: An explicit flux-form semi-Lagrangian shallow-water model on the sphere, Q. J. Roy. Meteor. Soc., 123, 2477–2498, 1997.

Lindzen, R. S. and Hou, A. Y.: Hadley circulations for zonally averaged heating centered off the equator, J. Atmos. Sci., 45, 2416–2427, 1988.

Lipat, B. R., Tselioudis, G., Grise, K. M., and Polvani, L. M.: CMIP5 models' shortwave cloud radiative response and climate sensitivity linked to the climatological Hadley cell extent, Geophys. Res. Lett., 44, 5739–5748, https://doi.org/10.1002/2017GL073151, 2017.

Neale, R. B. and Hoskins, B. J.: A standard test for AGCMs including their physical parameterizations. II: Results for the Met Office model, Atmos. Sci. Lett., 1, 108–114, https://doi.org/10.1006/asle.2000.0024, 2000.

Pauluis, O.: Boundary layer dynamics and cross-equatorial Hadley circulation, J. Atmos. Sci., 61, 1161–1173, 2004.

Quilfen, Y., Chapron, B., Bentamy, A., Gourrion, J., Elfouhaily, T., and Vandemark, D.: Global ERS-1/2 and NSCAT observations: Upwind/crosswind and upwind/downwind measurements, J. Geophys. Res., 104, 11459–11469, 1999.

Schneider, E. K.: Axially symmetric steady-state models of the basic state for instability and climate studies. Part II: Nonlinear calculations, J. Atmos. Sci., 34, 280–297, 1977.

Simmons, A. J. and Burridge, D. M.: An energy and angular-momentum conserving vertical finite-difference scheme and hybrid vertical coordinates, Mon. Weather Rev., 109, 758–766, 1981.

Taylor, K. E.: Summarizing multiple aspects of model performance in a single diagram, J. Geophys. Res., 106, 7183–7192, 2001.

Uppala, S. M., Kållberg, P. W., Simmons, A. J., Andrae, U., Bechtold, V. D. C., Fiorino, M., Gibson, J. K., Haseler, J., Hernandez, A., Kelly, G. A., Li, X., Onogi, K., Saarinen, S., Sokka, N., Allan, R. P., Andersson, E., Arpe, K., Balmaseda, M. A., Beljaars, A. C. M., Berg, L. V. D., Bidlot, J., Bormann, N., Caires, S., Chevallier, F., Dethof, A., Dragosavac, M., Fisher, M., Fuentes, M., Hagemann, S., Hólm, E., Hoskins, B. J., Isaksen, L., Janssen, P. A. E. M., Jenne, R., Mcnally, A. P., Mahfouf, J. F., Morcrette, J. J., Rayner, N. A., Saunders, R. W., Simon, P., Sterl, A., Trenberth, K. E., Untch, A., Vasiljevic, D., Viterbo, P., and Woollen, J.: The ERA-40 Re-Analysis, Q. J. Roy. Meteor. Soc., 131, 2961–3012, https://doi.org/10.1256/qj.04.176, 2005.

Walker, C. C. and Schneider, T.: Eddy influences on Hadley circulations: simulations with an idealized GCM, J. Atmos. Sci., 63, 3333–3350, 2006.

White, A. A., Hoskins, B. J., Roulstone, I., and Staniforth, A.: Consistent approximate models of the global atmosphere: shallow, deep, hydrostatic, quasi-hydrostatic, and non-hydrostatic, Q. J. Roy. Meteor. Soc., 131, 2081–2107, 2005.

Williamson, D. L., Olson, J. G., Hannay, C., Toniazzo, T., Taylor, M., and Yudin, V.: Energy considerations in the Community Atmosphere Model (CAM), J. Adv. Model. Earth Sy., 7, 1178–1188, https://doi.org/10.1002/2015MS000448, 2015.

More precisely, we used a pre-release of CESM2.1.1 (no. 20, 22 March 2019). In terms of the simulations presented in this paper, the differences with the full 2.1.1 release only affect the F2000 cases at f19 resolution, for which slightly different emission datasets are used to force the simulations. The impacts of this are of negligible consequence for the results discussed in this section.

Short summary

We show that ensuring global conservation of the angular (rotational) momentum (AM) of the atmosphere along the Earth's axis of rotation, which is a property of the governing equations, has important and beneficial consequences for the quality of the numerical simulation of the general circulation of the atmosphere. We discuss the causes of non-conservation in the FV dynamical core of the Community Atmosphere Model (CAM), propose remedies, and show their impact in correcting systematic biases.

We show that ensuring global conservation of the angular (rotational) momentum (AM) of the...

Geoscientific Model Development

An interactive open-access journal of the European Geosciences Union