Scaling and computation of smooth atmospheric motions

We introduce a general scaling of the inviscid Eulerian equations which is satisfied by all members of the set of adiabatic smooth stratified atmospheric motions. Then we categorize the members into mutually exclusive subsets. By applying the bounded derivative principle to each of the subsets, we determine the specific scaling satisfied by that subset. One subset is midlatitude motion which is hydrostatic and has equal horizontal length scales. Traditionally, the primitive equations have been used to describe these motions. However it is well known that the use of the primitive equations for a limited area forecast of these motions leads to an ill-posed initial-boundary value problem. We introduce an alternate system which accurately describes this type of motion and can be used to form a well-posed initid-boundary value problem. We prove that the new system can also be used for any adiabatic or diabatic smooth stratified flow. Finally, we present supporting numerical results.


Introduction
Typically, the hyperbolic systems which describe atmospheric and oceanographic motions contain solutions with more than one time scale.For example, the shallow-water equations govern two classes of motion with different time scales: "slow" Rossby-type motions and "fast" inertia-gravity motions.Meteorologically the main interest is in the first type of motions.Therefore several methods have been developed to filter out the fast waves.In the method called initialization one keeps the underlying equations but prepares the initial data so that the fast waves will not be excited.Alternatively the underlying equations are changed to a set called the reduced system which only allow the slow motions.
I Sponsored by the National Science Foundation.'The work of this author was supported in part by NSF under Grant ATM-A201207 and by the Office of Naval Research under contract N00014-83-KO422.
In the early days of numerical weather prediction, the second approach was used and the reduced system for large scale flow in the middle latitudes was called the quasi-geostrophic system (Phillips, 1956).However the quasi-geostrophic models overpredicted the westward movement of ultralong waves.Phillips (1963) showed that the assumption of small horizontal divergence used in the derivation of the quasi-geostrophic system is inappropriate for ultralong waves and asserted that this was the reason for the poor forecasts of those waves.To relax the assumption of small horizontal divergence, Smagorinsky ( 1963) developed a model based on the primitive equations, i.e., the Eulerian equations modified solely by the assumption of hydrostatic equilibrium.Currently most large scale weather models use the primitive equations with the nonlinear normal mode initialization method (Baer, 1977; Machenhauer, 1977).This can be considered as a partially reduced system in conjunction with initialization, i.e. the hydrostatic assumption eliminates the fast vertical sound waves and the nonlinear normal mode initialization procedure prepares the initial data so that the remaining fast waves which may have large amplitudes will not be excited.
Experience has shown that the primitive equations are not totally satisfactory for limited area models.Analytic results (Oliger and Sundstrom, 1978) now indicate that the use of the primitive equations in a limited area leads to an ill-posed problem.In recent years (Kreiss, 1980;Tadmor, 1982) a mathematical theory for symi metric hyperbolic partial differential equations with different time scales has been developed.The theory, which we briefly review in Section 2, provides a way to initialize the data for initialboundary problems and allows one to derive reduced systems which are automatically wellposed for the limited area problem.By applying the theory to the inviscid Eulerian equations, we have developed an accurate and well-posed limited area forecasting model.
The first step of the theory requires that we scale the system of equations for the motion of interest.
In Section 3, We introduce a general scaling of the inviscid Eulerian equations which is satisfied by all members of the set of adiabatic smooth stratified atmospheric motions.This allows us to consider the entire set of motions with a single scaled system.Then we categorize the members of the set into mutually exclusive subsets and apply the bounded derivative principle to determine the specific scaling for each subset.In Section 4, we consider the subset which includes midlatitude flows which are hydrostatic and have equal length scales.We prove that such motions must satisfy either the scaling introduced by Charney (1948) or the scaling discussed by Phillips (1963).Traditionally the primitive equations have been used to describe these flows.In section 5, we give a simple demonstration that the use of the primitive equations for limited area forecasts of these flows leads to an ill-posed problem.Unfortunately the bounded derivative principle can not be applied directly to the scaled Eulerian system for these motions because of the extreme skewing of the equations (Browning and Kreiss, 1984).However we introduce an alternate system which describes these motions accurately and which can be used to form a well-posed limited area model (Browning and Kreiss, 1984).In Sections 6 and 7, we show that the modified system can also be used for adiabatic hydrostatic motions with unequal horizontal length scales and adiabatic nonhydro-static motions.Diabatic effects are considered in Section 8. Section 9 contains supporting numerical results.

