**Model description paper**
13 Feb 2019

**Model description paper** | 13 Feb 2019

# FVM 1.0: a nonhydrostatic finite-volume dynamical core for the IFS

Christian Kühnlein Willem Deconinck Rupert Klein Sylvie Malardel Zbigniew P. Piotrowski Piotr K. Smolarkiewicz Joanna Szmelter and Nils P. Wedi

^{1},

^{1},

^{2},

^{1,3},

^{4},

^{1},

^{5},

^{1}

**Christian Kühnlein et al.**Christian Kühnlein Willem Deconinck Rupert Klein Sylvie Malardel Zbigniew P. Piotrowski Piotr K. Smolarkiewicz Joanna Szmelter and Nils P. Wedi

^{1},

^{1},

^{2},

^{1,3},

^{4},

^{1},

^{5},

^{1}

^{1}European Centre for Medium-Range Weather Forecasts, Reading, UK^{2}FB Mathematik und Informatik, Freie Universität Berlin, Berlin, Germany^{3}Laboratoire de l'Atmosphére et des Cyclones, Météo-France, La Reunion, France^{4}Institute of Meteorology and Water Management – National Research Institute, Warsaw, Poland^{5}Loughborough University, Leicestershire, LE11 3TU, UK

^{1}European Centre for Medium-Range Weather Forecasts, Reading, UK^{2}FB Mathematik und Informatik, Freie Universität Berlin, Berlin, Germany^{3}Laboratoire de l'Atmosphére et des Cyclones, Météo-France, La Reunion, France^{4}Institute of Meteorology and Water Management – National Research Institute, Warsaw, Poland^{5}Loughborough University, Leicestershire, LE11 3TU, UK

**Correspondence**: Christian Kühnlein (christian.kuehnlein@ecmwf.int)

**Correspondence**: Christian Kühnlein (christian.kuehnlein@ecmwf.int)

Received: 22 Sep 2018 – Discussion started: 11 Oct 2018 – Revised: 17 Jan 2019 – Accepted: 22 Jan 2019 – Published: 13 Feb 2019

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.

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^{1} and high efficiency
(Ritchie et al., 1995; Temperton et al., 2001). In global models, the ST method typically
uses a spherical harmonics representation in spectral space and (reduced,
i.e. quasi-uniform) Gaussian grids (Hortal and Simmons, 1991; Wedi et al., 2015).

At ECMWF, the first forecast model using the ST method became operational in 1983, and the technique is still successfully applied today with the efficient SISL integration of the hydrostatic primitive equations in the Integrated Forecasting System (IFS) (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^{2}, leaving the associated effects solely to the
explicit non-linear residual. Furthermore, in the constant-coefficient
semi-implicit scheme different (i.e. split) boundary conditions are applied
in the linear operator and the non-linear residual. Although still a research
issue, the constant-coefficient semi-implicit scheme may incur reduced
stability under more complex orography for future high-resolution forecasts
in the nonhydrostatic regime.

At ECMWF, the reign of the ST method with the SISL integrators still continues, but future challenges, especially with respect to HPC, nonhydrostatic modelling, and complex orography, can be foreseen. The IFS-FVM represents an alternative dynamical core formulation that can complement IFS-ST with regard to these issues. However, to make IFS-FVM a useful option for global medium-range weather forecasting at ECMWF, it needs to be shown that the model formulation can provide (at least) comparable solution quality to the established IFS. In particular, a fundamental scientific question is whether a second-order FV method on the co-located meshes employed in IFS-FVM can sustain the accuracy of the ST method of IFS-ST. Another important question concerns the computational efficiency of IFS-FVM. At ECMWF and generally in NWP, tight constraints exist with regard to the runtime of the forecast models on the employed supercomputers. Therefore, we will evaluate the basic efficiency in terms of the time to solution of the current IFS-FVM formulation relative to the operational hydrostatic IFS-ST and its nonhydrostatic extension. In the present paper, these issues are investigated using relevant atmospheric flow benchmarks such as those defined in the context of the Dynamical Core Model Intercomparison Project (DCMIP; Ullrich et al., 2017). The DCMIP-2016 benchmarks involve large-scale hydrostatic and small-scale nonhydrostatic flows on the sphere and also emphasize the interaction of the 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.

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.

## 2.1 Finite-Volume Module of the IFS

IFS-FVM solves the deep-atmosphere^{3}, nonhydrostatic, fully compressible equations with a
generalized height-based terrain-following vertical coordinate. Numerical
integration of the governing equations employs a centred two-time-level
semi-implicit scheme that provides unconditional stability in 3-D with
respect to the fast acoustic and buoyant modes, as well as slower rotational
modes (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

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}}\phantom{\rule{0.125em}{0ex}}{\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

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*_{0} ≡ 10^{5} 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:

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

**) in thermal wind balance, which in terms of the momentum Eq. (1b) is**

