# DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models

Paul A. Ullrich, Christiane Jablonowski, James Kent, Peter H. Lauritzen, Ramachandran Nair, Kevin A. Reed, Colin M. Zarzycki, David M. Hall, Don Dazlich, Ross Heikes, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Paul A. Ullrich, Christiane Jablonowski, James Kent, Peter H. Lauritzen, Ramachandran Nair, et al. "DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models." Geoscientific Model Development, 2017, 10, pp.4477-4509. 10.5194/gmd-10-4477-2017. insu-03727037

HAL Id: insu-03727037

https://insu.hal.science/insu-03727037

Submitted on 27 Jul 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models

Paul A. Ullrich¹, Christiane Jablonowski², James Kent³, Peter H. Lauritzen⁴, Ramachandran Nair⁴, Kevin A. Reed⁵,
Colin M. Zarzycki⁴, David M. Hall⁶, Don Dazlich⁷, Ross Heikes⁷, Celal Konor⁷, David Randall⁷, Thomas Dubos⁸,
Yann Meurdesoif⁸, Xi Chen⁹, Lucas Harris⁹, Christian Kühnlein¹⁰, Vivian Lee¹¹, Abdessamad Qaddouri¹¹,
Claude Girard¹¹, Marco Giorgetta¹², Daniel Reinert¹³, Joseph Klemp⁴, Sang-Hun Park¹⁴, William Skamarock⁴,
Hiroaki Miura¹⁵, Tomoki Ohno¹⁹, Ryuji Yoshida¹⁶, Robert Walko¹⁷, Alex Reinecke¹⁸, and Kevin Viner¹⁸

¹University of California, Davis, Davis, CA, USA

²University of Michigan, Ann Arbor, MI, USA

³University of South Wales, Pontypridd, Wales, UK

⁴National Center for Atmospheric Research, Boulder, CO, USA

⁵Stony Brook University, Stony Brook, NY, USA

⁶University of Colorado, Boulder, Boulder, CO, USA

⁷Colorado State University, Fort Collins, CO, USA

⁸Laboratoire de Météorologie Dynamique, Institut Pierre-Simon Laplace (IPSL), Paris, France

⁹Geophysical Fluid Dynamics Laboratory (GFDL), Princeton, NJ, USA

¹⁰European Center for Medium-Range Weather Forecasting (ECMWF), Reading, UK

¹¹Environment and Climate Change Canada (ECCC), Dorval, Québec, Canada

¹²Max Planck Institute for Meteorology, Hamburg, Germany

¹³Deutscher Wetterdienst (DWD), Offenbach am Main, Germany

¹⁴Yonsei University, Seoul, South Korea

¹⁵University of Tokyo, Bunkyo, Tokyo, Japan

¹⁶RIKEN AICS/Kobe University, Kobe, Japan

¹⁷University of Miami, Coral Gables, FL, USA

¹⁸Naval Research Laboratory, Monterey, CA, USA

¹⁹Japan Agency for Marine-Earth Science and Technology, Yokohama, Kanagawa, Japan

*Correspondence to:* Paul A. Ullrich (paullrich@ucdavis.edu)

Received: 6 May 2017 – Discussion started: 6 June 2017

Revised: 10 October 2017 – Accepted: 18 October 2017 – Published: 6 December 2017

**Abstract.** Atmospheric dynamical cores are a fundamental component of global atmospheric modeling systems and are responsible for capturing the dynamical behavior of the Earth's atmosphere via numerical integration of the Navier-Stokes equations. These systems have existed in one form or another for over half of a century, with the earliest discretizations having now evolved into a complex ecosystem of algorithms and computational strategies. In essence, no two dynamical cores are alike, and their individual successes suggest that no perfect model exists. To better understand modern dynamical cores, this paper aims to provide a com-

prehensive review of 11 non-hydrostatic dynamical cores, drawn from modeling centers and groups that participated in the 2016 Dynamical Core Model Intercomparison Project (DCMIP) workshop and summer school. This review includes a choice of model grid, variable placement, vertical coordinate, prognostic equations, temporal discretization, and the diffusion, stabilization, filters, and fixers employed by each system.
# 1 Introduction

The Dynamical Core Model Intercomparison Project (DCMIP) is an ongoing effort targeting the intercomparison of a fundamental component of global atmospheric modeling systems: the dynamical core. Although this component's role is simply to solve the equations of fluid motion governing atmospheric dynamics (the Navier-Stokes equations), there are numerous confounding factors and compromises that arise from making global simulations computationally feasible. These factors include the choice of model grid, variable placement, vertical coordinate, prognostic equations, representation of topography, numerical method, temporal discretization, physics-dynamics coupling frequency, and the manner in which artificial diffusion, stabilization, filters, and/or energy/mass fixers are applied.

To advance the intercomparison project and provide a unique educational opportunity for students, DCMIP hosted a multidisciplinary 2-week summer school and model intercomparison project at the National Center for Atmospheric Research (NCAR) in June 2016, that invited graduate students, postdocs, atmospheric modelers, expert lecturers, and computer specialists to create a stimulating, unique, and hands-on driven learning environment. The 2016 workshop and summer school followed from earlier DCMIP and dynamical core workshops (held in 2012 and 2008, respectively), and other model intercomparison efforts. Its goals were to provide an international forum for discussing outstanding issues in global atmospheric models and provide a unique training experience for the future generation of climate scientists. Special attention was paid to the role of simplified physical parameterizations, physics-dynamics coupling, non-hydrostatic atmospheric modeling, and variable-resolution global modeling. The summer school and model intercomparison project promoted active learning, innovation, discovery, mentorship, and the integration of science and education. Modeling groups were then invited to contribute model descriptions and results to the intercomparison effort for publication.

The summer school directly benefited its participants by providing a unique educational experience and an opportunity to interact with modeling teams from around the world. The workshop is expected to have further repercussions on the development of operational atmospheric modeling systems by allowing modeling groups to assess their models in the context of the global dynamical core ecosystem. Past and present intercomparison efforts have been leveraged by modeling groups to improve their own models, in turn leading to a positive impact on the quality of weather and climate simulations. The workshop component of DCMIP has also advanced our knowledge of (1) the relative behaviors exhibited by atmospheric dynamical cores, (2) differences that arise among mechanisms for coupling the physical parameterizations and dynamical core, and (3) the impacts of variable-resolution refinement regions and transition zones

in global atmospheric simulations. Notably, the use of idealized test cases to isolate specific phenomena gave us a unique opportunity to assess specific differences that arise due to the choice of dynamical core. Another important outcome of the workshop was the development of a standard test case suite and benchmark set of simulations that can be used for assessment of any future dynamical core. The test cases introduced in the 2016 workshop build on the previous DCMIP test case suites (Jablonowski et al., 2008; Ullrich et al., 2012) with tests that now incorporate simplified moist physics.

This paper is the first in a series of papers documenting the results of this workshop. Its purpose is two-fold: first, to review the multitude of technologies and techniques that have been developed for non-hydrostatic global atmospheric modeling; and second, to provide a mechanism to understand the differences that arise in the test cases of later papers in this series. For ease of reference, a list of mathematical symbols that are employed in this paper (and subsequent DCMIP papers) is given in Table 1. Section 2 then provides a brief overview of each of the participating models, along with a tabulation of relevant details about the dynamical core design. The body of this paper is dedicated to an overview of techniques available for building the infrastructure of a global dynamical core: Sect. 3 describes aspects of the horizontal discretization, including model grids and horizontal placement of prognostic variables; Sect. 4 describes the vertical placement of model variables and choice of vertical coordinates; Sect. 5 describes aspects of variable placement and prognosis; Sect. 6 describes diffusion, stabilization, filters, and fixers employed by these models; and Sect. 7 describes temporal discretizations. The summary and conclusions then follow in Sect. 8. Finally, Appendix A provides a comprehensive overview of the various forms the Navier-Stokes equations take in dynamical cores, and has been included as a resource for dynamical core developers.

# 2 Dynamical cores

This section provides a brief overview of key discretization choices, along with unique features or design specifications from participating dynamical cores. Further details on these choices can be found in subsequent sections. In total, simulation results and model descriptions have been submitted from 11 dynamical cores (see Table 2). The prognostic variables employed and horizontal discretizations for these dynamical cores are summarized in Table 3. The vertical staging of variables and vertical coordinate choice are summarized in Table 4. Principal options for diffusion, stabilization, filters, or fixers along with the temporal discretization for these models are summarized in Table 5. A brief description of each participant model follows, focused on the unique features and decisions underlying the model design.
<table><caption>**Table 1.** A standard list of symbols used throughout this paper and in the DCMIP.</caption><thead><tr><th>Symbol</th><th>Description</th></tr></thead><tbody><tr><td>λ</td><td>Longitude (in radians)</td></tr><tr><td>φ</td><td>Latitude (in radians)</td></tr><tr><td>z</td><td>Height with respect to mean sea level (set to zero)</td></tr><tr><td>s</td><td>Vertical model coordinate</td></tr><tr><td>p<sub>s</sub></td><td>Surface pressure (p<sub>s</sub> of moist air if q > 0)</td></tr><tr><td>Φ</td><td>Geopotential</td></tr><tr><td>Φ<sub>s</sub></td><td>Surface geopotential</td></tr><tr><td>z<sub>s</sub></td><td>Surface elevation with respect to mean sea level (set to zero)</td></tr><tr><td>u</td><td>Zonal wind velocity</td></tr><tr><td>v</td><td>Meridional wind velocity</td></tr><tr><td>w</td><td>Vertical wind velocity</td></tr><tr><td>ζ</td><td>GEM vertical coordinate velocity</td></tr><tr><td>u</td><td>3-D wind vector</td></tr><tr><td>u<sub>h</sub></td><td>Horizontal wind vector</td></tr><tr><td>v<sub>h</sub></td><td>Horizontal wind vector with covariant components</td></tr><tr><td>ω</td><td>Vertical pressure velocity</td></tr><tr><td>D</td><td>Divergence of the horizontal wind vector</td></tr><tr><td>ζ</td><td>Vertical component of relative vorticity</td></tr><tr><td>p</td><td>Pressure (pressure of moist air if q > 0)</td></tr><tr><td>e</td><td>Internal energy</td></tr><tr><td>ρ</td><td>Total air density</td></tr><tr><td>ρ<sub>d</sub></td><td>Dry air density</td></tr><tr><td>ρ<sub>s</sub></td><td>Pseudo-density</td></tr><tr><td>T</td><td>Temperature</td></tr><tr><td>T<sub>v</sub></td><td>Virtual temperature</td></tr><tr><td>θ</td><td>Potential temperature</td></tr><tr><td>θ<sub>v</sub></td><td>Virtual potential temperature</td></tr><tr><td>θ<sub>il</sub></td><td>Ice-liquid potential temperature</td></tr><tr><td>θ<sub>ρ</sub></td><td>Density potential temperature</td></tr><tr><td>q</td><td>Specific humidity</td></tr><tr><td>q<sub>v</sub></td><td>Water vapor mixing ratio</td></tr><tr><td>q<sub>c</sub></td><td>Cloud water mixing ratio</td></tr><tr><td>q<sub>r</sub></td><td>Rain water mixing ratio</td></tr><tr><td>q<sub>i</sub></td><td>General tracer mixing ratio</td></tr></tbody></table>

## 2.1 Accelerated Climate Model for Energy-Atmosphere (ACME-A)

The Accelerated Climate Model for Energy-Atmosphere (ACME-A) has much in common with the Community Atmosphere Spectral Element Model (CAM-SE) (Dennis et al., 2012) as both share a common origin in the High Order Method Modeling Environment (HOMME) (Taylor and Fournier, 2010). ACME-A employs both a hydrostatic model and an experimental non-hydrostatic compressible shallow-atmosphere model. Both variants are designed to be mass and energy conserving, with nearly optimal parallel scalability at large core counts. ACME-A is built upon an unstructured grid of quadrilateral elements arranged in a cubed-sphere configuration (Sect. 3.2), although unstructured, re-

gionally refined meshes with conforming edges may also be employed. The fluid equations are discretized using dimensional splitting, with a nodal fourth-order spectral element discretization in the horizontal and vertical floating Lagrangian levels in hybrid terrain-following pressure coordinates (Sect. 4.2.3). Vertical operators are based on the mimetic (mass- and energy-conserving) second-order finite difference discretization of Simmons and Burridge (1981). All fields are co-located in the horizontal, in the sense that they share the same fourth-order basis functions. Tracer transport is subcycled relative to the hydrodynamics, using the spectral element method, with tracer mass as the prognostic variable.

## 2.2 Colorado State University (CSU) model

The Colorado State University (CSU) model is a finite-volume model using an optimized geodesic grid (Heikes and Randall, 1995; Heikes et al., 2013) (Sect. 3.4), with height as the vertical coordinate. The model is based on the non-hydrostatic unified system of equations proposed by Arakawa and Konor (2009), which filters vertically propagating sound waves but allows the Lamb wave and does not require a reference state. The horizontal wind field is determined by predicting the vertical component of the vorticity and the divergence of the horizontal wind, and then solving a pair of two-dimensional Poisson equations for a stream function and velocity potential. Horizontal diffusion is included in the form of a fourth-order hyperviscosity operator applied on constant height surfaces ($\nabla_z^4$) that acts on the vorticity, divergence, potential temperature, and tracer (Sect. 6.2). The CSU model supports both third-order and fifth-order upstream-weighted, finite-volume advection schemes, with positivity preservation enforced via mass borrowing.

## 2.3 DYNAMICO

DYNAMICO is a mimetic finite-difference/finite-volume model using a geodesic grid (Sect. 3.4) and a floating vertical mass coordinate (Sect. 4.2.3). Although originally a hydrostatic model, it has been recently extended to solve the shallow-atmosphere non-hydrostatic Euler equations. DYNAMICO's design uniquely combines a representation of the prognostic and diagnostic fields following the ideas of discrete differential geometry (Dubos et al., 2015). It includes a novel Hamiltonian formulation of the equations of motion in non-Eulerian coordinates (Dubos and Tort, 2014) which is imitated at the discrete level using building blocks from the literature (Thuburn et al., 2009; Ringler et al., 2010) and (up to the addition of explicit diffusion) leads to an energy-conserving spatial discretization. It also incorporates a novel explicit-implicit splitting which results in a simple, efficient, and scalable implicit solver while allowing stable time steps close or identical to those of the hydrostatic solver). Horizontal diffusion is included via a fourth-order hyperviscos-
<table><caption>**Table 2.** Participating modeling centers and associated dynamical cores that have submitted a model description and/or simulation results.</caption><thead><tr><th>Short name</th><th>Long name</th><th>Modeling center or group</th></tr></thead><tbody><tr><td>ACME-A</td><td>Atmosphere model of the Accelerated Climate Model for Energy</td><td>Sandia National Laboratories and University of Colorado, Boulder, USA</td></tr><tr><td>CSU</td><td>Colorado State University Model</td><td>Colorado State University, USA</td></tr><tr><td>DYNAMICO</td><td>DYNAMical core on the ICOsahedron</td><td>Institut Pierre Simon Laplace (IPSL), France</td></tr><tr><td>FV<sup>3</sup></td><td>GFDL Finite-Volume Cubed-Sphere Dynamical Core</td><td>Geophysical Fluid Dynamics Laboratory, USA</td></tr><tr><td>FVM</td><td>Finite Volume Module of the Integrated Forecasting System</td><td>European Centre for Medium-Range Weather Forecasts</td></tr><tr><td>GEM</td><td>Global Environmental Multiscale model</td><td>Environment and Climate Change Canada, Canada</td></tr><tr><td>ICON</td><td>ICOsahedral Non-hydrostatic model</td><td>Max-Planck-Institut für Meteorologie, Germany</td></tr><tr><td>MPAS</td><td>Model for Prediction Across Scales</td><td>National Center for Atmospheric Research, USA</td></tr><tr><td>NICAM</td><td>Non-hydrostatic Icosahedral Atmospheric Model</td><td>AORI/JAMSTEC/AICS, Japan</td></tr><tr><td>OLAM</td><td>Ocean Land Atmosphere Model</td><td>Duke University/University of Miami, USA</td></tr><tr><td>Tempest</td><td>Tempest Non-hydrostatic Atmospheric Model</td><td>University of California, Davis, USA</td></tr></tbody></table>

<table><caption>**Table 3.** Details on the prognostic variables and horizontal discretization for participating dynamical cores. The equation set indicates whether a model is hydrostatic (H) or non-hydrostatic (NH), and whether the model presently supports the deep-atmosphere formulation (D). Only three numerical methods are represented among participating models, namely finite difference (FD), finite volume (FV), and spectral element (SE). More details on horizontal staggering can be found in Sect. 3.8.</caption><thead><tr><th>Short name</th><th>Equation set</th><th>Prognostic variables</th><th>Horizontal grid</th><th>Numerical method</th><th>Horizontal staggering</th></tr></thead><tbody><tr><td>ACME-A</td><td>H/NH</td><td>u<sub>h</sub>, w, ρ<sub>s</sub>, ρ<sub>s</sub>θ, Φ, ρ<sub>s</sub>q<sub>i</sub></td><td>Cubed sphere (Sect. 3.2)</td><td>SE</td><td>A grid</td></tr><tr><td>CSU</td><td>NH (unified)</td><td>ζ, D, w, ρ<sub>s</sub>, θ<sub>v</sub>, q<sub>i</sub></td><td>Geodesic (Sect. 3.4)</td><td>FV</td><td>Z grid</td></tr><tr><td>DYNAMICO</td><td>H/NH</td><td>v<sub>h</sub>, ρ<sub>s</sub>w, ρ<sub>s</sub>, ρ<sub>s</sub>θ<sub>v</sub>, Φ, ρ<sub>s</sub>q<sub>i</sub></td><td>Geodesic (Sect. 3.4)</td><td>FV</td><td>C grid</td></tr><tr><td>FV<sup>3</sup></td><td>NH</td><td>u<sub>h</sub>, w, ρ<sub>s</sub>, ρ<sub>s</sub>θ<sub>v</sub>, Φ, ρ<sub>s</sub>q<sub>i</sub></td><td>Cubed sphere (Sect. 3.2)</td><td>FV</td><td>D grid</td></tr><tr><td>FVM</td><td>NH (D)</td><td>ρ<sub>d</sub>, u<sub>h</sub>, w, θ', q<sub>i</sub></td><td>Octahedral (Sect. 3.6)</td><td>FV</td><td>A grid</td></tr><tr><td>GEM</td><td>NH</td><td>u<sub>h</sub>, w, ζ, T<sub>v</sub>, p, q<sub>i</sub></td><td>Yin-Yang (Sect. 3.7)</td><td>FD</td><td>C grid</td></tr><tr><td>ICON</td><td>NH (D)</td><td>u<sub>h</sub>, w, ρ, θ<sub>v</sub>, ρq<sub>i</sub></td><td>Icosahedral triangular (Sect. 3.3)</td><td>FV</td><td>C grid</td></tr><tr><td>MPAS</td><td>NH</td><td>ρ<sub>d</sub>u<sub>h</sub>, ρ<sub>d</sub>w, ρ<sub>d</sub>, ρ<sub>d</sub>θ<sub>v</sub>, ρ<sub>d</sub>q<sub>i</sub></td><td>CCVT (Sect. 3.5)</td><td>FV</td><td>C grid</td></tr><tr><td>NICAM</td><td>NH</td><td>ρu<sub>h</sub>, ρw, ρ, ρe, ρq<sub>i</sub></td><td>Geodesic (Sect. 3.4)</td><td>FV</td><td>A grid</td></tr><tr><td>OLAM</td><td>NH (D)</td><td>ρu<sub>h</sub>, ρw, ρ, ρθ<sub>il</sub>, ρq<sub>i</sub></td><td>Geodesic (Sect. 3.4)</td><td>FV</td><td>C grid</td></tr><tr><td>Tempest</td><td>NH</td><td>u<sub>h</sub>, w, ρ, ρθ<sub>v</sub>, ρq<sub>i</sub></td><td>Cubed sphere (Sect. 3.2)</td><td>SE</td><td>A grid</td></tr></tbody></table>

ity operator (Sect. 6.3). In addition, it features a conservative positive-definite transport scheme based on a slope-limited finite-volume approach (Dubey et al., 2015).

## 2.4 FV cubed (FV³)