Mathematical background
Consider a hyperbolic system equations of the form where m 1 is a natural number and E > 0 is a small constant.u = (ul, ..., u , ) ~ is a vector function depending on x = ( x , , x2, x,) and t.D = diag ( d l , ..., d,,) is a positive definite diagonal matrix which depends smoothly on x and t.The operator is a linear first order skew self-adjoint differential operator with matrix coefficients depending smoothly on x and E. Finally P I is a quasilinear first order skew self-adjoint differential operator whose coefficients depend smoothly on x, t , and u.
Problems of this kind have solutions on different time scales: fast scales of order E-', 6 < m, and a slow scale of order unity.We are only interested in the solutions which evolve on the slow time scge.In a number of papers (Kreiss, 1980;Browning et al., 1980;Browning and Kreiss, 1982) we have shown that these slow solutions can be obtained using the following simple principle.

Choose the initial data in such a way that at t = 0 time derivatives of the solution of order less than or equal to k are of order unity
We have proved that this principle leads to solutions which only vary on the slow scale on a time interval 0 < t < T with Tindependent of E. The size of Twill in general depend on k.The larger the number of derivatives of order unity at t = 0, the longer it takes for the fast waves to appear and to pollute the slow solution.The number of linearly independent constraints which the initial data have to satisfy so that the time derivatives are of order The constraints consist of a number of elliptic differential equations which the initial data have to satisfy.As long as the solution stays on the slow times because these constraints are satisfied whenever the time derivatives are of order unity.This enables us to derive reduced systems which describe the slow solutions accurately by replacing part of the original hyperbolic system by the elliptic equations.
adiabatic exponent, and k = (0, 0, 1)' is the unit vector in the vertical direction.The total differential operator d/dt is given by d a a a a a -= -+ y. v = -+ We assume that f is given by the tangent plane approximation = 2R @ , + (ylr) cos @ , , where 2R z s-' is the earth's angular speed, 0, is the latitude of the coordinate origin, and r z 10' m is the radius of the earth.This assumption clearly defines the role which the variation of the

