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 finite-volume 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 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 integration scheme of the spectral-transform IFS.

1 Introduction

Notwithstanding the achievements made over the last decades (Bauer et al., 2015), 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.

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, non-linear 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¹ 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 advances helped to sustain the performance of the ST method (for details see Wedi et al., 2013, 2015; Wedi, 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, 2015). 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 predictor–corrector 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 distributed-memory 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 (Smolarkiewicz et al., 2016, 2017; Kühnlein and Smolarkiewicz, 2017). 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 (Wedi et al., 2015), 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, 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 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 (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 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², 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 dynamical core with selected parameterizations 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 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 finite-volume 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.

2 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 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 is provided (see main text for description). Abbreviations are as follows: finite element (FE), finite difference (FD), spectral transform (ST), finite volume (FV), two time level (2-TL), semi-implicit (SI), iterative-centred implicit (ICI). A summary of variables is provided in Table D1.

2.1 Finite-Volume Module of the IFS

IFS-FVM solves the deep-atmosphere³, 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 (Smolarkiewicz et al., 2014, 2016). In Sect. 2.1.2, we extend the IFS-FVM semi-implicit 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) (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 (Sect. 2.1.2 and Appendix A). In addition, improved efficacy with the Eulerian NFT MPDATA advection is sought by rigorous implementation of variable time stepping. The unstructured horizontal spatial discretization 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 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 (Deconinck et al., 2017). 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.

2.1.1 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

$$\begin{array}{}\text{(1a)}& {\displaystyle }{\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\right)=\mathrm{0},\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\mathit{u}}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\mathit{u}\right)=\\ {\displaystyle }& {\displaystyle }\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}[-{\mathit{\theta }}_{\mathit{\rho }}\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}^{\prime }+\mathit{g}(\mathrm{1}-{\displaystyle \frac{\mathit{\vartheta }}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}\left({\mathit{\theta }}_{\mathrm{a}}+{\mathit{\theta }}^{\prime }\right))\\ \text{(1b)}& {\displaystyle }& {\displaystyle }-\mathit{f}\times (\mathit{u}-{\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}{\mathit{u}}_{\mathrm{a}})+{\mathcal{M}}^{\prime }+{\mathit{P}}^{\mathit{u}}],\text{(1c)}& {\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\theta }}^{\prime }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\theta }}^{\prime }\right)=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}\left[-{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\cdot \mathrm{\nabla }{\mathit{\theta }}_{\mathrm{a}}+P^{{\mathit{\theta }}^{\prime }}\right],\text{(1d)}& {\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{r}_k}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{r}_k\right)=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{P}^{r_k},r_k=r_{\mathrm{v}},r_{\mathrm{l}},r_{\mathrm{r}},r_{\mathrm{i}},r_{\mathrm{s}}\text{(1e)}& {\displaystyle }& {\displaystyle }{\mathit{\phi }}^{\prime }=c_{\mathrm{pd}}\left[{\left({\displaystyle \frac{R_{\mathrm{d}}}{p_{\mathrm{0}}}}{\mathit{\rho }}_{\mathrm{d}}\mathit{\theta }(\mathrm{1}+r_{\mathrm{v}}/\mathit{\epsilon })\right)}^{R_{\mathrm{d}}/c_{\mathrm{vd}}}-{\mathit{\pi }}_{\mathrm{a}}\right],\end{array}$$

which describe the conservation laws of dry mass (Eq. 1a), momentum (Eq. 1b), dry entropy (Eq. 1b), and water substance (Eq. 1b). A summary of variables and physical constants is provided in Tables D1 and D2, respectively. Dependent variables in Eqs. (1a)–(1b) are dry density ρ_d, three-dimensional physical velocity vector $\mathit{u}=[u,v,w{]}^{\mathrm{T}}$, potential temperature perturbation θ^′, and a modified Exner pressure perturbation ${\mathit{\phi }}^{\prime }\equiv c_{\mathrm{pd}}{\mathit{\pi }}^{\prime }$, as well as the water substance mixing ratios $r_k={\mathit{\rho }}_k/{\mathit{\rho }}_{\mathrm{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. 1b). An additional prognostic equation for cloud fraction Λ_a employed with the IFS cloud parameterization is implemented as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathrm{\Lambda }}_{\mathrm{a}}}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathrm{\Lambda }}_{\mathrm{a}}\right)=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{P}^{{\mathrm{\Lambda }}_{\mathrm{a}}}.\end{array}$$

Furthermore, the thermodynamic variables are related by the gas law (Eq. 1b); the Exner pressure is $\mathit{\pi }=(p/p_{\mathrm{0}}{)}^{R_{\mathrm{d}}/c_{\mathrm{pd}}}$ and the potential temperature is $\mathit{\theta }=T/\mathit{\pi }$, where p is the total pressure, p₀ ≡ 10⁵ Pa, and T is the (absolute) temperature. In Sect. 2.1.2, the fully compressible system (1a)–(1b) is augmented with the prognostic Helmholtz Eq. (17) derived from the advective form of the gas law (Eq. 1b) for the implicit (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 $\mathit{\vartheta }\equiv (\mathrm{1}+r_{\mathrm{v}}/\mathit{\epsilon })/(\mathrm{1}+r_t)$ (Emanuel, 1994) with $\mathit{\epsilon }=R_{\mathrm{d}}/R_{\mathrm{v}}$, and the total water mixing ratio r_t represents the sum over all the individual mixing ratios:

$$\begin{array}{}\text{(3)}& r_t=\sum _k{r}_k=r_{\mathrm{v}}+r_{\mathrm{l}}+r_{\mathrm{r}}+r_{\mathrm{i}}+r_{\mathrm{s}}.\end{array}$$

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. (1b) (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)–(1b) is the underlying perturbation form. The perturbations of potential temperature ${\mathit{\theta }}^{\prime }=\mathit{\theta }-{\mathit{\theta }}_{\mathrm{a}}$ and Exner pressure ${\mathit{\pi }}^{\prime }=\mathit{\pi }-{\mathit{\pi }}_{\mathrm{a}}$ correspond to deviations from an ambient state (denoted by the subscript “a”) that satisfies a general balanced subset of Eqs. (1a)–(1b). 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

$$\begin{array}{}\text{(4)}& \mathrm{0}=-{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}_{\mathrm{a}}-\mathit{g}-\mathit{f}\times {\mathit{u}}_{\mathrm{a}}+\mathcal{M}\left({\mathit{u}}_{\mathrm{a}}\right),\end{array}$$

where φ_a≡c_pdπ_a and ${\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}={\mathit{\theta }}_{\mathrm{a}}(\mathrm{1}+r_{\mathrm{va}}/\mathit{\epsilon })/(\mathrm{1}+r_{\mathrm{va}})$. The perturbation form (1a)–(1b) 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 (Smolarkiewicz et al., 2019). 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.

All governing equations 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 $\mathit{x}=[x,y,z{]}^{\mathrm{T}}$, the 3-D nabla operator ∇ with respect to x, the Jacobian 𝒢 of the coordinate transformations (i.e. the square root of the determinant of the metric tensor), a matrix of metric coefficients $\stackrel{\mathrm{̃}}{\mathbf{G}}$, its transpose ${\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}$, and the contravariant velocity $\mathit{v}=\dot{\mathit{x}}={\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}+{\mathit{v}}^g$ where ${\mathit{v}}^g\equiv \partial \mathit{x}/\partial 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 $\mathit{g}=[\mathrm{0},\mathrm{0},-g{]}^{\mathrm{T}}$, the Coriolis parameter f≡2Ω with Ω the angular velocity vector of the Earth's rotation, and 𝓜^′ subsumes the metric forces due to the curvature of the sphere:

$$\begin{array}{}\text{(5)}& {\mathcal{M}}^{\prime }(\mathit{u},{\mathit{u}}_{\mathrm{a}},{\mathit{\theta }}_{\mathit{\rho }}/{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}})=\mathcal{M}\left(\mathit{u}\right)-({\mathit{\theta }}_{\mathit{\rho }}/{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}})\mathcal{M}\left({\mathit{u}}_{\mathrm{a}}\right).\end{array}$$

