## Abstract

The dynamic core of the Ocean–Land–Atmosphere Model (OLAM), which is a new global model that is partly based on the Regional Atmospheric Modeling System (RAMS), is described and tested. OLAM adopts many features of its predecessor, but its dynamic core is new and incorporates a global geodesic grid with triangular mesh cells and a finite-volume discretization of the nonhydrostatic compressible Navier–Stokes equations. The spatial discretization of horizontal momentum is based on a C-staggered grid and uses a method that has not been previously applied in atmospheric modeling. The temporal discretization uses a unique form of time splitting that enforces consistency of advecting mass flux among all conservation equations. OLAM grid levels are horizontal, and topography is represented by the shaved-cell method. Aspects of the shaved-cell method that pertain to the OLAM discretization on the triangular mesh are described, and a method of conserving momentum in shaved cells on a C-staggered grid is presented.

The dynamic core was tested in simulations with multiple vertical model levels and significant vertical motion. The tests include an idealized global circulation simulation, a cold density current, and mountain-wave flow over an orographic barrier, all of which are well-known standard benchmark experiments. OLAM gave acceptable results in all tests, demonstrating that its dynamic core produces accurate and robust solutions.

## 1. Introduction

The Ocean–Land–Atmosphere Model (OLAM) is a new numerical simulation model that is partly based on the Regional Atmospheric Modeling System (RAMS). The principal motivation for developing OLAM was to replace the limited-area configuration and lateral boundaries of RAMS with a global domain, thereby greatly expanding the range of systems that can be simulated and studied. Many components of RAMS were incorporated into OLAM with little modification, which saved considerable development effort. These include initialization methods, I/O data formats, model configuration options, and parameterizations of various physical processes. The OLAM dynamic core, however, is new because the RAMS grid structure and quasi-Boussinesq formulation are not suitable for a global domain.

A companion paper (Walko and Avissar 2008, hereafter WA) describes several aspects of OLAM’s new dynamic core, including the global grid structure, which consists of spherical triangles, and the spatial and temporal discretization methods as applied to the shallow-water system of equations. The OLAM simulations of five shallow-water test cases described by Williamson et al. (1992) are also described in WA and it is shown that OLAM gives good results for all tests.

The purpose of this paper is to examine the formulation and performance of the OLAM dynamic core with multiple vertical levels, particularly where vertical motion plays a significant role. A central component of the formulation is the spatial and temporal discretization of the compressible nonhydrostatic Navier–Stokes equations, which we describe in detail. Some aspects of the discretization are unique, while certain others have appeared elsewhere but have not previously been tested in atmospheric model applications. This paper also describes how topography is represented in OLAM, which is by means of shaved grid cells (e.g., Adcroft et al. 1997; Marshall et al. 1997) on a height–coordinate grid with horizontal surfaces. We explain particular details of the shaved-cell method that follow directly from the spatial discretization method used in OLAM. We also discuss how momentum is conserved in the C-staggered-grid implementation of shaved coordinates, which to our knowledge has not been addressed elsewhere.

Three types of simulations are performed in order to test and validate the aspects of the OLAM dynamic core that are discussed in this paper. The first is the Held–Suarez (1994) experiment, which replicates a simplified form of the general circulation of the global atmosphere. The second is the two-dimensional density current simulation described by Straka et al. (1993). The third is a series of mountain-wave simulations described in Dudhia (1993) that provide a significant test of the shaved-cell method. These simulations include examples of both nonhydrostatic and quasi-hydrostatic vertical motion.

The remainder of this paper is organized as follows. Section 2 describes the full set of continuous equations that comprise the OLAM dynamic core, plus a general scalar conservation equation that applies to moisture variables and other prognosed fields in most OLAM model applications. Section 3 explains the finite-volume discretization of the governing equations, while section 4 describes the temporal discretization. Representation of topography is covered in section 5. Finally, results from the test simulations are given in section 6.

## 2. Governing equations

OLAM’s dynamic core solves the following conservation equations for momentum, potential temperature, mass, and general scalar quantities:

and the equation of state,

In these equations, **u** and **U** ≡ *ρ***u** are velocity and momentum vectors, respectively; and **g** and **Ω** are the earth’s gravity and angular velocity vectors, respectively. Subscript *i* represents a vector component in the *x*_{i} direction, *t* is time, *p* is pressure, and *θ* is potential temperature. Here *C*_{p} and *C _{υ}* are the specific heats of dry air at constant pressure and constant volume, respectively;

*R*

_{d}and

*R*are gas constants for dry air and water vapor, respectively; and

_{υ}*p*

_{0}is a reference pressure equal to 10

^{5}Pa. Total density,

*ρ*, is given by the sum of densities of dry air, water vapor, and liquid plus ice condensate:

Scalar variable *s* represents the specific density or concentration (relative to *ρ*) of any prognostic scalar quantity, such as various classes of ice and liquid hydrometeors and aerosols, and *s*_{d}, *s*_{υ}, and *s*_{c} are specific densities of dry air, water vapor, and total condensate, respectively. Here, *F*_{i}, *H*, *M*, and *Q* are forcing terms for momentum, internal energy, mass, and scalar fields, respectively. These terms represent processes such as radiative transfer, microphysical phase changes, surface fluxes, and/or optional nudging to observational data, as applicable to each equation.

The ice–liquid potential temperature, described in Tripoli and Cotton (1981) and Walko et al. (2000), is used in OLAM as the prognostic internal energy variable. It has the desirable property of being nearly constant in a parcel for processes of transport and internal phase change. It is empirically related to potential temperature by

where *q*_{lat} is the latent heat required to vaporize any liquid and ice water present and *T* is air temperature. Application of Eqs. (6) and (7) to Eq. (5) gives the following equation for pressure:

where

The coefficient *β* depends primarily on the specific densities of water in all three phases, although a weak dependence on air temperature arises from Eq. (7). Bryan and Fritsch (2004) used simulations of ascending and descending moist parcels to compare the performance of several different definitions of ice–liquid potential temperature, including Eq. (7) and alternatives that have greater accuracy and fewer simplifying assumptions. He concluded that temperature errors resulting from Eq. (7) were primarily due to the neglect of specific heats of liquid and ice hydrometeors. Walko et al. (2000) incorporated Eq. (7) into a combined system of equations for hydrometeor mass and energy budgets that included hydrometeor specific heats. We plan to carry out numerical tests identical to those in Bryan and Fritsch (2004) to evaluate whether the combined system eliminates most of the error that he found.

