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

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

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 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 . Section  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.  describes aspects of the horizontal discretization, including model grids and horizontal placement of prognostic variables; Sect.  describes the vertical placement of model variables and choice of vertical coordinates; Sect.  describes aspects of variable placement and prognosis; Sect.  describes diffusion, stabilization, filters, and fixers employed by these models; and Sect.  describes temporal discretizations. The summary and conclusions then follow in Sect. . Finally, Appendix  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.

A standard list of symbols used throughout this paper and in the DCMIP.

Symbol	Description
$\mathit{\lambda }$	Longitude (in radians)
$\mathit{\phi }$	Latitude (in radians)
$z$	Height with respect to mean sea level (set to zero)
$s$	Vertical model coordinate
$p_{\mathrm{s}}$	Surface pressure ($p_{\mathrm{s}}$ of moist air if $q>\mathrm{0}$)
$\mathrm{\Phi }$	Geopotential
${\mathrm{\Phi }}_{\mathrm{s}}$	Surface geopotential
$z_{\mathrm{s}}$	Surface elevation with respect to mean sea level (set to zero)
$u$	Zonal wind velocity
$v$	Meridional wind velocity
$w$	Vertical wind velocity
$\dot{\mathit{\zeta }}$	GEM vertical coordinate velocity
$\mathit{u}$	3-D wind vector
${\mathit{u}}_{\mathrm{h}}$	Horizontal wind vector
${\mathit{v}}_{\mathrm{h}}$	Horizontal wind vector with covariant components
$\mathit{\omega }$	Vertical pressure velocity
$D$	Divergence of the horizontal wind vector
$\mathit{\zeta }$	Vertical component of relative vorticity
$p$	Pressure (pressure of moist air if $q>\mathrm{0}$)
$e$	Internal energy
$\mathit{\rho }$	Total air density
${\mathit{\rho }}_{\mathrm{d}}$	Dry air density
${\mathit{\rho }}_{\mathrm{s}}$	Pseudo-density
$T$	Temperature
$T_{\mathrm{v}}$	Virtual temperature
$\mathit{\theta }$	Potential temperature
${\mathit{\theta }}_{\mathrm{v}}$	Virtual potential temperature
${\mathit{\theta }}_{\mathrm{il}}$	Ice–liquid potential temperature
${\mathit{\theta }}_{\mathit{\rho }}$	Density potential temperature
$q$	Specific humidity
$q_{\mathrm{v}}$	Water vapor mixing ratio
$q_{\mathrm{c}}$	Cloud water mixing ratio
$q_{\mathrm{r}}$	Rain water mixing ratio
$q_i$	General tracer mixing ratio

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 ). The prognostic variables employed and horizontal discretizations for these dynamical cores are summarized in Table . The vertical staggering of variables and vertical coordinate choice are summarized in Table . Principal options for diffusion, stabilization, filters, or fixers along with the temporal discretization for these models are summarized in Table . A brief description of each participant model follows, focused on the unique features and decisions underlying the model design.

Participating modeling centers and associated dynamical cores that have submitted a model description and/or simulation results.

Short name	Long name	Modeling center or group
ACME-A	Atmosphere model of the Accelerated	Sandia National Laboratories and
	Climate Model for Energy	University of Colorado, Boulder, USA
CSU	Colorado State University Model	Colorado State University, USA
DYNAMICO	DYNAMical core on the ICOsahedron	Institut Pierre Simon Laplace (IPSL), France
FV${}^{\mathrm{3}}$	GFDL Finite-Volume Cubed-Sphere Dynamical Core	Geophysical Fluid Dynamics Laboratory, USA
FVM	Finite Volume Module of the Integrated Forecasting System	European Centre for Medium-Range Weather Forecasts
GEM	Global Environmental Multiscale model	Environment and Climate Change Canada, Canada
ICON	ICOsahedral Non-hydrostatic model	Max-Planck-Institut für Meteorologie, Germany
MPAS	Model for Prediction Across Scales	National Center for Atmospheric Research, USA
NICAM	Non-hydrostatic Icosahedral Atmospheric Model	AORI/JAMSTEC/AICS, Japan
OLAM	Ocean Land Atmosphere Model	Duke University/University of Miami, USA
Tempest	Tempest Non-hydrostatic Atmospheric Model	University of California, Davis, USA

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

Short name	Equation set	Prognostic variables	Horizontal grid	Numerical	Horizontal
				method	staggering
ACME-A	H/NH	${\mathit{u}}_{\mathrm{h}}$, $w$, ${\mathit{\rho }}_{\mathrm{s}}$, ${\mathit{\rho }}_{\mathrm{s}}\mathit{\theta }$, $\mathrm{\Phi }$, ${\mathit{\rho }}_{\mathrm{s}}{q}_i$	Cubed sphere (Sect. )	SE	A grid
CSU	NH (unified)	$\mathit{\zeta }$, $D$, $w$, $p_{\mathrm{s}}$, ${\mathit{\theta }}_{\mathrm{v}}$, $q_i$	Geodesic (Sect. )	FV	Z grid
DYNAMICO	H/NH	${\mathit{v}}_{\mathrm{h}}$, ${\mathit{\rho }}_{\mathrm{s}}w$, ${\mathit{\rho }}_{\mathrm{s}}$, ${\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}$, $\mathrm{\Phi }$, ${\mathit{\rho }}_{\mathrm{s}}{q}_i$	Geodesic (Sect. )	FV	C grid
FV${}^{\mathrm{3}}$	NH	${\mathit{u}}_{\mathrm{h}}$, $w$, ${\mathit{\rho }}_{\mathrm{s}}$, ${\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}$, $\mathrm{\Phi }$, ${\mathit{\rho }}_{\mathrm{s}}{q}_i$	Cubed sphere (Sect. )	FV	D grid
FVM	NH (D)	${\mathit{\rho }}_{\mathrm{d}}$, ${\mathit{u}}_{\mathrm{h}}$, $w$, ${\mathit{\theta }}^{\prime }$, $q_i$	Octahedral (Sect. )	FV	A grid
GEM	NH	${\mathit{u}}_{\mathrm{h}}$, $w$, $\dot{\mathit{\zeta }}$, $T_{\mathrm{v}}$, $p$, $q_i$	Yin–Yang (Sect. )	FD	C grid
ICON	NH (D)	${\mathit{u}}_{\mathrm{h}}$, $w$, $\mathit{\rho }$, ${\mathit{\theta }}_{\mathrm{v}}$, $\mathit{\rho }{q}_i$	Icosahedral triangular (Sect. )	FV	C grid
MPAS	NH	${\mathit{\rho }}_{\mathrm{d}}{\mathit{u}}_{\mathrm{h}}$, ${\mathit{\rho }}_{\mathrm{d}}w$, ${\mathit{\rho }}_{\mathrm{d}}$, ${\mathit{\rho }}_{\mathrm{d}}{\mathit{\theta }}_{\mathrm{v}}$, ${\mathit{\rho }}_{\mathrm{d}}{q}_i$	CCVT (Sect. )	FV	C grid
NICAM	NH	$\mathit{\rho }{\mathit{u}}_{\mathrm{h}}$, $\mathit{\rho }w$, $\mathit{\rho }$, $\mathit{\rho }e$, $\mathit{\rho }{q}_i$	Geodesic (Sect. )	FV	A grid
OLAM	NH (D)	$\mathit{\rho }{\mathit{u}}_{\mathrm{h}}$, $\mathit{\rho }w$, $\mathit{\rho }$, $\mathit{\rho }{\mathit{\theta }}_{\mathrm{il}}$, $\mathit{\rho }{q}_i$	Geodesic (Sect. )	FV	C grid
Tempest	NH	${\mathit{u}}_{\mathrm{h}}$, $w$, $\mathit{\rho }$, $\mathit{\rho }{\mathit{\theta }}_{\mathrm{v}}$, $\mathit{\rho }{q}_i$	Cubed sphere (Sect. )	SE	A grid

Vertical staggering (detailed in Sect. ) and vertical coordinates (detailed in Sect. ) for participating dynamical cores.

Acronym	Vertical staggering	Vertical coordinate
ACME-A	Co-located	Floating mass (Sect. )
CSU	Lorenz	Fixed height
DYNAMICO	Lorenz	Floating mass (Sect. )
FV${}^{\mathrm{3}}$	Co-located	Floating mass (Sect. )
FVM	Co-located	Fixed height
GEM	Modified Charney–Phillips (Sect. )	Log pressure (Sect. )
ICON	Lorenz	Fixed height
MPAS	Lorenz	Fixed height
NICAM	Lorenz	Fixed height
OLAM	Lorenz	Fixed height with cut cells (Sect. )
Tempest	Lorenz	Fixed height

Principal options for diffusion, stabilization, filters, or fixers in participating dynamical cores (detailed in Sect. ) and temporal discretization (detailed in Sect. ).

Acronym	Principal options for diffusion,	Temporal discretization
	stabilization, filters, or fixers
ACME-A	Fourth-order horizontal hyperviscosity	KGU53
CSU	Fourth-order horizontal hyperviscosity	Third-order Adams–Bashforth (AB3)
DYNAMICO	Fourth-order horizontal hyperviscosity	ARK(2,3,2)
FV${}^{\mathrm{3}}$	Divergence damping, hyperviscosity	Forward–backward /semi-implicit
FVM	Monotonic limiting	Semi-implicit (Sect. )
GEM	Hyperviscosity	Semi-implicit (Sect. )
ICON	Divergence damping, Smagorinsky, hyperdiffusion	Predictor–corrector
MPAS	Smagorinsky, hyperdiffusion	Split-explicit
NICAM	3-D divergence damping, Smagorinsky, hyperviscosity	Split-explicit
OLAM	Divergence/vorticity damping	Second-order Adams–Bashforth, Lax–Wendroff (for tracers)
Tempest	Fourth-order horizontal hyperviscosity	ARS(2,3,2)

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) as both share a common origin in the High Order Method Modeling Environment (HOMME) . 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. ), although unstructured, regionally 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. ). Vertical operators are based on the mimetic (mass- and energy-conserving) second-order finite difference discretization of . 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 finite-volume model using an optimized geodesic grid (Sect. ), with height as the vertical coordinate. The model is based on the non-hydrostatic unified system of equations proposed by , 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 (${\mathrm{\nabla }}_z^{\mathrm{4}}$) that acts on the vorticity, divergence, potential temperature, and tracer (Sect. ). 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. ) and a floating vertical mass coordinate (Sect. ). 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 . It includes a novel Hamiltonian formulation of the equations of motion in non-Eulerian coordinates which is imitated at the discrete level using building blocks from the literature 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 hyperviscosity operator (Sect. ). In addition, it features a conservative positive-definite transport scheme based on a slope-limited finite-volume approach .

FV cubed (FV${}^{\mathrm{3}}$)

The GFDL Finite-Volume Cubed-Sphere Dynamical Core (FV${}^{\mathrm{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. ) 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.  below). The non-hydrostatic formulation extends the hydrostatic model described in 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 (see also Sect. ), 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 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. ). 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.

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 . FVM solves the non-hydrostatic Euler equations on an octahedral reduced Gaussian grid (Sect. ) with a height-based terrain-following vertical coordinate . 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 rotational modes (Sect. ). The associated 3-D Helmholtz problem is solved iteratively using a bespoke preconditioned generalized conjugate residual approach. The integration procedure uses the non-oscillatory, finite-volume MPDATA (multidimensional positive definite advection transport algorithm) advection scheme . The non-oscillatory (i.e., monotonic) MPDATA also provides sufficient dissipation/diffusion to stabilize the model, so no other explicit filtering mechanism is required (Sect. ). 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 .

Global Environmental Multiscale (GEM) model

The Global Environmental Multiscale (GEM) model is used for operational forecasting at Environment and Climate Change Canada. GEM solves the non-hydrostatic Euler equations on the Yin–Yang grid (Sect. ) with Arakawa C-grid staggering of prognostic variables. The vertical coordinate is a unique hybrid terrain-following coordinate of a log-hydrostatic-pressure type (Sect. ) and the vertical discretization is based on the Charney–Phillips grid (Sect. ). A two-time-level, semi-Lagrangian implicit time discretization is implemented as described in Sect. . 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 Schwarz-type domain decomposition method . 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. ).