The GFDL Finite-Volume Cubed-Sphere Dynamical Core (FV³, or sometimes written FV3) is a finite-volume model that solves the non-hydrostatic Euler equations on the equiangular gnomonic cubed-sphere grid (Sect. 3.2) with a floating Lagrangian vertical coordinate. The Lagrangian vertical coordinate deforms so that the flow is constrained to follow the Lagrangian surfaces, allowing vertical transport to be represented implicitly without additional advection terms (see Sect. 4.2.3 below). The non-hydrostatic formulation extends the hydrostatic model described in Lin (2004) by adding a prognostic vertical velocity and geometric height of each grid cell, which can then be used to compute density. The discretization is on the C–D grid as described by Lin

and Rood (1997) (see also Sect. 3.8), although the prognostic horizontal winds are stored in the native gnomonic local coordinate. All variables are 3-D cell-mean values, except for the horizontal winds, which are 2-D face-mean values on their respective staggerings; as a result, diagnostic vorticity is a 3-D cell-mean value. Fluxes are computed using the piece-wise parabolic method of Colella and Woodward (1984) with an optional monotonicity constraint; in non-hydrostatic applications, the monotonicity constraint is used primarily for tracer transport. Since divergence is effectively invisible to the solver, a 2-D divergence damping is applied to control numerical noise as divergent modes cascade to the grid scale (Sect. 6.4). Implicit viscosity is applied through the monotonicity constraint; if non-monotonic advection is used for the momentum and total air mass, a weak explicit hyperviscosity is applied for stability and to alleviate numerical noise. Explicit viscosity is applied every acoustic time step.
<table><caption>**Table 4.** Vertical staggering (detailed in Sect. 4.1) and vertical coordinates (detailed in Sect. 4.2) for participating dynamical cores.</caption><thead><tr><th>Acronym</th><th>Vertical staggering</th><th>Vertical coordinate</th></tr></thead><tbody><tr><td>ACME-A</td><td>Co-located</td><td>Floating mass (Sect. 4.2.3)</td></tr><tr><td>CSU</td><td>Lorenz</td><td>Fixed height</td></tr><tr><td>DYNAMICO</td><td>Lorenz</td><td>Floating mass (Sect. 4.2.3)</td></tr><tr><td>FV<sup>3</sup></td><td>Co-located</td><td>Floating mass (Sect. 4.2.3)</td></tr><tr><td>FVM</td><td>Co-located</td><td>Fixed height</td></tr><tr><td>GEM</td><td>Modified Charney-Phillips (Sect. 4.1)</td><td>Log pressure (Sect. 4.2.2)</td></tr><tr><td>ICON</td><td>Lorenz</td><td>Fixed height</td></tr><tr><td>MPAS</td><td>Lorenz</td><td>Fixed height</td></tr><tr><td>NICAM</td><td>Lorenz</td><td>Fixed height</td></tr><tr><td>OLAM</td><td>Lorenz</td><td>Fixed height with cut cells (Sect. 4.2.4)</td></tr><tr><td>Tempest</td><td>Lorenz</td><td>Fixed height</td></tr></tbody></table>

<table><caption>**Table 5.** Principal options for diffusion, stabilization, filters, or fixers in participating dynamical cores (detailed in Sect. 6) and temporal discretization (detailed in Sect. 7).</caption><thead><tr><th>Acronym</th><th>Principal options for diffusion, stabilization, filters, or fixers</th><th>Temporal discretization</th></tr></thead><tbody><tr><td>ACME-A</td><td>Fourth-order horizontal hyperviscosity</td><td>KGU53 (Guerra and Ullrich, 2016)</td></tr><tr><td>CSU</td><td>Fourth-order horizontal hyperviscosity</td><td>Third-order Adams-Bashforth (AB3)</td></tr><tr><td>DYNAMICO</td><td>Fourth-order horizontal hyperviscosity</td><td>ARK(2,3,2) (Giraldo et al., 2013)</td></tr><tr><td>FV<sup>3</sup></td><td>Divergence damping, hyperviscosity</td><td>Forward-backward (Lin and Rood, 1997)/semi-implicit</td></tr><tr><td>FVM</td><td>Monotonic limiting</td><td>Semi-implicit (Smolarkiewicz et al., 2014) (Sect. 7.2)</td></tr><tr><td>GEM</td><td>Hyperviscosity</td><td>Semi-implicit (Girard et al., 2014) (Sect. 7.3)</td></tr><tr><td>ICON</td><td>Divergence damping, Smagorinsky, hyperdiffusion</td><td>Predictor-corrector</td></tr><tr><td>MPAS</td><td>Smagorinsky, hyperdiffusion</td><td>Split-explicit (Klemp et al., 2007)</td></tr><tr><td>NICAM</td><td>3-D divergence damping, Smagorinsky, hyperviscosity</td><td>Split-explicit (Klemp et al., 2007)</td></tr><tr><td>OLAM</td><td>Divergence/vorticity damping</td><td>Second-order Adams-Bashforth, Lax-Wendroff (for tracers)</td></tr><tr><td>Tempest</td><td>Fourth-order horizontal hyperviscosity</td><td>ARS(2,3,2) (Ascher et al., 1997)</td></tr></tbody></table>

## 2.5 Finite-Volume Module (FVM) of the Integrated Forecasting System

The Finite-Volume Module (FVM) of the Integrated Forecasting System (IFS) is currently under development at ECMWF (Smolarkiewicz et al., 2016; Kühnlein and Smolarkiewicz, 2017; Smolarkiewicz et al., 2017). FVM solves the non-hydrostatic Euler equations on an octahedral reduced Gaussian grid (Sect. 3.6) with a height-based terrain- following vertical coordinate (Szmelter and Smolarkiewicz, 2010; Smolarkiewicz et al., 2016). The horizontal spatial discretization uses the median-dual finite-volume approach, combined with a structured-grid finite-difference method in the vertical. In both the horizontal and vertical discretizations, all variables are co-located. A centered two- time-level, semi-implicit integration scheme is employed with 3-D implicit treatment of acoustic, buoyant, and ro- tational modes (Smolarkiewicz et al., 2014) (Sect. 7.2). The associated 3-D Helmholtz problem is solved itera- tively using a bespoke preconditioned generalized conjugate residual approach. The integration procedure uses the non- oscillatory, finite-volume MPDATA (multidimensional posi-

tive definite advection transport algorithm) advection scheme (Smolarkiewicz and Szmelter, 2005; Kühnlein and Smolarkiewicz, 2017). The non-oscillatory (i.e., monotonic) MP- DATA also provides sufficient dissipation/diffusion to stabi- lize the model, so no other explicit filtering mechanism is required (Sect. 6.5). Note that the octahedral reduced Gaus- sian grid is also employed in the spectral-transform dynami- cal core of the presently operational IFS at ECMWF, which facilitates interoperability of the two formulations. However, FVM is not restricted to this grid and offers capabilities to- wards broad classes of meshes (Szmelter and Smolarkiewicz, 2010; Kühnlein et al., 2012; Deconinck et al., 2017).

## 2.6 Global Environmental Multiscale (GEM) model

The Global Environmental Multiscale (GEM) model (Gi- rard et al., 2014) is used for operational forecasting at Environment and Climate Change Canada. GEM solves the non-hydrostatic Euler equations on the Yin-Yang grid (Kageyama and Sato, 2004) (Sect. 3.7) with Arakawa C- grid staggering of prognostic variables. The vertical coordi- nate is a unique hybrid terrain-following coordinate of a log-
hydrostatic-pressure type (Sect. 4.2.2) and the vertical discretization is based on the Charney–Phillips grid (Sect. 4.1). A two-time-level, semi-Lagrangian implicit time discretization is implemented as described in Sect. 7.3. It gives rise to an iterative process where each step requires the solution of a linear system of equations that is reduced to a Helmholtz problem for one composite variable. For this problem, a direct solver is involved, using the Schwarzs-type domain decomposition method (Qaddouri et al., 2008). Semi-Lagrangian advection is also used for tracer transport. To eliminate numerical noise, an explicit hyperviscosity is employed for wind components and tracers via applications of the Laplacian operator, applied after the completion of the physics time step (Sect. 6.6).

## 2.7 ICOsahedral Non-hydrostatic (ICON) model

The ICOsahedral Non-hydrostatic (ICON) model (Zängl et al., 2015) is a finite-volume model that solves the non-hydrostatic Euler equations in 2-D vector-invariant form on an icosahedral (triangular) grid (Sect. 3.3) with Arakawa C-grid staggering, and further utilizing a smoothed terrain-following height-based Lorenz vertical discretization. Prognostic horizontal velocities are stored as normal wind components at the edge midpoints of full levels. Prognostic vertical velocity is stored at the circumcenters of the triangles on half levels. The discretization employs a two-time-level predictor-corrector scheme, which is explicit in all terms except for those describing the vertical propagation of sound waves. For stabilization of the divergence term on the triangular C grid, the divergence in a triangle is computed from modified normal wind components, resulting from a weighted average, including normal winds on edges of adjacent cells. Further divergence damping is applied to the normal wind at every substep. Rayleigh damping is applied to the vertical wind in layers close to the model top in order to avoid the reflection of gravity waves. The horizontal diffusion, which is applied at full model time steps, combines a flow-dependent Smagorinsky scheme with a background fourth-order Laplacian diffusion operator (Sect. 6.7). For tracer transport, a flux-form semi-Lagrangian scheme with monotone flux limiters is used, which leads to local mass conservation and consistency with the air motion. Specifically, the average air mass flux of the dynamical substeps is provided to the tracer transport to allow for mass-consistent transport. These numerical methods have been chosen for high numerical efficiency, and they rely on next-neighbor communication only, thus allowing massive parallelization.

## 2.8 Model for Prediction Across Scales (MPAS)

The Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012) is a finite-volume model that solves the non-hydrostatic Euler equations using an Arakawa C-grid staggering on a centroidal Voronoi tessellation mesh

(Sect. 3.5) and the mimetic TRISK discretization (Thuburn et al., 2009; Ringler et al., 2010). In the vertical, MPAS employs a Lorenz-type second-order nodal finite volume method with a smoothed terrain-following height coordinate. Advection is nominally third- to fourth-order and is handled in accordance with Skamarock and Gassmann (2011). The prognostic variables are dry air pseudo-density ($\tilde{\rho}_d$), dry momentum ($\tilde{\rho}_d u$), and a modified moist potential temperature. Integration in time is handled via the split-explicit method of Klemp et al. (2007). Various filters are available for controlling spurious oscillations, including Smagorinsky-type eddy viscosity, fourth-order hyperdiffusion, and 2-D and 3-D divergence damping operators (Sect. 6.8).

## 2.9 Non-hydrostatic ICOsahedral Atmospheric Model (NICAM)

Non-hydrostatic ICOsahedral Atmospheric Model (NICAM) is a finite-volume model that solves the non-hydrostatic Euler equations using a geodesic grid (Sect. 3.4) optimized with spring dynamics using the method of Tomita et al. (2002). A terrain-following height coordinate system is used in the vertical (Tomita and Satoh, 2004) with Lorenz staggering. Instead of temperature or potential temperature, total energy is prognosed following the method of Satoh (2002, 2003). All prognostic variables are collocated horizontally at the mass centroid of each hexagonal/pentagonal cell to mitigate accuracy reduction under cell averaging, which is required in converting cell-integrated quantities to point values at cell centroids. The use of cell centroids ensures quasi-second-order accuracy of the gradient and divergence operators of NICAM (Tomita et al., 2001). For integration in time, a two-stage Runge-Kutta scheme is usually employed because of low computational cost, although a three-stage Runge-Kutta scheme (Wicker and Skamarock, 2002) is available and recommended. The split-explicit time discretization is used for the horizontally propagating sound waves with the 3-D divergence damping term (Skamarock and Klemp, 1992) (Sect. 6.9). An implicit time discretization is adopted for the vertically propagating wave modes. A variant of the piece-wise linear transport scheme (Miura, 2007; Niwa et al., 2011) is used with a flux limiter of Thuburn (1997) for passive tracer transports.

## 2.10 Ocean-Land-Atmosphere Model (OLAM)

Ocean-Land-Atmosphere Model (OLAM) (Walko and Avissar, 2008a, b, 2011) is a finite-volume model that solves the deep-atmosphere non-hydrostatic Euler equations in momentum conservation form on a hexagonal Voronoi mesh (Sect. 3.4) with Arakawa C-grid staggering. The model supports optional local mesh refinement, which introduces some pentagons and heptagons to the grid. Height is the vertical coordinate, and a Lorenz vertical grid staggering is used. A unique feature of OLAM is that grid levels are horizontal
and intersect topography (Sect. 4.2.4). This avoids a number of well-documented errors associated with terrain-following grids and also eliminates the need for evaluation of coordinate transformation terms. Topography is represented as a smooth (non-stepped) surface by means of cut cells whose surfaces and volume are reduced according to the portion of each cell that is below ground. The OLAM cut-cell formulation conserves mass and momentum. Acoustic modes are solved explicitly in the horizontal, using time splitting and a second-order Lax–Wendroff method, and implicitly in the vertical. Tracer transport is second order in space and time, using the scheme of Miura (2007), with consistent fluxes obtained by time averaging over the acoustic time steps.

## 2.11 Tempest

The Tempest model (Ullrich, 2014a; Guerra and Ullrich, 2016) is an experimental test bed for high-performance numerical methods that solves the non-hydrostatic Euler equations on a cubed-sphere grid (Sect. 3.2) using a horizontally co-located spectral element discretization. In the vertical, Tempest uses an Eulerian finite-volume discretization with Lorenz staggering and terrain-following height coordinates. The implementation includes both fully explicit time integration and a horizontally explicit vertically implicit formulation that is solved with a third-order implicit-explicit additive Runge–Kutta scheme from Ascher et al. (1997). Fourth-order hyperviscosity is used in the horizontal to prevent a buildup of energy at the grid scale (Sect. 6.1). The model further provides an optional upwind-biased transport scheme in the vertical column. Tracer transport is performed using the spectral element method with the same time step as the hydrodynamics and using the tracer mass density as a prognostic variable. As with the hydrodynamics, tracer transport is performed explicitly in the horizontal and implicitly in the vertical.

## 3 Horizontal discretization and model grids

The horizontal discretization determines how the atmosphere, which consists of a set of approximately continuous fields, is mapped into a very limited and discrete computational space. The horizontal discretization essentially consists of two major choices: the model grid, which determines the density and connectivity of discrete regions (Staniforth and Thuburn, 2012), and the arrangement of prognostic and diagnostic variables around each grid region (Arakawa and Lamb, 1977). In order to meet demands for high computational efficiency and equal partitioning of computation across large parallel systems, modern dynamical cores have explored a number of options for model grids. The choice of model grid can be motivated by simplicity, as in the case of the latitude–longitude grid; by a desire to maintain a local Cartesian structure, as with the cubed-sphere grid; or

to support grid isotropy and homogeneity, as with many of the hexagonal or Voronoi grids that have been employed. The choice of grid may be further decided by the numerical method; for instance, finite element models that use tensor products to define basis functions require grids consisting entirely of quadrilaterals. Inevitably, a choice must be made, and the pros and cons of that choice will impact other decisions related to the model. To better understand the options that are available to dynamical core developers, we begin by reviewing many of the model grids that have been employed in global dynamical cores around the world. Then, in Sect. 3.8, we discuss the “staggering” of model variables, referring to the distribution of variables within and around each grid cell.

### 3.1 Latitude–longitude grid

The classic latitude–longitude grid is produced by subdividing the sphere along lines of constant latitude and longitude. The latitude–longitude grid has the benefits of being globally rectilinear, which simplifies data access and subdivision of computation across processors, and yields a vector basis that is locally orthogonal nearly everywhere. This structure accurately maintains purely zonal flows and simplifies data post-processing for visualization. Because of the convergence of grid lines near the poles, the operational use of this grid requires that the associated numerical scheme be resilient to arbitrarily small Courant numbers, or that polar filtering be employed to remove unstable computational modes (Lin, 2004). This grid is presently employed in many global models, including the UK Met Office New Dynamics and ENDGame dynamical cores (Davies et al., 2005; Wood et al., 2014). The latitude–longitude grid is also an option in the GEM model.

### 3.2 Cubed-sphere grid

The equiangular, gnomonic cubed-sphere grid (Sadourny, 1972; Ronchi et al., 1996; Putman and Lin, 2007) consists of six Cartesian patches arranged along the faces of a cube which is then inflated onto a spherical shell. More information on this choice of grid can be found in Ullrich (2014a). On the equiangular cubed-sphere grid, coordinates are given as $(\alpha, \beta, p)$, with central angles $\alpha, \beta \in [-\frac{\pi}{4}, \frac{\pi}{4}]$ and panel index $p$. The structure of this grid supports refinement through stretching (Schmidt, 1977; Harris et al., 2016) or nesting (Harris and Lin, 2013). The Cartesian structure of cubed-sphere grid panels is advantageous for numerical methods that are formulated in Cartesian coordinates or that utilize dimension splitting. Nonetheless, special treatment of the panel boundaries is often necessary since they represent coordinate discontinuities. This grid is depicted in Fig. 1a. Among the DCMIP2016 models, the cubed-sphere grid is employed by the ACME, $FV^3$, and Tempest dynamical cores.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) A cubed-sphere grid. (b) An icosahedral (triangular) grid with additional refinement over Europe, as indicated in red. (c) An icosahedral (hexagonal) grid.</figcaption></figure>

## 3.3 Icosahedral (triangular) grid

The icosahedral triangular grid is derived from the spherical icosahedron that consists of 20 equilateral spherical triangles, 30 great circle edges, and 12 vertices. These initial triangles are then subdivided repeatedly until the desired mean resolution is obtained. For a single subdivision, each edge is divided in *n* arcs of equal length, thus defining new vertices, which by proper connection to other new vertices result in *n*² triangles filling the original triangle. By construction, the new vertices share six triangles; thus, the refinement process brakes the initial isotropy of the icosahedron and results in non-equilateral triangles of different sizes. Among the DCMIP2016 models, the icosahedral (triangular) grid is employed operationally in the ICON dynamical core.

Several methods are available for subdividing the triangular regions. One such approach is implemented by the ICON grid generator, which allows an “arbitrary” subdivision factor *n* for the first refinement step only, the so-called root refinement. Typical choices are *n* = 2, 3, or 5. All additional *m* refinement steps use *n* = 2; i.e., they are bisection steps. A global grid resulting from a root division factor *n* and *m* bisections, denominated as *RnBm* grid, has $n_c = 20 \cdot n^2 \cdot 2^{2m}$ cells, $n_e = 3/2 \cdot n_c$ edges, and $n_v = 10 \cdot n^2 \cdot 2^{2m} + 2$ vertices. The anisotropy of global grids is reduced by the spring dynamics of Tomita et al. (2001). An example of such a grid is depicted in Fig. 1b. A discussion of the effective resolution of such grids is given in Dipankar et al. (2015). The ICON grid generator further allows for inset regional grids, produced by additional refinement steps that are only applied over a limited region or set of regions. The dynamical core then allows for either one-way or two-way coupling of the refined region to the parent model. The current operational numerical weather prediction of the Deutscher Wetterdienst (German Weather Service, DWD), for instance, uses a *R3B7* global grid with 2949120 cells and 13 km mean resolution in combination with a refined region over Europe at 6.5 km resolution.

## 3.4 Icosahedral (hexagonal) grid/geodesic grid

The icosahedral (hexagonal) grid, also commonly referred to as the geodesic grid, is most directly obtained by taking the dual to the icosahedral (triangular grid) – that is, by replacing grid nodes with spherical polygons. The resulting grid's cells are hexagonal, except for 12 pentagonal cells. Given an icosahedral–triangular mesh, vertices of the corresponding icosahedral–hexagonal mesh are then defined as either circumcenters or barycenters of triangles, leading to either a Voronoi mesh, used by DYNAMICO (see also Sect. 3.5), or a barycentric mesh, used by NICAM. A Voronoi mesh has the property that triangular edges are perpendicular to edges of hexagons/pentagons, facilitating the formulation of certain finite-difference and finite-volume numerical schemes. The resulting highly homogeneous and isotropic grid then appears analogous to the grid in Fig. 1c. Unlike the cubed-sphere and icosahedral (triangular) grids, grid cells on this geodesic grid are guaranteed to be edge neighbors (cells that share a given edge) if they are also node neighbors (cells that share a given node). Among the DCMIP2016 models, the geodesic grid is employed by the CSU, DYNAMICO, NICAM, and OLAM dynamical cores.

It is often useful to optimize icosahedral–hexagonal grids as well. DYNAMICO applies a number of iterations of Lloyd’s algorithm (Lloyd, 1982), following by replacing the vertices of the original triangular mesh by the centroid of hexagons/pentagons, then regenerating the icosahedral–hexagonal mesh. This improves the homogeneity of the grid (e.g., ratio of largest cell area to smallest cell area), but several thousand iterations can be required for a significant improvement.