## 3. Finite-volume discretization

Equations (1)–(4) are discretized using the finite-volume method (Hirsch 1990; Adcroft et al. 1997) by integrating over control volumes Ψ bounded by surface areas *σ* that represent individual grid cells. Control volume integrals of vector divergences are transformed to surface integrals by application of the Gaussian divergence theorem:

where *d*** σ** =

**n**

*dσ*,

**n**is the outward normal unit vector to the surface element, and

*dσ*is a differential element of surface area. In finite-volume form, Eqs. (1)–(4) are written as

In the C-grid stagger implementation, control volumes for scalar variables are prism cells, each having two horizontal and three vertical faces (Fig. 1a), and the momentum component normal to a cell face *j*, *U*_{j} = **U**_{j} · **n**_{j}, is defined and prognosed at each face. As in Wenneker et al. (2002, hereafter WSW02; Wenneker et al. 2003, hereafter WSW03), the control volume for any horizontal momentum component *U*_{i} at vertical face *i* is chosen to be the union of the two prism cells that are adjacent to that face (Fig. 1b). Each of these combined volumes has eight total faces: four horizontal (because the top or bottom surfaces of the joined prism cells are not exactly coplanar on the spherical earth) and four vertical. The control volume for any vertical momentum component *U*_{i} at horizontal face *i* consists of the upper half of the lower scalar prism cell and the lower half of the upper scalar prism cell that are adjacent to that face (Fig. 1c). These volumes have five total faces.

Volume integrals in Eqs. (11)–(14) are represented as products of the mean value in the grid cell (for which we use an overbar symbol) and the volume of the cell:

Surface integrals of fluxes in Eqs. (11)–(14) are discretized as a summation over individual faces of the control volumes, and over each face area, *σ*_{j}, the flux is partitioned into mean and subgrid scale (SGS) contributions:

With application of Eqs. (15) and (16), and the assumption that the pressure gradient commutes with the volume integral, Eqs. (11)–(14) may be written as

All areas *σ*_{j} and volumes Ψ in Eqs. (17)–(20) are precomputed and stored in memory during OLAM initialization. Reduction of their values in grid cells that intersect or are embedded within a topographic barrier is described in detail in section 5.

Advective fluxes *u*_{ij}*U*_{j}, Θ_{j}*U*_{j}, and *s*_{j}*U*_{j} are each the product of advecting mass flux *U*_{j} and an advected quantity. These products are evaluated at each face *j* of a scalar control volume, which is where *U*_{j} is defined. (The one exception is the control volume for vertical momentum, shown in Fig. 1c, to whose top and bottom faces *U*_{j} must be vertically interpolated.) Advected quantities *u*_{i} , Θ, and *s* are defined as means over their respective control volumes and must be interpolated to each face *j* according to the advection algorithm used; an added *j* subscript denotes their interpolated values. WSW02 present both first-order upstream and second-order centered interpolation methods, although their numerical results are based on the first-order scheme. OLAM uses the second-order scheme described in Crowley (1968) and Tremback et al. (1987) for momentum and scalar advection. Next, we will describe the spatial differencing for this scheme. The time differencing is described in section 4.

We first consider horizontal advection of horizontal momentum, which is the most complicated of the advection terms. WSW02 and WSW03 describe the algorithm for two-dimensional planar flow. Arrangement of prism cells in spherical earth geometry requires an extension of the algorithm that links horizontal and vertical velocity components in the interpolation that yields *u*_{ij}, so we provide details for the 3D scheme here. Figure 2a depicts a horizontal computational stencil for the horizontal momentum value defined at location u0. The shaded area is the control volume. Vertical velocity and scalar values at all numbered *w* locations and horizontal velocity values at all numbered *u* locations are required to evaluate the right-hand side of Eq. (17) for horizontal momentum. Eight mass flux values are responsible for transport across the u0 control volume faces. Four are horizontal momentum values defined and prognosed at u1, u2, u3, and u4 locations, and four are vertical momentum values located on the top and bottom faces of the control volume (at the horizontal coordinates of w1 and w2 but out of the plane of the figure.) Transported velocity *u*_{ij} interpolated to each control volume face is defined as the component of total velocity at location *j* that is parallel to the direction of u0, which is indicated by the arrow. This is always different from the direction of *U*_{j} at that location.

Using face u1 as an example, *u*_{ij} is obtained from linear interpolation between both upstream and downstream values, that is, between a value inside and a value outside the control volume. This interpolation has the following form:

where Courant number *C*, which can be positive or negative depending on *U*_{j}, provides the upwind bias of the Crowley scheme. The inside value is the horizontal velocity at u0, which is already parallel to the direction of *u*_{ij}. The outside value is a velocity component parallel to *u*_{ij}, which is linearly constructed from horizontal velocity values at u5 and u6, and vertical velocity values at w3, which are averaged between the top and bottom of the prism cell, according to

Constant coefficients *a*, *b*, and *c* are determined by solving the matrix equation:

where **N**() is a column vector that contains the (*x*, *y*, *z*) components of the unit vector aligned in the positive velocity direction for the given stagger point in parentheses. The unit vector components are defined relative to the earth Cartesian system. In spherical geometry, **N**(w3) has a nonzero projection onto **N**(u0). For specialized limited domain OLAM simulations in planar geometry, projection coefficient *c* is zero, and projection coefficients *a* and *b* are evaluated from a reduced 2 × 2 matrix equation. For efficiency, coefficients *a*, *b*, and *c* are computed once and stored in memory.

We next consider horizontal advection of scalars and of vertical velocity, for which the horizontal stencil and control volume are depicted in Fig. 2b. Advection across the control volume face at u1 is driven by mass flux *U*_{j} at that location, and the advected quantity is interpolated with upwind bias between values at w0 and w1 locations. In the case of advected vertical velocity, the direction of w1 vertical velocity is not exactly parallel to that at w0 due to the earth’s curvature, and horizontal velocities at u4 and u5 are not exactly perpendicular to w0. Thus, the velocity component outside the w0 control volume face that is parallel to w0 is linearly constructed from w1, u4, and u5 velocities in a manner similar to Eq. (22).