ICOsahedral Non-hydrostatic (ICON) model

The ICOsahedral Non-hydrostatic (ICON) model is a finite-volume model that solves the non-hydrostatic Euler equations in 2-D vector-invariant form on an icosahedral (triangular) grid (Sect. ) 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. ). 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) 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. ) and the mimetic TRiSK discretization . 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 . The prognostic variables are dry air pseudo-density (${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{d}}$), dry momentum (${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{d}}\mathit{u}$), and a modified moist potential temperature. Integration in time is handled via the split-explicit method of . 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. ).

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. ) optimized with spring dynamics using the method of . A terrain-following height coordinate system is used in the vertical with Lorenz staggering. Instead of temperature or potential temperature, total energy is prognosed following the method of . 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 . 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 is available and recommended. The split-explicit time discretization is used for the horizontally propagating sound waves with the 3-D divergence damping term (Sect. ). An implicit time discretization is adopted for the vertically propagating wave modes. A variant of the piecewise linear transport scheme is used with a flux limiter of for passive tracer transports.

Ocean–Land–Atmosphere Model (OLAM)

Ocean–Land–Atmosphere Model (OLAM) is a finite-volume model that solves the deep-atmosphere non-hydrostatic Euler equations in momentum conservation form on a hexagonal Voronoi mesh (Sect. ) 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. ). 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 , with consistent fluxes obtained by time averaging over the acoustic time steps.

Tempest

The Tempest model is an experimental test bed for high-performance numerical methods that solves the non-hydrostatic Euler equations on a cubed-sphere grid (Sect. ) 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 . Fourth-order hyperviscosity is used in the horizontal to prevent a buildup of energy at the grid scale (Sect. ). 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.

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 , and the arrangement of prognostic and diagnostic variables around each grid region . 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. , 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 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 . This grid is presently employed in many global models, including the UK Met Office New Dynamics and ENDGame dynamical cores . The latitude–longitude grid is also an option in the GEM model.

Cubed-sphere grid

The equiangular, gnomonic cubed-sphere grid 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 . On the equiangular cubed-sphere grid, coordinates are given as $(\mathit{\alpha },\mathit{\beta },p)$, with central angles $\mathit{\alpha },\mathit{\beta }\in [-\frac{\mathit{\pi }}{\mathrm{4}},\frac{\mathit{\pi }}{\mathrm{4}}]$ and panel index $p$. The structure of this grid supports refinement through stretching or nesting . 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. a. Among the DCMIP2016 models, the cubed-sphere grid is employed by the ACME, FV${}^{\mathrm{3}}$, and Tempest dynamical cores.

(a) A cubed-sphere grid. (b) An icosahedral (triangular) grid with additional refinement over Europe, as indicated in red. (c) An icosahedral (hexagonal) grid.