Details about the specification of the curvilinear space, as well as explicit expressions of the Coriolis term $-\mathit{f}\times \mathit{u}$, the metric forces 𝓜, and the gravity g under shallow- or deep-atmosphere equations, are given in Appendix C. Last but not least, the symbols ${\mathit{P}}^{\mathit{u}}=[P^u,P^v,P^w{]}^{\mathrm{T}}$, $P^{{\mathit{\theta }}^{\prime }}$, $P^{r_k}$ on the RHS of Eqs. (1a)–(1b) denote the respective forcings from physics parameterizations.

Note that additional RHS terms not explicitly provided in the governing Eqs. (1a)–(1b) 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) and Klemp et al. (2008).

2.1.2 Semi-implicit numerical integration

To facilitate a compact description of the integration scheme, each of the governing Eqs. (1a)–(1b), (2), (17) is accommodated in a generalized conservation law form:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\partial G\mathrm{\Psi }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{V}\mathrm{\Psi }\right)=G\left({\mathcal{R}}^{\mathrm{\Psi }}+P^{\mathrm{\Psi }}\right),\end{array}$$

in which Ψ denotes the prognostic model variable, V the advector, G a generalized density, P^Ψ represents the forcing from physics parameterization, and ℛ^Ψ the remaining right-hand side⁴; 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 $\mathrm{d}\mathrm{\Psi }/\mathrm{d}t={\mathcal{R}}^{\mathrm{\Psi }}+P^{\mathrm{\Psi }}$, where $\mathrm{d}/\mathrm{d}t=\partial /\partial t+\mathit{v}\cdot \mathrm{\nabla }$ represents the total derivative.

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 ℛ^Ψ. The rightmost column refers to the governing equation for each dependent variable Ψ.

Building on the earlier works by Smolarkiewicz et al. (2014, 2016), the two-time-level numerical integrators of IFS-FVM for Eq. (6) can be subsumed in the following template scheme:

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Psi }}_{\mathit{i}}^{n+\mathrm{1}}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},{\mathit{V}}^{n+\mathrm{1}/\mathrm{2}},G^n,G^{n+\mathrm{1}},\mathit{\delta }t)+b^{\mathrm{\Psi }}\mathit{\delta }t{\mathcal{R}}^{\mathrm{\Psi }}{|}_{\mathit{i}}^{n+\mathrm{1}}\\ \text{(7)}& {\displaystyle }& {\displaystyle }\equiv {\widehat{\mathrm{\Psi }}}_{\mathit{i}}+b^{\mathrm{\Psi }}\mathit{\delta }t{\mathcal{R}}^{\mathrm{\Psi }}{|}_{\mathit{i}}^{n+\mathrm{1}},\end{array}$$

where

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{̃}}{\mathrm{\Psi }}={\mathrm{\Psi }}^n+a^{\mathrm{\Psi }}\mathit{\delta }t{\mathcal{R}}^{\mathrm{\Psi }}{|}^n+\mathit{\delta }t{P}^{\mathrm{\Psi }}{|}^n,\end{array}$$

and ${\widehat{\mathrm{\Psi }}}_{\mathit{i}}={\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},{\mathit{V}}^{n+\mathrm{1}/\mathrm{2}},G^n,G^{n+\mathrm{1}},\mathit{\delta }t)$ in Eq. (7) symbolizes a flux-form Eulerian NFT advective transport scheme based on MPDATA, as described in Kühnlein and Smolarkiewicz (2017) and Appendix A of the present paper. The n, $n+\mathrm{1}/\mathrm{2}$, n+1 indices denote full and half-time levels, and $\mathit{\delta }t=t^{n+\mathrm{1}}-t^n$ is the time step increment of the semi-implicit dynamics. The vector index $\mathit{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 ℛ^Ψ at tⁿ and tⁿ⁺¹, respectively, are given in Table 2. Apart from the incorporation of the physics parameterization P^Ψ, the semi-implicit scheme (7) with weights $a^{\mathrm{\Psi }}\equiv b^{\mathrm{\Psi }}\equiv \mathrm{0.5}$ is fully congruent with the second-order trapezoidal-rule trajectory integral of the corresponding ordinary differential equation

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\mathrm{\Psi }}{\mathrm{d}t}}={\mathcal{R}}^{\mathrm{\Psi }}\end{array}$$

of Eq. (6) (Smolarkiewicz and Margolin, 1993). Due to the congruency of Eq. (7) with the ODE (Eq. 9), the advection operator ${\widehat{\mathrm{\Psi }}}_{\mathit{i}}$ may equally represent a second-order accurate semi-Lagrangian scheme (Smolarkiewicz et al., 2014). In terms of the coupling to physics parameterizations formally incorporated with first-order accuracy, the associated forcing $P^{\mathrm{\Psi }}{|}^n=P^{\mathrm{\Psi }}(t_{\mathrm{phys}},\mathrm{\Delta }{t}_{\mathrm{phys}})$ can optionally be evaluated with an equal or longer time step Δt_phys=N_sδt (with integer $N_{\mathrm{s}}=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }$) 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 ${\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},{\mathit{V}}^{n+\mathrm{1}/\mathrm{2}},G^n,G^{n+\mathrm{1}},\mathit{\delta }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 Kühnlein and Smolarkiewicz (2017). 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 ${\mathit{V}}^{n+\mathrm{1}/\mathrm{2}}$ at the intermediate time level $t^{n+\mathrm{1}/\mathrm{2}}$ provided as an 𝒪(δt²) estimate, which is explained below.

Given the preceding discussion, the semi-implicit solution procedure proceeds from an atmospheric state at tⁿ to a state at tⁿ⁺¹ as described in the following. The solution procedure commences with the integration of the mass continuity Eq. (1a) as