Subgrid-scale diffusive transport is evaluated based on *K* theory and gradients of velocity components or scalar values across each control volume face. Velocity components outside each face in the direction parallel to prognosed velocity inside the control volume are constructed in the same manner as for advection. Values of the turbulent mixing coefficient in OLAM reflect the physical magnitude of turbulence intensity and are not required to be artificially large for numerical damping because there is a small but sufficient amount of damping in the advection operator. The model may be run with a zero mixing coefficient.

The Coriolis force term in Eq. (17) is a source term for *U*_{i} but can arise only from wind components **u** that are orthogonal to *U*_{i}. The orthogonal components must therefore be linearly constructed and interpolated from neighboring grid cells. For example, in Fig. 2a, horizontal velocity values at u1, u2, u3, and u4, and vertical velocity values at w1 and w2 are used for constructing the velocity at u0 (including horizontal and vertical components), which is normal to u0. Weighting coefficients for this operation are obtained by inverting a 3 × 3 matrix of unit vector components.

Evaluation of the pressure gradient force term in Eq. (17) in a horizontal direction *x*_{i}, which is parallel to u0 in Fig. 2a, requires the pressure values at barycentric locations w1, w2, w3, w4, w5, and w6 unless the line between w1 and w2 is parallel to u0. A linear weighting of all six pressure values is determined by inversion of a matrix, as described in WSW02.

OLAM’s horizontally unstructured mesh requires one horizontal and one vertical index for each grid cell. Special bookkeeping procedures are required to store the mesh topology. The computational stencils in Fig. 2 form the basis of the bookkeeping. Each *u*, *w*, and *m* point has a unique global horizontal index value, and each keeps a list of its immediate neighbors, which it identifies by names u1, u2, w1, etc. These names are integer constants to which the value of the unique global index for that point are assigned. Thus, the value of the w3 index attached to w0 is the global index of the w3 point. These bookkeeping procedures provide considerable ease of programming of all physical processes. For example, no separate indexing is required for grid cells inside, outside, or on the boundary of a refined mesh area.

## 4. Time integration method

OLAM uses a time-split integration technique in which Eq. (20) and the turbulent diffusion, Coriolis force, and physics parameterization terms in Eqs. (17)–(19) are applied on the regular “long” model time step Δ*t*, and the remaining terms in Eqs. (17)–(19) (e.g., advective transport, pressure gradient force, and the gravity term) are evaluated on a shorter “small” time step Δ*τ*. Time-splitting enables the fast-moving acoustic waves to be represented in the dynamic system with reasonable efficiency. The following small time step terms are evaluated explicitly in time: horizontal advection of all quantities, vertical advection of horizontal momentum, and the horizontal pressure gradient force. This keeps the number of arithmetic operations per time step to a minimum, and avoids the need to solve a global elliptic equation.^{1} Vertical advection of mass, potential temperature, and vertical momentum, the vertical pressure gradient force, and the gravity term are all evaluated with a Crank–Nicolson implicit scheme so that fine vertical grid spacing does not excessively limit the small time step in connection with vertically propagating acoustic modes.

Evaluation of momentum and mass advection terms on the small time steps is novel and distinguishes our scheme from the time-splitting method used by Klemp and Wilhelmson (1978) and also from the third-order Runge–Kutta method (RK3) and Crowley time-splitting schemes presented in Wicker and Skamarock (2002). The recent scheme of Klemp et al. (2007) does include advection of a perturbation mass quantity, but not momentum, on the small time step. It might appear that evaluating momentum advection on small time steps is a poor choice that would decrease computational efficiency. However, for comparison, the RK3 time-splitting method of Wicker and Skamarock (2002) performs three advection evaluations per long time step and requires about 60% more small time steps to be performed for each long time step. Thus, the computational expense of the OLAM and RK3 time-split schemes appears roughly equivalent, assuming identical small and long time steps.

There were two motivations for developing and testing the new scheme. The first was to require that the mass flux, *U*_{j}, which is responsible for advective transport across control volume boundaries, be consistent (equal) between Eqs. (17)–(20), at least in the time average. The other schemes referenced above do not enforce this equality because they evaluate momentum advection on the long time step before all the *U*_{j} values that advect potential temperature (and possibly mass) on the corresponding small time steps are known. This leads to the possibility that individual grid cells with a zero long-term trend of pressure [i.e., with nondivergent time-averaged *θ*_{j}*U*_{j} values applied in Eq. (18)], could use *U*_{j} values in Eqs. (18) and/or (20) that have persistent nonzero divergence. One reason that this situation could arise is that small time step divergences can oscillate between positive and negative values at the frequency of and in response to the updates of the long time step tendency terms. We have observed this behavior in RAMS. If the *U*_{j} values used on the long time step for transporting momentum in Eq. (18) were taken from a single small time step, which happened to be one where the divergence differed from the long-term average, the inconsistency would arise. Although we have not documented a specific case where inconsistency of *U*_{j} leads to numerical problems or inaccuracies, consistency seems a good principle to enforce.

The second motivation for our scheme was the realization that the accuracy of a time-differencing scheme increases with smaller time steps. A similar point was made by Durran (1991) in his comparison of the third-order Adams–Bashforth scheme (AB3) with Runge–Kutta methods. He found that although the latter is more accurate for a given time step length, the lower cost of AB3 per time step enables it to be run at a smaller time step with greater overall accuracy and comparable cost compared with RK3. Thus, a simple time differencing scheme can be competitive with a more complex and accurate scheme by virtue of its smaller time step.

We found that the following scheme produces solutions that are stable for long and short time steps of comparable length to other compressible models such as RAMS, and are also accurate for a range of simulations, as demonstrated in WA and further demonstrated below. For simplicity, we return to the continuous equations, express long time step forcing terms as *G*, assume that Θ = *θ*, consider flow in two dimensions, and separate the momentum equation into horizontal and vertical parts, with *w* and *W* representing vertical velocity and momentum, respectively. Equations that advance momentum, mass, and potential temperature forward one small time step from *τ* to *τ* + Δ*τ* are written as