OLAM optimizes by applying the spring dynamics method of Tomita et al. (2001) to the dual triangular mesh prior to its mapping to the Voronoi mesh. When local mesh refinement is applied, which OLAM achieves in a series of one or more resolution-doubling steps, each spanning a transition zone that is three grid rows wide (Fig. 2), the equilibrium spring length is scaled to the target grid cell size in each refinement level and is varied incrementally across the transition zone. Spring dynamics is further modified by forcing
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Detail of one step of local mesh refinement used by the OLAM Voronoi mesh. The transition zone is constructed by explicit topological reconnection of the grid lines, which produces pairings of heptagons (red dots) and pentagons (blue dots) along the refinement perimeter.</figcaption></figure>

angles on the dual triangular mesh in the transition zone in order to move the triangle edges closer to the centers of the hexagon edges they intersect.

## 3.5 Constrained centroidal Voronoi tessellation (CCVT) meshes

Given a set of $N$ distinct points on the sphere $x_i$ (referred to as the generators, $1 \le i \le N$), the Voronoi tessellation (or the Voronoi diagram) associated with the generators is the set of polygons $\Omega_i$ consisting of all points that are closer (in the sense of great-circle distance) to $x_i$ than any other $x_j$ with $i \ne j$ (Okabe et al., 2009). For a given set of generators, this tiling is unique and completely covers the sphere, and thus can be employed in conjunction with many finite volume methods. However, for an arbitrary set of generators, it is easy to produce highly distorted polygons, particularly if the density of generators varies substantially. This has led to the development of the constrained centroidal Voronoi tessellation (CCVT) (Du et al., 2003), which imposes the additional requirement that the set of generators be coincident with the centroids of each polygon. Given a desired polygonal density function, several algorithms have been developed to generate CCVTs both in Cartesian and spherical geometry (i.e., for ocean basins or ice sheets) (Ringler et al., 2008). Figure 3 depicts one such CCVT grid that is compatible with the MPAS model. CCVT grids are often confused with deformations of the icosahedral (hexagonal) grid described in Sect. 3.4, since both typically contain a large number of hexagonal elements; however, CCVT grids are fundamentally constructed using a

<figure><img src="image_6.png" /><figcaption>**Figure 3.** A constrained centroidal Voronoi tessellation mesh with localized grid density that could be employed in the MPAS model.</figcaption></figure>

very different technique. Although hexagons are, by far, the most common polygon on CCVT grids, CCVT grids on the sphere will also include at least 12 pentagons and sometimes other polygons with more than six sides. Quadrilateral elements are theoretically possible but are never found in practice on the final grid due to this being a locally unstable solution of the underlying CCVT system of equations.

## 3.6 Octahedral reduced Gaussian grid

As with the classical reduced Gaussian grid of Hortal and Simmons (1991), the octahedral reduced Gaussian grid (Malardel et al., 2016; Smolarkiewicz et al., 2016) specifies the latitudes according to the roots of the Legendre polynomials. The two grids differ in the arrangement of the points along the latitudes, which follows a simple rule for the octahedral grid: starting with 20 points on the first latitude around the poles, 4 points are added with every latitude towards the Equator, whereby the spacing between points along the latitudes is uniform and there are no points at the Equator. The octahedral reduced Gaussian grid is suitable for transformations involving spherical harmonics and has been introduced for operational weather prediction with the spectral dynamical core of the IFS at ECMWF in 2016. Figure 4 depicts the octahedral reduced Gaussian grid nodes together with the edges of the primary mesh as applied in the context of the finite-volume discretization of FVM (Sect. 2.5).

## 3.7 Yin-Yang grid

The overset Yin-Yang grid (Kageyama and Sato, 2004) has two Cartesian grid components (subsets of a latitude-longitude grid) which are geometrically identical (see Fig. 5). These components are combined to cover a spherical surface with partial overlap along their borders. The Yin component covers the latitude-longitude region

$$ \left(-\frac{\pi}{4} - \delta_{\theta} \le \theta \le \frac{\pi}{4} + \delta_{\theta}\right) \cap \left(-\frac{3\pi}{4} - \delta_{\lambda} \le \lambda \le \frac{3\pi}{4} + \delta_{\lambda}\right), $$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Locations of the octahedral reduced Gaussian grid nodes (a), and the edges of the primary mesh connecting the nodes as applied with the finite-volume discretisation in FVM (b). A coarse octahedral grid with only 24 latitudes between each pole and the Equator ("O24") is used for illustration. The dual mesh resolution of the octahedral reduced Gaussian grid is about a factor of 2 finer at the poles than at the Equator; see Smolarkiewicz et al. (2016).</figcaption></figure>

(1)

where $\delta_{\lambda}$, $\delta_{\theta}$ are small buffers that are proportional to the respective grid spacings and are required to enforce a minimum overlap in the overset methodology. For instance, a common configuration employed by the GEM model for DCMIP fixes $\delta_{\theta} = 2^{\circ}$ and $\delta_{\lambda} = 3\delta_{\theta}$. The Yang component covers an analogous area but is rotated perpendicularly so as to cover the region of the sphere outside of the Yin grid. This grid is employed by the GEM model, utilizing a pair of local area models based with the numerics from the GEM latitude-longitude model.

## 3.8 Horizontal staggering

The horizontal placement of variables impacts a number of properties of the numerical method, including how energy and enstrophy conservation is managed, any computational modes that might arise due to differencing, dispersion properties, and the maximum stable time-step size for explicit time-stepping schemes (Randall, 1994; Ullrich, 2014b). The original four Arakawa grids (Arakawa and Lamb, 1977), denoted with letters A through D, were initially designed for rectilinear meshes but were later adapted for a variety of unstructured grids. Later, other grid types were added, including the Z grid, which used the vertical component of vorticity and the horizontal divergence in place of the velocity components (Randall, 1994), and the ZM grid, which extends the B grid to hexagons by placing the velocity at hexagonal nodes (Ringler and Randall, 2002). By interpreting "staggerings" to be analogous to a choice of finite element basis, new staggerings are under development in the context of mixed finite element methods (Cotter and Shipton, 2012). Among the models that participated in DCMIP, only four grids were represented: the A grid, which involves simple co-location of

all velocity components and scalar fields; the C grid, which places perpendicular velocity components on grid edges; the D grid, which places parallel velocity components on grid edges; and the Z grid, which co-locates the vorticity, divergence, and buoyancy variables (see Fig. 6).

Arguments in favor or against particular staggerings have generally emerged from linear analyses and typically in the absence of either implicit or explicit diffusion. In this context, the A grid tends to support large time-step sizes but produces unphysical phase speeds and negative group velocities at high wavenumbers, including a stationary $2\Delta x$ wavelength mode (even in the context of finite element methods); the C grid better represents short wave modes and does not support extraneous computational modes (as long as the number of horizontal faces is equal to twice the number of volumes) but typically has a more restrictive time step with explicit time-stepping schemes than the A grid; the D grid provides a better representation of vorticity but produces unphysical effects analogous to those on the A grid at high wavenumbers that must be controlled with divergence damping; finally, the Z grid yields optimal dispersion properties but requires the inversion of a Poisson problem at each time step to extract the velocity field from the divergence and vorticity.

Other specialized staggerings have been developed that couple horizontal staggering with the formulation of the time integrator. In the FV³ model, although velocities are stored in accordance with the D-grid arrangement, at the intermediate stages of the forward-backward time-stepping scheme, velocities are actually prognosed on the C grid. The intermediate velocities then act as a simplified Riemann solver: the intermediate stage velocities are time centered and can be used to compute the fluxes and advance the flux terms by a full acoustic time step. More details on this approach can be found in Lin and Rood (1997).
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The Yin-Yang grid is a combination of two limited-domain latitude-longitude grids assembled to provide complete coverage of the sphere.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Horizontal staggering options represented among DCMIP models, in this case depicted on a rectilinear grid and geodesic grid. Here, $\eta$ denotes the buoyancy variable.</figcaption></figure>

# 4 Vertical discretization

Because of the vast differences between horizontal and vertical scales in global simulations, most atmospheric models use dimension splitting in order to separate the horizontal discretization from the vertical discretization. In this section, design considerations related to the vertical column are discussed, including the staggering of prognostic and diagnostic variables, and the choice of vertical coordinate.

## 4.1 Vertical staggering

Along with the choice of prognostic variables, the vertical discretization of the equations of motion also allows for the staggered placement of prognostic variables. As with hydrostatic models, certain discretizations give rise to spurious computational modes that can contaminate the solution (Tokioka, 1978; Arakawa and Moorthi, 1988). The choice of vertical staggering may also impact many phys-

ically relevant properties of the model near the grid scale, such as the phase speed of Rossby waves (Thuburn and Woollings, 2005). Finally, the choice of vertical staggering can have impacts on the physics-dynamics coupling (Holdaway et al., 2013a, b). Taken altogether, these issues suggest care should be taken when selecting the discretization. Since co-located discretizations of the non-hydrostatic equations generally require some additional effort to control spurious computational modes, it is more common to employ either (a) a Lorenz-type staggering (Lorenz, 1960), which places horizontal velocity, buoyancy, and thermodynamic variables on model levels, and vertical velocity on model interfaces; or (b) a Charney-Phillips-type staggering (Charney and Phillips, 1953), which places horizontal velocity and buoyancy variables on model levels and vertical velocity and thermodynamic variables on model interfaces (see Fig. 7). These approaches can be further augmented as needed, for instance, by shifting the vertical velocity and thermodynamic
### 4.2.3 Floating Lagrangian coordinates (ACME-A, DYNAMICO, and FV³)

In the floating Lagrangian formulation (Starr, 1945; Lin, 2004), the vertical coordinate is chosen to represent an artificial tracer with monotonically increasing or decreasing mixing ratio $s$ in the vertical. The actual mixing ratio at initiation is arbitrary and can be constructed to be height-like (i.e., $s = z$) or mass-like, i.e.,

$$s = \int_{z}^{\infty} \rho_{0}g dz, \qquad (7)$$

in which case a 3-D reference density field $\rho_0$ can be imposed. Of primary importance is the fact that the vertical coordinate satisfies

$$\dot{s} = \frac{ds}{dt} = 0, \qquad (8)$$

which greatly simplifies the associated prognostic velocity and continuity equations. Floating Lagrangian coordinates are often paired with a vertical remapping operation that corrects for strong grid distortions that may occur after sufficiently long model integrations.

### 4.2.4 Cut cells in OLAM

A pure $z$ coordinate with horizontal grid levels is used in OLAM (Walko and Avissar, 2008b) in order to completely avoid topographic imprinting on the model grid levels (Fig. 8). This implies that grid levels intersect the topographic surface, leading to some grid cells being partially above and partly below the surface. The face areas of these so-called cut cells are reduced accordingly, which in turn regulates cell-to-cell flux transport in accordance with the kinematic constraint imposed by the topography. Cut-cell volumes are also reduced, and volumes and surface areas of all cells appear explicitly in the finite-volume formulation of the mass and momentum conservation equations. One or more methods are used to avoid the so-called small cell problem where volume to area ratios of cut cells are much smaller than those for full cells and therefore can lead to instability. The smallest cells are eliminated by adjusting topography slightly, which is usually justified by noting that local topographic sampling is approximate. In larger cut cells, volumes can be increased (without changing surface areas) which stabilizes the cell at the expense of slowing its response to advected transients. When either of the above adjustments is unacceptable for a particular application, a flux-balance method based partly on Berger and Helzel (2012) is used to stabilize small cut cells.

## 5 Prognostic equations and treatment of moisture

The Navier-Stokes equations that govern atmospheric motion can take on many forms, depending on the choice of

prognostic variables and coordinate system. A derivation of many forms of these equations can be found in Appendix A. The particular prognostic equations used by the model can impact the presence of computational modes, the accuracy of the model in representing the physical modes of the atmosphere (Thuburn and Woollings, 2005), and the ability of the model to conserve important invariants such as energy (Dubos and Tort, 2014). The remainder of this section gives some specific examples of prognostic equations used by the DCMIP models, including any special treatment of terms related to moist physics.

### 5.1 ACME-A

ACME-A presently solves the compressible shallow-atmosphere equations using a hybrid terrain-following pressure vertical coordinate $\eta$, similar to the model of Laprise (1992). The 2-D vector-invariant form of the prognostic horizontal velocity Eq. (A62) is employed, in conjunction with prognostic potential temperature (Eq. A57), pseudo-density (Eq. A55), and geopotential (Eq. A27). The vertical velocity equation is formulated analogous to that of GEM:

$$\frac{dw}{dt} = -g_c \left( 1 - \frac{\partial p}{\partial s} \right). \qquad (9)$$

### 5.2 CSU

The CSU model uses the vorticity divergence form of the equations of motion, as described in Sect. A10, discretized on the geodesic mesh with absolute vorticity and velocity divergence scalars stored at cell centers. The unified approximation of the equations of motion (Arakawa and Konor, 2009) is employed to avoid vertically propagating sound waves.

### 5.3 DYNAMICO

The prognostic equations employed by DYNAMICO are based on a Hamiltonian formulation (Dubos and Tort, 2014). The specific prognostic variables employed are pseudo-density $\rho_s$, mass-weighted tracers (potential temperature, water species), geopotential $\Phi$, horizontal covariant components of momentum, and mass-weighted vertical momentum $W = \rho_s g^{-2} d\Phi/dt = \rho_s g^{-1} w$. Prognostic equations are in flux form for mass (Eq. A55) and $W$ (Eq. A23), in advective form for $\Phi$ (Eq. A27), and in vector-invariant form for covariant horizontal momentum (Eq. A76).

### 5.4 FV³

The hydrostatic FV³ model uses a mass-based floating Lagrangian coordinate along with the shallow-atmosphere approximation (Lin, 2004). Prognostic equations include horizontal velocity in 2-D vector-invariant form (Eq. A38), pseudo-density (Eq. A55), and virtual potential temperature (Eq. A57). The non-hydrostatic model further incorporates
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) A terrain-following coordinate passing over rough topography. (b) A cut-cell coordinate used for representing the same topography.</figcaption></figure>

prognostic geopotential (Eq. A27) and vertical momentum (Eq. A37).

## 5.5 FVM

The FVM formulation is based on conservation laws for dry mass (Eq. 10a), momentum (Eq. 10b), and dry entropy (Eq. 10c) in Eulerian flux form, which are similar to Eq. (A9) for $ρ_d$, Eq. (A23), and Eq. (A13) for $θ$, respectively. Moreover, underlying the conservation laws in FVM is a perturbational form with respect to a balanced ambient state and a generalized curvilinear coordinate formulation in a geospherical framework. Following Smolarkiewicz et al. (2017), the FVM governing equations can concisely be written as

$$ \frac{\partial \mathcal{G} \rho_d}{\partial t} + \bar{\nabla} \cdot (\mathbf{v} \mathcal{G} \rho_d) = 0, \quad (10a) $$

$$ \frac{\partial \mathcal{G} \rho_d \mathbf{u}}{\partial t} + \bar{\nabla} \cdot (\mathbf{v} \mathcal{G} \rho_d \mathbf{u}) = \mathcal{G} \rho_d \left( -\theta_{\rho} \tilde{\mathbf{G}} \bar{\nabla} \phi' - \mathbf{k} g \left( \frac{\theta'}{\theta_a} + \varepsilon_b q'_v - q_c - q_r \right) - 2\Omega \times \left( \mathbf{u} - \frac{\theta_{\rho}}{\theta_{\rho a}} \mathbf{u}_a \right) + \mathcal{M}' \right], \quad (10b) $$

$$ \frac{\partial \mathcal{G} \rho_d \theta'}{\partial t} + \bar{\nabla} \cdot (\mathbf{v} \mathcal{G} \rho_d \theta') = - \mathcal{G} \rho_d \tilde{\mathbf{G}}^T \mathbf{u} \cdot \bar{\nabla} \theta_a, \quad (10c) $$

$$ \phi' = c_{pd} \left[ \left( \frac{R_d}{p_0} \rho_d \theta (1 + q_v / \varepsilon) \right)^{R_d / c_{vd}} - \pi_a \right]. \quad (10d) $$

Dependent variables in Eq. (10) are dry density $ρ_d$, 3-D physical velocity vector $\mathbf{u}$, potential temperature perturbation $θ'$, and a modified Exner pressure perturbation $\phi'$, with the thermodynamic variables related by the gas law Eq. (10d). Primes indicate perturbations with respect to the prescribed ambient state denoted by subscript “a”; see Prusa et al. (2008) and Smolarkiewicz et al. (2014) for discussions. The symbol $g$ in Eq. (10b) denotes the gravitational acceleration and $\varepsilon_b = 1/\varepsilon - 1$. As far as geometric aspects are concerned, the nabla operator $\bar{\nabla}$ represents the 3-D vector of

partial derivatives with respect to the curvilinear coordinates, along with the Jacobian $\tilde{\mathbf{G}}$, a matrix of metric coefficients $\tilde{\mathbf{G}}$, its transpose $\tilde{\mathbf{G}}^T$, and the contravariant velocity $\mathbf{v} = \tilde{\mathbf{G}}^T \mathbf{u}$, where a contribution from optional time dependency of the curvilinear coordinates is neglected for simplicity (Kühnlein and Smolarkiewicz, 2017). The symbol $\mathcal{M}' = \mathcal{M}'(\mathbf{u}, \mathbf{u}_a, \theta_{\rho} / \theta_{\rho a})$ in Eq. (10b) subsumes the metric forces in the spherical domain (Smolarkiewicz et al., 2017).

## 5.6 GEM

In GEM, the non-hydrostatic equations are written explicitly as deviations from hydrostatic balance represented by

$$ \mu = \frac{\partial p}{\partial s} - 1, \quad (11) $$

where $s$ (denoted $\pi$ in GEM documentation) is given by Eq. (3). In this case, the equations of GEM model (Girard et al., 2014) are concisely given by

$$ \frac{d\mathbf{u}_h}{dt} + f\mathbf{k} \times \mathbf{u}_h + R_d T_v \nabla_{\zeta} \log p + (1+\mu) \nabla_{\zeta} \Phi = 0, \quad (12) $$

$$ \frac{dw}{dt} - g_c \mu = 0, \quad (13) $$

$$ \frac{d}{dt} \log \left( \frac{\partial s}{\partial \zeta} \right) + \nabla_{\zeta} \cdot \mathbf{u}_h + \frac{\partial \dot{\zeta}}{\partial \zeta} = 0, \quad (14) $$

$$ \frac{d \log T_v}{dt} - \frac{R_d}{c_p} \frac{d \log p}{dt} = 0, \quad (15) $$

$$ \frac{\partial \Phi}{\partial s} + \frac{R_d T_v}{p} = 0, \quad (16) $$

$$ \frac{d\Phi}{dt} - g_c w = 0. \quad (17) $$

Here, $\nabla_{\zeta}$ denotes the horizontal gradient along $\zeta$ surfaces. With respect to the treatment of moisture in GEM, the cloud water and all non-gases are embedded in the total air density $ρ$, affecting the virtual temperature defined in Eq. (A7). Also, specific mass is used in GEM (not mixing ratio).
## 5.7 ICON

ICON solves a non-hydrostatic equation set based on Gassmann and Herzog (2008) using terrain-following z coordinates. The governing equations describe the mixture of a two-component system of dry air and water, where water is allowed to occur in all three phases, including precipitating constituents. Following Wacker et al. (2006), the barycentric (bc) velocity $\mathbf{u}_{bc} = \sum_k \rho_k \mathbf{u}_k / \sum_k \rho_k$ – that is, the mass-weighted sum of all constituent-specific velocities (including sedimentation velocities) – is used as a prognostic variable. In contrast to Gassmann and Herzog (2008), a vector-invariant form is only used for the horizontal velocity Eq. (A33), whereas the vertical velocity equation is solved in advective form. The pressure-gradient force is formulated according to Eq. (A20).

Additional prognostic variables include total air density (Eq. A10), virtual potential temperature (Eq. A57), and mass fractions $q_k = \rho_k / \rho$ of all constituents (except for dry air) for which the prognostic equation reads

$$ \frac{\partial \rho q_k}{\partial t} + \nabla \cdot (\rho q_k \mathbf{u}_{bc}) = -\nabla \cdot \mathbf{J}_k + \sigma_k, \quad (18) $$

with $\sigma_k$ describing sources/sinks due to phase changes, and $\mathbf{J}_k = \rho q_k (\mathbf{u}_k - \mathbf{u}_{bc})$ denoting diffusion fluxes, which account for the motion of constituents relative to the frame of reference set by $\mathbf{u}_{bc}$.

