DYNAMICO-1 . 0 , an icosahedral hydrostatic dynamical core designed for consistency and versatility

The design of the icosahedral dynamical core DYNAMICO is presented. DYNAMICO solves the multi-layer rotating shallow-water equations, a compressible variant of the same equivalent to a discretization of the hydrostatic primitive equations in a Lagrangian vertical coordinate, and the primitive equations in a hybrid mass-based vertical coordinate. The common Hamiltonian structure of these sets of equations is exploited to formulate energy-conserving spatial discretizations in a unified way. The horizontal mesh is a quasi-uniform icosahedral C-grid obtained by subdivision of a regular icosahedron. Control volumes for mass, tracers and entropy/potential temperature are the hexagonal cells of the Voronoi mesh to avoid the fast numerical modes of the triangular C-grid. The horizontal discretization is that of Ringler et al. (2010), whose discrete quasi-Hamiltonian structure is identified. The prognostic variables are arranged vertically on a Lorenz grid with all thermodynamical variables collocated with mass. The vertical discretization is obtained from the three-dimensional Hamiltonian formulation. Tracers are transported using a second-order finite-volume scheme with slope limiting for positivity. Explicit Runge–Kutta time integration is used for dynamics, and forward-in-time integration with horizontal/vertical splitting is used for tracers. Most of the model code is common to the three sets of equations solved, making it easier to develop and validate each piece of the model separately. Representative three-dimensional test cases are run and analyzed, showing correctness of the model. The design permits to consider several extensions in the near future, from higher-order transport to more general dynamics, especially deep-atmosphere and non-hydrostatic equations.


Introduction
In the last 2 decades, a number of groups have explored the potential of quasi-uniform grids for overcoming wellknown deficiencies of the latitude-longitude mesh applied to atmospheric general circulation modelling (Williamson, 2007).Particularly compelling has been the computational bottleneck created by the convergence of the meridians at the pole, which prevents efficient distribution of the computational load among many computers.Quasi-uniform grids have no such singular points and are free of this bottleneck.The first attempts at using quasi-uniform grids (Sadourny et al., 1968;Sadourny, 1972) failed at delivering important numerical properties that could be achieved on Cartesian longitude-latitude grids (Arakawa, 1966;Sadourny, 1975a, b;Arakawa and Lamb, 1981).For this reason the balance has been in favour of longitude-latitude grids until the recent advent of massively parallel computing that provided a strong incentive to revisit these grids.
Since one reason for using quasi-uniform grids is the capability of benefiting from the computing power of massively parallel supercomputers, many groups have set highresolution modelling as a primary target.For the dynamical core, which solves the fluid dynamical equations of motion, T. Dubos et al.: DYNAMICO-1.0, an icosahedral hydrostatic dynamical core this generally implies solving a non-hydrostatic set of equations.Indeed, the hydrostatic primitive equations commonly used in climate-oriented global circulation models (GCMs) assume that the modelled motions have horizontal scales much larger than the scale height, typically about 10 km on Earth.Some hydrostatic models on quasi-uniform grids have been developed but essentially as a milestone towards a nonhydrostatic model (Wan et al., 2013).
In fact in many areas of climate research high-resolution modelling can still be hydrostatic.For instance palaeoclimate modelling must sacrifice atmospheric resolution for simulation length, so that horizontal resolutions typical of the Coupled Model Intercomparison Project (CMIP)-style climate modelling are so far beyond reach, and would definitely qualify as high-resolution for multi-millenial-scale simulations.Similarly, three-dimensional modelling of giant planets is so far unexplored since resolving their small Rossby radius requires resolutions of a fraction of a degree.Modelling at Institut Pierre Simon Laplace (IPSL) focusses to a large extent on climate timescales and has diverse interests ranging from palaeo-climate to modern climate and planetology.When IPSL embarked in 2009 in an effort to develop a new dynamical core alongside Laboratoire de Météorologie Dynamique -Zoom (LMD-Z) (Hourdin et al., 2013), a medium-term goal was therefore set to focus on hydrostatic dynamics in order to best serve the IPSL community with increased efficiency and versatility.
By versatility we mean the ability to relax in the dynamical core certain classical assumptions that are accurate for the Earth atmosphere but not necessarily for planetary atmospheres, or may have small but interesting effects on Earth.For instance in LMD-Z it is possible to assume for dry air a non-ideal perfect gas with temperature-dependent thermal capacities and this feature is used to model Venus (Lebonnois et al., 2010).In a similar vein a parallel effort has been undertaken to relax the shallow-atmosphere approximation in LMD-Z and solve the deep-atmosphere quasihydrostatic equations (White and Bromley, 1995;Tort and Dubos, 2014a).Although this feature is not yet implemented in DYNAMICO, the same prognostic variables have been adopted in DYNAMICO as in the deep-atmosphere LMD-Z (Tort et al., 2014b), in order to facilitate upcoming generalizations of DYNAMICO, including generalizations to nonhydrostatic dynamics.
LMD-Z is a finite-difference dynamical core but the kinematic equations (transport of mass, entropy/potential temperature, chemical species) are discretized in flux form, leading to the exact discrete conservation of total mass, total entropy/potential temperature and species content.Upwindbiased reconstructions and slope limiters are used for the transport of species, which is consistent with mass transport and monotonic (Hourdin and Armengaud, 1999).Horizontal dynamics is discretized in vector-invariant form following the enstrophy-conserving scheme of Sadourny (1975b).Unlike the vast majority of hydrostatic dynamical cores, Sim-mons and Burridge (1981) is not used for vertical momentum transport and hydrostatic balance.Another discretization is used, which also exactly preserves energy (Hourdin, 1994).Due to this emphasis on exact discrete conservation properties in LMD-Z, a critical design goal of DYNAMICO was to have at least equivalent properties of conservation and consistency.
Pursuing both objectives of consistency and versatility (as defined above) implies that generic approaches must be found, rather than solutions tailored to a specific equation set.For instance the equivalence of mass and pressure, the proportionality of potential and internal energies are valid only for the hydrostatic primitive equations and cease to be valid in a deep-atmosphere geometry, or even in a shallowatmosphere geometry with a complete Coriolis force (Tort and Dubos, 2014a).The Bernoulli function appearing in the vector-invariant form of the equations of motion is the sum of kinetic energy and geopotential only if an ideal perfect gas (with temperature-independent thermal capacities) is assumed (Tort and Dubos, 2014b).The same assumption is required to have internal energy and enthalpy proportional to temperature, as in Simmons and Burridge (1981).For versatility the dynamical core should not critically rely on such accidental relationships.This raises the question of what assumptions can be made that are both common to all potential target equation sets and sufficient to obtain the desired consistency properties.The answer to this question that has emerged during the DYNAMICO project is that the Hamiltonian formulation of the equations of motion is a sufficient common structure from which discrete consistency can be obtained for all well-formed equation sets.This idea is not really new.In fact it has been advocated for some time now by Salmon (1983Salmon ( , 2004) ) who applied it to the Saint-Venant equations.However, the Hamiltonian approach has been applied only once to date to derive a full-fledged three-dimensional dynamical core, by Gassmann (2013).Gassmann (2013) uses the Hamiltonian formulation of the fully compressible equations in Eulerian coordinates.This Hamiltonian theory has been recently extended for compressible hydrostatic flows and for non-Eulerian vertical coordinates (Tort and Dubos, 2014b;Dubos and Tort, 2014) and serves as the basis to formulate the discretization of dynamics in DYNAMICO.
In addition to the above approach, building blocks for DY-NAMICO include a positive-definite finite-volume transport scheme (Lauritzen et al., 2014b) and finite-difference operators generalizing Sadourny's scheme to general unstructured spherical meshes.A partial generalization has been achieved by Bonaventura and Ringler (2005) but still lacked a discrete conservation of potential vorticity/potential enstrophy and exact discrete geostrophic equilibria, two properties tied together as discussed by Thuburn (2008).A full generalization was obtained later by Thuburn et al. (2009) and Ringler et al. (2010) assuming a Delaunay-Voronoi pair of primal and dual meshes, or more generally orthogonal primal and dual meshes (see Sect. 2).Thuburn et al. (2014) further generalize to a wide class of non-orthogonal dual meshes, targeting the cubed sphere which has a better balance between the degrees of freedom for mass and velocity, thus avoiding numerical modes present in triangular meshes and their dual (Gassmann, 2011;Weller et al., 2012).However the accuracy of finite differences on the cubed sphere is poor and a triangular-hexagonal grid yields much more accurate results for a number of degrees of freedom similar to that of the cubed sphere (Thuburn et al., 2014).On Delaunay-Voronoi meshes placing mass inside triangles leads to a branch of non-stationary numerical modes that must be controlled by a non-trivial amount of dissipation (Rípodas et al., 2009;Wan et al., 2013), while placing mass inside Voronoi domains leads to stationary numerical modes, which requires no or very little dissipation for stable integrations (Ringler et al., 2010;Skamarock et al., 2012;Gassmann, 2013).DYNAM-ICO follows the second option.
The present paper is organized as follows.Section 2 describes how the transport of mass, potential temperature and other tracers is handled by DYNAMICO.For this the grid and the discrete representation of scalar and vector fields are introduced.Mass fluxes through control volumes boundaries are provided by the dynamics, as described in Sect.3. Following the Hamiltonian approach, the primary quantity is the total energy, which is discretized first vertically then horizontally then yields the discrete expressions for the Bernoulli function and other quantities appearing in the curl-form equation of motion.Section 4 is devoted to energetic consistency.The discrete energy budget of DYNAMICO is derived, and the underlying Hamiltonian structure of the TRiSK scheme (Thuburn et al., 2009;Ringler et al., 2010) is identified.In Sect. 5 sample numerical results are presented, verifying the correctness of DYNAMICO and its ability to perform climate-style integrations.Our main contributions are summarized and discussed in Sect.6, and future work is outlined.

