- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications

Journal cover
Journal topic
**Geoscientific Model Development**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer-review process
- For authors
- For editors and referees
- EGU publications

- 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

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

**Correspondence**: Celal S. Konor (csk@atmos.colostate.edu)

Abstract

Back to toptop
We have used a normal-mode 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
normal-mode 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 normal-mode analyses with the Z, C, D, A, E and B grids generally confirm the conclusions of previous shallow-water studies for the cyclone-resolving scales (with low horizontal wavenumbers). We conclude that, aided by nonhydrostatic effects, the Z and C grids become overall more accurate for cloud-resolving resolutions (with high horizontal wavenumbers) than for the cyclone-resolving scales.

A companion paper, Part 2, discusses the impacts of the discretization on the
Rossby modes on a midlatitude *β* plane.

How to cite

Back to top
top
How to cite.

Konor, C. S. and Randall, D. A.: Impacts of the horizontal and vertical grids on the numerical solutions of the dynamical equations – Part 1: Nonhydrostatic inertia–gravity modes, Geosci. Model Dev., 11, 1753–1784, https://doi.org/10.5194/gmd-11-1753-2018, 2018.

1 Introduction

Back to toptop
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 finite-differenced 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 shallow-water 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 quasi-hydrostatic 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 three-dimensional anelastic system introduced by Lipps and Hemler (1982) because the curl of its pressure-gradient force vanishes, which gives an important simplicity to our normal-mode 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 middle-latitude 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 three-dimensional nonhydrostatic system to study the impact of the horizontal discretization on the dispersion of the waves has an important advantage over the two-dimensional shallow-water 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 shallow-water 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 normal-mode analyses are generally consistent with the shallow-water analyses discussed by previous authors for the cyclone-resolving scales. Here, we present a three-dimensional view of the dispersion of the modes by including an accurate assessment of the performance of the horizontal grids for the nonhydrostatic cloud-resolving scales.

For large horizontal grid spacings, the results of our analyses are also applicable to the quasi-hydrostatic systems. The normal-mode analyses presented by Arakawa and Konor (2009) show that the frequencies of the medium- and large-scale inertia–gravity waves, with horizontal wave lengths of approximately 60 km and larger, are nearly identical in the anelastic and quasi-hydrostatic systems. The dispersion of midlatitude Rossby waves is also nearly identical in the two systems, with the exception of the ultra-long waves. The only major difference is that the quasi-hydrostatic 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 normal-mode 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 second-order accurate finite-difference schemes for simplicity. The use of higher-order schemes would lead to minor quantitative differences in our results by reducing errors in the well-resolved 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, quasi-hydrostatic and
shallow-water 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.

2 Linearized anelastic equations with an isothermal basic state

Back to toptop
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.

The linearized horizontal momentum equation is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{v}}{\partial t}}=-f\mathit{k}\times \mathit{v}-{\mathrm{\nabla}}_{\mathrm{H}}P,\end{array}$$

where ** v** is the horizontal velocity,

By taking $\mathit{k}\cdot {\mathrm{\nabla}}_{\mathrm{H}}\times \left(\mathrm{1}\right)$, we obtain the linearized vertical vorticity equation as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial {\mathit{\omega}}_{z}}{\partial t}}=-fD,\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial D}{\partial t}}=f{\mathit{\omega}}_{z}-{\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}P,\end{array}$$

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

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial w}{\partial t}}=-\left({\displaystyle \frac{\partial}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)P+B,\end{array}$$

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

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial B}{\partial t}}=-{N}^{\mathrm{2}}w,\end{array}$$

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
three-dimensionally and temporally constant temperature of the isothermal
basic state. The anelastic continuity equation is

$$\begin{array}{}\text{(6)}& D+\left({\displaystyle \frac{\partial}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)w=\mathrm{0},\end{array}$$

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

$$\begin{array}{}\text{(7)}& \left({\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}}+{N}^{\mathrm{2}}\right)D=\left[{\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {z}^{\mathrm{2}}}}-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)}^{\mathrm{2}}\right]{\displaystyle \frac{\partial P}{\partial t}}.\end{array}$$

The three-dimensional 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

$$\begin{array}{ll}\text{(8)}& {\displaystyle}& {\displaystyle}{\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}P+\left[{\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {z}^{\mathrm{2}}}}-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)}^{\mathrm{2}}\right]P={\displaystyle}& {\displaystyle}f{\mathit{\omega}}_{z}+\left({\displaystyle \frac{\partial}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)B.\end{array}$$

Although the elliptic equation (Eq. 8) is not needed for the normal-mode 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 quasi-static, given by

$$\begin{array}{}\text{(9a)}& \mathrm{0}=f{\mathit{\omega}}_{z}-{\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}P\end{array}$$

and

$$\begin{array}{}\text{(9b)}& \mathrm{0}=-\left({\displaystyle \frac{\partial}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)P+B.\end{array}$$

By using Eq. (9a) in Eq. (10), we find that

$$\begin{array}{}\text{(9c)}& \mathrm{0}=-f\left({\displaystyle \frac{\partial}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right){\mathit{\omega}}_{z}+{\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}B\end{array}$$

is also satisfied for the steady balanced state. The majority of the discrete systems that we will discuss have steady balanced-state solutions corresponding to Eq. (9a)–(9c). However, we found some exceptions with the CD-grid discretization, depending on the time-integration 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
middle-latitude *f* plane, we seek solutions of Eqs. (2), (3)
and (7) of the form

$$\begin{array}{}\text{(10)}& \mathit{\varphi}\left(x,y,z,t\right)=\text{Re}\left\{\widehat{\mathrm{\Phi}}{e}^{\underset{\sim}{i}\left(kx+\mathrm{\ell}y+mz-\mathit{\nu}t\right)}\right\},\end{array}$$

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

$$\begin{array}{}\text{(11)}& m\equiv {\displaystyle \frac{\mathit{\pi}n}{{z}_{\text{T}}}}\phantom{\rule{0.25em}{0ex}}\text{for}\phantom{\rule{0.25em}{0ex}}n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots}\end{array}$$

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.

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

$$\begin{array}{}\text{(12)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{N}^{\mathrm{2}}\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}}.\end{array}$$

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 steady-state 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.

3 Horizontal discretization of the linear anelastic equations on different grids and discrete dispersion equations

Back to toptop
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 non-unique (or multiple) physical solutions that do not interact in the linear system. These non-unique 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 normal-mode 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 normal-mode analysis similar to the one used for the continuous system. Fig. 1 shows portions of the horizontal grids discussed in this paper.

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 shallow-water 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.

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

$$\begin{array}{ll}\text{(13)}& {\displaystyle}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{i,j}}{\partial t}}=-f{D}_{i,j},\text{(14)}& {\displaystyle}& {\displaystyle \frac{\partial {D}_{i,j}}{\partial t}}={\displaystyle}& {\displaystyle}f{\left({\mathit{\omega}}_{z}\right)}_{i,j}-{\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{1},j}+{P}_{i-\mathrm{1},j}+{P}_{i,j+\mathrm{1}}+{P}_{i,j-\mathrm{1}}-\mathrm{4}{P}_{i,j}\right),\\ \text{(15)}& {\displaystyle}& {\displaystyle}\left({\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}}+{N}^{\mathrm{2}}\right){D}_{i,j}=\left[{\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {z}^{\mathrm{2}}}}-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)}^{\mathrm{2}}\right]{\displaystyle \frac{\partial {P}_{i,j}}{\partial t}}.\end{array}$$

