## Introduction

A two-level pressure treatment is described that ensures that the anelastic mass-continuity equation is satisfied to round-off level. Satisfying mass continuity significantly increases model stability and accuracy. In the presence of modestly steep topography, the computational burden of the diagnostic elliptic pressure solver is equivalent to that of a direct solver.

The procedure described represents a recent development to the nonhydrostatic anelastic terrain-following coordinate model initially developed by Clark (1977). In that model, the diagnostic pressure was solved using dimension reduction, resulting in decoupled two-dimensional Helmholtz equations for each vertical grid level. A block-iteration algorithm was used in which the pressure terms associated with topographic gradients were taken to the right-hand side and treated explicitly; that is, these terms were not part of the direct solver. One of the drawbacks of this approach is the slow rate of convergence of the pressure solution for relatively steep topography (e.g., Schumann and Volkert 1984). One measure of convergence used by Clark (1977) is the mean standard deviation of the mass divergence. A potential advantage of this type of algorithm is that the major computational effort is spent on the vertical dimension reduction that occurs in the hyperbolic or forward-marching portion of the code. In nonvertically stretched coordinates, fast Fourier transforms can be used (Ogura 1969). Dimension reduction is an advantage in a distributed memory computational environment because no data exchange between nodes is required for the vertical dimension reduction. Internodal data transfer first occurs during the establishment of the two-dimensional Helmholtz equations. Before describing the new techniques, some of the model's developments will be briefly described.

Multiple levels of vertical and horizontal grid refinement were added by Clark and Farley (1984) and, more recently, vertical grid stretching with vertical refinement was added by Clark and Hall (1996). Clark et al. (1997) converted the code to shared-memory open multiprocessing parallelization using domain decomposition in the two horizontal dimensions. This domain decomposition architecture was recently converted from open multiprocessing software, designed for shared-memory systems (Open MP) to distributed-memory message passing interface (MPI). It is in part a result of the MPI conversion that led to the testing of some of the ideas presented in this paper. Another incentive was to find an efficient pressure solver that could be used for relatively steep topography.

The two-level pressure technique is a simple procedure and is best suited to block iterative methods. Instead of using only the latest pressure *p*^{υ} (where *υ* denotes the number of iterations) to update the velocities, both *p*^{υ} and *p*^{υ−1} are used as field variables in a manner ensuring that mass continuity is satisfied to round-off error. This adds a significant degree of stability to the model in the presence of relatively steep terrain, disallowing the mass divergence to contribute to kinetic energy production through the nonlinear advection terms. This procedure is similar in some respects to the Piacsek and Williams (1970) technique in which they added corrective terms to the nonlinear advection terms to reduce the effect of mass divergence. Using this two-level technique, the mean standard deviation of mass divergence no longer serves as a measure for the convergence of the pressure.

A procedure is described in which metric filtering can be used indirectly as part of the block iteration in cases for which the topography exceeds about 70% (35°), requiring multiple block iterations. In such cases, the filtering is progressively decreased to zero with each block iteration, resulting in a solution to the full or unfiltered elliptic equation. The majority of model applications require only a single block iteration to solve the diagnostic pressure equation, and so this technique is typically not used.

Simple experiments of flow over a two-dimensional ridge with both constant and vertically varying air densities are used to assess the two-level pressure technique. The flow involves breaking waves giving a strong dynamic response that should provide an evaluation of these techniques under severe conditions. The errors are also compared with those that would be expected using a numerical quadrature from Bernardet (1995). Comparisons are performed at both early and later times in the flow evolution. This approach allows evaluation during both the transient and the near-steady-state periods of the flow. It will be shown that the errors are modest for all of the cases evaluated.

## Dynamic model equations

*ϕ*(

*x*+ Δ

*x*/2) +

*ϕ*(

*x*− Δ

*x*/2)]/2 and

*δ*

_{x}

*ϕ*= [

*ϕ*(

*x*+ Δ

*x*/2) −

*ϕ*(

*x*− Δ

*x*/2)]/Δ

*x*}, the three-dimensional momentum equations are approximated aswhere

*λ*is the latitude (taken as zero in this paper);

*p*and