Scaling of the basic equations
We want to investigate the adiabatic smooth stratified motions of the atmosphere which are essentially hyperbolic, i.e. those which satisfy the inviscid Eulerian equations.As we have stated earlier, we locally introduce new dependent and independent variables so that the new dependent variables and their first derivatives are of order unity.We stress that after the scaling we assume that the new dependent variables and their first derivatives are of the same size, i.e. some are not smaller than others.We also assume that if the first derivatives are of order unity then the same is true of the derivatives for several orders higher.In particular this assumption excludes the occurence of jump discontinuites in the first derivatives.When such weak solutions are present dissipative forces cannot be neglected.
The adiabatic inviscid Eulerian equations in Cartesian coordinates x, y , and z directed eastward, northward, and upward, respectively, can be written as (e.g.Kasahara, 1974) where t is time, V = ( u, u, w)' is velocity, p is density, p is pressure, and s = pp-l'vis proportional to the reciprocal of the potential temperature.Also, f = f (y) is the Coriolis parameter, g = 9.8 m ss2 is the constant gravity acceleration, y = 1.4 is the Coriolis parameter plays in the equations to follow.
We shall now introduce dimensionless variables to identify the relative magnitude of all terms in the equations.We change independent variables via the relations where L , , L,, D, and Tare the representative scales along the x, y , z , and t axes, respectively.We also change dependent variables.For the velocity we introduce the relations (1 + S,P'/P,) (1 + S,p'/po)-''v = R, P;;"" S ~( Z ) (1 + S, s'), (3.5) where p'(z) = D,(ln p,), with Do = lo4 m.The factor of Do ensures that p' is of order unity.We have also scaled gravity as g = Gg' (G = 10 m s-,). (3.8) We scale the Coriolis parameter in the form We assume that (3.10) so that for each equation the dimensionless time derivative will be of the same order of magnitude as the horizontal advection terms.Typically there is some cancellation between the terms of (In so)z = s; ' (sJz = p;' (p,J, -y1 po'(po),, e.g. in an isothermal atmosphere (In so)z = -( 1y -l ) 0;' z -0.3 D;,.The factor of 10 Do on the right-hand side of the definition of S is to ensure that S is of order unity.For the lower atmosphere typically -1.4 < p' < -1.3 and -4 < S < -1 (Gates, 1960) Any variable in (3.12) which depends on Po also depends on R, in such a manner that it is the ratio P,R;' which is the important value.The ratio of the horizontal mean of the pressure over the horizontal mean of the density is approximately independent of height so that choosing surface values for Po and R, is no restriction.
To simplify the ensuing presentation, we will drop the O ( S , ) terms in (3.11)   = Do(~npo),.
We want to show that D < D o .Assume that D 9 Do.Then the variation of the background terms such as Z relative to the variation of the motion would be large, i.e., if we differentiate a background term with respect to z' the derivative would be larger than unity.This leads to an immediate contradiction since if we differentiate (3.13a) with respect to z' there would be a solitary large term, namely (10 S,)-l S, Z,, w'.This result is not surprising as it is well known that the largest equivalent depth of the atmosphere, which meteorologists call the external mode, is on the order of Do.We also want to show that S, < 1.If S, 1 then it follows from (3.13b) that yw:, + S2 S;I j w ' = 0, i.e. to first approximation w' = 0 in the case of no topography.Thus we did not use the right scaling for w.We summarize our results as In the remainder of the presentation we will drop the prime notation with the understanding that all variables are dimensionless.

Adiabatic hydrostatic motions with equal horizontal length scales
In the current and following two sections, we will consider adiabatic motions of the atmosphere which satisfy S, S, % 1.We now show that such motions are hydrostatic, i.e., that S, = 1.From Section 3, we know that S, < 1. Assume that S, 4 1.Then by (3.13e) to first approximation 4, = 0, i.e. we would not have made the correct vertical scaling.Thus necessarily S, = 1, the depth of the motion is 10 Km ( D = Do), and S, = 9,.To first approximation the hydrostatic relation holds.The fact that the second term needs to be included in (4.1) has also been shown by White (1977).
Next we want to show the well-known fact that the hydrostatic relation (4.1) is much better in the degree of approximation than quasi-geostrophy or semi-geostrophy, i.e. we want to show that S, S, % S, 9, .In order that there not be a solitary large term in the entropy equation, the following must hold : (10 SJ-1 s, < 1 which we can also write in the form s, 2 10-1 s,. By the definitions in (3.12),
(4) The motion is genuinely three dimensional.We want to show that such a motion must be quasi-geostrophic.Let us assume that this were not the case, i.e. that S, < 1 or S, < 1.If S, < 1 while S, % 1, then (3.13~) and (3.13d) require that we rescale u and v.The rescaling leads to a quasigeostrophic system.If S, < 1 while S, % 1 then to first approximation p , = p y = 0 which would say that we did not remove the horizontal mean of the pressure correctly.Thus we can assume that both S, and S, are less than or equal to one.Eq. (3.13b) can be written  If S, < I, then we can neglect the terms on the right-hand sides of (4.5a) and (4.5b), i.e. to first approximation the motion satisfies the barotropic vorticity equations.Since the barotropic equations do not contain any z derivatives, the motion in the plane z = C, is independent of the motion in the plane z = C, where C, and C, are any nonnegative constants with C, # C,.If the only three dimensional smooth solutions were solid body rotations, then such layered equations would be adequate.From observations we know that the motions are genuinely three dimensional.For initial conditions of a smooth solution which does not represent a solid body rotation, the layered system would yield z derivatives that would become large because of shearing and a bound on the z derivatives of the dependent variables would not exist.This shows that S , = 1.With this value for S,, eq. ( 4.2) requires that S, 2 lo-'.From the definition of S , which contradicts our assumption on the maximum size of the velocity.
We have proved that necessarily S,* 1 and S4 % 1.For the proper balance of terms we must also have that S , = S,.In combination these relations say that smooth solutions of the adiabatic system (3.13) which are in the midlatitudes, are in hydrostatic balance, and have equal horizontal length scales must be quasi-geostrophic.Of course it is well known that there are smooth extratropical motions which are hydrostatic and have equal length scales which are not quasi-geostrophic.There are only two assumptions in this section which can be altered to account for these motions.We could relax assumption 4 to permit the motion to be two dimensional.The other possibility is that the motions are driven by heating.We will discuss diabatic effects in Section 8.
If L B r, then (4.5b) would say that we did not scale u correctly so that L must be in the range lo6 < L < 10' .Let us now consider the case that L = 10' m.Then f, = O( 1) and from (4.1 I) we must have that S, = 1.From (3.4) and (4.2) lo-' 2 S , 2 10-1, i.e., S, == LO-'.Substituting this expression for S, into (4.10)leads to the value n = 2.By the same arguments as above we also find that T = lo6 s, U = 10 m s-,, W = lo-, m s-l, and q = 8 = lo-*.This is essentially the scaling discussed by Phillips (1963).
As we have noted earlier, we must have that (10 S , ) -, S, < 1.In both of the above scalings (10 S , ) -, S, = 1.For hydrostatic solutions if (10 S,)--' S, < 1 then to first approximation we could drop the lower term in (4.6a) and the resulting system could have extremely large exponential growth since the large lower order terms would not necessarily be skew symmetric.To avoid this possibility from now on we will assume that for hydrostatic solutions (10 S,)-l S, = 1.We now want to consider what must happen to smooth solutions of (4.6) as we approach the equator, i.e. f -, 0. For f = O(E) equations (4.6~) and (4.6d) require that we introduce the new variable # = E(' (Browning et al., 1980) so that we must reduce S, by one order of magnitude.But then the entropy equation (4.6a) requires that we make the change of variable w = EW'.For the case that n = 2 (4.1 1) becomes and the change of variable implies that 6 = O(E).Since f = O ( E ) and f, = 0(1), this means that we would not have scaled u correctly.When n = 1 the change of variable implies that 6 = O ( 2 ) which leads to the barotropic vorticity equations.Thus, for the adiabatic hydrostatic equal-length scale case, there are no smooth solutions at the equator unless Tellus 38A (1986).4 they are layered.It is conceivable that a smooth disturbance could cross the equatorial region in a sufficiently small time interval so that the growth of the derivatives due to shearing would be small.Although both (3.13) and the system we will introduce in the next section allow such a mathematical possibility, we believe that diabatic effects, which are discussed in Section 8, are more reasonable to describe smooth solutions near the equator.

Initial-boundary value problems for the scaled system
In this section, we want to consider initial boundary problems for the scaled system we derived in Section 4. As we have seen earlier q < E and therefore a reasonable assumption seems to be to use the primitive equations derived from (4.6) by setting q = 0. Oliger and Sundstrom (1978) observed that the initial-boundary value problem for the primitive equations is not well posed.A simple demonstration of this fact can be obtained by studying the initial-boundary value problem for the two dimensional constant coefficient version of --and U > 0 is a constant.We have dropped the undifferentiated terms that are not essential in reaching the conclusion we are seeking.We consider (5.1) in the domain 0 < x < o0,O < z < 2n, and t 2 0 with the solution being 2nperiodic in z.It would be more realistic to assume that w = 0 at z = 0, 1.But then we would develop s and w into sine series and the remaining variables in cosine series.
Thus the assumption of periodicity in z is no restriction.Fourier transforming (5.1) with respect to z gives us where the hat notation indicates the Fourier transform of the corresponding variable and m is the real dual variable of z.Eliminating i~ and f we need only consider the system For a first order hyperbolic system of this form, the number of positive or negative eigenvalues of the matrix L determines the number of characteristics entering or exiting the region.Thus there are two characteristics entering the region if and otherwise there is only one.This shows that the number of boundary conditions is determined by the vertical wave number and that one cannot give a fixed number of boundary conditions in physical space if the problem is to be well posed.
By setting v = 0 the speed of the vertical soundwaves is increased to infinity.Instead we will now slow these soundwaves down, i.e., we shall consider the modified system --_ - The choice of the new vertical sound speed is determined by the requirement that the solutions of the new system be smooth up to the boundary (Browning and Kreiss, 1984).The easiest way to see what value to choose for a is to note from (5.2b) that Applying the operator dldt to this expression we find that By neglecting the term of order E~ we obtain an elliptic equation for 4. In the smoothness argument it is crucial that the coefficients of the elliptic equation be of order unity in order to estimate 4 correctly.Thus we must choose We first consider the case that n = 1.We want to show that any smooth solution of system (4.6) with this scaling can be computed accurately by using (5.2) with a = c4.Let a subscript 0 on a dependent variable denote a given smooth solution of (4.6).
Then we write the solution of (5.2) in the form (5.4) Note that since v < lo-' e2, F is of order unity.
Substituting the expressions (5.4) into (5.2) we find that the perturbation variables with subscript 1 satisfy the linear system We have dropped a number of undifferentiated terms which have no influence on the arguments to follow.To further simplify the presentation we freeze the almost constant coefficients and consider the system dw1 E2 -+ # l z + s, =o.

dt
(5.8b) ( 5 .8 ~) (5.8d) (5.8e) (5.6a) (5.6b) ( 5 .6 ~) (5.6d) (5.6e) To complete the proof of our first assertion it will suffice to show that the solution of (5.6) with initial data equal to zero is small.To obtain the estimates we need we first symmetrize (5.6) by introducing the new variables and a similar estimate for all of its derivatives.
Finally we can use the elliptic equation for 4 that we derived earlier to obtain estimates of the same form for 9, and all of its derivatives.

If we were only interested in computing motions
with n = 2, then we could choose a = E6 and repeat the argument we used for the case n = 1.However in practice both motions can be present and the question arises as to what will happen if we also use the value a = c4 when n = 2.In that case a-I q = E' and by using essentially the same proof as for the case n = 1 we can show that (5.2) with a = 6 and n = 2 approximates solutions of (4.6) with L = 10' m up to an error term of order 2.
The modified system (5.2) has very nice mathematical properties.Since it is essentially a symmetric hyperbolic system, appropriate boundary conditions can be chosen for (5.2) so that the initial-boundary value problem will be well-posed and the solution will be smooth up to the boundary (Browning and Kreiss, 1984).For example on the western boundary one can choose to give s, u +

~' ' ~( y p ~p o ~) -~' ~
4, u, and w at inflow points and u + ~"~(yp~p;')-"~ 4 at outflow points.We note that the use of these same boundary conditions with (4.6) leads to solutions which are not smooth up to the boundary because of the extreme skewing of the system.Thus it is fortuitous that we can show that the smooth solutions of (4.6) can be computed by using (5.2), a system which is well understood mathematically.

Adiabatic hydrostatic motions with unequal horizontal length scales
In this section, we will again consider adiabatic motions of the atmosphere with S, S, B 1 but will change Assumption 1 of Section 4 to be (1) The motion has unequal horizontal length scales From Section 4, we know that the vertical depth scale of such a motion is 10 km so that S, = 1 and S, = s,.We rotate the original coordinate system about the z axis through the angle p where p is the angle between the x axis and the direction of the longer part of the motion.The corresponding transformation of independent variables is given by As is well known, the equations of motion are invariant under such a rotation and the scaled equations in the new coordinate system are the same as (3.13) except that x, y , u, and u are everywhere replaced by 2, Y, li, and 6, respectively.
Of course in the new coordinate system f must be considered to be a function of both x' and 9.In the rotated coordinate system the unequal length scales are represented by L , B L,. From (3.10) the magnitudes of the horizontal velocity components satisfy a similar relation, i.e., U % V. We relax Assumption (3) to be (3) The maximum transverse velocity is small compared with 100 m s-l ( V < IOOms-').
We will now show that this set of assumptions sufficiently determines the scaling of the motion to allow us to prove that the modified system we L , # L,.
introduced in the previous section can also be used for the motions considered in this section.
We want to show that the motion is geostrophic in at least one direction.We first assert that the scaling of the motion satisfies 9, % 1.Let us assume that S3 < 1.Then from (3.12) and the fact that U 3 V we have that S, < 9, < 1.Since we can not have a single large term in any equation necessarily S,, 9, < 1.From (4.5) if S, < 1 the motion would to first approximation satisfy the barotropic vorticity equations.Therefore by the same argument as in Section 4, S, = 1.With this value for S,, (4.2) requires that s, 2 lo-'.
From the definition of s', Po($, R,J-']"* 2 100 m s-1, which contradicts our assumption on the size of V.
Thus we have proved that necessarily s, B 1. If 9, < 1 then to first approximation q$= 0, i.e. we would not have made the correct scaling in the y'direction.
Therefore g4 > 1 and necessarily 9, = f, which means that to first approximation the motion must be geostrophic in at least one direction.
= S, = S,, E -~ = g3 = 9, (0 < a < b), and q -l = S, S, with E as in Section 4 and 0 < t] < lo-, Then in terms of E and q we can write the rotated system as We now define Note that there can be solutions of (6.1) such that S, = 1 so that the horizontal divergence is of order unity.For example this occurs when a = 0 and f f i A-.Q 0 (E).Also from (4.3) S, S4 = 2 n 2 T 2 and if we assume the same time scale as in Section 4, then a + b = 2 .
A specific skewed flow satisfying (6.1) with a = f and b = 3 was considered earlier by Dickinson (1968).
We can use a proof similar to the one in Section 5 to show that smooth hydrostatic solutions with unequal length scales can be computed with only an error of order t2 by using the modified system ds dt -+ F(z) w = 0, (6.2a) (6.2b) d4 s, -+ p0p;"yd + s, p'(z) wl = 0, dt du' If we rotate the dimensional version of (6.2) through the angle -p, then we obtain the dimensional version of system (5.2),i.e. the same modified system works for the motions described in Section 4 and for the motions discussed in this section.