Coefficients *A*, *B*, and *C* are constants approximately equal to ½, but may be independently increased slightly to provide damping of the solution (in practice, all are usually set to 0.55). Terms that contain these constants in their superscripts thus apply approximately midway through the small time step. Note that advecting mass flux terms *U*^{τ+AΔτ} and *W*^{τ+AΔτ} are identical in all four equations. Horizontal flux *U*^{τ+AΔτ} is evaluated explicitly using a simple forward extrapolation from the current and previous prognosed values:

Vertical flux *W*^{τ+AΔτ}, pressure, and density are interpolated in time between the current and future prognosed values:

Advected quantities *u*^{τ}, *w*^{τ}, *θ*^{τ} that appear on the right-hand side of Eqs. (24), (25), and (27) at the current small time step *τ* are evaluated using the second-order Crowley (1968) method. This is another feature that distinguishes our scheme from the methods of Klemp and Wilhelmson (1978), Wicker and Skamarock (2002), and Klemp et al. (2007), which all use the forward–backward scheme to integrate the acoustic terms. Horizontal and vertical advection are evaluated simultaneously, without including cross-derivative terms, even though the Crowley method is known to be potentially unstable in this form at higher Courant number. However, the solution of advective and acoustic modes on the same time step naturally keeps Courant numbers low for atmospheric flows. Although not implemented here, a monotone advection operator can easily be incorporated into the Crowley scheme. For example, the slope limiter method of Mesinger and Jovic (2003) is an efficient procedure for determining a representation of field values that is piecewise linear within each control volume and possibly discontinuous across control volume boundaries. Modification of the Crowley method to integrate over these piecewise elements produces a monotonic solution, with only modest added cost to the overall scheme even with its frequent evaluation on the small time step.

An implicit equation for Δ*W* is obtained by the following procedure. Equations (29)–(31) are substituted into Eqs. (25)–(27) to give

where the following terms represent explicitly evaluated contributions:

Pressure at the future small time step is obtained by linearizing the equation of state about the thermodynamic state that existed at time level *t*, at the beginning of the current series of small time steps:

Finally, Eqs. (39) and (33) are substituted into Eq. (32), which gives the implicit equation for the change in *W* from *τ* to *τ* + Δ*τ* :

Equation (40) is solved implicitly using a tridiagonal solver. The resulting value for Δ*W* is then back substituted into Eqs. (33), (34), and (39) to obtain updated values for density, potential temperature, and pressure. After these terms are computed, *U*^{τ+AΔτ} is determined using Eq. (24). The order of these operations is different from the other time-split schemes referenced above, in which horizontal velocity or momentum values are updated first, followed by vertical momentum and thermodynamic quantities.

On each long time step, time averages are computed of small time step values of *U*^{τ+AΔτ} and *W*^{τ+AΔτ}, and these are used in Eq. (4) to advect all other conserved scalar quantities. This extends mass flux consistency to transport of all scalars.

OLAM does not generally require special filtering techniques for numerical stability. The inherent damping in the Crowley advection scheme (e.g., Tremback et al. 1987) and the implicit solution of the vertical acoustic terms appear to provide adequate stability, even when horizontal and vertical turbulent mixing coefficients are zero. However, the upper boundary condition in OLAM is a rigid lid, which reflects upward-propagating gravity waves downward, back into the domain. Thus, in some simulations, a Rayleigh friction damping layer is activated in the upper model levels.

## 5. Representation of topography

Topography is represented in OLAM by the volume-fraction or shaved grid cell method (Hirt 1992; Steppeler et al. 2002; Adcroft et al. 1997; Bonaventura 2000) rather than the terrain-following coordinates used in RAMS and most other models (Fig. 3). In the shaved grid cell method, model levels remain horizontal in the presence of topography and intersect the ground surface. Grid cells that are completely below ground are excluded from numerical computations. Cells that are only partially submerged are modified in a manner that represents the kinematic effect of the topographic barrier.

An advantage of horizontal grid mesh surfaces is that horizontal gradient evaluation is simple and efficient because it avoids the transformation terms that are required for terrain-following coordinates. More importantly, horizontal gradient computation with terrain-following coordinates is subject to error because it is inextricably linked to the vertical gradient, which can be much larger (e.g., Sundqvist 1976; Mesinger et al. 1988; Mellor et al. 1994). This is particularly true of the pressure gradient computation. A common method for reducing this error has been to subtract the hydrostatic pressure field of a horizontally homogeneous reference state from the full pressure, and to compute the horizontal gradient of the remaining perturbation pressure field (e.g., Klemp and Wilhelmson 1978; Tomita and Satoh 2004). This approach is most effective in limited-area models where the reference temperature can be chosen close to the actual temperature everywhere in the domain. However, over a global domain, atmospheric temperature varies so widely with location (often spanning a 100-K range near sea level) that any chosen reference temperature must depart by more than 15% from actual temperature at some locations. This means that the perturbation pressure field will have a vertical gradient more than 15% as large in magnitude as the full pressure gradient at some locations, a value that is still orders of magnitude larger than the horizontal gradient. The total separation of horizontal and vertical gradient evaluation in shaved-coordinate models avoids this type of error, and allows full pressure to be used if desired.

A disadvantage of the height coordinates used with shaved grid cells is that high vertical resolution near the ground is more expensive to attain over a wide range of topographic height. Vertical grid levels are commonly spaced more closely near the ground in order to increase resolution, and with terrain-following coordinates, this is accomplished over all topography heights with only a few additional levels. To obtain high near-ground vertical resolution with height coordinates, model levels must be closely spaced over the range of topography height where high resolution is required. If this is, say, up to 4 or 5 km above sea level, the number of levels that must be used can be substantial. A modeler may choose to reduce this expense by requiring high resolution up to only 2 km or so, foregoing enhanced resolution over higher-terrain areas because they have relatively infrequent occurrence globally. The expense of higher vertical resolution is partially offset by the greater efficiency of horizontal gradient computation and the exclusion of underground cells from computation. Thus, a shaved coordinate model can run with more levels than a terrain-following model with equal computational expense.

A disadvantage of shaved grid cells is that representation of flow kinematics and surface fluxes at the ground is more complex than with terrain-following coordinates. Grid cells that intersect topography have modified properties and require special treatment. The method developed for the triangular mesh of OLAM is similar to other methods but differs in a number of details. We describe the scheme next.