Kinematics
In this section we describe how the transport of mass, potential temperature and other tracers is handled by DYNAM-ICO, using mass fluxes computed by the dynamics as described in Sect.3. We use bold face letters for vectors in three-dimensional physical space and for points on the unit sphere.Space-dependent fields are functions of a vector n on the unit sphere and a generalized vertical coordinate η.Especially the geopotential (n, η, t) is a dependent quantity.Using the dot notation for the Lagrangian (material) derivative, u = ṅ is an angular velocity tangent to the unit sphere , i.e. n • u = 0.The Eulerian position r of a fluid parcel in physical space is determined by the geopotential considered as a vertical Eulerian coordinate and n, i.e. r = r( , n).An expression for r( , n) is not needed to solve the transport equa-tions and needs to be specified only when dealing with the dynamics (see Sect. 3).Denoting ∂ α = ∂/∂α for α = n, η, t, the continuous flux-form budget for mass, potential temperature θ and tracer q are (1) where µ is the pseudo-density such that total mass is µd 2 ndη, = µθ, Q = µq, U = µu is the horizontal mass flux vector, W = µ η is the mass flux per unit surface through model layers η = cst.
The following subsections describe the grid, indexing conventions, the discrete mass and potential temperature budgets, and finally the positive-definite finite-volume scheme used for additional tracers.

Icosahedral-hexagonal grid, staggering and discrete objects
The mesh is based on a tessellation of the unit sphere (Sadourny et al., 1968).Each triangle has a global index v and each vertex has a global index i.Several points are associated with each index i or v. Mesh generation and smoothing is described in Appendix A and numerical stability issues arising in the calculation of spherical-geometric entities are raised and solved in Appendix B. By joining v points one obtains the hexagonal-pentagonal mesh, with control volumes indexed by i and vertices indexed by v. Mass will be associated with hexagonal control volumes and i points, so we will refer to this mesh as the primal mesh, while the triangular mesh will be referred to as dual.Additional quantities are associated with primal edges joining v points, and dual edges joining i points.Both types of edges are indexed by e.These notations follow Thuburn et al. (2009).Lorenz staggering is used in the vertical.Full vertical levels are indexed by k = 1. ..K. Interfaces between full levels are indexed by l = 1/2. ..L = K + 1/2.Following the spirit of discrete exterior calculus (DEC; see, e.g.Thuburn and Cotter, 2012), we associate to each scalar or vector field a discrete description reflecting the underlying differential-geometric object, i.e. 0-forms (scalar functions), 1-forms (vector fields with a curl), 2-forms (vector fields with a divergence) and 3-forms (scalar densities).Scalar densities include µ and .We describe them by discrete values µ ik , ik defined as their integral over the threedimensional control volumes (µ ik is in units of kg).Scalar functions include θ ik = ik /µ ik and specific volume α ik ; 2forms include the fluxes of mass and potential temperature.The horizontal mass flux vector U = µu is described by its integrals U ek over a vertical boundary between two hexagonal control volumes and the vertical mass flux per unit surface W = µ η by its integral W il over the pseudo-horizontal boundary between two adjacent control volumes located one above another (unit: kg s −1 ).Averages and finite differences are decorated with the location of the result, i.e. δ k lies at full levels, m l lies at interfaces, and ek is collocated with U ek (see Fig. 1).Especially, using the notations of Ringler et al. ( 2010) t ev v ek .
Operators δ i , δ e and δ v are discrete versions of the twodimensional div, grad and curl operators.They are mimetic in the sense that they satisfy for any A e , B i the discrete formulae e A e δ e B + i B i δ i A = 0, (4) Equation ( 4) is a discrete integration-by-parts formula and Eq. ( 5) imitates curl grad = 0 (Bonaventura and Ringler, 2005).Notice that the generic A, B used here are unrelated to quantities A (areas) and B (Bernoulli function) defined later.