By using

$$\begin{array}{}\text{(16)}& {\mathit{\varphi}}_{i,j}\left(z,t\right)=\text{Re}\left\{\widehat{\mathrm{\Phi}}{e}^{\underset{\sim}{i}\left(kdi+\mathrm{\ell}dj+mz-\mathit{\nu}t\right)}\right\}\end{array}$$

in Eqs. (14)–(16), and seeking nontrivial solutions, we can obtain the discrete dispersion relation as

$$\begin{array}{}\text{(17)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{N}^{\mathrm{2}}\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

where the factors

$$\begin{array}{}\text{(18)}& \mathit{\xi}\equiv {\displaystyle \frac{\mathrm{sin}\left(\frac{\mathrm{1}}{\mathrm{2}}kd\right)}{\frac{\mathrm{1}}{\mathrm{2}}kd}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\mathit{\eta}\equiv {\displaystyle \frac{\mathrm{sin}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\ell}d\right)}{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\ell}d}}\end{array}$$

arise from the finite-difference 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 higher-order 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 lower-tropospheric 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 cloud-resolving 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 Z-grid 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 Z-grid solutions maintain positive (or zero) group velocity, but near the SRHS the Z-grid 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.

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

We define the discrete vorticity and divergence in terms of the components of the horizontal velocity on the C grid (Fig. 1b) as

$$\begin{array}{ll}\text{(19a)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({u}_{i+\mathrm{1}/\mathrm{2},j}-{u}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}}+{v}_{i+\mathrm{1},j+\mathrm{1}/\mathrm{2}}-{v}_{i,j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{}\text{(19b)}& {D}_{i,j}\equiv {\displaystyle \frac{\mathrm{1}}{d}}\left({u}_{i+\mathrm{1}/\mathrm{2},j}-{u}_{i-\mathrm{1}/\mathrm{2},j}+{v}_{i,j+\mathrm{1}/\mathrm{2}}-{v}_{i,j-\mathrm{1}/\mathrm{2}}\right),\end{array}$$

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

$$\begin{array}{ll}\text{(20)}& {\displaystyle}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}}{\partial t}}={\displaystyle}& {\displaystyle}-f{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i,j}+{D}_{i+\mathrm{1},j}+{D}_{i,j+\mathrm{1}}+{D}_{i+\mathrm{1},j+\mathrm{1}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(21)}& {\displaystyle \frac{\partial {D}_{i,j}}{\partial t}}=& {\displaystyle}\phantom{\rule{0.25em}{0ex}}f{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.{\displaystyle}& {\displaystyle}\left.+{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right]\\ {\displaystyle}& {\displaystyle}-{\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{1},j}+{P}_{i-\mathrm{1},j}+{P}_{i,j+\mathrm{1}}+{P}_{i,j-\mathrm{1}}-\mathrm{4}{P}_{i,j}\right).\end{array}$$

By using Eq. (17) and

$$\begin{array}{ll}\text{(22)}& {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\left(z,t\right)={\displaystyle}& {\displaystyle}\text{Re}\left\{\widehat{\mathrm{\Phi}}{e}^{\underset{\sim}{i}\left[k\left(i+\frac{\mathrm{1}}{\mathrm{2}}\right)d+\mathrm{\ell}\left(j+\frac{\mathrm{1}}{\mathrm{2}}\right)d+mz-\mathit{\nu}t\right]}\right\}\end{array}$$

in Eqs. (21), (22) and (16), we obtain the discrete dispersion equation for the inertia–gravity waves on the C grid as

$$\begin{array}{}\text{(23)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{N}^{\mathrm{2}}\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+{\mathit{\mu}}^{\mathrm{2}}{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

where *ξ* and *η* are defined by Eq. (18) and

$$\begin{array}{}\text{(24)}& \mathit{\mu}\equiv \mathrm{cos}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}kd\right)\mathrm{cos}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\ell}d\right).\end{array}$$

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 higher-order Laplacian would slightly improve the errors for
the well-resolved 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 C-grid frequency solution is almost identical to
the Z-grid solution. The solutions differ for the modes with vertical integer
wavenumber *n*=1947, for which the C-grid solutions generate negative group
velocities ($\partial \left|\mathit{\nu}\right|/\partial {k}^{*}<\mathrm{0})$. This means
that the nonhydrostatic C-grid 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 C-grid 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 Z-grid 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.

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 finite-difference 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.

We define the discrete vorticity and divergence from the components of the horizontal velocity on the D grid (Fig. 1c) as

$$\begin{array}{ll}\text{(25a)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i,j}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({v}_{i-\mathrm{1}/\mathrm{2},j}-{v}_{i+\mathrm{1}/\mathrm{2},j}+{u}_{i,j-\mathrm{1}/\mathrm{2}}-{u}_{i,j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(25b)}& {\displaystyle}& {\displaystyle}{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({u}_{i+\mathrm{1},j+\mathrm{1}/\mathrm{2}}-{u}_{i,j+\mathrm{1}/\mathrm{2}}+{v}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}}-{v}_{i+\mathrm{1}/\mathrm{2},j}\right),\end{array}$$

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