The basic shaved-cell method is described in detail in Adcroft et al. 1997. It is a natural application of the finite-volume formulation of transport terms in Eqs. (17)–(20). Grid cell volumes Ψ and face areas *σ* appear explicitly in these terms. Their values are defined such that they represent only the portion that is above ground. Thus, for any grid face that is completely below ground, *σ* = 0, which has the effect of preventing transport across that face. Figure 4 depicts the properties of OLAM grid cells that are partially submerged in the ground. Topography height values are defined at locations that correspond to the vertices of the triangular mesh. Piecewise linear interpolation across each scalar control volume of heights from its three vertices defines a planar topographic surface. The portion of the grid cell volume and each of its five planar faces that lies above this plane is determined, and these values are assigned to Ψ and *σ*_{j}. Thus, the rate of flow through the grid cell surfaces is regulated to the correct amount for the above-ground space occupied by the cell.

In OLAM, the reduction of Ψ and *σ* by topography is assumed to be uniform throughout the control volume and over each face, as depicted in the bottom panel of Fig. 4. This is consistent with defining scalar value *s* at the cell barycenter rather than at the center of the unsubmerged portion. This is important so that *s* remains at the same geopotential as its lateral neighbors and horizontal and vertical gradients remain separable. This contrasts with, for example, the method used by Calhoun and LeVeque (2000) which centers scalar values and fluxes in the unsubmerged portion of the cell volume and surfaces.

The governing equations used in Adcroft et al. (1997) and in OLAM follow conservation principles by casting transport terms in conservative flux form. Shaved cells adhere to these principles in a straightforward way for scalar quantities because every grid cell face that is open to scalar transport lies between two scalar cells in which scalars are fully prognostic. [The shaved-cell method described by Steppeler et al. (2002) differs in this respect because transport, other than for mass, is cast in advective rather than in flux form.] Scalar cells that are not prognostic because they are fully submerged have *σ* = 0 on all faces, so they neither gain nor lose a scalar quantity. The situation is more complicated for momentum transport in a staggered grid system, and additional care is required to enforce conservation. We have not seen a full discussion of this topic in the literature so we provide the details here. We focus on the grid configuration used in OLAM, although the methods we describe are more general and apply to the more common quadrilateral meshes as well.

Consider vertical and horizontal cross sections through the grid, from point *A* to point *B* as depicted in Fig. 5. The topographic surface is represented by the dot–dash line, and scalar control volumes that are completely submerged are shaded gray. Scalar control volumes are labeled with *s*, horizontal momentum surfaces with *u*, and vertical momentum surfaces with *w*. Recall the choice of control volumes for horizontal and vertical momentum (Fig. 1). Thus, the control volume for horizontal momentum *U* at u3 is the combined control volumes for s3 and s4, and the control volume for vertical momentum *W* at w4 is the upper half of s4 and the lower half of s5. We assume for simplicity that horizontal flow normal to the cross section (i.e., across points labeled u0) is zero. Consider first the *U* momentum budget at u3. Its control volume has fluxes across surfaces at w3, w4, u2, u4, w1, and w2. None of these presents any problem because the first three represent exchanges with adjacent control volumes where *U* is fully prognosed (i.e., at u1 and u5) and the last three fluxes are zero. However, a problem occurs in prognosing momentum at u2 and u6. Control volume u2 exchanges momentum across surface u3 with control volume u4, while control volume u6 exchanges momentum across surface w4 with control volume u4. The problem is that momentum cannot be fully prognosed at u4 because one side is submerged. Thus, for example, there is no proper way to compute a pressure gradient force. It may seem that momentum and velocity at u4 should be set to zero because the entire face is below ground, but this violates momentum conservation and also produces noisy, poorly behaved solutions. It must be remembered that u4 is both a surface (whose area is zero), and a control volume, whose volume is greater than zero because of the contribution from s4. As an active control volume, u4 can and does gain and lose momentum via transport with adjacent control volumes. Although momentum in the u4 control volume cannot be prognosed, it can act as an intermediary between u2 and u6. Thus, u4 momentum is semiprognostic, being governed by transport terms alone. What momentum is received from either of its neighbors (u2 and u6) is subsequently donated to the other. Conservation is guaranteed because control volume u4 is identical in this case to control volume s4 and follows the same flux-conservation principles.

A similar issue arises with *w* control volumes, although the solution is somewhat different. Vertical momentum in control volume w3 can exchange momentum with control volume w1, and momentum in control volume w4 exchanges with both w2 and w5, yet w1, w2, and w5 cannot be fully prognostic. We solve this problem by combining w1 and w3 into a single control volume, and also combining w2 with w4 and w5 with w6. Control volume surfaces are correspondingly combined as well. Thus, the control volume for w4 includes all of surfaces u3 and u4 and half of surfaces u5 and u6.

Two further issues needed to be addressed in this formulation of shaved coordinates. One is the commonly discussed problem that the *σ*/Ψ ratio of a control volume can be much greater in partially submerged cells than in unsubmerged cells. Cell s4 in Fig. 5 illustrates an example, with volume reduced by more than 2/3 and surface areas u3 and w4 reduced by less than 1/3. Numerically, this is equivalent to reducing grid spacing, which can place severe Courant–Friedrichs–Lewy (CFL) stability restrictions on the time step. A common remedy used by Hirt (1992), Steppeler et al. (2002), and Bonaventura (2000) has been to increase Ψ without changing *σ* in order to reduce the surface to volume ratio to an acceptable level. While this approach is stable and produces accurate solutions in some cases, there are situations where the solution is unsatisfactory. We illustrate one example in Fig. 6, which depicts the passage of a cold front in a uniform flow (note the use of Ω to represent volumes in the figure). Grid cells in the lowest layer are mostly submerged, so their volumes and vertical face areas are much less than in the layer above (i.e., Ω_{1} ≪ Ω_{2} and *σ*_{1} ≪ *σ*_{2}). If Ω_{1} were artificially increased to guarantee stability, given that the top surface of control volume Ω_{1} is fully open, the temperature would decrease more slowly in Ω_{1} than in Ω_{2}, causing unrealistic retardation of the front in that layer. In this example, the temperature in Ω_{1} would eventually cool to the correct value, however. A more serious example occurs in a steady-state flow that will be described in the next section. For this case, volume reduction is unsuitable. An alternative remedy for small control volumes, used for example by Ye et al. (1999), is to combine partly submerged scalar control volumes with the unsubmerged or less submerged neighbor above. This approach is impractical in the present context because the widely different pressures in the two cells cannot be readily combined. Still another remedy, described in LeVeque (1996) and Calhoun and LeVeque (2000), implicitly balances inward and outward fluxes of small control volumes. This is accomplished by recognizing that transport traverses the entire volume in less than one time step, so the outward flux must be directly tied to the inward flux on the same time step. This approach, although not yet implemented in OLAM, is our chosen method for stabilizing values in small control volumes because it always produces stable and accurate solutions, although at somewhat greater computational expense than the other methods described. For the simulations with topography described below, we resorted to slightly raising topography to submerge the smallest control volumes and to reducing the time step where necessary to avoid instability.

