Methods for assessment of models 08 May 2018
Methods for assessment of models  08 May 2018
Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – Part 1: Nonhydrostatic inertia–gravity modes
 Department of Atmospheric Science, Colorado State University, Fort Collins, CO 80523, USA
 Department of Atmospheric Science, Colorado State University, Fort Collins, CO 80523, USA
Correspondence: Celal S. Konor (csk@atmos.colostate.edu)
Hide author detailsCorrespondence: Celal S. Konor (csk@atmos.colostate.edu)
We have used a normalmode analysis to investigate the impacts of the horizontal and vertical discretizations on the numerical solutions of the nonhydrostatic anelastic inertia–gravity modes on a midlatitude f plane. The dispersion equations are derived from the linearized anelastic equations that are discretized on the Z, C, D, CD, (DC), A, E and B horizontal grids, and on the L and CP vertical grids. The effects of both horizontal grid spacing and vertical wavenumber are analyzed, and the role of nonhydrostatic effects is discussed. We also compare the results of the normalmode analyses with numerical solutions obtained by running linearized numerical models based on the various horizontal grids. The sources and behaviors of the computational modes in the numerical simulations are also examined.
Our normalmode analyses with the Z, C, D, A, E and B grids generally confirm the conclusions of previous shallowwater studies for the cycloneresolving scales (with low horizontal wavenumbers). We conclude that, aided by nonhydrostatic effects, the Z and C grids become overall more accurate for cloudresolving resolutions (with high horizontal wavenumbers) than for the cycloneresolving scales.
A companion paper, Part 2, discusses the impacts of the discretization on the Rossby modes on a midlatitude β plane.
 Article
(2330 KB)  Companion paper

Supplement
(8718 KB)  BibTeX
 EndNote