[Figure omitted. See PDF]

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^{\mathrm{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=\mathrm{2}$, 3, or 5. All additional $m$ refinement steps use $n=\mathrm{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_{\mathrm{c}}=\mathrm{20}\cdot n^{\mathrm{2}}\cdot {\mathrm{2}}^{\mathrm{2}m}$ cells, $n_{\mathrm{e}}=\mathrm{3}/\mathrm{2}\cdot n_{\mathrm{c}}$ edges, and $n_{\mathrm{v}}=\mathrm{10}\cdot n^{\mathrm{2}}\cdot {\mathrm{2}}^{\mathrm{2}m}+\mathrm{2}$ vertices. The anisotropy of global grids is reduced by the spring dynamics of . An example of such a grid is depicted in Fig. b. A discussion of the effective resolution of such grids is given in . 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 $R\mathrm{3}B\mathrm{7}$ global grid with 2 949 120 cells and 13 $\mathrm{km}$ mean resolution in combination with a refined region over Europe at 6.5 $\mathrm{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. ), 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. c. 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 , 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 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. ), 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.

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.

[Figure omitted. See PDF]

Constrained centroidal Voronoi tessellation (CCVT) meshes

Given a set of $N$ distinct points on the sphere $x_i$ (referred to as the generators, $\mathrm{1}\le i\le N$), the Voronoi tessellation (or the Voronoi diagram) associated with the generators is the set of polygons ${\mathrm{\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$ . 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) , 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) . Figure  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. , 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.

A constrained centroidal Voronoi tessellation mesh with localized grid density that could be employed in the MPAS model.

[Figure omitted. See PDF]

Octahedral reduced Gaussian grid

As with the classical reduced Gaussian grid of , the octahedral reduced Gaussian grid 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  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. ).

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 .

[Figure omitted. See PDF]

Yin–Yang grid

The overset Yin–Yang grid has two Cartesian grid components (subsets of a latitude–longitude grid) which are geometrically identical (see Fig. ). These components are combined to cover a spherical surface with partial overlap along their borders. The Yin component covers the latitude–longitude region $${\displaystyle }{\displaystyle }\left(-,{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}},-,{\mathit{\delta }}_{\mathit{\theta }},\le ,\mathit{\theta },\le ,{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}},+,{\mathit{\delta }}_{\mathit{\theta }}\right)\cap \left(-,{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{4}}},-,{\mathit{\delta }}_{\mathit{\lambda }},\le ,\mathit{\lambda },\le ,{\displaystyle \frac{\mathrm{3}\mathit{\pi }}{\mathrm{4}}},+,{\mathit{\delta }}_{\mathit{\lambda }}\right),$$ where ${\mathit{\delta }}_{\mathit{\lambda }},{\mathit{\delta }}_{\mathit{\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 ${\mathit{\delta }}_{\mathit{\theta }}=\mathrm{2}$${}^{\circ }$ and ${\mathit{\delta }}_{\mathit{\lambda }}=\mathrm{3}{\mathit{\delta }}_{\mathit{\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.

The Yin–Yang grid is a combination of two limited-domain latitude–longitude grids assembled to provide complete coverage of the sphere.

[Figure omitted. See PDF]

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 . The original four Arakawa grids , 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 , and the ZM grid, which extends the B grid to hexagons by placing the velocity at hexagonal nodes . 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 . 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. ).

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 $\mathrm{2}\mathrm{\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.

Horizontal staggering options represented among DCMIP models, in this case depicted on a rectilinear grid and geodesic grid. Here, $\mathit{\eta }$ denotes the buoyancy variable.

[Figure omitted. See PDF]

Other specialized staggerings have been developed that couple horizontal staggering with the formulation of the time integrator. In the FV${}^{\mathrm{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 .

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.

(a) A Lorenz-type variable staggering for a model utilizing height coordinates, (b) a Charney–Phillips-type variable staggering for a model utilizing height coordinates, and (c) a modified Charney–Phillips-type staggering used in the GEM model that introduces a new near-surface level for vertical velocity and temperature.

[Figure omitted. See PDF]

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 . The choice of vertical staggering may also impact many physically relevant properties of the model near the grid scale, such as the phase speed of Rossby waves . Finally, the choice of vertical staggering can have impacts on the physics–dynamics coupling . 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 , which places horizontal velocity, buoyancy, and thermodynamic variables on model levels, and vertical velocity on model interfaces; or (b) a Charney–Phillips-type staggering , which places horizontal velocity and buoyancy variables on model levels and vertical velocity and thermodynamic variables on model interfaces (see Fig. ). These approaches can be further augmented as needed, for instance, by shifting the vertical velocity and thermodynamic variables from the bottom boundary to an intermediate level, as in the GEM model. Note that, in general, tracer variables are co-located with the buoyancy variable.

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,\mathit{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 , which is in terms of the altitude $z$ and takes the form $${\displaystyle }s(\mathit{x},z)=z_{\mathrm{top}}\left[{\displaystyle \frac{z-z_{\mathrm{s}}\left(\mathit{x}\right)}{z_{\mathrm{top}}-z_{\mathrm{s}}\left(\mathit{x}\right)}}\right],$$ where $\mathit{x}$ denotes the horizontal position, $z_{\mathrm{s}}\left(\mathit{x}\right)$ is the height of the topography at that position, and $z_{\mathrm{top}}$ denotes the height of the model top (typically independent of position). Analogous formulations are available for mass-based ($\mathit{\sigma }$ 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 or . 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., ). 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 are a generalization of pressure-based coordinates to non-hydrostatic models, with a vertical coordinate defined as the total gravity-weighted overhead mass, $${\displaystyle }s=\underset{z}{\overset{\mathrm{\infty }}{\int }}\mathit{\rho }g\mathrm{d}z.$$ Under this definition, $${\displaystyle \frac{\partial s}{\partial z}}=-\mathit{\rho }g.$$

GEM $\mathit{\zeta }$ coordinate

The vertical coordinate in the GEM model, denoted $\mathit{\zeta }$, is a hybrid terrain-following coordinate of a log-hydrostatic-pressure type. Taking $s$ (denoted $\mathit{\pi }$ in GEM documentation) as given in Eq. (), $\mathit{\zeta }$ is given by the relation $${\displaystyle }\log s=A\left(\mathit{\zeta }\right)+B\left(\mathit{\zeta }\right)\left[\log ,s,\left(,z_{\mathrm{s}},\right),-,{\mathit{\zeta }}_{\mathrm{s}}\right],$$ with $${\displaystyle }A\left(\mathit{\zeta }\right)=\mathit{\zeta },\text{and}B\left(\mathit{\zeta }\right)={\left({\displaystyle \frac{\mathit{\zeta }-{\mathit{\zeta }}_{\mathrm{top}}}{{\mathit{\zeta }}_{\mathrm{s}}-{\mathit{\zeta }}_{\mathrm{top}}}}\right)}^r.$$ Here, ${\mathit{\zeta }}_{\mathrm{s}}=\log \left({\mathrm{10}}^{\mathrm{5}}\right)$, ${\mathit{\zeta }}_{\mathrm{top}}=\log \left(s_{\mathrm{top}}\right)$, $s_{\mathrm{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${}^{\mathrm{3}}$)

In the floating Lagrangian formulation , 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., $${\displaystyle }s=\underset{z}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{\mathrm{0}}g\mathrm{d}z,$$ in which case a 3-D reference density field ${\mathit{\rho }}_{\mathrm{0}}$ can be imposed. Of primary importance is the fact that the vertical coordinate satisfies $${\displaystyle }\dot{s}={\displaystyle \frac{\mathrm{d}s}{\mathrm{d}t}}=\mathrm{0},$$ 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 in order to completely avoid topographic imprinting on the model grid levels (Fig. ). 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 is used to stabilize small cut cells.

(a) A terrain-following coordinate passing over rough topography. (b) A cut-cell coordinate used for representing the same topography.

[Figure omitted. See PDF]

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 . 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 , and the ability of the model to conserve important invariants such as energy . 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 shallow-atmosphere equations using a hybrid terrain-following pressure vertical coordinate $\mathit{\eta }$, similar to the model of . The 2-D vector-invariant form of the prognostic horizontal velocity Eq. () is employed, in conjunction with prognostic potential temperature (Eq. ), pseudo-density (Eq. ), and geopotential (Eq. ). The vertical velocity equation is formulated analogous to that of GEM: $${\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}}{\displaystyle }=-g_{\mathrm{c}}\left(\mathrm{1},-,{\displaystyle \frac{\partial p}{\partial s}}\right).$$

CSU

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

DYNAMICO

The prognostic equations employed by DYNAMICO are based on a Hamiltonian formulation . The specific prognostic variables employed are pseudo-density ${\mathit{\rho }}_{\mathrm{s}}$, mass-weighted tracers (potential temperature, water species), geopotential $\mathrm{\Phi }$, horizontal covariant components of momentum, and mass-weighted vertical momentum $W={\mathit{\rho }}_{\mathrm{s}}{g}^{-\mathrm{2}}d\mathrm{\Phi }/\mathrm{d}t={\mathit{\rho }}_{\mathrm{s}}{g}^{-\mathrm{1}}w$. Prognostic equations are in flux form for mass (Eq. ) and $W$ (Eq. ), in advective form for $\mathrm{\Phi }$ (Eq. ), and in vector-invariant form for covariant horizontal momentum (Eq. ).

FV${}^{\mathrm{3}}$

The hydrostatic FV${}^{\mathrm{3}}$ model uses a mass-based floating Lagrangian coordinate along with the shallow-atmosphere approximation . Prognostic equations include horizontal velocity in 2-D vector-invariant form (Eq. ), pseudo-density (Eq. ), and virtual potential temperature (Eq. ). The non-hydrostatic model further incorporates prognostic geopotential (Eq. ) and vertical momentum (Eq. ).

FVM

The FVM formulation is based on conservation laws for dry mass (Eq. ), momentum (Eq. ), and dry entropy (Eq. ) in Eulerian flux form, which are similar to Eq. () for ${\mathit{\rho }}_{\mathrm{d}}$, Eq. (), and Eq. () for $\mathit{\theta }$, 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 , the FVM governing equations can concisely be written as $$\begin{array}{cc}& {\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}}{\partial t}}+\stackrel{\mathrm{‾}}{\mathrm{\nabla }}\cdot \left(\mathit{v},\mathcal{G},{\mathit{\rho }}_{\mathrm{d}}\right)=\mathrm{0},{\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\mathit{u}}{\partial t}}+\stackrel{\mathrm{‾}}{\mathrm{\nabla }}\cdot \left(\mathit{v},\mathcal{G},{\mathit{\rho }}_{\mathrm{d}},\mathit{u}\right)\\ {\displaystyle }& {\displaystyle }=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\left(-,{\mathit{\theta }}_{\mathit{\rho }},\stackrel{\mathrm{̃}}{\mathbf{G}},\stackrel{\mathrm{‾}}{\mathrm{\nabla }},{\mathit{\varphi }}^{\prime },-,\mathit{k},,g,\left({\displaystyle \frac{{\mathit{\theta }}^{\prime }}{{\mathit{\theta }}_{\mathrm{a}}}},+,{\mathit{\epsilon }}_b,,q_{\mathrm{v}}^{\prime },-,q_{\mathrm{c}},-,q_{\mathrm{r}}\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..-,\mathrm{2},\mathbf{\Omega },\times ,\left(\mathit{u},-,{\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}}{{\mathit{\theta }}_{\mathit{\rho }a}}},{\mathit{u}}_a\right),+,{\mathcal{M}}^{\prime }\right],& {\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\theta }}^{\prime }}{\partial t}}+\stackrel{\mathrm{‾}}{\mathrm{\nabla }}\cdot \left(\mathit{v},\mathcal{G},{\mathit{\rho }}_{\mathrm{d}},,{\mathit{\theta }}^{\prime }\right)=-\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^T\mathit{u}\cdot \stackrel{\mathrm{‾}}{\mathrm{\nabla }}{\mathit{\theta }}_a,& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}^{\prime }=c_{\mathrm{pd}}\left[{\left({\displaystyle \frac{R_{\mathrm{d}}}{p_{\mathrm{0}}}},{\mathit{\rho }}_{\mathrm{d}},,\mathit{\theta },,(,\mathrm{1},+,q_{\mathrm{v}},/,\mathit{\epsilon },)\right)}^{R_{\mathrm{d}}/c_{\mathrm{vd}}},-,{\mathit{\pi }}_a\right].\end{array}$$ Dependent variables in Eq. () are dry density ${\mathit{\rho }}_{\mathrm{d}}$, 3-D physical velocity vector $\mathit{u}$, potential temperature perturbation ${\mathit{\theta }}^{\prime }$, and a modified Exner pressure perturbation ${\mathit{\varphi }}^{\prime }$, with the thermodynamic variables related by the gas law Eq. (). Primes indicate perturbations with respect to the prescribed ambient state denoted by subscript “a”; see and for discussions. The symbol $g$ in Eq. () denotes the gravitational acceleration and ${\mathit{\epsilon }}_b=\mathrm{1}/\mathit{\epsilon }-\mathrm{1}$. As far as geometric aspects are concerned, the nabla operator $\stackrel{\mathrm{‾}}{\mathrm{\nabla }}$ represents the 3-D vector of partial derivatives with respect to the curvilinear coordinates, along with the Jacobian $\mathcal{G}$, a matrix of metric coefficients $\stackrel{\mathrm{̃}}{\mathbf{G}}$, its transpose ${\stackrel{\mathrm{̃}}{\mathbf{G}}}^T$, and the contravariant velocity $\mathit{v}={\stackrel{\mathrm{̃}}{\mathbf{G}}}^T\mathit{u}$, where a contribution from optional time dependency of the curvilinear coordinates is neglected for simplicity . The symbol ${\mathcal{M}}^{\prime }={\mathcal{M}}^{\prime }(\mathit{u},{\mathit{u}}_{\mathrm{a}},{\mathit{\theta }}_{\mathit{\rho }}/{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}})$ in Eq. () 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 $$\mathit{\mu }={\displaystyle \frac{\partial p}{\partial s}}-\mathrm{1},$$ where $s$ (denoted $\mathit{\pi }$ in GEM documentation) is given by Eq. (). In this case, the equations of GEM model are concisely given by $$\begin{array}{}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{\mathit{u}}_{\mathrm{h}}}{\mathrm{d}t}}+f\mathit{k}\times {\mathit{u}}_{\mathrm{h}}+R_{\mathrm{d}}{T}_{\mathrm{v}}{\mathrm{\nabla }}_{\mathit{\zeta }}\log p+\left(\mathrm{1},+,\mathit{\mu }\right){\mathrm{\nabla }}_{\mathit{\zeta }}\mathrm{\Phi }=\mathrm{0},& {\displaystyle }& {\displaystyle \frac{\mathrm{d}w}{\mathrm{d}t}}-g_{\mathrm{c}}\mathit{\mu }=\mathrm{0},& {\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\log \left({\displaystyle \frac{\partial s}{\partial \mathit{\zeta }}}\right)+{\mathrm{\nabla }}_{\mathit{\zeta }}\cdot {\mathit{u}}_{\mathrm{h}}+{\displaystyle \frac{\partial \dot{\mathit{\zeta }}}{\partial \mathit{\zeta }}}=\mathrm{0},& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\log T_{\mathrm{v}}}{\mathrm{d}t}}-{\displaystyle \frac{R_{\mathrm{d}}}{c_{\mathrm{p}}}}{\displaystyle \frac{\mathrm{d}\log p}{\mathrm{d}t}}=\mathrm{0},& {\displaystyle }& {\displaystyle \frac{\partial \mathrm{\Phi }}{\partial s}}+{\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}}{p}}=\mathrm{0},& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}t}}-g_{\mathrm{c}}w=\mathrm{0}.\end{array}$$ Here, ${\mathrm{\nabla }}_{\mathit{\zeta }}$ denotes the horizontal gradient along $\mathit{\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 $\mathit{\rho }$, affecting the virtual temperature defined in Eq. (). Also, specific mass is used in GEM (not mixing ratio).

ICON

ICON solves a non-hydrostatic equation set based on 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 , the barycentric (bc) velocity ${\mathit{u}}_{\mathrm{bc}}={\sum }_k{\mathit{\rho }}_k{\mathit{u}}_k/{\sum }_k{\mathit{\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 , a vector-invariant form is only used for the horizontal velocity Eq. (), whereas the vertical velocity equation is solved in advective form. The pressure-gradient force is formulated according to Eq. ().

Additional prognostic variables include total air density (Eq. ), virtual potential temperature (Eq. ), and mass fractions $q_k={\mathit{\rho }}_k/\mathit{\rho }$ of all constituents (except for dry air) for which the prognostic equation reads $${\displaystyle \frac{\partial \mathit{\rho }{q}_k}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho },q_k,{\mathit{u}}_{\mathrm{bc}}\right)=-\mathrm{\nabla }\cdot {\mathbf{J}}_k+{\mathit{\sigma }}_k,$$ with ${\mathit{\sigma }}_k$ describing sources/sinks due to phase changes, and ${\mathbf{J}}_k=\mathit{\rho }{q}_k\left({\mathit{u}}_k,-,{\mathit{u}}_{\mathrm{bc}}\right)$ denoting diffusion fluxes, which account for the motion of constituents relative to the frame of reference set by ${\mathit{u}}_{\mathrm{bc}}$.

The specific heat capacities and ideal gas constant are approximated to be equal to their dry values $R^{\ast }\approx R_{\mathrm{d}}$, $c_{\mathrm{p}}^{\ast }\approx c_{\mathrm{pd}}$, $c_{\mathrm{v}}^{\ast }\approx c_{\mathrm{vd}}$. The model also uses a prognostic equation for Exner pressure to simplify the treatment of vertical sound wave propagation, given by $${\displaystyle \frac{\partial \mathit{\pi }}{\partial t}}+{\displaystyle \frac{R_{\mathrm{d}}}{c_{\mathrm{vd}}}}{\displaystyle \frac{\mathit{\pi }}{\mathit{\rho }{\mathit{\theta }}_{\mathrm{v}}}}\mathrm{\nabla }\cdot \left({\mathit{u}}_{\mathrm{bc}},\mathit{\rho },{\mathit{\theta }}_{\mathrm{v}}\right)=\widehat{Q},$$ where $\widehat{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_{\mathrm{n}}$, making use of reconstructed tangential velocity components $v_{\mathrm{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. () by setting the vertical component of ${\mathit{u}}_{\mathrm{bc}}$ to zero at the lower boundary: $w_{\mathrm{bc}}{|}_{\mathrm{sfc}}=\mathrm{0}$. In addition, only the vertical diffusion fluxes ${\mathbf{J}}_k$ of sedimentation constituents and the surface evaporation flux ${\mathbf{J}}_{\mathrm{v}}{|}_{\mathrm{sfc}}$ are taken into account. The counter flux of non-sedimentation constituents is discarded. Since in the given framework the continuity Eq. () is only valid if the constraint ${\sum }_k{\mathbf{J}}_k=\mathrm{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 , based on the formulation of . The MPAS model uses the momentum form of the update equations, as described in Sect. , with dry mass utilized for the density variable ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{\mathrm{s}}$. MPAS further evolves dry mass using a continuity equation of Eq. () and moist potential temperature following Eq. ().

NICAM

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

OLAM

OLAM solves the deep-atmosphere, fully compressible equations in mass- and momentum-conserving finite-volume form using Eqs. (), (), and (). Prognostic variables are the three components of momentum, ice–liquid potential temperature ${\mathit{\theta }}_{\mathrm{il}}$ , total density $\mathit{\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 $\mathit{\theta }$. The latter is in turn diagnosed from ${\mathit{\theta }}_{\mathrm{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 , 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.

Tempest

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

Diffusion, stabilization, filters, and fixers

Most dynamical cores implement specialized techniques for diffusion or stabilization (see Table ). 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 . 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 $\mathit{\rho }$, $\mathit{\theta }$, and tracer variables via repeated application of a scalar Laplacian . Vector hyperviscosity is also applied by decomposing the horizontal vector Laplacian into divergence damping and vorticity damping terms via the vector identity $${\displaystyle }{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{u}}_{\mathrm{h}}=\mathrm{\nabla }\mathrm{\nabla }\cdot {\mathit{u}}_{\mathrm{h}}+\mathrm{\nabla }\times \mathrm{\nabla }\times {\mathit{u}}_{\mathrm{h}}.$$ 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 .

CSU

The CSU model utilizes an explicit diffusion scheme that consists of fourth-order hyperdiffusion (${\mathrm{\nabla }}^{\mathrm{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_{\mathrm{diff}}$ Runge–Kutta time step in a forward-Euler manner, with $N_{\mathrm{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 $\mathit{\tau }$: given $\mathit{\tau }$, the hyperviscous coefficient that multiplies operator $D^p$ is ${\mathit{\delta }}^{\mathrm{2}p}{\mathit{\tau }}^{-}\mathrm{1}$, where ${\mathit{\delta }}^{-}\mathrm{2}$ is the largest eigenvalue of operator $D$. For DCMIP, DYNAMICO uses $p=\mathrm{2}$ (fourth-order hyperviscosity) for all filters.

FV${}^{\mathrm{3}}$

Explicit dissipation in FV${}^{\mathrm{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${}^{\mathrm{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${}^{\mathrm{3}}$ is constructed with a flexible-lid (constant-pressure) upper boundary that is effective at damping internal gravity wave modes; however, FV${}^{\mathrm{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 cut-off pressure.

FV${}^{\mathrm{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 Exner-function 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_{\mathrm{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 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 and a global mass fixer . Those approaches have been evaluated using chemical constituents such as ozone .

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 , and are summarized here. For damping, in the corrector step, a fourth-order divergence damping term $F_{\mathrm{d}}\left(\mathit{v}\right)$ 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: $${\displaystyle }{F}_{\mathrm{d}}\left(\mathit{v}\right)=-f_{\mathrm{d}}{\stackrel{\mathrm{‾}}{a_{\mathrm{c}}}}^{\mathrm{2}}\mathrm{\nabla }\stackrel{\mathrm{̃}}{\mathrm{\nabla }}\cdot \left\{\mathrm{\nabla },\left[\stackrel{\mathrm{̃}}{\mathrm{\nabla }},\cdot ,v,+,{\displaystyle \frac{\mathrm{\Delta }}{\mathrm{\Delta }z}},\left(w,-,{\stackrel{\mathrm{‾}}{{\stackrel{\mathrm{‾}}{w_{cc}}}^c}}^i\right)\right]\right\}.$$ $f_{\mathrm{d}}$ typically attains values between $\frac{\mathrm{1}}{\mathrm{1000}\mathrm{\Delta }t}$ and $\frac{\mathrm{1}}{\mathrm{250}\mathrm{\Delta }t}$, and $\stackrel{\mathrm{‾}}{a_{\mathrm{c}}}$ is the global mean cell area.

ICON also includes Rayleigh damping on $w$ following , 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_{\mathrm{D}\mathrm{2}}\left(v_n\right)$) and potential temperature ($F_{\mathrm{D}\mathrm{2}}\left(\mathit{\theta }\right)$) combined with a fourth-order background diffusion of velocity $F_{\mathrm{D}\mathrm{4}}\left(v_n\right)$, defined via $$\begin{array}{cccc}{\displaystyle }{F}_{\mathrm{D}\mathrm{2}}\left(v_n\right)& {\displaystyle }=\mathrm{4}{K}_{\mathrm{h}}{\stackrel{\mathrm{̃}}{\mathrm{\nabla }}}^{\mathrm{2}}\left(v_n\right),& {\displaystyle }{F}_{\mathrm{D}\mathrm{2}}\left(\mathit{\theta }\right)& {\displaystyle }=a_{\mathrm{c}}\stackrel{\mathrm{̃}}{\mathrm{\nabla }}\cdot \left(K_{\mathrm{h}},{\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{\mathrm{\Delta }n}}\right),\\ & {\displaystyle }{F}_{\mathrm{D}\mathrm{4}}\left(v_n\right)& {\displaystyle }=-k_{\mathrm{4}}{a}_{\mathrm{e}}^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathrm{\nabla }}}^{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathrm{\nabla }}}^{\mathrm{2}}\right(v_n\left)\right),\end{array}$$ where $a_{\mathrm{c}}$ and $a_{\mathrm{e}}$ denote the area associated with the cell and edge under consideration, respectively. An empirically determined offset of $\mathrm{0.75}{k}_{\mathrm{4}}{a}_{\mathrm{e}}$ is subtracted from $K_{\mathrm{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 $\mathrm{km}$) or finer because the advection of vertical wind speed has no implicit damping of small-scale structures. This term appears as $${\displaystyle }{F}_{\mathrm{D}}\left(w\right)=-k_w{a}_{\mathrm{c}}^{\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}\left({\mathrm{\nabla }}^{\mathrm{2}}\right(w\left)\right).$$

MPAS

The MPAS model applies fourth-order hyperdiffusion and Smagorinsky diffusion , as described in . When applied to the momentum, the Laplacian is evaluated as $${\displaystyle }{\mathrm{\nabla }}^{\mathrm{2}}{u}_i={\displaystyle \frac{\partial }{\partial x_i}}{\mathrm{\nabla }}_s\cdot \mathit{v}-{\displaystyle \frac{\partial \mathit{\eta }}{\partial x_j}},$$ where $u_i$ is the edge-normal velocity defined on cell edge $i$, $\mathit{\eta }$ is the vertical component of the relative vorticity, computed on vertices, and ${\mathrm{\nabla }}_s\cdot \mathit{v}$ is the horizontal divergence on $s$ surfaces, computed on edges. The evaluation of divergence and vorticity in this expression is described in . 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_{\mathrm{h}}$, defined in terms of Cartesian velocities $(u,v)$: $$K_{\mathrm{h}}=c_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\ell }}^{\mathrm{2}}\sqrt{(u_x-v_y{)}^{\mathrm{2}}+(u_y+v_x{)}^{\mathrm{2}}},$$ where $c_{\mathrm{s}}$ is a constant parameter and $\mathrm{\ell }$ is the grid scale. The diffusion operator then takes the form $\mathrm{\nabla }\cdot \left(K_{\mathrm{h}}\mathrm{\nabla }\mathit{\psi }\right)$ for a scalar field $\mathit{\psi }$.

NICAM

NICAM implements three types of diffusion: 3-D divergence damping, fourth-order horizontal hyperdiffusion, and sixth-order vertical hyperdiffusion, as described in . Specifically, the divergence damping term 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 and .

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 $\mathrm{km}$ of the domain, where the model top is 35 or 40 $\mathrm{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 . 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 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 $\mathrm{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 , quasi-hydrostatic , pseudo-incompressible , or unified approximation . 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 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, ). 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 , 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 .

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${}^{\mathrm{3}}$ 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 . The horizontal component of the pressure-gradient force is evaluated backward in time using the algorithm of ; 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 .

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 .

MPAS and NICAM use a split-explicit formulation 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 responsible 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 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.

The FVM semi-implicit method

A characteristic feature of the FVM (Sect. ) 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 and for dry dynamics, whereas it can be found in and for extension to moist-precipitating dynamics. Here, we provide a short outline of the solution procedure for the compressible Euler Eq. (). It employs the two-time-level, second-order-accurate template algorithm given as $$\begin{array}{cc}{\displaystyle }{\mathit{\psi }}_i^{n+\mathrm{1}}& {\displaystyle }={\mathcal{A}}_i({\stackrel{\mathrm{̃}}{\mathit{\psi }}}^n,{\mathit{v}}^{n+\mathrm{1}/\mathrm{2}},G^n,G^{n+\mathrm{1}})+\mathrm{0.5}\mathrm{\Delta }t{R}^{\mathit{\psi }}{|}_i^{n+\mathrm{1}},\\ & {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\psi }}}^n& {\displaystyle }\equiv {\mathit{\psi }}^n+\mathrm{0.5}\mathrm{\Delta }t{R}^{\mathit{\psi }}{|}^n,\end{array}$$ where $\mathit{\psi }$ represents the solution variable, $R^{\mathit{\psi }}$ is the respective right-hand side, $\mathcal{A}$ symbolizes the advective transport operator given by the non-oscillatory, finite-volume MPDATA scheme , 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. () can then be divided into three steps. First, the homogenous mass continuity Eq. () is integrated with $\mathit{\psi }\equiv {\mathit{\rho }}_{\mathrm{d}}$, $\mathit{v}\equiv \mathit{v}\mathcal{G}$, $G\equiv \mathcal{G}$, and $R^{{\mathit{\rho }}_{\mathrm{d}}}\equiv \mathrm{0}$ in Eq. (). Second, given already updated moisture variables , the thermodynamic (Eq. ) and momentum (Eq. ) equations enter Eq. () with $\mathit{\psi }=u,v,w,{\mathit{\theta }}^{\prime }$, $\mathit{v}\equiv \mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}$, $G\equiv \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}$, and the right-hand-side $R^{\mathit{\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. ) 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 ${\mathit{\varphi }}^{\prime }$ using an advective form of the equation of state (Eq. ). An $\mathcal{O}\left(\mathrm{\Delta }{t}^{\mathrm{2}}\right)$ integration of this equation with a Euler backward scheme, in the spirit of Eq. (), leads to a Helmholtz equation . The associated 3-D elliptic boundary value problem is solved iteratively using a preconditioned generalized conjugate residual approach; see for a recent overview and comprehensive list of references. Non-linear terms in $R^{\mathit{\psi }}{|}^{n+\mathrm{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 $${\displaystyle \frac{\mathrm{d}F}{\mathrm{d}t}}+G=\mathrm{0},$$ where $\mathrm{d}/\mathrm{d}t$ 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. (): $${\displaystyle \frac{F^{\mathrm{A}}-F^{\mathrm{D}}}{\mathrm{\Delta }t}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},+,\mathit{ϵ}\right)G^{\mathrm{A}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},-,\mathit{ϵ}\right)G^{\mathrm{D}}=\mathrm{0},$$ where “A” stands for the arrival position at model grid point $({\mathbf{r}}_{\mathrm{h}},\mathit{\zeta },t)$ and “D” for the departure position $({\mathbf{r}}_{\mathrm{h}}-\mathrm{\Delta }{\mathbf{r}}_{\mathrm{h}},\mathit{\zeta }-\mathrm{\Delta }\mathit{\zeta },t-\mathrm{\Delta }t)$ due to the displacements $\mathrm{\Delta }{\mathbf{r}}_{\mathrm{h}},\mathrm{\Delta }\mathit{\zeta }$ having occurred during the time step $\mathrm{\Delta }t$. $G$ is averaged between these two positions with a possible slight off-centering $\mathit{ϵ}$. The displacements are themselves calculated by solving, again using the Lagrangian method, the following equations: $${\displaystyle \frac{\mathrm{d}{\mathbf{r}}_{\mathrm{h}}}{\mathrm{d}t}}-{\mathit{u}}_{\mathrm{h}}=\mathrm{0};{\displaystyle \frac{\mathrm{d}\mathit{\zeta }}{\mathrm{d}t}}-\dot{\mathit{\zeta }}=\mathrm{0},$$ discretized in the same way (trapezoidal method) though without off-centering: $${\displaystyle \frac{\mathrm{\Delta }{\mathbf{r}}_{\mathrm{h}}}{\mathrm{\Delta }t}}-{\displaystyle \frac{{{\mathit{u}}_{\mathrm{h}}}^{\mathrm{A}}+{{\mathit{u}}_{\mathrm{h}}}^{\mathrm{D}}}{\mathrm{2}}}=\mathrm{0};{\displaystyle \frac{\mathrm{\Delta }\mathit{\zeta }}{\mathrm{\Delta }t}}-{\displaystyle \frac{{\dot{\mathit{\zeta }}}^{\mathrm{A}}+{\dot{\mathit{\zeta }}}^{\mathrm{D}}}{\mathrm{2}}}=\mathrm{0}.$$ The process is of course a multi-step iterative one, since both positions and velocities at departure positions (past time $t-\mathrm{\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. () is rewritten as follows: $${\displaystyle \frac{F^{\mathrm{A}}}{\mathit{\tau }}}+G^{\mathrm{A}}={\displaystyle \frac{F^{\mathrm{D}}}{\mathit{\tau }}}-\mathit{\beta }{G}^{\mathrm{D}}\equiv R^{\mathrm{D}},$$ where $\mathit{\tau }=\left(\mathrm{1},/,\mathrm{2},+,\mathit{ϵ}\right)\mathrm{\Delta }t$ and $\mathit{\beta }=\left(\mathrm{1},/,\mathrm{2},-,\mathit{ϵ}\right)/\left(\mathrm{1},/,\mathrm{2},+,\mathit{ϵ}\right)$. Secondly, separating linear from non-linear parts, we get $${\displaystyle }{L}^{\mathrm{A}}+N^{\mathrm{A}}=R^{\mathrm{D}},$$ with $$\begin{array}{cc}{\displaystyle }{L}^{\mathrm{A}}& {\displaystyle }={\left[{\displaystyle \frac{F^{\mathrm{A}}}{\mathit{\tau }}},+,G^{\mathrm{A}}\right]}_{\mathrm{linear}},\text{ and }\\ & {\displaystyle }{N}^{\mathrm{A}}& {\displaystyle }\equiv {\displaystyle \frac{F^{\mathrm{A}}}{\mathit{\tau }}}+G^{\mathrm{A}}-{\left[{\displaystyle \frac{F^{\mathrm{A}}}{\mathit{\tau }}},+,G^{\mathrm{A}}\right]}_{\mathrm{linear}}.\end{array}$$ Note that both $F$ and $G$ may require linearization. $L^{\mathrm{A}}$ may then be obtained if $N^{\mathrm{A}}$ is first guessed; once $L^{\mathrm{A}}$ is found, an estimate of $N^{\mathrm{A}}$ is obtained and $L^{\mathrm{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 . 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 . 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 , , , , , , , , and .

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 , conservation of invariants , or representation of the kinetic energy spectrum . 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.

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

Short Name	Code availability
ACME-A	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.${}^{†}$
CSU	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).

Short Name	Code availability
DYNAMICO	DYNAMICO is open source and available online from IPSL Forge (http://forge.ipsl.jussieu.fr/dynamico) or directly by request to Thomas Dubos ([email protected]). The release used for DCMIP2016 can be found via Zenodo (http://dx.doi.org/10.5281/zenodo.583718).
FV${}^{\mathrm{3}}$	FV${}^{\mathrm{3}}$ 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.
FVM	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”.${}^{\mathrm{a}}$
GEM	Due to licensing requirements, GEM model code is only available by request to Abdessamad Qaddouri ([email protected]) or Vivian Lee ([email protected]).
ICON	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 ([email protected]) or Marco Giorgetta ([email protected]).
MPAS	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).
NICAM	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.
OLAM	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://doi.org/10.5281/zenodo.582308).

Short name	Code availability
Tempest	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).

${}^{\mathrm{a}}$ In compliance with the GMD editorial requirements, this code has been made available to the topical editor in charge of the paper.

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 . The material derivative is used for quantities in the Lagrangian frame (following individual air parcels) and is given by $${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}={\displaystyle \frac{\partial }{\partial t}}+\mathit{u}\cdot \mathrm{\nabla },$$ where $\mathit{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 $${\displaystyle \frac{\mathrm{d}{q}_i}{\mathrm{d}t}}=\mathrm{0}.$$

A list of physical constants used in this document.

Constant	Description	Value
$a_{\mathrm{ref}}$	Radius of the Earth	$\mathrm{6.37122}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{m}$
${\mathrm{\Omega }}_{\mathrm{ref}}$	Rotational speed of the Earth	$\mathrm{7.292}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{1}}$
$g_{\mathrm{c}}$	Gravitational acceleration	$\mathrm{9.80616}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}$
$p_{\mathrm{0}}$	Reference pressure	$\mathrm{1000}\mathrm{hPa}$
$c_{\mathrm{pd}}$	Specific heat capacity of dry air at constant pressure	$\mathrm{1004.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$c_{\mathrm{pv}}$	Specific heat capacity of water vapor at constant pressure	$\mathrm{1930.0}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$c_{\mathrm{vd}}$	Specific heat capacity of dry air at constant volume	$\mathrm{717.5}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$c_{\mathrm{vv}}$	Specific heat capacity of water vapor at constant volume	$\mathrm{1460.0}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$R_{\mathrm{d}}$	Gas constant for dry air	$\mathrm{287.0}\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$R_{\mathrm{v}}$	Gas constant for water vapor	$\mathrm{461.5}$ $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$
$\mathit{\epsilon }$	Ratio of $R_{\mathrm{d}}$ to $R_{\mathrm{v}}$	$\mathrm{0.622}$
$M_{\mathrm{v}}$	Constant for virtual temperature conversion	$\mathrm{0.608}$
${\mathit{\rho }}_{\mathrm{water}}$	Reference density of water	1000 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$

Diagnostic relationships

The atmospheric fluid is assumed to be an ideal gas. For moist air, the ideal gas constant $R^{\ast }$, specific heat capacity at constant pressure $c_{\mathrm{p}}^{\ast }$, and specific heat capacity at constant volume $c_{\mathrm{v}}^{\ast }$ are given by $$\begin{array}{cccc}{\displaystyle }{R}^{\ast }& {\displaystyle }=R_{\mathrm{d}}+(R_{\mathrm{v}}-R_{\mathrm{d}})q,& {\displaystyle }{c}_{\mathrm{p}}^{\ast }=& {\displaystyle }{c}_{\mathrm{pd}}+(c_{\mathrm{pv}}-c_{\mathrm{pd}})q,\\ & {\displaystyle }{c}_{\mathrm{v}}^{\ast }& {\displaystyle }=c_{\mathrm{vd}}+(c_{\mathrm{vv}}-c_{\mathrm{vd}})q.\end{array}$$ Note that in many models, $R^{\ast }$, $c_{\mathrm{p}}^{\ast }$, and $c_{\mathrm{v}}^{\ast }$ are approximated by $R_{\mathrm{d}}$, $c_{\mathrm{pd}}$, and $c_{\mathrm{vd}}$, respectively. Dry air, water vapor, and moist air quantities all satisfy the linear relationship $R=c_{\mathrm{p}}-c_{\mathrm{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 ${\mathit{\rho }}_{\mathrm{d}}$, moist density $\mathit{\rho }$, pressure $p$, vapor pressure $e$, temperature $T$, virtual temperature $T_{\mathrm{v}}$, Exner pressure $\mathit{\pi }$, potential temperature $\mathit{\theta }$, and virtual potential temperature ${\mathit{\theta }}_{\mathrm{v}}$. Common ratios $\mathit{\kappa }=R^{\ast }/c_{\mathrm{p}}^{\ast }$, $\mathit{\epsilon }=R_{\mathrm{d}}/R_{\mathrm{v}}$, and $\mathit{\gamma }=c_{\mathrm{p}}^{\ast }/c_{\mathrm{v}}^{\ast }$ 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 ${\mathit{\theta }}_{\mathit{\rho }}$.

Relationships between key thermodynamic variables arise from the ideal gas law, along with definitions of Exner pressure, potential temperature, and virtual potential temperature: $$\begin{array}{cccccc}{\displaystyle }p=& {\displaystyle }\mathit{\rho }{R}_{\mathrm{d}}{T}_{\mathrm{v}},& {\displaystyle }\mathit{\pi }=& {\displaystyle }{\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^{\mathit{\kappa }},& {\displaystyle }\mathit{\theta }=& {\displaystyle }T{\left({\displaystyle \frac{p_{\mathrm{0}}}{p}}\right)}^{\mathit{\kappa }},\\ & {\displaystyle }{\mathit{\theta }}_{\mathrm{v}}=& {\displaystyle }{T}_{\mathrm{v}}{\left({\displaystyle \frac{p_{\mathrm{0}}}{p}}\right)}^{\mathit{\kappa }},\end{array}$$ which further give rise to $${\displaystyle }p={\displaystyle }{\left({\displaystyle \frac{\mathit{\rho }{R}_{\mathrm{d}}{\mathit{\theta }}_{\mathrm{v}}}{p_{\mathrm{0}}^{\mathit{\kappa }}}}\right)}^{\mathit{\gamma }},{\displaystyle }\mathit{\pi }={\displaystyle }{\left({\displaystyle \frac{\mathit{\rho }{R}_{\mathrm{d}}{\mathit{\theta }}_{\mathrm{v}}}{p_{\mathrm{0}}}}\right)}^{R^{\ast }/c_{\mathrm{v}}^{\ast }},{\displaystyle }\mathit{\theta }={\displaystyle \frac{T}{\mathit{\pi }}},{\displaystyle }{\mathit{\theta }}_{\mathrm{v}}={\displaystyle \frac{T_{\mathrm{v}}}{\mathit{\pi }}}.$$ Note that virtual temperature is typically written as $${\displaystyle }{T}_{\mathrm{v}}={\displaystyle }T\left(\mathrm{1},+,{\displaystyle \frac{(\mathrm{1}-\mathit{\epsilon })}{\mathit{\epsilon }}},q\right),$$ which arises from the relationship $${\displaystyle }{T}_{\mathrm{v}}={\displaystyle \frac{T}{\mathrm{1}-\frac{e}{p}(\mathrm{1}-\mathit{\epsilon })}},$$ upon applying $e/p=q/\mathit{\epsilon }$.

Prognostic equations for thermodynamic variables

Note that, as a consequence of Eq. (), the following simplifications can be applied: $${\displaystyle \frac{\mathrm{1}}{T_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{v}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{T}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}},{\displaystyle \frac{\mathrm{d}{R}^{\ast }}{\mathrm{d}t}}={\displaystyle }\mathrm{0},{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{p}}^{\ast }}{\mathrm{d}t}}={\displaystyle }\mathrm{0},{\displaystyle \frac{\mathrm{d}{c}_{\mathrm{v}}^{\ast }}{\mathrm{d}t}}={\displaystyle }\mathrm{0}.$$ Mass conservation is typically represented through the continuity equation, which can be written in the Lagrangian frame as $${\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}={\displaystyle }-\mathit{\rho }\mathrm{\nabla }\cdot \mathit{u},$$ or equivalently in the Eulerian frame: $${\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}={\displaystyle }-\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{u}\right).$$ Further prognostic relationships can be derived from the thermodynamic equation, including the diabatic heating rate $J$: $${\displaystyle \frac{\mathrm{1}}{T}}{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}-{\displaystyle \frac{\mathit{\kappa }}{p}}{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}={\displaystyle \frac{J}{T{c}_{\mathrm{p}}^{\ast }}},$$ which can be alternatively written as $${\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}={\displaystyle \frac{J\mathit{\theta }}{T{c}_{\mathrm{p}}^{\ast }}},\text{ or }{\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{v}}}{\mathrm{d}t}}={\displaystyle \frac{J{\mathit{\theta }}_{\mathrm{v}}}{T{c}_{\mathrm{p}}^{\ast }}}.$$ These equations can then be combined with Eq. () to obtain $${\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }{\mathit{\theta }}_{\mathrm{v}}\right)+\mathrm{\nabla }\cdot \left(\mathit{\rho }{\mathit{\theta }}_{\mathrm{v}}\mathit{u}\right)={\displaystyle \frac{J\mathit{\rho }{\mathit{\theta }}_{\mathrm{v}}}{T{c}_{\mathrm{p}}^{\ast }}},$$ or similarly for $\mathit{\theta }$. In conjunction with the material derivative of the ideal gas law, $${\displaystyle \frac{\mathrm{1}}{p}}{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{v}}}{\mathrm{d}t}},$$ the thermodynamic equation can be written in the form $${\displaystyle \frac{c_{\mathrm{v}}^{\ast }}{R^{\ast }{T}_{\mathrm{v}}}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{v}}}{\mathrm{d}t}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}={\displaystyle \frac{J}{T{R}^{\ast }}}.$$ Then, substituting Eq. () gives a prognostic equation for virtual temperature: $${\displaystyle \frac{c_{\mathrm{v}}^{\ast }}{R^{\ast }}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{v}}}{\mathrm{d}t}}+T_{\mathrm{v}}\mathrm{\nabla }\cdot \mathit{u}={\displaystyle \frac{J{T}_{\mathrm{v}}}{T{R}^{\ast }}}.$$ The prognostic equation for temperature is identical except with $T$ substituted for $T_{\mathrm{v}}$. An analogous equations for pressure can be obtained through a similar procedure: $${\displaystyle \frac{c_{\mathrm{v}}^{\ast }}{c_{\mathrm{p}}^{\ast }}}{\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}+p\mathrm{\nabla }\cdot \mathit{u}={\displaystyle \frac{Jp}{T{c}_{\mathrm{p}}^{\ast }}},$$ and similarly for Exner pressure: $${\displaystyle \frac{c_{\mathrm{v}}^{\ast }}{R^{\ast }}}{\displaystyle \frac{\mathrm{d}\mathit{\pi }}{\mathrm{d}t}}+\mathit{\pi }\mathrm{\nabla }\cdot \mathit{u}={\displaystyle \frac{J\mathit{\pi }}{T{c}_{\mathrm{p}}^{\ast }}}.$$

Momentum equations

In coordinate-invariant form, the prognostic velocity equations may be written in either the Lagrangian or Eulerian frame as $${\displaystyle \frac{\mathrm{d}\mathit{u}}{\mathrm{d}t}}={\displaystyle \frac{\partial \mathit{u}}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }\mathit{u}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }p-\mathrm{2}\mathbf{\Omega }\times \mathit{u}-\mathrm{\nabla }\mathrm{\Phi },$$ where $\mathbf{\Omega }$ denotes the planetary vorticity vector and $\mathrm{\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 $\mathit{u}\cdot \mathrm{\nabla }\mathit{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 $${\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }p=-c_{\mathrm{p}}^{\ast }\mathit{\theta }\left[\mathrm{\nabla },\mathit{\pi },-,\mathit{\pi },\ln ,\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right),\mathrm{\nabla },\mathit{\kappa }\right],$$ which follows from Eq. (). Note that often in non-hydrostatic models, $\mathit{\kappa }$ is assumed constant and the $\mathrm{\nabla }\mathit{\kappa }$ term neglected. A second form of Eq. () emerges on substituting the vector calculus identity $${\displaystyle }\mathit{u}\cdot \mathrm{\nabla }\mathit{u}={\displaystyle }\mathrm{\nabla }K+\mathit{\zeta }\times \mathit{u},$$ where $K=\frac{\mathrm{1}}{\mathrm{2}}(\mathit{u}\cdot \mathit{u})$ is the 3-D specific kinetic energy and $\mathit{\zeta }=\mathrm{\nabla }\times \mathit{u}$ is the 3-D relative vorticity vector. This gives rise to the 3-D vector-invariant form, $${\displaystyle \frac{\partial \mathit{u}}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }p-\mathrm{\nabla }(K+\mathrm{\Phi })-(\mathit{\zeta }+\mathrm{2}\mathbf{\Omega })\times \mathit{u}.$$ Because no gradients of vectors appear in this equation, it avoids derivatives of the coordinate basis that would arise from the momentum transport term $\mathit{u}\cdot \mathrm{\nabla }\mathit{u}$ in Eq. (). In conjunction with Eq. (), Eq. () also gives rise to the flux-form momentum equations, $${\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }\mathit{u}\right)=-\mathrm{\nabla }\cdot (\mathit{u}\otimes \mathit{u}+\mathcal{I}p)-\mathrm{2}\mathbf{\Omega }\times \left(\mathit{\rho }\mathit{u}\right)-\mathit{\rho }\mathrm{\nabla }\mathrm{\Phi },$$ where $\mathit{u}\otimes \mathit{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 $\mathrm{\Phi }$ and $\mathbf{\Omega }$. For deep-atmosphere models, one typically chooses $${\displaystyle }\mathrm{\Phi }=g_{\mathrm{c}}{a}^{\mathrm{2}}\left[{\displaystyle \frac{\mathrm{1}}{a}},-,{\displaystyle \frac{\mathrm{1}}{a+z}}\right],\text{ and }\mathbf{\Omega }=\mathrm{\Omega }(\mathit{k}\sin \mathit{\phi }+\mathit{j}\cos \mathit{\phi }),$$ where $g_{\mathrm{c}}$ is gravitational acceleration at the surface, $a$ is the radius of the planet, $\mathrm{\Omega }$ is the rotation rate (in s${}^{-\mathrm{1}}$), $\mathit{\phi }$ is the latitude, $\mathit{j}$ is the unit vector oriented in the meridional direction, and $\mathit{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=\mathrm{d}z/\mathrm{d}t$, $${\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}t}}={\displaystyle \frac{a^{\mathrm{2}}{g}_{\mathrm{c}}w}{(a+z{)}^{\mathrm{2}}}}.$$ For shallow-atmosphere models, the geopotential takes the simpler form $${\displaystyle }\mathrm{\Phi }=g_{\mathrm{c}}z,\text{ and }\mathbf{\Omega }=\mathrm{\Omega }\sin \mathit{\phi }\mathit{k},$$ where $z$ is the altitude above the surface. In this case, we write $\mathrm{2}\mathbf{\Omega }=f\mathit{k}$, where $f=\mathrm{2}\mathrm{\Omega }\sin \mathit{\phi }$ is the Coriolis parameter. The evolution equation for the shallow-atmosphere geopotential is then $${\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}t}}=g_{\mathrm{c}}w.$$

Orthogonal formulation

Under the orthogonal formulation, projection of a vector field $\mathit{b}$ onto its horizontal components is defined via $${\displaystyle }[\mathit{b}{]}_z=\mathit{b}-(\mathit{b}\cdot \mathit{k})\mathit{k}.$$ When applied to the velocity vector, this gives rise to the decomposition $${\displaystyle }\mathit{u}={\mathit{u}}_{\mathrm{h}}+w\mathit{k},$$ where $\mathit{k}=\mathrm{\nabla }z$ is the unit vector in the vertical direction and ${\mathit{u}}_{\mathrm{h}}=[\mathit{u}{]}_z$ (${\mathit{u}}_{\mathrm{h}}$ is aligned with surfaces of constant $z$). In the orthogonal formulation, the material derivative expands as $${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}={\displaystyle \frac{\partial }{\partial t}}+{\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }+w\mathit{k}\cdot \mathrm{\nabla }.$$ For the special case of the material derivative applied to scalars, this equation can also be written as $${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}={\displaystyle \frac{\partial }{\partial t}}+{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_z+w{\displaystyle \frac{\partial }{\partial z}}.$$ where ${\mathrm{\nabla }}_{z}b=[\mathrm{\nabla }b{]}_z$ denotes the gradient along constant $z$ surfaces. From here, the vector-invariant form velocity equation obtained by multiplying Eq. () by $\mathit{k}$ expands as $${\displaystyle \frac{\partial w}{\partial t}}={\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial z}}-{\displaystyle \frac{\partial }{\partial z}}(K+\mathrm{\Phi })-\left[(,\mathit{\zeta },+,\mathrm{2},\mathbf{\Omega },),\times ,\mathit{u}\right]\cdot \mathit{k},$$ which, from ${\mathit{u}}_{\mathrm{h}}=\mathit{u}-w\mathit{k}$, then gives rise to $${\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}={\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{z}p-{\mathrm{\nabla }}_z(K+\mathrm{\Phi })-{\left[(,\mathit{\zeta },+,\mathrm{2},\mathbf{\Omega },),\times ,\mathit{u}\right]}_z.$$ Due to its association with hydrostatic models, it is common to use the 2-D kinetic energy, $K_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}({\mathit{u}}_{\mathrm{h}}\cdot {\mathit{u}}_{\mathrm{h}})$. Decomposing the momentum transport term into horizontal and vertical components gives $${\displaystyle }\mathit{u}\cdot \mathrm{\nabla }\mathit{u}={\displaystyle }{\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{h}}+(\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{h}})\times \left(w\mathit{k}\right)+(\mathit{u}\cdot \mathrm{\nabla }w)\mathit{k}.$$ The first term in this expression admits the relationships $$\begin{array}{cc}& {\displaystyle }{\left[{\mathit{u}}_{\mathrm{h}},\cdot ,\mathrm{\nabla },{\mathit{u}}_{\mathrm{h}}\right]}_z& {\displaystyle }={\mathrm{\nabla }}_z{K}_{\mathrm{2}}+{\mathit{\zeta }}_{\mathrm{h}}\mathit{k}\times {\mathit{u}}_{\mathrm{h}},{\displaystyle }({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }{\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}& {\displaystyle }=-{\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})\\ {\displaystyle }& {\displaystyle }=-K_{\mathrm{2}}(\mathrm{\nabla }\cdot \mathit{k})-\frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{u}}_{\mathrm{h}},\cdot ,(,\mathrm{\nabla },\times ,{\mathit{u}}_t,)\right.\\ & {\displaystyle }& {\displaystyle }\left..+,(,\mathrm{\nabla },\times ,{\mathit{u}}_{\mathrm{h}},),\cdot ,{\mathit{u}}_t\right],\end{array}$$ where ${\mathit{\zeta }}_{\mathrm{h}}=(\mathrm{\nabla }\times \mathit{u})\cdot \mathit{k}=(\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}$ is the relative vorticity scalar and ${\mathit{u}}_t=\mathit{k}\times {\mathit{u}}_{\mathrm{h}}$. Note that this equation does incorporate metric terms associated with horizontal advection of $\mathit{k}$, which must be accounted for.

Thus, the vertical velocity equation, obtained by taking Eq. ()$\cdot \mathit{k}$, is $$\begin{array}{cc}{\displaystyle \frac{\partial w}{\partial t}}=& {\displaystyle }{\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-w{\displaystyle \frac{\partial w}{\partial z}}-{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{z}w-{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial z}}\\ & {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial z}}-(\mathrm{2}\mathbf{\Omega }\times {\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}.\end{array}$$ Then, subtracting Eq. () $\cdot \mathit{k}$ from Eq. () gives $$\begin{array}{cc}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}& {\displaystyle }=-w({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-w{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial z}}-{\mathrm{\nabla }}_z(K_{\mathrm{2}}+\mathrm{\Phi })\\ & {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{z}p-{\mathit{\zeta }}_{\mathrm{h}}\mathit{k}\times {\mathit{u}}_{\mathrm{h}}-[\mathrm{2}\mathbf{\Omega }\times \mathit{u}{]}_z.\end{array}$$ Note that, under the shallow-atmosphere approximation, the metric term ${\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})$ in Eq. () is set equal to zero in accordance with .

Arbitrary vertical coordinates

The dynamical equations are now formulated in terms of the vertical coordinate $s(t,\mathit{x},z)$ with $\partial s/\partial z\ne \mathrm{0}$ everywhere, i.e., following (hereafter K74). Since $\mathit{x}$ and $t$ are shared between the two coordinate systems, the chain rule can be applied to obtain expressions $$\begin{array}{cccc}{\displaystyle \frac{\partial }{\partial z}}=& {\displaystyle \frac{\partial s}{\partial z}}{\displaystyle \frac{\partial }{\partial s}},& {\displaystyle }{\mathrm{\nabla }}_s=& {\displaystyle }{\mathrm{\nabla }}_z+\left({\mathrm{\nabla }}_{s}z\right)(\mathit{k}\cdot \mathrm{\nabla }),\\ & {\displaystyle }{\left({\displaystyle \frac{\partial }{\partial t}}\right)}_s=& {\displaystyle \frac{\partial }{\partial t}}+{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s(\mathit{k}\cdot \mathrm{\nabla }),\end{array}$$ 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 $${\displaystyle \frac{\partial s}{\partial z}}={\displaystyle }{\left({\displaystyle \frac{\partial z}{\partial s}}\right)}^{-\mathrm{1}},{\displaystyle }{\mathrm{\nabla }}_{z}s={\displaystyle }-\left({\displaystyle \frac{\partial s}{\partial z}}\right){\mathrm{\nabla }}_{s}z,{\displaystyle \frac{\partial s}{\partial t}}={\displaystyle }-{\displaystyle \frac{\partial s}{\partial z}}{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s.$$ From here, Eq. () also gives rise to $${\displaystyle }{\mathrm{\nabla }}_z={\mathrm{\nabla }}_s-{\displaystyle \frac{\partial s}{\partial z}}\left({\mathrm{\nabla }}_{s}z\right){\displaystyle \frac{\partial }{\partial s}},$$ which can be used directly to rewrite Eq. () or Eq. () in terms of derivatives over $s$. Note that the operators ${\mathrm{\nabla }}_z$ and ${\mathrm{\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 $\mathit{k}\cdot {\mathrm{\nabla }}_{z}A=\mathrm{0}$ and $\mathit{k}\cdot {\mathrm{\nabla }}_{s}A=\mathrm{0}$ for any scalar field $A$.

From Eq. (), it can be shown that the 2-D divergence on $s$ surfaces (given by K74 Eq. 3.17) is $${\displaystyle }{\mathrm{\nabla }}_s\cdot {\mathit{u}}_{\mathrm{h}}={\mathrm{\nabla }}_z\cdot {\mathit{u}}_{\mathrm{h}}+\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left({\mathrm{\nabla }}_{s}z\right)\cdot \left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\right),$$ and that the 2-D curl is given by $${\displaystyle }{\mathrm{\nabla }}_s\times {\mathit{u}}_{\mathrm{h}}={\mathrm{\nabla }}_z\times {\mathit{u}}_{\mathrm{h}}+\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left({\mathrm{\nabla }}_{s}z\right)\times \left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\right),$$ where ${\mathrm{\nabla }}_z\times {\mathit{u}}_{\mathrm{h}}=\mathit{k}(\mathit{k}\cdot (\mathrm{\nabla }\times {\mathit{u}}_{\mathrm{h}}\left)\right)$. Notably, these expressions are valid for both shallow- and deep-atmosphere formulations.

The generalized velocity $\dot{s}$ following a fluid parcel is defined by $${\displaystyle }\dot{s}\equiv {\displaystyle \frac{\mathrm{d}s}{\mathrm{d}t}}={\displaystyle \frac{\partial s}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }s={\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{z}s+\left[w,-,{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\right]{\displaystyle \frac{\partial s}{\partial z}}.$$ Then, using Eqs. () and () to rewrite Eq. () gives an expression for the material derivative for scalars on $s$ surfaces: $$\begin{array}{cc}{\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=& {\displaystyle }{\left({\displaystyle \frac{\partial A}{\partial t}}\right)}_s-{\displaystyle \frac{\partial s}{\partial z}}{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s{\displaystyle \frac{\partial A}{\partial s}}\\ & {\displaystyle }+& {\displaystyle }{\mathit{u}}_{\mathrm{h}}\cdot \left[{\mathrm{\nabla }}_s,A,+,\left(,{\mathrm{\nabla }}_z,s,\right),{\displaystyle \frac{\partial A}{\partial s}}\right]+w{\displaystyle \frac{\partial s}{\partial z}}{\displaystyle \frac{\partial A}{\partial s}}& {\displaystyle }=& {\displaystyle }{\left({\displaystyle \frac{\partial A}{\partial t}}\right)}_s+{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{s}A+\dot{s}{\displaystyle \frac{\partial A}{\partial s}}.\end{array}$$ A similar expression arises for vectors, although in this case ${\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathbf{a}\ne {\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_z\mathbf{a}$ implies we cannot use the operator ${\mathrm{\nabla }}_s$ in the form (), and instead obtain $${\displaystyle \frac{\mathrm{d}\mathbf{a}}{\mathrm{d}t}}={\displaystyle }{\left({\displaystyle \frac{\partial \mathbf{a}}{\partial t}}\right)}_s+\left[{\mathit{u}}_{\mathrm{h}},\cdot ,\mathrm{\nabla },\mathbf{a},+,(,{\mathit{u}}_{\mathrm{h}},\cdot ,{\mathrm{\nabla }}_s,z,),(,\mathit{k},\cdot ,\mathrm{\nabla },\mathbf{a},)\right]+\dot{s}{\displaystyle \frac{\partial \mathbf{a}}{\partial s}}.$$

Conservation laws in arbitrary vertical coordinates

Using Eq. (), we observe that the 3-D divergence on the sphere takes the form $${\displaystyle }\mathrm{\nabla }\cdot \mathit{u}={\mathrm{\nabla }}_z\cdot {\mathit{u}}_{\mathrm{h}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\alpha }w\right),$$ where $\mathit{\alpha }=\mathrm{1}$ for shallow-atmosphere models and $\mathit{\alpha }=r^{\mathrm{2}}=(a+z{)}^{\mathrm{2}}$ for deep-atmosphere models. Using $w=\mathrm{d}z/\mathrm{d}t$, this last term also takes the form $${\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\alpha }w\right)={\displaystyle \frac{\partial w}{\partial z}}+{\displaystyle \frac{w}{\mathit{\alpha }}}{\displaystyle \frac{\partial \mathit{\alpha }}{\partial z}}={\displaystyle \frac{\partial w}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{\displaystyle \frac{\mathrm{d}\mathit{\alpha }}{\mathrm{d}t}}.$$ Using Eq. () to rewrite Eq. () gives rise to $${\displaystyle }w={\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s+{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{s}z+\dot{s}{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},$$ which is then differentiated to yield K74 Eq. (3.16), $${\displaystyle \frac{\partial w}{\partial z}}=\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left[{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}},{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},+,\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\right),\cdot ,\left(,{\mathrm{\nabla }}_s,z,\right)\right]+{\displaystyle \frac{\partial \dot{s}}{\partial s}}=\mathrm{0}.$$ Substituting this expression into continuity Eq. () and using Eqs. (), (), and () then leads to $$\begin{array}{cc}{\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[\mathit{\alpha },{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},\mathit{\rho }\right]+\mathit{\alpha }{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}}\mathit{\rho }[{\mathrm{\nabla }}_z\cdot {\mathit{u}}_{\mathrm{h}}\\ & {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\right)\cdot \left({\mathrm{\nabla }}_{s}z\right)]+\mathit{\alpha }{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}}\mathit{\rho }{\displaystyle \frac{\partial \dot{s}}{\partial s}}=\mathrm{0}.\end{array}$$ Defining the pseudo-density as $${\displaystyle }{\mathit{\rho }}_{\mathrm{s}}=\mathit{\alpha }{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}}\mathit{\rho },$$ and using Eq. () in the form $${\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{s}}}{\mathrm{d}t}}={\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{s}}}{\partial t}}\right)}_s+{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_s{\mathit{\rho }}_{\mathrm{s}}+\dot{s}{\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{s}}}{\partial s}},$$ along with Eq. (), leads to $${\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{s}}}{\partial t}}\right)}_s+{\mathrm{\nabla }}_s\cdot \left({\mathit{\rho }}_{\mathrm{s}}{\mathit{u}}_{\mathrm{h}}\right)+{\displaystyle \frac{\partial }{\partial s}}\left({\mathit{\rho }}_{\mathrm{s}}\dot{s}\right)=\mathrm{0}.$$ Hence, for any quantity that is conserved following a fluid parcel (i.e., $\mathrm{d}q/\mathrm{d}t=\mathrm{0}$), $${\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{s}}q}{\partial t}}\right)}_s+{\mathrm{\nabla }}_s\cdot \left({\mathit{\rho }}_{\mathrm{s}}q{\mathit{u}}_{\mathrm{h}}\right)+{\displaystyle \frac{\partial }{\partial s}}\left({\mathit{\rho }}_{\mathrm{s}}q\dot{s}\right)=\mathrm{0}.$$ In particular, the prognostic equation for virtual potential temperature (or equivalently for potential temperature) reads $${\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}}{\partial t}}\right)}_s+{\mathrm{\nabla }}_s\cdot \left({\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}{\mathit{u}}_{\mathrm{h}}\right)+{\displaystyle \frac{\partial }{\partial s}}\left({\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}\dot{s}\right)={\displaystyle \frac{J{\mathit{\rho }}_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{v}}}{c_{\mathrm{p}}^{\ast }}}.$$