The specific heat capacities and ideal gas constant are approximated to be equal to their dry values $R^* \approx R_d$, $c_p^* \approx c_{pd}$, $c_v^* \approx c_{vd}$. The model also uses a prognostic equation for Exner pressure to simplify the treatment of vertical sound wave propagation, given by

$$ \frac{\partial \pi}{\partial t} + \frac{R_d}{c_{vd}} \frac{\pi}{\rho \theta_v} \nabla \cdot (\mathbf{u}_{bc} \rho \theta_v) = \hat{Q}, \quad (19) $$

where $\hat{Q}$ is an appropriately formulated diabatic heat term. The horizontal uses a Arakawa C-grid formulation on the triangular grid to prognose horizontal velocities normal to triangle edges $v_n$, making use of reconstructed tangential velocity components $v_t$.

In the current implementation, the following simplifications are made with regard to the treatment of moisture: the atmospheric mass loss/gain due to precipitation/evaporation is neglected in the total mass continuity Eq. (A10) by setting the vertical component of $\mathbf{u}_{bc}$ to zero at the lower boundary: $w_{bc|sfc} = 0$. In addition, only the vertical diffusion fluxes $\mathbf{J}_k$ of sedimentation constituents and the surface evaporation flux $\mathbf{J}_v|sfc$ are taken into account. The counter flux of non-sedimentation constituents is discarded. Since in the given framework the continuity Eq. (A10) is only valid if the constraint $\sum_k \mathbf{J}_k = 0$ holds, it is (implicitly) assumed that a fictitious counter flux of dry air compensates for the considered vertical diffusion fluxes. As a consequence, ICON currently conserves the global integral of total air mass rather than dry air mass.

## 5.8 MPAS

The evolution equations used by MPAS are fully described in Skamarock et al. (2012), based on the formulation of Dutton (1986). The MPAS model uses the momentum form of the update equations, as described in Sect. A11, with dry mass utilized for the density variable $\tilde{\rho}_s$. MPAS further evolves dry mass using a continuity equation of Eq. (A10) and moist potential temperature following Eq. (A13).

## 5.9 NICAM

NICAM prognoses horizontal and vertical momentum analogous to the approach described in Sect. A11. It further evolves the density perturbation from the background reference state using Eq. (A10) and sensible heat part of internal energy. A detailed explanation of the evolution equations can be found in Satoh et al. (2008).

## 5.10 OLAM

OLAM solves the deep-atmosphere, fully compressible equations in mass- and momentum-conserving finite-volume form using Eqs. (A10), (A23), and (A13). Prognostic variables are the three components of momentum, ice-liquid potential temperature $\theta_{il}$ (Walko et al., 2000), total density $\rho$, specific density of total water, and specific bulk density and/or bulk number concentration of various scalar quantities including liquid and ice hydrometeors, aerosols, and trace gases. For DCMIP, the latter are limited to cloud and rain specific bulk densities. Water vapor density is diagnosed by subtracting bulk densities of all liquid and ice hydrometeors from the total water density, dry air density is diagnosed by subtracting total water density from total density, and pressure is diagnosed based on the equation of state and values of dry air density, water vapor density, and potential temperature $\theta$. The latter is in turn diagnosed from $\theta_{il}$ and from the latent heat required to convert any hydrometeors present to the vapor phase. Velocity components are diagnosed by dividing momentum components by total density.

Momentum is C-staggered in the horizontal and vertical (Lorenz vertical staggering is used), meaning that prognosed components live on the grid cell faces and are each normal to the respective face, and the pressure-gradient force is evaluated and applied at those locations. However, evaluation of advective and turbulent momentum transport (as well as the Coriolis force) involves a diagnostic reconstruction of the total momentum vector at the centers of scalar grid cells (Perot, 2000), and cell-to-cell flux of momentum is computed from that reconstruction using the same A-grid control volumes as for scalars. This arrangement is particularly convenient for the cut-cell formulation at the topographic surface where reduced cell face areas and volumes regulate momentum and scalar fluxes in an identical manner.
## 5.11 Tempest

Tempest is a shallow-atmosphere Eulerian model with terrain-following z coordinates with prognostic density (Eq. A55), virtual potential temperature (Eq. A57), and vector-invariant form for covariant horizontal velocity (Eq. A76) and vertical momentum (Eq. A32).

# 6 Diffusion, stabilization, filters, and fixers

Most dynamical cores implement specialized techniques for diffusion or stabilization (see Table 5). Diffusion is a numerical technique that removes spurious numerical noise from the simulation, where the numerical noise typically arises because of inaccuracies in the treatment of waves with wavelengths near the grid scale. Diffusion also includes mechanisms for damping vertically propagating internal gravity waves, such as model-top Rayleigh layers, which are fairly ubiquitous across models and hence not discussed in detail here. Stabilization is a numerical technique that prevents energy growth and allows the model to be run over long periods. Diffusion or stabilization options include physically motivated turbulence parameterizations, added viscosity or hyperviscosity terms with tunable coefficients, off-centering, or wave-mode filters. Since the discretization can also lead to an unphysical loss of mass or energy, mass or energy fixers are also employed to replace lost mass or energy to the system. A comprehensive overview of schemes for diffusion and stabilization schemes can be found in Jablonowski and Williamson (2011). In this section, we discuss some of the diffusion and stabilization strategies employed by the DCMIP suite of dynamical cores.

## 6.1 ACME-A/Tempest

In both ACME-A and Tempest, scalar hyperviscosity is employed for $\rho$, $\theta$, and tracer variables via repeated application of a scalar Laplacian (Dennis et al., 2012; Ullrich, 2014a). Vector hyperviscosity is also applied by decomposing the horizontal vector Laplacian into divergence damping and vorticity damping terms via the vector identity

$$ \nabla^2 \mathbf{u}_h = \nabla \nabla \cdot \mathbf{u}_h + \nabla \times \nabla \times \mathbf{u}_h. \quad (20) $$

Both viscosity operations are applied after the completion of all Runge-Kutta subcycles. Several limiter options are available for tracer transport, including a sign-preserving limiter and a monotone optimization base limiter described in Guba et al. (2014).

## 6.2 CSU

The CSU model utilizes an explicit diffusion scheme that consists of fourth-order hyperdiffusion ($\nabla^4$) applied to the vorticity, divergence, and potential temperature. The model

does not include any explicit diffusion in the vertical column. However, for the idealized DCMIP test cases, explicit diffusion was disabled.

## 6.3 DYNAMICO

In DYNAMICO, (hyper-)diffusive filters are used to eliminate spurious noise due to the energy-conservative centered discretization. Filters are applied every $N_{\text{diff}}$ Runge-Kutta time step in a forward-Euler manner, with $N_{\text{diff}}$ as large as allowed by stability. The scalar Laplacian is computed as the divergence of the gradient, and the vector Laplacian is decomposed into its divergent (grad div) and rotational (curl curl) parts. The strength of filtering is controlled by dissipation timescales $\tau$: given $\tau$, the hyperviscous coefficient that multiplies operator $D^p$ is $\delta^{2p}\tau^{-1}$, where $\delta^{-2}$ is the largest eigenvalue of operator $D$. For DCMIP, DYNAMICO uses $p=2$ (fourth-order hyperviscosity) for all filters.

## 6.4 FV³

Explicit dissipation in FV³ is applied separately to the divergence and to the horizontal fluxes in the governing equations. The D-grid discretization applies no direct implicit dissipation to the divergence, so divergence damping is an intrinsic part of the solver algorithm since otherwise there are no processes by which energy contained in the divergent modes is removed at the grid scale. FV³ has options for fourth-, sixth-, or eighth-order divergence damping; a second-order option is also available for use in idealized convergence tests, which can be applied in addition to the higher-order diffusion. The monotonicity constraint used in computing the fluxes in the momentum, thermodynamic, and mass continuity equations is sufficient to damp and stabilize the non-divergent component of the flow. The model also supports an option to apply hyperdiffusion to the fluxes in each of these equations, with the exception of the tracer transport, which always uses monotonic transport with no explicit diffusion. The hyperdiffusion is of the same order as but much smaller than the divergence damping. Both divergence damping and hyperdiffusion are applied along the Lagrangian surfaces and are recomputed every acoustic time step.

FV³ is constructed with a flexible-lid (constant-pressure) upper boundary that is effective at damping internal gravity wave modes; however, FV³ also applies second-order diffusion to all fields, except the tracers, to create a sponge layer, typically comprising the top two layers of the domain, to damp other signals reaching the top of the domain. An energy-conserving Rayleigh damping is also available, applied consistently to all three components of the winds, which is strongest in the top layer of the domain and becomes weaker with distance until it reaches a runtime-specified cutoff pressure.

FV³ has an option to restore lost energy by the adiabatic dynamics, in whole or a fraction thereof (decided by
a namelist option at runtime), by globally adding a Exner-function weighted potential temperature increment. This is only done before the physics is called and is not used in idealized simulations.

## 6.5 FVM

Within the dynamical core, FVM does not apply any explicit dissipation/diffusion. For the DCMIP test cases, the implicit regularization of the monotonic MPDATA provides sufficient dissipation/diffusion needed to remove excess energy from the finest scales and maintain model stability. An absorbing layer is also available for damping vertically propagating waves near the model top.

## 6.6 GEM

An explicit hyperviscosity in GEM is handled via applications of the Laplacian operator for both wind components and tracers. A vertical sponge layer, which uses a Laplacian operator, is employed on wind components and $T_v$ with a vertical modulation on the topmost levels. For stabilization purposes, the temporal discretization of GEM also uses an off-centering parameter. The quasi-monotone, semi-Lagrangian (QMSL) method (Bermejo and Staniforth, 1992) is used operationally to ensure tracer monotonicity for specific humidity and different hydrometeors. Other options are now available in GEM, including a mass-conserving monotonic scheme (Sørenson et al., 2013) and a global mass fixer (Bermejo and Conde, 2002). Those approaches have been evaluated using chemical constituents such as ozone (de Grandpré et al., 2016).

## 6.7 ICON

The ICON model employs damping and diffusion operators for numerical stabilization and dynamic closure. The details of this scheme appear in Sect. 2.4 and 2.5 of Zängl et al. (2015), and are summarized here. For damping, in the corrector step, a fourth-order divergence damping term $F_d(v)$ is applied in order to allow calling the (relatively) computationally expensive diffusion operator (see below) at the physics time steps without incurring numerical stability problems under extreme conditions:

$$F_d(v) = -f_d \bar{a}_c^2 \nabla \tilde{\nabla} \cdot \left\{ \nabla \left[ \tilde{\nabla} \cdot v + \frac{\Delta}{\Delta z} \left( w - \overline{w_{cc} c^i} \right) \right] \right\}. \quad (21)$$

$f_d$ typically attains values between $\frac{1}{1000 \Delta t}$ and $\frac{1}{250 \Delta t}$, and $\bar{a}_c$ is the global mean cell area.

ICON also includes Rayleigh damping on $w$ following Klemp et al. (2008), which serves to prevent unphysical reflections of gravity waves at the model top. The Rayleigh damping is restricted to a fixed number of levels below the model top, and the damping coefficient is given by a hyperbolic tangent.

The horizontal diffusion consists of a flow-dependent second-order Smagorinsky diffusion of velocity ($F_{D2}(v_n)$) and potential temperature ($F_{D2}(\theta)$) combined with a fourth-order background diffusion of velocity $F_{D4}(v_n)$, defined via

$$ \begin{aligned} F_{D2}(v_n) &= 4K_h \tilde{\nabla}^2(v_n), & F_{D2}(\theta) &= a_c \tilde{\nabla} \cdot \left( K_h \frac{\Delta\theta}{\Delta n} \right), \\ F_{D4}(v_n) &= -k_4 a_c^2 \tilde{\nabla}^2(\tilde{\nabla}^2(v_n)), & & \end{aligned} \quad (22) $$

where $a_c$ and $a_e$ denote the area associated with the cell and edge under consideration, respectively. An empirically determined offset of $0.75k_4a_e$ is subtracted from $K_h$ in order to avoid excessive diffusion under weakly disturbed conditions.

A fourth-order computational diffusion is also available for vertical wind speed $w$. This filter term is needed at resolutions of O(1 km) or finer because the advection of vertical wind speed has no implicit damping of small-scale structures. This term appears as

$$F_D(w) = -k_w a_c^2 \nabla^2(\nabla^2(w)). \quad (23)$$

## 6.8 MPAS

The MPAS model applies fourth-order hyperdiffusion and Smagorinsky diffusion (Smagorinsky, 1963), as described in Skamarock et al. (2012). When applied to the momentum, the Laplacian is evaluated as

$$\nabla^2 u_i = \frac{\partial}{\partial x_i} \nabla_s \cdot v - \frac{\partial \eta}{\partial x_j}, \quad (24)$$

where $u_i$ is the edge-normal velocity defined on cell edge $i$, $\eta$ is the vertical component of the relative vorticity, computed on vertices, and $\nabla_s \cdot v$ is the horizontal divergence on $s$ surfaces, computed on edges. The evaluation of divergence and vorticity in this expression is described in Ringler et al. (2010). The fourth-order hyperdiffusion operator is then computed by twice applying the above Laplacian operator to the momentum.

Smagorinsky diffusion, which is often applied in atmospheric models to parameterize turbulent processes, uses a second-order Laplacian and a physically motivated eddy viscosity $K_h$, defined in terms of Cartesian velocities $(u, v)$:

$$K_h = c_s^2 \ell^2 \sqrt{(u_x - v_y)^2 + (u_y + v_x)^2}, \quad (25)$$

where $c_s$ is a constant parameter and $\ell$ is the grid scale. The diffusion operator then takes the form $\nabla \cdot (K_h \nabla \psi)$ for a scalar field $\psi$.

## 6.9 NICAM

NICAM implements three types of diffusion: 3-D divergence damping, fourth-order horizontal hyperdiffusion, and sixth-order vertical hyperdiffusion, as described in Tomita and Satoh (2004). Specifically, the divergence damping term
(Skamarock and Klemp, 1992) aims to suppress instabilities that arise due to the time-splitting scheme and is applied to both horizontal and vertical velocities. The hyperdiffusion operators are applied to all prognostic variables. For tracer advection, upwinding is used to remove spurious oscillations, as described in Miura (2007) and Niwa et al. (2011).

## 6.10 OLAM

OLAM requires two types of artificial damping. In the upper layers of the model, vertical velocity and small-scale horizontal divergence are damped in order to attenuate gravity waves and thereby mitigate their reflection off the rigid top boundary of the domain. The damping layer is commonly applied in the uppermost 10 km of the domain, where the model top is 35 or 40 km above sea level. The damping rate is zero at the bottom of the damping layer and increases upward, usually linearly. Throughout the model domain, vertical vorticity is filtered horizontally at the smallest resolvable scale in order to control a spurious computational mode. This vertical vorticity mode is inherent in C-staggered momentum formulations on hexagonal meshes because horizontal velocities are more numerous than twice the number of scalar (mass) values and are thus underconstrained (Weller et al., 2012; Weller, 2012). The vorticity filter is constructed so as to have zero impact on divergence at any scale. Upwinding in the Lax–Wendroff formulation of the advection operator provides sufficient damping so that no other type of filtering is required.

# 7 Temporal discretizations

Temporal discretizations are important for capturing the discrete dynamical evolution of the global atmosphere. In the past two decades, a variety of new temporal discretizations have been developed, leaving behind the days when the leapfrog scheme was ubiquitous across models. This diversity is in part because of the demands of non-hydrostatic models: unlike their hydrostatic counterparts, non-hydrostatic atmospheric models must include a mechanism for dealing with vertically propagating sound waves. These waves are meteorologically insignificant, but with a vertical grid spacing of 100 m, a purely explicit temporal discretization of the unmodified fluid equations would require a time-step size on the order of 1 s or less. Consequently, sound waves are either filtered explicitly through the use of an alternative equation set or artificially slowed through the use of implicit temporal discretizations. Some commonly employed alternative equation sets include the anelastic (Ogura and Phillips, 1962), quasi-hydrostatic (Orlanski, 1981), pseudo-incompressible (Durran, 1989), or unified approximation (Arakawa and Konor, 2009). These filtered equation sets generally require that a global elliptic solve be performed as prognostic variables are updated. In

this section, we discuss the time-stepping schemes that have been employed across the DCMIP suite of models.

## 7.1 Mixed implicit–explicit, forward–backward, semi-implicit, and additive Runge–Kutta schemes

Implicit–explicit schemes are a broad category of time-integration schemes that divide the terms of the prognostic equations into a set of explicitly integrated terms and implicitly integrated terms. At the very least, terms associated with vertically propagating sound waves are included among the implicit terms. For the remaining terms, there is some freedom in choosing how to integrate terms associated with vertical advection and horizontally propagating sound waves. Semi-implicit schemes are one such class of schemes that typically incorporate horizontally propagating sound waves into the implicit solve and thus rely on a global Helmholtz-type solve. Additive Runge–Kutta schemes are another mechanism to ensure high-order temporal accuracy, and many such schemes have been described throughout the literature (see, for example, Weller et al., 2013; Ullrich and Jablonowski, 2012b). Several examples of these schemes can be found among the DCMIP models.

ACME-A and Tempest both use the ARS(2,3,2) scheme described in Ascher et al. (1997), with all horizontal and vertical advection terms treated explicitly and the remaining vertical terms, associated with sound wave propagation, treated implicitly. A number of different ARK schemes have been compared and contrasted in this framework, with significant implications for model performance and stability (Gardner et al., 2017).

CSU uses a semi-implicit time-integration scheme with third-order Adams–Bashforth scheme for explicit integration of the continuity equation, potential temperature equation, and terms related to advection. Since potential temperature is updated prior to the computation of the pressure-gradient force, this term can be thought of as implicit in time. The horizontal wind field is then predicted through integration of the vorticity and divergence of the horizontal wind and a multi-grid method applied to solve a pair of two-dimensional Poisson equations for the stream function and velocity potential, which are then differentiated to obtain the velocity field. Horizontal diffusion is then applied forward in time.

FV³ and its predecessors are integrated using a forward–backward integration for the Lagrangian dynamics. With the exception of the pressure-gradient force, all of the terms in the momentum, energy, and mass equations are expressible as fluxes and thus can be integrated using the explicit forward-in-time algorithm described by Lin and Rood (1997). The horizontal component of the pressure-gradient force is evaluated backward in time using the algorithm of Lin (1997); the non-hydrostatic component of the vertical pressure-gradient force is evaluated using a semi-implicit solver. This forward–backward time step is referred to as the “acoustic” time step, although the full solver is advanced on
each of these acoustic time steps. Physics tendencies are applied impulsively at prescribed intervals, consistent with the forward-in-time discretization; the physics time step is typically much longer than the acoustic time step.

DYNAMICO uses an additive Runge-Kutta time scheme with two Butcher tableaus, one explicit and one implicit. A Hamiltonian splitting decides which terms of the equations of motion are treated explicitly or implicitly (Dubos and Dubey, 2017). As a result, the implicit terms couple the vertical acceleration due to the pressure gradient and the adiabatic pressure change due to vertical displacements of fluid parcels. The resulting implicit problem reduces to independent, scalar, purely vertical, non-linear problems which are solved to machine precision in two Newton iterations involving one tridiagonal solve each. The overall time scheme has a HEVI (horizontally explicit, vertically implicit) structure. Currently, the second-order, three-stage ARK(2,3,2) scheme is used (Giraldo et al., 2013).

ICON consists of a two-time-level predictor-corrector scheme, which is explicit for all terms except for those describing the vertical propagation of sound waves. No time splitting is used with respect to sound waves, because the ratio between the speed of sound and the maximum wind speed in the mesosphere, which is in part covered by the vertical domain, can be close to 1. Instead, time splitting is employed to dynamics on the one hand and horizontal diffusion, tracer transport, and fast physics on the other hand. Typically, a full time step consists of four or five dynamical substeps in which a constant forcing originating from the slow physics is applied. Mass-consistent transport is achieved by passing time-averaged air-mass fluxes from the dynamical substeps to the transport scheme. The details of the predictor-corrector scheme, including measures to increase the numerical efficiency and to optimize the accuracy, are described in Sect. 2.4 of Zängl et al. (2015).

MPAS and NICAM use a split-explicit formulation (Klemp et al., 2007) consisting of an outer Runge-Kutta loop (typically RK3) and inner acoustic loop. At the beginning of each Runge-Kutta subcycle, tendencies are computed for each of the prognostic variables and stored for the duration of the subcycle. Several iterations of an acoustic loop are then performed with a time step much smaller than required for the Runge-Kutta subcycle. Within the acoustic loop, an implicit solve for vertically integrated sound waves is performed to avoid time-step constraints that may arise from vertically propagating sound waves.