*ϕ*are two different iteration levels of the pressure. Here

*ρ*is the Jacobian-weighted air density, where

*ρ*=

*ρ*

_{0}(

*z*)(1 −

*h*/

*H*) →

*ρ*

_{00}(

*ξ*) when

*h*→ 0, which means that for zero topography the right-hand side of (3) reduces to pf

*z*= −cos

*λρ*

_{0}

*δ*

_{ξ}(

*p*/

*ρ*

_{0}), which is a form in which the perturbation pressure does not contribute to vorticity production (Lipps and Hemler 1982). It is also a form that is, in general, energy conservative and reduces to the deep equations of Ogura and Phillips (1962) when

*θ*is expanded about a constant. In (1)–(3), FX, FY, and FZ represent the advection, Coriolis, buoyancy, and mixing terms. The pressure terms (pf) remain on the right-hand side for the purpose of illustrating the two-level pressure formulation.

*δ*

_{x}

*ρ*

^{x}

*u*

*δ*

_{y}

*ρ*

^{y}

*υ*

*δ*

_{ξ}

*ρ*

^{ξ}

*ω*

*ω*(=

*dξ*/

*dt*) and the Cartesian components of the air velocity is taken as

*z*

*F*

*ξ*

*h*

*H*

*h,*

*h*=

*h*(

*x,*

*y*) is the height of the topography,

*z*is the Cartesian height,

*H*is the fixed height of the model domain, and

*ξ*is the transformed vertical coordinate in the nonorthogonal system. In the finite-difference formulation,

*ξ*is taken as

*ξ*

_{k}

*k*

*ξ*

_{k}

*ξ*is a constant and

*k*is the vertical grid index. The boundary conditions on

*F*(

*ξ*) are taken as

*F*

*F*

*H*

*H.*

*ξ*) =

*ξ,*the above system reduces to the shearing coordinates as used in Gal-Chen and Somerville (1975), Clark (1977), and Schumann et al. (1987), among others. In the general systemBefore deriving the diagnostic pressure equation, (1)–(3) and (5) must be combined to form the

*ω*tendency equationwhereHere, FΩ includes the well-known Christoffel terms. However, because they are derived using the discretized equations, using the kinetic-energy-conservative form (Arakawa 1966), the nonlinear advection terms conserve both momentum and kinetic energy to

*O*(Δ

*t*

^{2}) (Vinokur 1974; Clark 1977). This is in contrast to the pioneering work of Gal-Chen and Somerville (1975) in which they discretized the

*u,*

*υ,*and

*ω*equations in a nonconservative manner in which neither momentum nor kinetic energy was conserved. A conservative form of the

*ω*tendency equation can be obtained by expanding the left-hand side of (10).

We combine (1), (2), (10), and (4) to obtain a diagnostic elliptic pressure equation. A tendency term for mass continuity appears in this equation, and, as is common procedure, the future value is put to zero. An accurate solution of this equation results in mass continuity [(4)] being satisfied at each grid point of the model. With much less rigor we can also ensure that mass continuity is satisfied if we treat all terms in (1), (2), and (10) that contain *ϕ* as being explicit and take them to the left-hand side. Now, when we form the diagnostic pressure equation, the explicit terms contain expressions using past iteration values of the pressure. If we solve this elliptic equation (which turns out to be the matrix block) accurately and continue to treat the pressure gradient terms containing *ϕ* explicitly when we update the velocity components, then again mass continuity will be satisfied. Separating the pressure into the two components *p* and *ϕ* and then applying a direct solver to invert the block for *p* represents one block iteration.

*Q*represents the divergence of the forcing terms FX, FY, FΩ, and the mass divergence tendency.

*p*and

*ϕ*in (1), (2), (3), and (10) it is easily shown that

*L*

*p,*

*ϕ*

*L*

_{0}

*p*

*L*

_{1}

*ϕ*

*Q,*

*p*is derived aswhere from (9)

*G*

*h*(

*x,*

*y*)/

*H.*

### Block-iteration solution

*G*(

*x,*

*y*), the left-hand side of (14) [or

*L*

_{0}(

*p*)] is separable and easily solved using direct methods. By putting

