FVM 1 . 0 : a nonhydrostatic finite-volume dynamical core formulation for IFS

We present a nonhydrostatic finite-volume global atmospheric model formulation for numerical weather prediction with the Integrated Forecasting System (IFS) at ECMWF, and compare it to the established operational spectral-transform formulation. The novel Finite-Volume Module of IFS (henceforth IFS-FVM) integrates the fully compressible equations using semi-implicit time stepping and non-oscillatory forward-in-time (NFT) Eulerian advection, whereas the spectral-transform IFS solves the hy5 drostatic primitive equations (optionally the fully compressible equations) using a semi-implicit semi-Lagrangian scheme. The IFS-FVM complements the spectral-transform counterpart by means of the finite-volume discretisation with a local communication footprint, fully conservative and monotone advective transport, all-scale deep-atmosphere fully compressible equations in a generalised height-based vertical coordinate, applicable on flexible meshes. Nevertheless, both the finite-volume and spectral-transform formulations can share the same quasi-uniform horizontal grid with co-located arrangement of variables, 10 geospherical longitude-latitude coordinates, and physical parametrisations, thereby facilitating their comparison, coexistence and combination in IFS. We highlight the advanced semi-implicit NFT finite-volume integration of the fully compressible equations of the novel IFSFVM considering comprehensive moist-precipitating dynamics with coupling to the IFS cloud parametrisation by means of a generic interface. These developments—including a new horizontal-vertical split NFT MPDATA advective transport scheme, 15 variable time stepping, effective preconditioning of the elliptic Helmholtz solver in the semi-implicit scheme, and a computationally efficient coding implementation—provide a basis for the efficacy of IFS-FVM and its application in global numerical weather prediction. Here, numerical experiments focus on relevant dry and moist-precipitating baroclinic instability at various resolutions. We show that the presented semi-implicit NFT finite-volume integration scheme on co-located meshes of IFS-FVM can provide highly competitive solution quality and computational performance to the proven semi-implicit semi-Lagrangian 20 integration scheme of the spectral-transform IFS. 1 Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-237 Manuscript under review for journal Geosci. Model Dev. Discussion started: 11 October 2018 c © Author(s) 2018. CC BY 4.0 License.


Introduction
Notwithstanding the achievements made over the last decades (Bauer et al., 2015), numerical weather prediction (NWP) faces the formidable challenge to resolve rather than parametrise essential small-scale forcings and circulations in the multiscale global flow-most notably processes associated with the surface, convective clouds, gravity waves and tropospherestratosphere interaction.While there is a need for advancement in many aspects of global NWP model infrastructures, prerequisite are the ability of the numerical model formulations to accurately predict atmospheric flows throughout the large-scale hydrostatic and small-scale nonhydrostatic regimes, and to run efficiently on emerging and future high-performance computing architectures.
The spectral-transform (ST)-also known as pseudo-spectral-method was introduced in NWP following the work by Eliasen et al. (1970) and Orszag (1970).As a model representation entirely in spectral (i.e.wavenumber) space is impractical for NWP, the ST method maps between spectral and grid-point space, in order to solve different parts of the governing equations in the space where the computations can be performed most efficiently.Typically, nonlinear terms and the physical parametrisations are computed in grid-point space1 .Horizontal derivatives are computed in spectral space with formally high accuracy, as are linear terms of the discretised governing model equations-in particular the constant-coefficient Helmholtz problem resulting from the semi-implicit time stepping (Robert et al., 1972) can be solved directly and accurately in spectral space.Facilitated by the ST method, the unconditional stability of the semi-implicit scheme combined with semi-Lagrangian (SL) advection in grid-point space permits very long time steps2 and high efficiency (Ritchie et al., 1995;Temperton et al., 2001).In global models, the ST method typically uses a spherical harmonics representation in spectral space and (reduced, i.e. quasi-uniform) Gaussian grids (Hortal and Simmons, 1991;Wedi et al., 2015).
At ECMWF, the first forecast model using the ST method became operational in 1983, and the technique is still successfully applied today with the efficient SISL integration of the hydrostatic primitive equations in the Integrated Forecasting System (IFS) (Wedi et al., 2015).Recent advancements helped to sustain the performance of the ST method (for details see Wedi et al. (2013); Wedi (2014); Wedi et al. (2015) and Section 2.2) and enabled real-time medium-range global weather forecasts at ECMWF with ≈ 9 km horizontal grid spacings in 2016 (Malardel et al., 2016).Furthermore, current research advanced the applicability of the ST method into the realm of global convection-permitting forecasts with kilometer-scale horizontal grid spacings (Wedi and Düben, 2017).While the viability of the ST method at ECMWF is ensured for the next decade, uncertainties particularly concerning the scalability of the quasi-global data-rich parallel communications in the STs exist in the longer term (Wedi et al., 2013(Wedi et al., , 2015)).In addition, these scalability issues could be exacerbated from the time when the increase in horizontal resolution makes the nonhydrostatic formulation based on the fully compressible equations (Bubnová et al., 1995;Bénard et al., 2010) essential, because the associated solution procedure in IFS requires-at least in its current implementation-a predictor-corrector approach in the semi-implicit integration scheme which, among others, involves a considerably larger number of STs per model time step (Bénard et al., 2010;Wedi et al., 2015).
The uncertainties concerning the ST method with regard to emerging and future high-performance computing (HPC) architectures is one of the main reasons for ECMWF and its European partners to look into alternative nonhydrostatic, all-scale global model formulations and discretisation schemes to be incorporated in IFS.With this objective in mind, the Finite-Volume Module of IFS (henceforth IFS-FVM) is currently under development at ECMWF (Smolarkiewicz et al., 2016;Kühnlein and Smolarkiewicz, 2017;Smolarkiewicz et al., 2017).An important property of the finite-volume (FV) method applied in IFS-FVM is a compact spatial discretisation stencil in "grid-point" space, associated with a local distributed-memory communication footprint that contrasts with the global character of the STs.In addition, advantages of the FV method are inherently local conservation and the ability to operate in complex, unstructured-mesh geometries.The lack of conservation is a common issue with standard SL schemes and a shortcoming in the current operational IFS.Conservation errors are presumed to contribute to significant (moisture and temperature) biases in the upper-troposphere lower-stratosphere region (Wedi et al., 2015), and an albeit small but systematic drift in air mass and tracer fields may affect the forecast quality at longer (sub-)seasonal forecast ranges (Thuburn, 2008;Diamantakis and Flemming, 2014).The ability of FV methods to operate in complex, unstructured mesh geometries is of high relevance to global NWP.In the global (spherical or spheroidal) domains, the FV technique provides ample freedom for implementing efficient quasi-uniform resolution meshes, which circumvent the stiffness associated with the converging meridians towards the poles of the classical regular longitude-latitude grids commonly employed with finite-difference (FD) discretisation methods (Prusa et al., 2008;Wood et al., 2014).Flexibility with respect to the mesh is also important for implementing variable and/or adaptive resolution in atmospheric modeling systems, where locally finer mesh spacings in sensitive regions (e.g.along coastlines or in mountainous areas) may provide an efficient way towards a more accurate representation of multi-scale interactions (Bacon et al., 2000;Weller et al., 2010;Kühnlein et al., 2012;Zarzycki et al., 2014).
By default, IFS-FVM employs 3D semi-implicit integrators for the nonhydrostatic fully compressible equations (Smolarkiewicz et al., 2014).The all-scale integrators in IFS-FVM are conceptually akin to the semi-implicit schemes in IFS-ST but more general.In both formulations, accurate and robust integration with large time steps is achieved by 3D implicit representation of fast acoustic and buoyant modes supported by the fully compressible equations.Furthermore, fully implicit representation of slow rotational modes is another common feature of both IFS-FVM and IFS-ST (Temperton, 2011;Smolarkiewicz et al., 2016).Although implicit time stepping is predominantly associated with the computational stability, there are indications for favourable balance and accuracy in multi-scale flows (Knoll et al., 2003;Dörnbrack et al., 2005;Wedi and Smolarkiewicz, 2006).In contrast to IFS-ST, IFS-FVM's semi-implicit schemes do not rely on constant coefficients of the operator that is represented implicitly, as is required with the spectral space representation (Bénard et al., 2010).In IFS-ST, the operator that is represented implicitly results from linearisation of the full nonlinear governing equations about a horizontal reference state, and the semi-implicit integration then treats the nonlinear residual (i.e.full nonlinear equations minus the linear operator) explicitly (see Section 2.2).Constant coefficients effectively exclude orographic forcing from the linear operator of the semi-implicit scheme3 , leaving the associated effects solely to the explicit nonlinear residual.Furthermore, in the constantcoefficient semi-implicit scheme different (i.e.split) boundary conditions are applied in the linear operator and the nonlinear residual.Albeit still a research issue, the constant-coefficient semi-implicit scheme may incur reduced stability under more complex orography for future high-resolution forecasts in the nonhydrostatic regime.
At ECMWF, the reign of the ST method with the SISL integrators still continues, but future challenges especially with respect to HPC, nonhydrostatic modeling and complex orography can be foreseen.IFS-FVM represents an alternative dynamical core formulation that can complement IFS-ST with regard to these issues.However, to make IFS-FVM a useful option for global medium-range weather forecasting at ECMWF, it needs to be shown that the model formulation can provide (at least) comparable solution quality to the established IFS.In particular, a fundamental scientific question is whether a second-order FV method on the co-located meshes employed in IFS-FVM can sustain the accuracy of the ST method of IFS-ST.Another important question concerns the computational efficiency of IFS-FVM.At ECMWF and generally in NWP, tight constraints exist with regard to the runtime of the forecast models on the employed supercomputers.Therefore, we will evaluate the basic efficiency in terms of time-to-solution of the current IFS-FVM formulation relative to the operational hydrostatic IFS-ST and its nonhydrostatic extension.In the present paper, these issues are investigated using relevant atmospheric flow benchmarks such as those defined in the context of the Dynamical Core Model Intercomparison Project (DCMIP, Ullrich et al., 2017).
The DCMIP-2016 benchmarks involve large-scale hydrostatic and small-scale nonhydrostatic flows on the sphere, and also emphasise the interaction of the dynamical core with selected parametrisations of sub-grid scale physical processes.In the present paper IFS-FVM is verified against the proven IFS-ST at ECMWF for the baroclinic instability benchmark in the hydrostatic regime, considering specific configurations and parametrisations of interest at ECMWF.IFS-FVM also participates in the wider DCMIP-2016 model intercomparison, and this includes the nonhydrostatic supercell test case (Zarzycki et al., 2018).
The paper is organised as follows.Section 2 addresses the IFS-FVM and IFS-ST model formulations, and juxtaposes their main formulation features.In particular, while Section 2.1 provides a description of the advanced semi-implicit finite-volume integration scheme of the novel IFS-FVM, Section 2.2 summarises briefly the established IFS-ST.Furthermore, Section 2.3 discusses some basic aspects of the coupling to physics parametrisations, and Section 2.4 describes the common octahedral reduced Gaussian grid applied at ECMWF.Having described the model formulations, the IFS-FVM and IFS-ST benchmark simulations are then compared in Section 3. Section 4 concludes the paper.