*x*
where *φ*_{a}≡*c*_{pd}*π*_{a} and
${\mathit{\theta}}_{\mathit{\rho}\mathrm{a}}={\mathit{\theta}}_{\mathrm{a}}\phantom{\rule{0.125em}{0ex}}(\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{\u0303}}{\mathbf{G}}$, its transpose
${\stackrel{\mathrm{\u0303}}{\mathbf{G}}}^{\mathrm{T}}$, and the contravariant velocity $\mathit{v}=\dot{\mathit{x}}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathbf{G}}}^{\mathrm{T}}\phantom{\rule{-0.125em}{0ex}}\mathit{u}\phantom{\rule{0.125em}{0ex}}+\phantom{\rule{0.125em}{0ex}}{\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

**≡2**

*f***Ω**with

**Ω**the angular velocity vector of the Earth's rotation, and 𝓜

^{′}subsumes the metric forces due to the curvature of the sphere:

Details about the specification of the curvilinear space, as well as explicit
expressions of the Coriolis term $-\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:

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

^{4}; see Table 2 for the respective specifications of Ψ,

**,**

*V**G*. The homogeneous mass continuity Eq. (1a) is a particular case of Eq. (6) and plays a fundamental role for the conservative advective transport of all other scalar variables. Note that because of the mass continuity Eq. (1a), the other scalar conservation laws (Eq. 6) are equivalent to the Lagrangian form $\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.

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:

where

and ${\widehat{\mathrm{\Psi}}}_{\mathit{i}}={\mathcal{A}}_{\mathit{i}}(\stackrel{\mathrm{\u0303}}{\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*^{n} and *t*^{n+1}, 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

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{\u0303}}{\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*^{2}) estimate, which is explained below.

Given the preceding discussion, the semi-implicit solution procedure proceeds
from an atmospheric state at *t*^{n} to a state at *t*^{n+1} as described in
the following. The solution procedure commences with the integration of the
mass continuity Eq. (1a) as

which straightforwardly returns the updated density
${\mathit{\rho}}_{\mathrm{d}\phantom{\rule{0.125em}{0ex}}\mathit{i}}^{n+\mathrm{1}}$. The 𝒪(*δ**t*^{2}) 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*^{n−1} and *t*^{n}; see Appendix A of
Kühnlein et al. (2012) for the procedure accounting for a variable time
step *δ**t*. In addition, the algorithm (10) defines the
advector ${\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}}{)}^{\u27c2}{|}^{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*^{n} 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*^{n+1} – the
scheme (7)–(8) for the thermodynamic
Eq. (1b) and momentum Eq. (1b)
is implemented as

where

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

*ϑ*

^{n+1}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

*θ*

^{0}for

*θ*

^{⋆}is provided as the explicit solution of full potential temperature:

whereas a first guess *u*^{0} for *u*^{⋆} is prescribed by linear
extrapolation of ** u** from the

*t*

^{n−1}and

*t*

^{n}time levels (again taking into account the variable time step

*δ*

*t*).

From Eqs. (11a)–(11b), closed-form
expressions for the discrete velocity update are derived by eliminating
Eq. (11a) for ${\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

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

**on the LHS and**

*u***L**

^{−1}its inverse, Eq. (14) can be symbolized as

and

respectively, where $\stackrel{\mathrm{\u02c7}}{\stackrel{\mathrm{\u02c7}}{\mathit{u}}}={\mathbf{L}}^{-\mathrm{1}}\widehat{\widehat{\mathit{u}}}$, and
$\mathbf{C}={\mathbf{L}}^{-\mathrm{1}}\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}\mathit{\delta}t\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathit{\rho}}^{\star}\stackrel{\mathrm{\u0303}}{\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

where $\mathit{\phi}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}{\mathit{\phi}}_{\mathrm{a}}+{\mathit{\phi}}^{\prime}$ and

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

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*

^{n+1}, which can be written compactly as

where the spatial grid index ** i** has been omitted again. The summation
ℓ in Eq. (20) is over the three divergence operators on the
RHS of Eq. (17). The symbolic notations ℒ(

*φ*

^{′}) 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