*ϕ*=

*p*

^{υ−1}, (14) is solved using dimension reduction (e.g., Ogura 1969; Haltiner and Williams 1980, p. 168) in the vertical or

*ξ*coordinate. This is achieved by first letting [

*G*

*p*] = 𝗧

^{−1}[

*r*] where the square brackets denote a column vector and 𝗧

^{−1}is the inverse eigenvector matrix of the vertical operator in (14). Multiplying (14) by 𝗧 results inwhere

*A*represents the vertical operator in (14). The formation of

*A*requires the application of Neumann boundary conditions on

*p*at the model surface and top using (10) consistent with

*ω*= 0. The order of

*A*is equal to the number of vertical levels, meaning that the boundary conditions on

*p*are treated implicitly. Because

**Λ**is the diagonal eigenvalue matrix, (15) represents a series of decoupled two-dimensional Helmholtz equations. They are solved using a combination of iterative and direct methods. For example, the larger eigenvalue equations require as few as three red–black Gauss–Seidel iterations (i.e., in which iterations progress through the grid in a checkerboard fashion, which is easily vectorized). Most of the computation cost of the pressure solver is in forming (15) and inverting [

*r*] back to

*p.*Both of these computations are performed on local nodes, and as such the method is well suited for the distributed memory architecture.

In the treatment of Clark (1977) and Schumann et al. (1987), block iterations were performed until |*p*^{υ} − *p*^{υ−1}| had acceptably converged. After the termination of the series of block iterations, *ϕ* was set equal to the current values of *p,* and the pressure-forcing terms, pf*x*(*p,* *p*), pf*y*(*p,* *p*), and pf*ω*(*p,* *p*) in (1)–(3) would be calculated. These terms would then be used when calculating the advanced-time-level velocity components. A simple measure of convergence was to monitor the root-mean-square mass divergence error. This measure has been the main test procedure used for this particular model since inception.

### Two-level pressure approach

A very simple modification to the procedures of Clark (1977) and Schumann et al. (1987) is to retain *ϕ* as a field variable and to use it to calculate pf*x*(*p,* *ϕ*), pf*y*(*p,* *ϕ*), and pf*ω*(*p,* *ϕ*). The reason for keeping *ϕ* as a field variable is that it can represent the previous time level of *p* instead of simply representing a last iteration during the current time step. This seemingly trivial modification to the block-iteration procedure assures that mass continuity is satisfied to round-off error, provided that one solves the two-dimensional Helmholtz equations to round off. The root-mean-square mass divergence diagnostic is now used to monitor the accuracy of the Helmholtz equation solutions. An immediate benefit of this approach is a notable increase in model stability. Using this two-level pressure technique, we have found that for most experiments we can use a single block iteration when solving the pressure (e.g., Clark et al. 1997, 2000).

### Block iteration with metric filtering

A practical technique used with terrain-following-coordinate models is to filter the coordinate transformation metrics (R. Meyer 1996, personal communication). This means that only the metric terms associated with the horizontal gradients of topography are filtered, that is, *J**J*^{13} and *J**J*^{23} but not *J**x*(*p,* *ϕ*), pf*y*(*p,* *ϕ*), and pf*ω*(*p,* *ϕ*) pressure gradient terms so that neither (4) nor (5) is directly affected by this approximation. Also, the treatment of diffusion terms is not affected, and so variables continue to mix along horizontal surfaces when appropriate. Although the model can be run using metric filtering, its discussion in this paper is confined to how it can be used as part of the pressure solver when multiple block iterations are required.

Metric filtering is used as a method to obtain a solution to (14) using progressively decreasing metric filtering with each block iteration. All terms in the prognostic equations, including the pressure gradient terms, use unfiltered values of *J**J*^{13} and *J**J*^{23}. Again, the procedure is seemingly trivial: one may start with metrics filtered using, say, five applications of a filter. After the completion of the first block iteration, these metrics are replaced with four applications of the filter and so on until the last block iteration uses the unmodified metrics. No attempt has been made to optimize this iteration procedure, but there certainly seems to be room for future exploration. The full block-iterative method with metric filtering provides us with a robust enough solver to treat grades in excess of 100% or 45° slopes.

