Adaptive mesh refinement techniques offer a flexible framework for future variable-resolution climate and weather models since they can focus their computational mesh on certain geographical areas or atmospheric events. Adaptive meshes can also be used to coarsen a latitude–longitude grid in polar regions. This allows for the so-called reduced grid setups.
A spherical, block-structured adaptive grid technique is applied to the Lin–Rood finite-volume dynamical core for weather and climate research. This hydrostatic dynamics package is based on a conservative and monotonic finite-volume discretization in flux form with vertically floating Lagrangian layers. The adaptive dynamical core is built upon a flexible latitude–longitude computational grid and tested in two- and three-dimensional model configurations. The discussion is focused on static mesh adaptations and reduced grids. The two-dimensional shallow water setup serves as an ideal testbed and allows the use of shallow water test cases like the advection of a cosine bell, moving vortices, a steady-state flow, the Rossby–Haurwitz wave or cross-polar flows. It is shown that reduced grid configurations are viable candidates for pure advection applications but should be used moderately in nonlinear simulations. In addition, static grid adaptations can be successfully used to resolve three-dimensional baroclinic waves in the storm-track region.
The use of variable grid resolutions has been discussed for regional atmospheric modelling over the past three decades (Fox-Rabinovitz et al. 1997). In particular, nested grids, stretched grids and dynamic adaptive mesh refinement (AMR) techniques have been suggested in the literature. Nested grids are widely used for local weather predictions and regional climate models to downscale a global large-scale simulation. Examples include the limited-area ‘Weather and Research Forecasting Model’ WRF (Skamarock et al. 2005) and the ‘Canadian Regional Climate Model’ CRCM (Caya & Laprise 1999). Typical grid ratios between 2 and 5 are chosen for the resolution jump (Denis et al. 2002). A fixed-size refined grid is then permanently embedded in a coarse-resolution general circulation model, which initializes the fine grid and periodically updates the lateral boundary conditions of the nested area. Most often, different sets of equations, physics approximations and numerical schemes are used in the two domains. Therefore, special care must be taken to minimize the consequent numerical and physical inconsistencies across the fine–coarse grid boundaries. In addition, nested grids might also be placed within the same model to allow for further refinement steps (Skamarock et al. 2005). Such an approach can also be viewed as a statically adaptive mesh application. It considerably improves the numerical consistency along the mesh interfaces due to the identical numerical setups in the nested domains. In addition, it allows for two-way mesh interactions that provide small-to-large scale feedbacks. In general, nested grid configurations make it possible to combine large-scale simulations with realistic meso-scale forecasts for selected domains. In contrast to dynamic AMR methods, the total number of grid points stays constant during a nested grid model run.
The goal of dynamic grid adaptations is to refine the mesh locally in advance of any important physical process that needs additional grid resolution, and to coarsen the grid if the additional resolution is no longer required. Dynamic AMR is therefore the most flexible variable-resolution global modelling technique that can readily vary the number of grid points as demanded by the adaptation criterion and the flow field. The grid interaction is automatically two-way interactive. Dynamically adaptive grids for atmospheric flows were first applied by Skamarock et al. (1989) and Skamarock & Klemp (1993) who discussed adaptive grid techniques for three-dimensional limited-area models in Cartesian geometry. More recently, Bacon et al. (2000) and Boybeyi et al. (2001) introduced the adaptive non-hydrostatic limited-area weather and dispersion model OMEGA for atmospheric transport and diffusion simulations. The model is based on unstructured, triangulated grids with rotated Cartesian coordinates that can be dynamically and statically adapted to features of interest. OMEGA has also been used as a regional hurricane forecasting system in spherical geometry (Gopalakrishnan et al. 2002).
In addition, a variety of statically and dynamically adaptive advection models and shallow water codes on the sphere have been proposed. A literature review of dynamically adaptive advection techniques is provided in Jablonowski et al. (2006) who applied a block-structured AMR method to a finite-volume advection algorithm. Most relevant to the advection studies presented here is the adaptive transport model by Hubbard & Nikiforakis (2003) who also used a block-structured finite-volume approach. Statically adaptive shallow water models on the sphere were developed by Ruge et al. (1995), Barros & Garcia (2004), Fournier et al. (2004), Giraldo & Warburton (2005), Ringler et al. (2008) and Weller & Weller (2008). Dynamically adaptive shallow water codes on the sphere have been discussed by Jablonowski (2004), Behrens et al. (2005), Jablonowski et al. (2006), Läuter et al. (2007), St-Cyr et al. (2008) and Heinze (2009). Behrens et al. (2005), Läuter et al. (2007) and Heinze (2009) utilized finite-element methods on unstructured triangulated meshes, whereas Jablonowski et al. (2006) and St-Cyr et al. (2008) introduced block-structured AMR techniques for finite-volume and high-order spectral-element methods. This paper is built upon the latter block-structured finite-volume approach and extends its applicability. In particular, the paper focuses on the use of adaptive blocks for two types of static mesh adaptations. These are the (i) so-called reduced grids that coarsen the latitude–longitude grid in longitudinal direction near polar regions and (ii) static refinements in storm-track regions for three-dimensional atmospheric flows. Reduced grids are attractive for computational performance reasons. This is particularly true if the Courant–Friedrichs–Lewy (CFL) stability condition in the polar regions limits the time step for the global domain. Reduced or ‘skipped’ grids have a long tradition in atmospheric general circulation modelling. They were first introduced in finite-difference models by Gates & Riegel (1962) and Kurihara (1965). The latter designed the so-called Kurihara grid with only four grid points next to the poles that offered an almost homogeneous grid distribution. But as noted by Shuman (1970) and Williamson & Browning (1973), the numerical design of reduced grids must be carefully chosen. Finite-difference models with curvilinear coordinates exhibit extreme curvature near the pole points. As a consequence, large truncation errors might be present in reduced grids. We will revisit this observation for finite-volume techniques.
This paper is organized as follows. Section 2 reviews the fundamental ideas behind the block-structured adaptive mesh methodology that is used in this paper. Both statically and dynamically adaptive meshes are supported, managed by a spherical mesh library for parallel computing architectures. The AMR technique is applied to a revised version of the finite-volume dynamical core, originally developed by Lin (2004) and Lin & Rood (1997). The model design and governing equations are briefly surveyed. The adaptive dynamical core is tested in two model configurations: the full three-dimensional hydrostatic dynamical core on the sphere and its corresponding single-level shallow water version. This shallow water setup serves as an ideal testbed for the two-dimensional adaptive mesh and reduced grid approach, which is discussed in §§3 and 4. In particular, several standard test cases from the Williamson et al. (1992) test suite, moving vortices (Nair & Jablonowski 2008) and a cross-polar flow test by McDonald & Bates (1989) are presented. In §5 the static adaptations in three-dimensional simulations are assessed using the idealized baroclinic wave test case by Jablonowski & Williamson (2006a,b). Section 6 provides a summary of the results and addresses future modelling challenges.
2. The adaptive finite-volume dynamical core
(a) Block-structured adaptive meshes and reduced grids on the sphere
The adaptive model design uses the spherical block-structured AMR grid library for parallel computing architectures by Oehmke & Stout (2001) and Oehmke (2004), which groups a latitude–longitude grid into horizontal, logically rectangular blocks. Adaptive blocks offer flexible non-uniform grids that have the potential to use the computing resources more efficiently than unstructured meshes or cell-based trees. In particular, the high potential for loop and cache optimization within a block, the reduced number of pointers to grid neighbours, and the quadtree-based block management with reduced parallel communication overhead render adaptive blocks highly competitive as also argued by Stout et al. (1997) and Colella et al. (2007). Adaptive blocks can be used for static and dynamic grid adaptations. Static adaptations are placed in user-defined regions of interest at the beginning of a forecast and are held fixed during the model run. Typical applications include static refinements in mountainous terrain or static refinements over selected geographical areas, e.g. for weather prediction applications. Dynamic adaptations can be used to track features of interest as they evolve. Examples are the life cycles of storm systems and cyclones or tracer transport. Adaptive blocks can also be used to define a special type of grid, the so-called reduced latitude–longitude grid. In a reduced grid, grid points are coarsened in longitudinal direction as the poles are approached. This effectively relaxes the strict time-step restrictions in a regular latitude–longitude grid if a CFL stability condition in the polar regions limits the time step for the global domain. In addition, it reduces the overall workload by reducing the number of grid points. From a high-level viewpoint, such a reduced mesh can be considered a statically adapted grid with one-dimensional mesh coarsenings.
Each adapted block is self-similar and contains an identical number of grid points. The number of grid points and initial blocks are user-defined and together determine the base resolution. In the experiments presented here 6×9 grid points per block in latitudinal × longitudinal direction are chosen. Furthermore, we select 6×8,12×16 or 24×32 initial blocks that correspond to the horizontal base resolutions 5°×5°,2.5°×2.5° and 1.25°×1.25°, respectively. In the event of refinement requests, a block is split into four, thereby doubling the horizontal resolution. Note that the vertical resolution is held fixed in the three-dimensional model runs. Coarsening requests reverse the two-dimensional refinement principle. Then four ‘children’ are coalesced into a single parent block, which reduces the grid resolution in each horizontal direction by a factor of 2. Neighbouring blocks can only differ by one refinement level. This leads to cascading refinement regions. More details on the block-structured adaptation technique, the finite-volume algorithm, the coarse–fine grid interface treatment including the flux corrections and interpolation techniques are provided in Jablonowski (2004), Jablonowski et al. (2006) and St-Cyr et al. (2008). The block-structured AMR algorithm ensures mass conservation. It is enforced via numerical flux corrections that average the fine grid fluxes and override the coarse grid flux at an adjacent grid interface. In this paper we focus on new aspects of the adapted blocks that address the reduced grid design and the performance of statically adapted grids in three-dimensional model runs.
The AMR grid library provides provisions for block-structured reduced grids in polar regions. Several of these reduced grid designs are explored here. As an example, two configurations are shown in figure 1a,b, which illustrates a reduced grid with base resolution 2.5°×2.5° plus one (rg1) and two (rg2) reduction levels. In addition, an adapted mesh with three static refinement levels is shown in figure 1c. The latter spans the resolutions 5°×5° to 0.625°×0.625° and is used for three-dimensional hydrostatic model simulations in §5. In the reduced grids, the grid reductions are placed at 75° and 60° (rg2) in both hemispheres. At each transition point the physical grid spacing Δx in longitudinal direction then changes abruptly by a factor of 2. Examples of the grid spacings are shown in table 1. The table also lists the grid distance near the poles, which determines the CFL numbers and maximum time steps in the advection examples in §3. This particular reduced grid design obeys the principle that the physical grid distances on the sphere do not greatly exceed the maximum grid spacing Δx≈278 km at the Equator. If wider grid distances are allowed in polar regions, the solution has a greater potential to degrade, especially in nonlinear simulations. The possible latitudinal positions of the reduction steps are restricted by the block-data design and determined by the initial number of blocks in both horizontal directions and the block size.
Each block has three additional cells, so-called ghost cells, in each horizontal direction that provide information from neighbouring blocks for the numerical scheme. The algorithm then loops over all blocks on the sphere before a Message Passing Interface (MPI) communication step on parallel computing architectures is triggered. This communication step updates the ghost cell region. The minimum block size that fulfils the parallel communication constraint is therefore 6×6 grid points per block. Note that the time step is the same in all blocks, which are treated as independent units. Depending on the particular adaptive or reduced grid interface, the ghost cell updates use conservative interpolation or averaging techniques, or just plain copies of the neighbouring data if the resolutions at an interface are the same. The parallel workload is managed by a load-balancing technique that can be called dynamically if required. It guarantees that an almost identical number of blocks is assigned to each processor on a parallel machine, which ensures the equal workload.
(b) The finite-volume dynamical core
We apply the block-structured adaptive mesh technique to a revised version of the finite-volume dynamical core, originally developed by Lin (2004) and Lin & Rood (1997). This hydrostatic dynamics package is based on a conservative two-dimensional finite-volume discretization in flux form and uses a floating Lagrangian coordinate in the vertical direction. The underlying advection scheme is upwind-biased and built upon the piecewise parabolic method (PPM) by Colella & Woodward (1984). This finite-volume advection algorithm uses a nonlinear monotonicity constraint that is a source of scale-selective numerical diffusion, especially near sharp gradients. Here the monotonicity option FFSL-4 from Lin & Rood (1996) is used. Other diffusion and filtering mechanisms include a horizontal second-order divergence damping mechanism, a fast Fourier transform (FFT) filter near the poles and a fourth-order Shapiro filter (Shapiro 1970) in the longitudinal direction polewards of 60°. The latter two stabilize the fast waves that originate from the pressure gradient terms and filter high-speed gravity waves at high latitudes. The filters are applied in tests of the full nonlinear shallow water and hydrostatic equations. They are neither needed nor applied in pure advection examples.
The PPM method is formally third-order accurate in each direction. Nevertheless, the order of accuracy reduces to second-order in two dimensions as shown by Jablonowski et al. (2006). The time-stepping scheme is explicit and stable for zonal and meridional Courant numbers CFL<1. This restriction arises since the semi-Lagrangian extension of the Lin & Rood (1996) advection algorithm is not used in the AMR model experiments. This keeps the width of the ghost cell regions small, but on the other hand requires small time steps if high wind speeds are present in polar regions. In this case the CFL condition is most restrictive due to the convergence of the meridians in the latitude–longitude grid.
From a high-level viewpoint, the three-dimensional model design can be viewed as vertically stacked shallow water layers. This provides an ideal framework for two-dimensional and three-dimensional studies since the two-dimensional shallow water model is equivalent to a single-level version of the three-dimensional design. The underlying hyperbolic shallow water system comprises the mass continuity equation and the momentum equation as shown in equations (2.1) and (2.2). Here the flux form of the mass conservation law and the vector-invariant form of the momentum equation are selected: 2.1 2.2 Here v is the spherical vector velocity with components u and v in the longitudinal (λ) and latitudinal (φ) directions, and Ωa=ζ+f denotes the absolute vorticity. The absolute vorticity is composed of the relative vorticity and the Coriolis parameter , where Ω is the angular velocity of the Earth. Furthermore, is the outward radial unit vector, ∇ and ∇· represent the horizontal gradient and divergence operators, and Φ=Φs+gh symbolizes the free surface geopotential, where Φs is the surface geopotential, h is the depth or mass of the fluid and g is gravity. In addition, D is the horizontal divergence and ν the divergence damping coefficient. A distinct advantage of this vector-invariant formulation is that the metric terms, which are singular at the poles in the chosen curvilinear spherical coordinate system, are hidden by the definition of the relative vorticity.
In three dimensions, the set of equations is very closely related to the shallow water system when replacing the height of the shallow water system with the hydrostatic pressure difference δp between two bounding Lagrangian surfaces in the vertical direction. Furthermore, the thermodynamic equation (2.4) in conservation form is added to the set: 2.3 2.4 2.5 Here Θ is the potential temperature and ρ denotes the density. The calculation of the pressure gradient terms (1/ρ)∇p+∇Φ is based on an integration over the pressure forces acting upon a finite volume. This approach constitutes the main difference between the shallow water system and the three-dimensional model setup. The underlying method has been proposed by Lin (1997, 2004) who defines the zonal and meridional components of the discretized pressure gradient terms as and .
In this primitive-equation formulation, the prognostic variables of the dynamical core are the wind components u and v, the potential temperature Θ and the hydrostatic pressure thickness δp. The geopotential Φ, on the other hand, is computed diagnostically via the numerical integration of the hydrostatic relation in pressure coordinates. It is important to note that the hydrostatic equation is the vertical coupling mechanism for the two-dimensional dynamical systems in each layer.
The pressure value at each Lagrangian surface can be directly derived when adding all the overlying pressure thicknesses within the vertical column. The pressure at the model top ptop is prescribed and set to 2.19 hPa in the current three-dimensional formulation with Nlev=26 vertical levels. There are a total of Nlev+1 Lagrangian surfaces that enclose Nlev Lagrangian layers. Each layer is allowed to float vertically as dictated by the hydrostatic flow and, as a consequence, the two Lagrangian surfaces bounding the finite volumes will deform over time. The displaced Lagrangian surfaces are then periodically mapped back to a fixed Eulerian reference system via monotonic and conservative interpolations. Lagrangian surfaces are non-penetrable material surfaces that do not allow transport processes across the boundaries. The periodic vertical remapping step essentially performs the function of a conservative vertical advection process as viewed from the fixed Eulerian reference frame. The lowermost Lagrangian surface coincides with the Earth’s surface.
Because of their equivalent designs, the two-dimensional shallow water version can be immediately extracted out of the full three-dimensional dynamical core with almost no changes of the source code. The main idea is to eliminate the influence of the thermodynamic equation in the momentum equation for the shallow water setup. The thermodynamic equation is linked to the momentum equation only via the computation of the geopotential gradient. Overall, three steps need to be performed. First, the number of vertical layers needs to be set to 1. Second, the potential temperature field Θ must be initialized with the constant 1 K and must not change during the shallow water run. Furthermore, the δp field now stands for the geopotential height h instead of a pressure thickness. Third, the gas constant for dry air Rd, the specific heat of dry air at constant pressure cp and need to be overwritten and set to unity. This guarantees that the discretized pressure gradient terms and in the three-dimensional momentum equations (Lin 1997, 2004) become identical to ∇h in the two-dimensional momentum equation. Additional details on the adaptive grid implementation of the finite-volume dynamical core can be found in Jablonowski (2004) and Jablonowski et al. (2006).
3. AMR and reduced grid advection tests
(a) Cosine bell advection test
The first test of the finite-volume advection scheme on reduced grids is based upon the solid body rotation of a cosine bell. This test is part of the standard shallow water test suite by Williamson et al. (1992) who describe the analytical initial conditions for the height h and the wind speeds (see test case 1). In brief, a cosine bell with peak amplitude of 1000 m is initially placed at (λc,φc)=(3π/2,0) and passively advected once around the sphere over a 12-day period. For the tests here, the flow orientation angle α=90° is selected that forces the cosine bell to cross both poles. This setup is therefore a stringent test of the reduction levels.
The cosine bell advection is tested on a 2.5°×2.5° base grid with one or two reduction levels positioned at 75° and 60° in both hemispheres. In addition, a single refinement level is overlaid in two experiments to dynamically track the path of the cosine bell. Refinements are triggered if one or more grid points within a block exceed the user-defined height threshold h≥53 m. Coarsenings are invoked if all grid points within a block fall below this value. Such an adaptation criterion has also been used in the cosine bell advection tests in Jablonowski et al. (2006) and St-Cyr et al. (2008) that serve as comparisons. Here the dynamic adaptations are primarily used to demonstrate their interplay with the reduced grid. In all experiments the time step is variable and matches the CFL number of 0.95.
Figure 2 shows snapshots of two reduced grid simulations with one and two reduction levels. Furthermore, a combined adapted and reduced grid (rg2) run is depicted. The composites display a daily time series of the cosine bell from day 1 through day 5 as it passes over the North Pole. The layout of the reduced blocks is indicated by the dotted grid lines. For clarity, no attempt is made to overlay the dynamic block adaptations in figure 2c. The figures show that the cosine bell passes over the North Pole without visible distortions or noise. The reduced grid rg1 has no negative impact on the quality of the simulations, which are almost indistinguishable from the 2.5° reference run (fig. 11a in Jablonowski et al. 2006). However, a minor degradation of the height field is apparent with two reduction levels (rg2). They slightly decrease the peak amplitude of the bell as it passes over the North Pole. The interplay between the reduced grid rg2 and dynamic adaptations is shown in figure 2c. It can be seen that the refinement region improves not only the shape but also the peak amplitude of the advected feature.
The tests are further quantified in table 2, which lists the normalized root mean square ℓ2 and infinity height error norms for the reduced, adapted and uniform grid runs after one revolution (12 days). The uniform-resolution simulations serve as reference runs. The norms are defined in Williamson et al. (1992) and Jablonowski et al. (2006) for the block-structured setup. The table also lists the final height of the cosine bell at day 12, the minimum and maximum number of blocks during the course of the simulation, the number of time steps for a 12-day run and CPU time on a single processor of a SUN Ultra 60 workstation. The latter three provide insight into the relative computational performance and workload of the model runs. They show that the reduced and adapted grid simulations exhibit speed-up factors between 6 and 16 when compared with their equivalent uniform-resolution runs. This is mainly attributable to the longer time steps permitted in the reduced and adapted grid runs. In these simulations the adaptive time step is determined by the smallest longitudinal grid spacing closest to the poles. Table 2 shows that all error measures in the 2.5° and 1.25° resolution groups are very similar. This confirms the observation that the reduced grids can be successfully used for the advection test. Here it is interesting to note that the adapted reduced grid simulations show slightly improved ℓ2 error measures despite the coarser longitudinal resolution at high latitudes. This effect is mainly attributable to the decrease in the total number of integration steps. The fewer invocations of the monotonicity constraint lessens the inherent numerical diffusion of the advection scheme and helps preserve the peak amplitude in the case of one reduction level. In the case of two reduction levels this time-stepping effect is overshadowed by the slight smoothing effect of the second grid reduction. The error measures in table 2 compare favourably to similar advection experiments documented in the literature. In particular, Rasch (1994) performed the cosine bell advection test with a reduced grid design and comparable base resolution (2.8°). The error measures for Rasch’s RG2.8 configuration with a monotonic upstream-biased transport scheme exceed the errors in table 2 by a factor of 1.5–2.
(b) Moving vortices on the sphere
A second more complex advection test confirms that reduced grids are viable candidates for tracer advection problems. The initial conditions are identical to those of Nair & Jablonowski (2008) and for brevity are not displayed here. In short, this advection test describes the roll-up of a passive tracer that is transported over the sphere over a 12-day time period. During this forecast period, the initially smooth passive tracer develops sharp gradients. They are depicted in figure 3a,b, which shows the initial field and analytical reference solution at day 12 after one revolution around the sphere. The advecting wind field is prescribed analytically. As before, the flow orientation angle is set to α=90°, which describes a translation over both pole points. Such a configuration challenges the reduced grid setups, especially at days 3 and 9 when the vortex centres cross both poles.
Figure 3c–f depicts snapshots of the tracer field at days 8 and 12, which are simulated with two reduced grid setups. These are the 2.5°×2.5° base resolution with two reduction levels at 75° and 60° (rg2) as shown in table 1, and the 1.25°×1.25° base resolution with four reductions at 82.5°,75°,67.5° and 60° in both hemispheres. Note that the polar grid spacing in the last configuration exceeds the maximum mesh spacing in the equatorial region with Δx≈139 km. In particular, the physical mesh distances Δx for the 1.25° base resolution and the four reduced grids are documented in table 3, which lists the mesh distances at the reduced grid transition points and at the position closest to the poles. Therefore, the configuration with four reductions tests the limits of the reduced grid in the block-structured AMR framework. The grid should not be reduced further or even be limited to three reduction levels in practice. The figures show that the reduced grid can be successfully employed for this more demanding advection test. The roll-up and solid body advection of the tracer are predicted reliably in both reduced grid simulations. The equivalent uniform grid runs at the resolutions 2.5°×2.5° and 1.25°×1.25° are visually indistinguishable from the reduced grid runs and therefore not shown here.
This observation is confirmed by global error statistics. The error measures of selected uniform and reduced grid simulations are listed in table 4. The table compares the normalized ℓ1,ℓ2 and error norms of the tracer distribution (Φ) after 12 days for different uniform and reduced grid simulations. In addition, it lists the constant length of the time step and the total number of time steps and blocks. They are relative measures of the computational cost when disregarding the overhead invoked by the adaptive mesh design of the reduced grids. The table shows that the error norms improve slightly in the reduced grid runs despite the coarsening of the polar grid spacing in longitudinal direction. This decrease in the error is again due to the reduced number of integration steps needed to cover the 12-day forecast period. As explained before, fewer integration steps invoke less numerical damping. Such damping is inherent in the monotonicity constraint of the finite-volume advection algorithm by Lin & Rood (1996). The reduced grid setups allow longer time steps since in this particular test the CFL stability condition is determined by the smallest grid spacing near the pole point. Here, all runs reflect identical CFL numbers of about 0.95. Doubling the longitudinal spacing allows for a doubling of the time step and speeds up the computation. In addition, the reduced number of blocks lessens the overall workload. The latter is partially compensated by the additional computational overhead that AMR invokes. The error measures in table 4 can also be compared to the dynamic AMR simulations by Nair & Jablonowski (2008) who list the errors for the flow orientation angle α=0° (table 3 in Nair & Jablonowski 2008). This flow follows a solid body rotation along the Equator where larger time steps can be used. They necessitate fewer integration steps that reduce the inherent numerical diffusion. Therefore, the error measures at day 12 with α=0° are lower at comparable resolutions and lie around ℓ2(Φ)=0.0226 and ℓ2(Φ)=0.0074 for the grid spacings 2.5° and 1.25°, respectively.
4. Test of the reduced grid in nonlinear shallow water simulations
In order to assess the performance of the reduced grid setups in the nonlinear shallow water system, three test scenarios are examined. These include the cross-polar rotating flow by McDonald & Bates (1989) as well as the standard shallow water tests 2 and 6 from the Williamson et al. (1992) test suite. The latter are the steady-state geostrophic flow and Rossby–Haurwitz wave.
(a) Cross-polar rotating high–low
First, the results of a uniform grid control run and reduced grid simulations are evaluated qualitatively using the McDonald & Bates (1989) cross-polar flow test. The test has also been used by Bates et al. (1990), Giraldo et al. (2002) and Nair et al. (2005). The initial condition consists of a geostrophically balanced flow pattern where the geopotential is given by 4.1 with gh0=5.768×104 m2 s−2. The Earth’s radius is denoted by a, the maximum wind speed is set to v0=20 m s−1. The wind components u and v can then be derived via the geostrophic relationship (see also Nair et al. 2005). This test does not have an analytical solution. It consists of two large waves, with the high wave in the west and the low wave in the east (in the Northern Hemisphere). The positions in the Southern Hemisphere are reversed. The wave rotates clockwise around the pole points and deforms slightly. After 5 days, the low and high waves exchange their positions and almost arrive at their starting locations at day 10. A cross-polar flow is well suited for the reduced grid assessments. The wind field reaches its maximum speed at the poles and exhibits strong gradients.
The results of the simulations in the Northern Hemisphere at day 10 are illustrated in figure 4, which displays the flow pattern of the uniform full grid (fg) control run (upper row) and the rg2 setup (bottom row). In particular, the geopotential height field and the wind components u and v (from left to right) are depicted. The reduced grid simulations with one (rg1, not shown) and two reduction levels (rg2) are visually almost indistinguishable from the full grid control run. There are no distortions in any of the large-scale flow patterns, which are still well resolved even in the rg2 case with only 36 grid points next to the poles. The time steps are Δt=240 s for the full grid and Δt=450 s for the reduced grid model runs. The time step cannot be reduced further in the rg2 simulations since the CFL condition at the poles no longer dominates the global time step. This purely qualitative assessment suggests that reduced grids are potentially viable for finite-volume simulations.
(b) Steady-state geostrophic flow
The aforementioned suggestion is only partly supported in the subsequent examples. The second test scenario is based on standard shallow water test 2 with flow orientation angle α=45°, which directs the maximum wind speeds of about 38.6 m s−1 across the mid-latitudes. The initial conditions are described in Williamson et al. (1992). In essence, the flow field consists of a large-scale steady-state geopotential wave in geostrophic balance, which is the analytical solution. The flow is expected to remain unchanged during the course of the integration and can be directly compared with its initial state. The polar regions exhibit a strong cross-polar flow with wind speeds of about 28 m s−1 at the 45° angle. This accentuates possible errors in the polar regions that result from the reduced grid coarsenings. It furthermore challenges the numerical scheme with maximum errors at the 45° angle. They are mainly due to the use of the latitude–longitude grid and the very nature of upwind-biased finite-volume schemes as argued in St-Cyr et al. (2008). Both reduced grid setups rg1 and rg2 have been tested as documented by the geopotential height and height error field in figure 5. The figure shows snapshots of the circulation at day 14 with overlaid blocks of the reduced grids. Visually, the contours of the geopotential height distributions in both reduced grid simulations are indistinguishable. In the case of one reduction level (rg1), a very moderate error pattern in the height field can be observed that is similar in full grid simulations (not shown). However, the errors increase considerably when adding a second reduction level (rg2). Here it is interesting to note that the errors of the reduced grid are concentrated in, but not confined to, the regions at high latitudes. Instead, the errors propagate and interact with other model errors over time so that the whole domain is affected. The time steps are Δt=200 s for the full grid and Δt=400 s for the rg1 and rg2 reduced grid simulations.
This qualitative assessment is further quantified in figure 6, which illustrates the normalized ℓ2 height error norms of the full grid and the reduced grid simulations. The temporal evolution of the normalized ℓ2 error norms confirms the significant increase in the solution error in the case of the rg2 setup, whereas error levels remain low for one reduction step. For comparison, Tolstykh (2002) also reports comparable error measures, with ℓ2(h)=4.3×10−4 at day 5 for uniform-grid finite-difference simulations at 2.5° resolution. In general, the error curves depict a linear increase in the solution error, which is overlaid by inertio-gravity wave oscillations. The latter are damped out over time. The sharp rise in the error level for rg2 is not tolerable for practical applications. Therefore, the rg2 design cannot be recommended for future nonlinear adaptive grid applications. This confirms the results by Lanser et al. (2000) who used this test with flow orientation angle α=90° to assess up to four reduction levels in a finite-volume approach. As in Lanser et al. (2000) we also find that the errors primarily arise in the velocity fields in polar regions due to the extreme curvature and metric terms in the spherical representation.
(c) Rossby–Haurwitz wave
The third assessment of the reduced grid design is based on shallow water test 6 that comprises a Rossby–Haurwitz wave with wavenumber 4. The initial conditions are described in Williamson et al. (1992). The test translates the symmetric wave pattern without change of shape from west to east. Figure 7 presents the geopotential height field at day 14 on a reduced grid with one reduction level as well as the time series of various normalized ℓ2(h) height error norms. In particular, the height error norms of the rg1 and rg2 reduced grid and fg uniform-grid control runs are compared. These error norms are computed with the help of the National Center for Atmospheric Research (NCAR) T511 spectral transform reference solution with quadratic transform grid (approx. 0.235° grid spacing) provided by the German Weather Service DWD. The solutions are available online at http://icon.enes.org/swm/stswm/node5.html as an archived network Common Data Form (NetCDF) dataset. Daily snapshots of the spectral transform simulation are provided. Here it is interesting to note that the ℓ2 height errors in the rg1 design almost overlay the full grid control run until the errors split after day 11. The overall solution error is therefore hardly affected by one reduction step, which is considered a viable option for nonlinear runs. In contrast, the rg2 model simulation again shows unacceptably large errors and cannot be recommended for future use.
The ℓ2 height error in both the uniform grid model run with ℓ2(h)=0.0033 and the rg1 simulation with ℓ2(h)=0.0036 at day 14 compare favourably to values published in the literature. These are listed for comparable resolutions in table 5. Besides from St-Cyr et al. (2008) these runs are compared with the NCAR spectral transform reference solution at the lower T213 triangular truncation (approx. 0.5625° grid spacing). Note that this reference solution has an uncertainty of about 0.0008 as assessed by Taylor et al. (1997) and Jakob-Chien et al. (1995).
Another aspect must be noted concerning the reduced grid setup. In test case 6 presented here, the time step Δt=360 s cannot be decreased in the two reduced grid simulations as seen before in the advection experiments or the polar flow shallow water problems. This is due to the fact that the time step is restricted not by the advection speed in polar regions, but by the gravity wave activity in mid-latitudes. Coarsening the polar regions therefore has no impact on the stability criterion and an identical time step must be used for all three simulations. This has implications with respect to the computational performance of the model runs. The clear speed-up of a reduced grid run in the advection experiment can no longer be achieved. On the contrary, choosing the rg1 reduced grid setup adds extra work at the reduced grid interfaces despite the 9 per cent reduction in the number of blocks. In Jablonowski (2004), the overhead was estimated to be approximately 27 per cent, which outweighs the advantages of the reduced number of blocks in the rg1 case. Therefore, it is not recommended to use a reduced grid if the global time step is not limited by the CFL restrictions in polar regions. However, it may be possible that future optimizations of the AMR library and the interpolation-averaging routines will reduce the overhead significantly. More detailed performance considerations are beyond the scope of this paper and a subject for future research.
In order to improve the accuracy of the zonal gradients in the immediate vicinity of the poles, special treatments could be introduced as suggested by Purser (1988). He pointed out that zonal derivatives can still be calculated with sufficient precision when high-frequency information is interpolated through an assumption of smoothness from nearest neighbours in the polar regions. None of these difficulties exist in spectral transform models, which need not calculate difference operators in grid point space. Instead, derivatives are evaluated in a spherical representation after applying Fourier and Legendre transforms. These so-called reduced Gaussian grids have been successfully used in spectral transform models for many years. Examples are discussed in Hortal & Simmons (1991) and Williamson & Rosinski (2000) who developed reduced grids for weather and climate applications. No significant loss of accuracy has been observed in comparison to full grid simulations. However, spectral transform models do not offer options for future adaptive or nested grids, which is the focus here.
5. Static refinements in the three-dimensional hydrostatic finite-volume dynamical core
The final assessment of the block-structured finite-volume design addresses the use of static adaptations in three-dimensional hydrostatic simulations. We place nested grids in a mid-latitudinal channel to test the influence of variable resolution on the evolution of baroclinic waves in the Northern Hemisphere. The baroclinic wave test has been proposed by Jablonowski & Williamson (2006a,b) (referred to as JW hereafter) who describe the initial analytical conditions. In short, the test starts from balanced and steady-state initial conditions with an overlaid perturbation in the zonal wind field. This perturbation triggers the growth of baroclinic waves over the course of 10 days. The balanced initial state is an analytical solution to the shallow atmosphere equations, which include both the hydrostatic primitive equations as well as the non-hydrostatic equation set. The flow comprises a zonally symmetric basic state with a jet in the mid-latitudes of each hemisphere and a quasi-realistic temperature distribution. Overall, the atmospheric conditions resemble the climatic state of a winter hemisphere reasonably well. The test design guarantees static, inertial and symmetric stability properties, but is unstable with respect to baroclinic or barotropic instability mechanisms. The baroclinic wave starts growing around day 4 and evolves rapidly thereafter, with explosive cyclogenesis at model day 8. The wave train breaks after day 9 and generates a full circulation in both hemispheres between days 20 and 30 depending on the model formulation.
All adaptive model runs use the base resolution 5°×5° that is statically refined with up to three refinement steps. The latter then bridges the resolutions 5°×5° to 0.625°×0.625° that have already been shown in figure 1c. Note that the finest resolution always covers the area between 30° and 75°N, which is the pre-selected storm-track region in the northern mid-latitudes. The refinements are confined to the horizontal directions so that the whole vertical column with 26 levels is refined in the event of refinement requests. Convergence studies by JW with uniform resolutions have shown that the representation of the baroclinic wave improves considerably with increasing horizontal resolution until the pattern converges at spatial resolutions around 1°–0.5°. It has also been observed that the flow is rather insensitive to the number of vertical levels, which justifies a purely horizontal refinement strategy. The same convergence behaviour is expected for the locally nested simulations.
The adaptive mesh designs are depicted in figure 8a–d by the overlaid adaptive blocks. The panels show snapshots of the surface pressure field at day 9 simulated with the nested grid resolution. In addition, uniform-resolution reference solutions with the same model are shown in figure 8e–h. The additional resolution in the nested grid simulations clearly helps intensify the baroclinic wave until the pattern converges at resolutions around 0.625°. Only minor deviations from the uniform grid reference solutions are visible. These are apparent in the regions south of 30° when transitioning to coarser resolutions. Overall, the adapted runs match the uniform grid simulations very closely. This is also confirmed in the time evolution of the minimum surface pressure depicted in figure 9. The figure shows the expected drop in minimum surface pressure at higher resolutions. The drop slows down considerably at the finest resolution, which indicates signs of convergence. The uniform and adapted model runs almost perfectly overlay each other.
A quantitative assessment of the adapted simulations is provided in figure 10. The figure displays the time evolution of the normalized ℓ2 and surface pressure error norms for both the nested and uniform model runs. The errors are computed with the help of the uniform high-resolution reference solution with 0.3125°×0.3125° grid spacing. The error norms are defined as in Williamson et al. (1992) and assess the global domain. Both the ℓ2 and errors show that the adapted runs start at a high error plateau before the baroclinic wave grows considerably after day 6. The initially flat error plateau is almost independent of the number of refinement steps. It is dominated by the 5°×5° errors of the coarse resolution in the Southern Hemisphere. During the simulation, the ℓ2 error cannot be significantly reduced below its early state at day 1, thereby indicating a lower limit for the mixed-resolution simulation. The nested runs can only slow down the error growth from day 6 onwards when the baroclinic wave contributes to the error pattern. At day 10 the errors in the nested simulations match the uniform runs very closely. However, the adapted run with three refinement levels shows slightly elevated ℓ2 error levels, partly due to the coarse grid error bound. This ℓ2 error gap at day 10 is further analysed below. In contrast, the error norms at day 10 almost perfectly overlay each other. They pick out the maximum error that lies in the refined region at day 10.
The errors displayed by the ℓ2 norm are further investigated in figure 11. The figure shows the time series of the contributions to the global error from (a) the Southern Hemisphere and (b) the refined region between 30° and 75°N. The ℓ2 errors in the Southern Hemisphere lie around an almost constant plateau in the adapted runs. This confirms the existence of the Southern Hemisphere error bound. The errors in the refined region (b) are more closely linked to the growing baroclinic wave. The aforementioned error gap at day 10 between the uniform-grid and refined 0.625° simulations becomes smaller without the contributions from the Southern Hemisphere. However, a perfect match in the refined region cannot be achieved. Figure 11b also displays that the errors in the adapted runs do not drop down to the error levels of the uniform runs during the early phases of the simulation. This suggests the existence of a second error bound before the growing baroclinic wave dominates the error pattern from day 6 onwards. This error bound is most likely linked to propagating truncation errors from the remaining coarse grid domain, wave reflections and refractions, and interpolation errors at fine–coarse grid boundaries that all interact with the flow in the fine-grid domain. Such errors are inherent in mixed-resolution simulations and must be kept small.
6. Summary and conclusions
A spherical, block-structured adaptive grid technique is assessed that has been applied to a revised version of the finite-volume dynamical core (Lin 2004) for weather and climate research. This hydrostatic dynamics package is based on a conservative and monotonic finite-volume discretization in flux form with vertically floating Lagrangian layers. It uses a latitude–longitude spherical grid that can be coarsened in longitudinal direction in polar regions to form a reduced grid. The adaptive model design is based upon the parallel AMR grid library by Oehmke (2004). The adaptive blocks can be used for static and dynamic grid adaptations. This paper focuses on the use of AMR techniques for reduced grids and static adaptations. Both the two-dimensional shallow water equations and the three-dimensional hydrostatic dynamical core are assessed. The shallow water setup serves as an ideal testbed for the simulations and allows the use of shallow water test cases like the advection of a cosine bell, moving vortices, a steady-state flow, the Rossby–Haurwitz wave or cross-polar flows.
The two-dimensional results show that reduced grid configurations are viable candidates for pure advection applications. The tracers are transported without visible distortions or noise across the reduced grid interfaces in polar regions. This extends the conclusion that reduced grids are suitable for passive advection (Williamson & Browning 1973) to finite-volume approaches. Several block-structured reduced grids are investigated that start from different base resolutions with up to four reduction levels. It is suggested that the physical resolution in longitudinal direction Δx in polar regions should not greatly exceed the maximum physical grid distance at the Equator. Otherwise, the transport algorithm might suffer from inaccuracies due to the enhanced spacing. It is also shown that reduced grid applications allow longer time steps if the simulation is bounded by a CFL restriction near the poles. Furthermore, reduced grids decrease the overall workload. This leads to computational speed-ups of the advection simulations by factors between 6 and 16 when compared with uniform resolution despite the computational overhead that AMR invokes. More detailed performance considerations and optimizations are beyond the scope of this paper. It is also shown that reduced grids should only be used moderately in nonlinear simulations. These simulations exhibit enhanced errors that are mainly generated in the velocity components in the polar regions due to the extreme curvature of curvilinear coordinates near the poles. The enhanced truncation errors propagate and interact with the flow in the whole model domain. Therefore, only very few reduction levels are acceptable for nonlinear simulations with a finite-volume approach in spherical geometry. The same applies to finite-difference models. Only spectral models do not exhibit this characteristic due to the very different and more accurate computation of the derivatives in spectral space. In addition, errors at fine–coarse grid interfaces arise in reduced grids. These interface errors are small but inherent in non-conforming AMR approaches with 2:1 or even higher jumps in the resolution. Smoothly varying unstructured grids or stretched grids can minimize such boundary effects.
This paper also demonstrates that static grid adaptations in three-dimensional model experiments can be successfully used to resolve growing baroclinic waves in the storm-track region. The baroclinic wave test by Jablonowski & Williamson (2006a,b) is used. Up to three refinement levels are tested in a mid-latitudinal periodic channel that span the resolutions 5°×5° to 0.625°×0.625° in the global domain. Global and regional error analyses of the adaptive model and uniform grid simulations are presented. They confirm that the nested domains reduce the errors considerably over the 10-day forecast period. The errors are mostly generated by the evolving baroclinic wave in the storm-track regime. The additional resolution clearly reduces both the global and regional error measures. The mixed-resolution regional errors compare very well with their corresponding uniform-resolution runs at the highest resolution. However, the mixed-resolution errors are bounded by contributions from the coarsest grid and additional errors triggered by wave reflections, refractions and interpolation at fine–coarse grid interfaces. Such errors are universal and inherent in mixed-resolution simulations (Ringler et al. 2008; Weller in press), especially for low-order (e.g. second-order) numerical schemes. The interface errors must be kept small and can be decreased or even hidden by switching to high-order numerical methods as demonstrated by Fournier et al. (2004) and St-Cyr et al. (2008) with spectral element schemes. The AMR interface effects are a subject of current and future research.
The work was partly supported by the Climate Change Prediction Program of the US Department of Energy (grant number DOE FG02 01 ER63248) and by the National Science Foundation (grants ATM-0723440 and ATM-0620065).
One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.
- © 2009 The Royal Society