$$\begin{array}{}\text{(10)}& {\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{n+\mathrm{1}}={\mathcal{A}}_{\mathit{i}}({\mathit{\rho }}_{\mathrm{d}}^n,(\mathit{v}\mathcal{G}{)}^{n+\mathrm{1}/\mathrm{2}},{\mathcal{G}}^n,{\mathcal{G}}^{n+\mathrm{1}},\mathit{\delta }t),\end{array}$$

which straightforwardly returns the updated density ${\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{n+\mathrm{1}}$. The 𝒪(δt²) estimate for the advector $(\mathit{v}\mathcal{G}{)}^{n+\mathrm{1}/\mathrm{2}}$ in Eq. (10) is implemented here by linear extrapolation of the advective velocities from the previous time levels tⁿ⁻¹ and tⁿ; 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 ${\mathit{V}}^{n+\mathrm{1}/\mathrm{2}}\equiv (\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{n+\mathrm{1}/\mathrm{2}}$ as face-normal mass fluxes to the dual cell $(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{⟂}{|}^{n+\mathrm{1}/\mathrm{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 Kühnlein and Smolarkiewicz (2017) for a discussion in the context of IFS-FVM.

Given the tendencies from physics parameterization $P^{{\mathit{\theta }}^{\prime }}$, P^u, P^v, $P^{r_k}$, $P^{{\mathrm{\Lambda }}_{\mathrm{a}}}$ formally evaluated at tⁿ and the advective transport of all scalar variables ${\widehat{{\mathit{\theta }}^{\prime }}}_{\mathit{i}}$, ${\widehat{\mathit{u}}}_{\mathit{i}}\equiv ({\widehat{u}}_{\mathit{i}},{\widehat{v}}_{\mathit{i}},{\widehat{w}}_{\mathit{i}})$, ${r_k}_{\mathit{i}}^{n+\mathrm{1}}\equiv {\widehat{r_k}}_{\mathit{i}}$, ${{\mathrm{\Lambda }}_{\mathrm{a}}}_{\mathit{i}}^{n+\mathrm{1}}\equiv {\widehat{{\mathrm{\Lambda }}_{\mathrm{a}}}}_{\mathit{i}}$ – with the advected water content mixing ratios ${\widehat{r_k}}_{\mathit{i}}$ and cloud fraction ${\widehat{{\mathrm{\Lambda }}_{\mathrm{a}}}}_{\mathit{i}}$ already representing the final solutions at tⁿ⁺¹ – the scheme (7)–(8) for the thermodynamic Eq. (1b) and momentum Eq. (1b) is implemented as

$$\begin{array}{ll}\text{(11a)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_{\mathit{i}}^{\prime }={\widehat{{\mathit{\theta }}^{\prime }}}_{\mathit{i}}-\mathrm{0.5}\mathit{\delta }t{\left[{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\cdot \mathrm{\nabla }{\mathit{\theta }}_{\mathrm{a}}\right]}_{\mathit{i}}{\displaystyle }& {\displaystyle }{\mathit{u}}_{\mathit{i}}={\widehat{\mathit{u}}}_{\mathit{i}}+\mathrm{0.5}\mathit{\delta }t[-{\mathit{\theta }}_{\mathit{\rho }}^{\star }\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}^{\prime }+\mathit{g}(\mathrm{1}-{\displaystyle \frac{{\mathit{\vartheta }}^{n+\mathrm{1}}}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}\left({\mathit{\theta }}_{\mathrm{a}}+{\mathit{\theta }}^{\prime }\right))\\ \text{(11b)}& {\displaystyle }& {\displaystyle }-\mathit{f}\times (\mathit{u}-{\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}^{\star }}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}{\mathit{u}}_{\mathrm{a}})+{\mathcal{M}}^{\prime }({\mathit{u}}^{\star },{\mathit{u}}_{\mathrm{a}},{\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}^{\star }}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}){]}_{\mathit{i}},\end{array}$$

where

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_{\mathit{\rho }}^{\star }={\mathit{\theta }}^{\star }{\mathit{\vartheta }}^{n+\mathrm{1}}\equiv {\displaystyle \frac{{\mathit{\theta }}^{\star }(\mathrm{1}+r_{\mathrm{v}}^{n+\mathrm{1}}/\mathit{\epsilon })}{(\mathrm{1}+r_t^{n+\mathrm{1}})}}.\end{array}$$

Due to the presence of non-linear terms, the discrete system (11a)–(11b) 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; Smolarkiewicz et al., 2014, 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 ϑⁿ⁺¹ defined in Eq. (12) that is composed of the already completed $r_k^{n+\mathrm{1}}$. At the beginning of the iterative execution of Eqs. (11a)–(11b), a first guess θ⁰ for θ^⋆ is provided as the explicit solution of full potential temperature:

$$\begin{array}{ll}{\displaystyle }{\mathit{\theta }}_{\mathit{i}}^{\mathrm{0}}={\widehat{\mathit{\theta }}}_{\mathit{i}}={\mathcal{A}}_{\mathit{i}}& {\displaystyle }({\mathit{\theta }}^n+\mathit{\delta }t{P}^{{\mathit{\theta }}^{\prime }}{|}^n,(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{⟂}{|}^{n+\mathrm{1}/\mathrm{2}},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^n,\\ \text{(13)}& {\displaystyle }& {\displaystyle }(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{n+\mathrm{1}},\mathit{\delta }t),\end{array}$$

whereas a first guess u⁰ for u^⋆ is prescribed by linear extrapolation of u from the tⁿ⁻¹ and tⁿ 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 ${\mathit{\theta }}_{\mathit{i}}^{\prime }$ 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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{u}+\mathrm{0.5}\mathit{\delta }t\mathit{f}\times \mathit{u}-(\mathrm{0.5}\mathit{\delta }t{)}^{\mathrm{2}}{\displaystyle \frac{\mathit{g}}{{\mathit{\theta }}_{\mathrm{a}}}}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\cdot \mathrm{\nabla }{\mathit{\theta }}_{\mathrm{a}}=\\ {\displaystyle }& {\displaystyle }\widehat{\mathit{u}}-\mathrm{0.5}\mathit{\delta }t[{\displaystyle \frac{\mathit{g}}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}\left({\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}-{\mathit{\vartheta }}^{n+\mathrm{1}}{\mathit{\theta }}_{\mathrm{a}}-{\mathit{\vartheta }}^{n+\mathrm{1}}\widehat{{\mathit{\theta }}^{\prime }}\right)-\mathit{f}\times {\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}^{\star }}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}{\mathit{u}}_{\mathrm{a}}\\ {\displaystyle }& {\displaystyle }-{\mathcal{M}}^{\prime }\left({\mathit{u}}^{\star },{\mathit{u}}_{\mathrm{a}},{\displaystyle \frac{{\mathit{\theta }}_{\mathit{\rho }}^{\star }}{{\mathit{\theta }}_{\mathit{\rho }\mathrm{a}}}}\right)]-\mathrm{0.5}\mathit{\delta }t{\mathit{\theta }}_{\mathit{\rho }}^{\star }\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}^{\prime }\\ \text{(14)}& {\displaystyle }& {\displaystyle }\equiv \widehat{\widehat{\mathit{u}}}-\mathrm{0.5}\mathit{\delta }t{\mathit{\theta }}_{\mathit{\rho }}^{\star }\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}^{\prime }.\end{array}$$

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 $\widehat{\widehat{\mathit{u}}}$, and the pressure gradient term with the lagged coefficient ${\mathit{\theta }}_{\mathit{\rho }}^{\star }$. Defining L as the linear operator acting on u on the LHS and L⁻¹ its inverse, Eq. (14) can be symbolized as

$$\begin{array}{}\text{(15)}& \mathbf{L}\mathit{u}=\widehat{\widehat{\mathit{u}}}-\mathrm{0.5}\mathit{\delta }t{\mathit{\theta }}_{\mathit{\rho }}^{\star }\stackrel{\mathrm{̃}}{\mathbf{G}}\mathrm{\nabla }{\mathit{\phi }}^{\prime },\end{array}$$

and

$$\begin{array}{}\text{(16)}& \mathit{u}=\stackrel{\mathrm{ˇ}}{\stackrel{\mathrm{ˇ}}{\mathit{u}}}-\mathbf{C}\mathrm{\nabla }{\mathit{\phi }}^{\prime },\end{array}$$

respectively, where $\stackrel{\mathrm{ˇ}}{\stackrel{\mathrm{ˇ}}{\mathit{u}}}={\mathbf{L}}^{-\mathrm{1}}\widehat{\widehat{\mathit{u}}}$, and $\mathbf{C}={\mathbf{L}}^{-\mathrm{1}}\mathrm{0.5}\mathit{\delta }t{\mathit{\theta }}_{\mathit{\rho }}^{\star }\stackrel{\mathrm{̃}}{\mathbf{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. 1b) 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)–(1b) with an auxiliary 3-D 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 (Eq. 1b) combined with the respective advective forms of the mass continuity Eq. (1a), thermodynamic Eq. (1b), and water vapour Eq. (1b) equations; see Smolarkiewicz et al. (2014, 2017) for further details. The governing equation for φ^′ is then implemented in conservation form consistent with Eqs. (1a)–(1b) and reads

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial \mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\phi }}^{\prime }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\phi }}^{\prime }\right)=\\ {\displaystyle }& {\displaystyle }\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}[-{\displaystyle \frac{R_{\mathrm{d}}}{c_{\mathrm{vd}}}}{\displaystyle \frac{\mathit{\phi }}{\mathcal{G}}}\mathrm{\nabla }\cdot (\mathcal{G}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u})-{\displaystyle \frac{\mathrm{1}}{\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}}}\mathrm{\nabla }\cdot (\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}{\mathit{\phi }}_{\mathrm{a}})\\ \text{(17)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{\phi }}_{\mathrm{a}}}{\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}}}\mathrm{\nabla }\cdot \left(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\right)+{\displaystyle \frac{R_{\mathrm{d}}}{c_{\mathrm{vd}}}}\mathit{\phi }\mathrm{\Pi }],\end{array}$$