In the discretization of the governing equations of atmosphere models, we try to maintain as many physically important properties of the continuous system as possible. These include mimetic properties, such as selected identities from vector calculus; linear properties, such as the stability of the discrete systems, discrete representations of the geostrophic and hydrostatic adjustment processes, and discrete wave dispersion that faithfully imitates the continuous dispersion; and nonlinear properties, such as conservation of mass, energy and enstrophy.
The purpose of this paper to discuss the maintenance of the linear properties of the finitedifferenced nonhydrostatic equations on selected horizontal and vertical grids, with a particular focus on wave dispersion. Although the linear properties of a discretized system should not be the only factor in selecting a grid, they must be a key consideration. There have been numerous published studies on the horizontal discretization of the shallowwater equations (e.g., Winninghoff, 1968; Arakawa and Lamb, 1977; Mesinger and Arakawa, 1976; Randall, 1994; Skamarock, 2008; Thuburn, 2008; Thuburn et al., 2009; Weller et al., 2012) and on the vertical discretization of the quasihydrostatic equations (e.g., Tokioka, 1978; Arakawa and Moorthi, 1988; Hollingsworth, 1995; Arakawa and Konor, 1996; Konor and Arakawa, 2000). There are also a few published studies of the vertical discretization of the nonhydrostatic equations (e.g., Girard et al., 2014; Thuburn and Woolings, 2004; Toy and Randall, 2009). We are not aware of any previous publications on the horizontal discretization of the nonhydrostatic equations.
Arakawa and Winninghoff (Winninghoff, 1968; Arakawa and Lamb, 1977; Arakawa, 1988) defined the A, B, C, D and E grids based on the horizontal staggering of the variables. More recently, Randall (1994) and Lin and Rood (1997) have, respectively, added the Z grid and the CD grid to this list. We have also examined a DC grid to aid our analysis of the CD grid (defined later). For vertical discretization, there are only two grids available, namely the Lorenz grid (Lorenz, 1960; hereafter the L grid) and the Charney–Phillips grid (Charney and Phillips, 1953; hereafter the CP grid). In this paper, we discuss the discretization of the nonhydrostatic equations on these various horizontal and vertical grids.
Our analysis is based on the threedimensional anelastic system introduced by Lipps and Hemler (1982) because the curl of its pressuregradient force vanishes, which gives an important simplicity to our normalmode analyses based on the vorticity and divergence equations. The anelastic system also excludes the physically insignificant acoustic waves and the Lamb wave, so that we can focus on the much more important inertia–gravity and Rossby waves. Although the anelastic system has inaccuracies, as pointed out by Davies et al. (2003), Arakawa and Konor (2009), Dukowicz (2013) and Dubos and Voitus (2014), a wide range of nonhydrostatic phenomena are accurately described by the anelastic system. Arakawa and Konor (2009) show that among the features that can be well simulated by the anelastic system are the frequency of the inertia–gravity modes for all horizontal scales and the frequency of the middlelatitude Rossby modes for medium and small horizontal scales. Since discretization errors are mostly confined near the shortest resolvable scales, our findings can be directly applied to all nonhydrostatic systems.
The use of a threedimensional nonhydrostatic system to study the impact of the horizontal discretization on the dispersion of the waves has an important advantage over the twodimensional shallowwater system, in that it allows us to directly assess the performance of the discretization as a function of the vertical wavenumber, vertical resolution or vertical grid spacing. With the shallowwater system, only an indirect assessment can be made, based on the values of the Rossby radius of deformation divided by the horizontal grid spacing, as discussed by Arakawa and Lamb (1977), Randall (1994) and others. The results of our anelastic normalmode analyses are generally consistent with the shallowwater analyses discussed by previous authors for the cycloneresolving scales. Here, we present a threedimensional view of the dispersion of the modes by including an accurate assessment of the performance of the horizontal grids for the nonhydrostatic cloudresolving scales.
For large horizontal grid spacings, the results of our analyses are also applicable to the quasihydrostatic systems. The normalmode analyses presented by Arakawa and Konor (2009) show that the frequencies of the medium and largescale inertia–gravity waves, with horizontal wave lengths of approximately 60 km and larger, are nearly identical in the anelastic and quasihydrostatic systems. The dispersion of midlatitude Rossby waves is also nearly identical in the two systems, with the exception of the ultralong waves. The only major difference is that the quasihydrostatic system includes the Lamb wave, while the anelastic system does not. We ignore the Lamb wave in our analyses.
The dispersion of inertia–gravity waves leads to geostrophic and hydrostatic adjustments in the nonhydrostatic system, which is key to the maintenance of the geostrophic and hydrostatic balances, respectively. Thus, the accurate simulation of wave dispersion, i.e., the phase and group velocities are essential to maintain the approximate balances and also bolster the stability of the discrete system. Since we focus on the impact of the discretizations on the frequency of the waves, rather than their amplitude, no explicit diffusion is included in our analysis.
The discretization process can give rise to computational modes. The sources and behaviors of the computational modes differ for each horizontal and vertical grid. While a normalmode analysis is an excellent method to study the modification of the physical modes with the discretization, it is not completely adequate to study the computational modes. This is one motivation for including in our study analyses of the numerical solutions of the linearized anelastic equations.
We have used secondorder accurate finitedifference schemes for simplicity. The use of higherorder schemes would lead to minor quantitative differences in our results by reducing errors in the wellresolved scales, but it would not change the main conclusions of our analysis, which deal with errors near the smallest resolvable scales.
This paper discusses the dispersion of the inertia–gravity waves on a midlatitude f plane. Section 2 presents the linearized anelastic equations. Section 3 discusses the horizontal discretization of these equations on the Z, C, D, CD, (DC), A, E and B grids, in that order. At the end of Sect. 3, we present plots of the discrete dispersion relations and a comparison of the performance of the grids in simulating the inertia–gravity waves. Vertical discretization on the L and CP grids is discussed in Sect. 4. Section 5 presents a summary table of the discrete dispersion relations, to facilitate comparisons. In Sect. 6, we present results from simulations of the inertia–gravity modes obtained by running linearized numerical models based on the various horizontal grids and discuss the sources and behaviors of the various computational modes. Section 7 analyses the nonhydrostatic effects on the performance of the grids by comparing the dispersion of the inertia–gravity modes with the nonhydrostatic, quasihydrostatic and shallowwater systems. A summary and conclusions are provided in Sect. 8. Additional details are given in the Supplement. Finally, Part 2 (Konor and Randall, 2018) discusses the dispersion of Rossby modes on a midlatitude β plane.
The linearization of the Lipps and Hemler (1983) anelastic equations is with respect to a resting isothermal basic state, following Arakawa and Konor (2009), Davies et al. (2003) and many others. All variables are weighted by the square root of the basic state density ${\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}\left(z\right)$, so that the dispersion relation does not have an imaginary part.
2.1 Basic equations
The linearized horizontal momentum equation is
where v is the horizontal velocity, t is time, f is the Coriolis parameter, k is the vertical unit vector, ∇_{H} is the horizontal del operator, P is the pressure. Although the definition of P is not necessary for this study, in the Lipps–Hemler anelastic system, it is given by $P\equiv {\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}{c}_{\mathrm{p}}{\mathit{\theta}}_{\mathrm{0}}{\mathit{\pi}}^{\prime}$, where c_{p} is the specific heat of dry air under constant pressure, θ_{0}=θ_{0}(z) is the basic state potential temperature, ${\mathit{\pi}}^{\prime}\equiv \mathit{\pi}{\mathit{\pi}}_{\mathrm{0}}$ is the perturbation Exner pressure, $\mathit{\pi}\equiv {\left(p/{p}_{\mathrm{00}}\right)}^{\mathit{\kappa}}$ is the Exner pressure, p is the pressure, p_{00} is a standard pressure, $\mathit{\kappa}\equiv R/{c}_{\mathrm{p}}$, R is the gas constant, π_{0}=π_{0}(z) is the basic state Exner pressure, $\mathit{\theta}\equiv T/\mathit{\pi}$ is the potential temperature, and T is the temperature.
By taking $\mathit{k}\cdot {\mathrm{\nabla}}_{\mathrm{H}}\times \left(\mathrm{1}\right)$, we obtain the linearized vertical vorticity equation as
where ${\mathit{\omega}}_{z}\equiv \mathit{k}\cdot {\mathrm{\nabla}}_{\mathrm{H}}\times \mathit{v}$ is the vertical vorticity and $D\equiv {\mathrm{\nabla}}_{\mathrm{H}}\cdot \mathit{v}$ is the divergence of the horizontal wind. By applying the divergence operator to Eq. (1), we obtain the linearized divergence equation as
where ${\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}$ is the horizontal Laplacian operator. To predict the horizontal velocity, we can directly integrate Eq. (1). An alternative approach is to first predict vorticity and divergence through Eqs. (2) and (3), respectively, and then obtain the horizontal momentum by using horizontal elliptic solvers. To study the linear modes, we will use the vorticity and divergence equations instead of the momentum equation.
The vertical momentum equation is
where w is the vertical velocity weighted by ${\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}$ and $B\equiv {\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}g{\mathit{\theta}}^{\prime}/{\mathit{\theta}}_{\mathrm{0}}$ is the buoyancy, g is the gravitational acceleration, and ${\mathit{\theta}}^{\prime}\equiv {\mathit{\theta}}^{\prime}{\mathit{\theta}}_{\mathrm{0}}$ is the perturbation potential temperature. Note that in Eq. (4) $\left(\mathrm{1}/\mathrm{2}\right)\left(\mathrm{1}/{\mathit{\rho}}_{\mathrm{0}}\right)\left(\partial {\mathit{\rho}}_{\mathrm{0}}/\partial z\right)$ appears as a result of the weighting of the pressure term by ${\mathit{\rho}}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}$ with the linearization. The thermodynamic equation is
where ${N}^{\mathrm{2}}\equiv \left(g/{\mathit{\theta}}_{\mathrm{0}}\right)\left(\partial {\mathit{\theta}}_{\mathrm{0}}/\partial z\right)$ is the square of the Brunt–Väisälä frequency associated with the basic state. Both $\left(\mathrm{1}/{\mathit{\rho}}_{\mathrm{0}}\right)\left(\partial {\mathit{\rho}}_{\mathrm{0}}/\partial z\right)\equiv \mathrm{1}/H$ and ${N}^{\mathrm{2}}=\left(g/{\mathit{\theta}}_{\mathrm{0}}\right)\left(\partial {\mathit{\theta}}_{\mathrm{0}}/\partial z\right)=g\mathit{\kappa}/H$ are constants for an isothermal atmosphere, where $H\equiv R{T}_{\mathrm{00}}/g$ is the scale height of the isothermal basic state, and T_{00} is the threedimensionally and temporally constant temperature of the isothermal basic state. The anelastic continuity equation is
where $D\equiv {\mathrm{\nabla}}_{\mathrm{H}}\cdot \mathit{v}$ is the divergence.
Equations (4)–(6) can be combined into a single equation by eliminating B and w, giving
The threedimensional elliptic equation that can be used to diagnose the pressure term can be obtained by eliminating $\partial D/\partial t$ between Eq. (3) and $\partial /\partial t$ of Eq. (6) and using Eq. (4). The result is
Although the elliptic equation (Eq. 8) is not needed for the normalmode analysis, the numerical integration of the linearized equations requires the use of a discrete version of Eq. (8).
These equations have solutions that are steady ($\partial /\partial t=\mathrm{0})$, balanced (i.e., D=0) and quasistatic, given by
and
By using Eq. (9a) in Eq. (10), we find that
is also satisfied for the steady balanced state. The majority of the discrete systems that we will discuss have steady balancedstate solutions corresponding to Eq. (9a)–(9c). However, we found some exceptions with the CDgrid discretization, depending on the timeintegration scheme. We discuss this point further later and give examples of the schemes that do or do not have steady balanced solutions.
To obtain the dispersion relation for the inertia–gravity modes on a middlelatitude f plane, we seek solutions of Eqs. (2), (3) and (7) of the form
where ϕ represents an arbitrary variable, Re indicates the real part of a complex function, $\widehat{\mathrm{\Phi}}$ is the complex amplitude, $\underset{\sim}{i}\equiv \sqrt{\mathrm{1}}$, k , ℓ and m are the wavenumbers of waves propagating in x, y and z directions, respectively, and ν is the frequency of the waves. Since we assume that the upper and lower boundaries are material surfaces (so that ${w}_{\text{S}}={w}_{\text{T}}=\mathrm{0})$, we define the vertical wavenumber m as
where n is the “integer wavenumber,” and z_{T} is the height of the upper boundary. The use of n is very convenient for the study of discretization because n is equal to the number of layers in a numerical model.
2.2 Dispersion of inertia–gravity waves on a midlatitude f plane
By applying Eq. (11) to ω_{z}, D and P in Eqs. (2), (3) and (7) and seeking nontrivial solutions, we obtain the dispersion relation as
Here, we have used $\left(\mathrm{1}/{\mathit{\rho}}_{\mathrm{0}}\right)\left(\partial {\mathit{\rho}}_{\mathrm{0}}/\partial z\right)\equiv \mathrm{1}/H$ and ${N}^{\mathrm{2}}=\left(g/{\mathit{\theta}}_{\mathrm{0}}\right)\left(\partial {\mathit{\theta}}_{\mathrm{0}}/\partial z\right)=g\mathit{\kappa}/H$. It should be noted that the dispersion equation (Eq. 13) is obtained by seeking nontrivial solutions for D among Eqs. (2), (3) and (7). If we seek nontrivial solutions for ω_{z} or P, we will find an additional mode with ν=0, which corresponds to the steadystate balanced solution. This comment applies to all discrete systems with the exception of some on the CD grid. We will briefly discuss the exceptions with the CD grid. According to Eq. (13), the frequency of inertia–gravity modes is bounded between the Brunt–Väisälä frequency (N) and the inertial frequency (f). For large horizontal and small vertical wavenumbers, the “gravitational term”, i.e., the first term in the numerator of Eq. (13), dominates the “inertial term”. The effect of rotation on these gravity modes is small.
The modes of the continuous system will be referred to as physical modes or true modes. The discrete solutions contain versions of these physical modes that are modified by discretization, and they may also contain additional computational modes (or solutions) that are purely the result of discretization. In some cases, the computational modes are dynamically inert and confined to the smallest resolvable scales. In other cases, the computational modes appear in the form of nonunique (or multiple) physical solutions that do not interact in the linear system. These nonunique solutions can also be dynamically inert on the smallest resolvable scales. The generation and behavior of computational modes in the horizontally discrete systems will be discussed in detail at the end this section.
We will try to avoid the repetition of results that have been previously presented by other authors. Many additional aspects of the A, B, C, D and E grids are discussed by Arakawa (2000, and the references therein). For additional discussion on the Z grid, the reader is referred to Randall (1994) and Chen (2015). Lin and Rood (1997) and Skamarock (2008) can be consulted for additional discussions on the CD grid. For recent applications of the C grid, the paper by Weller et al. (2012) and the references therein are recommended.
We use the vorticity and divergence equations in our normalmode analyses for all grids, since the inertia–gravity modes are primarily controlled by the vorticity, divergence and pressure (or mass), and the horizontal velocity is determined by the vorticity and divergence.
Next, we discuss the dispersion and modification of the physical modes in the discrete solutions through a normalmode analysis similar to the one used for the continuous system. Fig. 1 shows portions of the horizontal grids discussed in this paper.
3.1 Solutions on the Z grid
We choose to discuss the Z grid first because, as pointed out by Randall (1994), it yields better solutions than the A, B, C, D and E grids in the discretization of the shallowwater equations on an f plane. Instead of predicting the horizontal velocity, a model that uses the Z grid predicts the vorticity and divergence without staggering. This approach requires the use of elliptic solvers to obtain the horizontal velocity from the predicted vorticity and divergence. The elliptic solvers add 10–20 % to the computational cost of the dynamical core.
3.1.1 Discrete dispersion of the inertia–gravity waves on an f plane
The vertically and temporally continuous and horizontally discrete versions of vorticity and divergence equations given by Eqs. (2) and (3), respectively, and Eq. (7) can be written on the Z grid shown in Fig. 1a as
By using
in Eqs. (14)–(16), and seeking nontrivial solutions, we can obtain the discrete dispersion relation as
where the factors
arise from the finitedifference form of the Laplacian of P in the divergence equation (Eq. 15). As the grid spacing approaches zero ($kd,\phantom{\rule{0.125em}{0ex}}\mathrm{\ell}d\to \mathrm{0})$ and therefore ${\mathit{\xi}}^{\mathrm{2}}={\mathit{\eta}}^{\mathrm{2}}\to \mathrm{1}$ according to Eq. (18), the discrete dispersion relation Eq. (18) approaches the continuous dispersion equation (Eq. 13), which confirms that the discrete solution is consistent. This property is shared by all of the dispersion relations that we discuss in the paper. For the smallest resolvable horizontal scales, for which $kd=\phantom{\rule{0.125em}{0ex}}\mathrm{\ell}d=\mathit{\pi}$, the factors in the Laplacian yield ${\mathit{\xi}}^{\mathrm{2}}=\phantom{\rule{0.125em}{0ex}}{\mathit{\eta}}^{\mathrm{2}}=\mathrm{4}/{\mathit{\pi}}^{\mathrm{2}}\sim \mathrm{0.4}$. Since these factors are bounded between 0 and 1, all discrete modes and all variables move or oscillate with a single frequency that is between N and f. We conclude that the discrete inertia–gravity modes are purely physical; i.e., they correspond to solutions of the continuous system, although they are slightly distorted by discretization errors.
The use of a higherorder Laplacian in Eq. (15) would not change the form of the discrete dispersion (Eq. 18), but it would change the definitions of the factors defined in Eq. (18) such that errors would be slightly improved near and at the smallest resolvable horizontal scale (hereafter, SRHS). This comment applies to all of the discrete Laplacian in the dispersion relations discussed in the paper.
For each grid, we present plots of the absolute value of the frequency as a function of ${k}^{*}\equiv \sqrt{{k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}}$, as given by the dispersion relations obtained for the various schemes for the given integer vertical wavenumber, n. The true frequency given by Eq. (13) is also superimposed in these plots to facilitate assessment of the results. We use k=ℓ to make the plots, which means that the modes would be oriented diagonally on a square horizontal grid. The scale height of the basic state is assumed to be H=24 km (or ${N}^{\mathrm{2}}=\mathrm{1.16}\times {\mathrm{10}}^{\mathrm{4}}$ s^{−2}), which corresponds to a typical lowertropospheric stability. The height of the domain is z_{T}=80 km, and the Coriolis parameter is $f={\mathrm{10}}^{\mathrm{4}}$ s^{−1}. The smallest vertical grid spacing is given by $\mathit{\delta}z={z}_{\text{T}}/n$. We present the frequency plots for the four different horizontal grid spacings d=2 km, d=10 km, d=25 km and d=100 km. A grid spacing of 2 km is selected to represent the cloudresolving models. The grid spacings of 10 and 25 km are selected to be representative of mesoscale and global numerical weather prediction models. The grid spacing of 100 km is selected to be representative of climate models. In these plots, we show the frequency lines for 10 integer vertical wavenumbers ranging from n=1 to n=1280, and the corresponding vertical grid spacings are also indicated.
Figure 2 shows that, in the Zgrid solutions, the discrete frequencies associated with small horizontal wavenumbers are virtually identical to the true frequencies for all vertical wavenumbers and horizontal grid spacings. This is expected because the discretization errors must be negligible away from the SRHS. For horizontal wavenumbers that are near and at the SRHS (indicated by thin dashed vertical lines in the plots), the discrete frequencies deviate downward from the true frequencies, and the deviations become larger as the vertical scales of the modes become smaller.
The group velocity, $\partial \left\mathit{\nu}\right/\partial {k}^{*}$, is the velocity with which the wave energy moves. If the group velocity in the discrete system is substantially slowed down or reversed, relative to that of the continuous system, then wave energy can be spuriously trapped in the source regions, which may lead to instability. For brevity, we only qualitatively assess the simulated group velocity as implied by the frequency plots. The group velocity of the inertia–gravity waves in the true solutions is always positive (or zero) as shown by the thin lines in Fig. 2. The Zgrid solutions maintain positive (or zero) group velocity, but near the SRHS the Zgrid group velocity is generally smaller than the true group velocity.
These results generally support Randall's (1994) conclusions regarding the Z grid with the low horizontal wavenumbers. Our results suggest that the nonhydrostatic effects reduce the errors for the deep modes with high horizontal wavenumbers. This will be further discussed in Sect. 7.
3.2 Solutions on the C grid
The C grid is arguably the most commonly used grid in atmosphere models. The Mintz–Arakawa UCLA general circulation model (GCM; Arakawa, 2000), developed starting in about 1965, was the first GCM and possibly the first atmospheric model to use the C grid (Akio Arakawa, personal communication, 2017). With the C grid, the normal components of the horizontal velocity, i.e., u and v, are predicted at the meridional and latitudinal walls, respectively, as shown in Fig. 1b. The mass, vertical velocity and thermodynamic variables are predicted at the cell centers. Following the standard practice, the vorticity and divergence are defined at the corners and centers of the cells on the C grid, respectively (see Fig. 1b). Since the mass and divergence are predicted at the cell centers and the normal component of velocity is predicted at the cell walls, the C grid is well suited for controlling the divergent component of the velocity. Although the same is not true for the rotational component of the velocity, some improvements have been developed in recent years (e.g., Adcroft et al., 1999; Thuburn et al., 2009; Weller et al., 2012).
Discrete dispersion of the inertia–gravity waves on an f plane
We define the discrete vorticity and divergence in terms of the components of the horizontal velocity on the C grid (Fig. 1b) as
and
respectively. Then, the vertically and temporally continuous and horizontally discrete versions of the vorticity and divergence equations corresponding to Eqs. (2) and (3) can be written as
and
By using Eq. (17) and
in Eqs. (21), (22) and (16), we obtain the discrete dispersion equation for the inertia–gravity waves on the C grid as
where ξ and η are defined by Eq. (18) and
The factor μ^{2} appears in Eq. (24) as a result of the averaging of the vorticity and divergence to each other's grid points. As k (and/or ℓ) approaches ${k}_{max}={\mathrm{\ell}}_{max}=\mathit{\pi}/d$, corresponding to the SRHS, μ approaches zero, which reduces the effect of rotation. For ${k}_{max}={\mathrm{\ell}}_{max}=\mathit{\pi}/d$, so that μ=0, the rotation term completely vanishes. This means that the interactions between the divergence and vorticity are broken. For this case, the frequency obtained by Eq. (24) applies only to the divergence and buoyancy, while the frequency associated with the vorticity is arbitrary. We will return to this point in Sect. 6.
The use of a higherorder Laplacian would slightly improve the errors for the wellresolved modes, as discussed earlier, but because of the presence of the averaging represented by μ, there would be no improvement of errors at the SRHS. The same comment applies to averaging on all of the other grids discussed in this paper.
The dispersion relations for the C grid are plotted in Fig. 3. By comparing the four panels of the figure, the first thing we see is that the character of the solutions changes with the horizontal grid spacing. For the d=2 km case shown in Fig. 3a, the Cgrid frequency solution is almost identical to the Zgrid solution. The solutions differ for the modes with vertical integer wavenumber n=1947, for which the Cgrid solutions generate negative group velocities ($\partial \left\mathit{\nu}\right/\partial {k}^{*}<\mathrm{0})$. This means that the nonhydrostatic Cgrid models using a horizontal grid spacing of 2 km and a vertical grid spacing of δz=41 m or larger (with a domain depth of z_{T}=80 km) produce wave solutions nearly as accurate as those of the Z grid. There is one issue with the C grid that applies for all horizontal grid spacings: the vorticity and divergence decouple at the SRHS (indicated with double dashed lines in the figures) because this mode cannot recognize the effects of rotation. Figure 3b–d show that for d=10 km, d=25 km and d=100 km the group velocities are reversed for modes shallower than n=390 (δz=205 m) for 10 km horizontal grid spacing, n=156 (δz=512 m) for 25 km case and n=39 (δz=2051 m) for 100 km horizontal grid spacing.
As the horizontal grid spacing of a Cgrid model decreases, the solution approaches the true solution for an increasingly wide range of vertical scales. This is because the effects of rotation are unimportant for the shortest modes on sufficiently fine horizontal grids. As a result, the effect of the averaging factor μ^{2} becomes negligible in the inertial term, and therefore the difference between the C and Zgrid physical solutions virtually disappears. As the horizontal scale of the modes approaches the SRHS, the coupling of vorticity and divergence becomes weaker due to the averaging and finally disappears at the SRHS. We classify the mode with the SRHS as dynamically inert.
These results generally support conclusions of previous studies of the behavior of the C grid for the low horizontal wavenumbers (Arakawa and Lamb, 1977; Randall, 1994). However, the nonhydrostatic effects on deep modes with small horizontal scales contribute to the reduction of the errors with the C grid, as with the Z grid. The reasons for this are discussed in Sect. 7.
The discretization of momentum equation based on the C grid staggering on the hexagonal and triangular grids presents the challenge of dealing with the unbalanced number of mass and velocity points. Gassmann (2011) shows that the discrete system of linear equations can be closed for the hexagonal C grid yielding a similar discrete dispersion equation to the one for the Cartesian C grid given by Eq. (18). Achieving a similar closure for the triangular C grid appears to be more challenging than the hexagonal version, however.
3.3 Solutions on the D grid
We discuss the D grid here in preparation for a discussion of the CD grid. The D grid was originally proposed for the convenience of maintaining the geostrophic balance between the wind and mass fields. The Mintz–Arakawa model developed around 1963 used the D grid (Akio Arakawa, private communication, 2017), but we are not aware of any other atmospheric model that has used it since then. On the D grid, the longitudinal component of velocity u and the meridional component v are paired with the meridional and longitudinal finitedifference derivatives of P, respectively (Fig. 1c). However, the placement of the velocity components is consistent with the placement of the divergence at the corners of the grid cells. This leads to the most unattractive and problematic feature of the D grid: the mass and/or divergence has to be averaged in the horizontally discrete continuity equation.
On the D grid, the vertical velocity can be placed either at the centers with the mass and thermodynamic variables or at the corners with the divergence. In this paper, we primarily discuss the former option, in which the vertical velocity is placed at the centers, but we discuss the latter at the end of this subsection.
Discrete dispersion of the inertia–gravity waves on an f plane
We define the discrete vorticity and divergence from the components of the horizontal velocity on the D grid (Fig. 1c) as
and
respectively. The vertically and temporally continuous and horizontally discrete versions of vorticity and divergence equations that are, respectively, given by Eqs. (2) and (3) can be written as
and
where ${\stackrel{\mathrm{\u203e}}{P}}_{i+/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv \frac{\mathrm{1}}{\mathrm{4}}\left({P}_{i,j}+{P}_{i+\mathrm{1},j}+{P}_{i,j+\mathrm{1}}+{P}_{i+\mathrm{1},j+\mathrm{1}}\right)$. The discrete version of Eq. (7) is
Using Eqs. (17) and (23) in Eqs. (27)–(29), we obtain the discrete dispersion relation for the inertia–gravity modes on the D grid as
where ξ and η are defined by Eq. (18), and μ is defined by Eq. (25). As with the C grid, the inertial term is multiplied by the averaging factor μ^{2} in Eq. (30). The Laplacian terms are multiplied by μ^{2} because the discrete Laplacian in Eq. (28) involves averaging. As discussed below, the multiplication of the gravitational and inertial terms by μ^{2} in the numerator of Eq. (30) creates a stationary mode at the SRHS.
The frequency plots for the D grid are shown in Fig. 4. Since the inertia–gravity solutions are qualitatively similar for all horizontal grid spacings (unlike the solutions on the C grid), the following discussion applies to all grid spacings. The Dgrid solutions for all vertical scales converge to a dynamically inert solution in all variables at the SRHS, while the true solutions are oscillating modes (as indicated by the double dashed lines in Fig. 4). As mentioned above, the zero frequency at the SRHS is a result of the factor μ^{2} in the numerator of Eq. (30). Since all variables have zero frequency, with or without rotation, the impact of the dynamically inert modes can be more severe than with the C grid. As a consequence, the group velocities of the Dgrid solutions are badly reversed near the SRHS.
These results are in agreement with the conclusions of Arakawa and Lamb (1977) and Randall (1994) for the D grid.
The configuration of the D grid for which the vertical velocity is defined at the corners, where the divergence is predicted, yields a dispersion relation with a numerator identical to that of Eq. (30) but without μ^{2} in the denominator. The frequency plots for this case are qualitatively similar to those shown in Fig. 4; they are omitted here for brevity. A detailed discussion of this version of the D grid, including the corresponding frequency plots, can be found in the Supplement.
3.4 Solutions on the CD grid
The CD grid proposed by Lin and Rood (1997) was first used in a fluxform semiLagrangian dynamical core. It is now also used in Geophysical Fluid Dynamics Laboratory (GFDL)'s FiniteVolume CubedSphere Dynamical Core (FV3) and in the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM). FV3 will also be used in the nextgeneration numerical prediction model of the US National Centers for Environmental Prediction. The CD grid is built on a predictor–corrector timeintegration scheme, in which the C and Dgrid discretizations are used for the predictor and corrector steps, respectively (see Fig. 1d). Because of this, an analysis of the dispersion properties of the CD grid has to include temporal discretization as well as spatial discretization.
Lin and Rood (1997) asserted that the CD grid combines the C grid's superiority for the prediction of the divergent part of velocity with the D grid's superiority for the prediction of the rotational part. They claimed that the solutions of the CD grid should be close to those of the Z grid. Skamarock (2008) examined the normal modes of the linearized shallowwater equations discretized on the CD grid for the timecontinuous case and concluded that the linear properties of the CD grid resemble those of the D grid rather than the Z or C grids.
In this section, we examine the performance of the CD grid as inferred from a normalmode analysis using secondorder accurate differences. As mentioned earlier, the main conclusions of our analysis apply without change to systems discretized with higherorder schemes. Guided by the two papers mentioned above, we temporally and horizontally discretize the anelastic equations on the CD grid. For the predictor step, we use the horizontally discrete vorticity and divergence equations on the C grid, as given by Eqs. (21) and (22), respectively. Then we horizontally discretize the vertical momentum equation (Eq. 4), the buoyancy equation (Eq. 5) and the continuity equation (Eq. 6). For the corrector step, we use the horizontally discrete vorticity and divergence equations on the D grid, as given by Eqs. (27) and (28), respectively. We also horizontally discretize the vertical momentum equation (Eq. 4), and the buoyancy equation (Eq. 5) and continuity equation (Eq. 6) on the D grid. For the predictor and corrector steps, the temporal discretizations can be summarized in compact form as
and
respectively. In Eqs. (31a) and (32b), ϕ is an arbitrary variable, the superscripts (n) and (n+1) denote the time steps (the notational conflict with the vertical wavenumber should not be confusing), (*) denotes a provisional value, τ is the time step, and F and G represent the tendency terms. Equation (31a) indicates that the provisional values are obtained by advancing a half a time step. All variables denoted by (n) and (n+1) are defined on the D grid, and those denoted by (*) are defined on the C grid. If Eq. (31a) is applied to the vorticity and divergence equations on the C grid, the initial values with superscript (n) need to be averaged to the C grid from the D grid. It should be noted that, on the D grid, the provisional values of vorticity (and divergence) do not have to be averaged for the use in the tendency term of the divergence (and vorticity) equations. We seek solutions of these equations in the following forms:
and
where (n) is the time step counter, and (*) denotes the provisional value of the complex amplitude of the variables after half a time step of integration in the prediction step on the C grid.
Skipping a few derivation steps (see the Supplement), we can write the discrete equations that predict the complex amplitudes as follows:
Predictor step on the C grid:
Corrector step on the D grid:
where ξ and ηare given by Eq. (18), μ is given by Eq. (25).
In Eqs. (32) and (33), the averaging factor μ appears as a result of averaging of the vorticity and divergence from the D grid to the C grid. The provisionally predicted vorticity ${\widehat{\mathit{\omega}}}_{z}^{(*)}$ and divergence ${\widehat{D}}^{(*)}$ reside at the corners and the centers, respectively. They are used without averaging in the divergence equation (Eq. 39) and vorticity equation (Eq. 38) on the D grid. The pressures on the two grids, denoted by the subscripts C and D, are diagnosed at the cell centers with the buoyancy and vertical velocity, but these two pressures are, in principal, different from each other. Since the pressure is a diagnostic variable, no time stamp is included.
We have tested several variations of the predictor–corrector timeintegration algorithm outlined above. There are two main reasons for doing so. First, we do not have access to the details of how the Lin and Rood (1997) scheme is implemented, and second, we want to examine the extent to which the solutions depend on the timeintegration schemes. We hope that the results of our analysis will provide some guidance to future users of the CD grid. The construction of the five schemes is summarized in Table 1.
Scheme I is the simplest. Its equations are given by Eqs. (32)–(42), with the assumption that $\widehat{P}\equiv {\widehat{P}}_{\text{C}}={\widehat{P}}_{\text{D}}$. To derive the dispersion equations, we first eliminate the variables with (*) in the correctionstep equations by using the equations from the predictor step. Then we use
where $\widehat{\mathrm{\Phi}}$ represents ${\widehat{\mathit{\omega}}}_{z}$, $\widehat{D}$ and $\widehat{B}$, to obtain the discrete dispersion relation for complex frequencies as
where we use the following definitions to shorten the equation:
The (complex) dispersion equation (Eq. 45) is obtained by requiring nontrivial solutions for $\widehat{D}$. The nontrivial solutions for ${\widehat{\mathit{\omega}}}_{z}$ and $\widehat{B}$ produce an additional mode, satisfying ${e}^{\underset{\sim}{i}\mathit{\nu}\mathit{\tau}}\mathrm{1}=\mathrm{0}$, which corresponds to a steady and neutral mode associated with the steady balanced state. By using ${e}^{\underset{\sim}{i}\mathit{\nu}\mathit{\tau}}\mathrm{1}=\mathrm{0}$ in Eq. (45), where ν_{r} and ν_{i} are the frequency and the rate of amplification/decay of a mode, respectively, we obtain the dispersion relations that govern the frequency and the amplification factor as
and
respectively.
A Newton–Raphson iteration is used to find the frequency ν_{r} that satisfies Eqs. (47a) and (48b) simultaneously. As a check, we also use a simple search algorithm to find the frequency. The two methods give the virtually the same unique solution unless the time step τ is very large. Although we obtain real frequencies for moderately unstable amplification factors (${e}^{{\mathit{\nu}}_{i}\mathit{\tau}}>\mathrm{1})$, we limit the discussion here to nearneutral cases (${e}^{{\mathit{\nu}}_{i}\mathit{\tau}}\approx \mathrm{1})$. Before discussing the results, we mention that the dispersion relations given by Eqs. (47a) and (48b) yield a solution satisfying ν_{r}=0 for the SRHS (μ=0), which indicates that the discrete solutions include a dynamically inert mode. The solutions of Eq. (48b) with ν_{r}→0 and μ=0 show that this mode is very weakly damped (${e}^{{\mathit{\nu}}_{i}\mathit{\tau}}={\mathit{\sigma}}_{f}<\mathrm{1})$. Recall that ${\mathit{\sigma}}_{f}\equiv \mathrm{1}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\tau}}^{\mathrm{2}}{f}^{\mathrm{2}}$, so that the damping rate is considered negligible for the time steps typically used in atmosphere models.
The steady balanced solution can be obtained by assuming nondivergent motion (${\widehat{D}}^{\left(n+\mathrm{1}\right)}={\widehat{D}}^{\left(n\right)}=\mathrm{0})$ and a steady and quasistatic state (${\widehat{\mathit{\omega}}}_{z}^{\left(n+\mathrm{1}\right)}{\widehat{\mathit{\omega}}}_{z}^{\left(n\right)}=\mathrm{0}$ and ${\widehat{B}}^{\left(n+\mathrm{1}\right)}{\widehat{B}}^{\left(n\right)}=\mathrm{0})$. To satisfy Eq. (42), the vertical velocity must also vanish (${\widehat{w}}^{\left(n+\mathrm{1}\right)}={\widehat{w}}^{\left(n\right)}=\mathrm{0})$ in the balanced state. We obtain the unique balanced solution given by
These are the discrete equivalents of the continuous relations given by Eq. (9a)–(11c). We conclude that the balancedstate solution is maintained by Scheme I.
Results
The results presented in this section, as inferred from our normalmode analyses, demonstrate that the dispersion of the inertia–gravity waves produced on the CD grid with various timeintegration schemes is in all cases close to the one produced by the D grid. This is in agreement with the conclusions of Skamarock (2008). For some of the timeintegration choices (not discussed here), we could not uniquely determine a steady balancedstate solution. For some choices, we could not even find real frequencies.
The frequency plots for the inertia–gravity modes on the CD grid with Scheme I are shown in Fig. 5. The time steps τ used in computing these frequencies are 50, 120, 150 and 210 s for the 2, 10, 25 and 100 km cases, respectively. Schemes II, III, IV and V produce very similar frequency plots (not shown here but included in the Supplement). Schemes III and V pass the provisional values of divergence from the C grid to D grid but do not pass the provisional values of vorticity. Thus, the vorticity is solely predicted on the D grid with these schemes. All of these plots show a strong qualitative resemblance to those for the D grid shown in Fig. 4, and therefore the discussion presented for the Dgrid case can be applied to the CD grid as well.
Solutions on the DC grid
Our normalmode analysis shows that, in agreement with the results of Skamarock (2008), the CD grid behaves similarly to the D grid rather than the C grid. One possible explanation is that the dispersion of the inertia–gravity waves on the CD grid is almost entirely dominated by the corrector step. To test this hypothesis, we constructed a “DC grid” by reversing the integration sequence used in the CD grid. With the DC grid, the D and C grids are used in the predictor and corrector steps, respectively. One of the solutions that corresponds to Scheme II (discussed above) produces a Cgridlike discrete dispersion relation for the DCgrid sequence. A detailed discussion of the DC grid can be found in the Supplement.
We do not advocate the use of the DC grid because we do not think it has any advantage over the C grid, and the C grid is both simpler and computationally less expensive than the DC grid.
3.5 Solutions on the A grid
The A grid is a completely unstaggered grid on which the mass and the zonal and meridional wind components are predicted/diagnosed at the same grid points (see Fig. 1e). The quadrilateral A grid can be viewed as the superposition of four noninteracting C grids or Z grids (see the Supplement). We are not aware of any major model that currently uses the quadrilateral A grid, but the icosahedral A grid is used by the Nonhydrostatic ICosahedral Atmospheric Model (NICAM) (Satoh et al., 2008), and the FlowFollowing Icosahedral Model (FIM) and Nonhydrostatic Icosahedral Model (NIM) developed at the National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) (Lee and MacDonald, 2009).
Discrete dispersion of the inertia–gravity waves on an f plane
We define the discrete vorticity and divergence from the components of the horizontal velocity on the A grid shown in Fig. 1e as
and
respectively. The vertically and temporally continuous and horizontally discrete versions of the vorticity and divergence equations, given by Eqs. (2) and (3), respectively, can be written as
and
By including Eq. (16), which is also the discrete version of Eq. (7) on the A grid, we obtain the discrete dispersion relation for the inertia–gravity modes on the A grid (following derivations parallel to those for the Z, C and D grids) as
where
The dispersion relation for the A grid resembles that of the Z grid except that the factors ξ and η that arise with the Z grid are replaced by $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ and $\stackrel{\mathrm{\u0303}}{\mathit{\eta}}$, respectively. The factors $\stackrel{\mathrm{\u0303}}{\mathit{\xi}}$ and $\stackrel{\mathrm{\u0303}}{\mathit{\eta}}$ become zero for the SRHS ($kd=\mathrm{\ell}d=\mathit{\pi})$, which yields ν=f for all vertical scales. The modes with the SRHS oscillate as purely inertial modes. In the absence of rotation, these modes do not propagate or oscillate. Although it cannot be detected by the normalmode analysis, the A grid generates multiple, linearly decoupled (or nonunique) solutions. This point is discussed further in Sect. 3.8. We will show the plots of the frequencies obtained with the A grid in Sect. 3.9, along with the corresponding plots for the E and B grids.
3.6 Solutions on the E grid
The E grid can be viewed as the superposition of two C grids, in which the cell centers of one C grid are placed at the corners of a second C grid, as indicated by solid boundary lines in Fig. 1f. On the resulting E grid, the zonal and meridional components of the horizontal velocity are predicted at the same points; this can be seen as an advantage over the C grid. As noted by Arakawa and Lamb (1977), the E grid is identical to a B grid rotated by 45^{∘} (shown in Fig. 1g) with a grid spacing of $d/\sqrt{\mathrm{2}}$. The borders of the cells of the B grid are indicated by the dashed lines in Fig. 1f. We will see next that from the vorticity and divergence point of view the E grid can also be viewed as a superposition of two independent and noninteracting Z grids.
The E grid was used in operational models of the US National Centers for Environmental Prediction for more than three decades (Janjic, 1984). We are not aware of any model that currently uses the E grid.
Discrete dispersion of the inertia–gravity waves on an f plane
We first define the discrete vorticity and divergence from the components of the horizontal velocity for the integer (i,j) and halfinteger ($i+\mathrm{1}/\mathrm{2},\phantom{\rule{0.125em}{0ex}}j+\mathrm{1}/\mathrm{2})$ center points of the E grid shown in Fig. 1f. We than write the discrete versions of the divergence, vorticity and mass equations, respectively, Eqs. (2), (3) and (7), for the integer and halfinteger points as follows:
Equations for the centers (integer points) of E grid:
Equations for the centers (halfinteger points) of E grid:
By using Eq. (17) or Eq. (23) separately in Eqs. (54)–(56) and (58)–(60), we obtain the same discrete dispersion relation for the inertia–gravity modes at the integer and halfinteger points as
where the factors ξ and η are defined by Eq. (18). The dispersion equation (Eq. 62) is identical to that of the Z grid given by Eq. (18). The important difference with the E grid, however, is that there are two solutions: one for the integer points and the other for the halfinteger points. This can be seen by a comparison of the discrete equations for the integer points (Eqs. 54–56) with the corresponding equations for the halfinteger points (Eqs. 58–60). The discrete equations for the integer points use only information from integer points, and those for the halfinteger points use only information from the halfinteger points. Of course, this means that the solution includes computational modes arising from multiple (or nonunique) solutions. Further discussion of the computationalmode on the E grid is given in Sect. 3.8.
3.7 Solutions on the B grid
On the B grid, the mass and thermodynamic variables are predicted at the cell centers while the horizontal velocity components are predicted at the cell corners (see Fig. 1g). The fact that the two components of the horizontal velocity are predicted at the same points is convenient for the Coriolis terms of the momentum equations. In the linearized system considered here, the vorticity and divergence are both predicted at the cell centers.
A version of the Goddard Institute for Space Studies global circulation model used the B grid (Hansen et al., 1983). The B grid has been used in many ocean models and can perform as well as the C grid in ocean models because the Rossby radius of deformation can be smaller than the grid spacing in ocean applications (Randall, 1994; Arakawa and Lamb, 1977). This is now changing as ocean models move to higher horizontal resolutions (e.g., Adcroft et al., 2016).
Discrete dispersion of the inertia–gravity waves on an f plane
We define the discrete versions of the vorticity and divergence on the B grid shown in Fig. 1g as
and
respectively. The vorticity equation is identical to that of the Z grid (Eq. 14) and E grid (Eq. 56). The divergence equation for the B grid is
The discrete Laplacian in the divergence equation (Eq. 64) uses information from the cells sharing corners instead of walls, which distinguishes the B grid from of the all other grids considered here. The B grid also uses the same mass equation as the Zgrid equation (Eq. 16) and Egrid equation (Eq. 56), which are also used with C and A grids. By using Eq. (17) in the vorticity equation (Eq. 14), divergence equation (Eq. 64) and mass equation (Eq. 56), we obtain the discrete dispersion relation for the inertia–gravity modes of the B grid as
where the factors ξ and η are defined by Eq. (18). The Laplacian term ${\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}{d}^{\mathrm{2}}{\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}$ becomes zero for the smallest resolvable horizontal scales $\left(kd=\mathrm{\ell}d=\mathit{\pi}\right)$, which leads to ν=f for all vertical scales. As with the A grid, the modes with the smallest resolvable horizontal scales appear as purely inertial modes. There are also multiple solutions with the B grid, as well as computational modes, which will be discussed next.
3.8 Computational modes on the A, B and E grids
The most easily recognizable feature of the computational modes is that they are dynamically inert at the SRHS. Consider the Laplacian operators of the A and B grids used on the righthand sides of Eqs. (51) and (64), respectively. The Laplacian operator of the A grid yields a zero result for patterns (1), (2), (3) and (4) shown in Fig. 6, presuming that the red and blue colors correspond to same magnitudes with opposite signs (also see Fig. 1e). These patterns are caused by the existence of four independent solutions as demonstrated in the Supplement. Because these four patterns give a Laplacian of zero, they remain unchanged. They are the computational modes of the A grid. On the other hand, the Laplacian of the B grid gives zero results for pattern (3), and therefore this pattern is the computational mode of the B grid. The solutions on the red and blue grid networks are completely decoupled from each other on the A and B grids (also see Fig. 1g). For a rotating case, the A and B grids produce inertial oscillations at all grid points.
The E grid holds two decoupled solutions: one for the network of centers (say red cells) governed by Eqs. (52)–(56) and the other for the network of corners (say blue cells) governed by Eqs. (58)–(60). Pattern (5) is the computational mode of the E grid. These decoupled solutions are also indicated by different colors in Fig. 1f.
The A grid can produce four independent and noninteracting (nonunique) solutions as indicated with different colors in Fig. 1e, and the B and E grids can produce two independent and noninteracting (nonunique) solutions as indicated with different colors in Fig. 1f and g, respectively. To avoid the separation of solutions, a horizontal mixing process such as diffusion is needed, but such a process can have adverse side effects such as unrealistic dissipation of the smallscale physical modes.
The Z, C, D and CD grids do not have independent and noninteracting solutions like those of the A, B and E grids. The C, D and CD grids do have dynamically inert modes for the SRHS because of the averaging of the variables.
3.9 An illustrative discussion of the dispersion of modes obtained with the A, E and B grids
The dispersion relations of the A, E and Bgrid solutions are free from the averaging factor μ (see Eqs. 52, 62 and 65) because these grids do not require averaging in their divergence and mass equations. All these grids yield ν=f at their SRHS, and consequently their frequency plots resemble each other. The frequency plots of the inertia–gravity modes for the A, E and B grids (see Fig. 1e, f and g) are shown in Figs. 7, 8 and 9, respectively. The horizontal grid spacing $\sqrt{\mathrm{2}}d$ is used in the E grid plots to maintain the same density of cell centers as the other grids. The frequency of all vertical scales, particularly deep vertical modes, makes a very sharp change near the SRHS. This sharp change causes the group velocity to reverse severely near the SRHS for all vertical modes.
In light of the discussion in the previous subsection, we conclude that the frequencies of the modes obtained using the A, E and B grids are not unique. Solutions as shown in Fig. 6 at the neighboring grid points can be decoupled from each other although they satisfy the same dispersion relation.
The conclusions of our normalmode analysis of the nonhydrostatic anelastic inertia–gravity wave solutions obtained with the A, E and B grids are not substantially different from the conclusions of earlier studies of the normalmodes of the shallowwater equations as described by many authors, including Arakawa and Lamb (1977) and Randall (1994).
There are two vertical grids available for the vertical discretization of models based on the height and pressure vertical coordinates. These are the Lorenz grid (or L grid; Lorenz, 1960) and the Charney–Phillips (or CP grid; Charney and Phillip, 1953). The vertical grid used with isentropic coordinates is analogous to the CP grid (see Konor and Arakawa, 1997, 2000; Thuburn and Woolings, 2004; Toy and Randall, 2007). The difference between the L and CP grids is that the L grid predicts the thermodynamic variables (in this paper the buoyancy) at the model “layers”, which are denoted by the dashed lines in Fig. 10, along with the horizontal momentum or vorticity and divergence and mass variable or pressure. The CP grid, on the other hand, predicts the thermodynamic variables at the model “interfaces”, which are the solid lines in Fig. 10, along with the vertical momentum.
The Lorenz (1960) grid is the most commonly used vertical grid for atmospheric models. Notable examples of the nonhydrostatic models that use the L grid are the Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012) and NICAM (Satoh et al., 2008). Many studies (e.g., Arakawa and Moorthi, 1988; Hollingsworth, 1995; Arakawa and Konor, 1996) have raised issues with the L grid. All of these studies find evidence of a dynamically inert mode in the vertical structure of the potential temperature. Girard et al. (2014) report an unusually high number of occurrences of the 2 dz wave in the vertical structure of simulations with GEM (a nonhydrostatic model), which is consistent with a dynamically inert mode. Arakawa and Moorthi (1988) and Arakawa and Konor (1996) additionally find evidence of a spurious rapid growth of short baroclinic modes in the Lgrid models. We will not discuss this subject further because it is outside the scope of this paper.
The alternative to the L grid is the CP grid, which is free of these problems. Examples of nonhydrostatic models that use the CP grid are the UK Met Office models (Wood et al., 2014) and GEM (Girard et al., 2014). Toy and Randall (2009) also use the CP grid in their vertical discretization.
Our purpose in this section is to present a comparison of the L and CP grids with the same normalmode analysis that we used for the comparison of the horizontal grids, so that we can assess the relative impact of the horizontal and vertical grids on the solutions of the nonhydrostatic anelastic equations.
4.1 L grid
The vertically discrete, temporally and horizontally continuous versions of the vorticity equation (Eq. 2), divergence equation (Eq. 3), vertical momentum equation (Eq. 4), thermodynamic equation (Eq. 5) and continuity equation (Eq. 6) on the L grid shown in Fig. 10 can be written as
and
respectively. Recall that $\left(\mathrm{1}/{\mathit{\rho}}_{\mathrm{0}}\right)\left(\partial {\mathit{\rho}}_{\mathrm{0}}/\partial z\right)=\mathrm{1}/H$ for an isothermal atmosphere. By using
in Eqs. (66)–(70), we obtain the discrete dispersion relation for the inertia–gravity modes on the L grid as
where
The vertical averages of buoyancy in the vertical momentum equation (Eq. 68) and the vertical velocity in the thermodynamic equation (Eq. 69) both result in the appearance of the averaging factor ${\mathit{\mu}}_{z}^{\mathrm{2}}$ in the gravity term, which prevents the smallest vertical scale (mδz=π) from recognizing the stability because ${\mathit{\mu}}_{z}^{\mathrm{2}}=\mathrm{0}$ for this scale. This permits a dynamically inert solution in which the buoyancy is completely decoupled from the rest of the variables. The frequency obtained from Eq. (72) in this case applies to the variables other than the buoyancy. The term 1∕4H^{2} is multiplied by the averaging factor ${\mathit{\mu}}_{z}^{\mathrm{2}}$ because the pressure and vertical velocity are vertically averaged in Eqs. (68) and (70), respectively, which slightly modifies the dispersion.
4.2 CP grid
The vertically discrete, and temporally and horizontally continuous, versions of the vorticity, divergence and continuity equations on the CP grid shown in Fig. 10 are given by Eqs. (66), (67) and (70), and the vertical momentum and thermodynamic equations can be written as
and
respectively. Following the procedure for the Lgrid case, we obtain the discrete dispersion relation for the inertia–gravity modes as
where ${\mathit{\mu}}_{z}^{\mathrm{2}}$ and ζ^{2} are given by Eq. (72). The main difference between the dispersion equation for the CP grid and that for the L grid is that the averaging factor ${\mathit{\mu}}_{z}^{\mathrm{2}}$ does not appear in the gravity term of the CPgrid dispersion equation given by Eq. (76), and thus the dynamically inert mode obtained with the L grid does not exist on the CP grid.
Figure 11 shows the frequencies as functions of composite horizontal wavenumber of inertia–gravity modes obtained with the L and CP grids. The true frequencies are also shown in separate panels of the figure. The figure shows the results for two vertical integer wavenumbers (or number of layers), namely ${n}_{max}=\mathrm{320}$ and 80. We included additional frequency lines corresponding to more integer vertical wavenumbers than were used in the plots of Sect. 3 (indicated by thinner solid lines in the plots). In the Lgrid solutions shown in Fig. 11b and e, the frequency of the smallest vertical resolvable mode, identified by n_{max}, deviates greatly from the true frequency, which yields values equal to or less than inertial frequency. As the vertical scale approaches the smallest resolvable scale, the modes gradually lose their ability to recognize the effects of buoyancy. For the mode with the smallest scale, the buoyancy is completely decoupled from the wind field; for that mode, the buoyancy is dynamically inert. The group velocity of the mode is reversed for the short horizontal scales. In contrast, the frequency of the CPgrid solutions is generally close to the true frequency but slightly higher.
Because the frequency errors with the L grid increase with increasing horizontal wavenumber, we can expect that the adverse effect of the use of the L grid on the dispersion of the inertia–gravity waves will be more problematic in highresolution nonhydrostatic applications than in lowresolution quasihydrostatic ones.
Table 2 summarizes the discrete dispersion relations obtained for the horizontal and the vertical grids. The range of the horizontal and vertical wavenumbers normalized by the horizontal and vertical grid spacing is indicated for each dispersion equation.
In this section, we first discuss the construction of linear anelastic numerical models based on the various horizontal grids and then present analyses of simulations of inertia–gravity modes on a middlelatitude f plane. Our purpose is to confirm the results of the normalmode analyses of the discrete equations and to investigate issues that are not completely revealed by the normalmode analysis. In the main text, we discuss only results for the Z, C, D and CD grids. Additional discussion on the A, B and E grids is given in the Supplement.
6.1 Equations of the C, D and CDgrid models
We write the horizontally discrete version equations (Eqs. 2–6 and 8) on the C and D grids as
Equations of the Cgrid model:
where
and
Equations of the Dgrid model:
where
and
The C grid involves averaging twice. First, we must average D from the centers to the ω_{z} points at the cell corners for use in the vorticity equation, Eq. (76), and second, we must average ω_{z} from corners to the D points at the cell centers for use in the divergence equation, Eq. (77).
Consider the perturbation patterns consisting of vertical and horizontal stripes or a checkerboard for the buoyancy $\stackrel{\mathrm{\u0303}}{B}$ (or divergence D) on the C grid. Such a perturbation given to D at the cell centers cannot be recognized at the cell corners where ω_{z} is predicted. These patterns therefore behave as pure gravity modes rather than inertia–gravity modes. We interpret these as physical modes behaving badly.
The D grid requires four separate fourpoint averages in four equations. Two of these are the averages of D from corners to the ω_{z} and B points at cell centers. The third one is the averaging of ω_{z} from the centers to the D points at the corners. The fourth one is the averaging of $\stackrel{\mathrm{\u0303}}{B}$ to the corners, for use on the righthand side of the elliptic equation. In comparison to the C grid, the D grid requires two additional averages in the continuity equation and in the elliptic equation that determines the pressure, P. We can avoid averaging in the continuity equation, but then we have to use averaging in the thermodynamic equation. If the perturbation patterns mentioned above are given on the D grid, they cannot be recognized, and so they appear as stationary, dynamically inert patterns.
Now we write the temporally and horizontally discrete version equations (Eqs. 2–6 and 8) on the CD grids as
Predictor step (*) on the C grid:
where ${\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}$ is given by Eq. (80). In the above equations,
Corrector step (n+1) on the D grid:
where ${\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}$ is given by Eq. (85). In the above equations,
On the CD grid, the principal prognostic variables that are denoted by (n) and (n+1), reside on the D grid. The provisional values of variables that are defined on the C grid are denoted by (*). The variables ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n\right)}$ in Eq. (86) and ${\stackrel{\mathrm{\u203e}}{D}}_{i,j}^{\left(n\right)}$ in Eq. (87) are averaged to the grid points where they are used on the C grid at the beginning of the predictor step. At the correction step on the D grid, the provisional value of D^{(*)} is directly used in the prediction of ω_{z}.
The equations given above belong to our primary model We also built a secondary model, in which the singleunderlined terms are replaced by ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}^{(*)}$, the doubleunderlined terms are replaced by ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{(*)}$, and the tripleunderlined term is replaced by $\stackrel{\mathrm{\u203e}}{{D}_{i,j}^{(*)}}$. The primary and secondary models are built by closely following Schemes V and II, respectively. These were discussed in Sect. 3. The two models produce identical numerical solutions, as discussed in the next section.
6.2 Checking the dispersion equations through numerical integrations
We have integrated the models described above to check the normalmode solutions discussed in Sect. 3 and to illustrate the behavior of the computational modes associated with the different grids. For the time discretization of the all models other than the CDgrid model, we use a simple forward timeintegration scheme that produces virtually neutral (very weakly unstable) solutions with short time steps in all models. We also tested other timedifferencing schemes in a single test case. The trapezoidal scheme produced neutral solutions very similar to those obtained with the forward scheme. A forward–backward predictor–corrector scheme was excessively dissipative, and a secondorder Adams–Bashforth scheme was slightly dissipative.
We first present results from the standing oscillation simulations, to look for the frequencies obtained through the normalmode analyses. The standing oscillations appear stationary because they are produced by the waves that propagate with the same phase speed in all directions on the horizontal plane. In these simulations, the vertical structure is continuous, and the vertical wavenumber n is prescribed. The simulations start from the initial buoyancy fields shown in Fig. 12, for which we have selected the horizontal wavelengths of L=4 km, L=20 km, L=50 km and L=200 km in both the x and y directions. These wavelengths correspond to the SRHS of d= 2 km, d= 10 km, d= 25 km and d= 100 km, respectively, which were examined in connection with the normalmode analyses. Here, we discuss results for L=4 km and L=200 km. Similar discussions for L=20 km and L=50 km can be found in the Supplement. The vertical wavenumbers used in the simulations are n= 320, 640 and 1280 for L=4 km, and n= 80, 160 and 320 for L=200 km. The sensitivity of the numerical solutions to the grid spacing is examined by using three different grid spacings, which are $d=L/\mathrm{80}$, $d=L/\mathrm{4}$ and $d=L/\mathrm{2}$. The initial perturbation field for $\stackrel{\mathrm{\u0303}}{B}$ (or any other variable) is shown in Fig. 12 for these three horizontal resolutions. With the highest horizontal resolution ($d=L/\mathrm{80})$ shown in Fig. 12a, the initial waves are well resolved and appear smooth. For the grid spacing $d=L/\mathrm{4}$, the initial field is not completely resolved. The initial field is only recognized as ones, zeros and negative ones, which appears like alternating upsidedown and inverted “pyramids” in Fig. 12b. The grid points marked by plus signs (+) are at the positive and negative extremes, and at the wall centers and the corners of the (upsidedown and inverted) pyramids. For $d=L/\mathrm{2}$, which corresponds to the shortest horizontal grid spacing to resolve the initial perturbation, the pyramidlike structure is again visible with the grid points at the extremes of the (upsidedown and inverted) pyramids (Fig. 12c).
The numerical frequencies of the standing oscillations for the high horizontal resolution case ($d=L/\mathrm{80})$ simulated by Z, C, CD and Dgrid models are given in Table 3. For comparison, we include in the table the true frequency obtained from Eq. (11). The numerical frequencies obtained with the C, CD and Dgrid models are very close to the true frequencies. This is expected because of the use of high horizontal resolutions in the simulations. All variables oscillate with the same frequency, and there is a unique solution. The numerical frequencies obtained with the CD grid are identical to those obtained with the D grid.
We have repeated the same simulation using $d=L/\mathrm{4}$, which is half of the shortest spacing needed to resolve the initial perturbation. The numerical frequencies obtained with the Z, C, CD and Dgrid models for the same horizontal resolution are tabulated in Table 4, which shows that all variables at all the grid points oscillate with the same frequency in the Zgrid simulation. In the C, CD and Dgrid solutions, all variables are also oscillating with the same frequency, but some of the modes have frequencies lower than the inertial frequency (underlined numbers), which indicates that these modes cannot properly recognize the rotation. This is generally due to the averaging of the divergence and vorticity to each other's grid points. In the L=4 km case, the C and Z grids produce similar frequencies, but the CD and D grids produce much lower frequencies. In the L=200 km case, the C grid produces lower frequencies than the Z grid, sometimes even lower than the inertial frequencies. The CD and D grids produce lower frequencies than the C grid. Note that the CD and Dgrid frequencies are almost identical to each other.
Finally we tabulate results obtained by using the shortest possible horizontal grid spacing ($d=L/\mathrm{2})$ to resolve the initial perturbation shown Fig. 12c in Table 5. This is the horizontal grid spacing for which the errors due to finite differencing are the largest, and the computational modes impact the solutions most strongly. The D and CDgrid solutions do not oscillate (indicated by zeros) but the Z grid produces oscillations. Averaging of the initial buoyancy perturbation to the divergence points completely wipes out the wave on the D and CD grids. The checkerboard pattern is dynamically inert. Nonoscillating solutions are also obtained by starting from the vorticity perturbations (as indicated by zeros within parentheses in Table 5). In the Cgrid solution, the buoyancy and divergence recognize the initial perturbation and do produce oscillations, although the vorticity is decoupled from the divergence and buoyancy, due to the averaging of the divergence to the vorticity points. For the short horizontal scales, the frequencies of the buoyancy and divergence in the Cgrid solutions are very close to those of the Zgrid solutions. For the long horizontal and short vertical scales, the frequency with which the buoyancy and divergence oscillate in the Cgrid solutions is considerably smaller than that of the Zgrid solutions and the inertial frequency (indicated by underlined numbers in Table 5). If the initial perturbation is given to the vorticity instead of the buoyancy in the C, CD and Dgrid simulations, the solutions are nonoscillatory for all variables (indicated by zeros within parentheses in Table 5). The Zgrid solution produces frequencies that are close to the true frequencies even though the initial perturbation is poorly resolved.
6.3 A numerical simulation of wave propagation
We have also made simulations to demonstrate the behavior of the computational modes during the propagation of inertia–gravity modes with the seven grids under discussion. These simulations start from a positive bellshaped buoyancy perturbation placed in the middle of the horizontal domain with rapidly decaying amplitude away from the center. We modified the bellshaped perturbation by setting it to zero at every other grid point. We refer to this modification as “initial gridscale noise”. The horizontal domain is 280 by 280 grid points, and the horizontal grid spacing is d=50 m, except that for the E grid we use d=70.71 m. The radius of the perturbation, that is the distance from the peak of the perturbation to where the perturbation becomes zero, is 950 m. The perturbation is continuous in the vertical, with the integer vertical wavenumber n=320. Figure 13 shows the buoyancy field in the 124by124 wide cornerend portion of the horizontal domain, after 100 simulated minutes of integration. For reference, we show the Zgrid result obtained without the superimposed initial gridscale noise. No major difference can be seen between the Zgrid solutions started with the computational mode and without it in the portion of the domain shown in Fig. 13. In the solution with the gridscale noise, there is, however, a remnant of the initial noise, which takes the form of gridscale standing inertia–gravity oscillation at and near the center of the domain where the initial perturbation is strongest. This is not visible in Fig. 13 because it is outside of the plotted domain. The noise gradually subsides in time as the bulk of the initial perturbation propagates outward from the center of the domain. The Cgrid simulation result is very close to the Zgrid result. It is evident that the initial noise does not have a permanent effect on the solutions with the Z and C grids.
The CD and Dgrid results resemble the Cgrid result except that the noise is apparent near the center of the domain in both simulations (see Fig. 13). Since the noise with the SRHS cannot propagate or oscillate, it is a permanent feature of the CD and Dgrid simulations. By setting the initial perturbation at every other grid point to zero, the dynamically inert computational mode described by pattern (3) in Fig. 6 is put in place on the A, B and E grids. The other solution in these grids recognizes the initial Gaussian perturbation and yields a propagating wave that has similar features to those produced by the Z, C, CD and D grids.
6.4 Numerical simulations of waves that are forced by noisy heating patterns
To assess the performance of the Z, C, D and CD grids in the presence of noisy heating patterns, we performed numerical integrations with the linear models discussed in Sect. 6.1, by adding a buoyancy forcing (or heating) term ${\left(\partial \stackrel{\mathrm{\u0303}}{B}/\partial t\right)}_{\text{forcing}}$ to their buoyancy prediction equations. The forcing term is only added to the buoyancy equation at the correction step in the CDgrid model. The noisy forcing pattern, which is created using a random number generator and remains unchanged during the integrations, is confined to a circular region at the center of the domain, and it does not produce a net horizontally averaged heating. The domain and basic state characteristics of the “forced” models are identical to those described earlier in Sect. 6.2. The horizontal grid spacing and the vertical wavenumber are 50 m and n= 320, respectively. We use ${\left(\partial \stackrel{\mathrm{\u0303}}{B}/\partial t\right)}_{\text{forcing}}=\mp \mathrm{1}$ s^{−3} and use no initial perturbations in these simulations. The details are described in the Supplement.
In Fig. 14, we show the time evolutions of the domain maxima and minima of the buoyancy for the Z, C, D and CDgrid simulations. In the Z and Cgrid simulations, the amplitude of buoyancy perturbation does not increase, which implies that the divergence (or vertical velocity) responds efficiently to counter the forcing (heating). In the D and CDgrid simulations, the divergence response is not strong enough to prevent an increase of the amplitudes of the buoyancy perturbations. The D and CDgrid models may need strong explicit diffusion to prevent an excessive accumulation of noise near the source of the noise.
In summary, only the Z grid generates solutions that correspond to those of the continuous system, in which all variables interact with each other on all scales. The frequencies simulated with the Z grid are remarkably close to the true frequencies.
The numerical simulations obtained with the C grid are overall very close to those with the Z grid. However, the C grid generates solutions in which vorticity is decoupled from the divergence for the modes with the SRHS.
The D and CDgrid results are virtually identical to each other. All the variables, i.e., vorticity, divergence, mass and vertical velocity, for the modes with the SRHS are completely decoupled from each other on the D and CD grids.
The A, B and E grids have multiple solutions for all resolutions, and they produce dynamically inert solutions, in which the variables are decoupled from each other, for the SRHS (discussed in the Supplement).
Discretization errors have most often been studied using the shallowwater equations. This is justifiable for an assessment of the errors with the Lamb wave and quasihydrostatic modes. The dispersion of the inertia–gravity modes for the shallowwater system is analogous to that for the Lamb wave, except that the fluid depth is multiplied by $\mathit{\gamma}\equiv \mathrm{1}/\left(\mathrm{1}\mathit{\kappa}\right)$ in the Lamb wave solution. The shallowwater system can also be used to assess the discretization errors of the quasihydrostatic system because in that system the square of the frequency of the inertia–gravity is a linear function of (k^{2}+ℓ^{2}), as it is with the Lamb wave and shallowwater systems. A detailed discussion of the continuous and discrete dispersion equations can be found in the Supplement.
In this respect, the analogy between the shallowwater (also quasihydrostatic) and anelastic systems is weak. The nonhydrostatic effects of the anelastic system, as with the fully compressible, unified and semihydrostatic systems, limit the frequency of the inertia–gravity modes to the Brunt–Väisälä frequency, N, as can be seen from the red curves in Fig. 15. In contrast, the frequency of the quasihydrostatic modes (black lines) and the Lamb wave (dashed green lines) increases without bound with increasing horizontal wavenumber. It appears that the deeper modes are “more nonhydrostatic” than the shallower ones in terms of the separation of the anelastic and quasihydrostatic frequencies in Fig. 15.
The Z and C grids actually perform better with the nonhydrostatic systems than they do with the quasihydrostatic or shallowwater systems. Consider the discrete dispersion equations (Eqs. 18 and 24) for the Z and C grids, respectively. The finitedifference errors are represented by ξ^{2} and η^{2}, which are close to unity for small wavenumbers and approximately equal to 0.4 at the SRHS. This also means that ξ^{2}k^{2}+η^{2}ℓ^{2} is approximately $\mathrm{0.9}\times {k}^{*}$ for ξ≡η at the SRHS. For deep anelastic modes with large composite horizontal wavenumber (k^{*}), the dependence of the frequency on the horizontal wavenumber is negligible, and therefore the error in the frequency due to finite differencing is small for both the Z and C grids. As mentioned earlier, the averaging of the Coriolis term with the C grid does not significantly influence the solution for sufficiently high horizontal wavenumbers.
In conclusion, the Z and C grids perform better with the nonhydrostatic equations than the quasihydrostatic or shallowwater equations, particularly in terms of their ability to resolve deep modes with high horizontal wavenumbers. The D, CD and other grids do not share this benefit from the nonhydrostatic effects because their errors are dominated by the averaging or computational modes at the SRHS. Nonhydrostatic effects do not produce similar benefits for the vertical discretization because the benefits are mostly confined to the deep modes, which are already well resolved in both the nonhydrostatic and quasihydrostatic discrete systems.
We have discussed the horizontal and vertical discretizations of the threedimensional nonhydrostatic linearized anelastic equations on the A, B, C, CD, (DC), D, E and Z horizontal grids, and the L and CP vertical grids, with an emphasis on middlelatitude inertia–gravity waves. The use of a threedimensional nonhydrostatic system instead of the twodimensional shallowwater system allows us to directly assess the performance of the discretization as a function of the vertical wavenumber, vertical resolution or vertical grid spacing.
The Z grid yields the most accurate inertia–gravity wave dispersion among the seven horizontal grids considered in this paper. It has no computational or dynamically inert modes. Although the frequency and group velocity of inertia–gravity modes in the Zgrid solutions are lower than the true ones, the numerical frequency never goes below the inertial frequency and the group velocity never reverses sign.
The C grid produces mixed results for the inertia–gravity solutions. For cloudresolving applications with small horizontal grid spacing represented by d=2 km in our analyses, the accuracy of the physical modes is nearly identical to that of the physical modes with the Z grid. Of course, there is a dynamically inert mode that decouples the divergence and vorticity for the SRHS in the Cgrid solution, but its impact may not be severe because it disperses like a pure gravity mode. However, since the threedimensional enstrophy cascades to the SRHS in nonlinear cloudresolving models, the impact of the dynamically inert mode can be severe, and a parameterized dissipation is needed to dissipate the enstrophy accumulated at the SRHS. For mesoscale applications, i.e., horizontal grid spacings in the range 10 to 25 km, the inertia–gravity modes of the Cgrid solution behave like those of the Zgrid solutions if the vertical wavenumbers are equal to n=320 or smaller for 10 km grid spacings and n=156 or smaller for 25 km grid spacings. With a domain height of z_{T}=80 km, wavenumbers n=320 and n=156 correspond to δz=250 m and δz=512 m, respectively. The inertia–gravity modes in the Cgrid solutions with a typical climate model horizontal grid spacing of 100 km are as accurate as the Zgrid solutions for vertical wavenumbers n=39 and below, which corresponds to δz=2051 m and larger. For vertical wavenumbers higher than n=39, the reversal of group velocity takes place for the modes over a rather wide range of horizontal scales. This might lead to noise or instability.
The inertia–gravity mode solutions with the D and CD grids are almost identical. On these grids, the divergence and the mass (buoyancy, pressure and vertical velocity) are placed at different grid points. The vorticity is placed at the same grid points with the mass. With this staggering, not only are the vorticity and divergence averaged to each other's grid points, but also the divergence is averaged to mass (and vertical velocity) points and the pressure is averaged to divergence points. The result is large errors in the dispersion of the inertia–gravity modes for all vertical scales near the SRHS, for all horizontal grid spacings. At the SRHS, there are dynamically inert modes in all variables. The group velocity reverses and becomes large near the SRHS. In nonlinear models, the D and CD grids may need explicit diffusion to remove the noise in short horizontal scales that results from the computational mode and the reversal of the group velocity.
The dispersions of the inertia–gravity waves with the A, B and E grids are similar. All suffer from the existence of multiple (or nonunique) physical solutions that do not interact with each other (in linear models), as demonstrated by the numerical simulations. The existence of the multiple solutions with these grids also leads to the dynamically inert modes at the SRHS. To avoid the separation of two (or four) solutions with these grids, a horizontal mixing process is needed.
The conclusions of our normalmode analyses of the nonhydrostatic anelastic inertia–gravity waves with the Z, C, D, CD, A, E and B grids for the low horizontal wavenumbers is consistent with earlier the shallowwater analyses by various authors. The Z and C grids produce overall smaller errors with the nonhydrostatic anelastic system than with the quasihydrostatic system, particularly for the deep modes with high horizontal wavenumbers.
Part 2 (Konor and Randall, 2018) discusses the dispersion of the middlelatitude Rossby waves for all the horizontal grids.
Fortran codes that are used to compute and plot the frequencies for the CD grid will be provided by the corresponding author upon request. Related files can also be found in https://doi.org/10.5281/zenodo.1117930.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1117532018supplement.
The authors declare that they have no conflict of interest.
We are grateful to Bill Skamarock for his comments and suggestions to improve
the manuscript. We thank the reviewer, Almut Gaßmann, for her constructive and
helpful comments. We also thank an anonymous reviewer for the helpful
comments and for going through the lengthy derivations. This research was
supported by the National Science Foundation (NSF) under AGS1500187, the US Department of Energy Office of
Science DESC07050 (SciDAC), DESC00016273 (ACME) and DESC00016305
(CMDV).
Edited by: Paul Ullrich
Reviewed by: Almut Gaßmann and one anonymous referee
Adcroft, A., Hill, C., and Marshall, J.: A new treatment of the Coriolis terms in C grid models at both high and low resolutions, Mon. Weather Rev., 127, 1928–1936, 1999.
Adcroft, A., Hallberg, R., Griffies, S., and Dunne, J.: Nextgeneration ocean and ice models at GFDL: MOM6, SI2, and icebergs, available at: http://cosima.org.au/wpcontent/uploads/2016/06/MOM6overviewMay27.pdf, 2016.
Arakawa, A.: Finitedifference methods in climate modeling, in: Physically Based Modelling and Simulation of Climate and Climate Change, Part I, edited by: Schlesinger, M. E., Kluwer Academic Publisher, 79–168, 1988.
Arakawa, A.: A personal perspective on the early years of general circulation modeling at UCLA, General Circulation Model Development, Past, Present, and Future, edited by: Randall, D. A., Academic Press, 70, 1–65, 2000.
Arakawa, A. and Konor, C. S.: Vertical differencing of primitive equations based on the CharneyPhillips grid in hybrid σp vertical coordinates, Mon. Weather Rev., 124, 511–528, 1996.
Arakawa, A. and Konor, C. S.: Unification of the anelastic and quasihydrostatic systems of equations, Mon. Weather Rev., 137, 710–726, 2009.
Arakawa, A. and Lamb, V.: Computational design of the basic dynamical processes of the UCLA general circulation model, in: Methods in Computational Physics, edited by: Chang, J., Academic Press, 17, 173–265, 1977.
Arakawa, A. and Moorthi, S.: Baroclinic instability in vertically discrete systems, J. Atmos. Sci., 45, 1688–1707, 1988.
Charney, J. G. and Phillips, N. A.: Numerical integration of the quasigeostrophic equations for barotropic and simple baroclinic flows, J. Meteorology, 10, 71–99, 1953.
Chen, Q.: On staggering techniques and the nonstaggered Z grid scheme, Numer. Math., 135, 1–21, 2015.
Davies, T., Staniforth, A., Wood, N., and Thuburn, J.: Validity of anelastic and other equation sets as inferred from normalmode analysis, Q. J. Roy. Meteor. Soc., 129, 2761–2775, 2003.
Dubos, T. and Voitus, F.: A semihydrostatic theory of gravitydominated compressible flow, J. Atmos.Sci., 71, 4621–4638, 2014.
Dukowicz, J. K.: Evaluation of various approximations in atmosphere and ocean modeling based on an exact treatment of gravity wave dispersion, Mon. Weather Rev., 141, 4487–4506, 2013.
Gassmann, A.: Inspection of hexagonal and triangular C grid discretizations of the shallow water equations, J. Comput. Phys., 230, 2706–2721, 2011.
Girard, C., Plante, A., Desgagné, M., McTaggardCowan, R., Côté, J., Charron, M., Gravel, S., Lee, V., Patoine, A., Qaddouri, A., Roch, M., Spacek, L., Tanguay, M., Vaillancourt, P. A., and Zadra, A.: Staggered vertical discretization of the Canadian Environmental Multiscale (GEM) model using a coordinate of the loghydrostaticpressure type, Mon. Weather Rev., 142, 1183–1196, 2014.
Hansen, J., Russell, G., Rind, D., Stone, P., Lacis, A., Lebedeff, S., Ruedy, R., and Travis, L.: Efficient threedimensional global models for climate studies: Model I and II, Mon. Weather Rev., 111, 609–662, 1983.
Hollingsworth, A.: A spurious mode in the “Lorenz” arrangement of Φ and T which does not exist in the “CharneyPhillips” arrangement, ECMWF Tech. Memo. No. 211, 12 pp., 1995.
Janjic, Z. I.: Nonlinear advection schemes and energy cascade on semistaggered grids, Mon. Weather Rev., 112, 1234–1245, 1984.
Konor, C. S. and Arakawa, A.: Design of an atmospheric model based on a generalized vertical coordinate, Mon. Weather Rev., 125, 1649–1673, 1997.
Konor, C. S. and Arakawa, A.: Choice of a vertical grid in incorporating condensation heating into an isentropic vertical coordinate model, Mon. Weather Rev., 128, 3901–3910, 2000.
Konor, C. S. and Randall, D. A.: Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – Part 2: Quasigeostrophic Rossby modes, Geosci. Model Dev., 11, 1785–1797, https://doi.org/10.5194/gmd1117852018, 2018.
Lee, J.L. and MacDonald, A. E.: A finitevolume icosahedral shallowwater model on a local coordinate, Mon. Weather Rev., 128, 3901–3910, 2009.
Lin, S.J. and Rood, R. B.: Multidimensional flux form semiLagrangian transport, Mon. Weather Rev., 1237, 1422–1437, 1996.
Lipps, F. B. and Hemler, R. S.: A scale analysis of deep moist convection and some related numerical calculations, J. Atmos. Sci., 39, 2192–2210, 1982.
Lorenz, E. N.: Energy and numerical weather prediction, Tellus, 12, 364–373, 1960.
Mesinger, F. and Arakawa, A.: Numerical methods used in atmospheric models, GARP Publication Series 17, Part I, 64 pp., 1976.
Randall, D. A.: Geostrophic adjustment and the finitedifference shallowwater equations, Mon. Weather Rev., 122, 1371–1377, 1994.
Satoh, M., Matsuno, T., Tomita, H., Nasuno, T., Iga, S., and Muira, H.: Nonhydrostatic icosahedral atmospheric model (NICAM) for global cloud resolving simulations, J. Comput. Phys., 227, 3486–3514, 2008.
Skamarock, W. C.: A linear analysis of the NCAR CCSM finitevolume dynamical core, Mon. Weather Rev., 136, 2112–2119, 2008.
Skamarock, W. C., Klemp, J. B., Duda, M., Fowler, L. D., Parks, S.H., and Ringler, T. D.: A multiscale nonhydrostatic model using centroidal Voronio tesselations and C grid staggering, Mon. Weather Rev., 140, 3090–3105, 2012.
Thuburn, J.: Numerical wave propagation on the hexagonal C grid, J. Comput. Phys., 227, 5826–5858, 2008.
Thuburn, J. and Woolings, T. J.: Vertical discretization for compressible Euler equation atmospheric models giving optimal representation of normal modes, J. Comput. Phys., 203, 386–404, 2004.
Thuburn, J., Ringler, T. D., Skamarock, W. C., and Klemp, J. B.: Numerical representation of geostropichic modes on arbitrarily structured C grids, J. Comput. Phys., 228, 8321–8335, 2009.
Tokioka, T.: Some consideration on vertical differencing, J. Meteor. Soc. Jpn., 56, 89–111, 1978.
Toy, M. and Randall, D. A.: Comment on the article “Vertical discretizations for compressible Euler equation atmospheric models giving optimal representation of normal modes” by Thuburn and Woolings, J. Comput. Phys., 223, 82–88, 2007.
Toy, M. and Randall, D. A.: Design of a nonhydrostatic model based on a generalized vertical coordinate, Mon. Weather Rev., 137, 2305–2330, 2009.
Weller, H., Thuburn, J., and Cotter, C. J.: Computational modes and grid imprinting on five quasiuniform spherical C grids, Mon. Weather Rev., 140, 2734–2755, 2012.
Winninghoff, F. J.: On the adjustment toward a geostrophic balance in a simple primitive equation model with application to the problems of initialization and objective analysis, PhD dissertation, University of California, Los Angeles, 1968.
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 massconserving semiimplicit semiLagrangian discretization of the deepatmosphere global nonhydrostatic equations, Q. J. Roy. Meteor. Soc., 140, 1505–1520, 2014.
 Abstract
 Introduction
 Linearized anelastic equations with an isothermal basic state
 Horizontal discretization of the linear anelastic equations on different grids and discrete dispersion equations
 Vertical discretization of linear anelastic equations on the L and CP grids and discrete dispersion equations
 A summary of the discrete dispersion equations
 Numerical solutions
 The impact of the nonhydrostatic effects on the horizontal grid selection
 Summary and conclusions
 Code and data availability
 Competing interests
 Acknowledgements
 References
 Supplement
 Abstract
 Introduction
 Linearized anelastic equations with an isothermal basic state
 Horizontal discretization of the linear anelastic equations on different grids and discrete dispersion equations
 Vertical discretization of linear anelastic equations on the L and CP grids and discrete dispersion equations
 A summary of the discrete dispersion equations
 Numerical solutions
 The impact of the nonhydrostatic effects on the horizontal grid selection
 Summary and conclusions
 Code and data availability
 Competing interests
 Acknowledgements
 References
 Supplement