OLAM uses a unique temporal discretization that combines elements of the Adams-Bashforth (AB2) scheme and a Lax-Wendroff formulation for advected quantities. However, instead of extrapolating all prognostic tendencies forward to the half-future time level as in AB2, the horizontal momentum components alone (not their tendencies) are extrapolated in time at the cell boundaries where they reside. The extrapolated momentum provides the time-centered cell-to-cell total mass flux across the grid cell faces that is re-

sponsible for advective transport. Advection of all quantities, including all three velocity components that are诊断ally reconstructed at scalar cell centers, and advancement in time from the current to the future time level is based on the time- and space-centered Lax-Wendroff formulation. This scheme is horizontally explicit, but a trapezoidal-implicit formulation is used in the vertical for stable integration of vertically propagating sound waves. A byproduct of the implicit formulation is an implicit time-centered vertical momentum that joins the time-extrapolated horizontal momentum to form a complete set of mass fluxes for advection. The vertical momentum equation is solved first so that the time-centered vertical momentum is available for computing transport of horizontal momentum and all scalar quantities. A time-split scheme is most often used where momentum and potential temperature are updated more frequently than other scalar fields in order to accommodate horizontally propagating sound waves.

## 7.2 The FVM semi-implicit method

A characteristic feature of the FVM (Sect. 2.5) time-stepping scheme is the 3-D implicit treatment of the fast buoyant and acoustic modes, and the slow rotational modes. Therefore, the model time step is identical for all processes and typically selected with regard to the stability of the advective transport scheme; i.e., the time step is continuously adapted according to a given maximum advective Courant number permitted by the MPDATA scheme. A comprehensive discussion of the integration scheme can be found in Smolarkiewicz et al. (2014, 2016) and Kühnlein and Smolarkiewicz (2017) for dry dynamics, whereas it can be found in Kurowski et al. (2014) and Smolarkiewicz et al. (2017) for extension to moist-precipitating dynamics. Here, we provide a short outline of the solution procedure for the compressible Euler Eq. (10). It employs the two-time-level, second-order-accurate template algorithm given as

$$ \psi_i^{n+1} = A_i(\tilde{\psi}^n, v^{n+1/2}, G^n, G^{n+1}) + 0.5 \Delta t R^\psi |_i^{n+1}, \\ \tilde{\psi}^n = \psi^n + 0.5 \Delta t R^\psi |_i^n, \qquad (26) $$

where $\psi$ represents the solution variable, $R^\psi$ is the respective right-hand side, $A$ symbolizes the advective transport operator given by the non-oscillatory, finite-volume MP-DATA scheme (Smolarkiewicz and Szmelter, 2005; Kühnlein and Smolarkiewicz, 2017), and the spatial mesh vector index $i \equiv (k, i)$ denotes the position on the hybrid horizontally unstructured, vertically structured computational mesh.

The solution procedure of Eq. (10) can then be divided into three steps. First, the homogenous mass continuity Eq. (10a) is integrated with $\psi \equiv \rho_d$, $v \equiv vG$, $G \equiv G$, and $R^\rho \equiv 0$ in Eq. (26). Second, given already updated moisture variables (Smolarkiewicz et al., 2017), the thermodynamic (Eq. 10c) and momentum (Eq. 10b) equations enter Eq. (26) with $\psi = u, v, w, \theta', v \equiv vG\rho_d, G \equiv G\rho_d$, and the right-hand-side $R^\psi$
which generally depends on all these prognostic variables. The high degree of implicitness in the representation of the right-hand-side forcings is achieved by inverting the overall discrete system (Eq. 26) to obtain closed-form expressions for the velocity updates; this procedure is facilitated by the co-located arrangement of all variables on the computational mesh. Retained on the right-hand side of the derived closed-form velocity expressions is the pressure-gradient term. The subsequent third step in the solution procedure is to formulate an implicit boundary value problem for the pressure variable $\phi'$ using an advective form of the equation of state (Eq. 10d). An $O(\Delta t^2)$ integration of this equation with a Euler backward scheme, in the spirit of Eq. (26), leads to a Helmholtz equation (Smolarkiewicz et al., 2014). The associated 3-D elliptic boundary value problem is solved iteratively using a preconditioned generalized conjugate residual approach; see Smolarkiewicz and Szmelter (2011) for a recent overview and comprehensive list of references. Non-linear terms in $R^{\psi|n+1}$ and the solution-dependent coefficients of the Helmholtz equation are lagged behind and executed in an outer iteration.

## 7.3 A semi-Lagrangian implicit time discretization in the GEM model

GEM differs from the approaches above by using a semi-Lagrangian advection. Any model equations, prognostic or diagnostic, are written in the form

$$ \frac{dF}{dt} + G = 0, \qquad (27) $$

where $d/dt$ is the Lagrangian derivative, $F$ contains the terms subject to this operator, and $G$ contains the remaining terms. The semi-Lagrangian approach consists in the following space-time discretization of Eq. (27):

$$ \frac{F^A - F^D}{\Delta t} + \left(\frac{1}{2} + \epsilon\right) G^A + \left(\frac{1}{2} - \epsilon\right) G^D = 0, \qquad (28) $$

where “A” stands for the arrival position at model grid point $(\mathbf{r}_h, \zeta, t)$ and “D” for the departure position $(\mathbf{r}_h - \Delta\mathbf{r}_h, \zeta - \Delta\zeta, t - \Delta t)$ due to the displacements $\Delta\mathbf{r}_h, \Delta\zeta$ having occurred during the time step $\Delta t$. $G$ is averaged between these two positions with a possible slight off-centering $\epsilon$. The displacements are themselves calculated by solving, again using the Lagrangian method, the following equations:

$$ \frac{d\mathbf{r}_h}{dt} - \mathbf{u}_h = 0; \quad \frac{d\zeta}{dt} - \dot{\zeta} = 0, \qquad (29) $$

discretized in the same way (trapezoidal method) though without off-centering:

$$ \frac{\Delta \mathbf{r}_h}{\Delta t} - \frac{\mathbf{u}_h^A + \mathbf{u}_h^D}{2} = 0; \quad \frac{\Delta \zeta}{\Delta t} - \frac{\dot{\zeta}^A + \dot{\zeta}^D}{2} = 0. \qquad (30) $$

The process is of course a multi-step iterative one, since both positions and velocities at departure positions (past time

$t - \Delta t$) are unknown as well as, of course, the velocities at arrival positions (time $t$). Once a first estimate of the departure positions is obtained, the model equations are solved to obtain a first estimate of the velocities at time $t$. The model equations must be solved simultaneously and this is only possible for the linear part $L$ which becomes a matrix inversion problem. Hence, a suitable linearization is considered. The unknown (arrival) linear $L$ and non-linear $N$ parts are then separated from the known (first departure estimate) remaining $R$ part. Thus, first separating space-time, Eq. (28) is rewritten as follows:

$$ \frac{F^A}{\tau} + G^A = \frac{F^D}{\tau} - \beta G^D \equiv R^D, \qquad (31) $$

where $\tau = (1/2 + \epsilon) \Delta t$ and $\beta = (1/2 - \epsilon)/(1/2 + \epsilon)$. Secondly, separating linear from non-linear parts, we get

$$ L^A + N^A = R^D, \qquad (32) $$

with

$$ L^A = \left[ \frac{F^A}{\tau} + G^A \right]_{\text{linear}}, \quad \text{and} \\ N^A = \frac{F^A}{\tau} + G^A - \left[ \frac{F^A}{\tau} + G^A \right]_{\text{linear}}. \qquad (33) $$

Note that both $F$ and $G$ may require linearization. $L^A$ may then be obtained if $N^A$ is first guessed; once $L^A$ is found, an estimate of $N^A$ is obtained and $L^A$ is recalculated. This is called the non-linear iteration process (one iteration is usually sufficient). The overall process is then repeated once, starting from a new estimate of the departure positions.

There are two intensive calculation sections in this process: the so-called semi-Lagrangian calculations (twice estimating departure positions, twice interpolating right-hand side $R$ on departure positions) and solution of the linear system (four times). Each time, the linear system is reduced to a Helmholtz problem for one composite variable. For this problem, a direct solver is involved using the Schwarz-type domain decomposition method on a Yin-Yang grid (Qaddouri et al., 2008). The composite variable solution is then used to update the prognostic variables (back substitution). At the end of the time step, the static halo region of both panels of the Yin-Yang grid is updated (Qaddouri and Lee, 2011). All required interpolations throughout the semi-Lagrangian process and between Yin and Yang grids are cubic interpolations.

## 8 Summary and conclusions

As discussed earlier, this paper represents the first in a series of papers documenting the results from the 2016 Dynamical Core Model Intercomparison Project workshop and summer school. In this paper, we have provided a description of the differences and similarities between participating
models, including the choice of computational grid, horizon-
tal staggering, vertical staggering, vertical coordinates, prog-
nostic equations, choice of diffusion, stabilization, filters and
fixers, and temporal discretization. The literature on dynam-
ical core development is vast, with origins that go back over
half a century. Consequently, the models discussed in this
paper only represent a sample of the many dynamical cores
that have been developed for general circulation modeling.
Some of the models that have not been discussed include
Fox-Rabinovitz et al. (1997), Prusa et al. (2008), Nair et al.
(2009), Baba et al. (2010), Donner et al. (2011), Ullrich and
Jablonowski (2012a), Gassmann (2013), Wood et al. (2014),
and Doyle (2014).

The vast diversity within the modern dynamical core ecosystem suggests that there is no consensus on a single approach that is intrinsically superior to other options. Choices made in the dynamical core confer advantages that include parallel scalability (Dennis et al., 2012), conservation of invariants (Thuburn, 2008), or representation of the kinetic energy spectrum (Skamarock, 2004). The repercussions that emerge from these choices can then be explored in the context of idealized test cases, such as the ones that have been proposed as part of DCMIP. The remaining papers in this series investigate how the models described in this paper are able to simulate three idealized test cases, each of which incorporates simplified model physics: a moist baroclinic wave, an idealized tropical cyclone, and a splitting supercell storm on a small planet. Where appropriate, metrics have been included that may be indicative of model performance. These tests can also be used for future dynamical core development to identify where a new dynamical core diverges from a suite of modern models.

*Code availability.* Information on the availability of source code for the models featured in this paper is tabulated below.

<table>
   <thead>
    <tr>
     <td>
      Short Name
     </td>
     <td>
      Code availability
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      ACME-A
     </td>
     <td>
      ACME, including ACME-A, is under active development funded by the US Department of Energy. ACME version 1.0 is scheduled to be publicly released under an open-source license in 2018 but is not available at present.†
     </td>
    </tr>
    <tr>
     <td>
      CSU
     </td>
     <td>
      CSU model source code is available under the Berkeley Software Distribution (BSD) three-clause license. The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.580099).
     </td>
    </tr>
   </tbody>
  </table>