where $\mathit{\phi }={\mathit{\phi }}_{\mathrm{a}}+{\mathit{\phi }}^{\prime }$ and

$$\begin{array}{}\text{(18)}& \mathrm{\Pi }=\left({\displaystyle \frac{P^{{\mathit{\theta }}^{\prime }}}{\mathit{\theta }}}+{\displaystyle \frac{P^{r_{\mathrm{v}}}/\mathit{\epsilon }}{\mathrm{1}+r_{\mathrm{v}}/\mathit{\epsilon }}}\right).\end{array}$$

We interpret Eq. (17) in terms of the generalized transport Eq. (6) for $\mathrm{\Psi }\equiv {\mathit{\phi }}^{\prime }$ (see Table 2), while setting

$$\begin{array}{}\text{(19)}& P^{{\mathit{\phi }}^{\prime }}=(R_{\mathrm{d}}/c_{\mathrm{vd}})\mathit{\phi }\mathrm{\Pi },\end{array}$$

and ${\mathcal{R}}^{{\mathit{\phi }}^{\prime }}$ as the remaining RHS consisting of the three divergence operators. Importantly, the $P^{{\mathit{\phi }}^{\prime }}$ given by Eq. (19) describes the pressure adjustment from the physics parameterization tendencies $P^{{\mathit{\theta }}^{\prime }}$ and $P^{r_{\mathrm{v}}}$. We have implemented the integration of Eq. (17) using the template semi-implicit scheme (7) with parameters $a^{{\mathit{\phi }}^{\prime }}=(\mathrm{1}-\mathit{\alpha })$ and $b^{{\mathit{\phi }}^{\prime }}=\mathit{\alpha }$, where $\mathit{\alpha }\in [\mathrm{0.5},\mathrm{1.0}]$. In the limit α=1.0, the template Eq. (7) represents the first-order backward Euler scheme with regard to ${\mathcal{R}}^{{\mathit{\phi }}^{\prime }}$. For α=0.5, the template Eq. (7) becomes the second-order trapezoidal scheme; in practice we may use weak off-centring (e.g. α=0.51) for regularization. For any specification of α, coupling with the solution procedure of the fully compressible Eqs. (1a)–(1b) is implemented through Eq. (16) that enters into the three occurrences of u on the RHS of Eq. (17) in ${\mathcal{R}}^{{\mathit{\phi }}^{\prime }}{|}^{n+\mathrm{1}}$. Furthermore, as indicated in Eq. (7), the (explicit) forward Euler scheme is used with regard to the forcing $P^{{\mathit{\phi }}^{\prime }}$. Reorganizing terms finally yields the elliptic Helmholtz equation for the pressure perturbation variable φ^′ at the future time level tⁿ⁺¹, which can be written compactly as

$$\begin{array}{ll}{\displaystyle }\mathrm{0}& {\displaystyle }=-\sum _{\mathrm{\ell }=\mathrm{1}}^{\mathrm{3}}\left({\displaystyle \frac{A_{\mathrm{\ell }}^{\star }}{{\mathit{\zeta }}_{\mathrm{\ell }}}}\mathrm{\nabla }\cdot {\mathit{\zeta }}_{\mathrm{\ell }}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}(\stackrel{\mathrm{ˇ}}{\stackrel{\mathrm{ˇ}}{\mathit{u}}}-\mathbf{C}\mathrm{\nabla }{\mathit{\phi }}^{\prime })\right)-B^{\star }({\mathit{\phi }}^{\prime }-\widehat{{\mathit{\phi }}^{\prime }})\\ \text{(20)}& {\displaystyle }& {\displaystyle }\equiv \mathcal{L}\left({\mathit{\phi }}^{\prime }\right)-R,\end{array}$$

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 ℒ(φ^′) and R refer to the implicit and explicit parts of the equation, respectively. The coefficients $A_{\mathrm{\ell }}^{\star }$, ζ_ℓ, and B^⋆ are defined accordingly as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{A}_{\mathrm{1}}^{\star }=\mathrm{1},A_{\mathrm{2}}^{\star }=A_{\mathrm{3}}^{\star }={\displaystyle \frac{c_{\mathrm{vd}}}{R_{\mathrm{d}}}}{\displaystyle \frac{{\mathit{\phi }}_{\mathrm{a}}}{\mathit{\phi }}},{\mathit{\zeta }}_{\mathrm{1}}=\mathcal{G},{\mathit{\zeta }}_{\mathrm{2}}=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{\mathit{\phi }}_{\mathrm{a}},\\ \text{(21)}& {\displaystyle }& {\displaystyle }{\mathit{\zeta }}_{\mathrm{3}}=\mathcal{G}{\mathit{\rho }}_{\mathrm{d}},B^{\star }={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }\mathit{\delta }t}}{\displaystyle \frac{c_{\mathrm{vd}}}{R_{\mathrm{d}}\mathit{\phi }}},\end{array}$$

and $\widehat{{\mathit{\phi }}^{\prime }}=\mathcal{A}({\mathit{\phi }}^{\prime n}+(\mathrm{1}-\mathit{\alpha })\mathit{\delta }t{\mathcal{R}}^{{\mathit{\phi }}^{\prime }}{|}^n+\mathit{\delta }t{P}^{{\mathit{\phi }}^{\prime }}{|}^n,(\mathit{v}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{⟂}{|}^{n+\mathrm{1}/\mathrm{2}}$, $(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^n,(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{n+\mathrm{1}},\mathit{\delta }t)$. The Helmholtz Eq. (20) extends the implementations in Smolarkiewicz et al. (2016, 2017) and Kühnlein and Smolarkiewicz (2017) 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 ℒ and enables rapid convergence of the GCR solver. The preconditioner 𝒫≈ℒ discards the off-diagonal entries of the matrix ${\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathbf{C}$. Inversion of 𝒫 then utilizes a weighted line Jacobi method explained in Appendix B. In the GCR solver and the preconditioner 𝒫, the coefficients $A_{\mathrm{2}}^{\star }$, $A_{\mathrm{3}}^{\star }$, and B^⋆ in Eq. (21) depend on $\mathit{\phi }={\mathit{\phi }}_{\mathrm{a}}+{\mathit{\phi }}^{\prime }$, 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 𝒪(δt²) 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³) 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; Kühnlein and Smolarkiewicz, 2017). 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)–(1b), and thus the semi-implicit model time step δt can be selected according to the stability of the advective transport scheme ${\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},{\mathit{V}}^{n+\mathrm{1}/\mathrm{2}},G^n,G^{n+\mathrm{1}},\mathit{\delta }t)$ in Eq. (7).