Numerical derivations of the kinetic energy and associated pressure truncation errors are presented in appendix A. Appendix B presents the numerical form of the pressure truncation errors for both the Clark (1977) and a particular Bernardet (1995) formulation under the simplified conditions when *f* = 0 and *ρ*_{0}(*z*) = constant.

## Experimental design

A simple three-dimensional experiment of stable air flowing over an east–west-oriented two-dimensional ridge was chosen to test the performance of the two-level pressure approach. Two cases were chosen. The first used a constant-with-height wind of *V*_{0} = 30 m s^{−1} directed normal to the ridge and with *ρ*(*z*) = 1. The second case allowed *ρ* to vary with height and used a constant-with-height wind of 15 m s^{−1}. The potential temperature *θ* was taken as a constant between *z* = 0 and 1 km. Above 1 km, a constant Brunt–Väisälä frequency of *N* = 0.01 s^{−1} was applied. Above *z* = 14 km, *N* = 0.02 was used. The topography was taken as a Gaussian-shaped ridge with a height of 6 and 3 km for the first and second cases, respectively. A half-width of 30 km was used. There is nothing special about these values other than that both experimental setups provided solutions with intense lee-slope winds resulting from gravity wave breaking. The horizontal spacing of the model was 3.7 km, and so the ridge is well resolved. The domain was taken as 444 km long in both the *x* and *y* directions.

This experimental framework results in breaking gravity waves causing the mean flow to amplify through reflective resonance (Peltier and Clark 1979). The mean flow speeds are amplified to over 100 and 45 m s^{−1} for the first and second cases, respectively. These flows result in strong local pressure gradients that should provide a good test of the pressure solution techniques being assessed. Table 1 describes the basic experiments. Experiments X1 and X5 represent the basic two-level pressure case (using a single block iteration) for constant and variable density, respectively. One can compare X1 (X5) with X2 (X4) to determine differences resulting from using the two-level approach with a single block iteration versus solutions with a fully converged pressure. The purpose of X3 is to assess the effect of using metric filtering as part of the block-iteration process. Experiments X2 and X4 are essentially the one-level approach of Clark (1977) because the two pressures have converged, which is why such a large value was chosen for *N*_{BLOCK}. To validate that X4 had converged, experiment X6 was also performed in which the two-level approach was replaced with the usual one-level approach; that is, only the most recent calculated pressure was used in (1)–(3).

## Results

### Basic dynamics

*θ*for experiments X1 and X2. Even with

*ρ*(

*z*) = 1, the flow results in breaking hydrostatic gravity waves that amplify through reflection, leading to intense winds of over 100 m s

^{−1}(e.g., Peltier and Clark 1979). There is also evidence of downstream nonhydrostatic resonant lee waves between

*z*= 0 and 6 km. The difference between experiments X1 and X2 at either the early (

*t*= 10 min) or late (

*t*= 60 min) time is negligible. Figure 2 shows the horizontal velocity for X1 at

*t*= 10 and 60 min and the differences between X2 and X1 at the same times. This figure shows

*υ*attaining maximum of 95 and −49 m s

^{−1}at

*t*= 60 min, whereas the maximum differences are ±1.6 m s

^{−1}. The maximum differences at

*t*= 60 min occur primarily in the wave-breaking region, where nonlinear effects are largest. It is well known that nonlinear solutions with slightly different initial conditions diverge with time. The differences between X2 and X1 shown in Fig. 2 are not only small but are within the range one might expect to result from predictability limits caused by nonlinear effects. For this reason, errors associated with the kinetic energy budget are used to assess the relative magnitude of the two-level pressure technique. The pressure truncation errors associated with the numerical formulation of the equations can be used as a baseline measure because these (as described in appendixes A and B) are an inherent part of the terrain-following system. The erroneous rate of kinetic energy production caused by the contribution from the two-level pressure technique can be estimated using the converged pressure from experiment X2 (

*p**) and the two pressures from X1 (

*p,*

*ϕ*) aswhere

*ũ,*

*υ̃*

*w̃*are defined in (A1). The results from (16) will be shown in the next section.