2-D vector-invariant form

The prognostic equations utilizing horizontal kinetic energy $K_{\mathrm{2}}$ in place of $K$ are derived by applying Eq. () to Eq. (), yielding $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial w}{\partial t}}\right)}_s=& {\displaystyle }{\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-{\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }w\\ {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left\{\left[{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,-,w\right],{\displaystyle \frac{\partial w}{\partial s}}\right.\\ & {\displaystyle }& {\displaystyle }\left..-,{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial s}},-,{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\displaystyle \frac{\partial p}{\partial s}}\right\}-(\mathrm{2}\mathbf{\Omega }\times {\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}.\end{array}$$ Similarly, from Eq. (), $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }-w({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})+\left[{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,-,w\right]\left({\displaystyle \frac{\partial s}{\partial z}}\right){\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\\ & {\displaystyle }& {\displaystyle }-{\mathit{\zeta }}_{\mathrm{h}}\mathit{k}\times {\mathit{u}}_{\mathrm{h}}-{\mathrm{\nabla }}_z(K_{\mathrm{2}}+\mathrm{\Phi })-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{z}p-[\mathrm{2}\mathbf{\Omega }\times \mathit{u}{]}_z.\end{array}$$ 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. () into Eq. () then gives $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial w}{\partial t}}\right)}_s=& {\displaystyle }{\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-{\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{s}w-\dot{s}{\displaystyle \frac{\partial w}{\partial s}}\\ & {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left[-,{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial s}},-,{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\displaystyle \frac{\partial p}{\partial s}}\right]-(\mathrm{2}\mathbf{\Omega }\times {\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}.\end{array}$$ Similarly, substituting Eq. () into Eq. () and using the identity $$\begin{array}{cc}{\displaystyle }({\mathit{u}}_{\mathrm{h}}\cdot {\mathrm{\nabla }}_{s}z)\left({\displaystyle \frac{\partial s}{\partial z}}\right){\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}& {\displaystyle }=({\mathrm{\nabla }}_s\times {\mathit{u}}_{\mathrm{h}})\times {\mathit{u}}_{\mathrm{h}}\\ & {\displaystyle }& {\displaystyle }-{\mathit{\zeta }}_{\mathrm{h}}\mathit{k}\times {\mathit{u}}_{\mathrm{h}}+\left({\mathrm{\nabla }}_{s}z\right){\displaystyle \frac{\partial K_{\mathrm{2}}}{\partial z}}\end{array}$$ then gives $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }-w({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-{\mathrm{\nabla }}_s{K}_{\mathrm{2}}-{\mathit{\zeta }}_s\mathit{k}\times {\mathit{u}}_{\mathrm{h}}-\dot{s}{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}\\ & {\displaystyle }& {\displaystyle }-{\mathrm{\nabla }}_z\mathrm{\Phi }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{z}p-[\mathrm{2}\mathbf{\Omega }\times \mathit{u}{]}_z,\end{array}$$ where $${\displaystyle }{\mathrm{\nabla }}_s\times {\mathit{u}}_{\mathrm{h}}=\mathit{k}{\mathit{\zeta }}_s,\text{ and }{\mathit{\zeta }}_s=\mathit{k}\cdot ({\mathrm{\nabla }}_s\times {\mathit{u}}_{\mathrm{h}}).$$ In this case, the vertical advection terms are removed when $\dot{s}=\mathrm{0}$, i.e., the vertical coordinate is advected with the 3-D wind $\mathit{u}$.

Note that under the shallow-atmosphere approximation, the first metric terms (those that include $({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})$) in Eqs. ()–() are typically dropped.

3-D vector-invariant form

From Eqs. () and (), the evolution equation for the 3-D velocity vector takes the form $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial \mathit{u}}{\partial t}}\right)}_s& {\displaystyle }={\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s(\mathit{k}\cdot \mathrm{\nabla }\mathit{u})-\mathrm{\nabla }(K+\mathrm{\Phi })\\ & {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\mathrm{\nabla }p-(\mathit{\zeta }+\mathrm{2}\mathbf{\Omega })\times \mathit{u}.\end{array}$$ Then, taking the dot product of this expression with $\mathit{k}$ gives $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial w}{\partial t}}\right)}_s=& {\displaystyle }\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left[{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,{\displaystyle \frac{\partial w}{\partial s}},-,{\displaystyle \frac{\partial }{\partial s}},(,K,+,\mathrm{\Phi },),-,{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\displaystyle \frac{\partial p}{\partial s}}\right]\\ & {\displaystyle }& {\displaystyle }-\left[(,\mathit{\zeta },+,\mathrm{2},\mathbf{\Omega },),\times ,\mathit{u}\right]\cdot \mathit{k},\end{array}$$ where we have used $\mathit{k}\cdot (\mathit{k}\cdot \mathrm{\nabla }\mathit{u})=\mathit{k}\cdot \mathrm{\nabla }w$. Similarly, the prognostic equation for horizontal velocity from Eq. () is reformulated as $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\left({\displaystyle \frac{\partial s}{\partial z}}\right){\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}}-{\mathrm{\nabla }}_z(K+\mathrm{\Phi })\\ & {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{z}p-{\left[(,\mathit{\zeta },+,\mathrm{2},\mathbf{\Omega },),\times ,\mathit{u}\right]}_z.\end{array}$$ Note that the vorticity term in this expression can be simplified further using $${\displaystyle }\left[\right(\mathit{\zeta }+\mathrm{2}\mathbf{\Omega })\times \mathit{u}{]}_z=-({\mathit{\zeta }}_{\mathrm{h}}+\mathit{k}\cdot \mathrm{2}\mathbf{\Omega }\left)\right({\mathit{u}}_{\mathrm{h}}\times \mathit{k})-w\mathit{k}\times (\mathit{\zeta }+\mathrm{2}\mathbf{\Omega }),$$ and $${\displaystyle }-\mathit{k}\times \mathit{\zeta }={\displaystyle }\mathit{k}\cdot \mathrm{\nabla }\mathit{u}-\mathrm{\nabla }(\mathit{k}\cdot \mathit{u})+\mathit{u}\cdot \mathrm{\nabla }\mathit{k}={\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial z}}-{\mathrm{\nabla }}_{z}w+{\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k}.$$