Discrete mass, potential temperature and tracer budgets
The discrete mass and potential temperature budgets are written in flux form: where we omit certain indices when there is no ambiguity (e.g. in Eq. 7 we omit the index e of θ * ek and U ek since operator δ i is always applied to quantities located on edges) and θ * ek , θ * il are values of θ reconstructed at interfaces between control volumes.Currently simple centred averages are used: Either a Lagrangian vertical coordinate or a mass-based vertical coordinate can be used.In the former case W = 0. Notice that if W = 0 and θ ik = θ k is initially uniform, it will remain so at later times for adiabatic motion.This corresponds to using an isentropic/isopycnal vertical coordinate.In the latter case (mass-based vertical coordinate) only the column-integrated mass M i is prognostic, while µ ik is diagnosed from M i : with a l , b l predefined profiles satisfying a = 0, b = 0 at the top and a = 1, b = 0 at the bottom.Then summing Eq. ( 6) over k and using no-flux top and bottom boundary conditions for W provides a prognostic equation for M i : 6) complemented by boundary conditions W = 0 at top and bottom is a diagnostic equation for W il .Equations ( 6)-( 7) are marched in time together with the dynamics using a Runge-Kutta time scheme with a time step τ (see Sect. 3).On the other hand, the additional tracers q are weakly coupled to the dynamics and can be stepped forward with a larger time step t = N transport τ with 1/N transport larger than the maximum Mach number in the flow.To this end, using simple bookkeeping, the dynamics provides timeintegrated fluxes U ek ,W il (both in units of kg) such that where δ t is a finite difference over N transport full Runge-Kutta time steps.Then Eq. ( 3) is discretized using horizontal-vertical splitting (Easter, 1993;Hourdin and Armengaud, 1999): where ik for m = 1, 2 are intermediate values, and q (m) are pointwise values of the tracer reconstructed from Q (m) and µ (m) (see below).The reconstruction operators satisfy the consistency principle that q (m) = 1 whenever 0) ; i.e. the tracer budget is consistent with the mass budget.
The vertical reconstruction is one-dimensional, piecewise linear, slope limited, and identical to Van Leer's scheme I (Van Leer, 1977;Hourdin and Armengaud, 1999).The horizontal advection scheme is identical to SLFV-SL of Lauritzen et al. ( 2012) and is detailed in Dubey et al. (2015).It relies on cellwise-linear reconstructions of q.For this a gradient is estimated in each cell using nearby values (Satoh et al., 2008) and limited to maintain positivity (Dukowicz and Kodis, 1987).The position at which the reconstructed value is evaluated is determined in a semi-Lagrangian fashion (Miura, 2007).

Dynamics
We now turn to the discretization of the momentum budget.A Hamiltonian formulation of the hydrostatic primitive equations in a generalized vertical coordinate is used (Dubos and Tort, 2014).From this formulation the energy budget is obtained invoking only integration by parts, a structure easy to reproduce at the discrete level in order to conserve energy.Before arriving, at the end of this section, at the fully discrete three-dimensional equations, we start from the Hamiltonian of the hydrostatic primitive equations.Introducing a vertical discretization (of the Hamiltonian) produces (the Hamiltonian of) a compressible multi-layer Saint-Venant model.The Boussinesq approximation, enforced by a Lagrange multiplier, yields a standard multi-layer Saint-Venant model.Finally, the horizontal discretization is described.

Continuous Hamiltonian
An ideal perfect gas with pα = RT and constant C p = R/κ is assumed where p is pressure, α specific volume and T temperature.Then where π is the Exner function and θ potential temperature.Note that, letting U (α, θ ) be specific internal energy, ∂U/∂α = −p, ∂U/∂θ = π , U + αp − θ π = 0. We work within the shallow-atmosphere and spherical geopotential approximation, so that gravity g is a constant, the elementary volume is The primitive equations are generated by the Hamiltonian: where f (n, η) = f d 2 n with the unit sphere and v = a 2 (u + n × ) is prognostic (Dubos and Tort, 2014).In Eq. ( 12) H is a functional of the three-dimensional fields µ, v, , and The terms in the integral are kinetic, internal and potential energy.The last term in Eq. ( 12) represents the work of pressure p ∞ exerted at the top η = 1 of the computational domain and sets the upper boundary condition p = p ∞ .Discretizing Hamiltonian Eq. ( 12) in the vertical direction yields a multi-layer Hamiltonian (Bokhove, 2002): where In order to reduce Eq. ( 13) to a multi-layer shallow-water system, the Boussinesq approximation is made by introducing into Eq.( 13) Lagrange multipliers λ k enforcing µ k = a 2 ρ r δ k g : where θ k is now the non-dimensional buoyancy of each layer.
Notice that the last term can be omitted (p ∞ = 0).Indeed,

Fully discrete Hamiltonian
We now discretize horizontally the Hamiltonians (Eqs.14, 13, 12).In addition to the kinematic degrees of freedom µ ik , ik we need to discretize the velocity degrees of freedom.Since we shall need the curl of v, it is a 1-form in the nomenclature of discrete differention geometry.Hence we describe v by the discrete integrals v ek = e v(n, η k ) • dl (unit: m 2 s −1 ), where e is a triangular edge.An approximation of H is then given by where R e = a 2 e ( × n) • dl is the planetary contribution to v e , d e is the (angular) length of triangular edge e and A ie is an (angular) area associated with edge e and to a cell i to which its belongs, with A ie = 0 if e is not part of the boundary of cell i. u ek is a first-order estimate of the component of u along e .In planar geometry, A ie = 1 4 l e d e is a consistent formula for A ie because it satisfies A i = e A ie (Ringler et al., 2010).It is therefore also consistent in spherical geometry, with A i e A ie .Letting A ie = 1 4 l e d e simplifies somewhat the kinetic energy term: Comparing Eqs. ( 15) and ( 13) it is clear that Eq. ( 15) is also a valid horizontal discretization of Eq. ( 13).Regarding Eq. ( 14), a discretization of the kinetic energy part is simply K as above.The other terms are discretized in a straightforward way: with λ ik the pointwise value of λ (0-form).

