TIME AND SPACE DISCRETIZATIONS

The equation sets are a prescription for advancing the meteorological fields in time for one time step. They furthermore determine the choice of the prognostic variables. For example, the thermodynamic energy can be represented by a temperature or potential temperature equation.

Alternatively, some models predict the evolution of moist thermodynamic quantities, like the virtual potential temperature, which is e.g. implemented in the Finite-Volume dynamical core (Lin, 2004) at the National Center for Atmospheric Research (NCAR) or the Geophysical Fluid Dynamics Laboratory (GFDL). Another alternative is a conservation equation for the total energy as documented in Satoh et al. (2008). However, the latter can create ambiguity with respect to the inclusion or non-inclusion of moisture-related contributions to the total energy. Once the equation set has been chosen, it needs to be determined on what kind of horizontal and vertical grids the fields are represented, and how each term in the set will be discretized. The advective terms can be treated from an Eulerian or a Lagrangian point of view either in advective or in flux form. In the Lagrangian point of view the set is supplemented by a trajectory algorithm. The semi-Lagrangian method is reviewed in Staniforth and Côté (1991). The remaining terms are responsible for the atmospheric oscillations and can be treated either implicitly or explicitly to preserve stability of the model with the chosen time step. Note that the discretization errors in time and in space are connected via the fluid equations.

6.3.1 Horizontal discretization

Traditionally, the global spectral-transform method on latitude-longitude (or Gaussian/reduced-Gaussian) computational grids has been the dominant choice for weather and climate GCMs for solving the horizontal components of the equations on the sphere. Newly-developed global models today are built upon quasi-uniform grids that have recently been reviewed by Staniforth and

Thuburn (2012). The quasi-uniform grids avoid the convergence of the meridians, the so-called pole problem, that has plagued latitude-longitude grids for so long. This convergence leads to very small grid spacings in the longitudinal direction, which either necessitate very short time steps or the application of polar filters, mostly in the form of Fourier filters. The latter become increasingly difficult to apply on parallel computing architectures since they require enhanced communication of the data along latitudes and thereby limit the parallel scalability.

Today, there is a general move towards local spatial discretization methods that avoid global communication on parallel hardware architectures. Examples are finite-difference (Qaddouri and Lee, 2011; Wood et al. 2014), finite-element (Kritsikis and Dubos, 2014; Melvin et al. 2014; Cotter and Thuburn, 2014; Staniforth et al. 2013; Cotter and Shipton, 2012), finite volume (Thuburn et al.

2014; Lee, 2013; Ullrich and Jablonowski, 2012b, Szmelter and Smolarkiewicz, 2010), spectral element (Giraldo and Rosmond, 2004; Taylor and Fournier, 2010; Dennis et al. 2012) and

discontinuous Galerkin (Nair et al. 2009; Bao et al. 2015) methods that utilize almost uniform polyhedral grids on the sphere. The expected future computer architectures with fast computations and relatively slow memory movement favour the development of high-order numerical

approximations by making extra local computations almost free. One also has to take into account that the anticipated amount of memory per processing element will be smaller than today. In a finite-element discretization one has the choice of reducing the element size (h-convergence) or increasing the discretization order in each element (p-convergence). At least for smooth problems the global error can then be minimized with a gain in efficiency but without forgetting that numerical spatial errors have to balance the errors that come from all other sources. A good introduction to these modern numerical methods can be found in Lauritzen et al. (2011).

Current dynamical core developments are on cubed-sphere meshes (Adcroft et al. 2004; Ullrich and Jablonowski, 2012b; Dennis et al. 2012), triangular grids (Walko and Avissar, 2008; Lee, 2013; Zängl et al. 2015; Satoh et al. 2014), hexagonal meshes (Lee and McDonald, 2009;

Skamarock et al. 2012; Gassmann, 2013) and the Yin-Yang composite grid system (Qaddouri and Lee, 2011; Qaddouri et al. 2012; Sakamoto et al. 2013). In the longer term one might use mesh-free representations but these approaches are not competitive today (Flyer et al. 2015).