<table>
   <thead>
    <tr>
     <td>
      Short Name
     </td>
     <td>
      Code availability
     </td>
    </tr>
   </thead>
   <tbody>
    <tr>
     <td>
      DYNAMICO
     </td>
     <td>
      DYNAMICO is open source and available online from IPSL Forge (http://forge.ipsl.jussieu.fr/dynamico) or directly by request to Thomas Dubos (dubos@lmd.polytechnique.fr). The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.583718).
     </td>
    </tr>
    <tr>
     <td>
      FV3
     </td>
     <td>
      FV3 model source code is available through the GFDL Virtual Lab (https://vlab.ncep.noaa.gov/web/fv3gfs). Access requires users to create a Virtual Lab account.
     </td>
    </tr>
    <tr>
     <td>
      FVM
     </td>
     <td>
      Model codes developed at ECMWF, including the IFS and FVM, are intellectual property of ECMWF and its member states. Although the FVM code is not publicly available at present, it is expected that FVM will be available in the near future under the OpenIFS license (http://www.ecmwf.int/en/research/projects/openifs). The repo tag for the version of FVM that applies for DCMIP is “v0.1”.a
     </td>
    </tr>
    <tr>
     <td>
      GEM
     </td>
     <td>
      Due to licensing requirements, GEM model code is only available by request to Abdessamad Qaddouri ( Abdessamad.qaddouri@canada.ca) or Vi-vian Lee ( vivian.lee2@canada.ca).
     </td>
    </tr>
    <tr>
     <td>
      ICON
     </td>
     <td>
      ICON is freely available to the scientific community for non-commercial research under an institutional license issued by project partners DWD and MPI-M. Because of the restrictions of this license, access to the code is only available by request to Günther Zängl ( guenther.zaengl@dwd.de) or Marco Giorgetta ( marco.giorgetta@mpimet.mpg.de).
     </td>
    </tr>
    <tr>
     <td>
      MPAS
     </td>
     <td>
      MPAS is an open-source model available under the BSD three-clause license via GitHub (https://github.com/MPAS-Dev/MPAS-Release). The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.583316).
     </td>
    </tr>
    <tr>
     <td>
      NICAM
     </td>
     <td>
      NICAM source code is available under the BSD two-clause license via Zenodo (http://dx.doi.org/10.5281/zenodo.580128). Further information on collaborating with the NICAM team can be found at http://nicam.jp/hiki/?Research+Collaborations.
     </td>
    </tr>
    <tr>
     <td>
      OLAM
     </td>
     <td>
      OLAM is open source and available online via SourceForge (https://sourceforge.net/projects/olam-model/). The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.582308).
     </td>
    </tr>
   </tbody>
  </table>
<table><thead><tr><td>Short name</td><td>Code availability</td></tr></thead><tbody><tr><td>Tempest</td><td>Tempest source code is available under the Lesser GNU Public License on GitHub (https://github.com/paullric/tempestmodel). The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.579649).</td></tr></tbody></table>

<footnote>ª In compliance with the GMD editorial requirements, this code has been made available to the topical editor in charge of the paper.</footnote>
## Appendix A: Moist non-hydrostatic equation sets

In this Appendix, we provide a detailed derivation of the fluid equations utilized by non-hydrostatic models. The physical constants which are used throughout this document are given in Table A1. The material derivative is used for quantities in the Lagrangian frame (following individual air parcels) and is given by

$$ \frac{d}{dt} = \frac{\partial}{\partial t} + \mathbf{u} \cdot \nabla, \qquad (A1) $$

where $\mathbf{u}$ denotes the 3-D vector velocity. Note that tracer variables $q_i$, including specific humidity $q$, in the absence of sources and sinks satisfy the simple Lagrangian relationship

$$ \frac{dq_i}{dt} = 0. \qquad (A2) $$

## A1 Diagnostic relationships

The atmospheric fluid is assumed to be an ideal gas. For moist air, the ideal gas constant $R^*$, specific heat capacity at constant pressure $c_p^*$, and specific heat capacity at constant volume $c_v^*$ are given by

$$ \begin{aligned} R^* &= R_d + (R_v - R_d)q, & c_p^* &= c_{pd} + (c_{pv} - c_{pd})q, \\ c_v^* &= c_{vd} + (c_{vv} - c_{vd})q. \end{aligned} \qquad (A3) $$

Note that in many models, $R^*$, $c_p^*$, and $c_v^*$ are approximated by $R_d$, $c_{pd}$, and $c_{vd}$, respectively. Dry air, water vapor, and moist air quantities all satisfy the linear relationship $R = c_p - c_v$. For a two-fluid system (dry air plus water vapor), two independent variables plus the specific humidity $q$ are needed to describe the thermodynamic state of the system. Key thermodynamic variables include dry air density $\rho_d$, moist density $\rho$, pressure $p$, vapor pressure $e$, temperature $T$, virtual temperature $T_v$, Exner pressure $\pi$, potential temperature $\theta$, and virtual potential temperature $\theta_v$. Common ratios $\kappa = R^*/c_p^*$, $\varepsilon = R_d/R_v$, and $\gamma = c_p^*/c_v^*$ are adopted here. Note that as additional water species are added (cloud water, rain water, etc.) additional independent variables are needed to capture the thermodynamic effects of these species, and the “virtual” quantities need to be modified accordingly, for instance, through the adoption of density potential temperature $\theta_p$.

Relationships between key thermodynamic variables arise from the ideal gas law, along with definitions of Exner pressure, potential temperature, and virtual potential tempera-

ture:

$$ \begin{aligned} p &= \rho R_d T_v, & \pi &= \left(\frac{p}{p_0}\right)^{\kappa}, & \theta &= T \left(\frac{p_0}{p}\right)^{\kappa}, \\ \theta_v &= T_v \left(\frac{p_0}{p}\right)^{\kappa}, \end{aligned} \qquad (A4) $$

which further give rise to

$$ \begin{aligned} p &= \left( \frac{\rho R_d \theta_v}{p_0^{\kappa}} \right)^{\gamma}, & \pi &= \left( \frac{\rho R_d \theta_v}{p_0} \right)^{R^*/c_v^*}, & \theta &= \frac{T}{\pi}, & \theta_v &= \frac{T_v}{\pi}. \end{aligned} \qquad (A5) $$

Note that virtual temperature is typically written as

$$ T_v = T \left( 1 + \frac{(1 - \varepsilon)}{\varepsilon} q \right), \qquad (A6) $$

which arises from the relationship

$$ T_v = \frac{T}{1 - \frac{e}{p}(1 - \varepsilon)}, \qquad (A7) $$

upon applying $e/p = q/\varepsilon$.

## A2 Prognostic equations for thermodynamic variables

Note that, as a consequence of Eq. (A2), the following simplifications can be applied:

$$ \frac{1}{T_v} \frac{dT_v}{dt} = \frac{1}{T} \frac{dT}{dt}, \frac{dR^*}{dt} = 0, \frac{dc_p^*}{dt} = 0, \frac{dc_v^*}{dt} = 0. \qquad (A8) $$

Mass conservation is typically represented through the continuity equation, which can be written in the Lagrangian frame as

$$ \frac{d\rho}{dt} = -\rho \nabla \cdot \mathbf{u}, \qquad (A9) $$

or equivalently in the Eulerian frame:

$$ \frac{\partial \rho}{\partial t} = -\nabla \cdot (\rho \mathbf{u}). \qquad (A10) $$

Further prognostic relationships can be derived from the thermodynamic equation, including the diabatic heating rate $J$:

$$ \frac{1}{T} \frac{dT}{dt} - \frac{\kappa}{p} \frac{dp}{dt} = \frac{J}{T c_p^*}, \qquad (A11) $$

which can be alternatively written as

$$ \frac{d\theta}{dt} = \frac{J\theta}{Tc_p^*}, \text{ or } \frac{d\theta_v}{dt} = \frac{J\theta_v}{Tc_p^*}. \qquad (A12) $$

These equations can then be combined with Eq. (A9) to obtain

$$ \frac{\partial}{\partial t}(\rho \theta_v) + \nabla \cdot (\rho \theta_v \mathbf{u}) = \frac{J \rho \theta_v}{T c_p^*}, \qquad (A13) $$
<table><caption>**Table A1.** A list of physical constants used in this document.</caption><thead><tr><th>Constant</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td>a<sub>ref</sub></td><td>Radius of the Earth</td><td>6.37122 × 10<sup>6</sup> m</td></tr><tr><td>Ω<sub>ref</sub></td><td>Rotational speed of the Earth</td><td>7.292 × 10<sup>-5</sup> s<sup>-1</sup></td></tr><tr><td>g<sub>c</sub></td><td>Gravitational acceleration</td><td>9.80616 m s<sup>-2</sup></td></tr><tr><td>p<sub>0</sub></td><td>Reference pressure</td><td>1000 hPa</td></tr><tr><td>c<sub>pd</sub></td><td>Specific heat capacity of dry air at constant pressure</td><td>1004.5 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>c<sub>pv</sub></td><td>Specific heat capacity of water vapor at constant pressure</td><td>1930.0 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>c<sub>vd</sub></td><td>Specific heat capacity of dry air at constant volume</td><td>717.5 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>c<sub>vv</sub></td><td>Specific heat capacity of water vapor at constant volume</td><td>1460.0 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>R<sub>d</sub></td><td>Gas constant for dry air</td><td>287.0 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>R<sub>v</sub></td><td>Gas constant for water vapor</td><td>461.5 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>ε</td><td>Ratio of R<sub>d</sub> to R<sub>v</sub></td><td>0.622</td></tr><tr><td>M<sub>v</sub></td><td>Constant for virtual temperature conversion</td><td>0.608</td></tr><tr><td>ρ<sub>water</sub></td><td>Reference density of water</td><td>1000 kg m<sup>-3</sup></td></tr></tbody></table>

or similarly for $\theta$. In conjunction with the material derivative of the ideal gas law,

$$ \frac{1}{p} \frac{dp}{dt} = \frac{1}{\rho} \frac{d\rho}{dt} + \frac{1}{T_v} \frac{dT_v}{dt}, \qquad (A14) $$

the thermodynamic equation can be written in the form

$$ \frac{c_v^*}{R^* T_v} \frac{dT_v}{dt} - \frac{1}{\rho} \frac{d\rho}{dt} = \frac{J}{T R^*}. \qquad (A15) $$

Then, substituting Eq. (A9) gives a prognostic equation for virtual temperature:

$$ \frac{c_v^*}{R^*} \frac{dT_v}{dt} + T_v \nabla \cdot \mathbf{u} = \frac{J T_v}{T R^*}. \qquad (A16) $$

The prognostic equation for temperature is identical except with $T$ substituted for $T_v$. An analogous equations for pressure can be obtained through a similar procedure:

$$ \frac{c_v^*}{c_p^*} \frac{dp}{dt} + p \nabla \cdot \mathbf{u} = \frac{Jp}{T c_p^*}, \qquad (A17) $$

and similarly for Exner pressure:

$$ \frac{c_v^*}{R^*} \frac{d\pi}{dt} + \pi \nabla \cdot \mathbf{u} = \frac{J\pi}{T c_p^*}. \qquad (A18) $$

## A3 Momentum equations

In coordinate-invariant form, the prognostic velocity equations may be written in either the Lagrangian or Eulerian frame as

$$ \frac{du}{dt} = \frac{\partial u}{\partial t} + u \cdot \nabla u = -\frac{1}{\rho} \nabla p - 2\Omega \times u - \nabla \Phi, \qquad (A19) $$

where $\Omega$ denotes the planetary vorticity vector and $\Phi$ is the geopotential function. The three terms on the right-hand side

of this expression correspond to pressure gradient, Coriolis, and gravitational force, respectively. In Eulerian form, one must be careful with the treatment of the momentum advection term $\mathbf{u} \cdot \nabla \mathbf{u}$, since in an arbitrary coordinate frame this term will give rise to Christoffel symbols associated with derivatives of the vector basis. Note that it is common to rewrite the pressure-gradient force using the relationship

$$ -\frac{1}{\rho} \nabla p = -c_p^* \theta \left[ \nabla \pi - \pi \ln \left( \frac{p}{p_0} \right) \nabla \kappa \right], \qquad (A20) $$

which follows from Eq. (A4). Note that often in non-hydrostatic models, $\kappa$ is assumed constant and the $\nabla \kappa$ term neglected. A second form of Eq. (A19) emerges on substituting the vector calculus identity

$$ \mathbf{u} \cdot \nabla \mathbf{u} = \nabla K + \zeta \times \mathbf{u}, \qquad (A21) $$

where $K = \frac{1}{2}(\mathbf{u} \cdot \mathbf{u})$ is the 3-D specific kinetic energy and $\zeta = \nabla \times \mathbf{u}$ is the 3-D relative vorticity vector. This gives rise to the 3-D vector-invariant form,

$$ \frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho} \nabla p - \nabla (K + \Phi) - (\zeta + 2\Omega) \times \mathbf{u}. \qquad (A22) $$

Because no gradients of vectors appear in this equation, it avoids derivatives of the coordinate basis that would arise from the momentum transport term $\mathbf{u} \cdot \nabla \mathbf{u}$ in Eq. (A19). In conjunction with Eq. (A9), Eq. (A19) also gives rise to the flux-form momentum equations,

$$ \frac{\partial}{\partial t}(\rho \mathbf{u}) = -\nabla \cdot (\mathbf{u} \otimes \mathbf{u} + \mathcal{I} p) - 2\Omega \times (\rho \mathbf{u}) - \rho \nabla \Phi, \qquad (A23) $$

where $\mathbf{u} \otimes \mathbf{u}$ denotes the outer product and $\mathcal{I}$ is the identity matrix.

The equations above still provide some flexibility with regard to the choice of $\Phi$ and $\Omega$. For deep-atmosphere models, one typically chooses

$$ \Phi = g_c a^2 \left[ \frac{1}{a} - \frac{1}{a+z} \right], \quad \text{and} \quad \Omega = \Omega(k \sin \varphi + j \cos \varphi), $$
$$ (A24) \qquad \text{which, from } \mathbf{u}_h = \mathbf{u} - \mathbf{w}\mathbf{k}, \text{ then gives rise to} $$

where $g_c$ is gravitational acceleration at the surface, $a$ is the radius of the planet, $\Omega$ is the rotation rate (in s⁻¹), $\varphi$ is the latitude, $j$ is the unit vector oriented in the meridional direction, and $k$ is the unit vector oriented in the vertical direction. For models that do not utilize a height-based vertical coordinate, the geopotential is generally treated as a prognostic variable, with an evolution equation that emerges from the definition $w = dz/dt$,

$$ \frac{d\Phi}{dt} = \frac{a^2 g_c w}{(a+z)^2}. \qquad (A25) $$

For shallow-atmosphere models, the geopotential takes the simpler form

$$ \Phi = g_c z, \quad \text{and} \quad \Omega = \Omega \sin \varphi k, \qquad (A26) $$

where $z$ is the altitude above the surface. In this case, we write $2\Omega = fk$, where $f = 2\Omega \sin\varphi$ is the Coriolis parameter. The evolution equation for the shallow-atmosphere geopotential is then

$$ \frac{d\Phi}{dt} = g_c w. \qquad (A27) $$

## A4 Orthogonal formulation

Under the orthogonal formulation, projection of a vector field $b$ onto its horizontal components is defined via

$$ [b]_z = b - (b \cdot k)k. \qquad (A28) $$

When applied to the velocity vector, this gives rise to the decomposition

$$ \mathbf{u} = \mathbf{u}_h + \mathbf{w}\mathbf{k}, \qquad (A29) $$

where $k = \nabla z$ is the unit vector in the vertical direction and $\mathbf{u}_h = [\mathbf{u}]_z$ ($\mathbf{u}_h$ is aligned with surfaces of constant $z$). In the orthogonal formulation, the material derivative expands as

$$ \frac{d}{dt} = \frac{\partial}{\partial t} + \mathbf{u}_h \cdot \nabla + \mathbf{w}\mathbf{k} \cdot \nabla. \qquad (A30) $$

For the special case of the material derivative applied to scalars, this equation can also be written as

$$ \frac{d}{dt} = \frac{\partial}{\partial t} + \mathbf{u}_h \cdot \nabla_z + \mathbf{w} \frac{\partial}{\partial z}. \qquad (A31) $$

where $\nabla_z b = [\nabla b]_z$ denotes the gradient along constant $z$ surfaces. From here, the vector-invariant form velocity equation obtained by multiplying Eq. (A22) by $k$ expands as

$$ \frac{\partial w}{\partial t} = -\frac{1}{\rho} \frac{\partial p}{\partial z} - \frac{\partial}{\partial z} (K + \Phi) - [(\zeta + 2\Omega) \times \mathbf{u}] \cdot \mathbf{k}, \qquad (A32) $$

$$ \frac{\partial \mathbf{u}_h}{\partial t} = -\frac{1}{\rho} \nabla_z p - \nabla_z (K + \Phi) - [(\zeta + 2\Omega) \times \mathbf{u}]_z. \qquad (A33) $$

Due to its association with hydrostatic models, it is common to use the 2-D kinetic energy, $K_2 = \frac{1}{2}(\mathbf{u}_h \cdot \mathbf{u}_h)$. Decomposing the momentum transport term into horizontal and vertical components gives

$$ \mathbf{u} \cdot \nabla \mathbf{u} = \mathbf{u}_h \cdot \nabla \mathbf{u}_h + (\nabla \times \mathbf{u}_h) \times (\mathbf{w}\mathbf{k}) + (\mathbf{u} \cdot \nabla \mathbf{w})\mathbf{k}. \qquad (A34) $$

The first term in this expression admits the relationships

$$ [\mathbf{u}_h \cdot \nabla \mathbf{u}_h]_z = \nabla_z K_2 + \zeta_h \mathbf{k} \times \mathbf{u}_h, \qquad (A35) \\ (\mathbf{u}_h \cdot \nabla \mathbf{u}_h) \cdot \mathbf{k} = -\mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k}) $$

$$ = -K_2(\nabla \cdot \mathbf{k}) - \frac{1}{2}[\mathbf{u}_h \cdot (\nabla \times \mathbf{u}_t) + (\nabla \times \mathbf{u}_h) \cdot \mathbf{u}_t], \qquad (A36) $$

where $\zeta_h = (\nabla \times \mathbf{u}) \cdot \mathbf{k} = (\nabla \times \mathbf{u}_h) \cdot \mathbf{k}$ is the relative vorticity scalar and $\mathbf{u}_t = \mathbf{k} \times \mathbf{u}_h$. Note that this equation does incorporate metric terms associated with horizontal advection of $k$, which must be accounted for.

Thus, the vertical velocity equation, obtained by taking Eq. (A19)·$k$, is

$$ \frac{\partial w}{\partial t} = \mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k}) - w \frac{\partial w}{\partial z} - \mathbf{u}_h \cdot \nabla_z w - \frac{\partial \Phi}{\partial z} - \frac{1}{\rho} \frac{\partial p}{\partial z} - (2\Omega \times \mathbf{u}_h) \cdot \mathbf{k}. \qquad (A37) $$

Then, subtracting Eq. (A37)·$k$ from Eq. (A19) gives

$$ \frac{\partial \mathbf{u}_h}{\partial t} = -w(\mathbf{u}_h \cdot \nabla \mathbf{k}) - w \frac{\partial \mathbf{u}_h}{\partial z} - \nabla_z (K_2 + \Phi) - \frac{1}{\rho} \nabla_z p - \zeta_h \mathbf{k} \times \mathbf{u}_h - [2\Omega \times \mathbf{u}]_z. \qquad (A38) $$

Note that, under the shallow-atmosphere approximation, the metric term $\mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k})$ in Eq. (A37) is set equal to zero in accordance with Phillips (1966).

## A5 Arbitrary vertical coordinates

The dynamical equations are now formulated in terms of the vertical coordinate $s(t, x, z)$ with $\partial s/\partial z \neq 0$ everywhere, i.e., following Kasahara (1974) (hereafter K74). Since $x$ and $t$ are shared between the two coordinate systems, the chain rule can be applied to obtain expressions

$$ \frac{\partial}{\partial z} = \frac{\partial s}{\partial z} \frac{\partial}{\partial s}, \qquad \nabla_s = \nabla_z + (\nabla_s z)(\mathbf{k} \cdot \nabla), \\ \left(\frac{\partial}{\partial t}\right)_s = \frac{\partial}{\partial t} + \left(\frac{\partial z}{\partial t}\right)_s (\mathbf{k} \cdot \nabla), \qquad (A39) $$

which correspond to derivatives in the vertical, in the horizontal, and in time. This final expression is used to describe
the rate of change of a quantity on s surfaces. These operators then yield the useful identities

$$ \frac{\partial s}{\partial z} = \left(\frac{\partial z}{\partial s}\right)^{-1} \nabla_z s = -\left(\frac{\partial s}{\partial z}\right) \nabla_s z, \quad \frac{\partial s}{\partial t} = -\left(\frac{\partial s}{\partial z}\right) \left(\frac{\partial z}{\partial t}\right)_s. \qquad (A40) $$

From here, Eq. (A39) also gives rise to

$$ \nabla_z = \nabla_s - \frac{\partial s}{\partial z} (\nabla_s z) \frac{\partial}{\partial s}, \qquad (A41) $$

which can be used directly to rewrite Eq. (A33) or Eq. (A38) in terms of derivatives over s. Note that the operators $\nabla_z$ and $\nabla_s$ are usually introduced in the context of 2-D flows; however, the construction described here has the advantage of working seamlessly in a 3-D context, while admitting the properties $k \cdot \nabla_z A = 0$ and $k \cdot \nabla_s A = 0$ for any scalar field A.

From Eq. (A39), it can be shown that the 2-D divergence on s surfaces (given by K74 Eq. 3.17) is

$$ \nabla_s \cdot \mathbf{u}_h = \nabla_z \cdot \mathbf{u}_h + \left(\frac{\partial s}{\partial z}\right) (\nabla_s z) \cdot \left(\frac{\partial \mathbf{u}_h}{\partial s}\right), \qquad (A42) $$

and that the 2-D curl is given by

$$ \nabla_s \times \mathbf{u}_h = \nabla_z \times \mathbf{u}_h + \left(\frac{\partial s}{\partial z}\right) (\nabla_s z) \times \left(\frac{\partial \mathbf{u}_h}{\partial s}\right), \qquad (A43) $$

where $\nabla_z \times \mathbf{u}_h = k(k \cdot (\nabla \times \mathbf{u}_h))$. Notably, these expressions are valid for both shallow- and deep-atmosphere formulations.

The generalized velocity $\dot{s}$ following a fluid parcel is defined by

$$ \dot{s} \equiv \frac{ds}{dt} = \frac{\partial s}{\partial t} + \mathbf{u} \cdot \nabla s = \mathbf{u}_h \cdot \nabla_z s + \left[ w - \left(\frac{\partial z}{\partial t}\right)_s \right] \frac{\partial s}{\partial z}. \qquad (A44) $$

Then, using Eqs. (A39) and (A44) to rewrite Eq. (A31) gives an expression for the material derivative for scalars on s surfaces:

$$ \frac{dA}{dt} = \left(\frac{\partial A}{\partial t}\right)_s - \frac{\partial s}{\partial z} \left(\frac{\partial z}{\partial t}\right)_s \frac{\partial A}{\partial s} + \mathbf{u}_h \cdot \left[ \nabla_s A + (\nabla_z s) \frac{\partial A}{\partial s} \right] + w \frac{\partial s}{\partial z} \frac{\partial A}{\partial s} \qquad (A45) $$

$$ = \left(\frac{\partial A}{\partial t}\right)_s + \mathbf{u}_h \cdot \nabla_s A + \dot{s} \frac{\partial A}{\partial s}. \qquad (A46) $$

A similar expression arises for vectors, although in this case $\mathbf{u}_h \cdot \nabla \mathbf{a} \neq \mathbf{u}_h \cdot \nabla_z \mathbf{a}$ implies we cannot use the operator $\nabla_s$ in the form (A39), and instead obtain

$$ \frac{d\mathbf{a}}{dt} = \left(\frac{\partial \mathbf{a}}{\partial t}\right)_s + \left[\mathbf{u}_h \cdot \nabla \mathbf{a} + (\mathbf{u}_h \cdot \nabla_s z)(\mathbf{k} \cdot \nabla \mathbf{a})\right] + \dot{s} \frac{\partial \mathbf{a}}{\partial s}. \qquad (A47) $$

## A6 Conservation laws in arbitrary vertical coordinates

Using Eq. (A42), we observe that the 3-D divergence on the sphere takes the form

$$ \nabla \cdot \mathbf{u} = \nabla_z \cdot \mathbf{u}_h + \frac{1}{\alpha} \frac{\partial}{\partial z}(\alpha w), \qquad (A48) $$

where $\alpha = 1$ for shallow-atmosphere models and $\alpha = r^2 = (a+z)^2$ for deep-atmosphere models. Using $w = dz/dt$, this last term also takes the form

$$ \frac{1}{\alpha} \frac{\partial}{\partial z}(\alpha w) = \frac{\partial w}{\partial z} + \frac{w}{\alpha} \frac{\partial \alpha}{\partial z} = \frac{\partial w}{\partial z} + \frac{1}{\alpha} \frac{d\alpha}{dt}. \qquad (A49) $$

Using Eq. (A40) to rewrite Eq. (A44) gives rise to

$$ w = \left(\frac{\partial z}{\partial t}\right)_s + \mathbf{u}_h \cdot \nabla_s z + \dot{s} \left(\frac{\partial s}{\partial z}\right)^{-1}, \qquad (A50) $$

which is then differentiated to yield K74 Eq. (3.16),

$$ \frac{\partial w}{\partial z} = \left(\frac{\partial s}{\partial z}\right) \left[ \frac{d}{dt} \left(\frac{\partial s}{\partial z}\right)^{-1} + \left(\frac{\partial \mathbf{u}_h}{\partial s}\right) \cdot (\nabla_s z) \right] + \frac{\partial \dot{s}}{\partial s} = 0. \qquad (A51) $$

Substituting this expression into continuity Eq. (A9) and using Eqs. (A48), (A49), and (A51) then leads to

$$ \frac{d}{dt} \left[ \alpha \left( \frac{\partial s}{\partial z} \right)^{-1} \rho \right] + \alpha \left( \frac{\partial s}{\partial z} \right)^{-1} \rho \left[ \nabla_z \cdot \mathbf{u}_h + \left( \frac{\partial s}{\partial z} \right) \left( \frac{\partial \mathbf{u}_h}{\partial s} \right) \cdot (\nabla_s z) \right] + \alpha \left( \frac{\partial s}{\partial z} \right)^{-1} \rho \frac{\partial \dot{s}}{\partial s} = 0. \qquad (A52) $$

Defining the pseudo-density as

$$ \rho_s = \alpha \left(\frac{\partial s}{\partial z}\right)^{-1} \rho, \qquad (A53) $$

and using Eq. (A46) in the form

$$ \frac{d\rho_s}{dt} = \left(\frac{\partial \rho_s}{\partial t}\right)_s + \mathbf{u}_h \cdot \nabla_s \rho_s + \dot{s} \frac{\partial \rho_s}{\partial s}, \qquad (A54) $$

along with Eq. (A42), leads to

$$ \left(\frac{\partial \rho_s}{\partial t}\right)_s + \nabla_s \cdot (\rho_s \mathbf{u}_h) + \frac{\partial}{\partial s}(\rho_s \dot{s}) = 0. \qquad (A55) $$

Hence, for any quantity that is conserved following a fluid parcel (i.e., $dq/dt = 0$),

$$ \left(\frac{\partial \rho_s q}{\partial t}\right)_s + \nabla_s \cdot (\rho_s q \mathbf{u}_h) + \frac{\partial}{\partial s}(\rho_s q \dot{s}) = 0. \qquad (A56) $$

In particular, the prognostic equation for virtual potential temperature (or equivalently for potential temperature) reads

$$ \left(\frac{\partial \rho_s \theta_v}{\partial t}\right)_s + \nabla_s \cdot (\rho_s \theta_v \mathbf{u}_h) + \frac{\partial}{\partial s}(\rho_s \theta_v \dot{s}) = \frac{J \rho_s \theta_v}{c_p^*}. \qquad (A57) $$
## A7 2-D vector-invariant form

The prognostic equations utilizing horizontal kinetic energy $K_2$ in place of $K$ are derived by applying Eq. (A39) to Eq. (A37), yielding

$$ \left(\frac{\partial w}{\partial t}\right)_s = \mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k}) - \mathbf{u}_h \cdot \nabla w + \left(\frac{\partial s}{\partial z}\right) \left\{ \left[ \left(\frac{\partial z}{\partial t}\right)_s - w \right] \frac{\partial w}{\partial s} - \frac{\partial \Phi}{\partial s} - \frac{1}{\rho} \frac{\partial p}{\partial s} \right\} - (2\Omega \times \mathbf{u}_h) \cdot \mathbf{k}. \quad (A58) $$

Similarly, from Eq. (A38),

$$ \left(\frac{\partial \mathbf{u}_h}{\partial t}\right)_s = -w(\mathbf{u}_h \cdot \nabla \mathbf{k}) + \left[\left(\frac{\partial z}{\partial t}\right)_s - w\right] \left(\frac{\partial s}{\partial z}\right) \frac{\partial \mathbf{u}_h}{\partial s} - \zeta_h \mathbf{k} \times \mathbf{u}_h - \nabla_z (K_2 + \Phi) - \frac{1}{\rho} \nabla_z p - [2\Omega \times \mathbf{u}]_z. \quad (A59) $$

Observe that both of these equations simplify when $w = (\partial z/\partial t)_s$, i.e., model levels are advected with the vertical wind.

An alternative form of these equation can similarly be obtained in terms of $\dot{s}$. Substituting Eq. (A44) into Eq. (A58) then gives

$$ \left(\frac{\partial w}{\partial t}\right)_s = \mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k}) - \mathbf{u}_h \cdot \nabla_s w - \dot{s} \frac{\partial w}{\partial s} + \left(\frac{\partial s}{\partial z}\right) \left[ -\frac{\partial \Phi}{\partial s} - \frac{1}{\rho} \frac{\partial p}{\partial s} \right] - (2\Omega \times \mathbf{u}_h) \cdot \mathbf{k}. \quad (A60) $$