2.1.3 Spatial discretization

The discretization framework of IFS-FVM combines a structured-grid FD–FV method in the vertical with an unstructured⁵ FV approach in the horizontal (Smolarkiewicz et al., 2016). 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 spatial discretization of the applied MPDATA schemes in IFS-FVM, symbolized by the operator 𝒜_i in Eq. (7), is described in Kühnlein and Smolarkiewicz (2017).

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

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 𝒱_i and 𝒱_j, respectively.

For a differentiable vector field A, the Gauss divergence theorem ${\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot \mathit{A}={\int }_{\partial \mathrm{\Omega }}\mathit{A}\cdot \mathit{n}$ applied over the control volume surrounding node $\mathit{i}=(k,i)$ in 3-D computational space is given as

$$\begin{array}{}\text{(22)}& (\mathrm{\nabla }\cdot \mathit{A}{)}_{\mathit{i}}={\displaystyle \frac{\mathrm{1}}{{\mathcal{V}}_i}}\sum _{j=\mathrm{1}}^{l\left(i\right)}{A}_{k,j}^{⟂}{S}_j+{\displaystyle \frac{A_{k+\mathrm{1}/\mathrm{2},i}^z-A_{k-\mathrm{1}/\mathrm{2},i}^z}{\mathit{\delta }z}},\end{array}$$

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⁶. The geometric quantities S_j, 𝒱_i, and δz of the computational space are independent of height. The fields $A_{k,j}^{⟂}$ and $A_{k+\mathrm{1}/\mathrm{2},i}^z$ are interpreted as the mean normal components of the vector A at the horizontal and vertical cell faces, respectively. Elementary approximations are given as $A_{k,j}^{⟂}=\mathrm{0.5}{\mathit{n}}_j\cdot ({\mathit{A}}_{k,i}+{\mathit{A}}_{k,j})$ in the horizontal and $A_{k+\mathrm{1}/\mathrm{2},i}^z=\mathrm{0.5}(A_{k+\mathrm{1},i}^z+A_{k,i}^z)$ 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}^{⟂}=\mathrm{0.5}{\mathit{n}}_j\cdot {\mathit{A}}_{k,j}$ and $A_{k+\mathrm{1}/\mathrm{2},i}^z=\mathrm{0.5}{A}_{k+\mathrm{1},i}^z$. 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

$$\begin{array}{ll}{\displaystyle }(\mathrm{\nabla }\mathrm{\Psi }{)}_{k,i}=(& {\displaystyle \frac{\mathrm{1}}{{\mathcal{V}}_i}}\sum _{j=\mathrm{1}}^{l\left(i\right)}\mathrm{0.5}{\mathrm{\Psi }}_{k,j}{S}_j^x,{\displaystyle \frac{\mathrm{1}}{{\mathcal{V}}_i}}\sum _{j=\mathrm{1}}^{l\left(i\right)}\mathrm{0.5}{\mathrm{\Psi }}_{k,j}{S}_j^y,\\ \text{(23)}& {\displaystyle }& {\displaystyle }\mathrm{0.5}{\displaystyle \frac{{\mathrm{\Psi }}_{k+\mathrm{1},i}-{\mathrm{\Psi }}_{k-\mathrm{1},i}}{\mathit{\delta }z}}),\end{array}$$

where $S_j^x$ and $S_j^y$ denote the x and y components of the oriented surface element S_j=S_jn_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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left({\displaystyle \frac{\mathrm{1}}{\mathcal{G}}}\mathrm{\nabla }\cdot \left(\mathcal{G}{\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\right)\right)}_{k,i}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{0.5}}{{\mathcal{G}}_{k,i}{\mathcal{V}}_i}}\sum _{j=\mathrm{1}}^{l\left(i\right)}{\mathcal{G}}_{k,j}({\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{1}}u{)}_{k,j}{S}_j^x+{\mathcal{G}}_{k,j}({\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{2}}v{)}_{k,j}{S}_j^y\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{0.5}}{{\mathcal{G}}_{k,i}\mathit{\delta }z}}({\mathcal{G}}_{k+\mathrm{1},i}{\left({\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{3}}u+{\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{3}}v+{\stackrel{\mathrm{̃}}{G}}_{\mathrm{3}}^{\mathrm{3}}w\right)}_{k+\mathrm{1},i}\\ \text{(24)}& {\displaystyle }& {\displaystyle }-{\mathcal{G}}_{k-\mathrm{1},i}{\left({\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{3}}u+{\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{3}}v+{\stackrel{\mathrm{̃}}{G}}_{\mathrm{3}}^{\mathrm{3}}w\right)}_{k-\mathrm{1},i}).\end{array}$$

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 distributed-memory 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 equal-region horizontal domain decomposition of IFS. The nearest-neighbour communications enabling the FV stencil operations in the parallel horizontal domain are performed on overlap regions (halos) between partitions, as well as a “periodic overlap” for the east–west boundary on the globe. 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-per-byte 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.

2.2 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⁷ 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 finite-element (FE) approach to discretize the integral operator in the vertical direction (Untch and Hortal, 2004)⁸. 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) with the Lorenz staggering. Prognostic variables of the HPEs and other main formulation features of IFS-ST are given in Table 1. Wedi et al. (2015) provide a recent overview of IFS-ST and a comprehensive list of references, while the official IFS documentation and changes with model cycles can be found on the ECMWF website (https://www.ecmwf.int/, last access: 6 February 2019).

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

$$\begin{array}{}\text{(25)}& X_{\mathrm{A}}^{n+\mathrm{1}}=X_{\mathrm{D}}^n+{\mathcal{N}}_{\mathrm{M}}^{n+\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathcal{L}}_{\mathrm{D}}^n+{\mathcal{L}}_{\mathrm{A}}^{n+\mathrm{1}}\right).\end{array}$$

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 ℒ symbolizes the linear operator and 𝒩 the non-linear residual $\mathcal{N}=\mathcal{M}-\mathcal{L}$, with ℳ being the full non-linear model. The operator ℒ results from the linearization of ℳ 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 ${\mathcal{L}}_{\mathrm{A}}^{n+\mathrm{1}}$ with ${\mathcal{L}}_{\mathrm{A}}^n$. Generally, $X_{\mathrm{D}}^n$ 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 compare the computational performance of the HPE and NH formulations of IFS-ST against IFS-FVM.

2.3 Some aspects of physics–dynamics coupling

The IFS physics parameterization package at ECMWF is applied in the same configuration throughout the medium-range, 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.

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.