Icosahedral and hexagonal grids are often used as optimized grids to enhance the regularity of the grid spacing even further. This can be accomplished via a spring dynamics algorithm (Iga and Tomita, 2014) or other iterative techniques (Miura and Kimoto, 2005; Heikes et al. 2013). The more uniform grid point distributions lessen potential grid imprinting issues. These arise when the

computational grid becomes visible as a flow feature, which is undesirable from a physical

viewpoint. A variety of such grid imprinting signatures has been analyzed by Lauritzen et al. (2010) and Weller et al. (2012). Experience furthermore shows that high-order numerical methods help reduce grid imprinting (Ullrich and Jablonowski, 2012b).

Concerning the horizontal grid staggering technique, the most dominant choice today is the Arakawa C-grid. It offers beneficial gravity wave propagation and dispersion characteristics, but also contains some computational modes on quasi-uniform grids as for example discussed by Thuburn and Staniforth (2004), Thuburn (2008b), Thuburn et al. (2009), Gassmann (2011b) and Weller (2012). As outlined in Staniforth and Thuburn (2012) the computational modes are wave modes that are supported by the numerical discretization, but have no analog among the waves that are represented by the continuous equations. The computational modes usually appear at or near the grid scale and can be the cause for numerical noise. Weller et al. (2012) found that such computational modes are most easily controlled on a hexagonal grid via a diffusive advection algorithm.

Although spectral transform methods are being predicted to be phased out, the current spectral model at the European Centre for Medium-Range Weather Forecasts (ECMWF), the Integrated Forecasting System (IFS), in its operational form (albeit hydrostatic) is the benchmark to beat, and it is not clear that any of the new developments are ready to replace it. In addition, there is no pole problem in the reduced Gaussian grid, at least for some time to come, and the IFS has also already been extended to a non-hydrostatic framework (Wedi et al. 2009). Recently, there have been some powerful developments in the use of spectral transform methods, and this is likely to extend the lifetime and use of spectral transform methods for perhaps another decade (Wedi et al.

2013; Wedi, 2014). In addition, hybrid approaches are under investigation that can combine multiple (e.g. soundproof and compressible) equation sets and discretization methods (Smolarkiewicz et al. 2014).

6.3.2 Vertical discretization

In the vertical there are generally four choices, which are a height-based (e.g. Gal-Chen and Sommerville, 1975; Schär et al. 2002; Klemp, 2011), pressure-based (e.g. Phillips, 1957; Simmons and Burridge, 1981), hydrostatic mass-based (e.g. Laprise, 1992; Bubnová et al. 1995, Côté et al.

1998; Wood and Staniforth, 2003) or isentropic-based (e.g. Hsu and Arakawa, 1990; Toy and Randall, 2009; Bleck et al. 2010) generalized coordinate with various arrangements of the variables in the vertical direction. All four choices are mostly employed as terrain-following coordinates. However, a pure height coordinate is suitable for models with shaved-cell or cut-cell

approaches as outlined in Section 6.3.5. The pressure-based vertical coordinates are typically selected for hydrostatic models. The other three choices are mostly used in combination with non-hydrostatic equation sets. The exception is the isentropic-based approach, which has been implemented in both hydrostatic and non-hydrostatic dynamical cores. Most new dynamical core developments today are either built upon the height-based or mass-based vertical coordinate, and emphasis is put on a smooth transition between the terrain-following levels and the levels of constant height (Klemp, 2011).

Concerning the vertical staggering of the prognostic variables, Lorenz staggering (co-located variables at full model levels except the vertical velocity is represented at model interfaces, see Lorenz (1960)) has been favoured in the past because conservation is easier to implement with this arrangement. Charney-Phillips staggering (both the vertical velocity and the thermodynamic variable are located at the model interface) is preferred for computational reasons as e.g.

extensively discussed by Thuburn and Woollings (2005), Girard et al. (2014) and Arakawa and Konor (1996). The choice of the vertical staggering seems to be less important if higher-order numerical discretizations in the vertical are used. However, care must still be taken to control computational modes (Staniforth and Wood, 2005).

