DCMIP2016: A Review of Non-hydrostatic Dynamical Core Design and Intercomparison of Participating Models

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 Naviér-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 5 aims to provide a comprehensive review of eleven 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 choice of model grid, variable placement, vertical coordinate, prognostic equations, temporal discretization, and the diffusion, stabilization, filters and fixers employed by each system.


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

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 staggering 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.  Virtual potential temperature θ il Ice-liquid potential temperature θ ρ Density potential temperature q Specific humidity q v Water vapor mixing ratio q c Cloud water mixing ratio q r Rain water mixing ratio q i General tracer mixing ratio 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 shallowatmosphere 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 cubedsphere 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.

Colorado State University (CSU) model
The Colorado State University (CSU) model is a finitevolume 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 (∇ 4 z ) 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.

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. DY-NAMICO's design uniquely combines a representation of the prognostic and diagnostic fields following the ideas of discrete differential geometry . 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 energyconserving 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-    ity operator (Sect. 6.3). In addition, it features a conservative positive-definite transport scheme based on a slope-limited finite-volume approach .

FV cubed (FV 3 )
The GFDL Finite-Volume Cubed-Sphere Dynamical Core (FV 3 , 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 piecewise 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.  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 terrainfollowing 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 twotime-level, semi-implicit integration scheme is employed with 3-D implicit treatment of acoustic, buoyant, and rotational modes (Smolarkiewicz et al., 2014) (Sect. 7.2). The associated 3-D Helmholtz problem is solved iteratively using a bespoke preconditioned generalized conjugate residual approach. The integration procedure uses the nonoscillatory, 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 stabilize the model, so no other explicit filtering mechanism is required (Sect. 6.5). Note that the octahedral reduced Gaussian grid is also employed in the spectral-transform dynamical 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 towards broad classes of meshes (Szmelter and Smolarkiewicz, 2010;Kühnlein et al., 2012;Deconinck et al., 2017).

Global Environmental Multiscale (GEM) model
The 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 Schwarztype 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).