Similarly, substituting Eq. (A44) into Eq. (A59) and using the identity

$$ (\mathbf{u}_h \cdot \nabla_s z) \left(\frac{\partial s}{\partial z}\right) \frac{\partial \mathbf{u}_h}{\partial s} = (\nabla_s \times \mathbf{u}_h) \times \mathbf{u}_h - \zeta_h \mathbf{k} \times \mathbf{u}_h + (\nabla_s z) \frac{\partial K_2}{\partial z} \quad (A61) $$

then gives

$$ \left(\frac{\partial \mathbf{u}_h}{\partial t}\right)_s = -w(\mathbf{u}_h \cdot \nabla \mathbf{k}) - \nabla_s K_2 - \zeta_s \mathbf{k} \times \mathbf{u}_h - \dot{s} \frac{\partial \mathbf{u}_h}{\partial s} - \nabla_z \Phi - \frac{1}{\rho} \nabla_z p - [2\Omega \times \mathbf{u}]_z, \quad (A62) $$

where

$$ \nabla_s \times \mathbf{u}_h = \mathbf{k} \zeta_s, \quad \text{and} \quad \zeta_s = \mathbf{k} \cdot (\nabla_s \times \mathbf{u}_h). \quad (A63) $$

In this case, the vertical advection terms are removed when $\dot{s} = 0$, i.e., the vertical coordinate is advected with the 3-D wind $\mathbf{u}$.

Note that under the shallow-atmosphere approximation, the first metric terms (those that include $(\mathbf{u}_h \cdot \nabla \mathbf{k})$) in Eqs. (A58)–(A62) are typically dropped.

## A8 3-D vector-invariant form

From Eqs. (A22) and (A39), the evolution equation for the 3-D velocity vector takes the form

$$ \left(\frac{\partial \mathbf{u}}{\partial t}\right)_s = \left(\frac{\partial z}{\partial t}\right)_s (\mathbf{k} \cdot \nabla \mathbf{u}) - \nabla(K + \Phi) - \frac{1}{\rho} \nabla p - (\zeta + 2\Omega) \times \mathbf{u}. \quad (A64) $$

Then, taking the dot product of this expression with $\mathbf{k}$ gives

$$ \left(\frac{\partial w}{\partial t}\right)_s = \left(\frac{\partial s}{\partial z}\right) \left[ \left(\frac{\partial z}{\partial t}\right)_s \frac{\partial w}{\partial s} - \frac{\partial}{\partial s}(K + \Phi) - \frac{1}{\rho} \frac{\partial p}{\partial s} \right] - \left[ (\zeta + 2\Omega) \times \mathbf{u} \right] \cdot \mathbf{k}, \quad (A65) $$

where we have used $\mathbf{k} \cdot (\mathbf{k} \cdot \nabla \mathbf{u}) = \mathbf{k} \cdot \nabla w$. Similarly, the prognostic equation for horizontal velocity from Eq. (A33) is reformulated as

$$ \left(\frac{\partial \mathbf{u}_h}{\partial t}\right)_s = \left(\frac{\partial z}{\partial t}\right)_s \left(\frac{\partial s}{\partial z}\right) \frac{\partial \mathbf{u}_h}{\partial s} - \nabla_z (K + \Phi) - \frac{1}{\rho} \nabla_z p - \left[ (\zeta + 2\Omega) \times \mathbf{u} \right]_z. \quad (A66) $$

Note that the vorticity term in this expression can be simplified further using

$$ [(\zeta + 2\Omega) \times \mathbf{u}]_z = - (\zeta_h + \mathbf{k} \cdot 2\Omega)(\mathbf{u}_h \times \mathbf{k}) - w \mathbf{k} \times (\zeta + 2\Omega), \quad (A67) $$

and

$$ -\mathbf{k} \times \zeta = \mathbf{k} \cdot \nabla \mathbf{u} - \nabla(\mathbf{k} \cdot \mathbf{u}) + \mathbf{u} \cdot \nabla \mathbf{k} = \frac{\partial \mathbf{u}_h}{\partial z} - \nabla_z w + \mathbf{u}_h \cdot \nabla \mathbf{k}. \quad (A68) $$

## A9 Covariant component formulation

In conjunction with Eq. (A41), the horizontal momentum equation (in 2-D vector-invariant form as Eqs. A59 or A62, or in 3-D vector-invariant form as Eq. A66) with an arbitrary vertical coordinate gives rise to a two-term pressure gradient. This can be avoided by prognosing the covariant components of the velocity in place of the physical velocity components. We define a horizontal covariance operator by

$$ [\mathbf{b}]_s = [\mathbf{b}]_z + (\nabla_s z)(\mathbf{k} \cdot \mathbf{b}). \quad (A69) $$

Applying this operator to the horizontal velocity gives

$$ \mathbf{v}_h = [\mathbf{u}]_s = \mathbf{u}_h + (\nabla_s z)w. \quad (A70) $$

For a time-dependent $s$ coordinate, we obtain the identity

$$ \left[ \frac{\partial}{\partial t} (\nabla_s z) \right]_s = \nabla_s \left( \frac{\partial z}{\partial t} \right)_s, \quad (A71) $$
and so can write

$$ \left(\frac{\partial \mathbf{v}_h}{\partial t}\right)_s = \left(\frac{\partial \mathbf{u}_h}{\partial t}\right)_s + (\nabla_s z)\left(\frac{\partial w}{\partial t}\right)_s + w \nabla_s \left(\frac{\partial z}{\partial t}\right)_s . \quad (\text{A72}) $$

Then, using Eqs. (A72), (A65), and (A66) and identity

$$ \begin{aligned} \left(\frac{\partial z}{\partial t}\right)_s \left(\frac{\partial s}{\partial z}\right) \left[\frac{\partial \mathbf{u}_h}{\partial s} + (\nabla_s z) \frac{\partial w}{\partial s}\right] + w \nabla_s \left(\frac{\partial z}{\partial t}\right)_s & \\ = \left(\frac{\partial s}{\partial z}\right) \left(\frac{\partial z}{\partial t}\right)_s \left\{\frac{\partial \mathbf{v}_h}{\partial s} - \nabla_s \left[\left(\frac{\partial s}{\partial z}\right)^{-1} w\right]\right\} & \\ + \nabla_s \left[\left(\frac{\partial z}{\partial t}\right)_s w\right], & \end{aligned} \quad (\text{A73}) $$

gives

$$ \begin{aligned} \left(\frac{\partial \mathbf{v}_h}{\partial t}\right)_s = & -\nabla_s \left[ K - w \left(\frac{\partial z}{\partial t}\right)_s + \Phi \right] - \frac{1}{\rho} \nabla_s p \\ & - \left[ (\zeta + 2\Omega) \times \mathbf{u} \right]_s \\ & + \left(\frac{\partial s}{\partial z}\right) \left(\frac{\partial z}{\partial t}\right)_s \left\{ \frac{\partial \mathbf{v}_h}{\partial s} - \nabla_s \left[ \left(\frac{\partial s}{\partial z}\right)^{-1} w \right] \right\}. \end{aligned} \quad (\text{A74}) $$

Finally, we can expand the vorticity term and hence obtain

$$ \begin{aligned} \left(\frac{\partial \mathbf{v}_h}{\partial t}\right)_s = & -\nabla_s \left[ K - w \left(\frac{\partial z}{\partial t}\right)_s + \Phi \right] - \frac{1}{\rho} \nabla_s p - [2\Omega \times \mathbf{u}]_s \\ & + [\mathbf{k} \cdot \nabla_s \times \mathbf{v}_h] (\mathbf{u}_h \times \mathbf{k}) \\ & - \dot{s} \left\{ \frac{\partial \mathbf{v}_h}{\partial s} - \nabla_s \left[ \left(\frac{\partial s}{\partial z}\right)^{-1} w \right] \right\} - \left(\frac{\partial z}{\partial t}\right)_s \mathbf{u}_h \cdot \nabla \mathbf{k}. \end{aligned} \quad (\text{A76}) $$

## A10 Vorticity divergence form

The vorticity divergence form of the dynamical equations in an arbitrary vertical coordinate predicts the absolute vorticity ($\zeta_h^*$) and velocity divergence ($D$) given by

$$ \zeta_h^* = (\nabla_s \times \mathbf{u}_h + 2\Omega) \cdot \mathbf{k}, \quad (\text{A77}) $$

and

$$ D = \nabla_s \cdot \mathbf{u}_h, \quad (\text{A78}) $$

respectively, instead of the horizontal velocity. The horizontal velocity can be obtained from the stream function $\psi$ and the velocity potential $\chi$ following

$$ \mathbf{u}_h = \mathbf{k} \times \nabla_s \psi + \nabla_s \chi. \quad (\text{A79}) $$

By using Eq. (A79) in Eqs. (A77) and (A78), we obtain the elliptic equations that diagnose the stream function and velocity potential from the predicted velocity and divergence as

$$ \nabla_s^2 = \zeta_h^* - 2\Omega \cdot \mathbf{k}, \quad \text{and} \quad \nabla_s^2 \chi = D, \quad (\text{A80}) $$

respectively.

By taking the material derivative (Eq. A46) of Eq. (A77) and using horizontal momentum Eqs. (A38), (A79), and (A80), the absolute vorticity prediction equation emerges:

$$ \begin{aligned} \left(\frac{\partial \zeta_h^*}{\partial t}\right)_s - J_s(\zeta_h^*, \psi) + \nabla_s \cdot (\zeta_h^* \nabla_s \chi) + \nabla_s \cdot \left(\dot{s} \frac{\partial}{\partial s} \nabla_s \psi\right) & \\ + \mathbf{k} \cdot \nabla_s \times \left(\dot{s} \frac{\partial}{\partial s} \nabla_s \chi\right) + J_s(\rho^{-1}, p) = 0, & \end{aligned} \quad (\text{A81}) $$

where $J_s(a, b) = \mathbf{k} \cdot \nabla_s \times (a \nabla_s b)$ is the Jacobian operator. It can also be shown that $\dot{s}$ relates to the vertical velocity $w$ through

$$ \dot{s} = \left(\frac{\partial s}{\partial z}\right) (w - w_c), \quad (\text{A82}) $$

where

$$ w_c = \left(\frac{\partial z}{\partial t}\right)_s + (\mathbf{k} \times \nabla_s \psi + \nabla_s \chi) \cdot (\nabla_s z). \quad (\text{A83}) $$

By taking the material derivative of Eq. (A78) and using Eqs. (A38), (A79), and (A80), we can obtain the divergence prediction equation

$$ \begin{aligned} \left(\frac{\partial D}{\partial t}\right)_s - J_s(\zeta_h^*, \chi) - \nabla_s \cdot (\zeta_h^* \nabla_s \psi) + \nabla_s \cdot \left(\dot{s} \frac{\partial}{\partial s} \nabla_s \chi\right) & \\ + \left(\mathbf{k} \times \frac{\partial}{\partial s} \nabla_s \psi\right) \cdot \nabla_s \dot{s} & \\ + \nabla_s \cdot (\nabla_s K_2 + g \nabla_s z) + \nabla_s \cdot \left(\frac{1}{\rho} \nabla_s p\right) = 0, & \end{aligned} \quad (\text{A84}) $$

where $K_2$ can be reformulated in terms of stream function and velocity potential as

$$ K_2 = \frac{1}{2} \left[ \nabla_s \cdot (\psi \nabla_s \psi) - \psi \nabla_s^2 \psi + \nabla_s \cdot (\chi \nabla_s \chi) - \chi \nabla_s^2 \chi \right] + J_s(\psi, \chi). \quad (\text{A86}) $$

## A11 Momentum form

The momentum form of the prognostic equations emerges by combining the prognostic velocity equations with a continuity equation. Essentially, any of the continuity equations can be chosen, as long as the mass field represented by the equation is everywhere non-zero. However, the most common options are moist pseudo-density (Ullrich and Jablonowski, 2012a) or dry pseudo-density (Skamarock et al., 2012). Here, we denote our density variable by $\tilde{\rho}_s$ and assume no external sources or sinks of $\tilde{\rho}$. Multiplying Eq. (A60) by $\tilde{\rho}_s$ and using
Eq. (A55) gives

$$ \left( \frac{\partial \tilde{\rho}_s w}{\partial t} \right)_s = \tilde{\rho}_s \mathbf{u}_h \cdot (\mathbf{u}_h \cdot \nabla \mathbf{k}) - \nabla_s \cdot (\tilde{\rho}_s \mathbf{u}_h w) - \frac{\partial}{\partial s} (\tilde{\rho}_s \dot{s} w) \\ + \tilde{\rho}_s \left( \frac{\partial s}{\partial z} \right) \left[ -\frac{\partial \Phi}{\partial s} - \frac{1}{\rho} \frac{\partial p}{\partial s} \right] - (2\Omega \times \tilde{\rho}_s \mathbf{u}_h) \cdot \mathbf{k}. \quad (A87) $$

Similarly, from Eq. (A62), we have

$$ \left( \frac{\partial \tilde{\rho}_s \mathbf{u}_h}{\partial t} \right)_s = - \tilde{\rho}_s w (\mathbf{u}_h \cdot \nabla \mathbf{k}) - \tilde{\rho}_s \nabla_s K_2 - \zeta_s \mathbf{k} \times \tilde{\rho}_s \mathbf{u}_h \\ - \mathbf{u}_h \nabla_s \cdot (\tilde{\rho}_s \mathbf{u}_h) - \frac{\partial}{\partial s} (\tilde{\rho}_s \dot{s} \mathbf{u}_h) \\ - \tilde{\rho}_s \left( \nabla_z \Phi + \frac{1}{\rho} \nabla_z p \right) - [2\Omega \times \tilde{\rho}_s \mathbf{u}]_z. \quad (A88) $$
*Author contributions.* Text in this paper describing individual models was provided by the respective modeling teams. Final composition and development of Appendix A was performed by PAU.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* DCMIP2016 is sponsored by the National Center for Atmospheric Research Computational Information Systems Laboratory, the Department of Energy Office of Science (award no. DE-SC0016015), the National Science Foundation (award no. 1629819), the National Aeronautics and Space Administration (award no. NNX16AK51G), the National Oceanic and Atmospheric Administration Great Lakes Environmental Research Laboratory (award no. NA12OAR4320071), the Office of Naval Research, and CU Boulder Research Computing. This work was made possible with support from our student and postdoctoral participants: Sabina Abba Omar, Scott Bachman, Amanda Back, Tobias Bauer, Vinicius Capistrano, Spencer Clark, Ross Dixon, Christopher Eldred, Robert Fajber, Jared Ferguson, Emily Foshee, Ariane Frassoni, Alexander Goldstein, Jorge Guerra, Chasity Henson, Adam Herrington, Tsung-Lin Hsieh, Dave Lee, Theodore Letcher, Weiwei Li, Laura Mazzaro, Maximo Menchaca, Jonathan Meyer, Farshid Nazari, John O'Brien, Bjarke Tobias Olsen, Hossein Parishani, Charles Pelletier, Thomas Rackow, Kabir Rasouli, Cameron Rencurrel, Koichi Sakaguchi, Gökhan Sever, James Shaw, Konrad Simon, Abhishek Srivastava, Nicholas Szapiro, Kazushi Takemura, Pushp Raj Tiwari, Chii-Yun Tsai, Richard Urata, Karin van der Wiel, Lei Wang, Eric Wolf, Zheng Wu, Haiyang Yu, Sungduk Yu, and Jiawei Zhuang. We would also like to thank Rich Loft, Cecilia Banger, Kathryn Peczkowicz, and Rory Kelly (NCAR); Perla Dinger, Carmen Ho, and Gina Skyberg (UC Davis); and Kristi Hansen (University of Michigan) for administrative support during the workshop and summer school.

Edited by: Julia Hargreaves

Reviewed by: Hilary Weller and Thomas Melvin

## References

Arakawa, A. and Konor, C. S.: Unification of the anelastic and quasi-hydrostatic systems of equations, Mon. Weather Rev., 137, 710–726, 2009.

Arakawa, A. and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, Meth. Comput. Phys., 17, 173–265, 1977.

Arakawa, A. and Moorthi, S.: Baroclinic instability in vertically discrete systems, J. Atmos. Sci., 45, 1688–1708, 1988.

Ascher, U. M., Ruuth, S. J., and Spiteri, R. J.: Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Appl. Numer. Math., 25, 151–167, 1997.

Baba, Y., Takahashi, K., Sugimura, T., and Goto, K.: Dynamical core of an atmospheric general circulation model on a yin-yang grid, Mon. Weather Rev., 138, 3988–4005, 2010.

Berger, M. and Helzel, C.: A simplified h-box method for embedded boundary grids, SIAM J. Sci. Comput., 34, A861–A888, 2012.

Bermejo, R. and Conde, J.: A conservative quasi-monotone semi-Lagrangian scheme, Mon. Weather Rev., 130, 423–430, https://doi.org/10.1175/1520-0493(2002)130<0423:ACQMSL>2.0.CO;2, 2002.

Bermejo, R. and Staniforth, A.: The conversion of semi-Lagrangian advection schemes to quasi-monotone schemes, Mon. Weather Rev., 120, 2622–2631, https://doi.org/10.1175/1520-0493(1992)120<2622:TCOSLA>2.0.CO;2, 1992.

Charney, J. G. and Phillips, N.: Numerical integration of the quasi-geostrophic equations for barotropic and simple baroclinic flows, J. Meteor., 10, 71–99, 1953.

Colella, P. and Woodward, P. R.: The Piecewise Parabolic Method (PPM) for Gas-Dynamical Simulations, J. Comput. Phys., 54, 174–201, 1984.

Cotter, C. J. and Shipton, J.: Mixed finite elements for numerical weather prediction, J. Comput. Phys., 231, 7076–7091, 2012.

Davies, T., Cullen, M., Malcolm, A., Mawson, M., Staniforth, A., White, A., and Wood, N.: A new dynamical core for the Met Office’s global and regional modelling of the atmosphere, Q. J. Roy. Meteorol. Soc., 131, 1759–1782, 2005.

Deconinck, W., Bauer, P., Diamantakis, M., Hamrud, M., Kühnlein, C., Maciel, P., Mengaldo, G., Quintino, T., Raoult, B., Smolarkiewicz, P. K., and Wedi, N. P.: Atlas: A library for numerical weather prediction and climate modelling, Comput. Phys. Commun., 220, 188–204, https://doi.org/10.1016/j.cpc.2017.07.006, 2017.

de Grandpré, J., Tanguay, M., Qaddouri, A., Zerroukat, M., and McLinden, C.: Semi-Lagrangian Advection of Stratospheric Ozone on a Yin-Yang Grid System, Mon. Weather Rev., 144, 1035–1050, https://doi.org/10.1175/MWR-D-15-0142.1, 2016.

Dennis, J. M., Edwards, J., Evans, K. J., Guba, O., Lauritzen, P. H., Mirin, A. A., St-Cyr, A., Taylor, M. A., and Worley, P. H.: CAM-SE: A scalable spectral element dynamical core for the Community Atmosphere Model, Int. J. High Perf. Comput. Appl., 26, 74–89, 2012.

Dipankar, A., Stevens, B., Heinze, R., Moseley, C., Zängl, G., Giorgetta, M., and Brdar, S.: Large eddy simulation using the general circulation model ICON, J. Adv. Model. Earth Sys., 7, 963–986, 2015.

Donner, L. J., Wyman, B. L., Hemler, R. S., Horowitz, L. W., Ming, Y., Zhao, M., Golaz, J. C., Ginoux, P., Lin, S. J., Schwarzkopf, M. D., and Austin, J.: The dynamical core, physical parameterizations, and basic simulation characteristics of the atmospheric component AM3 of the GFDL global coupled model CM3, J. Climate, 24, 3484–3519, 2011.

Doyle, J. D.: A next generation atmospheric prediction system for the Navy, Tech. rep., DTIC Document, available at: http://www.dtic.mil/docs/citations/ADA616720, Naval Research Laboratory, Monterey, CA, 2014.

Du, Q., Gunzburger, M. D., and Ju, L.: Constrained centroidal Voronoi tessellations for surfaces, SIAM J. Sci. Comput., 24, 1488–1506, 2003.

Dubey, S., Dubos, T., Hourdin, F., and Upadhyaya, H. C.: On the inter-comparison of two tracer transport schemes on icosahedral grids, Appl. Math. Model., 39, 4828–4847, https://doi.org/10.1016/j.apm.2015.04.015, 2015.
Dubos, T. and Meurdesoif, Y.: DYNAMICO Model, https://doi.org/10.5281/zenodo.583718, last access: 28 November 2017.