$$\begin{array}{ll}\text{(26)}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{i,j}}{\partial t}}=& {\displaystyle}\phantom{\rule{0.25em}{0ex}}-f{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{D}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.{\displaystyle}& {\displaystyle}\left.+{D}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(27)}& {\displaystyle}& {\displaystyle \frac{\partial}{\partial t}}{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}={\displaystyle}& {\displaystyle}f{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j+\mathrm{1}}+{\left({\mathit{\omega}}_{z}\right)}_{i,j+\mathrm{1}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j}+{\left({\mathit{\omega}}_{z}\right)}_{i,j}\right]\\ {\displaystyle}& {\displaystyle}-{\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({\stackrel{\mathrm{\u203e}}{P}}_{i+\mathrm{3}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\stackrel{\mathrm{\u203e}}{P}}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{\stackrel{\mathrm{\u203e}}{P}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{3}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}\left.+{\stackrel{\mathrm{\u203e}}{P}}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-\mathrm{4}{\stackrel{\mathrm{\u203e}}{P}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

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

$$\begin{array}{ll}\text{(28)}& {\displaystyle}\left({\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}}+{N}^{\mathrm{2}}\right)& {\displaystyle}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{D}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.\right.{\displaystyle}& {\displaystyle}\left.+{D}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right)]\\ {\displaystyle}& {\displaystyle}-\left[{\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {z}^{\mathrm{2}}}}-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right)}^{\mathrm{2}}\right]{\displaystyle \frac{\partial}{\partial t}}{P}_{i,j}=\mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(29)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{\mathit{\mu}}^{\mathrm{2}}{N}^{\mathrm{2}}\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+{\mathit{\mu}}^{\mathrm{2}}{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{{\mathit{\mu}}^{\mathrm{2}}\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

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 D-grid 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 D-grid 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.

The CD grid proposed by Lin and Rood (1997) was first used in a flux-form semi-Lagrangian dynamical core. It is now also used in Geophysical Fluid Dynamics Laboratory (GFDL)'s Finite-Volume Cubed-Sphere Dynamical Core (FV3) and in the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM). FV3 will also be used in the next-generation numerical prediction model of the US National Centers for Environmental Prediction. The CD grid is built on a predictor–corrector time-integration scheme, in which the C- and D-grid 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 shallow-water equations discretized on the CD grid for the time-continuous 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 normal-mode analysis using second-order accurate differences. As mentioned earlier, the main conclusions of our analysis apply without change to systems discretized with higher-order 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

$$\begin{array}{}\text{(30a)}& {\mathit{\varphi}}^{(*)}={\mathit{\varphi}}^{\left(n\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}{F}^{\left(n\right)},\end{array}$$

and

$$\begin{array}{}\text{(30b)}& {\mathit{\varphi}}^{(n+\mathrm{1})}={\mathit{\varphi}}^{\left(n\right)}+\mathit{\tau}{G}^{(*)},\end{array}$$

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:

$$\begin{array}{ll}\text{(31a)}& {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{i,j}^{(*)}\left(z\right)\equiv \text{Re}\left\{{\widehat{\mathrm{\Phi}}}^{(*)}{e}^{\underset{\sim}{i}\left(kdi+\mathrm{\ell}dj+mz\right)}\right\}{\displaystyle}& {\displaystyle}\text{and}\\ {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{(*)}\left(z\right)=\text{Re}\left\{{\widehat{\mathrm{\Phi}}}^{(*)}{e}^{\underset{\sim}{i}\left[k\left(i+\mathrm{1}/\mathrm{2}\right)d+\mathrm{\ell}\left(j+\mathrm{1}/\mathrm{2}\right)d+mz\right]}\right\},\end{array}$$

and

$$\begin{array}{ll}\text{(31b)}& {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{i,j}^{\left(n\right)}\left(z\right)\equiv \text{Re}\left\{{\widehat{\mathrm{\Phi}}}^{\left(n\right)}{e}^{\underset{\sim}{i}\left(kdi+\mathrm{\ell}dj+mz\right)}\right\}{\displaystyle}& {\displaystyle}\text{and}\\ {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n\right)}\left(z\right)=\text{Re}\left\{{\widehat{\mathrm{\Phi}}}^{\left(n\right)}{e}^{\underset{\sim}{i}\left[k\left(i+\mathrm{1}/\mathrm{2}\right)d+\mathrm{\ell}\left(j+\mathrm{1}/\mathrm{2}\right)d+mz\right]}\right\},\end{array}$$

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:

$$\begin{array}{}\text{(32)}& {\displaystyle}& {\displaystyle}{\widehat{\mathit{\omega}}}_{z}^{(*)}=\mathit{\mu}{\widehat{\mathit{\omega}}}_{z}^{\left(n\right)}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}f{\widehat{D}}^{\left(n\right)},\text{(33)}& {\displaystyle}& {\displaystyle}{\widehat{D}}^{(*)}=\mathit{\mu}{\widehat{D}}^{\left(n\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}\left(f{\widehat{\mathit{\omega}}}_{z}^{\left(n\right)}+{L}^{\mathrm{2}}{\widehat{P}}_{\text{C}}\right),\text{(34)}& {\displaystyle}& {\displaystyle}{\widehat{w}}^{(*)}={\widehat{w}}^{\left(n\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}\left[-\left(\underset{\sim}{i}m+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){\widehat{P}}_{\text{C}}+{\widehat{B}}^{\left(n\right)}\right],\text{(35)}& {\displaystyle}& {\displaystyle}{\widehat{B}}^{(*)}={\widehat{B}}^{\left(n\right)}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}{N}^{\mathrm{2}}{\widehat{w}}^{\left(n\right)},\text{(36)}& {\displaystyle}& {\displaystyle}{\widehat{D}}^{(*)}+\left(\underset{\sim}{i}m-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){\widehat{w}}^{(*)}=\mathrm{0},\text{(37)}& {\displaystyle}& {\displaystyle}\text{In(33),}\phantom{\rule{0.25em}{0ex}}{L}^{\mathrm{2}}\equiv {\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(38)}& {\displaystyle}& {\displaystyle}{\widehat{\mathit{\omega}}}_{z}^{\left(n+\mathrm{1}\right)}={\widehat{\mathit{\omega}}}_{z}^{\left(n\right)}-\mathit{\tau}f{\widehat{D}}^{(*)},\text{(39)}& {\displaystyle}& {\displaystyle}{\widehat{D}}^{\left(n+\mathrm{1}\right)}={\widehat{D}}^{\left(n\right)}+\mathit{\tau}\left(f{\widehat{\mathit{\omega}}}_{z}^{(*)}+\mathit{\mu}{L}^{\mathrm{2}}{\widehat{P}}_{\text{D}}\right),\text{(40)}& {\displaystyle}& {\displaystyle}{\widehat{w}}^{\left(n+\mathrm{1}\right)}={\widehat{w}}^{\left(n\right)}+\mathit{\tau}\left[-\left(\underset{\sim}{i}m+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){\widehat{P}}_{\text{D}}+{\widehat{B}}^{(*)}\right],\text{(41)}& {\displaystyle}& {\displaystyle}{\widehat{B}}^{\left(n+\mathrm{1}\right)}={\widehat{B}}^{\left(n\right)}-\mathit{\tau}{N}^{\mathrm{2}}{\widehat{w}}^{(*)},\text{(42)}& {\displaystyle}& {\displaystyle}\mathit{\mu}{\widehat{D}}^{\left(n+\mathrm{1}\right)}+\left(\underset{\sim}{i}m-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){\widehat{w}}^{\left(n+\mathrm{1}\right)}=\mathrm{0},\end{array}$$

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 time-integration 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 time-integration 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 correction-step equations by using the equations from the predictor step. Then we use

$$\begin{array}{}\text{(43)}& {\widehat{\mathrm{\Phi}}}^{\left(n+\mathrm{1}\right)}={e}^{-\underset{\sim}{i}\mathit{\nu}\mathit{\tau}}{\widehat{\mathrm{\Phi}}}^{\left(n\right)},\end{array}$$

where $\widehat{\mathrm{\Phi}}$ represents ${\widehat{\mathit{\omega}}}_{z}$, $\widehat{D}$ and $\widehat{B}$, to obtain the discrete dispersion relation for complex frequencies as

$$\begin{array}{ll}\text{(44)}& {\displaystyle}& {\displaystyle}{\mathit{\mu}}^{\mathrm{2}}\left[{\left({e}^{-\underset{\sim}{i}\mathit{\nu}\mathit{\tau}}-{\mathit{\sigma}}_{N}\right)}^{\mathrm{2}}+{\mathit{\tau}}^{\mathrm{2}}{N}^{\mathrm{2}}\right]{L}^{\mathrm{2}}{\displaystyle}& {\displaystyle}+\left[{\left({e}^{-\underset{\sim}{i}\mathit{\nu}\mathit{\tau}}-{\mathit{\sigma}}_{f}\right)}^{\mathrm{2}}+{\mathit{\mu}}^{\mathrm{2}}{\mathit{\tau}}^{\mathrm{2}}{f}^{\mathrm{2}}\right]{\mathit{\sigma}}_{m}^{\mathrm{2}}=\mathrm{0},\end{array}$$

where we use the following definitions to shorten the equation:

$$\begin{array}{ll}\text{(45)}& {\displaystyle}& {\displaystyle}{\mathit{\sigma}}_{N}\equiv \mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\tau}}^{\mathrm{2}}{N}^{\mathrm{2}},{\displaystyle}& {\displaystyle}{\mathit{\sigma}}_{f}\equiv \mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\tau}}^{\mathrm{2}}{f}^{\mathrm{2}}\\ {\displaystyle}& {\displaystyle}\text{and}\phantom{\rule{0.25em}{0ex}}{\mathit{\sigma}}_{m}^{\mathrm{2}}\equiv {m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}.\end{array}$$

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

$$\begin{array}{ll}\text{(46a)}& {\displaystyle}& {\displaystyle}{e}^{\mathrm{2}{\mathit{\nu}}_{i}\mathit{\tau}}\left({\mathit{\mu}}^{\mathrm{2}}{L}^{\mathrm{2}}+{\mathit{\sigma}}_{m}^{\mathrm{2}}\right)\mathrm{cos}\left(\mathrm{2}{\mathit{\nu}}_{r}\mathit{\tau}\right){\displaystyle}& {\displaystyle}-\mathrm{2}{e}^{{\mathit{\nu}}_{i}\mathit{\tau}}\left({\mathit{\mu}}^{\mathrm{2}}{\mathit{\sigma}}_{N}{L}^{\mathrm{2}}+{\mathit{\sigma}}_{f}{\mathit{\sigma}}_{m}^{\mathrm{2}}\right)\mathrm{cos}\left({\mathit{\nu}}_{r}\mathit{\tau}\right)\\ {\displaystyle}& {\displaystyle}+\left({\mathit{\sigma}}_{N}^{\mathrm{2}}+{\mathit{\tau}}^{\mathrm{2}}{N}^{\mathrm{2}}\right){\mathit{\mu}}^{\mathrm{2}}{L}^{\mathrm{2}}+\left({\mathit{\sigma}}_{f}^{\mathrm{2}}+{\mathit{\tau}}^{\mathrm{2}}{\mathit{\mu}}^{\mathrm{2}}{f}^{\mathrm{2}}\right){\mathit{\sigma}}_{m}^{\mathrm{2}}=\mathrm{0},\end{array}$$

and

$$\begin{array}{}\text{(46b)}& {\displaystyle}{e}^{{\mathit{\nu}}_{i}\mathit{\tau}}={\displaystyle \frac{\mathrm{2}\left({\mathit{\mu}}^{\mathrm{2}}{\mathit{\sigma}}_{N}{L}^{\mathrm{2}}+{\mathit{\sigma}}_{f}{\mathit{\sigma}}_{m}^{\mathrm{2}}\right)}{\left({\mathit{\mu}}^{\mathrm{2}}{L}^{\mathrm{2}}+{\mathit{\sigma}}_{m}^{\mathrm{2}}\right)}}{\displaystyle \frac{\mathrm{sin}\left({\mathit{\nu}}_{r}\mathit{\tau}\right)}{\mathrm{sin}\left(\mathrm{2}{\mathit{\nu}}_{r}\mathit{\tau}\right)}},\end{array}$$

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 near-neutral 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 quasi-static 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

$$\begin{array}{ll}\text{(47)}& {\displaystyle}& {\displaystyle}f{\widehat{\mathit{\omega}}}_{z}=-{L}^{\mathrm{2}}\widehat{P},{\displaystyle}& {\displaystyle}f\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){\widehat{\mathit{\omega}}}_{z}=\left(\underset{\sim}{i}m-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){L}^{\mathrm{2}}\widehat{B}\\ {\displaystyle}& {\displaystyle}\text{and}\phantom{\rule{0.25em}{0ex}}\left(\underset{\sim}{i}m+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right)\widehat{P}=\widehat{B}.\end{array}$$