IFS model formulations
The IFS comprises a comprehensive model infrastructure to perform data assimilation and to run deterministic and probabilistic global weather forecasts with various ranges and resolutions, supplemented with pre-and post-processing capabilities.The dynamical core lies at the heart of the NWP model infrastructure.
Table 1.Summary of the main formulation features of IFS-FVM and IFS-ST.For IFS-ST, information about the hydrostatic formulation and its nonhydrostatic extension are provided (see main text for description).Abbreviations are: finite-element (FE), finite-difference (FD), spectral-transform (ST), finite-volume (FV), semi-implicit (SI), two-time-level (2-TL).A summary of variables is provided in Table D1.

Finite-volume module of IFS
IFS-FVM solves the deep-atmosphere4 , non-hydrostatic, fully compressible equations with a generalised height-based terrainfollowing vertical coordinate.Numerical integration of the governing equations employs a centred two-time-level semi-implicit scheme that provides unconditional stability in 3D with respect to the fast acoustic and buoyant modes, as well as slower rotational modes (Smolarkiewicz et al., 2014(Smolarkiewicz et al., , 2016)).In Section 2.1.2,we extend the IFS-FVM semi-implicit integration to comprehensive moist-precipitating dynamics coupled to the IFS cloud physics parametrisations-this generalises the simplified moist-precipitating dynamics with different cloud physics coupling described in Smolarkiewicz et al. (2017).The IFS-FVM semi-implicit integration is combined with non-oscillatory forward-in-time (NFT) Eulerian advection based on MPDATA (Multidimensional Positive Definite Advection Transport Algorithm) (Kühnlein and Smolarkiewicz, 2017).In the present work, new efficient horizontal-vertical split NFT advective transport schemes based on MPDATA are developed and applied (Section 2.1.2and Appendix A).In addition, improved efficacy with the Eulerian NFT MPDATA advection is sought by rigorous implementation of variable time stepping.As reviewed in Section 2.1.3,the unstructured horizontal spatial discretisation uses the median-dual FV approach of Szmelter and Smolarkiewicz (2010), combined with a structured-grid FD/FV approach in the vertical direction (Smolarkiewicz et al., 2016).In both the horizontal and the vertical discretisation all dependent variables are co-located.Improved efficacy further results from advanced preconditioning of the elliptic Helmholtz solver in the semi-implicit scheme addressed in Section 2.1.2and Appendix B, and an efficient coding implementation with some aspects indicated in Section 2.1.3.To facilitate interoperability between the different numerical methods in IFS, IFS-FVM uses the median-dual FV mesh about the nodes of the octahedral reduced Gaussian grid.However, the IFS-FVM numerical formulation is not restricted to this grid and offers capabilities towards broad classes of meshes including adaptivity (Szmelter and Smolarkiewicz, 2010;Kühnlein et al., 2012).The IFS-FVM employs flexible parallel data structures provided by ECMWF's Atlas library (Deconinck et al., 2017).The main model formulation features of IFS-FVM are summarised in Table 1, and shown alongside the corresponding IFS-ST properties discussed below in Section 2.2.