The second remaining issue concerns the use of second- or higher-order advection schemes in connection with shaved grid cells. We found some instances of weak oscillatory behavior in scalar and momentum fields adjacent to the surface. These would be prevented by using a monotone advection scheme, but since one has not yet been implemented in OLAM, we instead increased the upwind bias in the advective algorithm in the lowest two model levels. This successfully removed the oscillations in the solution.

## 6. Test simulations

This section presents the results from three types of simulations, each of which tests aspects of the dynamic core described in this paper. In contrast to tests of the shallow-water equations on the sphere that were presented in our companion paper (WA), vertical motion plays an essential role in the simulations described here. The first test, proposed by Held and Suarez (1994, hereafter HS), is an idealized simulation of the general circulation of the global atmosphere, where vertical motion is quasi-hydrostatic, even though it is modeled here with the nonhydrostatic compressible equations. The second test is the density current simulation of Straka et al. (1993, hereafter S93), which is performed in a 2D vertical–horizontal domain with a limited-area version of OLAM. The third test is a set of six simulations described in Dudhia (1993) of flow over a topographic barrier for different values of the barrier width. These are also conducted with a quasi-2D limited-area domain. These tests are an essential step in determining whether all elements of the OLAM dynamic core, some of which are unique and have not previously been applied in atmospheric modeling, can accurately represent atmospheric dynamic processes. As in WA, these tests are designed to isolate fundamental dynamic processes such as advection, geostrophic adjustment, and wave propagation in an idealized setting. Thus, the tests exclude radiative transfer, all phases of water, and surface fluxes.

### a. Held–Suarez simulation

The HS simulation is run with multiple vertical levels for a deep atmosphere approximating that of the earth. Standard model parameterizations of moisture and radiative transfer are absent, but the global temperature field is forced by Newtonian relaxation toward a specified profile that is constant in time and longitude, and dependent on height and latitude. The profile temperature is isothermal (200 K) above a surface roughly corresponding to the earth’s tropopause with a height near the 200-mb level at the equator and the 400-mb level at the Poles. At lower heights, the specified temperature is symmetric across the equator, and has a surface (1000 mb) value of 315 K at the equator and 255 K at the Poles. Full details of the forcing are given in HS and in Fox-Rabinovitz et al. (1997, hereafter FR97).

The goal of the HS experiment is to generate the model climate, defined as a long-term average of the model state, under the imposed forcing. The imposed baroclinic forcing generates circulations that, under the influence of the earth’s geometry and rotation, develop into flow patterns that resemble midlatitude jet streams and baroclinic waves, the Hadley circulation, and other features of the earth’s atmosphere. These circulations transport heat vertically and latitudinally, which maintains the tropical troposphere somewhat cooler and the polar troposphere somewhat warmer than the specified profile. Thus, the Newtonian relaxation serves to add heat at low latitudes and remove it at high latitudes, similar to the radiation balance of the earth. The effects of surface friction are approximated by removing momentum through Newtonian relaxation below the 700-mb level.

The initial conditions of the HS simulation are arbitrary. We chose zero wind and horizontally homogeneous temperature. The OLAM simulation is run with *N* = 26, which is a horizontal resolution comparable to that used by HS and FR97. Spacing between OLAM vertical levels expands geometrically with height in order to approximate the constant increment between pressure levels used in HS. However, since pressure increments become smaller in the upper levels, a total of 26 levels are used to achieve pressures near zero at the model top. The simulation is carried out for 1200 days and the model climate is fairly well established after the first 200 days. Model statistics are obtained by averaging over the intervening 1000 days and also over longitude. The latter is done by binning triangular grid cells into latitude bands in 2° increments, while linearly weighting each cell between the two latitude bands that bound it. This results in some latitudinal smoothing of the zonally averaged fields.

Figure 7 shows four results from the OLAM simulation. The mean zonal wind (Fig. 7a) contains westerly jets of 30 m s^{−1} located near 40° latitude and 250 mb in both hemispheres. Easterly flow occurs at all heights near the equator, with a stronger easterly jets at the highest levels over the equator and near the surface at about 15°N/S latitude. Weak easterly flow occurs in the polar latitudes. The zonally averaged model atmosphere is cooler at the equator and warmer at the Poles than the imposed temperature field at all heights. The lowest temperature (Fig. 7b) occurs at the equatorial tropopause and is below the imposed value of 200 K. The strongest meridional wind (Fig. 7c) is an equatorward flow near the surface and 15° of latitude, corresponding to the trade wind belts of the Hadley circulation. Return (poleward) flow occurs above 800 mb between 0° and 30° latitude. A weaker, reversed circulation, corresponding to the Ferrel cell, occurs in the midlatitudes. The maximum potential temperature variance (Fig. 7d) exceeds 35 K^{2} and occurs near 40° latitude and 750 mb.

The major features of these OLAM results are very similar to HS and FR97. Locations and magnitudes of the zonal jet streams and meridional circulations are nearly identical. Temperatures are also very similar, although HS obtain a slightly lower equatorial tropopause temperature, and the maximum temperature variance in OLAM falls a little short of the 40 K^{2} peak produced in HS. This is due in part to interpolation from OLAM grid cells to latitude bands for the plot, and also possibly due to damping in the OLAM advection operator, which is less scale selective than the eighth-order Shapiro filter used in HS.

### b. Density current simulation