These are the discrete equivalents of the continuous relations given by Eq. (9a)–(11c). We conclude that the balanced-state solution is maintained by Scheme I.

The results presented in this section, as inferred from our normal-mode analyses, demonstrate that the dispersion of the inertia–gravity waves produced on the CD grid with various time-integration 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 time-integration choices (not discussed here), we could not uniquely determine a steady balanced-state 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 D-grid case can be applied to the
CD grid as well.

Our normal-mode 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 C-grid-like discrete dispersion relation for the DC-grid 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.

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 non-interacting 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 Flow-Following 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).

We define the discrete vorticity and divergence from the components of the horizontal velocity on the A grid shown in Fig. 1e as

$$\begin{array}{}\text{(48a)}& {\left({\mathit{\omega}}_{z}\right)}_{i,j}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({v}_{i+\mathrm{1},j}-{v}_{i-\mathrm{1},j}-{u}_{i,j+\mathrm{1}}+{u}_{i,j-\mathrm{1}}\right),\end{array}$$

and

$$\begin{array}{}\text{(48b)}& {D}_{i,j}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({u}_{i+\mathrm{1},j}-{u}_{i-\mathrm{1},j}+{v}_{i,j+\mathrm{1}}-{v}_{i,j-\mathrm{1}}\right),\end{array}$$

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

$$\begin{array}{}\text{(49)}& {\displaystyle \frac{\partial}{\partial t}}{\left({\mathit{\omega}}_{z}\right)}_{i,j}=-f{D}_{i,j},\end{array}$$

and

$$\begin{array}{ll}\text{(50)}& {\displaystyle}& {\displaystyle \frac{\partial}{\partial t}}{D}_{i,j}={\displaystyle}& {\displaystyle}f{\left({\mathit{\omega}}_{z}\right)}_{i,j}-{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{2},j}+{P}_{i-\mathrm{2},j}+{P}_{i,j+\mathrm{2}}+{P}_{i,j-\mathrm{2}}-\mathrm{4}{P}_{i,j}\right).\end{array}$$

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