Governing equations
Building on the formulation of moist-precipitating dynamics described in Smolarkiewicz et al. (2017), the fully compressible equations considered in IFS-FVM are given as which describe the conservation laws of dry mass (1a), momentum (1b), dry entropy (1c) and water substance (1d).Dependent variables in (1) are dry density ρ d , three-dimensional physical velocity vector u ≡ (u, v, w), potential temperature perturbation θ ′ , and a modified Exner pressure perturbation ϕ ′ ≡ c pd π ′ , as well as the water substance mixing ratios r k = ρ k /ρ d (i.e. the ratio of density of the individual water substance category ρ k to the density of dry air ρ d )-with the current cloud parametrisation of IFS (Forbes et al., 2010), five categories for water substance are considered (vapor r v , liquid r l , rain r r , ice r i , snow r s ), each described by the respective PDE (1d).An additional prognostic equation for cloud fraction Λ a employed with the IFS cloud parametrisation is implemented as Furthermore, the thermodynamic variables are related by the gas law (1e)-the Exner pressure is π = (p/p 0 ) R d /c pd and the potential temperature is θ = T /π, where p is the total pressure, p 0 ≡ 10 5 Pa, and T is the (absolute) temperature.A quantity which appears in various rhs terms of the momentum equation (1b) is the density potential temperature θ ρ = θ ϑ, where ϑ ≡ , 1994) with ε = R d /R v and the total water mixing ratio r t represents the sum over all the individual mixing ratios The multiplying factor ϑ appears explicitly in the buoyancy term of (1b), in order to expose the potential temperature perturbation θ ′ for implicit coupling to the thermodynamic equation (1c) (see Section 2.1.2).Note that a high-order approximation of the buoyancy term applied in Kurowski et al. (2014); Smolarkiewicz et al. (2017) for consistency with soundproof models at mesoscales, has been replaced here by the full, unabbreviated form essential for the accurate representation of moistprecipitating dynamics at planetary scales (Section 3).
Another important aspect of the governing equations ( 1) is the underlying perturbation form.The perturbations of potential temperature θ ′ = θ−θ a and Exner pressure π ′ = π−π a correspond to deviations from an ambient state (denoted by the subscript "a") that satisfies a general balanced subset of (1).A straightforward example applied in this paper is a stationary atmosphere θ a (x), u a (x), π a (x), r va (x) in thermal wind balance, which in terms of the momentum equation ( 1b) is where ϕ a ≡ c pd π a and θ ρa = θ a (1 + r va /ε)/(1 + r va ) .The perturbation form ( 1) is analytically equivalent to the corresponding unperturbed form of the fully compressible equations but has favourable properties for numerical integration; see Section 2.1.2and Prusa et al. (2008); Smolarkiewicz et al. (2014).Moreover, there are other analytically equivalent perturbation forms of the fully compressible equations implemented in IFS-FVM, which may differ in the degree of implicitness permitted in the integration (Smolarkiewicz et al., 2018).
The governing equations in (1)-( 2) are formulated with respect to generalised curvilinear coordinates embedded in a geospherical framework.At the most elementary level, the generalised curvilinear coordinate formulation can be used to implement fixed terrain-following levels with appropriate boundary conditions, but the model formulation optionally permits quite general moving meshes in the vertical and the horizontal directions.Symbols associated with the geometric aspects of the model are the transformed curvilinear coordinates x ≡ (x, y, z), the 3D nabla operator ∇ with respect to x, the Jacobian G of the coordinate transformations (i.e. the square root of the determinant of the metric tensor), a matrix of metric coefficients G, its transpose G T , and the contravariant velocity v = ẋ = G T u + v g where v g ≡ ∂x/∂t is the mesh velocity; see Prusa and Smolarkiewicz (2003); Kühnlein et al. (2012) for extended discussion.Following Smolarkiewicz et al. (2017), further symbols on the rhs of the momentum equation (1b) denote the gravity vector g = (0, 0, −g), the Coriolis parameter f ≡ 2Ω with Ω the angular velocity vector of the Earth's rotation, and M M M ′ subsumes the metric forces due to the curvature of the sphere Explicit expressions for M M M are given in Eq. (C3) of Appendix C, that also provides details about the specification of the curvilinear coordinate space in IFS-FVM.A summary of variables and physical constants is provided in Tables D1 and D2, respectively.Last but not least, the symbols P P P u ≡ (P u , P v , P w ), P θ ′ , P r k on the rhs of (1) denote the respective forcings from physical parametrisations.
Note that additional rhs terms not explicitly provided in the governing equations (1) may describe Rayleigh-type damping and/or Laplacian diffusion, applied especially to model wave-absorbing layers at the domain boundaries, see e.g.Prusa and Smolarkiewicz (2003); Klemp et al. (2008).