A variant of experiment X1 was also attempted with the two-level pressure technique turned off (i.e., using the approach of Clark 1977) with a single block iteration and was found to be numerically unstable. That particular experiment stopped after 57 time steps. These results indicate that the two-level pressure technique significantly improves numerical stability by ensuring that mass continuity is satisfied.

### Pressure truncation analysis

The kinetic energy source/sink terms described in appendix B are associated with truncation errors resulting from the numerical approximation of the pressure gradient terms. As discussed earlier, these terms are always zero for the analytical equations but are nonzero for the numerically approximated equations where topography is present. In the case of nonzero topography, Γ in (A3) has finite values, and it is of interest, then, to assess its magnitude to see whether it might represent a significant source of energy and to see whether it differs substantially among different numerical formulations of the equations. Because numerical instability is typically dominated by the nonlinear terms, it seems unlikely that these linear pressure truncation error terms would do much more than lead to systematic errors.

Figure 3 shows the second-order pressure truncation errors Γ for experiments X1 and X2 calculated using (A3) and the two-level errors Δke using (16). We see some differences between the Γ errors for the two-level pressure case using a single block iteration (X1) and the two-level case using 10 block iterations (X2) but only at the early time of *t* = 10 min when the flow is still relatively transient. As the flow approaches steady state at *t* = 60 min these errors are nearly identical. The Δke errors are between 5 and 10 times as small as the Γ errors at both times. The spatial distribution of the Δke differs considerably from that of the Γ errors because they result from both the transient behavior of the flow as well as from grid deformations resulting from the coordinate transformation, whereas the Γ errors mostly reflect effects from grid deformations.

A comparison between X3 and X2 (not shown) showed negligible differences, indicating that block iterations with metric filtering gave effectively the same results as without for this case. Other experiments used block iterations with metric filtering successfully for much steeper topography (about 50°), obtaining effectively converged solutions.

Figure 4 shows the pressure truncation errors for X1B calculated a posteriori using (A3) and *p* from X1. At *t* = 10 and 60 min, the errors resulting from the Bernardet (1995) formulation (described in appendix A) are very similar to those for X1 and X2. This particular experiment suggests there would be little advantage in converting the horizontal pressure gradient forcing to (B1), which is cast in a nonconservative form with respect to the conservation of momentum.

At *t* = 60 min, the errors shown in Fig. 3 are near their maximum amplitudes for the experiment. A rough estimate of their impact can be made by considering a trajectory passing through the positive region near the surface for X1 at *t* = 60 min. Assuming a width of 10 km, a mean speed of 100 m s^{−1}, and an average error of 0.06 W m^{−3}, the error in the wind speed would be about 3.5 m s^{−1}. The source/sink region is narrow and we would expect a similar deceleration to follow as a result of the negative (nearly equal) contribution. This error is still relatively small because of resolving the mountain half-width with eight grid points. At lower resolutions, modelers might consider evaluating these errors because they could seriously affect mountain–valley and/or flow separation calculations.

### Variable-density case

Experiments X4 and X5 were performed to demonstrate the nature of the pressure truncation error terms for vertically varying air density. The experiments were modified by reducing the mountain height to 3 km and the ambient wind speed to 15 m s^{−1}. The flow resulting from experiments X1–X3 with constant density was unrealistically intense. The current experiment is much closer to what one might expect in typical windstorm simulations.

Figure 5 shows *θ* for X4 at *t* = 10 and 60 min, showing the wave-breaking nature of the flow. Differences in *θ* between X5 and X4 are negligible and are not shown. A comparison between X6 and X4 showed negligible differences, indicating that X4 represents a converged solution. Figure 6 shows Γ for X4 and X5 and Δke using the converged *p** from X4 and (*p,* *ϕ*) from X5. We see that the amplitudes of the Γ errors are now much smaller when compared with those obtained in experiments X1 and X2. In this case, the regions Γ form narrow patches near the surface and are only one or two grid points wide in the direction of the flow. An estimated speed error resulting from these pressure truncation errors is about 1.4 m s^{−1}. This estimate assumed an average error of 0.01 W m^{−3} over a distance of about 4 km where the wind speed is about 45 m s^{−1}, which corresponds to the strong positive region of error at *t* = 60 min. This represents an error of about 3.1% that is similar to the earlier estimates of 3.5%. As in the previous section, the Δke errors are again between 5 and 10 times as small as the Γ errors, which are an inherent part of the terrain-following coordinates.