Adiabatic nonhydrostatic motions
In this section, we will consider the remaining subset of adiabatic smooth stratified solutions of the atmosphere, i.e. nonhydrostatic motions with s, s, Q 1.
(7.4) Eq. ( 7.4) requires that T Q ( IODi/Po)"2 = 10' s so that the time scale of adiabatic and hydrostatic motions must be less than two minutes.If we require that max (U, V) Q 100 m s-', then max ( L l , L 2 ) = T max (U, V) < lo4 m, i.e. the longest length scale of the motion is less than 10Km.
We now want to know whether or not the motion is geostrophic.It suffices to consider the case that L , 2 L,.The g4 s 1 if and only if L , L;I 9 ( 2 n T)-' 2 10'.
If L , = L,, we have an immediate contradiction.
Then necessarily s4 Q 1 and since S, < 9, the motion cannot be geostrophic, i.e. adiabatic nonhydrostatic motions with equal length scales are not geostrophic.If L , > L, then the ratio of the length scales has to be much greater than 100 in order for the motion to be geostrophic in at least one direction.
In order that the motion does not become 2dimensional in the equal-length scale case, the following must hold: for these motions as the system has only a small amount of skewing.

Diabatic motions
In this section, we want to discuss the role of heating in the determination of smooth solutions of the atmosphere described by the inviscid Eulerian equations.Heating appears as two forcing terms added to (3.1 1): one in the pressure equation and one in the entropy equation.We write the forcing term in the entropy equation as ds dt -+ (10 S,)-' S, S' W = Ho S'h, (8.1) where h(x, y, z, t) is a specified smooth function.If the heating is of order unity or less, i.e., if H, < 1 then it is easily seen that the arguments of the previous sections are unchanged.However, if H, 9 1, then the arguments of the previous sections involving the entropy equation are no longer valid.
But then we must have that H, = (10 S,)-I S, and to first approximation ~= h , (8.2) i.e., to first approximation the vertical velocity w is equal to the forcing term h.Eq. ( 8.2) should only be considered to be valid above the surface of the earth.At the surface if the heating is large (H, 9 1) (8.2) leads to an apparent contradiction since w = 0 at z = 0.The apparent contradiction is resolved by realizing that we have not taken the boundary layer into account.The boundary layer adds additional terms to (8.1) that can balance the large heating term near the surface.
If one restricts oneself to motions where (8.2) is valid, then one only need to solve the simple reduced system -+ S , , , + S , f u = O Thus in essence one is only solving the reduced system of the shallow water equations (Browning and Kreiss, 1982).The entropy can be obtained from the vertical velocity equation dt ax ay az h, + uh, + vh, + s, hh, + s, s5[#z + S,(O.l f# + gs)l = 0 after the solution of (8.3) is obtained.In the next section we will show that the balance indicated in (8.2) is present in small scale anelastic models.This means that for three dimensional small scale models one need only to solve two dimensional elliptic equations rather than a three dimensional elliptic equation for 4. Also, when H , 9 1, to correctly compute the balance between w and h in (8.1) in a standard anelastic model requires considerable numerical accuracy.If instead one uses (8.3) then no large terms appear and less accuracy is required to compute the same motion.
In practice the heating is computed by a parameterization process, i.e., h also depends on the dependent variables.If Ho < 1 and if h was only a function of the independent and the undifferentiated dependent variables our previous arguments would be unaffected.Of course if the parameterization also involves derivatives of the dependent variables, the system would have to be analyzed taking into account the form of the heating.If H, 9 1 and h is also a function of the dependent variables we would not be able to use (8.3) and again the system with heating would have to be analyzed.

Numerical results
The analysis of the previous sections suggests that the approximate system should be a viable candidate for a limited area forecasting model.To see how a model based on the approximate system behaves in an actual forecast situation, we have made a number of real data forecasts using the ECMWF FGGE Level 111-b analyzed data set of January 13, 1979 which has a grid increment of 1.875O (n/96 radians) in latitude (8) and longitude (A).The data set contains the height, the horizontal components of velocity u and v, and vertical velocity w on 15 constant pressure surfaces (1O00, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20, and 10 mb) at OO:OO, 06:00, 12:OO, and 18:OO Greenwich Mean Time (GMT).For a forecasting window we chose from 30° to 60° latitude and from 260° to 290° longitude, i.e. the eastern portion of the United States.We chose this area because of the high density of observing stations and relatively flat terrain.
For a static check of the data we computed the horizontal divergence 8 = (r cos  recompute u and u.This is a reversible process in the continuum and should be approximately rever--sible for a numerical method applied to smooth solutions.However the error was quite large.Thus we knew that the data contained energy in the high to the global data 7 times.The number 7 was chosen so that the composite stencil of the filter when centered on any grid point of the window would not contain the north pole.This simply avoided additional coding and was adequate for our purposes, i.e. after the filtering the process mentioned above was essentially reversible.Since the data set was given on constant pressure surfaces, we interpolated the data to constant height surfaces at 1 Km intervals up to 18 Km by linear interpolation in u, v, and log p.We obtained the vertical component of velocity from If the inviscid momentum equations accurately describe the large scale dynamics of the atmosphere (which is generally believed to be true) and if there were no error in the observations or analysis, then system (9.1)should give a perfect forecast.We approximated (9.la) at the interior grid points of the window by replacing each partial derivative by its corresponding centered second-order difference approximation and averaged 4 in the lower order term & in time.We assumed that the vorticity was known at inflow points on the boundary.At an outflow point of the boundary we approximated the normal derivative of C by a stable approximation (Browning and Kreiss, 1982).We also assumed the normal velocity component was known on all boundaries which was sufficient to allow us to solve for u and t, from c and 6 using the Helmholtz relations.We solved the elliptic equations for u and v by using the direct methods of Swarztrauber and Sweet (1975).As might be expected from the scalings of Sections 4 and 6, the terms in (9.la) containing derivatives of p can be neglected without affecting these errors.
Also for this case the forcing term f 8 was more important than the terms involving W. This is to be expected for motions described in Section 4 and motions of Section 6 where a > 0. There is a question of how much of the error is due to observational and analysis error and how much is due to the numerical approximations.We have also run the model with twice the number of grid points with essentially no change in the error.Thus we believe the error in the forecast is due to observational and analysis errors.We call this experiment the control run because we cannot expect for any model using the data to do any better than the control experiment unless the heating is the dominant force in a short term large scale forecast.Then if an accurate parameterization of the heating existed the model might be able to do better than (9.1).However we shall see that while the heating is important to the forecast of the large scale, it is not dominant.
Our second forecast used a spherical coordinate version of the modified system introduced in Section 5. We approximated (9.2) at the interior grid points of the window by replacing each partial derivative by its corresponding centered second-order difference approximation.At inflow points on the western boundary we assumed that s, u + (ypopi1)-1/2 #, u, and w were known and updated the outflow Although these errors are quite close to the errors of the control run, we also wanted to attempt to include heating in the model as heating is present in the atmosphere.We computed the "atmospheric heating" from the equation CH = S, + r-l[(cos 6)-l rig, + 6So1 + WS, + fW.We then ran the modified model with the forcing term S'H added to the right-hand side of (9.2a) and the term -ypo S'H added to the right side of (9.2b).Fig. 4 shows the (a) zonal and (b) meridional components of the wind computed by the "diabatic" model at the same height and time as in We also developed a standard two-dimensional anelastic model to show that the balance between the vertical velocity and heating indicated in Section 8 is present.The two dimensional anelastic version of (7.5) is given by HCy, z, t) is the rate of heating per unit mass, c,, is the specific heat for dry air at constant pressure, and T,,(z) is the vertical profile of the horizontal mean of the temperature.We obtained eq.(9.3d) from the continuity equation instead of the pressure equation.Dropping the time derivative term in the continuity equation results in the same error as dropping the time derivative term in the pressure equation.However the continuity equation is simpler than the pressure equation when heating is present.
We now derive the elliptic equation for 4. For simplicity we will assume an isothermal mean state with pressure at the surface given by the standard surface pressure (1.01325 x lo5 kg m-' s-*).Then j, P; and f a r e constants.Differentiating (9.= I,, + + @'+ s') + F( + [gs, + Kijs -@'w)2 (9.5) We consider system (9.3) with (9.5) in the region R = ( ( y , z, f ) : O Q y Q Y, 0 Q z Q Z, f 2 0).We assume that all dependent variables are Y-periodic.At the ground (z = 0) the correct physical boundary condition is w = 0.At z = Z we also assume w = 0.Although this is not physically correct, we will choose the heating so that the disturbance is contained in the lower part of R and the upper boundary condition will pose no problem.
Then at the bottom and top of R we have from (9.3c) fbz + 2fb == --gs, (9.6) which provides the necessary boundary conditions for (9.5).
To discretize the initial-boundary value problem, we choose spatial increments Ay = Y/J and Az = Z / K , where J and K are natural numbers.The temporal increment At is determined from the stability requirement.We approximated each of the partial derivatives in (9.3) by its corresponding second-order centered finite difference approximation.Assuming the grid function s, v, and w, are known at two time levels we can solve (9.5)-(9.6)for ( using the direct methods of Swarztrauber and Sweet (1975) if we replace each of the partial derivatives on the right-side of (9.5) by its corresponding centered second-order finite difference approximation.Then we can solve the finite difference approximations of (9.3) for s, v, and w at the next time level.We are then in a position to be able to repeat the entire procedure.The values of the parameters for the model are Ho is the amplitude of H and is equivalent to a heating rate of 31.25' per hour for a short period.Such a rate of heating is common in small scale anelastic models (Hsiao-ming Hsu, personal communication).H I is chosen so that the heating is in a fairly narrow band and so that the total integral of heating is zero.H2 is selected so that the majority of the heating is in the lower part of R tailing off smoothly to zero at the bottom and top.Finally H, is such that the heating starts off smoothly from 0 at time r = 0 and peaks at 3 h.This allows us to start the model at rest, i.e., all variables can initially be set to zero.

CP To
Contour values for these four figures are given in units of lO-'s-I.Note that the advective term (Fig. 6a) is two orders of magnitude smaller than either the lower order term (Fig. 6b) or the heating term (Fig. 6c), while the sum of the lower order term and the heating term (Fig. 6d) is the same size as the advective term, Of course this means that we can replace w in (9.3d) by -(kP To)-' H. Then in the two dimensional case we really only need to compute a one dimensional elliptic equation for 4.
Similar savings will occur in three dimensions.

Tellus
equal to zero.Using standard L 2 energy estimates for the initial value problem we obtain in every finite time interval 0 < t < 7 the estimate a (5.8a) estimate for all space and time derivatives.From (5.7) this gives us for the unsymmetrized variables s,, u,, and u I the estimate IIS,II + lIUIlI + I l u 1 l l g 0 ( ~-" ) and a similar estimate for all of their derivatives.Using the differential equation (5.6a), we can then obtain for w, the estimate I I ~, I I G 0(&3-n) + s, s, @w + po p;' [yd + s, j ( z ) wl = 0, dt and d are as in (3.13).
S , = D -' T W = l .Since T < 100 the motion must satisfy the constraint D Q 10, W. If we choose D = Do, then W would have to have a value of 100 m s-,.Thus the motions must have a depth scale D such that D < Do.From (3.4) and (4.2) we have that S, = lo-' and we can use (5.2)

Fig. 1 .
Fig. 1.The eastern United States observed (a) zonal and (b) meridional components of the wind at 9 km and OO:oO GMT for January 13, 1979.Contour values are given in units of m s-'.

Fig. 2 .
Fig. 2. The eastern United States observed (a) zonal and (b) meridional components of the wind at 9 km and 18:OO GMT for January 13, 1979.Contour values are given in units of m s-I.

Fig. 3 .
Fig. 3.The eastern United States (a) zonal and (b) meridional components of the wind at 9 km and 18:OO GMT for January 13, 1979 computed by the model based on the forced barotropic vorticity system (9.1).Contour values are given in units of m s-l.
Richardson's equation.To obtain boundary data at intermediate times, we applied quadratic interpolation in time.Figs. 1 and 2 are the resultant (a) zonal and (b) meridional components of the wind at 9 Km and OO:OO GMT and 18:OO GMT, respectively.The contour intervals are given in units of m s-l.As a dynamic test of the data, we ran a forecast starting at OO:OO GMT using a forced barotropic vorticity model.The model was based on the system of equations c, + r q c o s e)-1 u ~A + vc81 + we, + BC + rl[(COS e)-l wA t,, -w8u,i + (cos e r " p y @$pA-pAp8)+ f B + r -' f $ v = ~, (9.la) (9.lb) 6 = (r cos e)-l bA + (COS 8 v)$i, where C = (r cos @-I [-(cos 8 u ) ~ + 2111.The bar notation indicates that the corresponding quantity was obtained from the data.
Fig.3shows the (a) zonal and (b) meridional components of the wind computed by the forced barotropic model at the same height and time as in Fig.2.The contour intervals are given in units of m s-I.As can be seen from the figures, the forecast is not perfect, but quite good.As a Tellus 38A (1986), 4 quantitative measure of the accuracy of the forecast we used the relative I, errors of u and v which had the values e(u) = ~l i i l l ~l l l i i -1411, = 0.12, e(u) = l l f i l l ~l l l f ivll, = 0.28.
variable u -(yp, pi1)-l12 # by the stable approximation mentioned above.At outflow points on the western boundary we assumed that u + ( y p , ~; l ) -" ~ 4 was known and updated the outflow variables s, u -( y p , p;')-lf2 #, u, and w by the same stable approximation.The remaining boundaries were treated analogously.For a we chose the value a = [(rA6)-l &I2.This ensures that all the terms in the elliptic equation for pressure derived in Section 5 are the same size as required for smoothness of the solution up to the boundary.The highest frequencies present in the model are given by v= (rcos OM)-' II + (rA@-l u + (&)-l w f (ypo/po)1'2 [ ( r cos BAA)-' + (rA6)-2 + ~( A Z ) -~] ~/ ~.The time step At must satisfy the CFL constraint max IvlAtG 1.To obtain an approximate value for At we neglected the term arising from the vertical advection and assumed a westerly wind with a speed of 100 m s-'.In this case v becomes v = (rM)-l 1224 + (6ypo/po)1'21 which gives an upper bound on the time step of approximately 225 s.We chose a time step of 180 s.The relative I, errors at 9 Km and 18:OO GMT for the model based on the adiabatic modified system were e( u) = 0.1 5, e(u) = 0.28.

Fig. 4 .
Fig. 4. The eastern United States (a) zonal and (b) meridional components of the wind at 9 km and 18:OO GMT for January 13, 1979 computed by the diabatic model based on the modified system (9.2).Contour values are given in units of rn s-'.

Fig. 2 .
Fig. 2. The contour intervals are given in units of m s-l.The relative f2 errors were e(u) = 0.13, e(v) = 0.28, which are essentially the same as the control run.We also developed a standard two-dimensional anelastic model to show that the balance between the vertical velocity and heating indicated in Section 8 is present.The two dimensional anelastic version of (7.5) is given by

Fig. 5 .
Fig. 5.The (a) horizontal and (b) vertical velocity components at three hours computed by the two dimensional anelastic model based on system (9.3).Contour values are given in units of m s-'.
Figs. 5a and 5b are contour plots of the horizontal and vertical velocity components c and w. respectively, at three hours.Contour values are given in units of m s-'.The typical signature of a thunderstorm can be seen in these figures.Figs.6a-d are contour plots at three hours of the following terms of the entropy equation: