FVM 1.0: a nonhydrostatic finite-volume dynamical core for the 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 the 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 hydrostatic 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 discretization with a local low-volume communication footprint, fully conservative and monotone advective transport, all-scale deep-atmosphere fully compressible equations in a generalized height-based vertical coordinate, and flexible horizontal meshes. Nevertheless, both the finite-volume and spectral-transform formulations can share the same quasi-uniform horizontal grid with co-located arrangement of variables, geospherical longitude– latitude coordinates, and physics parameterizations, thereby facilitating their comparison, coexistence, and combination in the IFS. We highlight the advanced semi-implicit NFT finitevolume integration of the fully compressible equations of IFS-FVM considering comprehensive moist-precipitating dynamics with coupling to the IFS cloud parameterization by means of a generic interface. These developments – including a new horizontal–vertical split NFT MPDATA advective transport scheme, variable time stepping, effective preconditioning of the elliptic Helmholtz solver in the semi-implicit scheme, and a computationally efficient implementation of the median-dual finite-volume approach – 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 semiimplicit NFT finite-volume integration scheme on co-located meshes of IFS-FVM can provide highly competitive solution quality and computational performance to the proven semiimplicit semi-Lagrangian integration scheme of the spectraltransform IFS.

Abstract. 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 the 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 hydrostatic 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 discretization with a local low-volume communication footprint, fully conservative and monotone advective transport, all-scale deep-atmosphere fully compressible equations in a generalized height-based vertical coordinate, and flexible horizontal meshes. Nevertheless, both the finite-volume and spectral-transform formulations can share the same quasi-uniform horizontal grid with co-located arrangement of variables, geospherical longitudelatitude coordinates, and physics parameterizations, thereby facilitating their comparison, coexistence, and combination in the IFS.
We highlight the advanced semi-implicit NFT finitevolume integration of the fully compressible equations of IFS-FVM considering comprehensive moist-precipitating dynamics with coupling to the IFS cloud parameterization by means of a generic interface. These developments -including a new horizontal-vertical split NFT MPDATA advective transport scheme, variable time stepping, effective precondi-tioning of the elliptic Helmholtz solver in the semi-implicit scheme, and a computationally efficient implementation of the median-dual finite-volume approach -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 semiimplicit NFT finite-volume integration scheme on co-located meshes of IFS-FVM can provide highly competitive solution quality and computational performance to the proven semiimplicit semi-Lagrangian integration scheme of the spectraltransform IFS.