Figure 7 shows the extreme values of *u,* *υ,* *w,* and form drag versus time for experiments X4 and X5. The extremes are shown for X4 (heavy lines), whereas differences between X4 minus X5 multiplied by a factor of 25 are shown as thin lines. The differences for *υ,* *w,* and form drag reach their maximum amplitudes during the early transient stage and asymptotically approach zero as the solutions approach steady state. Because the experimental design is a 2D–3D experiment with mean flow in the *y* direction, the values of *u* become significant only after round-off error excitation and nonlinear interactions lead to the development of the *u* component of flow. At about *t* = 50 min we see a significant increase in the magnitudes of *u* as a result of gravity wave breaking and developing nonlinearities.

## Conclusions

One case chosen to evaluate the two-level pressure technique was that of constant-density stable air flowing over a 6-km-high Gaussian-shaped ridge. This flow evolved into extremely severe winds, probably well beyond anything that might occur naturally. The ambient winds of 30 m s^{−1} were amplified to over 110 m s^{−1} through the effects of wave breaking. A second case was chosen in which air density varies in the vertical. The mountain height was reduced to 3 km and the ambient flow to 15 m s^{−1}. This case produced strong winds of about 45 m s^{−1}, more in line with naturally observed cases.

All comparisons indicate that the two-level pressure technique added considerable numerical stability and accuracy to the model. The increased numerical stability was argued to be a direct result of the two-level pressure technique, guaranteeing that the mass-continuity constraint is satisfied to effectively minimize round-off error irrespective of how many block iterations are performed. This has the effect of minimizing some of the important nonlinear energy sources that can lead to numerical instability in the presence of finite-amplitude mass divergence in the anelastic system. The typical application of this kind of anelastic model requires only a single block iteration. This means that the cost of solving the diagnostic pressure equation is effectively that of using a direct solver. This type of procedure is particularly well suited to the MPI environment because the majority of the computational effort is performed on the vertical dimension reduction that is parallelized in a manner equivalent to the hyperbolic (or forward marching) portion of the code. The only special data transfer required for the pressure solution is confined to the treatment of the horizontal Helmholtz equations. The two-level pressure approach would seem counterproductive for iterative schemes because its application would require iteration to an acceptable level of convergence on the matrix block rather than on the full operator.

The magnitude and distribution of the errors were evaluated and shown, for the constant-density case, to be relatively small because the mountain was well resolved. In this calculation a local error estimate of about 3.5 m s^{−1} for a wind speed of 100 m s^{−1} was estimated. The error estimate for the variable-density case was 1.4 m s^{−1} for a maximum wind of about 45 m s^{−1}, representing a 3.1% error. The errors resulting from the two-level pressure technique were found to be between 5 and 10 times as small as the second-order numerical errors Γ in the terrain-following system. Although the Γ errors were relatively small for the well-resolved topographic cases presented, these errors are large enough that numerical modelers studying mountain–valley circulations and/or flow separation in complex topography might consider evaluating their impact, particularly when the mountains are poorly resolved. The pressure truncation error analysis presented in appendix A might be particularly appropriate when mountains are approximated using a minimal number of grid cells. Also, the differences in Γ errors between model formulations could be sufficient to account for the level of differences found in the model comparison study of Doyle et al. (2000). A simple method of evaluating the Γ errors was presented in appendix A.

Numerical models attempt to conserve quantities, such as momentum and kinetic energy, for purposes of accuracy of the model budgets. However, stability is no guarantee of accuracy. Numerous derivations of the pressure gradient terms have been proposed using Taylor series analysis in Cartesian coordinates to determine the most accurate formulation of the gradient term to apply in the transformed system without evaluating the effect on the kinetic energy budget (or other moments of the system). The final error for the kinetic energy budget is a result of all the coupled terms, which include the horizontal and vertical pressure gradients, the numerical formulation of the relation between the Cartesian and transformed velocity components, mass continuity, and the time-stepping procedure. A study that might prove worthwhile would be to use the method described in appendix A to evaluate the conservative properties of schemes such as those proposed by Mahrer (1984) and others. Another study that might be worthwhile is to explore improved formulations of the horizontal pressure gradients that might lead to reduced Γ errors. Such an analysis must, of course, be performed in the transformed system of coordinates.