and $\widehat{{\mathit{\phi}}^{\prime}}=\mathcal{A}({\mathit{\phi}}^{\prime n}+(\mathrm{1}-\mathit{\alpha})\phantom{\rule{0.125em}{0ex}}\mathit{\delta}t\phantom{\rule{0.125em}{0ex}}{\mathcal{R}}^{{\mathit{\phi}}^{\prime}}{|}^{n}+\mathit{\delta}t{P}^{{\mathit{\phi}}^{\prime}}{|}^{n},(\mathit{v}\mathcal{G}{\mathit{\rho}}_{\mathrm{d}}{)}^{\u27c2}{|}^{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{\u0303}}{\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*^{2})
integration of Eq. (17) with the first-order backward Euler
scheme *α*=1 maintains second-order accuracy of all variables except
*φ*^{′} over a single time step because *φ*^{′} enters the pressure
gradient term of the momentum equation with the factor 0.5*δ**t*, hence
resulting in an *O*(*δ**t*^{3}) integral. The accumulation of first-order
errors in *φ*^{′} can be mitigated by solving Eq. (17) with
the weakly off-centred trapezoidal scheme using *α*=0.51; see
Benacchio et al. (2014) for alternative design and pertinent discussion. All
IFS-FVM results presented in Sect. 3 were obtained with the
backward Euler scheme *α*=1, as this has been the default in the
earlier implementations. However, using the extended-range forecast of the
baroclinic instability benchmark, we will demonstrate that either *α*=1 or *α*=0.51 shows the same close agreement with IFS-ST. Moreover,
both choices *α*=1 or *α*=0.51 provide essentially identical
computational performance.

Overall, the 3-D implicit scheme with respect to *φ*^{′} permits time
steps equivalent to soundproof models
(Kurowski et al., 2014; Smolarkiewicz et al., 2014; 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{\u0303}}{\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^{5} 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*.

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

where *l*(*i*) numbers the edges connecting node (*k*,*i*) with its horizontal
neighbours (*k*,*j*), and *S*_{j} refers both to the face per se and its surface
area^{6}. The geometric
quantities *S*_{j}, 𝒱_{i}, and *δ**z* of the computational space
are independent of height. The fields ${A}_{k,j}^{\u27c2}$ 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}^{\u27c2}=\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}{\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}^{\u27c2}=\mathrm{0.5}\phantom{\rule{0.125em}{0ex}}{\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

where ${S}_{j}^{x}$ and ${S}_{j}^{y}$ denote the *x* and *y* components of the oriented
surface element **S**_{j}=*S*_{j}*n*_{j}. Given Eq. (22) and
Eq. (23) in the computational space, they are augmented
with metrics of the curvilinear coordinate framework to obtain the
respective physical divergence and gradient appearing in the governing
equations of Sect. 2.1.1 and
2.1.2, respectively. For illustration, one example for the
revised implementation of the velocity divergence in the first term on the
RHS of Eq. (17) at the node (*k*,*i*) is

The details about the specification of the curvilinear coordinate framework under shallow- and deep-atmosphere equations are described in Appendix C.

For IFS-FVM, the mesh generation and mesh data structures, as well as the nearest-neighbour 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^{7} and
buoyant modes, and the 3-D SL advection
(Ritchie et al., 1995; Temperton et al., 2001; Hortal, 2002). Therefore, the constant time
step of the SISL integration in IFS-ST can be selected according to optimal
efficacy rather than stability. The ST method, which is applied along model
levels in the horizontal, is combined with a finite-element (FE) approach to
discretize the integral operator in the vertical direction
(Untch and Hortal, 2004)^{8}. In 2002, this vertical FE scheme of Untch and Hortal (2004)
with a co-located arrangement of prognostic variables replaced the former FD
scheme of Simmons and Burridge (1981) 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

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.

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*^{n} and arrival
point at *t*^{n+1}. Thereby, the tendencies at *t*^{n+1} use a provisional
first guess from the explicit dynamics as input, and the final tendencies are
applied at the arrival point of the SL trajectory; see Wedi (1999) for
details^{9}. The basic approach of the physics–dynamics
coupling in the IFS uses sequential splitting of tendencies, i.e. the various
processes are integrated one after another, and the updated tendencies are
used as input to the subsequent process; see Beljaars et al. (2018)
for discussion. More details about the IFS physics parameterizations can be
found in the general IFS documentation.

The physics–dynamics coupling in IFS-FVM differs from IFS-ST. As explained in
Sect. 2.1.2, the current implementation of IFS-FVM incorporates
the tendencies from physics parameterizations *P*^{Ψ} by means of a
first-order coupling at *t*^{n}. Therefore, the fields that enter the
parameterizations are from *t*^{n}, and there is no averaging between tendencies
from *t*^{n} and *t*^{n+1}. While the incorporation of the tendencies from
physics parameterization deviates from IFS-ST, the sequential splitting
between the various IFS physics parameterizations is kept exactly the same.
Incorporating the physics parameterizations with first order at *t*^{n} is
motivated by the generally smaller time steps *δ**t* in IFS-FVM than in
IFS-ST and the desire to implement straightforward options for subcycling of
the dynamics (see below). In addition, numerical experimentation with IFS-FVM
so far has shown favourable results in terms of the incorporation of the
physics at *t*^{n}. Nevertheless, different forms of coupling IFS-FVM to the
physics parameterizations will be explored in the future.

The IFS-FVM code has its own interface to the IFS physics parameterizations.
Among others, it involves conversion between IFS-FVM's variables and those
employed in IFS-ST (see Table 1)^{10}, interpolations to vertical interfaces, and the provision of local
quantities describing the mesh geometry, but also a number of technical
aspects due to some differences in the computational design of IFS-FVM and
IFS-ST. Some of these operations in the interface may be removed at a later
stage when the IFS physics parameterization package becomes more harmonized
with IFS-FVM. However, generally the coupling is facilitated by common
features of IFS-FVM and IFS-ST, such as 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}}:$

where

The physics tendency *P*^{Ψ} is evaluated with the physics time step
Δ*t*_{phys} and is then reused for the *N*_{s} subcycling
steps with *δ**t*. 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^{11}, the current
operational IFS-ST evaluates the physics parameterizations at every
semi-implicit time step *δ**t*. However, due to the Eulerian versus SL
advection, IFS-FVM uses considerably smaller time steps than
IFS-ST^{12}, and
therefore physics–dynamics coupling may use some form of subcycling as
indicated above or other approaches, such as parallel splitting, to remain
competitive. The cost per time step of the full IFS physics parameterization
package can be up to 40 % of the forecast model. For the idealized DCMIP
experiments considered in the present work, we focus on the parameterization
of clouds and stratiform precipitation incorporated by means of the general
IFS-FVM and IFS-ST interfaces described above.

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

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^{13}. For the comparison of computational efficiency in
Sect. 3.3, we employ the (O1280,TCo1279) horizontal grid,
corresponding to about 9 km spacing, with 137 stretched vertical levels,
again with a similar distribution in IFS-FVM and IFS-ST. With regard to the
time step size, IFS-ST uses constant increments *δ**t* of 1800, 1200, and
450 s for the TCo159, TCo319, and TCo1279 grids, respectively. In contrast,
IFS-FVM generally applies variable time stepping that targets a maximum
horizontal advective Courant number of 0.95 – the actual maximum 3-D Courant
number can be significantly larger than 1 as this is permitted by the
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.

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.

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