Semi-implicit numerical integration
To facilitate a compact description of the integration scheme, each of the governing equations of the system (1) is accommodated in a generalised conservation law form in which Ψ denotes the prognostic model variable, V the advector, G a generalised density, P Ψ represents the forcing from physical parametrisation and R Ψ the remaining right-hand-side;5 see Table 2 for the respective specifications of Ψ, V, G.The homogeneous mass continuity equation ( 1a) is a particular case of ( 6) and plays a fundamental role for the conservative advective transport of all other scalar variables.Note that because of the mass continuity equation (1a), the other scalar conservation laws ( 6) are equivalent to the Lagrangian form dΨ/dt = R Ψ + P Ψ , where d/dt = ∂/∂t + v • ∇ represents the total derivative.
Building on the earlier works by Smolarkiewicz et al. (2014Smolarkiewicz et al. ( , 2016)), the two-time-level numerical integrators of IFS-FVM for ( 6) can be subsumed in the following template scheme where The n, n+1/2, n+1 indices denote full and half time levels, and δt = t n+1 −t n the time step increment of the semi-implicit dynamics.The vector index i = (k, i) marks the node positions k and i of the, respectively, vertical and horizontal computational mesh, thereby revealing the 3D co-located spatial arrangement of all dependent variables underlying IFS-FVM's discretisation.
The definitions of the weights a Ψ and b Ψ for the incorporation of the right-hand-side terms R Ψ at t n and t n+1 , respectively, are given in Table 2. Apart from the incorporation of the physical parametrisation P Ψ , the semi-implicit scheme (7) with weights a Ψ ≡ b Ψ ≡ 0.5 is fully congruent with the second-order trapezoidal-rule trajectory integral of the corresponding ordinary differential equation of ( 6) (Smolarkiewicz and Margolin, 1993).In this paper and standardly in IFS-FVM, the operator 7) represents a flux-form Eulerian NFT advection scheme based on MPDATA methods (Smolarkiewicz and Szmelter, 2005;Kühnlein and Smolarkiewicz, 2017), but it could equally be a SL scheme (Smolarkiewicz et al., 2014).The basic MPDATA is second-order accurate given the advector V n+1/2 at the intermediate time level t n+1/2 is provided as an O(δt 2 ) estimate.
Note that due to the congruency of ( 7) with the ODE (9), the operator Ψ i may equally represent second-order semi-Lagrangian advection (Smolarkiewicz et al., 2014).In terms of the coupling to physics parametrisations, formally incorporated with firstorder accuracy, the associated forcing P Ψ | n = P Ψ (t phys , ∆t phys ) can optionally be evaluated with an equal or longer time step ∆t phys = N s δt (with integer N s = 1, 2, 3, ..) than δt applied in (7); see Section 2.3 about physics-dynamics coupling.Given the preceding discussion, the semi-implicit solution procedure of the governing system (1) proceeds from an atmospheric state at t n to a state at t n+1 as described in the following.The solution procedure commences with the integration of the mass continuity equation (1a) as which straightforwardly returns the updated density ρ n+1 d i .The O(δt 2 ) estimate for the advector (vG) n+1/2 in (10) is implemented here by linear extrapolation of the advective velocities from the previous time levels t n−1 and t n ; see Appendix A of Kühnlein et al. (2012) for the procedure accounting for a variable time step δt.In addition, the algorithm ( 10 port of all other scalar variables (see Table 2), a common approach to enable mass-compatible and monotonic solutions; see Smolarkiewicz et al. (2016) and Kühnlein and Smolarkiewicz (2017) for discussion in the context of IFS-FVM.
Given the tendencies from physical parametrisation P θ ′ , P u , P v , P r k , P Λa formally evaluated at t n , and the advective transport of all scalar variables θ ′ i , ≡ Λ ai -with the advected water content mixing ratios r ki and cloud fraction Λ ai representing already the final solutions at t n+1 -the scheme ( 7)-( 8) for the thermodynamic (1c) and momentum (1b) equations is implemented as where Due to the presence of nonlinear terms, the discrete system ( 11) is executed iteratively (Smolarkiewicz and Dörnbrack, 2008;Smolarkiewicz et al., 2014), and lagged quantities from the previous iteration are simply denoted by the superscript ⋆ .With the respective prognostic model variables θ ′ and u in (11) the n + 1 time level index has been dropped, but the n + 1 time level index is retained with the coefficient ϑ n+1 defined in (12) that is composed of the already updated r n+1 k .At the beginning of the iterative execution of (11), a first guess θ 0 for θ ⋆ is provided as the explicit solution of full potential temperature whereas a first guess u 0 for u ⋆ is prescribed by linear extrapolation of u from the t n−1 and t n time levels (again, taking into account the variable time step δt).
From (11), closed-form expressions for the discrete velocity update are derived by eliminating (11a) for θ ′ i in the buoyancy term of (11b)-thereby implementing 3D fully implicit treatment of buoyant modes in a moist-precipitating atmosphere-and gather the terms with linear dependence on u i on the lhs, which results in As all terms are co-located, the spatial mesh vector index i has been omitted in ( 14).The rhs of ( 14) is composed of all explicitly known terms, summarised as u, and the pressure gradient term with the lagged coefficient θ ⋆ ρ .Defining L as the linear operator acting on u on the lhs, and L −1 its inverse, ( 14) can be symbolised, respectively as where ǔ = L −1 u, and C = L −1 0.5 δt θ ⋆ ρ G denotes a 3 × 3 matrix of known coefficients.The solution algorithm presented up to ( 16), still requires to specify the Exner pressure perturbation ϕ ′ .A straightforward computation may employ the gas law (1e), using the updated variables ρ d , θ, and r v .However, the resulting 3D explicit acoustic integration is subject to very small time steps in order to maintain numerical stability, and inefficient for NWP.Therefore, a final step in IFS-FVM's numerical solution procedure is to augment the fully compressible equations ( 1) with an auxiliary 3D implicit boundary value problem for the pressure perturbation variable ϕ ′ (Smolarkiewicz et al., 2014).The formulation of this implicit boundary value problem originates from the advective (or Lagrangian) form of the gas law (1e), combined with the respective advective forms of the mass continuity (1a), thermodynamic (1c) and water vapour (1d) equations, see Smolarkiewicz et al. (2014Smolarkiewicz et al. ( , 2017) ) for further details.The governing equation for ϕ ′ is then implemented in conservation form consistent with (1), and reads where ϕ = ϕ a + ϕ ′ and Implicit coupling with the solution procedure of (1) is performed by inserting (16) into the three occurrences of u on the rhs of (17).Interpreting (17) in terms of the generalised transport equation ( 6) for Ψ ≡ ϕ ′ (cf.Table 2), while setting P ϕ ′ ≡ (R d /c vd ) ϕ Π and R ϕ ′ the remaining rhs consisting of the three divergence operators, the integration of ( 17) can be performed using a ϕ ′ = 0 and b ϕ ′ = 1.0 in the template scheme (7), i.e.Euler backward scheme with regard to R ϕ ′ .This is combined with a Euler forward scheme with regard to P ϕ ′ implemented as , δt).Reorganising terms finally yields the elliptic Helmholtz equation for the future pressure perturbation variable ϕ ′ , which can be written as (Smolarkiewicz et al., 2014(Smolarkiewicz et al., , 2017) ) where again due to the co-location of variables and terms, the spatial grid index i has been omitted.The summation ℓ in ( 19) is over the three divergence operators on the rhs of (17), respectively (7), while the coefficients A ⋆ ℓ , B ⋆ and ζ ℓ are defined accordingly.In IFS-FVM, the 3D boundary value problem ( 19) is solved iteratively using a nonsymmetric preconditioned Generalised Conjugate Residual (GCR) approach-see Smolarkiewicz and Szmelter (2011) for a recent discussion and Smolarkiewicz and Margolin (2000); Smolarkiewicz et al. (2004) for tutorials-while the dependence of the coefficients A ⋆ ℓ and B ⋆ on ϕ ′ in ( 19) is lagged behind.The preconditioning of GCR as applied in this paper employs a preconditioning operator P, that simplifies the linear operator of ( 19) by discarding the off-diagonal entries of the matrix G T C. Inversion of P then utilizes a weighted line Jacobi method, see Appendix B for details.As far as the adiabatic dynamics is concerned, the O(δt 2 ) integration of (17) with the Euler scheme maintains second-order accuracy of all variables except ϕ ′ over a single time step, because ϕ ′ enters the pressure gradient term of the momentum equation with the factor 0.5δt, hence resulting in an O(δt 3 ) integral. 6An improved first guess for ϕ ′ entering the elliptic problem ( 19) can be achieved by employing a Euler forward scheme based on (17).The 3D implicit scheme with respect to acoustic modes permits time steps equivalent to soundproof models (Kurowski et al., 2014;Smolarkiewicz et al., 2014;Kühnlein and Smolarkiewicz, 2017).Together with the full 3D implicit incorporation of buoyant modes described above, the semi-implicit integration is unconditionally stable with respect to all waves supported by the fully compressible equations (1), and thus the semi-implicit model time step δt can be selected according to the stability of the advective transport scheme A i in (7).

Spatial discretization
The discretisation framework of IFS-FVM combines a structured-grid FD method in the vertical with an unstructured7 FV approach in the horizontal (Smolarkiewicz et al., 2016).The FV discretisation and differentiation on the spherical surfaces adopts the median-dual approach described in Szmelter and Smolarkiewicz (2010).All dependent variables are co-located in the nodes in 3D.The consistent spatial discretisation of the applied MPDATA schemes in IFS-FVM, symbolised by the operator A i in (7), is described in Kühnlein and Smolarkiewicz (2017).
The schematic in Fig. 1 illustrates an arbitrary unstructured mesh on a 2D horizontal plane.The median-dual FV approach defines the control volume containing the node i by connecting the barycentres of polygonal mesh cells encompassing the node i with the midpoints of the edges originating in the node i.All geometric elements such as cell volume, cell face area, and face normals are evaluated from vector calculus in the computational space, i.e. in terms of x and y coordinates (see Section 2.1.1)on a zonally-periodic horizontal plane (Szmelter and Smolarkiewicz, 2010).The FV approach in the horizontal is combined with standard second-order FDs on the vertical structured grid with independent coordinate z. Figure 1.Schematic of the median-dual mesh in 2D.The edge connecting nodes i and j of the primary polygonal mesh pierces, precisely in the edge centre, the face Sj shared by computational dual cells surrounding nodes i and j.Open circles represent geometrical barycentres of the primary mesh, solid black lines mark the primary mesh and grey lines indicate dual cells with control volumes Vi and Vj, respectively.
For a differentiable vector field A, the Gauss divergence theorem Ω ∇ • A = ∂Ω A • n applied over the control volume surrounding node i = (k, i) in 3D is given as where l(i) numbers edges connecting node (k, i) with its horizontal neighbours (k, j), and S j refers both to the face per se and its surface area8 -the geometric quantities S j , V i and δz are independent of height.The fields A ⊥ k,j and A z k+1/2,i are interpreted as the mean normal components of the vector A at the horizontal and vertical cell faces, respectively.Basic approximations applied in IFS-FVM are given as A ⊥ k,j = 0.5 n j • [A k,i + A k,j ] in the horizontal, where n j is the respective mean outward unit normal to the face S j , and the vertical A z k+1/2,i = 0.5 A z k,i + A z k+1,i .For a differentiable scalar field Ψ, the 3D nabla operator ∇Ψ can also be interpreted in terms of the Gauss divergence theorem, leading to where S x j and S y j denote the x and y components of the oriented surface element S j = S j n j (Szmelter and Smolarkiewicz, 2010;Smolarkiewicz et al., 2016).Given ( 20) and ( 21) in the computational space, they are augmented with metrics of the curvilinear framework to obtain the respective physical divergence and gradient appearing in the governing equations of the previous Sections 2.1.1 and 2.1.2,respectively; details about the specification of the curvilinear space in IFS-FVM are reviewed in Appendix C.
For IFS-FVM, the mesh generation and mesh data structures as well as the nearest-neighbour distributed-memory communication, are handled by ECMWF's Atlas library, comprehensively described in Deconinck et al. (2017).Atlas is also designed to make use of specific programming models to support accelerators, although these have not been explored with IFS-FVM.For the quasi-uniform octahedral reduced Gaussian grid (see Section 2.4) the parallelisation adopts the equal-regions horizontal domain decomposition of IFS (Smolarkiewicz et al., 2016).20) and ( 21), outer loops over nodes of the co-located grid (re)compute the required quantities on edges "on the fly".The code based on the resulting longer, fused node-based loops perform a larger overall number of computations as quantities on the edges are computed more than once, but it is overall more efficient as many intermediate memory stores are avoided and there is a greater chance for data in cache to be reused.Due to the underlying FD/FV discretisation leading to stencil computations, IFS-FVM is naturally cache and memory-bandwidth bound.The node-based programming improves flop-per-byte ratio and relaxes the dependency on relatively (as compared to pure computation) slow memory access.As a consequence, it provides a significant gain in efficiency as compared to the hybrid edge-and node-based programming applied in Smolarkiewicz et al. (2016).In addition, the fused loops and data locality in the node-based code effectively enabled the performance of the sharedmemory parallelisation with OpenMP, as the threads may operate on more local memory regions and avoid thread conflicts and cache trashing.An overall speed-up of IFS-FVM by a factor of ∼3-4 has been found on ECMWF's Cray XC40 by converting from the hybrid edge-and node-based to the purely node-based code.

Spectral-transform IFS
The spectral-transform IFS (simply denoted as IFS-ST in this paper) operational at ECMWF is based on the hydrostatic primitive equations (HPEs) formulated in a hybrid η-pressure terrain-following vertical coordinate (Simmons and Burridge, 1981).The HPEs are integrated numerically using the efficient centred two-time-level SISL scheme, that permits long time steps due to the unconditional stability provided by the fully implicit treatment of the fast acoustic 9 and buoyant modes, and the 3D semi-Lagrangian advection (Ritchie et al., 1995;Temperton et al., 2001;Hortal, 2002).Therefore, the global, constant time step of the SISL integration in IFS-ST can be selected according to the optimal efficacy, rather than stability.The ST method, which is applied along model levels in the horizontal, is combined with a finite-element (FE) approach to discretize the integral operator in the vertical direction (Untch and Hortal, 2004)  10 .In 2002, this vertical FE scheme of Untch and Hortal (2004) with a co-located arrangement of prognostic variables replaced the former finite-difference scheme of Simmons and Burridge (1981) with the Lorenz staggering.Prognostic variables of the HPEs and other main formulation features of IFS-ST are given in Table 1.
The IFS-ST uses a discrete spherical-harmonics representation of the spectral space (Wedi et al., 2013).At every time step, the model fields are transposed between grid-point, Fourier and spherical-harmonics representation.General concerns about the computational efficiency of the Legendre transforms between Fourier modes and spherical-harmonics (Williamson, 2007), could be mitigated by adopting a fast Legendre transform (FLT, Wedi et al., 2013), which is employed together with fast Fourier transforms (FFTs).Furthermore, the increasing importance of nonlinearities of the rhs forcing terms in the governing 9 The HPEs analytically filter internal acoustic modes but support the external Lamb mode. 10 The FE implementation of the discrete vertical integral operator is based on the Galerkin method using cubic B-splines as basis functions.
equations and aliasing at higher resolutions stimulated the adoption of a cubic truncation of the spherical harmonics in the ST method (Wedi, 2014).The cubic versus the former linear truncation basically samples the highest wavenumber with four instead of two grid points.Special treatment is required near the poles with the reduced grids where the ratio approaches that of the linear truncation (Wedi, 2014).The extra sampling of a particular spectral resolution with a relatively larger grid size in the cubic truncation led to substantial further improvement in efficiency and accuracy of the IFS-ST at ECMWF (Wedi et al., 2015).The cubic truncation in combination with the octahedral reduced Gaussian grid went operational at ECMWF in 2016.
The nomenclature for this grid configuration is defined as "TCo" (for "Triangular-Cubic" truncation and "octahedral" grid) followed by the number of waves in spectral space (Malardel et al., 2016).The octahedral reduced Gaussian grid (or simply "octahedral grid"), which is also employed with IFS-FVM, is reviewed in Section 2.4.
The semi-Lagrangian advection scheme in IFS-ST is based on Ritchie et al. (1995).The SL trajectories are computed with an iterative algorithm.At each iteration of the algorithm, the wind at the midpoint in time and space is re-evaluated using the second order time-extrapolating algorithm SETTLS (Stable Extrapolation Two-Time-Level Scheme, Hortal, 2002).The two-time-level semi-implicit integration of IFS-ST follows the template Here, X represents the prognostic variables, the subscripts A, M and D denote, respectively, the arrival, mid-and departure points, while n and n + 1 again refer to the current and future time step.In addition, the operator L symbolises the linear is interpolated at the departure point by quasi-cubic interpolation.A horizontal quasi-monotonic interpolation is used for the horizontal wind, the temperature and the surface pressure and a 3D quasi-monotonic limiter is used for the specific humidity.
All the other water content variables are estimated at the departure point using linear interpolation and thus no monotonic filter is needed.The SL scheme in IFS-ST is applied to specific (per unit mass) variables whose equations are written in advective form.It is then intrinsically non-conservative (Malardel and Ricard, 2015), but global mass-fixers have been developed for the atmospheric composition applications (Diamantakis and Flemming, 2014;Diamantakis and Augusti-Panareda, 2017).Notably, with the cubic truncation and the octahedral reduced Gaussian grid, global mass conservation is nearly exact in IFS, see Wedi et al. (2015).
The IFS also includes various research options which are not yet applied in the operational configuration.Most notably, the nonhydrostatic (NH) formulation based on the fully compressible equations (Bubnová et al., 1995;Bénard et al., 2010), which has been made available to ECMWF by Météo-France and the Aladin Consortium.The HPE and NH formulations of IFS-ST employ the same SISL integrators, but the NH extension requires a predictor-corrector approach-the so-called iterative-centred-implicit scheme (Bénard et al., 2010)-for stability in global configurations.Prognostic variables and main characteristics of the NH formulation are also given in Table 1.Furthermore, currently the NH formulation is restricted to the vertical FD scheme by Bubnová et al. (1995).The HPE and NH formulations of IFS-ST standardly use the shallow- atmosphere approximation.In Section 3, we will compare the computational performance of the HPE and NH formulations of IFS-ST against IFS-FVM.

Some aspects about physics-dynamics coupling
The IFS physics parametrisation package at ECMWF is applied in the same configuration throughout the medium-range, sub-seasonal and seasonal forecasting systems.The physics package includes parametrisation of radiation, moist convection, clouds and stratiform precipitation, surface processes, sub-grid scale turbulence, as well as orographic and non-orographic gravity wave drag.
In IFS-ST, the physics-dynamics coupling employs the SLAVEPP (Semi-Lagrangian Averaging of Physical Parametrisation) scheme (Wedi, 1999).SLAVEPP targets second-order accuracy by averaging tendencies from selected physics parametrisations along the SL trajectory, at the midpoint in space-time between the departure point at t n and arrival point at t n+1 .11Thereby, the tendencies at t n+1 use a provisional first guess from the explicit dynamics as input, and the final tendencies are applied at the arrival point of the SL trajectory, see Wedi (1999) for details.The basic approach of the physics-dynamics coupling in IFS uses sequential splitting (or fractional stepping), i.e. the various processes are integrated one after another, and the updated tendencies are used as input to the subsequent process, see Beljaars et al. (2018) for discussion.A more comprehensive explanation of the IFS physics parametrisations can be found in the general IFS documentation.
The physics-dynamics coupling in IFS-FVM differs from IFS-ST.As explained in Section 2.1.2,the current implementation of IFS-FVM incorporates the tendencies from physics parametrisations P Ψ by means of a first-order coupling at t n .Therefore, the fields that enter the parametrisations are from t n , and there is no averaging between tendencies from t n and t n+1 .While the incorporation of the tendencies from physics parametrisation deviates from IFS-ST, the sequential splitting between the various IFS physics parametrisations is kept exactly the same.Incorporating the physics parametrisations with first order at t n is motivated by the generally smaller time steps δt in IFS-FVM than in IFS-ST, and the desire to implement straightforward options for subcycling of the dynamics (see below).In addition, numerical experimentation with IFS-FVM so far has shown favourable results in terms of the incorporation of the physics at t n .Nevertheless, different forms of coupling IFS-FVM to the physical parametrizations will be explored in the future.
The IFS-FVM code has its own interface to the IFS physics parametrisations.Among others, it involves conversion between IFS-FVM's variables and those employed in IFS-ST (see Table 1) 12 , interpolations to vertical interfaces, and provision of local quantities describing the mesh geometry, but also a number of technical aspects due to some differences in the computational design of IFS-FVM and IFS-ST.Some of these operations in the interface may be removed at a later stage when IFS physics parametrisation package becomes more harmonised with IFS-FVM.However, generally the coupling is facilitated by common features of IFS-FVM and IFS-ST, such as longitude-latitude coordinates, the octahedral grid, and co-located arrangement of variables.The IFS-FVM interface to the physics parametrisations also includes an option for subcycling of the dynamics.The template scheme (7) for one physics time step from t N to t N + ∆t phys ≡ t N + N s δt, can be written as ℓ = 1, N s : where The physics tendency P Ψ is evaluated with the physics time step ∆t phys and is then reused for the N s subcycling steps with δt.
IFS-FVM has been coded rigorously with a variable time stepping capability (Kühnlein and Smolarkiewicz, 2017).In case of the physics-dynamics coupling with subcycling (23), we adapt the semi-implicit time step δt only every N s time steps with the corresponding physics time step given as ∆t phys = N s δt.Apart from radiation13 , the current operational IFS-ST evaluates the physics parametrisations at every semi-implicit time step δt.However, due to the Eulerian versus SL advection, IFS-FVM uses considerably smaller time steps than IFS-ST,14 and therefore physics-dynamics coupling may use some form of subcycling as indicated above, or other approaches such as parallel splitting, to remain competitive.The cost per time step of the full IFS physics parametrisation package can be up to 40 % of the forecast model.For the idealized DCMIP experiments considered in the present work, we focus on the parametrisation of clouds and stratiform precipitation, incorporated by means of the general IFS-FVM and IFS-ST interfaces described above.

Octahedral reduced Gaussian grid
As with the classical reduced Gaussian grid of Hortal and Simmons (1991), the octahedral reduced Gaussian grid (or simply "octahedral grid") (Wedi et al., 2015;Malardel et al., 2016;Smolarkiewicz et al., 2016) specifies the latitudes according to the roots of the Legendre polynomials.The two grids differ in the arrangement of the points along the latitudes, which follows a simple rule for the octahedral grid: starting with ∼ 20 points on the first latitude around the poles, four points are added with every latitude towards the equator, whereby the spacing between points along the individual latitudes is uniform and there are no points at the equator.The octahedral grid is suitable for transformations involving spherical harmonics, and has been introduced for operational weather prediction with IFS at ECMWF in 2016.Figure 2 depicts the 'O24' octahedral grid nodes, together with the edges of the primary mesh and dual resolution as applied in IFS-FVM.Compared to the classical reduced Gaussian grid of Hortal and Simmons (1991), the octahedral grid provides a much more uniform dual mesh resolution in the FV context (Malardel et al., 2016).Negligible grid imprinting in IFS-FVM with the octahedral grid will be shown by means of numerical experiments in Section 3. In contrast, IFS-FVM generally applies variable time stepping that targets a maximum horizontal advective Courant number of 0.95-the actual maximum 3D Courant number is significantly larger than 1 as this is permitted by the horizontal-vertical split NFT advection in IFS-FVM (Appendix A).All IFS-FVM and IFS-ST results presented in this section were obtained without any explicit diffusion or regularization.
For the dry and moist configurations, the baroclinic instability evolution starts from two zonal jet flows in the mid-latitudes of each global hemisphere, that are in thermal wind balance with the meridional temperature gradient.A local zonal velocity perturbation in the form of a simple exponential bell (tapered to zero in the vertical) excites the instability, leading to eastward propagating Rossby modes (Ullrich et al., 2016).Here, we apply the triggering zonal velocity perturbation in both the northern and southern global hemispheres.This dual triggering departs from (Ullrich et al., 2016), where the perturbation was only applied in the northern hemisphere, but it permits a clean evaluation of kinetic energy spectra relatively early in the baroclinic wave evolution, enables study of the solution symmetry about the equator and is more relevant to real weather.After an initial period of linear growth, the instability enters the nonlinear stage from 6-7 days of simulation time.Our analysis will focus on simulation results at day 10 and day 15.

Results for dry simulations
Figures 3 and 4 present the results for the dry adiabatic simulations at day 10, when a large amplitude baroclinic wave has developed and formed sharp fronts in the lower troposphere.Generally very close agreement is found between the finitevolume (IFS-FVM) and spectral-transform (IFS-ST) solutions.This is emphasised by the difference plots in the bottom row of the horizontal and vertical cross sections in Figs. 3 and 4, respectively.The solutions show identical phase propagation and amplitude of the baroclinic wave throughout the entire vertical depth of the simulation domain, which applies to the very coarse (O160,TCo159) and coarse (O320,TCo319) grids.Where present, differences between IFS-FVM and IFS-ST become smaller with the higher resolution.Figure 5 compares the pressure field much later into the nonlinear baroclinic instability 30 evolution at day 15, for the finer of the two grids.This depiction shows close agreement between IFS-FVM and IFS-ST also at this later stage.There are no signs of any significant grid imprinting, and both dynamical core formulations provide visibly symmetric solutions around the equator with the octahedral grid.Kinetic energy spectra evaluated at day 15 in Fig. 6

Simulation results for moist-precipitating configuration with IFS cloud parametrisation
Next we present results for the moist-precipitating baroclinic instability with coupling to the IFS cloud parametrisation.Fig- ure 7 shows the instantaneous large-scale precipitation rate at the surface 16 for the (O160,TCo159) and (O320,TCo319) grids at day 10.For any of these grids, both model formulations show five rainbands with essentially identical phase, as emphasized by the overlay with the 0.5 mm/day black contour line of the corresponding other model formulation.The elongated rainbands are associated with the lifting along sharp frontal zones.Precipitation amounts are overall similar but somewhat higher local values exist for IFS-FVM, particularly in the two easternmost rainbands when looking at the (O160,TCo159) grid.Figure 8 is analogous to Fig. 7 but for day 15.As can be expected, the spread between the different model formulations becomes larger.
However, there is still reasonably close agreement, especially for the higher-resolution grid (O320,TCo319) in the right column of Fig. 8. Here, the location of the easternmost frontal zone and associated rainband agrees closely considering the late stage of the baroclinic instability evolution.Figure 9 supplements the precipitation plots with the corresponding pressure field on day 15.In addition to to the standard configurations of IFS-FVM and IFS-ST where the physics parametrisation is evaluated every dynamics time step N s = 1, Fig. 9 also provides the IFS-FVM result with subcyling (middle panel) where the parametrisations are evaluated every N s = 3 semi-implicit time steps δt; see 2.3 for discussion of the physics-dynamics coupling.Again, the pressure fields of all three simulations resemble each other closely, often even in the location and magnitude of smaller structures, while the modified physics-dynamics coupling frequency N s = 3 to the cloud parametrisation seems to have only a small impact on the solution.Furthermore, none of the simulations shows significant grid imprinting in the pressure fields, but the solution symmetry about the equator is broken in both IFS-FVM and IFS-ST as a result of the incorporation of the 16 The precipitation rate represents the liquid and rain (excluding ice and snow) sedimentation flux at the surface.cloud parametrisation (in contrast to the dry results shown before in Fig. 5).The analysis of the simulations is supplemented in Compared to the spectra of the dry simulations in Fig. 6, IFS-FVM and IFS-ST consistently show a considerably larger kinetic energy in the scales smaller than wave number ≈ 120 in the mid-troposphere at about 500 hPa.Overall, the presented consistent results of IFS-FVM and IFS-ST attest to the quality of the presented dry and moist-precipitating FV formulations along with the coupling to the IFS physics parametrisation.The proposed scheme implements mass-compatible second-order Strang-splitting as explained in the following.The overall semi-implicit integration of the fully compressible equations (1) proceeds exactly as explained in Section 2.1.2,but with the 3D NFT advection operator A i split into purely horizontal A xy i and vertical A z i schemes, respectively.For each model time step δt, these are applied in the sequence A z i → A xy i → A z i using half time steps in the two vertical sweeps and the full time step in the horizontal part.Specifically, the split scheme commences with the integration of the mass continuity equation (1a) as , G [2] , δt) , d , (v z G) n+1/2 , G [2] , G n+1 , 0.5δt) , which provides the updated densities ρ d , ρ n+1 d and accumulates normal mass fluxes (v z Gρ d ) [1] , (v ⊥ h Gρ d ) [2] , (v z Gρ d ) [3]  for the three sub-steps.For compatibility with mass continuity, these quantities are then all employed in the subsequent advective transport of scalar variables Ψ (8) as i ( Ψ, (v z Gρ d ) [1] , (Gρ d ) n , (Gρ d ) [1] , 0.5δt) , (A2)
In (A1) and (A2), the implementation of the horizontal advection transport A xy follows the horizontal part of the unstructuredmesh FV MPDATA of Kühnlein and Smolarkiewicz (2017).The vertical scheme A z is a corresponding 1D structured-grid MPDATA.Results from numerical experimentation relevant to NWP show that the presented horizontally-vertically split NFT scheme based on MPDATA can be considerably more efficient than the standard fully multidimensional (unsplit) MPDATA of Kühnlein and Smolarkiewicz (2017).This is particularly due to the integration of the vertical parts A z i with δt/2 each, which mitigates the vertical stability restriction while not adding any significant computational cost18 .Overall, the horizontalvertical splitting of A i can enable a more than twice larger time step in the integration than the unsplit formulation.In addition, the split scheme facilitates the application of higher-order, e.g.Waruszewski et al. (2018), and/or flux-form semi-Lagrangian advective transport in the vertical.While a detailed presentation and analysis will be provided in a future publication, results so far indicate comparable solution quality of the split versus unsplit schemes for global atmospheric flow benchmarks.All IFS-FVM results presented in this paper were obtained using the split scheme (A1)-(A2) for A i in ( 7)-(8).
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.In IFS-FVM, there two alternative implementations of the NFT advection transport operator A i .Standard formulations of MPDATA for integrating the fully compressible equations in the horizontally-unstructured vertically-structured discretisation framework of IFS-FVM are described inKühnlein and Smolarkiewicz (2017).These MPDATA formulations are fully multidimensional (i.e.unsplit), equipped with non-oscillatory enhancement(Smolarkiewicz and Grabowski, 1990;Smolarkiewicz and Szmelter, 2005), and qualify for implicit large-eddy simulation (ILES) of high-Reynolds number atmospheric flows, e.g.Domaradzki et al. (2003);Piotrowski et al. (2009);Smolarkiewicz et al. (2013).In addition to the standard, fully multidimensional advection schemes for A i (7), a more efficient horizontal-vertical split implementation of MPDATA in IFS-FVM has been developed and it is outlined in Appendix A. All IFS-FVM results presented in Section 3 of this paper were obtained with this horizontal-vertical split scheme.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.In the present work, the programming of the discrete differential operators (20) and (21) in the IFS-FVM modern Fortran code has been augmented from earlier implementations in Smolarkiewicz et al. (2016).While Smolarkiewicz et al. (2016) used a hybrid edge-and node-based programming, already different from the edge-based codes described in Szmelter and Smolarkiewicz (2010), the current IFS-FVM programming represents effectively a node-based implementation.When evaluating ( operator and N the nonlinear residual N = M − L, with M being the full nonlinear model.The operator L results from linearisation of M with regard to a horizontally-homogeneous reference state.Before the final semi-implicit solution with (22) at each time step, the explicit guess of the future model state X is computed by replacing L n+1 A with L n A .Generally, X n D Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.We use two different sizes of the octahedral reduced Gaussian grid for comparing the solution quality of IFS-FVM and IFS-ST.Considered are very coarse (O160,TCo159) and coarse (O320,TCo319) grids by current NWP standards, corresponding to about 64 and 32 km nominal horizontal grid spacings, respectively; see Section 2.2 and 2.4 for the nomenclature that defines the grids.For the two different horizontal grid sizes, both IFS-FVM and IFS-ST employ 60 stretched vertical levels with a similar height distribution.The vertical spacing in terms of height between model levels ranges from about 12 meters near the ground to 4 km near to model top located at 48 km.For the comparison of computational efficiency in Section 3.3, we employ the (O1280,TCo1279) grid, corresponding to about 9 km spacing, with 137 stretched vertical levels.With regard to the time step size, IFS-ST uses constant increments δt of 1800, 1200, and 450 s for the TCo159, TCo319, and TCo1279 grids, respectively.
reveal a strikingly similar distribution of variance across wave numbers from 1 to ∼ 200.As all the simulations are without any explicit diffusion, the IFS-FVM spectra at the high wave numbers attest to the implicit scale-selective regularisation with artificial viscosity Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.provided by the nonoscillatory finite-volume MPDATA advection 15 .The slope of the IFS-ST and IFS-FVM spectra with respect to wavenumber l is somewhat shallower than l −3 at large scales.This is consistent with results from other dynamical cores for this baroclinic instability test case studied in the context of the High-Impact Weather Prediction Project (HIWPP,Whitaker, 2014).The spectra of IFS-ST feature some accumulation of energy near the truncation scale of the triangular-cubic grid, corresponding to four times the grid spacing.This increased energy near the truncation scale does not grow and is to a 5 large extent controlled by the spectral filtering at every time step of the nonlinear terms with the triangular-cubic grid, among other mechanisms of implicit dissipation such as the semi-Lagrangian interpolation.

Figure 3 .
Figure3.Baroclinic instability at day 10: the first two rows show pressure on the lowest full level (hPa) obtained with IFS-FVM and IFS-ST, while the last row depicts the corresponding difference between the solutions.The left and right columns are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.15The unsplit NFT MPDATA advection(Kühnlein and Smolarkiewicz, 2017) features lower implicit diffusion near the grid scale than the split advection applied here.

Figure 4 .
Figure 4. Dry baroclinic instability at day 10: the first two rows show meridional wind v (m s −1 ) along a zonal-height cross section at 50 • N obtained with IFS-FVM and IFS-ST, while the last row depicts the difference of v (m s −1 ) between their solutions.The left and right columns are for the (O160,TCo159) and intermediate (O320,TCo319) horizontal grids, respectively.

Fig. 10
Fig. 10 with the time series of the minimum near-surface pressure (left panel) and the area-integrated precipitation rate (right panel).The temporal evolution of these two quantities is close between IFS-FVM and IFS-ST.Particularly, the minimum nearsurface pressure agrees almost exactly at day 15, although small differences occur over the course of the simulation.The onset and subsequent increase of the precipitation matches well in IFS-FVM and IFS-ST, and the later variations in the precipitation rate are similar, with no systematic under-or overestimation.Kinetic energy spectra evaluated at day 15 are shown in Fig. 11.

Figure 7 .
Figure 7. Moist-precipitating baroclinic instability at day 10: surface precipitation rate (mm/day) obtained with IFS-FVM and IFS-ST coupled to the same IFS cloud microphysics parametrisation.The upper (lower) row shows shaded contours from the IFS-FVM (IFS-ST) simulations, overlaid by the IFS-ST (IFS-FVM) black contour line of 0.5 mm/day.The left and right columns are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

Figure 8 .
Figure 8. Moist-precipitating baroclinic instability at day 15: surface precipitation rate (mm/day) obtained with IFS-FVM and IFS-ST coupled to the same IFS cloud microphysics parametrisation.The upper (lower) row shows shaded contours from the IFS-FVM (IFS-ST) simulations, overlaid by the IFS-ST (IFS-FVM) black contour line of 0.5 mm/day.The left and right columns are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

Figure 9 .
Figure9.Moist-precipitating baroclinic instability at day 15: pressure on the lowest full level (hPa) obtained with IFS-FVM and IFS-ST coupled to the same IFS cloud microphysics parametrisation.The IFS-FVM results in the first and second row employed the standard coupling at every time step Ns = 1 or subcycling of dynamics for three time steps Ns = 3, respectively.The simulations were performed with the (O320,TCo319) grid.

Figure 10 .
Figure 10.Moist-precipitating baroclinic instability: Time series of minimum pressure on the lowest full level (left) and area-integrated rain rate (right).The blue and red lines correspond to the IFS-FVM and IFS-ST results, respectively.

Figure 11 .
Figure 11.Moist-precipitating baroclinic instability at day 15: Kinetic energy spectra obtained with IFS-FVM and IFS-ST using the (O320,TCo319) grid.The blue vertical line indicates the spatial scale corresponding to four times the nominal grid spacing of the O320 octahedral grid, which also represents the cubic truncation scale with TCo319 applied in IFS-ST.The spectra are shown on model levels near the surface and at ∼ 500 hPa.

3. 3
Computational efficiencyComputational efficiency of NWP models is crucial.For current HPC architectures and model resolutions, the operational IFS-ST at ECMWF represents one of the most efficient dynamical core formulations for global NWP.Although the novel IFS-FVM is considered for future applications in the nonhydrostatic regime and potentially different requirements in terms of the computing hardware, an important basic question concerns its performance for current NWP configurations.Of interest is the 5 relative performance of IFS-FVM to both the hydrostatic IFS-ST and its nonhydrostatic extension (see Section 2.2).In order to Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.modates to IFS-FVM's unstructured horizontal discretisation enabling broad classes of meshes over the surface of the Earth's sphere/spheroid, and the structured grid in the (stiff) vertical direction.

Table 2 .
Specification of prognostic model variables and corresponding parameters in the template scheme (7)-(8).Columns represent the dependent variable Ψ, the advector V, a generalised density G, and a Ψ , b Ψ the weights for the incorporation of the rhs forcings R Ψ .

Table D2 .
List of physical constants Specific heat capacity of dry air at constant pressure 1004.5 J kg −1 K −1 c vd Specific heat capacity of dry air at constant volume 717.5 J kg−1 K −1 R d Gas constant for dry air 287.0 J kg −1 K −1 R v Gas constant for water vapor 461.5 J kg −1 K −1 ε Ratio of R d to R v 0.622 supported in part by funding received from the European Research Council under the European Union's Seventh Framework Programme FP7/2012/ERC (grant agreement no.320375), and the Horizon 2020 Research and Innovation Programme ESCAPE (grant agreement no.671627) and ESIWACE (grant agreement no.675191).Rupert Klein thanks ECMWF for support under their Fellowship program, and acknowledges partial support of his contributions by Deutsche Forschungsgemeinschaft, Grant CRC 1114 "Scaling Cascades in Complex Systems", Project A02.Zbigniew Piotrowski acknowledges partial support from the "Numerical weather prediction for sustainable Europe" 10 project, carried out within the FIRST TEAM programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund. Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-237Manuscript under review for journal Geosci.Model Dev. Discussion started: 11 October 2018 c Author(s) 2018.CC BY 4.0 License.