Discrete equations of motion
We now write the equations of motion corresponding to the discrete Hamiltonians.First, mass fluxes must be computed for use by kinematics.They are computed as U ek is therefore a centred estimate of the mass flux across the face orthogonal to edge e .Next hydrostatic balance is expressed as ∂H /∂ il = 0 or equivalently H = 0, where H is induced by arbitrary, independent variations of only.For the compressible Hamiltonian Eq. ( 15) this yields Therefore, a 2 A i δ l p i + gµ i l = 0 with the upper boundary condition p iK = p ∞ +gµ iL /(2a 2 A i ).These are discrete versions for ∂ η p + µg = 0 and p(η = 1) = p ∞ .p ik can be determined starting from the top level.Alternatively one can define a pressure p * il at layer interfaces by When η is mass-based, one finds from Eq. ( 9) that p * il = a l p s i + C l with surface pressure p s i = p ∞ +gM i and C l = gb l +(1−a l )p ∞ ; i.e. the usual way to diagnose the vertical pressure profile from surface pressure is recovered.Once p ik has been determined, the specific volume α ik = α(p ik , θ ik ) follows.The geopotential is obtained by integrating: starting from the ground, where s i is the time-independent surface geopotential.
On the other hand, for the incompressible Hamiltonian Eq. ( 16), geopotential il is obtained by enforcing the constraint ∂H /∂λ ik = 0, i.e.Eq. ( 18) but with specific volume α ik = 1/ρ r independent from pressure.Furthermore, Therefore, λ ik satisfies the same equations as p ik but with (1 − θ ik ) µ ik instead of µ ik , which shows that θ ik acts indeed as a buoyancy θ = (ρ r − ρ) /ρ r .The Lagrange multipliers λ ik enforcing the incompressibility constraint are to be interpreted as the pressure at full model levels, a typical outcome within the Boussinesq approximation (Holm et al., 2002).Finally, the horizontal momentum balance is written in vector-invariant form.When W = 0, where and the ⊥ operator is defined in Ringler et al. (2010) through antisymmetric weights w ee = −w e e : (q k U k ) ⊥ e = e w ee q ee U e where q ee = q * e k + q * ek 2 (20) with q * ek a value of potential vorticity reconstructed at e points from values at v points q vk = δ v v k /µ v , where µ vk is µ integrated over triangular control volumes defined as an area-weighted sum of neighbouring µ ik (Ringler et al., 2010).Currently, a centred average q * ek = q k e is used but other reconstructions, including upwind-biased reconstructions, could be used as well (Ringler et al., 2010).The weights w ee are obtained by Thuburn et al. (2009), Eq. ( 33) as a function of the ratios R iv = A iv /A i satisfying v R iv = 1, i.e.
v A iv = A i .Using the compressible Hamiltonian Eq. ( 15) one finds where is an approximation of kinetic energy 1 2 a 2 u • u.Therefore geopotential at full levels is defined as a centred average of il and Exner pressure is diagnosed in each control volume using the equation of state.Because p ik = p(α ik , θ ik ), Eq. ( 21) simplifies to π ik = c p (p ik /p r ) κ .In practice π ik and α ik are both diagnosed from p ik , θ ik when solving the hydrostatic balance.
On the other hand, using the incompressible Hamiltonian Eq. ( 16) yields As already mentioned, changing p ∞ only modifies the upper boundary condition and only adds a constant to λ ik .Since only δ e B k is important for dynamics, the value of p ∞ is arbitrary and can be set to 0. Now if θ ik = θ k is horizontally uniform, θ * ek = θ k and and Eq. ( 19) takes the expected form for a multi-layer shallow-water model.In the more general case where θ ik is not uniform, Eq. ( 19) is a discretization of the vectorinvariant form of Ripa equations (Ripa, 1993).
When W = 0 an additional term takes into account vertical momentum transport: where v * el is a value of v e reconstructed at interfaces.Here a centred average v * el = v e l is used.The above discretization does not possess particular conservation properties and other equally accurate formulae could be explored.

Time marching
After spatial discretization one obtains a large set of ordinary algebraic equations where y = (M i , ik , v ik ) with a mass-based coordinate and y = (µ ik , ik , v ik ) with a Lagrangian coordinate.Geopotential ik is diagnosed from y when computing the trends f (y) (details below).Equation ( 26) is advanced in time using a scheme of Runge-Kutta type.Temporal stability is limited by the external mode, which propagates at the speed of sound c.For a p stage scheme, about (p/C max ) × cT /δx evaluations of f are necessary to simulate a time T with resolution δx where C max is the maximum time step allowed to integrate the differential equation dx/dt = ix.It is therefore desirable to maximize the effective Courant number C eff = C max /p.The design goals of the time scheme are to be fully explicit for simplicity, second-order accurate and with a favourable effective Courant number for efficiency.Two-stage Runge-Kutta schemes of the order of 2 are unconditionally unstable for imaginary eigenvalues and ruled out.All explicit p step RK schemes of order p are equivalent for linear equations; p = 3 and p = 4 yield C (3) Kinnmark and Gray (1984b) provided p stage Runge-Kutta schemes with optimal C max = p − 1 and an order of 2 for odd p (referred to as RK2.p below).Third-and fourthorder accuracy are achievable at a small price in terms of stability, i.e.C max = (p − 1) 2  − 1 (Kinnmark and Gray, 1984a).Hence, for p = 4 and p = 5 optimal schemes are RK4 and RK2.5, the latter having C eff = 0.8, about 13 % larger than C RK4 eff .Currently, the following scheme y n −→ T. Dubos et al.: DYNAMICO-1.0, an icosahedral hydrostatic dynamical core y n+1 is implemented for RK4: where τ ≤ 2 √ 2δx/c is the time step and y n y(nτ ).This is a low-storage scheme since the same memory space can be used for y 1 , y 2 , y 3 and y n+1 .It is also very easy to implement.It is fourth-order accurate for linear equations but only second-order accurate for non-linear equations.A similar sequence is used for RK2.5.
Furthermore, the last step is similar to an Euler step; hence, so that the time-integrated mass fluxes expected by the transport scheme are simply U ek = τ U 3 ek , W il = τ W 3 il or their sum over N transport successive time steps (see Sect. 3).