Kinetic energy spectra evaluated at day 15 in Fig. 7 reveal a
strikingly similar distribution of variance across wavenumbers from 1 to
∼200. As all the simulations are without any explicit diffusion, the
IFS-FVM spectra at the high wavenumbers attest to the implicit
scale-selective regularization with artificial viscosity provided by the
non-oscillatory finite-volume MPDATA advection^{14}. The slope of the
IFS-ST and IFS-FVM spectra with respect to wavenumber *l* is somewhat
shallower than *l*^{−3} at large scales. This is consistent with results from
other dynamical cores for this baroclinic instability test case studied in
the context of the High-Impact Weather Prediction Project
(HIWPP; Whitaker, 2014). The spectra of IFS-ST feature some accumulation of
energy near the scale of the triangular-cubic truncation, corresponding to
4 times the grid spacing. This increased energy at these scales does not
grow and is to a large extent controlled by the spectral filtering of the
non-linear terms at every time step, among other mechanisms of implicit
dissipation such as the semi-Lagrangian interpolation.

## 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^{15} for the
(O160,TCo159) and (O320,TCo319) grids at day 10. For any of these grids, both
model formulations show five rainbands with essentially identical phase, as
emphasized by the overlay with the 0.5 mm day^{−1} black contour line of the corresponding other model formulation.
The elongated rainbands are associated with the lifting along sharp frontal
zones. Precipitation amounts are overall similar but somewhat higher local
values exist for IFS-FVM, particularly in the two easternmost rainbands when
looking at the (O160,TCo159) grid. Figure 9 is analogous
to Fig. 8 but for day 15. As can be expected, the spread
between the different model formulations becomes larger. However, there is
still reasonably close agreement, especially for the higher-resolution grid
(O320,TCo319) in the right column of Fig. 9. Here, the
locations of the easternmost frontal zone and associated rainband agree
closely considering the late stage of the baroclinic instability evolution.
Figure 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.

## 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^{16}. Importantly, while IFS-ST used the constant time
step of 450 s, IFS-FVM employed variable time stepping according to the
maximum permitted advective Courant number; therefore, in order to obtain
realistic numbers for IFS-FVM, the timings were evaluated between day 10 and
15 in the fully non-linear stage of the baroclinic instability evolution.

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.

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

Model codes developed at ECMWF are the intellectual property of ECMWF and its member states, and therefore the IFS code is not publicly available. Access to a reduced version of the IFS code may be obtained from ECMWF under an OpenIFS licence (see http://www.ecmwf.int/en/research/projects/openifs for further information, last access: 6 February 2019).

The model output data can be downloaded from https://doi.org/10.5281/zenodo.1445597 (Kühnlein et al., 2018).

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

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}}^{\u27c2}\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{\u0303}}{\mathrm{\Psi}}$
(Eq. 8) as

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^{17}. 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).

The bespoke preconditioner solves for the solution error *e* of the pressure
perturbation variable *φ*^{′}:

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

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

where 𝒟 is the diagonal coefficient of 𝒫_{h},
specified as

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