This experiment, described by S93, is fundamentally two-dimensional and is performed on a domain that is 6.4 km deep and 38.4 km long. To simulate this case with OLAM’s triangular mesh, a section of grid cells was set up similar to that depicted in Fig. 5. A long channel consisting of a single row of alternating equilateral triangles and multiple vertical levels is configured with *σ* = 0 on grid cell faces that coincide with the channel walls (e.g., at locations labeled u0). (This experiment has no topography, so gray shading in Fig. 5 should be ignored for this description.) Horizontal grid spacing for this case is defined as the incremental distance along the channel between successive triangle barycenters. The ends of the channel are normally open in performing this experiment (e.g., S93), but open boundary conditions are somewhat difficult to configure with the triangular mesh so we instead used cyclic conditions between the channel ends. This detail is unimportant because the channel is too long for the end boundary to significantly influence the solution over the simulation period.

The air in the channel is initially motionless and isentropic, with ** θ** = 300 K, except for a cold perturbation centered around a point that is 3 km above the ground and midway along the length of the domain. The cold perturbation is elliptical in shape with maximum horizontal and vertical extents of 8 and 4 km, respectively. Temperature in the center is 15 K lower than ambient temperature at the 3-km level. Subgrid turbulent mixing is performed with a constant imposed mixing coefficient equal to 75 m

^{2}s

^{−1}. Bilateral symmetry in this experiment may be exploited by simulating only half of the domain, but we have not done this. The simulation is conducted for 900 s, during which time the cold perturbation descends to the surface and spreads laterally along the channel (in both directions). The flow becomes highly nonlinear, and strong wind shear at the edge of the advancing cold air pool generates a series of Kelvin–Helmholtz waves. The simulation is conducted with different grid spacings, and the finer grids resolve more of the small-scale turbulent motions. The fixed mixing coefficient prevents turbulent energy from cascading indefinitely to finer scales. The simulation is thus an excellent test of a model’s ability to converge on an accurate nonlinear solution as resolution is increased.

Figure 8 shows contours of the OLAM-simulated potential temperature field at simulation times *t* = 0, 300, 600, and 900 s using Δ*x* = Δ*z* = 50 m. The initial distribution of potential temperature is seen in the top panel. By 300 s, the cold air has descended to the surface and is spreading laterally, with a sharp frontal boundary near *x* = 4 km. Flow at this time is still mostly laminar. At 600 s, shear instability has spawned a coherent eddy circulation centered near *x* = 5 km, and a second eddy near *x* = 9 km is beginning to form. At 900 s, a total of three distinct eddies have formed, and the cold front has advanced to *x* = 15.5 km. Close visual comparison of these results with the reference solution presented in S93 shows that both are essentially identical.

Following the procedure described in S93, self-convergence properties of OLAM relative to grid spacing are evaluated by repeating the density current simulation at different spatial resolutions and comparing their results. A high-resolution reference solution is generated using Δ*x* = Δ*z* = 25 m. All simulations, including the reference simulation, use the same time step Δ*t* = 0.03125 s. Differences between the reference simulation and each lower resolution simulation are obtained in terms of the rms difference in the potential temperature field, also termed the *L*_{2} norm for *θ*. Self-convergence properties relative to temporal resolution are similarly evaluated by repeating the simulation with different time steps and fixed spatial resolution. A separate series is run for grid spacings of 200, 100, and 50 m. These results are plotted in Fig. 9, which may be compared directly with Fig. 4 in S93. It can be seen that spatial convergence error in OLAM is very close to that obtained in S93 with a nearly quadratic decrease rate at the highest resolution. Temporal convergence is close to first order with errors very close to S93 in the case of the 200-m grid spacing. The 100- and 50-m grid spacings produce higher temporal convergence errors in OLAM, in contrast to the results of S93.

### c. Flow over a topographic barrier

This set of six test simulations is described in Dudhia (1993) and simulates the generation and propagation of inertial gravity waves in a deep stable current that traverses an orographic barrier with a bell-shaped cross section. The upstream boundary condition consists of velocity *u* = 10 m s^{−1} at all heights and a potential temperature field in which Brunt–Väisälä frequency *N* = 10^{−2} s^{−1} is uniform with height. An effective Coriolis parameter *f* = 10^{−4} s^{−1} that acts only on horizontal motion is set by aligning the vertical domain axis with the earth’s axis and reducing the earth angular velocity appropriately. Five of these simulations use a barrier height *h* of 400 m, and a barrier half-width *a* equal to 100 m, 1 km, 10 km, 100 km, and 1000 km, respectively. Each of these simulations is run for a period equal to 21.6 *a/u*, which is sufficient for the solution to reach a quasi-steady state. The sixth simulation uses a barrier with height and half-width equal to 1 km. This last case is nonlinear, with *Nh/u* = 1, and produces breaking waves. It is run for 1 h, by which time the maximum tilt of isentropic surfaces has reached the vertical.

Although these six cases are, like the density current simulations of the previous section, essentially two-dimensional flows, those with wider barriers produce a significant crosswind component (parallel to the barrier axis) due to the Coriolis force. Thus, confinement of the flow between two impermeable channel walls is not suitable for the present experiments. Instead, a limited-area, planar configuration of the OLAM grid was developed that provides freedom of movement in both horizontal directions. A section of the grid that shows its full width and a portion of its length is shown in Fig. 10. Cyclic lateral boundary conditions are applied across the sides of the channel, giving the channel an effective width half that occupied by the grid cells, or twice that used for the density current experiments of the previous section. Cyclic boundary conditions are also applied between the ends of the channel, as in the previous section, because open boundary conditions are somewhat complex to implement with the triangular mesh. However, wave disturbances propagate sufficient distance during the orographic simulations to exit one end of the channel and reenter the other, so we use a somewhat longer channel than Dudhia, and in the region near both ends, we apply Rayleigh damping that essentially removes the disturbance before it cycles back into the domain. Rayleigh friction is also applied at heights greater than 12 km in order to prevent gravity wave reflection at the top boundary. Turbulent diffusion and surface drag are zero for these simulations.

Horizontal grid spacing, defined for the channel as in the density current experiment, is set to *a*/5, as in Dudhia. Vertical grid spacing is 100 m at the lowest model level and expands by a factor of 1.05 at each successive level. The 400-m-high topographic barrier penetrates nearly four model levels, while the 1000-m barrier penetrates more than eight levels. These simulations are therefore a significant test of the shaved grid representation of topography.