Recap: computation of trends in a mass coordinate
At the beginning of this computation v ek , M i , ik are known.Cell-integrated mass µ ik and potential temperature θ ik are diagnosed using Eqs.( 9) and (8).Pressure p ik follows from hydrostatic balance (see Sect. 3.3), then Exner pressure and specific volume π ik , α ik .Geopotential is obtained bottomup using Eq. ( 18), then the Bernoulli function (Eqs. 22 and 23).
From µ ik , v ek horizontal mass fluxes U ek are obtained then, by vertical integration, ∂M i /∂t.Then ∂µ i /∂t is obtained and injected into the mass budget (Eq. 1) to compute the vertical mass flux W il by a top-down integration.The potential temperature fluxes and trend are then computed using Eqs.( 7) and (8).Finally, the velocity trend is computed following Eq.( 25).
From µ ik , v ek horizontal mass fluxes U ek are obtained then ∂µ i /∂t.The trends of potential temperature and velocity are finally computed using Eq. ( 7) with W il = 0 and Eq. ( 19).

Filters
Centred schemes need stabilization to counteract the generation of grid-scale features in the flow.Linear sources of gridscale noise, e.g.dispersive numerical errors, may be handled by filters, e.g.upwinding or hyperviscosity.Other sources are genuinely non-linear, e.g. the downward cascade of energy or enstrophy.Here we handle these sources through hyperviscosity as well, rather than with a proper turbulence model, e.g.Smagorinsky (1963), following a widespread although disputable practice (see Gassmann, 2013).
For this purpose hyper-diffusion is applied every N diff time step in a forward-Euler manner: v ek := (28) where the exponent p is 1 or 2, the dissipation timescales τ θ , τ ω , τ δ serve to adjust the strength of filtering, the length scales L θ , L ω , L δ are such that L −2 θ , L −2 ω , L −2 δ are the largest eigenvalue of the horizontal dissipation operators D θ , D ω , D δ defined as These positive-definite operators correspond to diffusing a scalar, vorticity and divergence.Notice, however, that filtering with p > 1, although it damps grid-scale noise, typically neither removes oscillations entirely nor guarantees positivity of the filtered field (see e.g.Jiménez, 2006).
are precomputed by applying D θ , D ω , D δ many times in sequence on random data, so that their largest eigenvalue is given by ratio of the norm of two successive iterates.This process converges very quickly and in practice 20 iterations are sufficient.The dissipation timescales and the exponents can be set to different values for θ, ω, δ.N diff is determined as the largest integer that ensures stability, i.e. such that N diff τ be smaller than all three dissipation timescales.

Conservation and stability
In addition to its aesthetic appeal, discrete conservation of energy has practical consequences in terms of numerical stability, which we discuss here using arguments similar to Geosci.Model Dev., 8, 3131-3150, 2015 www.geosci-model-dev.net/8/3131/2015/energy-Casimir stability theory (Arnold, 1965).Indeed, if a dynamical system conserves a convex integral quantity, then any state of the system which is a minimum of that quantity is necessarily a stable steady state.For instance the states of rest of the shallow-water equations minimize a linear combination of total energy and mass.Each additional conserved integral quantity widens the family of steady states that can be proven to be stable.In the discussion below we assume that the discrete equations of motion conserve total energy.The additional conserved quantities then depend on the vertical coordinate used.Assuming a Lagrangian vertical coordinate, the additional integral quantities conserved by the discrete equations of motion are, for each layer, the horizontally integrated mass and potential temperature i µ ik , i ik , which form a subset of the Casimir invariants of the continuous equations (Dubos and Tort, 2014) The above reasoning shows that linearization of the discrete equations of motion around a steady state making H convex yields linear evolution equations with purely imaginary eigenvalues.Forward integration in time is then linearly stable provided the relevant Courant-Friedrichs-Lewy condition is satisfied.In particular, it is not necessary for linear stability that the time-marching scheme conserves energy.
With a mass-based vertical coordinate, the exchange of mass between layers reduces the set of discrete Casimir invariants to total mass and potential temperature i M i , ik ik .Considering the linear combination H = H − i M i −π ik ik one finds the condition ∂H /∂ ik = π .It is impossible to satisfy both hydrostatic balance and a uniform Exner pressure; hence, no feasible state minimizes H . On the other hand, if cell-integrated entropy S ik is prognosed instead of potential temperature, one finds that isothermal states of rest minimize H = H − i M i − T ik S ik (not shown).
We now proceed to derive the discrete energy budgets corresponding to a Lagrangian and a mass-based vertical coordinate.In these calculations only the adiabatic terms are considered, and the effect of the hyperviscous filters is omitted.

Lagrangian vertical coordinate
When W = 0, the continuous-time energy budget reads: Using the discrete integration-by-parts formula (4) and the antisymmetry property w ee +w e e = 0, one finds dH /dt = 0.More generally, similar calculations yield the temporal evolution of an arbitrary quantity F (µ ik , ik , v ek , il , λ ik ): Equations ( 29)-( 32) imitate at the discrete level the Hamiltonian formulations obtained in Dubos and Tort (2014).Discrete conservation of energy then appears as a consequence of the antisymmetry of the brackets {F, H } µ , {F, H } , {F, H } v , the formulation of hydrostatic balance as ∂H /∂ il = 0, and, in the incompressible case, of the constraint ∂H /∂λ ik = 0.The antisymmetry of {F, H } µ , {F, H } is equivalent to the discrete integrationby-parts formula (4), itself equivalent to the discretization of the horizontal div and grad operators being compatible (see e.g.Taylor and Fournier, 2010).The antisymmetry of {F, H } v results from w ee = −w ee and q ee = q e e (Ringler et al., 2010).

One-layer shallow-water equations
In the simplest case of a single layer without topography ( s = 0), the incompressible Hamiltonian Eq. ( 14) with = where i = gh i is the geopotential at the "top" of the model and we have taken into account the constraint µ i = A i h i , where h i is interpreted as the thickness of the fluid layer.
Hamiltonian H is precisely the one considered in Ringler et al. (2010).The discrete equations of motion also reduce to their energy-conserving scheme (not shown).Equation ( 29) reduces to This is a discrete imitation of the shallow-water Poisson bracket.Had we used the enstrophy-conserving scheme of Ringler et al. (2010) instead of the energy-conserving scheme, {F, H } v would have been: This discrete bracket is not antisymmetric.Comparing Eqs. ( 32) and ( 34), one sees that the energy-conserving bracket ( 32) is the antisymmetrization of Eq. ( 34); i.e.