Introduction
Notwithstanding the achievements made over the last decades , numerical weather prediction (NWP) faces the formidable challenge of resolving rather than parameterizing essential small-scale forcings and circulations in the multi-scale global flow -most notably processes associated with the surface, convective clouds, gravity waves, and troposphere-stratosphere interaction. While there is a need for advancement in many aspects of global NWP model infrastructures, prerequisites 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 (HPC) architectures.
Published by Copernicus Publications on behalf of the European Geosciences Union. The spectral-transform (ST) -also known as pseudospectral -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 physics parameterizations are computed in grid-point space. Horizontal derivatives are computed in spectral space with formally high accuracy, as are linear terms of the discretized 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 steps 1 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) . Recent advances helped to sustain the performance of the ST method (for details see Wedi et al., 2013Wedi et al., , 2015Wedi, 2014 and Sect. 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 kilometre-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 concerning the scalability of the non-local high-volume parallel communications in the STs exist in the longer term (Wedi et al., 2013. 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. This is because the associated solution procedure in the IFS requires -at least in its current implementation -a predictorcorrector approach in the semi-implicit integration scheme that involves a considerably larger number of STs per model time step (Bénard et al., 2010;Wedi et al., 2015). Furthermore, the large overlap regions between parallel distributedmemory partitions required in the SL scheme of IFS-ST can also be an issue in the longer term.
The uncertainties concerning the SISL integration based on the ST method with regard to emerging and future 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 discretization schemes to be incorporated in the IFS. With this objective in mind, the Finite-Volume Module of the IFS (henceforth IFS-FVM) is under development at ECMWF ). An important property of the finite-volume (FV) method applied in IFS-FVM is a compact spatial discretization stencil in "grid-point" space, associated with a distributed-memory communication footprint that is predominantly local and performed using thin overlap regions with the nearest neighbours, in contrast to the non-local high-volume communications required in IFS-ST. 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 , and a 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, unstructuredmesh 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 that circumvent the polar anisotropy of the classical regular longitude-latitude grids commonly employed with finite-difference (FD) discretization 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 modelling systems, in which locally finer mesh spacings in sensitive regions (e.g. storm tracks) 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 3-D semi-implicit integrators for the nonhydrostatic fully compressible equations . 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 3-D 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 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 linearization of the full non-linear governing equations about a horizontal reference state, and the semi-implicit integration then treats the non-linear residual (i.e. full non-linear equations minus the linear operator) explicitly (see Sect. 2.2). Constant coefficients effectively exclude orographic forcing from the linear operator of the semi-implicit scheme 2 , leaving the associated effects solely to the explicit non-linear residual. Furthermore, in the constant-coefficient semi-implicit scheme different (i.e. split) boundary conditions are applied in the linear operator and the non-linear residual. Although 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 modelling, and complex orography, can be foreseen. The 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 the 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 emphasize the interaction of the dy-2 Including orography in the implicit part involves multiplications which are standardly performed in grid-point space in the context of the ST method. In principle, one could carry out the necessary multiplications in spectral space but this is usually avoided because of computational complexity. namical core with selected parameterizations of sub-gridscale 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 parameterizations 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 organized as follows. Section 2 addresses the IFS-FVM and IFS-ST model formulations and juxtaposes their main formulation features. In particular, while Sect. 2.1 provides a description of the advanced semi-implicit finitevolume integration scheme of the novel IFS-FVM, Sect. 2.2 briefly summarizes the established IFS-ST. Furthermore, Sect. 2.3 discusses some basic aspects of the coupling to physics parameterizations, and Sect. 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 Sect. 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 preprocessing and postprocessing capabilities. The dynamical core lies at the heart of the NWP model infrastructure.

Finite-Volume Module of the IFS
IFS-FVM solves the deep-atmosphere 3 , nonhydrostatic, fully compressible equations with a generalized height-based terrain-following vertical coordinate. Numerical integration of the governing equations employs a centred two-time-level semi-implicit scheme that provides unconditional stability in 3-D with respect to the fast acoustic and buoyant modes, as well as slower rotational modes . In Sect. 2.1.2, we extend the IFS-FVM semiimplicit integration to comprehensive moist-precipitating dynamics coupled to the IFS cloud physics parameterizations -this generalizes 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)    (2010) combined with a structured-grid FD-FV approach in the vertical direction . In Sect. 2.1.3, we present a revised computationally efficient implementation of the median-dual FV approach. High efficacy also results from effective preconditioning of the Krylov-subspace solver for the Helmholtz problem arising in the semi-implicit time stepping of the fully compressible equations addressed in Sect. 2.1.2 and Appendix B. To facilitate interoperability between the different numerical methods in the IFS, IFS-FVM uses the median-dual FV mesh built 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 . The main model formulation features of IFS-FVM are summarized in Table 1 and shown alongside the corresponding IFS-ST properties discussed below in Sect. 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 asTS1 which describe the conservation laws of dry mass (Eq. 1a), momentum (Eq. 1b), dry entropy (Eq. 1c), and water substance (Eq. 1d). A summary of variables and physical constants is provided in Tables D1 and D2, respectively. Dependent variables in Eqs. (1a)-(1e) are dry density ρ d , threedimensional physical velocity vector u = [u, v, w] T , 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 the density of the individual water substance category ρ k to the density of dry air ρ d ); with the current cloud parameterization of the IFS (Forbes et al., 2010), five categories for water substance are considered (vapour r v , liquid r l , rain r r , ice r i , snow r s ), each described by the respective PDE (Eq. 1d). An additional prognostic equation for cloud fraction a employed with the IFS cloud parameterization is implemented as Furthermore, the thermodynamic variables are related by the gas law (Eq. 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. In Sect. 2.1.2, the fully compressible system (1a)-(1e)TS2 is augmented with the prognostic Helmholtz Eq. (17) derived from the advective form of the gas law (Eq. 1e) for the im-plicit (rather than explicit) solution with respect to the Exner pressure perturbation ϕ .
A quantity which appears in various right-hand-side (RHS) terms of the momentum Eq. (1b) is the density potential temperature θ ρ = θ ϑ, where ϑ ≡ (1 + r v /ε)/(1 + r t ) (Emanuel, 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 Eq. (1b) in order to expose the potential temperature perturbation θ for implicit coupling to the thermodynamic Eq. (1c) (see Sect. 2.1.2). Note that a high-order approximation of the buoyancy term applied in Kurowski et al. (2014) and 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 moist-precipitating dynamics at planetary scales (Sect. 3).
Another important aspect of the governing Eqs. (1a)-(1e) 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 Eqs. (1a)-(1e). A straightforward example applied in Sect. 3 of this paper is a stationary atmosphere u a (x), θ a (x), π a (x), r va (x) in thermal wind balance, which in terms of the momentum Eq. (1b) is where ϕ a ≡ c pd π a and θ ρa = θ a (1 + r va /ε)/(1 + r va ). The perturbation form (1a)-(1e) TS3 is analytically equivalent to the fully compressible equations for full variables but has favourable properties for numerical integration; see Sect. 2.1.2 and Prusa et al. (2008) and 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 . The optimal specification of ambient states for NWP and climate modelling is ongoing research. In general, ambient states can be as simple as stationary vertical profiles in hydrostatic balance, but bespoke definitions designed for a particular class of applications can substantially benefit the accuracy and robustness of the model integration.
The governing equations in Eqs. (1a)-(2) TS4 are formulated with respect to generalized curvilinear coordinates embedded in a geospherical framework. At the most elementary level, the generalized 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] T , the 3-D 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) and Kühnlein et al. (2012) for extended discussion. Following Smolarkiewicz et al. (2017), further symbols on the RHS of the momentum Eq. (1b) denote the gravity vector g = [0, 0, −g] T , the Coriolis parameter f ≡ 2 with the angular velocity vector of the Earth's rotation, and M subsumes the metric forces due to the curvature of the sphere: Details about the specification of the curvilinear space, as well as explicit expressions of the Coriolis term −f × u, the metric forces M, and the gravity g under shallow-or deepatmosphere equations, are given in Appendix C. Last but not least, the symbols P u = [P u , P v , P w ] T , P θ , P r k on the RHS of Eqs. (1a)-(1e) denote the respective forcings from physics parameterizations.
Note that additional RHS terms not explicitly provided in the governing Eqs. (1a)-(1e) may describe Rayleightype damping and/or Laplacian diffusion applied especially to model wave-absorbing layers at the domain boundaries; see e.g. Prusa and Smolarkiewicz (2003) and 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 (1a)-(1e)TS5 is accommodated in a generalized conservation law form: in which denotes the prognostic model variable, V the advector, G a generalized density, P represents the forcing from physics parameterization, and R the remaining righthand side 4 ; see Table 2 for the respective specifications of , V , G. The homogeneous mass continuity Eq. (1a) is a particular case of Eq. (6) and plays a fundamental role for the conservative advective transport of all other scalar variables. Note that because of the mass continuity Eq. (1a), the other scalar conservation laws (Eq. 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 Eq. (6) can be 4 As an example, R θ ≡ − G T u · ∇θ a when considering Eq. (1c).

6
C. Kühnlein et al.: FVM 1.0: a nonhydrostatic finite-volume dynamical core for the IFS subsumed in the following template scheme: and i = A i ( , V n+1/2 , G n , G n+1 , δt) in Eq. (7) symbolizes a flux-form Eulerian NFT advective transport scheme based on MPDATA, as described in  and Appendix A of the present paper. The n, n + 1/2, n + 1 indices denote full and half-time levels, and δt = t n+1 − t n is the time step increment of the semi-implicit dynamics. The vector index i = (k, i) marks the node positions k and i of the vertical and horizontal computational mesh, respectively, thereby revealing the 3-D co-located spatial arrangement of all dependent variables underlying IFS-FVM's discretization. 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 physics parameterization P , the semiimplicit 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 Eq. (6) (Smolarkiewicz and Margolin, 1993). Due to the congruency of Eq. (7) with the ODE (Eq. 9), the advection operator i may equally represent a second-order accurate semi-Lagrangian scheme . In terms of the coupling to physics parameterizations formally incorporated with first-order 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 Eq. (7); see Sect. 2.3 about physics-dynamics coupling.
There are two alternative implementations of the NFT advective transport scheme A i ( , V n+1/2 , G n , G n+1 , δt) based on MPDATA. First, the standard MPDATA formulations for integrating the fully compressible equations in the horizontally unstructured vertically structured discretization framework of IFS-FVM are provided in . These schemes 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), and Smolarkiewicz et al. (2013). Secondly, a more efficient horizontal-vertical split advective transport scheme based on MPDATA has been developed and is outlined in Appendix A. All IFS-FVM results presented in Sect. 3 were obtained with this horizontal-vertical split scheme. Note that the basic unsplit and horizontal-vertical split MPDATA schemes are second-order accurate given the advector V n+1/2 at the intermediate time level t n+1/2 provided as an O(δt 2 ) estimate, which is explained below.
Given the preceding discussion, the semi-implicit solution procedure of the governing system (1a)-(1e) TS6 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 Eq. (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 Eq. (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) defines the advector V n+1/2 ≡ (vGρ d ) n+1/2 as face-normal mass fluxes to the dual cell (vGρ d ) ⊥ | n+1/2 for the advective transport of all other scalar variables (see Table 2), a common approach to enable mass-compatible and monotonic solutions; see Smolarkiewicz et al. (2016) and  for a discussion in the context of IFS-FVM.
Given the tendencies from physics parameterization P θ , P u , P v , P r k , P a formally evaluated at t n and the advective transport of all scalar variables θ ≡ ai -with the advected water content mixing ratios r k i and cloud fraction ai already representing the final solutions at t n+1 -the scheme (7)-(8) for the thermodynamic Eq. (1c) and momentum Eq. (1b) is implemented as where .
Due to the presence of non-linear terms, the discrete system (11a)-(11b) TS7 is executed iteratively (Smolarkiewicz and Dörnbrack, 2008;Smolarkiewicz et al., 2014), with lagged quantities from the previous iteration denoted by the superscript . Typically, it uses one corrective iteration, which results in a predictor-corrector approach. The predictor of Eqs. (11a)-(11b) alone is second-order accurate, and the predictor-corrector already closely approximates the corresponding trapezoidal integral (Smolarkiewicz and Szmelter, 2009; 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 generalized density G, and a , b the weights for the incorporation of the RHS forcings R . The rightmost column refers to the governing equation for each dependent variable . Smolarkiewicz et al., 2014Smolarkiewicz et al., , 2019. With the respective prognostic model variables θ and u in Eqs. (11a)-(11b) 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 Eq. (12) that is composed of the already completed r n+1 k . At the beginning of the iterative execution of Eqs. (11a)-(11b), 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 Eqs. (11a)-(11b), closed-form expressions for the discrete velocity update are derived by eliminating Eq. (11a) for θ i in the buoyancy term of Eq. (11b) -thereby implementing 3-D fully implicit treatment of buoyant modes in a moist-precipitating atmosphere -and gathering the terms with linear dependence on u i on the left-hand side (LHS), which results in As all terms are co-located, the spatial mesh vector index i has been omitted in Eq. (14). The RHS of Eq. (14) is composed of all explicitly known terms, summarized 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, Eq. (14) can be symbolized as and respectively, 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 Eq. (16) still requires the Exner pressure perturbation ϕ to be specified. A straightforward computation may employ the gas law (Eq. 1e) using the updated variables ρ d , θ , and r v . However, the resulting 3-D explicit acoustic integration is subject to very small time steps (in order to maintain numerical stability) and thus inefficient for NWP. Therefore, a final step in IFS-FVM's numerical solution procedure is to augment the fully compressible Eqs. (1a)-(1e) with an auxiliary 3-D implicit boundary value problem for the pressure perturbation variable ϕ . The formulation of this implicit boundary value problem originates from the advective (or Lagrangian) form of the gas law (Eq. 1e) combined with the respective advective forms of the mass continuity Eq. (1a), thermodynamic Eq. (1c), and water vapour Eq. (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 Eqs. (1a)-(1e) and reads where ϕ = ϕ a + ϕ and We interpret Eq. (17) in terms of the generalized transport Eq. (6) for ≡ ϕ (see Table 2), while setting and R ϕ as the remaining RHS consisting of the three divergence operators. Importantly, the P ϕ given by Eq. (19) describes the pressure adjustment from the physics parameterization tendencies P θ and P r v . We have implemented the integration of Eq. (17) using the template semi-implicit scheme (7) with parameters a ϕ = (1−α) and b ϕ = α, where α ∈ [0.5, 1.0]. In the limit α = 1.0, the template Eq. (7) represents the first-order backward Euler scheme with regard to R ϕ . For α = 0.5, the template Eq. (7) becomes the secondorder trapezoidal scheme; in practice we may use weak offcentring (e.g. α = 0.51) for regularization. For any specification of α, coupling with the solution procedure of the fully compressible Eqs. (1a)-(1e) is implemented through Eq. (16) that enters into the three occurrences of u on the RHS of Eq. (17) in R ϕ | n+1 . Furthermore, as indicated in Eq. (7), the (explicit) forward Euler scheme is used with regard to the forcing P ϕ . Reorganizing terms finally yields the elliptic Helmholtz equation for the pressure perturbation variable ϕ at the future time level t n+1 , which can be written compactly as where the spatial grid index i has been omitted again. The summation in Eq. (20) is over the three divergence operators on the RHS of Eq. (17). The symbolic notations L(ϕ ) and R refer to the implicit and explicit parts of the equation, respectively. The coefficients A , ζ , and B are defined accordingly as Smolarkiewicz et al. (2016Smolarkiewicz et al. ( , 2017 and  that all used the backward Euler scheme α = 1. The 3-D elliptic boundary value problem (Eq. 20) is solved using a nonsymmetric preconditioned generalized conjugate residual (GCR) approach; see Smolarkiewicz and Szmelter (2011) for a recent discussion and Smolarkiewicz and Margolin (2000) and Smolarkiewicz et al. (2004) for tutorials. The crux of the preconditioning is the direct inversion of the vertical part of the problem that dramatically reduces the condition number of the full linear operator L and enables rapid convergence of the GCR solver. The preconditioner P ≈ L discards the off-diagonal entries of the matrix G T C. Inversion of P then utilizes a weighted line Jacobi method explained in Appendix B. In the GCR solver and the preconditioner P, the coefficients A 2 , A 3 , and B in Eq. (21) depend on ϕ = ϕ a +ϕ , lagged behind in the outer iteration of Eq. (11b). Importantly, in the predictor step of the outer iteration, an improved first guess for ϕ , can be achieved by employing the explicit forward Euler scheme for Eq. (17), obtained by setting α = 0 in Eq. (7).
As far as the adiabatic dynamics is concerned, the O(δt 2 ) integration of Eq. (17) with the first-order backward Euler scheme α = 1 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. The accumulation of first-order errors in ϕ can be mitigated by solving Eq. (17) with the weakly off-centred trapezoidal scheme using α = 0.51; see Benacchio et al. (2014) for alternative design and pertinent discussion. All IFS-FVM results presented in Sect. 3 were obtained with the backward Euler scheme α = 1, as this has been the default in the earlier implementations. However, using the extended-range forecast of the baroclinic instability benchmark, we will demonstrate that either α = 1 or α = 0.51 shows the same close agreement with IFS-ST. Moreover, both choices α = 1 or α = 0.51 provide essentially identical computational performance.
Overall, the 3-D implicit scheme with respect to ϕ permits time steps equivalent to soundproof models (Kurowski et al., 2014;Smolarkiewicz et al., 2014;. Together with the full 3-D implicit incorporation of buoyant and rotational modes described above, the semi-implicit integration is unconditionally stable with respect to all waves supported by the fully compressible Eqs. (1a)-(1e), and thus the semi-implicit model time step δt can be selected according to the stability of the advective transport scheme A i ( , V n+1/2 , G n , G n+1 , δt) in Eq. (7).

Spatial discretization
The discretization framework of IFS-FVM combines a structured-grid FD-FV method in the vertical with an unstructured 5 FV approach in the horizontal . The FV discretization and differentiation on the spherical surfaces adopt the median-dual approach described in Szmelter and Smolarkiewicz (2010). All dependent variables are co-located in the nodes in 3-D. The consistent spa- Figure 1. Schematic of the median-dual mesh in 2-D. The edge connecting nodes i and j of the primary polygonal mesh pierces, precisely in the edge centre, the face S j 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 blue lines indicate dual cells with control volumes V i and V j , respectively. tial discretization of the applied MPDATA schemes in IFS-FVM, symbolized by the operator A i in Eq. (7), is described in .
The schematic in Fig. 1 illustrates an arbitrary unstructured mesh on a 2-D 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 V i , cell face area S j , and face normals n j are evaluated from vector calculus in the computational space, i.e. in terms of x and y coordinates (see Sect. 2.1.1) on a zonally periodic horizontal plane (Szmelter and Smolarkiewicz, 2010). The unstructured FV discretization in the horizontal is combined with the standard second-order FD-FV method on the vertical structured grid with an independent coordinate z.
For a differentiable vector field A, the Gauss divergence theorem ∇ · A = ∂ A · n applied over the control volume surrounding node i = (k, i) in 3-D computational space is given as where l(i) numbers the edges connecting node (k, i) with its horizontal neighbours (k, j ), and S j refers both to the face per se and its surface area 6 . The geometric quantities S j , V i , and δz of the computational space 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. Elementary approximations are 6 Note that in IFS-FVM, S j and V i have dimensions of length and area, and the actual face areas and volumes of prismatic cells are S j δz and V i δz, respectively, in computational space  given as A ⊥ k,j = 0.5 n j · (A k,i + A k,j ) in the horizontal and A z k+1/2,i = 0.5(A z k+1,i + A z k,i ) in the vertical (Szmelter and Smolarkiewicz, 2010;Smolarkiewicz et al., 2016). However, when applied in Eq. (22) without non-linear operations involved at the faces, cancellation of the node value A k,i occurs, and this is exploited here to obtain a more efficient implementation of the median-dual approach by simply using A ⊥ k,j = 0.5 n j ·A k,j and A z k+1/2,i = 0.5A z k+1,i . Similarly, applying the same cancellation of the node value k,i of a differentiable scalar field , the 3-D nabla operator ∇ in computational space interpreted in terms of the Gauss divergence theorem is 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 . Given Eq. (22) and Eq. (23) in the computational space, they are augmented with metrics of the curvilinear coordinate framework to obtain the respective physical divergence and gradient appearing in the governing equations of Sect. 2.1.1 and 2.1.2, respectively. For illustration, one example for the revised implementation of the velocity divergence in the first term on the RHS of Eq. (17) at the node (k, i) is The details about the specification of the curvilinear coordinate framework under shallow-and deep-atmosphere equations are described in Appendix C. For IFS-FVM, the mesh generation and mesh data structures, as well as the nearest-neighbour distributedmemory communication using MPI, are handled by ECMWF's Atlas library, comprehensively described in Deconinck et al. (2017). Atlas is also designed to make use of specific programming paradigms to support accelerators, although these have not yet been explored with IFS-FVM. For the quasi-uniform octahedral reduced Gaussian grid (Sect. 2.4), the parallelization of Atlas adopts the equalregion horizontal domain decomposition of IFS. The nearestneighbour communications enabling the FV stencil operations in the parallel horizontal domain are performed on overlap regions (halos) between partitions, as well as a "peri-10 C. Kühnlein et al.: FVM 1.0: a nonhydrostatic finite-volume dynamical core for the IFS odic overlap" for the east-west boundary on the spherical-spheroidalCE1 Earth. With the median-dual FV approach presented above, an overlap region of one element is typically used.
In the present work, the programming of the discrete differential operators (22) and (23) in the IFS-FVM modern Fortran code has been comprehensively revised 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 a purely node-based implementation. When evaluating Eqs. (22) and (23), 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 performs 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 discretization leading to stencil computations, IFS-FVM is naturally cache and memory-bandwidth bound. The node-based programming improves the flop-perbyte ratio and relaxes the dependency on relatively (compared to pure computation) slow memory access. As a consequence, it provides a significant gain in efficiency 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 shared-memory parallelization 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.
We note that the GCR solver for the Helmholtz problem (20) uses global communications in the computation of inner products (characteristic of Krylov-subspace methods) and residual error norms. This has been implemented for the present paper such that summation is first done locally for each distributed-memory task, and the resulting single numbers are then reduced over all tasks.

Spectral-transform IFS
The spectral-transform IFS (denoted as IFS-ST in this paper) that is operational at ECMWF is based on the hydrostatic primitive equations (HPEs). The HPEs are formulated in a hybrid sigma-pressure terrain-following vertical coordinate following Simmons and Burridge (1981). The HPEs are integrated numerically using the efficient centred two-time-level SISL scheme, which permits long time steps due to the unconditional stability provided by the fully implicit treatment of the fast acoustic 7 and buoyant modes, and the 3-D SL advection (Ritchie et al., 1995;Temperton et al., 2001;Hortal, 2002). Therefore, the constant time step of the SISL integration in IFS-ST can be selected according to optimal efficacy rather than stability. The ST method, which is applied along model levels in the horizontal, is combined with a finiteelement (FE) approach to discretize the integral operator in the vertical direction (Untch and Hortal, 2004) 8 . In 2002, this vertical FE scheme of Untch and Hortal (2004) with a co-located arrangement of prognostic variables replaced the former FD scheme of Simmons and Burridge (1981)  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 gridpoint, 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 non-linearities in the RHS forcing terms in the governing 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 for which it always approaches 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 the efficiency and accuracy of the IFS-ST at ECMWF . 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, which is also employed with IFS-FVM, is reviewed in Sect. 2.4.
The semi-Lagrangian advection scheme in IFS-ST is based on Ritchie et al. (1995). The SL trajectories are com-puted 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 the arrival, middle, and departure points, respectively, and n and n+1 again refer to the current and future time step. In addition, the operator L symbolizes the linear operator and N the non-linear residual N = M−L, with M being the full non-linear model. The operator L results from the linearization of M with regard to a horizontally homogeneous reference state. Before the final semi-implicit solution with Eq. (25) 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 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 3-D 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 (ICI) scheme (Bénard et al., 2010) -for stability in global configurations. The ICI scheme in IFS-ST requires recomputation of the semi-Lagrangian trajectory and interpolations in the corrector step, which is not needed in the predictor-corrector approach of IFS-FVM (Sect. 2.1.2). The prognostic variables and main characteristics of the NH formulation are also given in Table 1. Furthermore, the NH formulation is currently restricted to the vertical FD scheme by Bubnová et al. (1995). By default, the HPE and NH formulations of IFS-ST use the shallow-atmosphere approximation. In Sect. 3, we will com-pare the computational performance of the HPE and NH formulations of IFS-ST against IFS-FVM.

Some aspects of physics-dynamics coupling
The IFS physics parameterization package at ECMWF is applied in the same configuration throughout the mediumrange, sub-seasonal, and seasonal forecasting systems. The physics package includes parameterization of radiation, moist convection, clouds and stratiform precipitation, surface processes, sub-grid-scale turbulence, and orographic and non-orographic gravity wave drag.
In IFS-ST, the physics-dynamics coupling employs the SLAVEPP (semi-Lagrangian averaging of physical parameterization) scheme (Wedi, 1999). SLAVEPP targets secondorder accuracy by averaging tendencies from selected physics parameterizations along the SL trajectory at the midpoint in space-time between the departure point at t n and arrival point 9 at t n+1 CE2 . Thereby, 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 the IFS uses sequential splitting of tendencies, 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. More details about the IFS physics parameterizations can be found in the general IFS documentation.
The physics-dynamics coupling in IFS-FVM differs from IFS-ST. As explained in Sect. 2.1.2, the current implementation of IFS-FVM incorporates the tendencies from physics parameterizations P by means of a first-order coupling at t n . Therefore, the fields that enter the parameterizations 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 parameterization deviates from IFS-ST, the sequential splitting between the various IFS physics parameterizations is kept exactly the same. Incorporating the physics parameterizations 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 physics parameterizations will be explored in the future.
The IFS-FVM code has its own interface to the IFS physics parameterizations. Among others, it involves conversion between IFS-FVM's variables and those employed 9 SLAVEPP is applied to tendencies from radiation, moist convection, and the cloud scheme, whereas tendencies from turbulence and gravity wave drag parameterizations are incorporated with first order at t n+1 (Wedi, 1999).  Table 1) 10 , interpolations to vertical interfaces, and the 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 the IFS physics parameterization package becomes more harmonized with IFS-FVM. However, generally the coupling is facilitated by common features of IFS-FVM and IFS-ST, such as longitudelatitude coordinates, the octahedral reduced Gaussian grid, and the co-located arrangement of dependent variables. The IFS-FVM interface to the physics parameterizations 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 : 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. The IFS-FVM has been coded rigorously with a variable time stepping capability . In the case of the physics-dynamics coupling with subcycling Eq. (26), 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 radiation 11 , the current operational IFS-ST evaluates the physics parameterizations 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 12 , 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 parameterization package can be up to 40 % of the forecast model. For the idealized DCMIP experiments considered in the present work, we focus on the parameterization of clouds and stratiform precipitation incorporated by means of the general IFS-FVM and IFS-ST interfaces described above.
10 Examples are conversions between mixing ratios r k and specific water content variables q k or between quantities in the heightversus pressure-based coordinate systems. 11 In the current high-resolution deterministic IFS forecasts on the O1280 grid, the radiation scheme is called every hour, compared to the semi-implicit model time step δt = 450 s, and is also run on the coarser O400 grid. 12 With the current formulations, the time step δt in IFS-FVM is typically about a factor of 6-7 smaller than in IFS-ST.

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") 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 a minimum of ∼ 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. Figure 2 depicts the O24 octahedral grid nodes, together with the corresponding edges of the primary mesh as applied in IFS-FVM. Also shown is the dual mesh spacing of the O1280 grid. 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 Sect. 3.

Benchmark simulation results
We study the solution quality and also the computational efficiency of IFS-FVM in comparison to the established IFS-ST. The hydrostatic IFS-ST represents a proven formulation for global medium-range NWP at ECMWF, and the aim is to reproduce its results with the novel IFS-FVM.
Baroclinic instability represents a common and relevant test problem to evaluate the performance of global NWP models in the large-scale hydrostatic regime. The underlying processes are fundamental to the life cycle (i.e. formation to decay) of high-and low-pressure systems in the mid-latitude "storm tracks" of the Earth's atmosphere. Here, we adopt the experimental set-up for a baroclinic wave life cycle used in the 2016 edition of DCMIP following Ullrich et al. (2014) and the documentation available in Ullrich et al. (2016). Note that the IFS is a highly optimized single-application model for real-time numerical weather prediction, and adapting the code to idealized configurations is not straightforward. Because of this and also a particular interest in requirements and applications at ECMWF, we consider configurations and physics parameterizations that depart from the test case specifications of DCMIP in Ullrich et al. (2016).
To study the accuracy of the novel nonhydrostatic IFS-FVM based on the finite-volume discretization, we verify its solution quality against IFS-ST based on the spectraltransform approach. This comparison is performed in Sect. 3.1 for dry adiabatic simulations of the baroclinic instability, i.e. dynamical core only, and then in Sect. 3.2 un- Figure 2. The locations of the octahedral reduced Gaussian grid nodes are shown in (a), using for illustration the very coarse O24 grid example with 24 latitudes between the pole and equator. Panel (b) depicts the associated edges of the primary mesh connecting the nodes as applied in the context of the FV discretization of IFS-FVM. Panel (c) provides the local spacing of the FV dual mesh for the O1280 grid corresponding to the highest-resolution deterministic forecast model currently at ECMWF. der consideration of moist-precipitating processes that involve coupling to the prognostic single-moment bulk microphysics parameterization of the operational IFS at ECMWF (Forbes et al., 2010). The computational efficiency of IFS-FVM alongside the hydrostatic and nonhydrostatic IFS-ST is studied in Sect. 3.3. Note that we use the nonhydrostatic IFS-ST only when looking at the efficiency of the various dynamical core formulations. In terms of the solution quality for the baroclinic instability, nonhydrostatic effects are entirely negligible at the considered coarse resolutions (see next paragraph), and therefore only results with the hydrostatic IFS-ST are shown and analysed.
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 Sect. 2.2 and 2.4 for the nomenclature that defines the grids. For the two horizontal grid sizes, both IFS-FVM and IFS-ST employ 60 stretched vertical levels. The height coordinate of IFS-FVM (Appendix C) is specified exactly according to the computed height of the hybrid sigma-pressure coordinate of IFS-ST at the initial time of the simulation. The lowest full level is located at a height of about 10 m, and the vertical spacing between model levels ranges from about 12 m near the ground to 4 km near the model top located at 48 km. In contrast to the height coordinate levels in IFS-FVM, the hybrid sigma-pressure coordinate levels in IFS-ST change with time, but overall the vertical spacing remains quite similar in both models over the course of the 15-day simulations. The similar spacing is particularly true for the vertical levels near the surface, where the terrain-following character of the hybrid sigmapressure coordinate dominates. Note that, as the maximum local difference in height of the lowest full levels is found to be smaller than 1 m over the 15-day baroclinic instability simulation, output fields of IFS-FVM and IFS-ST may be straightforwardly compared at this level without using inter-polation 13 . For the comparison of computational efficiency in Sect. 3.3, we employ the (O1280,TCo1279) horizontal grid, corresponding to about 9 km spacing, with 137 stretched vertical levels, again with a similar distribution in IFS-FVM and IFS-ST. 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. In contrast, IFS-FVM generally applies variable time stepping that targets a maximum horizontal advective Courant number of 0.95 -the actual maximum 3-D Courant number can be significantly larger than 1 as this is permitted by the horizontalvertical split NFT advection in IFS-FVM (Appendix A). As explained in Sect. 2.1.2, we have implemented the integration of Eq. (17) using an off-centred variant of the template semi-implicit scheme (7). All IFS-FVM results presented in this paper used the backward Euler scheme α = 1, as this has been the default so far in previous implementations. However, one exception to this is the middle panel in Fig. 5, which shows, for comparison to the default α = 1, the corresponding result obtained with the weakly off-centred trapezoidal scheme using α = 0.51.
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 midlatitudes of each global hemisphere that are in thermal wind balance with the meridional temperature gradient. The definition of the balanced initial state is given by analytical functions provided in Ullrich et al. (2014). Here, this balanced zonal flow is also used to specify the ambient state variables u a (y, z), θ a (y, z), π a (y, z) of the governing fully compressible Eqs. (1a)-(1e) that fulfil Eq. (4). 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), in which 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 the study of the solution symmetry about the equator, and is more relevant to real weather. Nevertheless, for reference we provide one illustration in Fig. 6 for the set-up in which the triggering of the baroclinic instability is applied in the Northern Hemisphere only. After an initial period of linear growth, the instability enters the non-linear 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 horizontal cross section of nearsurface pressure and the zonal height cross section of meridional wind, respectively, for the dry adiabatic simulations at day 10. At this stage, a large-amplitude baroclinic wave has developed and formed sharp fronts in the lower troposphere. Generally very close agreement is found between the finite-volume (IFS-FVM) and spectral-transform (IFS-ST) solutions. This is emphasized by the difference plots in the bottom row of the horizontal and vertical cross sec-tions 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 non-linear baroclinic instability evolution at day 15 for the finer of the two grids. This depiction shows close agreement between IFS-FVM (using the default α = 1 as well as the α = 0.51) and IFS-ST also at this later stage. The various dynamical core formulations provide visibly symmetric solutions around the equator, and there are no signs of any significant grid imprinting at the scale of the wavenumberfour irregularities of the octahedral grid. The latter is corroborated by Fig. 6 that shows the analogous result to Fig. 5 (again only α = 1 for IFS-FVM), but in contrast to all other plots, the trigger of the baroclinic instability is applied in the Northern Hemisphere only.
Kinetic energy spectra evaluated at day 15 in Fig. 7 reveal a strikingly similar distribution of variance across wavenumbers from 1 to ∼ 200. As all the simulations are without any explicit diffusion, the IFS-FVM spectra at the high wavenumbers attest to the implicit scale-selective regularization with artificial viscosity provided by the non-oscillatory finite-volume MPDATA advection 14 . 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 scale of the triangular-cubic truncation, corresponding to 4 times the grid spacing. This increased energy at these scales does not grow and is to a large extent controlled by the spectral filtering of the non-linear terms at every time step, among other mechanisms of implicit dissipation such as the semi-Lagrangian interpolation.

Simulation results for moist-precipitating configuration with the IFS cloud parameterization
Next we present results for the moist-precipitating baroclinic instability with coupling to the IFS cloud parameterization. Figure 8 shows the instantaneous large-scale precipitation rate at the surface 15 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 −1 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 9 is analogous to Fig. 8 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. 9. Here, the locations of the easternmost frontal zone and associated rainband agree closely considering the late stage of the baroclinic instability evolution. Figure  ysis of the simulations is supplemented in Fig. 11 with the time series of the minimum near-surface pressure (panel a) and the area-integrated precipitation rate (panel b). The temporal evolution of these two quantities is close between IFS-FVM and IFS-ST. Particularly, the minimum near-surface pressure agrees almost exactly at day 15, although small differences occur over the course of the simulation. The onset and subsequent increase in precipitation matches well in IFS-FVM and IFS-ST, and the later variations in the precipitation rate are similar, with no systematic underestimation or overestimation. Kinetic energy spectra evaluated at day 15 are shown in Fig. 12. Compared to the spectra of the dry simulations in Fig. 7, IFS-FVM and IFS-ST consistently show a considerably larger kinetic energy in the scales smaller than wavenumber ≈ 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 moistprecipitating FV formulations along with the coupling to the IFS physics parameterization.

Computational efficiency
The computational 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. The IFS-  FVM is envisaged for future applications in the nonhydrostatic regime running on future HPC architectures, but its computational performance on the current HPC facility at ECMWF sheds some light on its potential; see also Kühnlein and Smolarkiewicz (2019). Of interest is the relative performance of IFS-FVM to both the hydrostatic IFS-ST and its nonhydrostatic extension (see Sect. 2.2). In order to emphasize elementary aspects of the dynamical cores, here we as-sess the efficiency of the dry formulations only; i.e. no tracers or moisture variables, no physics parameterizations or coupled models, no I/O, and minimal diagnostics. Furthermore, we use the baroclinic instability benchmark of Sect. 3.1 but configured similar to HRES-ECMWF's highest-resolution deterministic forecast model. Figure 13 highlights the runtimes of the three different dynamical core formulations. The HRES forecast model configuration at ECMWF is currently based on the O1280/TCo1279 horizontal grid (corresponding to about 9 km grid spacing; Fig. 2c) with 137 vertical levels and run using 350 nodes (about one electrical group) on ECMWF's Cray XC40 supercomputer 16 . Importantly, while IFS-ST used the constant time step of 450 s, IFS-FVM employed variable time stepping according to the maximum permitted advective Courant number; therefore, in order to obtain realistic numbers for IFS-FVM, the timings were evaluated between day 10 and 15 in the fully non-linear stage of the baroclinic instability evolution. 16 Each node on this supercomputer consists of two Intel Xeon EP E5-2695 v4 "Broadwell" processors, each with 18 cores, which for the 350 compute nodes employed results in a total of 12 600 cores. Here, a hybrid MPI-OpenMP parallelization with six threads was used by all three dynamical cores. Figure 13 reveals that the time to solution that can be achieved with IFS-FVM for this configuration is only about twice as large as the operational hydrostatic IFS-ST and compares favourably with the nonhydrostatic IFS-ST. Although these performance measures merely represent snapshots at the current state of development, they highlight the potential of the numerical integration schemes applied in IFS-FVM to become competitive with state-of-the-art operational global weather forecasting models. Important aspects are that the FV method offers the prospect of better scalability and efficiency with respect to future HPC and that IFS-FVM employs substantially smaller time steps (again, about a factor 6-7 smaller compared to IFS-ST), which can be beneficial for accuracy. Further significant efficiency improvements of the IFS-FVM dynamical core are in preparation, and the work will be extended to physics-dynamics coupling with the smaller time steps.

Conclusions
Supporting substantially higher resolution in global NWP may ultimately demand local numerical discretizations to solve the governing nonhydrostatic equations in NWP models in a computationally efficient manner. The IFS-FVM successfully implements such a discretization and thus complements the operational hydrostatic IFS-ST and its nonhydrostatic extension at ECMWF. At the same time, the IFS-FVM introduces several useful new features into the IFS, such as conservative and monotone advective transport, deepatmosphere all-scale governing equations, and fully flexible unstructured FV meshes with optional variable resolution or meshes defined about the nodes of the operational octahedral grid.
The paper highlighted the semi-implicit NFT finitevolume integration of the fully compressible equations of the novel IFS-FVM considering comprehensive moistprecipitating dynamics with coupling to the IFS cloud parameterization by means of a generic interface applicable for coupling to the full IFS physics parameterization package. Developments such as the new horizontal-vertical directionally split NFT advective transport scheme based on MP-DATA, variable time stepping, effective preconditioning of the Krylov-subspace solver for the elliptic Helmholtz problem arising in the semi-implicit scheme, and an efficient node-based implementation of the median-dual FV approach provide a basis for the overall efficacy of IFS-FVM and application in global NWP at ECMWF.
The IFS-ST is applied successfully for operational forecasting at ECMWF and is therefore considered an appropriate reference model. It was shown that the presented semi-implicit NFT finite-volume integration scheme on colocated meshes can achieve comparable solutions to the proven spectral-transform IFS-ST. Here, the study focused on medium-and extended-range simulation of the dry and  moist-precipitating baroclinic instability benchmark at various resolutions. While the baroclinic instability benchmark aims at global atmospheric dynamics in the hydrostatic regime, referenced supplementary studies with IFS-FVM emphasize non-orographic and orographic flows in the nonhydrostatic regime. In addition to solution quality, we have demonstrated highly competitive computational efficiency of the presented semi-implicit NFT finite-volume integration of IFS-FVM in comparison to the semi-implicit semi-Lagrangian integration of IFS-ST.
Common aspects of the finite-volume and spectraltransform model formulations are the octahedral reduced Gaussian grid, the co-location of variables, the geospherical framework, and the physics parameterizations. Sharing these properties facilitates the comparison of the different discretizations and physics-dynamics coupling. Moreover, it provides numerous benefits for the general IFS model infrastructure, data assimilation, and ensemble system. Ongoing work advances IFS-FVM to full-physics global mediumrange NWP at convection-resolving resolutions .
Code availability. Model codes developed at ECMWF are the intellectual property of ECMWF and its member states, and therefore the IFS code is not publicly available. Access to a reduced version of the IFS code may be obtained from ECMWF under an OpenIFS licence (see http://www.ecmwf.int/en/research/projects/openifs for further information, last access: 6 February 2019).
Data availability. The model output data can be downloaded from https://doi.org/10.5281/zenodo.1445597 . We consider the advection operator A i in the two-time-level semi-implicit integration scheme (7) to be directionally split in the horizontal and vertical directions. This splitting is motivated by the observation that NWP models typically have a larger restriction on the time step in the vertical than the horizontal direction. For example, in the current operational configuration of the IFS run at TCo1279/L137 (≈ 9 km horizontal grid spacing and 137 stretched vertical levels), the advective Courant numbers are up to a factor of 2 larger in the vertical than in the horizontal direction. The horizontal-vertical splitting also accommodates IFS-FVM's unstructured horizontal discretization, enabling broad classes of meshes over the surface of the Earth's sphere-spheroid CE3 , and the structured grid in the (stiff) vertical direction.
The proposed scheme implements mass-compatible second-order Strang splitting as explained in the following. The overall semi-implicit integration of the fully compressible Eqs. (1a)-(1e) proceeds exactly as explained in Sect. 2.1.2, but with the 3-D 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 Eq. (1a) as , (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 (Eq. 8) asTS9 In Eqs. (A1) and (A2), the implementation of the horizontal advection transport A xy follows the horizontal part of the unstructured-mesh FV MPDATA of . The vertical scheme A z is a corresponding 1-D 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 . 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 cost 17 . Overall, the horizontal-vertical 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 higherorder, 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 a 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 Eqs. (7)-(8).

Appendix B: Weighted line Jacobi preconditioner
The bespoke preconditioner solves for the solution error e of the pressure perturbation variable ϕ : wherer denotes the residual error of Eq. (20). The preconditioning operator P is then decomposed into vertical and horizontal parts (Smolarkiewicz and Margolin, 2000), and the residual problem is solved iteratively according to P z (e µ+1 ) + P h (e µ ) −r = 0, where µ numbers the iterations, of which there are typically two. The vertical part P z is inverted directly with a tridiagonal algorithm. The horizontal part P h is lagged behind, except for its diagonal entries. The actual implementation is given as P z (e µ+1 ) + P h (e µ ) + D(e µ+1 − e µ ) −r = 0 , where D is the diagonal coefficient of P h , specified as with B 11 and B 22 referring to the diagonal entries of G T C. Subsequently, Eq. (B3) is executed as TS10 e µ+1 = ω D − P z −1 De µ −P h (e µ )+r +(1 − ω) e µ (B5) with the weight ω = 0.7. 17 Compared to the unsplit scheme, the particular horizontalvertical splitting also does not incur any additional parallel communication in the context of the horizontal domain decomposition of IFS-FVM.
Remarks from the language copy-editor CE1 According to our house standards, we use an en dash rather than a hyphen to indicate pairings (or in this case a dual option). Should this perhaps be "spherical or spheroidal" or "spherical and/or spheroidal"? CE2 A footnote cannot come after a full stop. Please check and confirm this alternative. CE3 We don't use slashes in this way according to our house standards. Should it be "or" or "and/or" instead here?
Remarks from the typesetter TS1 Please note that the font size of Eq. (1d) has been decreased to make it fit into one line. Furthermore, please check page break. If Eq. (1b) should not be split over two columns, please let me know.

TS8
I have contacted our image processing department to see if the figure quality can be improved. Please note that if the current versions do not meet your expectations we would kindly ask you to provide vector graphics in which our house standards have been applied (addition of panel labels, deleted hyphens during the copy-editing). Thank you very much in advance.

TS9
Please note that the font size was increased but the equation needed to be split.

TS10
Please note that this change will have to be approved by the editor. Could you please provide a short statement regarding this change that will be forwarded by us to the editor? Thank you very much in advance.

TS11
Please note that we already used the command "tilde" instead of "widetilde" for these instances. "Widetilde" has only been used for G.