$$\begin{array}{}\text{(51)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{N}^{\mathrm{2}}\left({\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({\stackrel{\mathrm{\u0303}}{\mathit{\xi}}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\stackrel{\mathrm{\u0303}}{\mathit{\eta}}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

where

$$\begin{array}{}\text{(52)}& {\displaystyle}\stackrel{\mathrm{\u0303}}{\mathit{\xi}}\equiv {\displaystyle \frac{\mathrm{sin}\left(kd\right)}{kd}}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}\stackrel{\mathrm{\u0303}}{\mathit{\eta}}\equiv {\displaystyle \frac{\mathrm{sin}\left(\mathrm{\ell}d\right)}{\mathrm{\ell}d}}.\end{array}$$

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
normal-mode analysis, the A grid generates multiple, linearly decoupled (or
non-unique) 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.

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
non-interacting 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.

We first define the discrete vorticity and divergence from the components of
the horizontal velocity for the integer (*i*,*j*) and half-integer
($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 half-integer points as follows:

$$\begin{array}{ll}\text{(53a)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i,j}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({v}_{i+\mathrm{1}/\mathrm{2},j}-{u}_{i,j+\mathrm{1}/\mathrm{2}}-{v}_{i-\mathrm{1}/\mathrm{2},j}+{u}_{i,j-\mathrm{1}/\mathrm{2}}\right),\\ \text{(53b)}& {\displaystyle}& {\displaystyle}{D}_{i,j}\equiv {\displaystyle \frac{\mathrm{1}}{d}}\left({u}_{i+\mathrm{1}/\mathrm{2},j}-{v}_{i,j-\mathrm{1}/\mathrm{2}}+{v}_{i,j+\mathrm{1}/\mathrm{2}}-{u}_{i-\mathrm{1}/\mathrm{2},j}\right),\end{array}$$

$$\begin{array}{ll}\text{(54)}& {\displaystyle}& {\displaystyle \frac{\partial}{\partial t}}({\mathit{\omega}}_{z}{)}_{i,j}=-f{D}_{i,j}\text{(55)}& {\displaystyle}& {\displaystyle \frac{\partial {D}_{i,j}}{\partial t}}={\displaystyle}& {\displaystyle}f{\left({\mathit{\omega}}_{z}\right)}_{i,j}-{\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{1},j}+{P}_{i-\mathrm{1},j}+{P}_{i,j+\mathrm{1}}+{P}_{i,j-\mathrm{1}}-\mathrm{4}{P}_{i,j}\right),\\ \text{(56)}& {\displaystyle}& {\displaystyle}\left({\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}}+{N}^{\mathrm{2}}\right){D}_{i,j}=-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){\displaystyle \frac{\partial {P}_{i,j}}{\partial t}},\end{array}$$

$$\begin{array}{ll}\text{(57a)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({v}_{i+\mathrm{1},j+\mathrm{1}/\mathrm{2}}-{u}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}}-{v}_{i,j+\mathrm{1}/\mathrm{2}}+{u}_{i+\mathrm{1}/\mathrm{2},j}\right),\\ \text{(57b)}& {\displaystyle}& {\displaystyle}{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{d}}\left({u}_{i+\mathrm{1},j+\mathrm{1}/\mathrm{2}}-{v}_{i+\mathrm{1}/\mathrm{2},j}+{v}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}}-{u}_{i,j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

$$\begin{array}{}\text{(58)}& {\displaystyle}{\displaystyle \frac{\partial}{\partial t}}{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}=-f{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}},\end{array}$$

$$\begin{array}{ll}\text{(59)}& {\displaystyle \frac{\partial {D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}}{\partial t}}=& {\displaystyle}\phantom{\rule{0.25em}{0ex}}f{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}{\displaystyle}& {\displaystyle}-{\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{3}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{P}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}+{P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{3}/\mathrm{2}}+{P}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\\ {\displaystyle}& {\displaystyle}\left.-\mathrm{4}{P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

$$\begin{array}{ll}\text{(60)}& {\displaystyle}& {\displaystyle}\left({\displaystyle \frac{{\partial}^{\mathrm{2}}}{\partial {t}^{\mathrm{2}}}}+{N}^{\mathrm{2}}\right){D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}={\displaystyle}& {\displaystyle}-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){\displaystyle \frac{\partial {P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}}{\partial t}}.\end{array}$$

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 half-integer points as

$$\begin{array}{}\text{(61)}& {\mathit{\nu}}^{\mathrm{2}}={\displaystyle \frac{{N}^{\mathrm{2}}\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({\mathit{\xi}}^{\mathrm{2}}{k}^{\mathrm{2}}+{\mathit{\eta}}^{\mathrm{2}}{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

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
half-integer 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 half-integer points (Eqs. 58–60). The discrete equations
for the integer points use only information from integer points, and those
for the half-integer points use only information from the half-integer
points. Of course, this means that the solution includes computational modes
arising from multiple (or non-unique) solutions. Further discussion of the
computational-mode on the E grid is given in Sect. 3.8.

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

We define the discrete versions of the vorticity and divergence on the B grid shown in Fig. 1g as

$$\begin{array}{ll}\text{(62a)}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i,j}\equiv & {\displaystyle}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({v}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{v}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.{\displaystyle}& {\displaystyle}\left.-{v}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-{v}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right)\\ {\displaystyle}& {\displaystyle}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({u}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{u}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}\left.-{u}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}-{u}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(62b)}& {\displaystyle}{D}_{i,j}\equiv & {\displaystyle}\phantom{\rule{0.25em}{0ex}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({u}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{u}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.{\displaystyle}& {\displaystyle}\left.-{u}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-{u}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right)\\ {\displaystyle}& {\displaystyle}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}d}}\left({v}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{v}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}\left.-{v}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}-{v}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right),\end{array}$$

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

$$\begin{array}{ll}\text{(63)}& {\displaystyle \frac{\partial}{\partial t}}{D}_{i,j}=& {\displaystyle}f{\left({\mathit{\omega}}_{z}\right)}_{i,j}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{1},j+\mathrm{1}}+{P}_{i+\mathrm{1},j-\mathrm{1}}\right.{\displaystyle}& {\displaystyle}\left.+{P}_{i-\mathrm{1},j+\mathrm{1}}+{P}_{i-\mathrm{1},j-\mathrm{1}}-\mathrm{4}{P}_{i,j}\right).\end{array}$$

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 Z-grid equation (Eq. 16) and E-grid 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

$$\begin{array}{ll}\text{(64)}& {\displaystyle}& {\displaystyle}{\mathit{\nu}}^{\mathrm{2}}={\displaystyle}& {\displaystyle \frac{{N}^{\mathrm{2}}\left({\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}}\right)+{f}^{\mathrm{2}}\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({\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}}\right)+\left({m}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

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.

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 right-hand 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 non-interacting (non-unique) solutions as indicated with different colors in Fig. 1e, and the B and E grids can produce two independent and non-interacting (non-unique) 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 small-scale physical modes.

The Z, C, D and CD grids do not have independent and non-interacting 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.

The dispersion relations of the A-, E- and B-grid 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 normal-mode 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 normal-modes of the shallow-water equations as described by many authors, including Arakawa and Lamb (1977) and Randall (1994).

4 Vertical discretization of linear anelastic equations on the L and CP grids and discrete dispersion equations

Back to toptop
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 L-grid 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 normal-mode 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.

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

$$\begin{array}{}\text{(65)}& {\displaystyle}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{k}}{\partial t}}=-f{D}_{k},\text{(66)}& {\displaystyle}& {\displaystyle \frac{\partial {D}_{k}}{\partial t}}=f{\left({\mathit{\omega}}_{z}\right)}_{k}-{\mathrm{\nabla}}_{\mathrm{H}}^{\mathrm{2}}{P}_{k},\end{array}$$