Covariant component formulation

In conjunction with Eq. (), the horizontal momentum equation (in 2-D vector-invariant form as Eqs.  or , or in 3-D vector-invariant form as Eq. ) 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 $${\displaystyle }[\mathit{b}{]}_s\equiv [\mathit{b}{]}_z+\left({\mathrm{\nabla }}_{s}z\right)(\mathit{k}\cdot \mathit{b}).$$ Applying this operator to the horizontal velocity gives $${\displaystyle }{\mathit{v}}_{\mathrm{h}}\equiv [\mathit{u}{]}_s={\mathit{u}}_{\mathrm{h}}+({\mathrm{\nabla }}_{s}z)w.$$ For a time-dependent $s$ coordinate, we obtain the identity $${\displaystyle }{\left[{\displaystyle \frac{\partial }{\partial t}},\left(,{\mathrm{\nabla }}_s,z,\right)\right]}_s={\displaystyle }{\mathrm{\nabla }}_s{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,$$ and so can write $${\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial t}}\right)}_s={\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial t}}\right)}_s+\left({\mathrm{\nabla }}_{s}z\right){\left({\displaystyle \frac{\partial w}{\partial t}}\right)}_s+w{\mathrm{\nabla }}_s{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s.$$ Then, using Eqs. (), (), and () and identity $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left[{\displaystyle \frac{\partial {\mathit{u}}_{\mathrm{h}}}{\partial s}},+,\left(,{\mathrm{\nabla }}_s,z,\right),{\displaystyle \frac{\partial w}{\partial s}}\right]+w{\mathrm{\nabla }}_s{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \frac{\partial s}{\partial z}}\right){\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\left\{{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial s}},-,{\mathrm{\nabla }}_s,\left[{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},w\right]\right\}\\ & {\displaystyle }& {\displaystyle }+{\mathrm{\nabla }}_s\left[{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,w\right],\end{array}$$ gives $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }-{\mathrm{\nabla }}_s\left[K,-,w,{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,+,\mathrm{\Phi }\right]-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{s}p\\ & {\displaystyle }& {\displaystyle }-{\left[(,\mathit{\zeta },+,\mathrm{2},\mathbf{\Omega },),\times ,\mathit{u}\right]}_s& {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\partial s}{\partial z}}\right){\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s\left\{{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial s}},-,{\mathrm{\nabla }}_s,\left[{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},w\right]\right\}.\end{array}$$ Finally, we can expand the vorticity term and hence obtain $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }-{\mathrm{\nabla }}_s\left[K,-,w,{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s,+,\mathrm{\Phi }\right]-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\mathrm{\nabla }}_{s}p-{\left[\mathrm{2},\mathbf{\Omega },\times ,\mathit{u}\right]}_s\\ {\displaystyle }& {\displaystyle }+\left[\mathit{k},\cdot ,{\mathrm{\nabla }}_s,\times ,{\mathit{v}}_{\mathrm{h}}\right]({\mathit{u}}_{\mathrm{h}}\times \mathit{k})\\ & {\displaystyle }& {\displaystyle }-\dot{s}\left\{{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial s}},-,{\mathrm{\nabla }}_s,\left[{\left({\displaystyle \frac{\partial s}{\partial z}}\right)}^{-\mathrm{1}},w\right]\right\}-{\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s{\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k}.\end{array}$$

Vorticity divergence form

The vorticity divergence form of the dynamical equations in an arbitrary vertical coordinate predicts the absolute vorticity (${\mathit{\zeta }}_{\mathrm{h}}^{\ast }$) and velocity divergence ($D$) given by $${\displaystyle }{\mathit{\zeta }}_{\mathrm{h}}^{\ast }=({\mathrm{\nabla }}_s\times {\mathit{u}}_{\mathrm{h}}+\mathrm{2}\mathbf{\Omega })\cdot \mathit{k},$$ and $${\displaystyle }D\equiv {\mathrm{\nabla }}_s\cdot {\mathit{u}}_{\mathrm{h}},$$ respectively, instead of the horizontal velocity. The horizontal velocity can be obtained from the stream function $\mathit{\psi }$ and the velocity potential $\mathit{\chi }$ following $${\displaystyle }{\mathit{u}}_{\mathrm{h}}=\mathit{k}\times {\mathrm{\nabla }}_s\mathit{\psi }+{\mathrm{\nabla }}_s\mathit{\chi }.$$ By using Eq. () in Eqs. () and (), we obtain the elliptic equations that diagnose the stream function and velocity potential from the predicted velocity and divergence as $${\displaystyle }{\mathrm{\nabla }}_s^{\mathrm{2}}={\mathit{\zeta }}_{\mathrm{h}}^{\ast }-\mathrm{2}\mathbf{\Omega }\cdot \mathit{k},\text{ and }{\mathrm{\nabla }}_s^{\mathrm{2}}\mathit{\chi }=D,$$ respectively.

By taking the material derivative (Eq. ) of Eq. () and using horizontal momentum Eqs. (), (), and (), the absolute vorticity prediction equation emerges: $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\mathit{\zeta }}_{\mathrm{h}}^{\ast }}{\partial t}}\right)}_s& {\displaystyle }-J_s({\mathit{\zeta }}_{\mathrm{h}}^{\ast },\mathit{\psi })+{\mathrm{\nabla }}_s\cdot \left({\mathit{\zeta }}_{\mathrm{h}}^{\ast }{\mathrm{\nabla }}_s\mathit{\chi }\right)+{\mathrm{\nabla }}_s\cdot \left(\dot{s},{\displaystyle \frac{\partial }{\partial s}},{\mathrm{\nabla }}_s,\mathit{\psi }\right)\\ & {\displaystyle }& {\displaystyle }+\mathit{k}\cdot {\mathrm{\nabla }}_s\times \left(\dot{s},{\displaystyle \frac{\partial }{\partial s}},{\mathrm{\nabla }}_s,\mathit{\chi }\right)+J_s({\mathit{\rho }}^{-\mathrm{1}},p)=\mathrm{0},\end{array}$$ where $J_s(a,b)=\mathit{k}\cdot {\mathrm{\nabla }}_s\times \left(a{\mathrm{\nabla }}_{s}b\right)$ is the Jacobian operator. It can also be shown that $\dot{s}$ relates to the vertical velocity $w$ through $${\displaystyle }\dot{s}=\left({\displaystyle \frac{\partial s}{\partial z}}\right)(w-w_{\mathrm{c}}),$$ where $${\displaystyle }{w}_{\mathrm{c}}\equiv {\left({\displaystyle \frac{\partial z}{\partial t}}\right)}_s+\left(\mathit{k},\times ,{\mathrm{\nabla }}_s,\mathit{\psi },+,{\mathrm{\nabla }}_s,\mathit{\chi }\right)\cdot \left({\mathrm{\nabla }}_{s}z\right).$$