with the weight *ω*=0.7.

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 d*x*_{c}=*r*cos *ϕ* d*λ*,
d*y*_{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{\u0303}}{x}=a\mathit{\lambda}$, $\stackrel{\mathrm{\u0303}}{y}=a\mathit{\varphi}$, $\stackrel{\mathrm{\u0303}}{z}={z}_{\mathrm{c}}$
($\stackrel{\mathrm{\u0303}}{x},\stackrel{\mathrm{\u0303}}{y},\stackrel{\mathrm{\u0303}}{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

where $\mathcal{F}(\stackrel{\mathrm{\u0303}}{t},\stackrel{\mathrm{\u0303}}{\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{\u0303}}{x}$, $y\equiv \stackrel{\mathrm{\u0303}}{y}$, combined with a height-based terrain-following vertical coordinate of the general form $z=z(\stackrel{\mathrm{\u0303}}{x},\stackrel{\mathrm{\u0303}}{y},\stackrel{\mathrm{\u0303}}{z})$. The most straightforward specification of the mapping $z=z(\stackrel{\mathrm{\u0303}}{x},\stackrel{\mathrm{\u0303}}{y},\stackrel{\mathrm{\u0303}}{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{\u0303}}{x},\stackrel{\mathrm{\u0303}}{y},\stackrel{\mathrm{\u0303}}{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{\u0303}}{\mathbf{G}}$ employed in the formalism of Sect. 2.1.1 is given as

where

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

The inverse metrics $\partial z/\partial \stackrel{\mathrm{\u0303}}{x}$, $\partial z/\partial \stackrel{\mathrm{\u0303}}{y}$ and $\partial z/\partial \stackrel{\mathrm{\u0303}}{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

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

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

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.

CK did most of the developments and numerical experiments presented in the paper. CK and PKS are the main developers of IFS-FVM. WD is the main developer of the Atlas library employed by IFS-FVM. RK contributed to developments of the advection scheme and time stepping. SM set IFS-ST up for the test cases and performed experiments for the comparison to IFS-FVM. ZPP contributed to the efficient coding of IFS-FVM. JS contributed to the preconditioning of the elliptic solver. NPW provided support with regard to the general IFS framework.

The authors declare that they have no conflict of interest.

We would like to thank two reviewers for their useful comments. Helpful
discussions with the physics parameterization team at ECMWF are gratefully
acknowledged. This work was supported in part by funding received from the
European Research Council under the European Union's Seventh Framework
Programme FP7/2012/ERC (grant agreement no. 320375), the Horizon 2020
Research and Innovation Programme ESCAPE (grant agreement no. 671627), and
ESIWACE (grant agreement no. 675191). Rupert Klein thanks ECMWF for support
under their fellowship programme and acknowledges partial support of his
contributions by the Deutsche Forschungsgemeinschaft, grant CRC 1114 “Scaling
Cascades in Complex Systems”, project A02. Zbigniew Piotrowski acknowledges
partial support from the “Numerical weather prediction for sustainable
Europe” project, carried out within the FIRST TEAM programme of the
Foundation for Polish Science co-financed by the European Union under the
European Regional Development Fund.

Edited by:
Paul Ullrich

Reviewed by: two anonymous referees

Bacon, D. P., Ahmad, N. N., Boybeyi, Z., Dunn, T. J., Hall, M. S., Lee, P. C. S., Sarma, R. A., Turner, M. D., Waight, K. T., Young, S. H., and Zack, J. W.: A dynamically adapting weather and dispersion model: the operational multiscale environment model with grid adaptivity (OMEGA), Mon. Weather Rev., 128, 2044–2076, 2000. a

Bauer, P., Thorpe, A., and Brunet, G.: The quiet revolution of numerical weather prediction, Nature, 525, 47–55, 2015. a

Beljaars, A., Balsamo, G., Bechthold, P., Bozzo, A., Forbes, R., Hogan, R. J., Köhler, M., Morcrette, J. J., Tompkins, A., Viterbo, P., and Wedi, N. P.: The numerics of physical parameterization in the ECMWF model, Front. Earth Sci., 6, 1–18, https://doi.org/10.3389/feart.2018.00137, 2018. a

Benacchio, T., O'Neill, W. P., and Klein, R.: A blended soundproof-to-compressible numerical model for small- to mesoscale atmospheric dynamics, Mon. Weather Rev., 142, 4416–4438, 2014. a

Bénard, P., Vivoda, J., Mašek, J., Smolíková, P., Yessad, K., Smith, C., Brožková, R., and Geleyn, J.-F.: Dynamical kernel of the Aladin-NH spectral limited-area model: Revised formulation and sensitivity experiments, Q. J. Roy. Meteor. Soc., 136, 155–169, 2010. a, b, c, d, e

Bubnová, R., Hello, G., Bénard, P., and Geleyn, J.-F.: Integration of the fully elastic equations cast in the hydrostatic pressure terrain-following coordinate in the framework of the ARPEGE/Aladin NWP system, Mon. Weather Rev., 123, 515–535, 1995. a, b, c

Cossette, J.-F., Smolarkiewicz, P. K., and Charbonneau, P.: The Monge-Ampère trajectory correction for semi-Lagrangian schemes, J. Comput. Phys., 274, 208–229, 2014. a

Deconinck, W., Bauer, P., Diamantakis, M., Hamrud, M., Kühnlein, C., Maciel, P., Mengaldo, G., Quintino, T., Raoult, B., Smolarkiewicz, P. K., and Wedi, N. P.: Atlas: A library for numerical weather prediction and climate modelling, Comput. Phys. Commun., 220, 188–204, 2017. a, b

Diamantakis, M. and Augusti-Panareda, A.: A positive definite tracer mass fixer for high resolution weather and atmospheric composition forecasts, ECMWF Technical Memorandum, 819, available at: https://www.ecmwf.int/en/elibrary/17914-positive-definite-tracer-mass-fixer-high-resolution-weather-and-atmospheric (last access: 7 February 2019), 2017. a

Diamantakis, M. and Flemming, J.: Global mass fixer algorithms for conservative tracer transport in the ECMWF model, Geosci. Model Dev., 7, 965–979, https://doi.org/10.5194/gmd-7-965-2014, 2014. a, b

Diamantakis, M. and Magnusson, L.: Sensitivity of the ECMWF model to semi-Lagrangian departure point iterations, Mon. Weather Rev., 144, 3233–3250, 2016. a

Domaradzki, J. A., Xiao, Z., and Smolarkiewicz, P. K.: Effective eddy viscosities in implicit large eddy simulations of turbulent flows, Phys. Fluids, 15, 3890–3893, 2003. a

Dörnbrack, A., Doyle, J. D., Lane, T. P., Sharman, R. D., and Smolarkiewicz, P. K.: On physical realizability and uncertainty of numerical solutions, Atmos. Sci. Lett., 6, 118–122, 2005. a

Eliasen, E., Machenbauer, B., and Rasmussen, E.: On a numerical method for integration of the hydrodynamical equations with a spectral representation of the horizontal fields, Tech. Rep. 2, Institute of Theoretical Meteorology, University of Copenhagen, 1970. a

Emanuel, K. A.: Atmospheric convection, Oxford University Press on Demand, 580 pp., New York, Oxford, 1994. a

Forbes, R. M., Tompkins, A. M., and Untch, A.: A new prognostic bulk microphysics scheme for the IFS, ECMWF Technical Memorandum, 649, available at: https://www.ecmwf.int/en/elibrary/9441-new-prognostic-bulk-microphysics-scheme-ifs (last access: 7 February 2019) 2010. a, b

Gal-Chen, T. and Somerville, R. C. J.: On the use of a coordinate transformation for the solution of the Navier-Stokes equations, J. Comput. Phys., 17, 209–228, 1975. a

Hortal, M.: The development and testing of a new two-time-level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model, Q. J. Roy. Meteor. Soc., 128, 1671–1687, 2002. a, b

Hortal, M. and Simmons, A.: Use of reduced Gaussian grids in spectral models, Mon. Weather Rev., 119, 1057–1074, 1991. a, b, c

Klemp, J. B.: A terrain-following coordinate with smoothed coordinate surfaces, Mon. Weather Rev., 139, 2163–2169, 2011. a

Klemp, J. B., Dudhia, J., and Hassiotis, A. D.: An upper gravity-wave absorbing layer for NWP applications, Mon. Weather Rev., 136, 3987–4004, 2008. a

Knoll, D. A., Chacon, L., Margolin, L. G., and Mousseau, V. A.: On balanced approximations for time integration of multiple time scale systems, J. Comput. Phys., 185, 583–611, 2003. a

Kühnlein, C. and Smolarkiewicz, P. K.: An unstructured-mesh finite-volume MPDATA for compressible atmospheric dynamics, J. Comput. Phys., 334, 16–30, 2017. a, b, c, d, e, f, g, h, i, j, k, l

Kühnlein, C. and Smolarkiewicz, P. K.: A nonhydrostatic finite-volume option for the IFS, ECMWF Newsletter, 158, 30–36, https://doi.org/10.21957/sd92ack6p3, 2019. a, b

Kühnlein, C., Smolarkiewicz, P. K., and Dörnbrack, A.: Modelling atmospheric flows with adaptive moving meshes, J. Comput. Phys., 231, 2741–2763, 2012. a, b, c, d, e, f, g

Kühnlein, C., Deconinck, W., Klein, R., Malardel, S., Piotrowski, Z., Smolarkiewicz, P., Szmelter, J., and Wedi, N.: FVM 1.0: A nonhydrostatic finite-volume dynamical core formulation for IFS, Zenodo, https://doi.org/10.5281/zenodo.1445597, 2018. a

Kurowski, M. J., Grabowski, W. W., and Smolarkiewicz, P. K.: Anelastic and compressible simulation of moist deep convection, J. Atmos. Sci., 71, 3767–3787, 2014. a, b

Malardel, S. and Ricard, D.: An alternative cell-averaged departure point reconstruction for pointwise semi-Lagrangian transport schemes, Q. J. Roy. Meteor. Soc., 141, 2114–2126, 2015. a

Malardel, S., Wedi, N. P., Deconinck, W., Diamantakis, M., Kühnlein, C., Mozdzynski, G., Hamrud, M., and Smolarkiewicz, P. K.: A new grid for the IFS, ECMWF Newsletter, 146, 23–28, 2016. a, b, c, d

Orszag, S. A.: Transform method for the calculation of vector-coupled sums: application to the spectral form of the vorticity equation, J. Atmos. Sci., 27, 890–895, 1970. a

Piotrowski, Z. P., Smolarkiewicz, P. K., Malinowski, S. P., and Wyszogrodzki, A. A.: On numerical realizability of thermal convection, J. Comput. Phys., 228, 6268–6290, 2009. a

Prusa, J. M. and Smolarkiewicz, P. K.: An all-scale anelastic model for geophysical flows: dynamic grid deformation, J. Comput. Phys., 190, 601–622, 2003. a, b, c, d, e, f, g

Prusa, J. M., Smolarkiewicz, P. K., and Wyszogrodzki, A. A.: EULAG, a computational model for multiscale flows, Comput. Fluids, 37, 1193–1207, 2008. a, b

Ritchie, H., Temperton, C., Simmons, A., Hortal, M., Davies, T., Dent, D., and Hamrud, M.: Implementation of the semi-Lagrangian method in a high-resolution version of the ECMWF forecast model, Mon. Weather Rev., 123, 489–514, 1995. a, b, c

Robert, A., Henderson, J., and Turnbull, C.: An implicit time integration scheme for baroclinic models of the atmosphere, Mon. Weather Rev., 100, 329–335, 1972. a

Schär, C., Leuenberger, D., Fuhrer, O., Lüthi, D., and Girard, C.: A new terrain-following vertical coordinate formulation for atmospheric prediction models, Mon. Weather Rev., 130, 2459–2480, 2002. a

Simmons, A. J. and Burridge, D. M.: An energy and angular-momentum conserving vertical finite-difference scheme and hybrid vertical coordinates, Mon. Weather Rev., 109, 758–766, 1981. a, b

Smolarkiewicz, P. and Szmelter, J.: A nonhydrostatic unstructured-mesh soundproof model for simulation of internal gravity waves, Acta Geophys., 59, 1109–1134, 2011. a

Smolarkiewicz, P. K. and Dörnbrack, A.: Conservative integrals of adiabatic Durran's equations, Int. J. Numer. Methods Fluids, 56, 1513–1519, 2008. a

Smolarkiewicz, P. K. and Grabowski, W.: The multidimensional positive definite advection transport algorithm: nonoscillatory option, J. Comput. Phys., 86, 355–375, 1990. a

Smolarkiewicz, P. K. and Margolin, L. G.: On forward-in-time differencing for fluids – Extension to a curvilinear framework, Mon. Weather Rev., 121, 1847–1859, 1993. a

Smolarkiewicz, P. K. and Margolin, L.: Variational methods for elliptic problems in fluid models, in: ECMWF Proceedings, Workshop on developments in numerical methods for very high resolution global models, Reading, UK, 137–159, 2000. a, b

Smolarkiewicz, P. K. and Pudykiewicz, J. A.: A class of semi-Lagrangian approximations for fluids, J. Atmos. Sci., 49, 2082–2096, 1992. a

Smolarkiewicz, P. K. and Szmelter, J.: MPDATA: an edge-based unstructured-grid formulation, J. Comput. Phys., 206, 624–649, 2005. a

Smolarkiewicz, P. K. and Szmelter, J.: Iterated upwind schemes for gas dynamics, J. Comput. Phys., 228, 33–54, 2009. a

Smolarkiewicz, P. K., Temperton, C., Thomas, S. J., and Wyszogrodzki, A. A.: Spectral Preconditioners for nonhydrostatic atmospheric models: extreme applications, in: Proceedings of the ECMWF seminar series on recent developments in numerical methods for atmospheric and ocean modelling, Reading, UK, 203–220, 2004. a

Smolarkiewicz, P. K., Szmelter, J., and Wyszogrodzki, A. A.: An unstructured-mesh atmospheric model for nonhydrostatic dynamics, J. Comput. Phys., 254, 184–199, 2013. a

Smolarkiewicz, P. K., Kühnlein, C., and Wedi, N. P.: A consistent framework for discrete integrations of soundproof and compressible PDEs of atmospheric dynamics, J. Comput. Phys., 263, 185–205, 2014. a, b, c, d, e, f, g, h, i, j

Smolarkiewicz, P. K., Deconinck, W., Hamrud, M., Kühnlein, C., Mozdzynski, G., Szmelter, J., and Wedi, N. P.: A finite-volume module for simulating global all-scale atmospheric flows, J. Comput. Phys., 314, 287–304, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Smolarkiewicz, P. K., Kühnlein, C., and Grabowski, W.: A finite-volume module for cloud-resolving simulations of global atmospheric flows, J. Comput. Phys., 341, 208–229, 2017. a, b, c, d, e, f, g

Smolarkiewicz, P. K., Kühnlein, C., and Wedi, N. P.: Semi-implicit integrations of perturbation equations for all-scale atmospheric dynamics, J. Comput. Phys., 376, 145–159, 2019. a, b

Szmelter, J. and Smolarkiewicz, P. K.: An edge-based unstructured mesh discretisation in geospherical framework, J. Comput. Phys., 229, 4980–4995, 2010. a, b, c, d, e, f, g

Temperton, C.: Treatment of the Coriolis terms in semi-Lagrangian spectral models, Atmos. Ocean, 35, 293–302, 2011. a

Temperton, C., Hortal, M., and Simmons, A.: A two-time-level semi-Lagrangian global spectral model, Q. J. Roy. Meteor. Soc., 127, 111–127, 2001. a, b

Thuburn, J.: Some conservation issues for the dynamical cores of NWP and climate models, J. Comput. Phys., 227, 3715–3730, 2008. a

Ullrich, P. A., Melvin, T., Jablonowski, C., and Staniforth, A.: A proposed baroclinic wave test case for deep- and shallow-atmosphere dynamical cores, Q. J. Roy. Meteor. Soc., 140, 1590–1602, 2014. a, b

Ullrich, P. A., Jablonowski, C., Reed, K. A., Zarzycki, C., Lauritzen, P. H., Nair, R. D., Kent, J., and Verlet-Banide, A.: Dynamical Core Model Intercomparison Project (DCMIP2016) Test Case Document, available at: https://github.com/ClimateGlobalChange/DCMIP2016 (last access: 6 February 2019), 2016. a, b, c, d

Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Zarzycki, C. M., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Dubos, T., Meurdesoif, Y., Chen, X., Harris, L., Kühnlein, C., Lee, V., Qaddouri, A., Girard, C., Giorgetta, M., Reinert, D., Klemp, J., Park, S.-H., Skamarock, W., Miura, H., Ohno, T., Yoshida, R., Walko, R., Reinecke, A., and Viner, K.: DCMIP2016: a review of non-hydrostatic dynamical core design and intercomparison of participating models, Geosci. Model Dev., 10, 4477–4509, https://doi.org/10.5194/gmd-10-4477-2017, 2017. a

Untch, A. and Hortal, M.: A finite-element scheme for the vertical discretization of the semi-Lagrangian version of the ECMWF forecast model, Q. J. Roy. Meteor. Soc., 130, 1505–1530, 2004. a, b

Waruszewski, M., Kühnlein, C., Pawlowska, H., and Smolarkiewicz, P. K.: MPDATA: Third-order accuracy for variable flows, J. Comput. Phys., 359, 361–379, 2018. a

Wedi, N. and Düben, P.: Extreme scaling for global weather forecasts at O(1km) horizontal resolution, Geophys. Res. Abstr., EGU2017-8671, EGU General Assembly 2017, Vienna, Austria, 2017. a

Wedi, N. P.: The numerical coupling of the physical parametrizations to the dynamical equations in a forecast model, ECMWF Technical Memorandum, 274, available at: https://www.weather.gov/media/sti/nggps/HIWPP_idealized_tests-v8 revised 05212015.pdf (last access: 7 February 2019) 1999. a, b, c

Wedi, N. P.: Increasing horizontal resolution in numerical weather prediction and climate simulations: illusion or panacea?, Philos. T. R. Soc. A, 372, 20130289, https://doi.org/10.1098/rsta.2013.0289, 2014. a, b, c

Wedi, N. P. and Smolarkiewicz, P. K.: Direct numerical simulation of the Plumb McEwan laboratory analog of the QBO, J. Atmos. Sci., 63, 3226–3252, 2006. a

Wedi, N. P., Hamrud, M., and Mozdzynski, G.: A fast spherical harmonics transform for global NWP and climate models, Mon. Weather Rev., 141, 3450–3461, 2013. a, b, c, d

Wedi, N. P., Bauer, P., Deconinck, W., Diamantakis, M., Hamrud, M., Kühnlein, C., Malardel, S., Mogensen, K., Mozdzynski, G., and Smolarkiewicz, P. K.: The modelling infrastructure of the Integrated Forecasting System: Recent advances and future challenges, ECMWF Technical Memorandum, 760, 1–48, 2015. a, b, c, d, e, f, g, h, i, j

Weller, H., Ringler, T., Piggott, M., and Wood, N.: Challenges facing adaptive mesh modeling of the atmosphere and ocean, B. Am. Meteorol. Soc., 91, 105–108, 2010. a

Whitaker, J.: HIWPP non-hydrostatic dynamical core tests: Results from idealized test cases, Tech. rep., available at: https://www.ecmwf.int/en/elibrary/13020-numerical-coupling-physical-parametrizations-dynamical-equations-forecast-model (last access: 7 February 2019), 2014. a

Williamson, D. L.: The evolution of dynamical cores for global atmospheric models, J. Meteorol. Soc. Jpn., 85, 241–269, 2007. a

Wood, N., Staniforth, A., White, A., Allen, T., Diamantakis, M., Gross, M., Melvin, T., Smith, C., Vosper, S., Zerroukat, M., and Thuburn, J.: An inherently mass-conserving semi-implicit semi-Lagrangian discretization of the deep-atmosphere global non-hydrostatic equations, Q. J. Roy. Meteor. Soc., 140, 1505–1520, 2014. a

Zarzycki, C. M., Jablonowski, C., and Taylor, M. A.: Using variable-resolution meshes to model tropical cyclones in the Community Atmosphere Model, Mon. Weather Rev., 142, 1221–1239, 2014. a

Zarzycki, C. M., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Ullrich, P. A., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Chen, X., Harris, L., Giorgetta, M., Reinert, D., Kühnlein, C., Walko, R., Lee, V., Qaddouri, A., Tanguay, M., Miura, H., Ohno, T., Yoshida, R., Park, S.-H., Klemp, J., and Skamarock, W.: DCMIP2016: The Splitting Supercell Test Case, Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd-2018-156, in review, 2018. a

^{1}

The SL schemes are subject to a topological realizability condition based on the Lipschitz number which is related to the flow deformation (Smolarkiewicz and Pudykiewicz, 1992; Cossette et al., 2014). However, in NWP this condition is typically much less restrictive than the advective CFL stability condition of Eulerian schemes; see e.g. Diamantakis and Magnusson (2016).

^{2}

Including orography in the implicit part involves multiplications which are standardly performed in grid-point space in the context of the ST method. In principle, one could carry out the necessary multiplications in spectral space but this is usually avoided because of computational complexity.

^{3}

The shallow-atmosphere equations,
the default in IFS-ST, are available by means of a simple switch *γ*
(Appendix C).

^{4}

As an example, ${\mathcal{R}}^{{\mathit{\theta}}^{\prime}}\equiv -{\stackrel{\mathrm{\u0303}}{\mathbf{G}}}^{\mathrm{T}}\mathit{u}\cdot \mathrm{\nabla}{\mathit{\theta}}_{\mathrm{a}}$ when considering Eq. (1b).

^{5}

Note that although the presented IFS-FVM formulation assumes an unstructured mesh with indirect addressing in the horizontal, the model may exploit structured or semi-structured grids on future HPC architectures.

^{6}

Note that in IFS-FVM, *S*_{j} and 𝒱_{i} have
dimensions of length and area, and the actual face areas and volumes of
prismatic cells are *S*_{j}*δ**z* and 𝒱_{i}*δ**z*, respectively, in computational space (Smolarkiewicz et al., 2016).

^{7}

The HPEs analytically filter internal acoustic modes but support the external Lamb mode.

^{8}

The FE implementation of the discrete vertical integral operator is based on the Galerkin method using cubic B splines as basis functions.

^{9}

SLAVEPP is applied to tendencies from radiation, moist
convection, and the cloud scheme, whereas tendencies from turbulence and
gravity wave drag parameterizations are incorporated with first order at
*t*^{n+1} (Wedi, 1999).

^{10}

Examples are
conversions between mixing ratios *r*_{k} and specific water content variables
*q*_{k} or between quantities in the height- versus pressure-based coordinate
systems.

^{11}

In the current
high-resolution deterministic IFS forecasts on the O1280 grid, the radiation
scheme is called every hour, compared to the semi-implicit model time step
*δ**t*=450 s, and is also run on the coarser O400 grid.

^{12}

With the current formulations, the time step *δ**t* in
IFS-FVM is typically about a factor of 6–7 smaller than in IFS-ST.

^{13}

For instance, the 1 m height difference corresponds to about 0.1 hPa near the surface, which is negligible with regard to the subsequent analysis.

^{14}

The unsplit NFT MPDATA advection (Kühnlein and Smolarkiewicz, 2017) features lower implicit diffusion near the grid scale than the split advection applied here.

^{15}

The precipitation rate represents the liquid and rain (excluding ice and snow) sedimentation flux at the surface.

^{16}

Each node on this supercomputer consists of two Intel Xeon EP E5-2695 v4 “Broadwell” processors, each with 18 cores, which for the 350 compute nodes employed results in a total of 12 600 cores. Here, a hybrid MPI–OpenMP parallelization with six threads was used by all three dynamical cores.

^{17}

Compared to the unsplit scheme, the particular horizontal–vertical splitting also does not incur any additional parallel communication in the context of the horizontal domain decomposition of IFS-FVM.

- Abstract
- Introduction
- IFS model formulations
- Benchmark simulation results
- Conclusions
- Code availability
- Data availability
- Appendix A: Horizontal–vertical splitting of the NFT advective transport
- Appendix B: Weighted line Jacobi preconditioner
- Appendix C: Geospherical framework, generalized terrain-following vertical coordinate, shallow- and deep-atmosphere equations
- Appendix D: Summary of variables and physical constants
- Author contributions
- Competing interests
- Acknowledgements
- References

- Abstract
- Introduction
- IFS model formulations
- Benchmark simulation results
- Conclusions
- Code availability
- Data availability
- Appendix A: Horizontal–vertical splitting of the NFT advective transport
- Appendix B: Weighted line Jacobi preconditioner
- Appendix C: Geospherical framework, generalized terrain-following vertical coordinate, shallow- and deep-atmosphere equations
- Appendix D: Summary of variables and physical constants
- Author contributions
- Competing interests
- Acknowledgements
- References