$$\begin{array}{ll}\text{(67)}& {\displaystyle \frac{\partial {w}_{k+\mathrm{1}/\mathrm{2}}}{\partial t}}=& {\displaystyle}-{\displaystyle \frac{{P}_{k+\mathrm{1}}-{P}_{k}}{\mathit{\delta}z}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right){\displaystyle \frac{{P}_{k+\mathrm{1}}+{P}_{k}}{\mathrm{2}}}{\displaystyle}& {\displaystyle}+{\displaystyle \frac{{B}_{k}+{B}_{k+\mathrm{1}}}{\mathrm{2}}},\end{array}$$

$$\begin{array}{}\text{(68)}& {\displaystyle}{\displaystyle \frac{\partial {B}_{k}}{\partial t}}=-{N}^{\mathrm{2}}{\displaystyle \frac{{w}_{k+\mathrm{1}/\mathrm{2}}+{w}_{k-\mathrm{1}/\mathrm{2}}}{\mathrm{2}}},\end{array}$$

and

$$\begin{array}{ll}\text{(69)}& {\displaystyle}{D}_{k}& {\displaystyle}+{\displaystyle \frac{{w}_{k+\mathrm{1}/\mathrm{2}}-{w}_{k-\mathrm{1}/\mathrm{2}}}{\mathit{\delta}z}}{\displaystyle}& {\displaystyle}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right){\displaystyle \frac{{w}_{k+\mathrm{1}/\mathrm{2}}+{w}_{k-\mathrm{1}/\mathrm{2}}}{\mathrm{2}}}=\mathrm{0},\end{array}$$

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

$$\begin{array}{ll}\text{(70)}& {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{k}\left(x,y,t\right)=\text{Re}\left\{\widehat{\mathrm{\Phi}}{e}^{\underset{\sim}{i}\left(kx+\mathrm{\ell}y+m\mathit{\delta}zk-\mathit{\nu}t\right)}\right\}{\displaystyle}& {\displaystyle}\text{and}\\ {\displaystyle}& {\displaystyle}{\mathit{\varphi}}_{k+\mathrm{1}/\mathrm{2}}\left(x,y,t\right)=\text{Re}\left\{\widehat{\mathrm{\Phi}}{e}^{\underset{\sim}{i}\left[kx+\mathrm{\ell}y+m\mathit{\delta}z\left(k+\mathrm{1}/\mathrm{2}\right)-\mathit{\nu}t\right]}\right\}\end{array}$$

in Eqs. (66)–(70), we obtain the discrete dispersion relation for the inertia–gravity modes on the L grid as

$$\begin{array}{}\text{(71)}& {\mathit{\nu}}^{\mathrm{2}}=\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{{\mathit{\mu}}_{z}^{\mathrm{2}}{N}^{\mathrm{2}}\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({\mathit{\zeta}}^{\mathrm{2}}{m}^{\mathrm{2}}+{\mathit{\mu}}_{z}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({\mathit{\zeta}}^{\mathrm{2}}{m}^{\mathrm{2}}+{\mathit{\mu}}_{z}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

where

$$\begin{array}{}\text{(72)}& \mathit{\zeta}\equiv {\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{1}}{\mathrm{2}}m\mathit{\delta}z}}\mathrm{sin}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m\mathit{\delta}z\right)\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}{\mathit{\mu}}_{z}\equiv \mathrm{cos}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}m\mathit{\delta}z\right).\end{array}$$

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∕4*H*^{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.

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

$$\begin{array}{ll}\text{(73)}& {\displaystyle}& {\displaystyle \frac{\partial {w}_{k+\mathrm{1}/\mathrm{2}}}{\partial t}}={\displaystyle}& {\displaystyle}-{\displaystyle \frac{{P}_{k+\mathrm{1}}-{P}_{k}}{\mathit{\delta}z}}+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\rho}}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\rho}}_{\mathrm{0}}}{\partial z}}\right){\displaystyle \frac{{P}_{k+\mathrm{1}}+{P}_{k}}{\mathrm{2}}}+{B}_{k+\mathrm{1}/\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(74)}& {\displaystyle \frac{\partial {B}_{k+\mathrm{1}/\mathrm{2}}}{\partial t}}=-{N}^{\mathrm{2}}{w}_{k+\mathrm{1}/\mathrm{2}},\end{array}$$

respectively. Following the procedure for the L-grid case, we obtain the discrete dispersion relation for the inertia–gravity modes as

$$\begin{array}{}\text{(75)}& {\mathit{\nu}}^{\mathrm{2}}=\phantom{\rule{0.125em}{0ex}}{\displaystyle \frac{{N}^{\mathrm{2}}\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+{f}^{\mathrm{2}}\left({\mathit{\zeta}}^{\mathrm{2}}{m}^{\mathrm{2}}+{\mathit{\mu}}_{z}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}{\left({k}^{\mathrm{2}}+{\mathrm{\ell}}^{\mathrm{2}}\right)+\left({\mathit{\zeta}}^{\mathrm{2}}{m}^{\mathrm{2}}+{\mathit{\mu}}_{z}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}\right)}},\end{array}$$

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 CP-grid 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
L-grid 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 CP-grid 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 high-resolution nonhydrostatic applications than in low-resolution quasi-hydrostatic ones.

5 A summary of the discrete dispersion equations

Back to toptop
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.

6 Numerical solutions

Back to toptop
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 middle-latitude
*f* plane. Our purpose is to confirm the results of the normal-mode analyses
of the discrete equations and to investigate issues that are not completely
revealed by the normal-mode 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.

We write the horizontally discrete version equations (Eqs. 2–6 and 8) on the C and D grids as

$$\begin{array}{}\text{(76)}& {\displaystyle}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}}{\partial t}}=-f{\stackrel{\mathrm{\u203e}}{D}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}},\text{(77)}& {\displaystyle}& {\displaystyle \frac{\partial {D}_{i,j}}{\partial t}}=f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}-{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j},\text{(78)}& {\displaystyle}& {\displaystyle \frac{\partial {\stackrel{\mathrm{\u0303}}{B}}_{i,j}}{\partial t}}={N}^{\mathrm{2}}{D}_{i,j},\text{(79)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){P}_{i,j}=f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}+{\stackrel{\mathrm{\u0303}}{B}}_{i,j},\end{array}$$

where

$$\begin{array}{ll}\text{(80)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{1},j}+{P}_{i-\mathrm{1},j}+{P}_{i,j+\mathrm{1}}+{P}_{i,j-\mathrm{1}}-\mathrm{4}{P}_{i,j}\right),\\ {\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u0303}}{B}}_{i,j}\equiv \left({\displaystyle \frac{\partial}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){B}_{i,j},\\ {\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{D}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i,j}+{D}_{i+\mathrm{1},j}+{D}_{i,j+\mathrm{1}}+{D}_{i+\mathrm{1},j+\mathrm{1}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(81)}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}\equiv & {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.{\displaystyle}& {\displaystyle}\left.+{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right].\end{array}$$

$$\begin{array}{ll}\text{(82)}& {\displaystyle}& {\displaystyle \frac{\partial {\left({\mathit{\omega}}_{z}\right)}_{i,j}}{\partial t}}=-f{\stackrel{\mathrm{\u203e}}{D}}_{i,j},\text{(83)}& {\displaystyle}& {\displaystyle \frac{\partial {D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}}{\partial t}}={\displaystyle}& {\displaystyle}f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}},\\ \text{(84)}& {\displaystyle}& {\displaystyle \frac{\partial {\stackrel{\mathrm{\u0303}}{B}}_{i,j}}{\partial t}}={N}^{\mathrm{2}}{D}_{i,j},\text{(85)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}={\displaystyle}& {\displaystyle}f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\stackrel{\mathrm{\u0303}}{B}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}},\end{array}$$