If centred finite-difference discretizations are used in the vertical direction, they formally degrade to first-order accuracy for the typical non-equidistant (stretched) vertical grids with smaller level spacings near the ground. Traditionally, continuous mappings are used to minimize vertical discretisation errors. Second-order finite-volume techniques (Ullrich and Jablonowski, 2012a, 2012b) or even higher-order finite-element techniques (Untch and Hortal, 2004) have been

explored for the vertical discretization. Such higher-order methods lead to a considerable increase in accuracy and reduce numerical noise. Alternatively, Lin (2004) implemented a floating

Lagrangian vertical coordinate that traps the flow field within the 2D horizontal layers for about 5-10 dynamics time steps. The flow fields are then periodically remapped to a reference grid, which captures the vertical advection. The floating Lagrangian method has gained increased popularity in recent years (Chen et al. 2013; Dubos and Tort, 2014; Penner et al. 2007), and has also been introduced into the Spectral-Element dynamical core at NCAR (Dennis, 2012). Vertically adaptive meshes and vertical nesting (McTaggart-Cowan et al. 2011) are other active areas of research.

6.3.3 Advection

The Lagrangian form for the advection terms allows for point-wise or flux-form semi-Lagrangian (SL) methods. The mass-conservative flux-form SL discretization has recently regained popularity for the advection component of GCMs due to its projected computational efficiency gains

especially in the presence of many tracers (Erath et al. 2012; Wong et al. 2013; Ullrich et al. 2014;

Wood et al. 2014). SL discretizations allow long, efficient time steps with Courant-Friedrichs-Lewy (CFL) numbers greater than one. However, it is likely that future computer architectures will favour Eulerian or SL methods with CFLs not much larger than 1. This will keep the computations local and reduce the parallel communication costs, thereby off-setting the costs of a shorter time step.

6.3.4 Temporal discretization

The minimization of global communication will favor a splitting between slow and fast waves, which can then be treated differently from the viewpoint of the time discretization scheme. A popular split is HEVI that stands for horizontal explicit and vertical implicit. Active research is also underway to explore other IMplicit/EXplicit (IMEX) combinations and non-traditional time integrators such as semi-implicit predictor corrector schemes or exponential time integration methods for stiff problems (Bao et al. 2015; Durran and Blossey, 2012, 2013; Garcia et al. 2014; Giraldo et al. 2013; Ullrich and Jablonowski, 2012a; Weller et al. 2013a; Clancy and Pudykiewicz, 2013a, 2013b). Research also continues into 3D implicit solvers suitable for massively parallel implementation. These typically feature hierarchical, multi-scale designs using multi-grid techniques (Heikes et al. 2013;

Müller and Scheichl, 2014). Looking further ahead “Parallel-in-time” time integrators would exploit untapped parallelism in GCMs, utilize the newest computing architectures more efficiently and thereby reduce the wall-clock execution time of GCMs (Haut and Wingate, 2014; Haut et al. 2014).

New efficient multi-scale time integration algorithms are an active area of research considering the

large impact of the chosen time-step on the overall timeliness of product delivery in NWP (Weinan et al. 2007).

6.3.5 Orography

Orography is most often considered via terrain-following vertical coordinates but this can lead to large discretization errors, instabilities or noise in the presence of steep terrain (Li et al. 2014;

Weller and Shahrokhi, 2014). The trend towards high horizontal resolutions will further steepen the mountain slopes and worsen this problem. More research is needed to assess the pros and cons of alternative approaches and to include them efficiently in global models such as the cut-cell (Good et al. 2014; Steppeler et al. 2011), shaved-cell (Yamazaki and Satomura, 2008) or the combined-grid approach (Yamazaki and Satomura, 2012, 2010). Moreover, the increasing, yet still only partially resolved, orographic gravity wave activity has profound implications on the global circulation. In particular, their non-local effect is increasingly more difficult to describe by columnar physical parameterizations of gravity wave drag as these are designed to switch off with increasing resolution. It can be shown that this problem is a function of increasing resolved topographic slopes. A recent review of orographic gravity-wave drag can be found in Teixeira (2014).

6.4 TREND TOWARDS HIGH-RESOLUTION, VARIABLE-RESOLUTION AND