By taking the material derivative of Eq. () and using Eqs. (), (), and (), we can obtain the divergence prediction equation $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial D}{\partial t}}\right)}_s& {\displaystyle }-J_s({\mathit{\zeta }}_{\mathrm{h}}^{\ast },\mathit{\chi })-{\mathrm{\nabla }}_s\cdot \left({\mathit{\zeta }}_{\mathrm{h}}^{\ast }{\mathrm{\nabla }}_s\mathit{\psi }\right)+{\mathrm{\nabla }}_s\cdot \left(\dot{s},{\displaystyle \frac{\partial }{\partial s}},{\mathrm{\nabla }}_s,\mathit{\chi }\right)\\ & {\displaystyle }& {\displaystyle }+\left(\mathit{k},\times ,{\displaystyle \frac{\partial }{\partial s}},{\mathrm{\nabla }}_s,\mathit{\psi }\right)\cdot {\mathrm{\nabla }}_s\dot{s}& {\displaystyle }& {\displaystyle }+{\mathrm{\nabla }}_s\cdot ({\mathrm{\nabla }}_s{K}_{\mathrm{2}}+g{\mathrm{\nabla }}_{s}z)+{\mathrm{\nabla }}_s\cdot \left({\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\mathrm{\nabla }}_s,p\right)=\mathrm{0},\end{array}$$ where $K_{\mathrm{2}}$ can be reformulated in terms of stream function and velocity potential as $$\begin{array}{cc}{\displaystyle }{K}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\mathrm{\nabla }}_s,\cdot ,\left(,\mathit{\psi },{\mathrm{\nabla }}_s,\mathit{\psi },\right),-,\mathit{\psi },{\mathrm{\nabla }}_s^{\mathrm{2}},\mathit{\psi },+,{\mathrm{\nabla }}_s,\cdot ,\left(,\mathit{\chi },{\mathrm{\nabla }}_s,\mathit{\chi },\right),-,\mathit{\chi },{\mathrm{\nabla }}_s^{\mathrm{2}},\mathit{\chi }\right]\\ & {\displaystyle }& {\displaystyle }+J_s(\mathit{\psi },\mathit{\chi }).\end{array}$$

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 or dry pseudo-density . Here, we denote our density variable by ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s$ and assume no external sources or sinks of $\stackrel{\mathrm{̃}}{\mathit{\rho }}$. Multiplying Eq. () by ${\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s$ and using Eq. () gives $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{s}w}{\partial t}}\right)}_s=& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}}\cdot ({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-{\mathrm{\nabla }}_s\cdot \left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}}w\right)-{\displaystyle \frac{\partial }{\partial s}}\left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s\dot{s}w\right)\\ & {\displaystyle }& {\displaystyle }+{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s\left({\displaystyle \frac{\partial s}{\partial z}}\right)\left[-,{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial s}},-,{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\displaystyle \frac{\partial p}{\partial s}}\right]-(\mathrm{2}\mathbf{\Omega }\times {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}})\cdot \mathit{k}.\end{array}$$ Similarly, from Eq. (), we have $$\begin{array}{cc}{\displaystyle }{\left({\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}}}{\partial t}}\right)}_s=& {\displaystyle }-{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_{s}w({\mathit{u}}_{\mathrm{h}}\cdot \mathrm{\nabla }\mathit{k})-{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathrm{\nabla }}_s{K}_{\mathrm{2}}-{\mathit{\zeta }}_s\mathit{k}\times {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}}\\ {\displaystyle }& {\displaystyle }-{\mathit{u}}_{\mathrm{h}}{\mathrm{\nabla }}_s\cdot \left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s{\mathit{u}}_{\mathrm{h}}\right)-{\displaystyle \frac{\partial }{\partial s}}\left({\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s\dot{s}{\mathit{u}}_{\mathrm{h}}\right)\\ & {\displaystyle }& {\displaystyle }-{\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s\left({\mathrm{\nabla }}_z,\mathrm{\Phi },+,{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}},{\mathrm{\nabla }}_z,p\right)-[\mathrm{2}\mathbf{\Omega }\times {\stackrel{\mathrm{̃}}{\mathit{\rho }}}_s\mathit{u}{]}_z.\end{array}$$

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.

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, Abhishekh 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 Banner, 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