ICOsahedral Non-hydrostatic (ICON) model
The ICOsahedral Non-hydrostatic (ICON) model  is a finite-volume model that solves the nonhydrostatic 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 terrainfollowing 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.

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 Cgrid 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 ( ρ d ), dry momentum ( ρ 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).

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 (2002Satoh ( , 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-secondorder 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 piecewise linear transport scheme (Miura, 2007;Niwa et al., 2011) is used with a flux limiter of Thuburn (1997) for passive tracer transports.

Ocean-Land-Atmosphere Model (OLAM)
Ocean-Land-Atmosphere Model (OLAM) (Walko andAvissar, 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.

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

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

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 (α, β, p), with central angles α, β ∈ [− π 4 , π 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 cubedsphere 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.

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 2 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 · n 2 · 2 2m cells, n e = 3/2 · n c edges, and n v = 10 · n 2 · 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 2 949 120 cells and 13 km mean resolution in combination with a refined region over Europe at 6.5 km resolution.

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

Constrained centroidal Voronoi tessellation (CCVT) meshes
Given a set of N distinct points on the sphere x i (referred to as the generators, 1 ≤ i ≤ N), the Voronoi tessellation (or the Voronoi diagram) associated with the generators is the set of polygons 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 = 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 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.

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

Yin-Yang grid
The overset Yin-Yang grid (Kageyama and Sato, 2004) has two Cartesian grid components (subsets of a latitudelongitude 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 (1) where δ λ , δ θ 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 δ θ = 2 • and δ λ = 3δ θ . 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 latitudelongitude model.

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 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 3 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).
. Horizontal staggering options represented among DCMIP models, in this case depicted on a rectilinear grid and geodesic grid.
Here, η denotes the buoyancy variable.

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.

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 Interface

Vertical coordinates
In the context of dimension splitting, the "horizontal" typically refers to either the contravariant basis, which is perpendicular to the vertical, or the covariant basis, which is directed along coordinate (e.g., terrain-following) surfaces. In contrast, the vertical dimension is strictly aligned with the radial vector pointing from the center of the Earth. Vertical position is typically labeled using an arbitrary function s(t, x, z) that is monotonic in z, so that model interfaces are equally spaced with respect to s. Typically, s is chosen so that the Earth's surface (the bottom boundary of the atmosphere) is a coordinate surface, allowing for easy specification of boundary conditions for the prognostic equations; this leads to the so-called "terrain-following" family of vertical coordinates. Perhaps the most common terrain-following coordinate is from Gal-Chen and Somerville (1975), which is in terms of the altitude z and takes the form where x denotes the horizontal position, z s (x) is the height of the topography at that position, and z top denotes the height of the model top (typically independent of position). Analogous formulations are available for mass-based (σ coordinates) and entropy-based vertical coordinates. Because the sharp variations in the coordinate surfaces are preserved far above a rough lower boundary, new coordinate formulations have been proposed that smooth coordinate surfaces, such as Schär et al. (2002) or Klemp (2011). All models in this paper except for OLAM use some variant of terrain-following coordinates, although work on developing modern cut-cell, embedded boundary and immersed boundary representations is ongoing (e.g., Lock et al., 2012). Note that time-dependent vertical coordinates are allowed and are typically referred to as "floating" coordinates. Several examples of vertical coordinates are now given.

Mass-based coordinates
Mass-based coordinates (Laprise, 1992) are a generalization of pressure-based coordinates to non-hydrostatic models, with a vertical coordinate defined as the total gravityweighted overhead mass, Under this definition,

GEM ζ coordinate
The vertical coordinate in the GEM model, denoted ζ , is a hybrid terrain-following coordinate of a log-hydrostaticpressure type. Taking s (denoted π in GEM documentation) as given in Eq.
(3), ζ is given by the relation with Here, ζ s = log(10 5 ), ζ top = log(s top ), s top is the coordinate value at the uppermost interface, and r is a variable exponent providing added freedom for adjusting the thickness of model layers over high terrain.

Floating Lagrangian coordinates (ACME-A, DYNAMICO, and FV 3 )
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., in which case a 3-D reference density field ρ 0 can be imposed. Of primary importance is the fact that the vertical coordinate satisfieṡ 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.

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.

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.

ACME-A
ACME-A presently solves the compressible shallowatmosphere equations using a hybrid terrain-following pressure vertical coordinate η, 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:

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.

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

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 Dependent variables in Eq. (10) are dry density ρ d , 3-D physical velocity vector u, potential temperature perturbation θ , and a modified Exner pressure perturbation φ , 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 ε b = 1/ε − 1. As far as geometric aspects are concerned, the nabla operator ∇ represents the 3-D vector of partial derivatives with respect to the curvilinear coordinates, along with the Jacobian G, a matrix of metric coefficients G, its transpose G T , and the contravariant velocity v = G T u, where a contribution from optional time dependency of the curvilinear coordinates is neglected for simplicity . The symbol M = M (u, u a , θ ρ /θ ρ a ) in Eq. (10b) subsumes the metric forces in the spherical domain .

GEM
In GEM, the non-hydrostatic equations are written explicitly as deviations from hydrostatic balance represented by where s (denoted π in GEM documentation) is given by Eq.
(3). In this case, the equations of GEM model (Girard et al., 2014) are concisely given by Here, ∇ ζ denotes the horizontal gradient along ζ 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).

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 u bc = k ρ k u k / k ρ 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 = ρ k /ρ of all constituents (except for dry air) for which the prognostic equation reads with σ k describing sources/sinks due to phase changes, and J k = ρq k (u k − u bc ) denoting diffusion fluxes, which account for the motion of constituents relative to the frame of reference set by u bc . The specific heat capacities and ideal gas constant are approximated to be equal to their dry values R * ≈ R d , c * p ≈ c pd , c * v ≈ c vd . The model also uses a prognostic equation for Exner pressure to simplify the treatment of vertical sound wave propagation, given by whereQ 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 u bc to zero at the lower boundary: w bc | sfc = 0. In addition, only the vertical diffusion fluxes J k of sedimentation constituents and the surface evaporation flux J v | sfc are taken into account. The counter flux of nonsedimentation constituents is discarded. Since in the given framework the continuity Eq. (A10) is only valid if the constraint k 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.

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 ρ s . MPAS further evolves dry mass using a continuity equation of Eq. (A10) and moist potential temperature following Eq. (A13).

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

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 θ il (Walko et al., 2000), total density ρ, 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 θ . The latter is in turn diagnosed from θ 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.

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.

ACME-A/Tempest
In both ACME-A and Tempest, scalar hyperviscosity is employed for ρ, θ, 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 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).

CSU
The CSU model utilizes an explicit diffusion scheme that consists of fourth-order hyperdiffusion (∇ 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.

DYNAMICO
In DYNAMICO, (hyper-)diffusive filters are used to eliminate spurious noise due to the energy-conservative centered discretization. Filters are applied every N diff Runge-Kutta time step in a forward-Euler manner, with N 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 τ : given τ , the hyperviscous coefficient that multiplies operator D p is δ 2p τ − 1, where δ − 2 is the largest eigenvalue of operator D. For DCMIP, DYNAMICO uses p = 2 (fourth-order hyperviscosity) for all filters.

FV 3
Explicit dissipation in FV 3 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 3 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 3 is constructed with a flexible-lid (constant-pressure) upper boundary that is effective at damping internal gravity wave modes; however, FV 3 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 3 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 Exnerfunction weighted potential temperature increment. This is only done before the physics is called and is not used in idealized simulations.

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.

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

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 typically attains values between 1 1000 t and 1 250 t , and 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 (θ )) combined with a fourthorder background diffusion of velocity F D4 (v n ), defined via 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 4 a 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

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 where u i is the edge-normal velocity defined on cell edge i, η is the vertical component of the relative vorticity, computed on vertices, and ∇ s · 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): where c s is a constant parameter and is the grid scale. The diffusion operator then takes the form ∇·(K h ∇ψ) for a scalar field ψ.

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 4494 P. A. Ullrich et al.: DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison (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).

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

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 nonhydrostatic 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 andPhillips, 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.

Mixed implicit-explicit, forward-backward, semi-implicit, and additive Runge-Kutta schemes
Implicit-explicit schemes are a broad category of timeintegration 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 3 and its predecessors are integrated using a forwardbackward 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 predictorcorrector 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 cellto-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 diagnostically 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 trapezoidalimplicit 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.

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. (2014Smolarkiewicz et al. ( , 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 moistprecipitating 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 where ψ represents the solution variable, R ψ 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 ≡ (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 ψ ≡ ρ d , v ≡ vG, G ≡ G, and R ρ d ≡ 0 in Eq. (26). Second, given already updated moisture variables , the thermodynamic (Eq. 10c) and momentum (Eq. 10b) equations enter Eq. (26) with ψ = u, v, w, θ , v ≡ vGρ d , G ≡ Gρ d , and the right-hand-side R ψ 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 closedform 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 φ using an advective form of the equation of state (Eq. 10d). An O( t 2 ) integration of this equation with a Euler backward scheme, in the spirit of Eq. (26), leads to a Helmholtz equation . 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 ψ | n+1 and the solution-dependent coefficients of the Helmholtz equation are lagged behind and executed in an outer iteration.

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 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): where "A" stands for the arrival position at model grid point (r h , ζ, t) and "D" for the departure position (r h − r h , ζ − ζ, t − t) due to the displacements r h , ζ having occurred during the time step t. G is averaged between these two positions with a possible slight off-centering . The displacements are themselves calculated by solving, again using the Lagrangian method, the following equations: discretized in the same way (trapezoidal method) though without off-centering: The process is of course a multi-step iterative one, since both positions and velocities at departure positions (past time t − 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: where τ = (1/2 + ) t and β = (1/2 − ) / (1/2 + ). Secondly, separating linear from non-linear parts, we get with 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 Schwarztype 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.

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, horizontal staggering, vertical staggering, vertical coordinates, prognostic equations, choice of diffusion, stabilization, filters and fixers, and temporal discretization. The literature on dynamical 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 andJablonowski (2012a), Gassmann (2013), Wood et al. (2014), andDoyle (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.   or similarly for θ. In conjunction with the material derivative of the ideal gas law, the thermodynamic equation can be written in the form Then, substituting Eq. (A9) gives a prognostic equation for virtual temperature: 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: and similarly for Exner pressure:

A3 Momentum equations
In coordinate-invariant form, the prognostic velocity equations may be written in either the Lagrangian or Eulerian frame as where denotes the planetary vorticity vector and 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 u · ∇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 which follows from Eq. (A4). Note that often in nonhydrostatic models, κ is assumed constant and the ∇κ term neglected. A second form of Eq. (A19) emerges on substituting the vector calculus identity where K = 1 2 (u · u) is the 3-D specific kinetic energy and ζ = ∇ × u is the 3-D relative vorticity vector. This gives rise to the 3-D vector-invariant form, Because no gradients of vectors appear in this equation, it avoids derivatives of the coordinate basis that would arise from the momentum transport term u · ∇u in Eq. (A19). In conjunction with Eq. (A9), Eq. (A19) also gives rise to the flux-form momentum equations, where u ⊗ u denotes the outer product and I is the identity matrix.
The equations above still provide some flexibility with regard to the choice of and . For deep-atmosphere models, one typically chooses = g c a 2 1 a − 1 a + z , and = (k sin ϕ + j cos ϕ), where g c is gravitational acceleration at the surface, a is the radius of the planet, is the rotation rate (in s −1 ), ϕ 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, For shallow-atmosphere models, the geopotential takes the simpler form = g c z, and = sin ϕk, where z is the altitude above the surface. In this case, we write 2 = f k, where f = 2 sin ϕ is the Coriolis parameter. The evolution equation for the shallow-atmosphere geopotential is then

A4 Orthogonal formulation
Under the orthogonal formulation, projection of a vector field b onto its horizontal components is defined via When applied to the velocity vector, this gives rise to the decomposition where k = ∇z is the unit vector in the vertical direction and u h = [u] z (u h is aligned with surfaces of constant z). In the orthogonal formulation, the material derivative expands as For the special case of the material derivative applied to scalars, this equation can also be written as where ∇ z b = [∇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 which, from u h = u − wk, then gives rise to Due to its association with hydrostatic models, it is common to use the 2-D kinetic energy, K 2 = 1 2 (u h · u h ). Decomposing the momentum transport term into horizontal and vertical components gives The first term in this expression admits the relationships where ζ h = (∇ × u) · k = (∇ × u h ) · k is the relative vorticity scalar and u t = k × 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 Then, subtracting Eq. (A37) ·k from Eq. (A19) gives Note that, under the shallow-atmosphere approximation, the metric term u h · (u h · ∇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 ∂s/∂z = 0 everywhere, i.e., following Kasahara (1974) From here, Eq. (A39) also gives rise to which can be used directly to rewrite Eq. (A33) or Eq. (A38) in terms of derivatives over s. Note that the operators ∇ z and ∇ 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 · ∇ z A = 0 and k · ∇ 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 and that the 2-D curl is given by where ∇ z × u h = k(k · (∇ × u h )). Notably, these expressions are valid for both shallow-and deep-atmosphere formulations.
The generalized velocityṡ following a fluid parcel is defined bẏ

A6 Conservation laws in arbitrary vertical coordinates
Using Eq. (A42), we observe that the 3-D divergence on the sphere takes the form where α = 1 for shallow-atmosphere models and α = r 2 = (a + z) 2 for deep-atmosphere models. Using w = dz/dt, this last term also takes the form Hence, for any quantity that is conserved following a fluid parcel (i.e., dq/dt = 0), ∂ρ s q ∂t s + ∇ s · (ρ s qu h ) + ∂ ∂s (ρ s qṡ) = 0.
In particular, the prognostic equation for virtual potential temperature (or equivalently for potential temperature) reads Observe that both of these equations simplify when w = (∂z/∂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ṡ. Substituting Eq. (A44) into Eq. (A58) then gives ∂w ∂t s =u h · (u h · ∇k) − u h · ∇ s w −ṡ ∂w ∂s Similarly, substituting Eq. (A44) into Eq. (A59) and using the identity (u h · ∇ s z) ∂s ∂z where ∇ s × u h = kζ s , and ζ s = k · (∇ s × u h ).
In this case, the vertical advection terms are removed wheṅ s = 0, i.e., the vertical coordinate is advected with the 3-D wind u.
Note that under the shallow-atmosphere approximation, the first metric terms (those that include (u h · ∇k)) in Eqs. (A58)-(A62) are typically dropped.
where we have used k ·(k ·∇u) = k ·∇w. Similarly, the prognostic equation for horizontal velocity from Eq. (A33) is reformulated as Note that the vorticity term in this expression can be simplified further using and −k × ζ = k · ∇u − ∇(k · u) + u · ∇k = ∂u h ∂z − ∇ z w + u h · ∇k. (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 (A76)

A10 Vorticity divergence form
The vorticity divergence form of the dynamical equations in an arbitrary vertical coordinate predicts the absolute vorticity (ζ * h ) and velocity divergence (D) given by and respectively, instead of the horizontal velocity. The horizontal velocity can be obtained from the stream function ψ and the velocity potential χ following u h = k × ∇ s ψ + ∇ s χ .

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 ρ s and assume no external sources or sinks of ρ. Multiplying Eq. (A60) by ρ s and using