{F, H }
In the limit of the linearized shallow-water equations on the f sphere (Thuburn et al., 2009), both brackets (32)-( 34) reduce to where f is the constant value of the Coriolis parameter and h is the background fluid layer thickness; i.e. h e = h + h e , h e h.
In Ringler et al. (2010), the energy-conserving discretizations of the mass flux, kinetic energy and Coriolis term were devised by choosing a certain form and stencil for each of them with undetermined coefficients, deriving the energy budget, and choosing the undetermined coefficients in such a way that all contributions cancel out.In hindsight this delicate task could have been avoided by following the approach used here, inspired by Gassmann (2013) and advocated since some time already by Salmon (1983Salmon ( , 2004)): discretizing the energy and the brackets, instead of the equations of motion themselves.The critical part is to discretize the brackets.Starting from the linearized bracket (35) implicitly derived in Thuburn et al. (2009), a straightforward non-linear generalization is Eq. ( 34), which can be antisymmetrized to yield Eq. ( 32).From this point of view all the critical building blocks of Ringler et al. (2010) were already obtained in Thuburn et al. (2009).The present approach generalizes this scheme to three-dimensional equations in a generalized vertical coordinate, exploiting recent advances in the relevant Hamiltonian formulation (Dubos and Tort, 2014).

Mass-based vertical coordinate
When a mass-based coordinate is used instead of a Lagrangian vertical coordinate, additional terms proportional to the vertical mass flux W il appear in the equations of motion and in the energy budget.These terms cancel each other for the continuous equations but not necessarily for the discrete equations.It is possible to obtain a cancellation by imitating at the discrete level a relationship between the functional derivatives of H due to invariance under a vertical relabelling (remapping) (Dubos and Tort, 2014).This strategy has been recently implemented in a longitude-latitude deepatmosphere quasi-hydrostatic dynamical core (Tort et al., 2014b).Tort et al. (2014b) estimate the numerical heat source due to the vertical transport terms as less than 10 −3 W m −2 in idealized climate experiments (Held and Suarez, 1994).Hence, cancelling this very small numerical heat source is not yet implemented in DYNAMICO and energy is not exactly conserved when a mass-based vertical coordinate is used.
So far we see no indication that this would damage longduration simulations (see numerical results in Sect.5) but in the future strict energy conservation may be offered as an option, together with the choice to prognose entropy instead of potential temperature.

Relation with Gassmann (2013)
At this point some important differences with respect to the approach of Gassmann (2013) can be highlighted.Firstly, since the vertical coordinate is non-Eulerian, geopotential depends on time and appears as an argument of the Hamiltonian.It therefore produces additional terms in the energy budget that vanish as shown in Sect.4.2.On the other hand vertical momentum is not prognostic, since the equations are hydrostatic.
Second, Gassmann (2013) prognoses contravariant momentum components while we prognose v ek , which are equivalent to covariant velocity components.Indeed, the latter appear as the preferred prognostic variables in the Euler-Lagrange equations of motion (Tort and Dubos, 2014b) and their Hamiltonian formulation in a general vertical coordinate (Dubos and Tort, 2014).An immediate advantage of prognosing v ek is that vorticity is trivially and naturally obtained along the lines of DEC.Furthermore, the horizontal mass flux appears in the mass and tracer budgets through its contravariant components, and the functional derivatives of the Hamiltonian with respect to covariant momentum  components are directly the contravariant mass flux components.This translates directly at the discrete level into Eq.( 17).As a result the discrete Poisson bracket is almost trivial, with the exception of the Coriolis part.Hence, averaging only occurs where it is unavoidable due to mesh staggering, i.e. in the discrete Coriolis term and in the Hamiltonian.
In a Eulerian formulation of the non-hydrostatic Euler equations, prognosing covariant components would have the drawback that the no-flux lower boundary condition involves a linear combination of all three covariant components, which on a staggered mesh may be difficult to discretize in an energy-conserving way.This may be a reason why Gassmann (2013) prognoses rather contravariant components.However, with a mass-based or Lagrangian coordinate and hydrostatic equations, this is not an issue since the no-flux lower boundary condition only enters the mass budget.

Results
In this section, the correctness of DYNAMICO is checked using a few idealized test cases.Horizontal resolutions of M = 32, 40, 64, 80, 128, 160 are used.Defining mean horizontal resolution δx as the distance between hexagon centres on a regular planar hexagonal mesh covering the surface of the Earth with 10 M 2 cells, or equivalently the length of the edges of triangles in a regular triangular mesh with 20 M 2 triangles and the same surface, these values translate into δx 280, 220, 140, 110, 70, 55 km.
Since our horizontal advection scheme is very similar to one scheme studied by Mittal and Skamarock (2010), we do not show two-dimensional results and focus on a threedimensional test case of the Dynamical Core Model Intercomparison Project (DCMIP) suite (Kent et al., 2014).The correctness of the three-dimensional dynamics solver is checked using the dry baroclinic instability set-up of Jablonowski and Williamson (2006).Finally, the forceddissipated set-up defined by Held and Suarez (1994) is carried out to demonstrate the suitability of DYNAMICO for climate type simulations.

Transport by a prescribed Hadley-like meridional circulation
This test case consist of a single layer of tracer, which deforms over the duration of simulation.The flow field is prescribed so that the deformed filament returns to its initial position in the end of simulation.We used resolutions M×N z = 40×30, 80×60, 160×120.The hybrid coefficients are computed so that the model levels are initially uniformly spaced.
Figure 2 shows the tracer profile at t = 12 h and t = 24 h for horizontal resolutions M = 80 and M = 160.At t = 24 h the tracer field should ideally be independent of latitude, so any latitudinal dependance results from numerical errors.Figure 2a shows that for coarse resolution the scheme is diffusive and the final profile is quite diffused particularly at the  2b show that the increasing resolution decreases the diffusive nature of the advection scheme.Moreover, the slope limiter successfully avoids the generation of spurious oscillations in the numerical solution.
Table 1 shows the global error norms for different horizontal and vertical resolutions.
As expected from two-dimensional test cases (Lauritzen et al., 2014b), our transport scheme is more diffusive than finite-volume schemes on essentially Cartesian meshes such as those presented in Kent et al. (2014).Sample solutions (their Fig. 6) and error norms (their Table 6) they present indicate that our scheme achieves, at resolution δx, an accuracy similar to these schemes at resolution 2δx.