## Acknowledgments

I am particularly grateful to Henryk Modzelewski at the University of British Columbia and Richard Farley at the South Dakota School of Mines and Technology for their assistance with the MPI conversion. I also would like to acknowledge the generous assistance of Morwenna Griffiths and Michael Reeder at Monash. This work was supported in part by a Victorian Partnership for Advanced Computing expertise grant. This work was in part also supported by the Forest Renewal British Columbia Contract PAR 02002-16 as part of the Firestorm project at UBC. I thank Bob Sharman, Duane Stevens, and Todd Lane for useful discussions that helped to improve the paper.

## REFERENCES

Arakawa, A. 1966. Computational design for long term integration of the equations of fluid motion: Two-dimensional incompressible flow. Part I.

*J. Comput. Phys.*1:119–143.Bernardet, P. 1995. The pressure term in the anelastic model: A symmetric elliptic solver for an Arakawa C grid in generalized coordinates.

*Mon. Wea. Rev.*123:2474–2490.Clark, T. L. 1977. A small-scale dynamic model using a terrain-following coordinate transformation.

*J. Comput. Phys.*24:186–214.Clark, T. L., and W. R. Farley. 1984. Severe downslope windstorm calculations in two and three spatial dimensions using anelastic interactive grid nesting: A possible mechanism for gustiness.

*J. Atmos. Sci.*41:329–350.Clark, T. L., and W. D. Hall. 1991. Multi-domain simulations of the time dependent Navier–Stokes equations: Benchmark error analyses of some nesting procedures.

*J. Comput. Phys.*92:456–481.Clark, T. L., and W. D. Hall. 1996. The design of smooth, conservative vertical grids for interactive grid nesting with stretching.

*J. Appl. Meteor.*35:1040–1046.Clark, T. L., , T. Keller, , J. Coen, , P. Neilley, , H-M. Hsu, , and W. D. Hall. 1997. Terrain-induced turbulence over Lantau Island: 7 June 1994 Tropical Storm Russ case study.

*J. Atmos. Sci.*54:1795–1814.Clark, T. L., , W. D. Hall, , R. M. Kerr, , D. Middleton, , L. Radke, , F. M. Ralph, , P. J. Neiman, , and D. Levinson. 2000. Origins of aircraft-damaging clear-air turbulence during the 9 December 1992 Colorado downslope windstorm: Numerical simulations and comparison with observations.

*J. Atmos. Sci.*57:1105–1131.Doyle, J. D. Coauthors,. 2000. An intercomparison of model predicted wave breaking for the 11 January 1972 Boulder windstorm.

*Mon. Wea. Rev.*128:901–914.Gal-Chen, T., and R. C. J. Somerville. 1975. On the use of a coordinate transformation for the solution of the Navier–Stokes equations.

*J. Comput. Phys.*17:209–228.Haltiner, G. J., and R. T. Williams. 1980.

*Numerical Prediction and Dynamic Meteorology*. John Wiley and Sons, 477 pp.Lipps, F. B., and R. S. Hemler. 1982. A scale analysis of deep moist convection and some related numerical problems.

*J. Atmos. Sci.*39:2192–2210.Mahrer, Y. 1984. An improved numerical approximation of the horizontal gradients in a terrain-following coordinate system.

*Mon. Wea. Rev.*112:918–922.Ogura, M. 1969. A direct solution of Poisson's equation by dimension reduction method.

*J. Meteor. Soc. Japan*47:319–323.Ogura, M., and N. A. Phillips. 1962. Scale analysis of deep and shallow convection in the atmosphere.

*J. Atmos. Sci.*19:173–179.Peltier, W. R., and T. L. Clark. 1979. The evolution and stability of finite-amplitude mountain waves. Part II: Surface wave drag and severe downslope windstorms.