In IFS-ST, the physics–dynamics coupling employs the SLAVEPP (semi-Lagrangian averaging of physical parameterization) scheme (Wedi, 1999). SLAVEPP targets second-order accuracy by averaging tendencies from selected physics parameterizations along the SL trajectory at the midpoint in space–time between the departure point at tⁿ and arrival point at tⁿ⁺¹. Thereby, the tendencies at tⁿ⁺¹ 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ⁿ. Therefore, the fields that enter the parameterizations are from tⁿ, and there is no averaging between tendencies from tⁿ and tⁿ⁺¹. 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ⁿ 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ⁿ. 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 in IFS-ST (see Table 1)¹⁰, 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 longitude–latitude 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+\mathrm{\Delta }{t}_{\mathrm{phys}}\equiv t^N+N_{\mathrm{s}}\mathit{\delta }t$ can be written as $\mathrm{\ell }=\mathrm{1},N_{\mathrm{s}}:$

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Psi }}_{\mathit{i}}(t^N+\mathrm{\ell }\mathit{\delta }t)& {\displaystyle }={\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},\mathit{V}(t^N+(\mathrm{\ell }\mathit{\delta }t-\mathrm{0.5})),\\ {\displaystyle }& {\displaystyle }G(t^N+(\mathrm{\ell }-\mathrm{1}\left)\mathit{\delta }t\right),G(t^N+\mathrm{\ell }\mathit{\delta }t),\mathit{\delta }t)\\ \text{(26)}& {\displaystyle }& {\displaystyle }+b^{\mathrm{\Psi }}\mathit{\delta }t{\mathcal{R}}_{\mathit{i}}^{\mathrm{\Psi }}(t^N+\mathrm{\ell }\mathit{\delta }t),\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Psi }}& {\displaystyle }=\mathrm{\Psi }(t^N+(\mathrm{\ell }-\mathrm{1}\left)\mathit{\delta }t\right)+a^{\mathrm{\Psi }}\mathit{\delta }t{\mathcal{R}}^{\mathrm{\Psi }}(t^N+(\mathrm{\ell }-\mathrm{1}\left)\mathit{\delta }t\right)\\ \text{(27)}& {\displaystyle }& {\displaystyle }+\mathit{\delta }t{P}^{\mathrm{\Psi }}(t^N,\mathrm{\Delta }{t}_{\mathrm{phys}}).\end{array}$$

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 (Kühnlein and Smolarkiewicz, 2017). 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¹¹, 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¹², 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.

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

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 spectral-transform 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 under 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 sigma–pressure 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 interpolation¹³. 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 horizontal–vertical 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.

Figure 3 Dry baroclinic instability at day 10: panels (a)–(d) show pressure on the lowest full level (hPa) obtained with IFS-FVM and IFS-ST, while (e) and (f) depict the corresponding difference between the solutions. Panels (a), (c), (e) and (b), (d), (f) are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

Figure 4 Dry baroclinic instability at day 10: panels (a)–(d) show meridional wind v (m s⁻¹) along a zonal height cross section at 50^∘ N obtained with IFS-FVM and IFS-ST, while (e) and (f) depict the difference of v (m s⁻¹) between their solutions. Panels (a), (c), (e) and (b), (d), (f) are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

Figure 5 Dry baroclinic instability at day 15: pressure on the lowest full level (hPa) obtained with IFS-FVM (default backward Euler scheme with α=1 or trapezoidal scheme with weak off-centring α=0.51 in the integration of Eq. (17) for the Exner pressure perturbation φ^′) and IFS-ST using the (O320,TCo319) grid.

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. 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)–(1b) 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.

Figure 6 Dry baroclinic instability at day 15 when the triggering of the baroclinic instability was applied in the Northern Hemisphere only: pressure on the lowest full level (hPa) obtained with IFS-FVM and IFS-ST using the (O320,TCo319) grid.

3.1 Results for dry simulations

Figures 3 and 4 present the horizontal cross section of near-surface 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 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 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 wavenumber-four 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.

Figure 7 Dry 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 4 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.

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¹⁴. The slope of the IFS-ST and IFS-FVM spectra with respect to wavenumber l is somewhat shallower than l⁻³ 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.

Figure 8 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 parameterization. The upper (lower) panels show shaded contours from the IFS-FVM (IFS-ST) simulations, overlaid by the IFS-ST (IFS-FVM) black contour line of 0.5 mm day⁻¹. Panels (a), (c) and (b), (d) are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

Figure 9 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 parameterization. The upper (lower) panels show shaded contours from the IFS-FVM (IFS-ST) simulations, overlaid by the IFS-ST (IFS-FVM) black contour line of 0.5 mm day⁻¹. Panels (a), (c) and (b), (d) are for the (O160,TCo159) and (O320,TCo319) horizontal grids, respectively.

3.2 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¹⁵ 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 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 10 supplements the precipitation plots with the corresponding pressure field on day 15. In addition to the standard configurations of IFS-FVM and IFS-ST, for which the physics parameterization is evaluated every dynamic time step N_s=1, Fig. 10 also provides the IFS-FVM result with subcyling (middle panel) whereby the parameterizations are evaluated every N_s=3 semi-implicit time step δt; see Sect. 2.3 for a 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 parameterization seems to have only a small impact on the solution. Furthermore, none of the simulations show 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 cloud parameterization (in contrast to the dry results shown before in Fig. 5). The analysis 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 moist-precipitating FV formulations along with the coupling to the IFS physics parameterization.

Figure 10 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 parameterization. The IFS-FVM results in (a) and (b) employed the standard coupling at every time step N_s=1 or subcycling of dynamics for three time steps N_s=3, respectively. The simulations were performed with the (O320,TCo319) grid.

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

Figure 12 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 4 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 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 assess 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¹⁶. 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.

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.

Figure 13 Elapsed time to run 1 day of the dry baroclinic instability benchmark similar to the current HRES configuration at ECMWF, i.e. the three different models – the nonhydrostatic IFS-FVM designated as FV(NH), the hydrostatic IFS-ST designated as ST(H), and the nonhydrostatic IFS-ST designated as ST(NH) – are set up for the O1280/TCo1279 horizontal grid (corresponding to about 9 km grid spacing) with 137 stretched vertical levels and employ 350 nodes of ECMWF's Cray XC40.

4 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, deep-atmosphere 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 finite-volume integration of the fully compressible equations of the novel IFS-FVM considering comprehensive moist-precipitating 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 MPDATA, 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 co-located 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 spectral-transform 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 medium-range NWP at convection-resolving resolutions (Kühnlein and Smolarkiewicz, 2019).

Appendix A: Horizontal–vertical splitting of the NFT advective transport

We consider the advection operator 𝒜_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 global meshes, 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)–(1b) proceeds exactly as explained in Sect. 2.1.2, but with the 3-D NFT advection operator 𝒜_i split into purely horizontal ${\mathcal{A}}_{\mathit{i}}^{xy}$ and vertical ${\mathcal{A}}_{\mathit{i}}^z$ schemes, respectively. For each model time step δt, these are applied in the sequence ${\mathcal{A}}_{\mathit{i}}^z\to {\mathcal{A}}_{\mathit{i}}^{xy}\to {\mathcal{A}}_{\mathit{i}}^z$ 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