The steady-state vertical motion fields produced in these OLAM simulations are shown in Figs. 11a–f, and correspond to Figs. 4 and 5 in Dudhia (1993). All results are essentially identical between the two studies in terms of their characterization of each different wave regime and their principal characteristics such as wavelength and phase tilt. Differences are noticeable in some of the finer details. Figure 11a, with the 100-m barrier half-width, represents a case where stationary gravity waves cannot form because flow traverses the barrier at a frequency *u/a* of 10^{−1}s^{−1}, which is much faster than *N*. Strong vertical motion occurs only at low levels comparable to the barrier height. Our results differ from Dudhia in that here, flow separation occurs near the top of the barrier, with subsidence displaced farther downstream and shallow return flow occurring on the lee slope of the barrier. Figure 11b, with the 1-km barrier half-width, shows a case where *u/a* is comparable to *N*, which produces stationary gravity waves of comparable horizontal and vertical wavelengths. Wave energy propagates upward and downwind relative to the barrier. For the even wider barrier of 10 km used in Fig. 11c, *u/a* is an order of magnitude less than *N*, which leads to hydrostatic stationary gravity waves with nearly vertical propagation of energy. With the 100-km-wide barrier used in Fig. 11d, *u/a* is comparable in magnitude to *f*, so the Coriolis force becomes a significant term in the horizontal momentum equation. Inertial oscillations are initiated by the barrier and advect downstream, producing a stationary pattern with alternating regions of subgeostrophic and supergeostrophic flow and, by continuity, upward and downward motion. In Fig. 11e, with the 1000-km barrier half-width, *u/a* is an order of magnitude less than *f*, producing a disturbance that is of low Rossby number. Steady-state flow over the barrier is mainly characterized by conservation of angular momentum (potential vorticity) wherein air columns develop anticyclonic relative vorticity while passing over the barrier because of their reduced height. Air travels to the left (into the plane of the channel) as it ascends the left side of the barrier, and turns to the right as it descends the other side, with crosswind velocity close to 3 m s^{−1}. Dudhia (1993) discusses these cases in more detail and points out their qualitative agreement with linear wave theory. Our purpose here is only to demonstrate that OLAM’s dynamic core produces qualitatively similar results.

In section 5, we discussed the severe CFL limitation that can arise from small shaved grid cells and a simple remedy used by some modelers to relax this limitation, namely, to restore some or all of the cell’s original volume. The inertial gravity wave case shown in Fig. 11d is one situation where this remedy produces a significantly different and incorrect solution. The stationary waves in the lee of the barrier represent a balance of forces in the momentum equation, primarily between the Coriolis force and the downstream advection of momentum. If the volume of shaved grid cells is set to a value larger than their actual unsubmerged volume, the momentum advective tendency is reduced, which shifts the solution to a different equilibrium state. We found in our simulations that if full volume is used for the cells that intersect topography, the phase lines become distorted and make a sudden upstream shift in the lowest few model levels. From a Lagrangian perspective, an air parcel passing through grid cells with reduced surface areas but full volume travels more slowly than a more elevated parcel that is in unshaved cells, while their inertial oscillations have identical frequencies. Thus, the parcel in shaved cells oscillates with a shorter wavelength, leading to the retardation of phase lines near the ground. Thus, it is necessary to resort to an alternate means of avoiding the CFL limitation, such as the implicit flux balance that was described earlier.

## 7. Summary and conclusions

The Ocean–Land–Atmosphere Model (OLAM) is a new global model based on RAMS that is designed to represent all scales of atmospheric flows, including mesoscale and microscale systems typically simulated in limited-area models. Accordingly, its dynamic core has been developed for nonhydrostatic compressible flow, with efficient representation of acoustic modes using numerical time splitting. This paper described the details of the OLAM dynamic core and demonstrated its performance in standard test simulations, with a particular focus on flows with significant vertical motion in a multilevel atmosphere.

The time-split temporal discretization method in OLAM is unique, and enforces mass-flux consistency of advective transport terms in all equations. The spatial discretization method, first presented in WSW02 and WSW03, is implemented and tested here for the first time in an atmospheric model, with modifications for a rotating, spherical earth with parameterized turbulent flow. Representation of topography, which is via the shaved grid cell method with height coordinates, also contains unique elements that relate to OLAM’s triangular mesh and to momentum conservation principles in a C-staggered grid. Specific details of these and other aspects of the OLAM dynamic core have been presented.

The following test simulations were performed in order to validate the new dynamic core: 1) the Held–Suarez (1994) experiment, which simulates an idealized general circulation of the global atmosphere, 2) the density current experiment described by S93, and 3) a series of mountain-wave flow experiments described by Dudhia (1993). OLAM simulation results were compared with previous studies and showed that OLAM successfully represented each case with solutions comparable to other models.

The main purpose of developing OLAM was to produce a global-scale model that has the capability to simulate specific regions at very high resolutions so that interactions and feedbacks between all atmospheric scales can be simulated and better studied. For example, such a modeling capability is particularly useful for the simulation of severe storms (e.g., hurricanes) and/or interactions between regional and global climate. Given that objective, it is essential to demonstrate that this new model is capable of performing the typical tests carried out to evaluate the performance of global scale, mesoscale, or microscale models independently. The tests performed here and in our companion paper (WA) contribute to this demonstration. Another essential capability is to run large simulations on distributed memory massively parallel computers, and to this end, an (Message Passing Interface) MPI-based domain decomposition of OLAM was recently implemented and is undergoing testing and evaluation.

## Acknowledgments

The development of OLAM was supported by the Edmund T. Pratt Jr., School of Engineering, Duke University, by NASA Grant NAG5-13781, and by the Gordon and Betty Moore Foundation. This paper benefited significantly from the anonymous reviewers’ comments.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Robert L. Walko, Department of Civil and Environmental Engineering, Duke University, P.O. Box 90287, Durham, NC 27708-0287. Email: robert.walko@duke.edu

^{1}

Explicit differencing does, however, usually limit the time step to smaller values than implicit methods. We point out, on the other hand, that longer time steps can have their own drawbacks. We give the example of rain falling in a model that uses 15-min time steps. The fall distance per time step is a few kilometers, which is at least a few model levels. This could prevent raindrops from participating in microphysical processes, such as evaporation, in the layers skipped during the time step. We thus favor a method that uses smaller time steps and fewer operations per time step.