*J. Atmos. Sci.*36:1498–1529.Piacsek, S. A., and G. P. Williams. 1970. Conservative properties of convection difference schemes.

*J. Comput. Phys.*6:392–405.Schuman, F. G. 1962. Numerical experiments with the primitive equations.

*Proc. Int. Symp. on Numerical Weather Prediction,*Tokyo, Japan, Meteorological Society of Japan, 85–107.Schumann, U., and H. Volkert. 1984. Three-dimensional mass- and momentum-consistent Helmholtz-equation in terrain-following coordinates.

*Efficient Solutions of Elliptic Systems,*W. Hackbusch, Ed., Notes on Numerical Fluid Mechanics, Vol. 10, Vieweg, 109–131.Schumann, U., , T. Hauf, , H. Holler, , H. Schmidt, , and H. Volkert. 1987. A mesoscale model for the simulation of turbulence, clouds and flow over mountains: Formulation and validation examples.

*Beitr. Phys. Atmos.*60:413–446.Sharman, R. D., , T. L. Keller, , and M. G. Wurtele. 1988. Incompressible and anelastic flow simulations on numerically generated grids.

*Mon. Wea. Rev.*116:1124–1136.Vinokur, M. 1974. Conservation equations of gas dynamics in curvilinear coordinate systems.

*J. Comput. Phys.*14:105–125.

## APPENDIX A

### Pressure Sources and Sinks of Kinetic Energy

*ũ,*

*υ̃*

*w̃,*respectively, where

*ũ,*

*υ̃*

*w̃*are defined aswhere

*τ*is the time index (

*t*=

*τ*Δ

*t*). The resulting kinetic energy equation iswhere

*q*

^{2}=

*ρ*

^{x}

*u*

^{2}

^{x}+

*ρ*

^{y}

*υ*

^{2}

^{y}+

*ρ*

^{ξ}

*w*

^{2}

^{ξ}. The first bracketed term on the right-hand side of (A2) represents energy transport, the second term represents the kinetic-to-convective conversion term, where

*θ*∗ =

*θ*′/

*θ*

It is important to note that Γ → 0 as Δ*x,* Δ*y,* and Δ*ξ* → 0. ADV (ke) in (A2) represents the nonlinear advection and mixing terms. Using (A3), different numerical formulations of the equations can be assessed. Forms for both the Bernardet (1995) and the current formulation are derived for the simplified conditions of *f* = 0 and *ρ*_{0}(*z*) = constant in appendix B. This analysis is different from previous papers, such as Mahrer (1984), in which accuracy is assumed to result from a precise matching of the gradient terms to Cartesian space interpolations for which the effects of other contributing terms are not considered.

*D*has finite values. The form of these terms iswhere Δ

_{u}= Δ

*t*

^{2}/2

*δ*

_{tt}

*u.*

Probably the main reason for the increased stability of the model when we use the two-level pressure scheme is that *D* approaches zero, to round-off error limits, and the source terms in (A4) effectively vanish.

## APPENDIX B

### Numerical Forms of the Pressure Truncation Errors

*ω*and Cartesian

*u*and

*w*[which is identical to (5) but for constant density], andfor the anelastic mass continuity, where constant density and zero Coriolis force were assumed. The corresponding kinetic energy equation is found asThe Bernardet error terms were also calculated using (A4) but using the velocity pressure quadrature leading to (B5). To round-off level, the errors were found to be identical to those presented in (B5). Formulations with significantly reduced errors can be obtained using a fully generalized transformation in which the grids are numerically generated (e.g., Sharman et al. 1988).

*p*=

*p*−

*ϕ.*

Experimental setup: *N*_{BLOCK} represents the number of block iterations; α is the metric-filtering reduction factor [e.g., *h _{x}* → (α

*h*)]. Experiments X2 and X4 did not use any metric filtering during the block iterations. X1B is identical to X1 except the pressure truncation error analysis uses the Bernardet numerical forms. Each experiment was run out to 60 min using a time step of 10 s

_{x}^{*}

The National Center for Atmospheric Research is sponsored by the National Science Foundation.