$$\begin{array}{ll}{\displaystyle }{\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{\left[\mathrm{1}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^z({\mathit{\rho }}_{\mathrm{d}}^n,(v^z\mathcal{G}{)}^{n+\mathrm{1}/\mathrm{2}},{\mathcal{G}}^n,{\mathcal{G}}^{\left[\mathrm{1}\right]},\mathrm{0.5}\mathit{\delta }t),\\ {\displaystyle }{\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{\left[\mathrm{2}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^{xy}({\mathit{\rho }}_{\mathrm{d}}^{\left[\mathrm{1}\right]},({\mathit{v}}_{\mathrm{h}}\mathcal{G}{)}^{n+\mathrm{1}/\mathrm{2}},{\mathcal{G}}^{\left[\mathrm{1}\right]},{\mathcal{G}}^{\left[\mathrm{2}\right]},\mathit{\delta }t),\\ \text{(A1)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{n+\mathrm{1}}={\mathit{\rho }}_{\mathrm{d}\mathit{i}}^{\left[\mathrm{3}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^z({\mathit{\rho }}_{\mathrm{d}}^{\left[\mathrm{2}\right]},(v^z\mathcal{G}{)}^{n+\mathrm{1}/\mathrm{2}},{\mathcal{G}}^{\left[\mathrm{2}\right]},{\mathcal{G}}^{n+\mathrm{1}},\mathrm{0.5}\mathit{\delta }t),\end{array}$$

which provides the updated densities ${\mathit{\rho }}_{\mathrm{d}}^{\left[\mathrm{1}\right]}$, ${\mathit{\rho }}_{\mathrm{d}}^{\left[\mathrm{2}\right]}$, ${\mathit{\rho }}_{\mathrm{d}}^{n+\mathrm{1}}$ and accumulates normal mass fluxes (v^z𝒢ρ_d)^[1], $(v_{\mathrm{h}}^{⟂}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{2}\right]}$, (v^z𝒢ρ_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 $\stackrel{\mathrm{̃}}{\mathrm{\Psi }}$ (Eq. 8) as

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Psi }}_{\mathit{i}}^{\left[\mathrm{1}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^z(\stackrel{\mathrm{̃}}{\mathrm{\Psi }},(v^z\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{1}\right]},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^n,(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{1}\right]},\mathrm{0.5}\mathit{\delta }t),\\ {\displaystyle }{\mathrm{\Psi }}_{\mathit{i}}^{\left[\mathrm{2}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^{xy}({\mathrm{\Psi }}^{\left[\mathrm{1}\right]},(v_{\mathrm{h}}^{⟂}\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{2}\right]},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{1}\right]},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{2}\right]},\mathit{\delta }t),\\ \text{(A2)}& {\displaystyle }{\mathrm{\Psi }}_{\mathit{i}}^{\left[\mathrm{3}\right]}& {\displaystyle }={\mathcal{A}}_{\mathit{i}}^z({\mathrm{\Psi }}^{\left[\mathrm{2}\right]},(v^z\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{3}\right]},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{\left[\mathrm{2}\right]},(\mathcal{G}{\mathit{\rho }}_{\mathrm{d}}{)}^{n+\mathrm{1}},\mathrm{0.5}\mathit{\delta }t),\end{array}$$

where ${\widehat{\mathrm{\Psi }}}_{\mathit{i}}\equiv {\mathrm{\Psi }}_{\mathit{i}}^{\left[\mathrm{3}\right]}$. In Eqs. (A1) and (A2), the implementation of the horizontal advection transport 𝒜^xy follows the horizontal part of the unstructured-mesh FV MPDATA of Kühnlein and Smolarkiewicz (2017). The vertical scheme 𝒜^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 Kühnlein and Smolarkiewicz (2017). This is particularly due to the integration of the vertical parts ${\mathcal{A}}_{\mathit{i}}^z$ with δt∕2 each, which mitigates the vertical stability restriction while not adding any significant computational cost¹⁷. Overall, the horizontal–vertical splitting of 𝒜_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 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 𝒜_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 φ^′:

$$\begin{array}{}\text{(B1)}& \mathcal{P}\left(e\right)=\widehat{r},\end{array}$$

where $\widehat{r}$ denotes the residual error of Eq. (20). The preconditioning operator 𝒫 is then decomposed into vertical and horizontal parts (Smolarkiewicz and Margolin, 2000), and the residual problem is solved iteratively according to

$$\begin{array}{}\text{(B2)}& {\mathcal{P}}_z\left(e^{\mathit{\mu }+\mathrm{1}}\right)+{\mathcal{P}}_{\mathrm{h}}\left(e^{\mathit{\mu }}\right)-\widehat{r}=\mathrm{0},\end{array}$$

where μ numbers the iterations, of which there are typically two. The vertical part 𝒫_z is inverted directly with a tridiagonal algorithm. The horizontal part 𝒫_h is lagged behind, except for its diagonal entries. The actual implementation is given as

$$\begin{array}{}\text{(B3)}& {\mathcal{P}}_z\left(e^{\mathit{\mu }+\mathrm{1}}\right)+{\mathcal{P}}_{\mathrm{h}}\left(e^{\mathit{\mu }}\right)+\mathcal{D}(e^{\mathit{\mu }+\mathrm{1}}-e^{\mathit{\mu }})-\widehat{r}=\mathrm{0},\end{array}$$

where 𝒟 is the diagonal coefficient of 𝒫_h, specified as

$$\begin{array}{}\text{(B4)}& {\mathcal{D}}_{k,i}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathcal{V}}_{\mathrm{i}}}}\sum _{\mathrm{\ell }=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{{A_{\mathrm{\ell }}^{\star }}_{k,i}}{{{\mathit{\zeta }}_{\mathrm{\ell }}}_{k,i}}}\sum _{j=\mathrm{1}}^{l\left(i\right)}{\displaystyle \frac{{{\mathit{\zeta }}_{\mathrm{\ell }}}_{k,j}}{{\mathcal{V}}_j}}\left({\mathcal{B}}_{k,j}^{\mathrm{11}}{S_j^x}^{\mathrm{2}}+{\mathcal{B}}_{k,j}^{\mathrm{22}}{S_j^y}^{\mathrm{2}}\right),\end{array}$$

with ℬ¹¹ and ℬ²² referring to the diagonal entries of ${\stackrel{\mathrm{̃}}{\mathbf{G}}}^{\mathrm{T}}\mathbf{C}$. Subsequently, Eq. (B3) is executed as

$$\begin{array}{}\text{(B5)}& e^{\mathit{\mu }+\mathrm{1}}=\mathit{\omega }{\left[\mathcal{D}+{\mathcal{P}}_z\right]}^{-\mathrm{1}}(\mathcal{D}{e}^{\mathit{\mu }}-{\mathcal{P}}_{\mathrm{h}}(e^{\mathit{\mu }})+\widehat{r})+\left(\mathrm{1}-\mathit{\omega }\right)e^{\mathit{\mu }}\end{array}$$

with the weight ω=0.7.

Appendix C: Geospherical framework, generalized terrain-following vertical coordinate, shallow- and deep-atmosphere equations

IFS-FVM's ability to accommodate complex mesh geometries results from two aspects of its formulation: the horizontal unstructured-mesh FV discretization and generalized curvilinear coordinate mappings embedded in a geospherical framework (Prusa and Smolarkiewicz, 2003; Szmelter and Smolarkiewicz, 2010).

In the geospherical curvilinear coordinate framework of Prusa and Smolarkiewicz (2003), the vector $\mathit{u}=[u,v,w{]}^{\mathrm{T}}$ represents the physical velocity with zonal, meridional, and vertical components aligned at every point of the spherical shell with axes of a local Cartesian frame (marked with the superscript “c”) tangent to the lower surface (r=a); here r is the radial component of the vector radius, and a is the radius of the sphere. Relations between the local Cartesian and the geospherical frame are therefore dx_c=rcos ϕ dλ, dy_c=r dϕ, and $z_{\mathrm{c}}=r-a$, where λ and ϕ denote longitude and latitude, respectively, in radians.