Baroclinic instability
The baroclinic instability benchmark of Jablonowski and Williamson (2006) is extensively used to test the response of three-dimensional atmospheric models to a controlled, evolving instability.The initial state for this test case is the sum of a steady-state, baroclinically unstable, zonally symmetric solution of the hydrostatic primitive equation and of a localized zonal wind perturbation triggering the instability in a deterministic and reproducible manner.
Even without the overlaid zonal wind perturbation, the initial state would not be perfectly zonally symmetric because the icosahedral grid, as other quasi-uniform grids, is not zonally symmetric.Therefore, the initial state possesses, in addition to the explicit perturbation, numerical deviations from zonal symmetry.This initial error, as well as truncation errors made at each time step by the numerical scheme, is not homogeneous but reflects the non-homogeneity of the grid.It nevertheless has the same symmetry as the grid, here wave number 5 symmetry.Due to the dynamical instability of the initial flow, the initial error is expected to trigger a wave number 5 mode of instability (provided such an unstable mode with that zonal wave number exists).Depending on the amplitude of the initial truncation error, this mode can become visible, a case of grid imprinting (Lauritzen et al., 2010).
Figure 3 presents results obtained at resolutions M = 32, 64, 128 (mean resolution 280, 140, 70 km) using 30 hybrid vertical levels and fourth-order filters (p = 2 in Eqs.27, 28).Dissipation time and time step are set to τ = 6, 3, 1.5 h and δt = 600, 300, 150 s, respectively.The right column shows the temperature field at pressure level 850 hPa at day 9.At this day the baroclinic wave is well developed.The wave crest is reasonably sharp at M = 32, and becomes sharper at a higher resolution.The simulated temperature field is qualitatively similar to those obtained at comparable resolutions by other models (e.g.Jablonowski and Williamson, 2006, Figs. 6, 7).
The left column shows surface pressure at day 12, after the baroclinic wave has broken, letting time for grid imprinting to develop.Grid imprinting in the Southern Hemisphere, measured quantitatively as in Lauritzen et al. (2010) as the root mean square departure of surface pressure from its unperturbed value of 1000 hPa, exceeds 0.5 hPa at day 9 at M = 32, at day 11 at M = 64 and at day 13 at M = 128.Comparing with Fig. 12 of Lauritzen et al. (2010), these values are in the low end of icosahedral models.Held and Suarez (1994) propose a benchmark to evaluate the statistically steady states produced by the dynamical cores used in climate models.Detailed radiative, turbulence and moist convective parametrization are replaced with very simple forcing and dissipation.The simple forcing and dissipation are designed in terms of a simple relaxation of the temperature field to a zonally symmetric state and Rayleigh damping of low-level winds to represent the Figure 4 presents statistics obtained when using horizontal resolution of 280 km (M = 32) and dissipation time τ = 6 h.The model is stable for longer dissipation times (τ = 24 h) but smaller values produce smoother fields.Statistics obtained at resolution 140 km (M = 64) with τ = 3 h are presented in Fig. 5. First-order statistics (panels ab) are close to those presented in Held and Suarez (1994) and present very little sensitivity to resolution.Second-order statistics are slightly more sensitive to resolution and increase slightly from M = 32 to M = 64.Temperature variance at M = 64 is close to that presented in Held and Suarez (1994) and slightly smaller than that obtained by Wan et al. (2013) on a triangular icosahedral grid at comparable resolution R2B5.

Contributions
A number of building blocks of DYNAMICO are either directly found in the literature or are adaptations of standard methods: explicit Runge-Kutta time stepping, mimetic horizontal finite-difference operators (Bonaventura and Ringler, 2005;Thuburn et al., 2009;Ringler et al., 2010), piecewiselinear slope-limited finite-volume reconstruction (Dukowicz and Kodis, 1987;Tomita et al., 2001), swept-area calculation of scalar fluxes (Miura, 2007), and directionally split time integration of three-dimensional transport (e.g.Hourdin and Armengaud, 1999).It is therefore useful to highlight the two specific contributions brought forward, in our opinion, in the design of DYNAMICO, and that can be of broader applicability for model design.
The first contribution is to separate kinematics from dynamics as strictly as possible.This separation means that the transport equations for mass, scalars and entropy use no information regarding the specific momentum equation being solved.This includes the equation of state as well as any metric information, which is factored into the prognosed T. Dubos et al.: DYNAMICO-1.0, an icosahedral hydrostatic dynamical core degrees of freedom and into the quantities derived from them (especially the mass flux).Metric information is not used to prognose tracer, mass and potential temperature.It is confined in a few operations computing the mass flux, Bernoulli function and Exner function from the prognostic variables.This formulation is in line with more general lines of thought known as physics-preserving discretizations (Koren et al., 2014) and discrete differential geometry (Thuburn and Cotter, 2012).Similarly, while we use the exact same hybrid vertical coordinate as most hydrostatic primitive equation models, we insist that it should be considered as mass based rather than pressure based.Indeed, the coincidence (up to time-independent multiplicative and additive factors) of mass and pressure is a peculiarity of the traditional shallowatmosphere hydrostatic equations with a pressure top boundary condition.Recognizing the fundamentally kinematic definition of the hybrid coordinate in terms of mass rather than pressure emphasizes its relevance for solving other equation sets, especially non-hydrostatic (Laprise, 1992).
The second contribution is to combine this kinematicsdynamics separation with a Hamiltonian formulation of the equations of motion to achieve energetic consistency.This approach extends the work of Gassmann (2013) to hydrostatic equations of motion and non-Eulerian vertical coordinates.This extension relies itself on a recent corresponding extension of the Hamiltonian theory of atmospheric fluid motion (Tort and Dubos, 2014b;Dubos and Tort, 2014).The Hamiltonian approach further confines the equationdependent parts of the numerical scheme to a single quantity, the total energy of the system expressed in terms of the prognostic variables and, in the case of hydrostatic equations, geopotential.The latter is a pseudo-prognostic variable that is an argument of the Hamiltonian but is diagnosed at each time step by enforcing the hydrostatic constraint, found to be simply the condition that the derivatives of the Hamiltonian with respect to geopotential degrees of freedom vanish.This variational formulation of hydrostatic balance was first identified in the context of the deep-atmosphere quasi-hydrostatic equations (Tort et al., 2014b) then generalized (Dubos and Tort, 2014) and applied to DYNAMICO within the shallowatmosphere approximation.Ultimately, the choice of a specific equation set boils down to choosing and discretizing the Hamiltonian, without changing the general structure of the algorithm computing the tendencies.
These two advances yield our design goals: consistency and versatility.The desired ability to solve different equation sets is currently limited to the hydrostatic primitive equations and the multi-layer Saint-Venant or Ripa equations, but little work is required to solve other similar equations like the recently derived non-traditional spherical shallow-water equations (Tort et al., 2014a).Whichever set of equations needs to be solved in the future, including the fully compressible Euler equations, energetic consistency is guaranteed if the general approach followed here and in Tort et al. (2014b) is applied.Furthermore, this approach is not limited to finite-difference schemes but can be extended to finite element schemes.
We would also like to emphasize what the Hamiltonian approach does not achieve.Good numerical dispersion crucially depends on grid staggering (for finite differences) or on the finite element spaces used to represent the various quantities.It is entirely possible to design an energy-conserving schemes with disastrous numerical dispersion properties.Other properties, such as exact geostrophic equilibria or a discrete potential vorticity budget, come in addition to the antisymmetry of the discrete Poisson bracket, as discussed in Sect. 4 (see also Cotter and Thuburn, 2014).However, the Hamiltonian formulation provides a divide-and-conquer strategy by allowing for the easy transfer of these additional properties to new sets of equations once they have been obtained for a specific one.