Dubos, T. and Tort, M.: Equations of Atmospheric Motion in Non-Eulerian Vertical Coordinates: Vector-Invariant Form and Quasi-Hamiltonian Formulation, Mon. Weather Rev., 142, 3860–3880, 2014.

Dubos, T., Dubey, S., Tort, M., Mittal, R., Meurdesoif, Y., and Hourdin, F.: DYNAMICO-1.0, an icosahedral hydrostatic dynamical core designed for consistency and versatility, Geosci. Model Dev., 8, 3131–3150, https://doi.org/10.5194/gmd-8-3131-2015, 2015.

Durran, D. R.: Improving the anelastic approximation, J. Atmos. Sci., 46, 1453–1461, 1989.

Dutton, J. A.: The ceaseless wind: An introduction to the theory of atmospheric motion, Dover Publications, 1986.

Fox-Rabinovitz, M. S., Stenchikov, G. L., Suarez, M. J., and Takacs, L. L.: A finite-difference GCM dynamical core with a variable-resolution stretched grid, Mon. Weather Rev., 125, 2943–2968, 1997.

Gal-Chen, T. and Somerville, R.: Numerical solution of the Navier-Stokes equations with topography, J. Comput. Phys., 17, 276–310, https://doi.org/10.1016/0021-9991(75)90054-6, 1975.

Gardner, D. J., Guerra, J. E., Hamon, F. P., Reynolds, D. R., Ullrich, P. A., and Woodward, C. S.: Implicit-explicit (IMEX) Runge-Kutta methods for non-hydrostatic atmospheric models, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2017-285, in review, 2017.

Gassmann, A.: A global hexagonal C-grid non-hydrostatic dynamical core (ICON-IAP) designed for energetic consistency, Q. J. Roy. Meteorol. Soc., 139, 152–175, 2013.

Gassmann, A. and Herzog, H. J.: Towards a consistent numerical compressible non-hydrostatic model using generalized Hamiltonian tools, Q. J. Roy. Meteorol. Soc., 134, 1597–1613, https://doi.org/10.1002/qj.297, 2008.

Giraldo, F. X., Kelly, J. F., and Constantinescu, E.: Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA), SIAM J. Sci. Comput., 35, B1162–B1194, 2013.

Girard, C., Plante, A., Desgagné, M., McTaggart-Cowan, R., Côté, J., Charron, M., Gravel, S., Lee, V., Patoine, A., Qaddouri, A., Roch, M., Spacek, L., Tanguay, M., Vaillancourt, P., and Zadra, A.: Staggered Vertical Discretization of the Canadian Environmental Multiscale (GEM) model using a coordinate of the log-hydrostatic-pressure type, Mon. Weather Rev., 142, 1183–1196, 2014.

Guba, O., Taylor, M., and St-Cyr, A.: Optimization-based limiters for the spectral element method, J. Comput. Phys., 267, 176–195, 2014.

Guerra, J. E. and Ullrich, P. A.: A high-order staggered finite-element vertical discretization for non-hydrostatic atmospheric models, Geosci. Model Dev., 9, 2007–2029, https://doi.org/10.5194/gmd-9-2007-2016, 2016.

Harris, L. M. and Lin, S.-J.: A Two-Way Nested Global-Regional Dynamical Core on the Cubed-Sphere Grid, Mon. Weather Rev., 141, 283–306, 2013.

Harris, L. M., Lin, S.-J., and Tu, C.: High resolution climate simulations using GFDL HiRAM with a stretched global grid,

J. Climate, 29, 4293–4314, https://doi.org/10.1175/JCLI-D-15-0389.s1, 2016.

Heikes, R. P. and Randall, D. A.: Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests, Mon. Weather Rev., 123, 1862–1880, 1995.

Heikes, R. P., Randall, D. A., and Konor, C. S.: Optimized icosahedral grids: Performance of finite-difference operators and multi-grid solver, Mon. Weather Rev., 141, 4450–4469, 2013.

Holdaway, D., Thuburn, J., and Wood, N.: Comparison of Lorenz and Charney–Phillips vertical discretisations for dynamics-boundary layer coupling. Part I: Steady states, Q. J. Roy. Meteorol. Soc., 139, 1073–1086, 2013a.

Holdaway, D., Thuburn, J., and Wood, N.: Comparison of Lorenz and Charney–Phillips vertical discretisations for dynamics-boundary layer coupling. Part II: Transients, Q. J. Roy. Meteorol. Soc., 139, 1087–1098, 2013b.

Hortal, M. and Simmons, A.: Use of reduced Gaussian grids in spectral models, Mon. Weather Rev., 119, 1057–1074, 1991.

Jablonowski, C. and Williamson, D. L.: The pros and cons of diffusion, filters and fixers in atmospheric general circulation models, in: Numerical Techniques for Global Atmospheric Models, 381–493, Springer, 2011.

Jablonowski, C., Lauritzen, P. H., Taylor, M. A., and Nair, R. D.: Idealized test cases for the dynamical cores of atmospheric general circulation models: a proposal for the NCAR ASP 2008 summer colloquium, Tech. rep., 2008.

Kageyama, A. and Sato, T.: The Yin-Yang grid: An overset grid in spherical geometry, Geochem. Geophys. Geosyst., 5, Q09005, https://doi.org/10.1029/2004GC000734, 2004.

Kasahara, A.: Various vertical coordinate systems used for numerical weather prediction, Mon. Weather Rev., 102, 509–522, 1974.

Klemp, J. B.: A terrain-following coordinate with smoothed coordinate surfaces, Mon. Weather Rev., 139, 2163–2169, 2011.

Klemp, J. B., Skamarock, W. C., and Dudhia, J.: Conservative split-explicit time integration methods for the compressible nonhydrostatic equations, Mon. Weather Rev., 135, 2897–2913, 2007.

Klemp, J. B., Dudhia, J., and Hassiotis, A. D.: An upper gravity-wave absorbing layer for NWP applications, Mon. Weather Rev., 136, 3987–4004, https://doi.org/10.1175/2008MWR2596.1, 2008.

Kühnlein, C. and Smolarkiewicz, P. K.: An unstructured-mesh finite-volume MPDATA for compressible atmospheric dynamics, J. Comput. Phys., 334, 16–30, 2017.

Kühnlein, C., Smolarkiewicz, P. K., and Dörnbrack, A.: Modelling atmospheric flows with adaptive moving meshes, J. Comput. Phys., 231, 2741–2763, 2012.

Kurowski, M. J., Grabowski, W. W., and Smolarkiewicz, P. K.: Anelastic and compressible simulation of moist deep convection, J. Atmos. Sci., 71, 3767–3787, 2014.

Laprise, R.: The Euler equations of motion with hydrostatic pressure as an independent variable, Mon. Weather Rev., 120, 197–207, 1992.

Lin, S.-J.: A finite-volume integration method for computing pressure gradient force in general vertical coordinates, Q. J. Roy. Meteorol. 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.: An explicit flux-form semi-Lagrangian shallow-water model on the sphere., Q. J. Roy. Meteorol. Soc., 123, 2477–2498, 1997.

Lloyd, S.: Least squares quantization in PCM, IEEE Trans. Info. Theory, 28, 129–137, 1982.

Lock, S.-J., Bitzer, H.-W., Coals, A., Gadian, A., and Mobbs, S.: Demonstration of a cut-cell representation of 3D orography for studies of atmospheric flows over very steep hills, Mon. Weather Rev., 140, 411–424, 2012.

Lorenz, E. N.: Energy and numerical weather prediction, Tellus, 12, 364–373, 1960.

Malardel, S., Wedi, N., Deconinck, W., Diamantakis, M., Kühnlein, C., Mozdzynski, G., Hamrud, M., and Smolarkiewicz, P.: A new grid for the IFS, ECMWF Newsletter, 146, 23–28, 2016.

Miura, H.: An upwind-biased conservative advection scheme for spherical hexagonal–pentagonal grids, Mon. Weather Rev., 135, 4038–4044, 2007.

Nair, R., Choi, H.-W., and Tufo, H.: Computational aspects of a scalable high-order discontinuous Galerkin atmospheric dynamical core, Computers & Fluids, 38, 309–319, 2009.

Niwa, Y., Tomita, H., Satoh, M., and Imasu, R.: A three-dimensional icosahedral grid advection scheme preserving monotonicity and consistency with continuity for atmospheric tracer transport, J. Meteor. Soc. Japan. Ser. II, 89, 255–268, 2011.

Ogura, Y. and Phillips, N. A.: Scale analysis of deep and shallow convection in the atmosphere, J. Atmos. Sci., 19, 173–179, 1962.

Okabe, A., Boots, B., Sugihara, K., and Chiu, S. N.: Spatial tesselations: concepts and applications of Voronoi diagrams, vol. 501, John Wiley & Sons, 2009.

Orlanski, I.: The quasi-hydrostatic approximation, J. Atmos. Sci., 38, 572–582, 1981.

Perot, B.: Conservation properties of unstructured staggered mesh schemes, J. Comput. Phys., 159, 58–89, 2000.

Phillips, N. A.: The equations of motion for a shallow rotating atmosphere and the “traditional approximation”, J. Atmos. Sci., 23, 626–628, 1966.

Prusa, J. M., Smolarkiewicz, P. K., and Wyszogrodzki, A. A.: EU-LAG, a computational model for multiscale flows, Computers & Fluids, 37, 1193–1207, 2008.

Putman, W. M. and Lin, S.-J.: Finite-volume transport on various cubed-sphere grids, J. Comput. Phys., 227, 55–78, 2007.

Qaddouri, A. and Lee, V.: The Canadian Global Environmental Multiscale model on the Yin-Yang grid system, Q. J. Roy. Meteorol. Soc., 137, 1913–1926, 2011.

Qaddouri, A., Laayouni, L., Loisel, S., Côté, J., and Gander, M.: Optimized Schwarz methods with an overset grid for the shallow-water equations: Preliminary results, Appl. Numer. Math., 58, 459–471, 2008.

Randall, D. A.: Geostrophic adjustment and the finite-difference shallow-water equations, Mon. Weather Rev., 122, 1371–1377, 1994.

Randall, D., Dazlich, D., and Heikes, R.: CSU Model, https://doi.org/10.5281/zenodo.580099, last access: 28 November 2017.

Ringler, T. D. and Randall, D. A.: The ZM grid: an alternative to the Z grid, Mon. Weather Rev., 130, 1411–1422, 2002.

Ringler, T. D., Ju, L., and Gunzburger, M.: A multiresolution method for climate system modeling: application of spherical centroidal Voronoi tessellations, Ocean Dyn., 58, 475–498, 2008.

Ringler, T. D., Thuburn, J., Klemp, J. B., and Skamarock, W. C.: A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids, J. Comput. Phys., 229, 3065–3090, 2010.

Ronchi, C., Ianoco, R., and Paolucci, P. S.: The ‘Cubed Sphere’: a new method for the solution of partial differential equations in spherical geometry, J. Comput. Phys., 124, 93–114, https://doi.org/10.1006/jcph.1996.0047, 1996.

Sadourny, R.: Conservative finite-difference approximations of the primitive equation on quasi-uniform spherical grids, Mon. Weather Rev., 100, 136–144, https://doi.org/10.1175/1520-0493(1972)100<0136:CFAOTP>2.3.CO;2, 1972.

Satoh, M.: Conservative scheme for the compressible nonhydrostatic models with the horizontally explicit and vertically implicit time integration scheme, Mon. Weather Rev., 130, 1227–1245, 2002.

Satoh, M.: Conservative scheme for a compressible nonhydrostatic model with moist processes, Mon. Weather Rev., 131, 1033–1050, 2003.

Satoh, M., Matsuno, T., Tomita, H., Miura, H., Nasuno, T., and Iga, S.-I.: Nonhydrostatic icosahedral atmospheric model (NICAM) for global cloud resolving simulations, J. Comput. Phys., 227, 3486–3514, 2008.

Schär, C., Leuenberger, D., Fuhrer, O., Lüthi, D., and Girard, C.: A new terrain-following vertical coordinate formulation for atmospheric prediction models, Mon. Weather Rev., 130, 2459–2480, 2002.

Schmidt, F.: Variable fine mesh in spectral global models, Beitr. Phys. Atmos., 50, 211–217, 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.

Skamarock, W. C.: Evaluating mesoscale NWP models using kinetic energy spectra, Mon. Weather Rev., 132, 3019–3032, 2004.

Skamarock, W. C. and Gassmann, A.: Conservative transport schemes for spherical geodesic grids: High-order flux operators for ODE-based time integration, Mon. Weather Rev., 139, 2962–2975, 2011.

Skamarock, W. C. and Klemp, J. B.: The stability of time-split numerical methods for the hydrostatic and the nonhydrostatic elastic equations, Mon. Weather Rev., 120, 2109–2127, 1992.

Skamarock, W. C., Klemp, J. B., Duda, M. G., Fowler, L. D., Park, S.-H., and Ringler, T. D.: A multiscale nonhydrostatic atmospheric model using centroidal Voronoi tessellations and C-grid staggering, Mon. Weather Rev., 140, 3090–3105, 2012.

Skamarock, W., Klemp, J., and Park, S.-H.: MPAS Model, https://doi.org/10.5281/zenodo.583316, last access: 28 November 2017.

Smagorinsky, J.: General circulation experiments with the primitive equations: I. The basic experiment, Mon. Weather Rev., 91, 99–164, 1963.

Smolarkiewicz, P. K. and Szmelter, J.: MPDATA: an edge-based unstructured-grid formulation, J. Comput. Phys., 206, 624–649, 2005.

Smolarkiewicz, P. K. and Szmelter, J.: A nonhydrostatic unstructured-mesh soundproof model for simulation of internal gravity waves, Acta Geophys., 59, 1109–1134, 2011.
Smolarkiewicz, P. K., Kühnlein, C., and Wedi, N. P.: A consistent framework for discrete integrations of soundproof and compressible PDEs of atmospheric dynamics, J. Comput. Phys., 263, 185–205, 2014.

Smolarkiewicz, P. K., Deconinck, W., Hamrud, M., Kühnlein, C., Mozdzynski, G., Szmelter, J., and Wedi, N. P.: A finite-volume module for simulating global all-scale atmospheric flows, J. Comput. Phys., 314, 287–304, 2016.

Smolarkiewicz, P. K., Kühnlein, C., and Grabowski, W.: A finite-volume module for cloud-resolving simulations of global atmospheric flows, J. Comput. Phys., 341, 208–229, 2017.

Sørensen, B., Kaas, E., and Korsholm, U. S.: A mass-conserving and multi-tracer efficient transport scheme in the online integrated Enviro-HIRLAM model, Geosci. Model Dev., 6, 1029–1042, https://doi.org/10.5194/gmd-6-1029-2013, 2013.

Staniforth, A. and Thuburn, J.: Horizontal grids for global weather and climate prediction models: a review, Q. J. Roy. Meteorol. Soc., 138, 1–26, 2012.

Starr, V. P.: A quasi-Lagrangian system of hydrodynamical equations, J. Meteorol., 2, 227–237, 1945.

Szmelter, J. and Smolarkiewicz, P. K.: An edge-based unstructured mesh discretisation in geospherical framework, J. Comput. Phys., 229, 4980–4995, 2010.

Taylor, M. A. and Fournier, A.: A compatible and conservative spectral element method on unstructured grids, J. Comput. Phys., 229, 5879–5895, 2010.

Thuburn, J.: A PV-based shallow-water model on a hexagonal-icosahedral grid, Mon. Weather Rev., 125, 2328–2347, 1997.

Thuburn, J.: Some conservation issues for the dynamical cores of NWP and climate models, J. Comput. Phys., 227, 3715–3730, 2008.

Thuburn, J. and Woollings, T.: Vertical discretizations for compressible Euler equation atmospheric models giving optimal representation of normal modes, J. Comput. Phys., 203, 386–404, 2005.

Thuburn, J., Ringler, T. D., Skamarock, W. C., and Klemp, J. B.: Numerical representation of geostrophic modes on arbitrarily structured C-grids, J. Comput. Phys., 228, 8321–8335, 2009.

Tokioka, T.: Some considerations on vertical differencing, Meteorol. Soc. Jpn. Journal, 56, 98–111, 1978.

Tomita, H. and Satoh, M.: A new dynamical framework of nonhydrostatic global model using the icosahedral grid, Fluid Dynam. Res., 34, 357–400, 2004.

Tomita, H., Tsugawa, M., Satoh, M., and Goto, K.: Shallow water model on a modified icosahedral geodesic grid by using spring dynamics, J. Comput. Phys., 174, 579–613, 2001.

Tomita, H., Satoh, M., and Goto, K.: An optimization of the icosahedral grid modified by spring dynamics, J. Comput. Phys., 183, 307–331, 2002.

Ullrich, P. A.: A global finite-element shallow-water model supporting continuous and discontinuous elements, Geosci. Model Dev., 7, 3017–3035, https://doi.org/10.5194/gmd-7-3017-2014, 2014a.

Ullrich, P. A.: Understanding the treatment of waves in atmospheric models. Part 1: The shortest resolved waves of the 1D linearized shallow-water equations, Q. J. Roy. Meteorol. Soc., 140, 1426–1440, 2014b.

Ullrich, P.: Tempest Model, https://doi.org/10.5281/zenodo.579649, last access: 28 November 2017.

Ullrich, P. A. and Jablonowski, C.: MCore: A non-hydrostatic atmospheric dynamical core utilizing high-order finite-volume methods, J. Comput. Phys., 231, 5078–5108, 2012a.

Ullrich, P. A. and Jablonowski, C.: Operator-split Runge-Kutta-Rosenbrock methods for nonhydrostatic atmospheric models, Mon. Weather Rev., 140, 1257–1284, 2012b.

Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R. D., and Taylor, M. A.: Dynamical core model intercomparison project (DCMIP) test case document, Tech. rep., 2012.

Wacker, U., Frisius, T., and Herbert, F.: Evaporation and Precipitation Surface Effects in Local Mass Continuity Laws of Moist Air, J. Atmos. Sci., 63, 2642–2652, 2006.

Walko, R.: OLAM Model, https://doi.org/10.5281/zenodo.582308, last access: 28 November 2017.

Walko, R. L. and Avissar, R.: The ocean–land–atmosphere model (OLAM). Part I: Shallow-water tests, Mon. Weather Rev., 136, 4033–4044, 2008a.

Walko, R. L. and Avissar, R.: The Ocean–Land–Atmosphere Model (OLAM). Part II: Formulation and tests of the nonhydrostatic dynamic core, Mon. Weather Rev., 136, 4045–4062, 2008b.

Walko, R. L. and Avissar, R.: A direct method for constructing refined regions in unstructured conforming triangular–hexagonal computational grids: Application to OLAM, Mon. Weather Rev., 139, 3923–3937, 2011.

Walko, R. L., Cotton, W. R., Feingold, G., and Stevens, B.: Efficient computation of vapor and heat diffusion between hydrometeors in a numerical model, Atmos. Res., 53, 171–183, 2000.

Weller, H.: Controlling the computational modes of the arbitrarily structured C grid, Mon. Weather Rev., 140, 3220–3234, 2012.

Weller, H., Thuburn, J., and Cotter, C. J.: Computational modes and grid imprinting on five quasi-uniform spherical C grids, Mon. Weather Rev., 140, 2734–2755, 2012.

Weller, H., Lock, S.-J., and Wood, N.: Runge–Kutta IMEX schemes for the horizontally explicit/vertically implicit (HEVI) solution of wave equations, J. Comput. Phys., 252, 365–381, 2013.

Wicker, L. J. and Skamarock, W. C.: Time-splitting methods for elastic models using forward time schemes, Mon. Weather Rev., 130, 2088–2097, 2002.

Wood, N., Staniforth, A., White, A., Allen, T., Diamantakis, M., Gross, M., Melvin, T., Smith, C., Vosper, S., Zerroukat, M., and Thuburn, J.: An inherently mass-conserving semi-implicit semi-Lagrangian discretization of the deep-atmosphere global non-hydrostatic equations, Q. J. Roy. Meteorol. Soc., 140, 1505–1520, 2014.

Zängl, G., Reinert, D., Rípodas, P., and Baldauf, M.: The ICON (ICOsahedral Non-hydrostatic) modelling framework of DWD and MPI-M: Description of the non-hydrostatic dynamical core, Q. J. Roy. Meteorol. Soc., 141, 563–579, https://doi.org/10.1002/qj.2378, 2015.