Consistent with Prusa and Smolarkiewicz (2003) (but not with their notation), we define a set of geospherical coordinates of the physical space S_p as $\stackrel{\mathrm{̃}}{x}=a\mathit{\lambda }$, $\stackrel{\mathrm{̃}}{y}=a\mathit{\varphi }$, $\stackrel{\mathrm{̃}}{z}=z_{\mathrm{c}}$ ($\stackrel{\mathrm{̃}}{x},\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z}$ in units of metres). The latter are related to the curvilinear coordinates $\mathit{x}=[x,y,z{]}^{\mathrm{T}}$ of the computational space S_t (see Sect. 2.1.1) by the general transformation

$$\begin{array}{}\text{(C1)}& (t,\mathit{x})=\left(\stackrel{\mathrm{̃}}{t},\mathcal{F}(\stackrel{\mathrm{̃}}{t},\stackrel{\mathrm{̃}}{\mathit{x}})\right),\end{array}$$

where $\mathcal{F}(\stackrel{\mathrm{̃}}{t},\stackrel{\mathrm{̃}}{\mathit{x}})$ represents a bijective map between the physical and computational systems (Prusa and Smolarkiewicz, 2003; Kühnlein et al., 2012). A default mapping in IFS-FVM uses no stretching with respect to the horizontal positions of the unstructured computational mesh $x\equiv \stackrel{\mathrm{̃}}{x}$, $y\equiv \stackrel{\mathrm{̃}}{y}$, combined with a height-based terrain-following vertical coordinate of the general form $z=z(\stackrel{\mathrm{̃}}{x},\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$. The most straightforward specification of the mapping $z=z(\stackrel{\mathrm{̃}}{x},\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$ is a basic terrain-following vertical coordinate given by means of analytical functions; e.g. Gal-Chen and Somerville (1975). However, the implemented general form $z=z(\stackrel{\mathrm{̃}}{x},\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$ admits variable vertical stretching with horizontal location, implicit–explicit smoothing of coordinate levels (Schär et al., 2002; Klemp, 2011), and hybrid specifications. Further note that in the numerical experiments of Sect. 3, IFS-FVM employed the vertical levels defined by the height of the hybrid sigma–pressure coordinate levels of IFS-ST. Overall, under the described coordinate mappings, the 3×3 coefficient matrix $\stackrel{\mathrm{̃}}{\mathbf{G}}$ employed in the formalism of Sect. 2.1.1 is given as

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{̃}}{\mathbf{G}}& {\displaystyle }=\left[\begin{array}{ccc}{\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{1}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{2}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{1}}^{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{1}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{2}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{2}}^{\mathrm{3}}\\ {\stackrel{\mathrm{̃}}{G}}_{\mathrm{3}}^{\mathrm{1}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{3}}^{\mathrm{2}}& {\stackrel{\mathrm{̃}}{G}}_{\mathrm{3}}^{\mathrm{3}}\end{array}\right]=\\ \text{(C2)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{ccc}\left(\mathrm{\Gamma }\cos \right(\stackrel{\mathrm{̃}}{y}/a){)}^{-\mathrm{1}}& \mathrm{0}& \left(\mathrm{\Gamma }\cos \right(\stackrel{\mathrm{̃}}{y}/a){)}^{-\mathrm{1}}\partial z/\partial \stackrel{\mathrm{̃}}{x}\\ \mathrm{0}& {\mathrm{\Gamma }}^{-\mathrm{1}}& {\mathrm{\Gamma }}^{-\mathrm{1}}\partial z/\partial \stackrel{\mathrm{̃}}{y}\\ \mathrm{0}& \mathrm{0}& \partial z/\partial \stackrel{\mathrm{̃}}{z}\end{array}\right],\end{array}$$

where

$$\begin{array}{}\text{(C3)}& \mathrm{\Gamma }=\mathrm{1}+\mathit{\gamma }\stackrel{\mathrm{̃}}{z}/a,\end{array}$$

with γ=0 and γ=1 for the shallow- and deep-atmosphere form of the governing Eqs. (1a)–(1b), respectively, and the indices 1, 2, and 3 correspond to x, y, and z components. The corresponding Jacobian of the coordinate mappings 𝒢, which appears in the system (1a)–(1b) as well as other equations, is

$$\begin{array}{}\text{(C4)}& \mathcal{G}={\mathrm{\Gamma }}^{\mathrm{2}}\cos (\stackrel{\mathrm{̃}}{y}/a)(\partial z/\partial \stackrel{\mathrm{̃}}{z}{)}^{-\mathrm{1}}.\end{array}$$

The inverse metrics $\partial z/\partial \stackrel{\mathrm{̃}}{x}$, $\partial z/\partial \stackrel{\mathrm{̃}}{y}$ and $\partial z/\partial \stackrel{\mathrm{̃}}{z}$ in Eqs. (C2) and (C4) are computed consistently with the FV discretization of Sect. 2.1.3 and using the Kronecker-delta identity; e.g. Kühnlein et al. (2012).

Furthermore, in the momentum Eq. (1b), the components of the Coriolis acceleration are

$$\begin{array}{ll}{\displaystyle }-\mathit{f}\times \mathit{u}=[& {\displaystyle }v{f}_{\mathrm{0}}\sin (\stackrel{\mathrm{̃}}{y}/a)-\mathit{\gamma }w{f}_{\mathrm{0}}\cos (\stackrel{\mathrm{̃}}{y}/a),\\ \text{(C5)}& {\displaystyle }& {\displaystyle }-u{f}_{\mathrm{0}}\sin (\stackrel{\mathrm{̃}}{y}/a),\mathit{\gamma }u{f}_{\mathrm{0}}\cos (\stackrel{\mathrm{̃}}{y}/a){]}^{\mathrm{T}},\end{array}$$

where $f_{\mathrm{0}}=\mathrm{2}\left|\mathbf{\Omega }\right|$, and the metric forcings due to the curvature of the sphere (i.e. component-wise Christoffel terms associated with the advective derivative of the physical velocity) are

$$\begin{array}{ll}{\displaystyle }\mathcal{M}\left(\mathit{u}\right)=\left(\mathrm{\Gamma }a{)}^{-\mathrm{1}}\right[& {\displaystyle }\tan (\stackrel{\mathrm{̃}}{y}/a)uv-\mathit{\gamma }uw,-\tan (\stackrel{\mathrm{̃}}{y}/a)uu\\ \text{(C6)}& {\displaystyle }& {\displaystyle }-\mathit{\gamma }vw,\mathit{\gamma }(uu+vv){]}^{\mathrm{T}}.\end{array}$$

The buoyancy term of Eq. (1b) contains the gravitational acceleration g, which is given as

$$\begin{array}{}\text{(C7)}& g=g_{\mathrm{0}}{\mathrm{\Gamma }}^{-\mathrm{2}}.\end{array}$$

Although not applied in the present paper, the optional time dependence of the generalized curvilinear coordinates enters through the mesh velocity v^g, as indicated in Sect. 2.1.1; see Prusa and Smolarkiewicz (2003) and Kühnlein et al. (2012) for further discussion.

Appendix D: Summary of variables and physical constants

Table D1 List of variables.

Table D2 List of physical constants.