Outlook for DYNAMICO
A Lagrangian vertical coordinate is currently available as an option.In the absence of the vertical remapping that must necessarily take place occasionally in order to prevent Lagrangian surfaces to fold or cross each other, this option can not be used over meaningful time intervals.However, it is convenient for development purposes since it allows one to investigate separately issues related to the vertical and horizontal discretizations.Nevertheless, a future implementation of vertical remapping would be a useful addition.There is room for improvement on other points.In particular, it may be worth improving the accuracy of the transport scheme, especially for water vapour and other chemically or radiatively active species.Regarding potential temperature, Skamarock and Gassmann (2011) have found that a third-order transport scheme for the potential temperature could significantly reduce phase errors in the propagation of baroclinic waves.
Whether more accurate transport of potential temperature is beneficial for climate modelling remains to be determined.
The Hamiltonian framework leaves a complete freedom with respect to the choice of a discrete Hamiltonian.Here the simplest possible second-order accurate approximation is used, but other forms may yield additional properties, such as a more accurate computation of the geopotential.Ongoing work suggests that it is possible to design a Hamiltonian such that certain hydrostatic equilibria are exactly preserved in the presence of arbitrary topography.Such a property is sometimes achieved by finite-volume schemes (Botta et al., 2004;Audusse et al., 2004), and its absence is one manifestation of the so-called pressure-gradient force error (Gary, 1973).
DYNAMICO is stabilized by (bi)harmonic operators of which we refer as filters rather than dissipation.Indeed, they are numerical devices aimed at stabilizing the model rather than physically based turbulence models such as nonlinear viscosity (Smagorinsky, 1963).Turbulence models induce a well-defined dissipation rate of resolved kinetic energy that should enter as a positive source term in the entropy budget in order to close the energy budget.Emulating this process in a discrete model can however prove difficult (Gassmann, 2013).Indeed, in order to convert into heat the kinetic energy destroyed by filters, one needs to recast their contribution to the energy budget as a positive-definite sum.Whether this can be done in DYNAMICO is left for future investigation.
Coupling DYNAMICO to the LMD-Z terrestrial physics package suite is ongoing.For planetary applications, it will be important to also check the discrete angular momentum budget (Lebonnois et al., 2012;Lauritzen et al., 2014a).
In the near future DYNAMICO should become able to solve richer, quasi-hydrostatic equations (White and Bromley, 1995;Tort and Dubos, 2014a) and to take into account deviations of the geopotential from spherical geometry (Tort and Dubos, 2014b).Extension to fully compressible Euler equations is the next step and should leave its general structure unchanged (Dubos and Tort, 2014).
demonstrating the sensitivity of Eq. (B2) to round-off error.Equation (B2) becomes useless when h ∼ ε/ h, which would happen with single-precision calculations at resolutions of about 1/1000 the planetary radius, i.e. 6 km on Earth.Conversely, Eq. ( B3) is stable and determines the position of p within round-off error.
Regarding the spherical centre of mass (Eq.A1), an exact Gauss formula exists for polygonal control volumes (not shown).Again this formula has large cancellation errors and yields C i with a round-off error O(ε/ h).For a spherical triangle, the planar centre of mass (equal-weight barycentre) projected onto the unit sphere yields a third-order accurate estimate of the true centre of mass.Therefore, subdividing a polygon into spherical triangles and forming an area-weighted sum of their barycentres yields a second-order accurate estimate of C i .This accuracy is sufficient for our purposes.An accurate and stable alternative is to decompose polygons into triangles and quadrangles, map the unit square to a spherical quadrangle or triangle use highorder Gauss-Legendre quadrature to evaluate Eq. (A1).
Finally, computing the area A of a spherical polygon should not be done using the simple but again unstable defect formula.Instead we decompose polygons into triangles and use l'Huillier formula: where A is the desired triangular area, 2s = a + b + c and a, b, c are the geodesic lengths of the sides of the triangle, computed as dist(p, q) = sin −1 p × q .Formula (B4) reduces for small triangles to the planar Henon formula, which demonstrates its stability.

Figure 1 .
Figure 1.Staggering and location of key prognostic and diagnostic variables.
but it would be possible to use more accurate, possibly upwind biased, reconstructions as in finite-volume advection schemes.Indeed, as shown byGassmann (2013), conserving discrete energy is possible when upwinding the advection of θ provided the same reconstructed values are reused in the curl-form momentum equation (see Sect. 4.2).

Figure 2 .
Figure 2. Latitude-altitude plot of advected tracer profile at the mid-time (t = 12 h, left) and in the end of the simulation (t = 24 h, right) for Hadley-like meridional circulation test case.Ideally, tracer isolines would be horizontal at the final time (contours separated by 0.1).

0, an icosahedral hydrostatic dynamical core changing
p ∞ only adds a constant to λ ik and does not change the motion (see Sect. 3.3).
. Stationary points of the pseudo-energy H = H − k k i µ ik − k π k i ik are such that ∂H /∂v ek = 0 (state of rest), ∂H /∂ ik = π ik = π k and ∂H /∂µ ik = i k = k .In the absence of topography, uniform i k and π ik in each layer are achieved if θ ik , µ ik , il do not depend on the horizontal position i.Such states of rest are stable provided H is convex.

Table 1 .
Global error norms for Hadley-like meridional circulation test case.Horizontal resolution is defined as 2R where 3 √ 3/2 10M 2 R 2 = 4πa 2 is the radius of the 10M 2 perfect and identical hexagons that would be needed to cover the surface 4πa 2 .