where

$$\begin{array}{ll}\text{(86)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{{d}^{\mathrm{2}}}}\left({P}_{i+\mathrm{3}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{P}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{3}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}+{P}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\left.-\mathrm{4}{P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right),\end{array}$$

$${\stackrel{\mathrm{\u0303}}{B}}_{i,j}\equiv \left({\displaystyle \frac{\partial}{\partial z}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}H}}\right){B}_{i,j},$$

$$\begin{array}{ll}{\displaystyle}{\stackrel{\mathrm{\u203e}}{D}}_{i,j}\equiv & {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{D}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}\left.+{D}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{D}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right),\end{array}$$

and

$$\begin{array}{ll}\text{(87)}& {\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\equiv {\displaystyle}& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j+\mathrm{1}}+{\left({\mathit{\omega}}_{z}\right)}_{i,j+\mathrm{1}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j}+{\left({\mathit{\omega}}_{z}\right)}_{i,j}\right].\end{array}$$

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 four-point 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 right-hand 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

$$\begin{array}{ll}\text{(88)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{(*)}={\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n\right)}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}f{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n\right)},\end{array}$$

$$\begin{array}{}\text{(89)}& {\displaystyle}{\displaystyle}{D}_{i,j}^{(*)}={\stackrel{\mathrm{\u203e}}{D}}_{i,j}^{\left(n\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}\left[\underset{\mathrm{\xaf}}{f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}^{(*)}}-{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}\right],\end{array}$$

$$\begin{array}{}\text{(90)}& {\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u0303}}{B}}_{i,j}^{(*)}={\stackrel{\mathrm{\u0303}}{B}}_{i,j}^{\left(n\right)}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\tau}{N}^{\mathrm{2}}{\stackrel{\mathrm{\u203e}}{D}}_{i,j}^{\left(n\right)},\text{(91)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){P}_{i,j}=\underset{\mathrm{\xaf}}{f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}^{(*)}}+{\stackrel{\mathrm{\u0303}}{B}}_{i,j}^{(*)},\end{array}$$

where ${\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i,j}$ is given by Eq. (80). In the above equations,

$$\begin{array}{ll}{\displaystyle}{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}^{(*)}\equiv & {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}{\left.{\left({\mathit{\omega}}_{z}\right)}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}\right]}^{(*)},\end{array}$$

$$\begin{array}{ll}{\displaystyle}{\stackrel{\mathrm{\u203e}}{D}}_{i,j}^{\left(n\right)}\equiv & {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({D}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}+{D}_{i+\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right.\\ {\displaystyle}& {\displaystyle}{\left.+{D}_{i-\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}+{D}_{i-\mathrm{1}/\mathrm{2},j-\mathrm{1}/\mathrm{2}}\right)}^{\left(n\right)}.\end{array}$$

$$\begin{array}{ll}\text{(92)}& {\displaystyle}& {\displaystyle}{\left({\mathit{\omega}}_{z}\right)}_{i,j}^{\left(n+\mathrm{1}\right)}={\left({\mathit{\omega}}_{z}\right)}_{i,j}^{\left(n\right)}-\mathit{\tau}f{D}_{i,j}^{(*)},\text{(93)}& {\displaystyle}& {\displaystyle}{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n+\mathrm{1}\right)}={\displaystyle}& {\displaystyle}{D}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n\right)}+\mathit{\tau}{\left[\underset{\mathrm{\xaf}}{\underset{\mathrm{\xaf}}{f{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}^{\left(n+\mathrm{1}\right)}}}-\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)\right]}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}},\\ \text{(94)}& {\displaystyle}& {\displaystyle}{\stackrel{\mathrm{\u0303}}{B}}_{i,j}^{\left(n+\mathrm{1}\right)}={\stackrel{\mathrm{\u0303}}{B}}_{i,j}^{\left(n\right)}+\mathit{\tau}{N}^{\mathrm{2}}\underset{\mathrm{\xaf}}{\underset{\mathrm{\xaf}}{\underset{\mathrm{\xaf}}{{\stackrel{\mathrm{\u203e}}{D}}_{i,j}^{\left(n+\mathrm{1}\right)}}}},\text{(95)}& {\displaystyle}& {\displaystyle}{\left({\stackrel{\mathrm{\u0303}}{\mathrm{\nabla}}}_{H}^{\mathrm{2}}P\right)}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}-\left({m}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}^{\mathrm{2}}}}\right){P}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}={\displaystyle}& {\displaystyle}f\underset{\mathrm{\xaf}}{\underset{\mathrm{\xaf}}{{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n+\mathrm{1}\right)}}}+{\stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u0303}}{B}}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n+\mathrm{1}\right)},\end{array}$$

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,

$$\begin{array}{ll}{\displaystyle}{\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n+\mathrm{1}\right)}\equiv & {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\left({\mathit{\omega}}_{z}\right)}_{i,j}+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j}+{\left({\mathit{\omega}}_{z}\right)}_{i,j+\mathrm{1}}\right.\\ {\displaystyle}& {\displaystyle}{\left.+{\left({\mathit{\omega}}_{z}\right)}_{i+\mathrm{1},j+\mathrm{1}}\right]}^{\left(n+\mathrm{1}\right)},\end{array}$$

$$}{\displaystyle}{\stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u0303}}{B}}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{\left(n+\mathrm{1}\right)}\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\stackrel{\mathrm{\u0303}}{B}}_{i,j}+{\stackrel{\mathrm{\u0303}}{B}}_{i+\mathrm{1},j}+{\stackrel{\mathrm{\u0303}}{B}}_{i,j+\mathrm{1}}+{\stackrel{\mathrm{\u0303}}{B}}_{i+\mathrm{1},j+\mathrm{1}}\right)}^{\left(n+\mathrm{1}\right)}.$$

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 single-underlined terms are replaced by ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i,j}^{(*)}$, the double-underlined terms are replaced by ${\stackrel{\mathrm{\u203e}}{\left({\mathit{\omega}}_{z}\right)}}_{i+\mathrm{1}/\mathrm{2},j+\mathrm{1}/\mathrm{2}}^{(*)}$, and the triple-underlined 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.

We have integrated the models described above to check the normal-mode 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 CD-grid model, we use a simple forward time-integration scheme that produces virtually neutral (very weakly unstable) solutions with short time steps in all models. We also tested other time-differencing 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 second-order Adams–Bashforth scheme was slightly dissipative.

We first present results from the standing oscillation simulations, to look
for the frequencies obtained through the normal-mode 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 normal-mode 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 upside-down 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 (upside-down and inverted) pyramids. For $d=L/\mathrm{2}$, which corresponds to
the shortest horizontal grid spacing to resolve the initial perturbation, the
pyramid-like structure is again visible with the grid points at the extremes
of the (upside-down 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 D-grid 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 D-grid 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 D-grid 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
Z-grid simulation. In the C-, CD- and D-grid 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 D-grid 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 CD-grid 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. Non-oscillating solutions are also obtained by starting from the vorticity perturbations (as indicated by zeros within parentheses in Table 5). In the C-grid 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 C-grid solutions are very close to those of the Z-grid solutions. For the long horizontal and short vertical scales, the frequency with which the buoyancy and divergence oscillate in the C-grid solutions is considerably smaller than that of the Z-grid 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 D-grid simulations, the solutions are non-oscillatory for all variables (indicated by zeros within parentheses in Table 5). The Z-grid solution produces frequencies that are close to the true frequencies even though the initial perturbation is poorly resolved.

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
bell-shaped buoyancy perturbation placed in the middle of the horizontal
domain with rapidly decaying amplitude away from the center. We modified the
bell-shaped perturbation by setting it to zero at every other grid point. We
refer to this modification as “initial grid-scale 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 124-by-124 wide corner-end portion of the
horizontal domain, after 100 simulated minutes of integration. For reference,
we show the Z-grid result obtained without the superimposed initial
grid-scale noise. No major difference can be seen between the Z-grid
solutions started with the computational mode and without it in the portion
of the domain shown in Fig. 13. In the solution with the grid-scale noise,
there is, however, a remnant of the initial noise, which takes the form of
grid-scale 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 C-grid simulation result is very close to
the Z-grid 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 D-grid results resemble the C-grid 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 D-grid 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.

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 CD-grid 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 CD-grid simulations. In the Z- and C-grid 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 CD-grid simulations, the divergence response is not strong enough to prevent an increase of the amplitudes of the buoyancy perturbations. The D- and CD-grid 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 CD-grid 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).

7 The impact of the nonhydrostatic effects on the horizontal grid selection

Back to toptop
Discretization errors have most often been studied using the shallow-water
equations. This is justifiable for an assessment of the errors with the Lamb
wave and quasi-hydrostatic modes. The dispersion of the inertia–gravity modes
for the shallow-water 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 shallow-water system can also be
used to assess the discretization errors of the quasi-hydrostatic 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 shallow-water systems. A detailed discussion of the continuous
and discrete dispersion equations can be found in the Supplement.

In this respect, the analogy between the shallow-water (also
quasi-hydrostatic) and anelastic systems is weak. The nonhydrostatic effects
of the anelastic system, as with the fully compressible, unified and
semi-hydrostatic 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 quasi-hydrostatic 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 quasi-hydrostatic frequencies in Fig. 15.

The Z and C grids actually perform better with the nonhydrostatic systems
than they do with the quasi-hydrostatic or shallow-water systems. Consider
the discrete dispersion equations (Eqs. 18 and 24) for the Z and
C grids, respectively. The finite-difference 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 quasi-hydrostatic or shallow-water 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 quasi-hydrostatic discrete systems.

8 Summary and conclusions

Back to toptop
We have discussed the horizontal and vertical discretizations of the three-dimensional 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 middle-latitude inertia–gravity waves. The use of a three-dimensional nonhydrostatic system instead of the two-dimensional shallow-water 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 Z-grid 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
cloud-resolving 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 C-grid solution, but its impact may not be severe because it
disperses like a pure gravity mode. However, since the three-dimensional
enstrophy cascades to the SRHS in nonlinear cloud-resolving 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 C-grid solution behave like
those of the Z-grid 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 C-grid solutions
with a typical climate model horizontal grid spacing of 100 km are as
accurate as the Z-grid 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 non-unique) 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 normal-mode 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 shallow-water analyses by various authors. The Z and C grids produce overall smaller errors with the nonhydrostatic anelastic system than with the quasi-hydrostatic system, particularly for the deep modes with high horizontal wavenumbers.

Part 2 (Konor and Randall, 2018) discusses the dispersion of the middle-latitude Rossby waves for all the horizontal grids.

Code and data availability

Back to toptop
Code and data availability.

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.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-11-1753-2018-supplement.

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements

Back to toptop
Acknowledgements.

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 AGS-1500187, the US Department of Energy Office of
Science DE-SC07050 (SciDAC), DE-SC00016273 (ACME) and DE-SC00016305
(CMDV).

Edited by: Paul Ullrich

Reviewed by: Almut Gaßmann and one anonymous referee

References

Back to toptop
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.: Next-generation ocean and ice models at GFDL: MOM6, SI2, and icebergs, available at: http://cosima.org.au/wp-content/uploads/2016/06/MOM6-overview-May-27.pdf, 2016.

Arakawa, A.: Finite-difference 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 Charney-Phillips grid in hybrid *σ*-*p* vertical
coordinates, Mon. Weather Rev., 124, 511–528, 1996.

Arakawa, A. and Konor, C. S.: Unification of the anelastic and quasi-hydrostatic 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 quasi-geostrophic equations for barotropic and simple baroclinic flows, J. Meteorology, 10, 71–99, 1953.

Chen, Q.: On staggering techniques and the non-staggered 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 normal-mode analysis, Q. J. Roy. Meteor. Soc., 129, 2761–2775, 2003.

Dubos, T. and Voitus, F.: A semihydrostatic theory of gravity-dominated 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., McTaggard-Cowan, 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 log-hydrostatic-pressure 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 three-dimensional 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 “Charney-Phillips” arrangement, ECMWF
Tech. Memo. No. 211, 12 pp., 1995.

Janjic, Z. I.: Nonlinear advection schemes and energy cascade on semi-staggered 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: Quasi-geostrophic Rossby modes, Geosci. Model Dev., 11, 1785–1797, https://doi.org/10.5194/gmd-11-1785-2018, 2018.

Lee, J.-L. and MacDonald, A. E.: A finite-volume icosahedral shallow-water model on a local coordinate, Mon. Weather Rev., 128, 3901–3910, 2009.

Lin, S.-J. and Rood, R. B.: Multidimensional flux form semi-Lagrangian 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 finite-difference shallow-water 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 finite-volume 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 quasi-uniform 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 mass-conserving semi-implicit semi-Lagrangian discretization of the deep-atmosphere global non-hydrostatic equations, Q. J. Roy. Meteor. Soc., 140, 1505–1520, 2014.

Short summary

We have discussed the discretizations of the three-dimensional nonhydrostatic linearized anelastic equations on the A, B, C, CD, (DC), D, E and Z horizontal grids, and on the L and CP vertical grids, with an emphasis on midlatitude inertia–gravity waves. The Z and C grids show the most accurate dispersion among the seven horizontal grids. The inertia–gravity mode solutions with the D and CD grids are almost identical. The A, B and E grids suffer from the multiple (or non-unique) physical modes.

We have discussed the discretizations of the three-dimensional nonhydrostatic linearized...

Geoscientific Model Development

An interactive open-access journal of the European Geosciences Union