Pressure Gradient in a Sigma Coordinate

A method is proposed to calculate the pressure gradient force in a sigma­ p coordinate over mountains. The method uses a local reference state to calculate the pressure perturbation, which is then used to calculate the pres­ sure gradient force. In addition to mathematical formulas, the author also presents the numerical results, showing the proposed method in�eed produces very accurate results. Another advantage to the method is that there is no need to calculat� any field beneath the ground.

It is well known that the errors involved in the pressure gradient term in a sigma coordinate can be quite large over mountain areas. Several different schemes have been developed in order to reduce this error (Kurihara, 1968;Corby et al., 1972;Phillips, 1974;Mahrer, 1984;Mihailovic and Janjic, 1986;Mesinger et al., 1988;Kuo et al., 1988;Juang and Ogura, 1990;Simmons and Chen, 1991;Sun et al., 1991;Danard et al., 1993). One popular use of such methods has been to define a horizontally-invariant reference state, (i.e., a universal reference state), such that the pressure gradient can be divided into two terms, naming the reference state and the perturbation, respectively. Since the reference state is horizontally invariant, (although it may change with time), the pressure gradient force is also zero. One needs only to evaluate the pressure gradient force of the perturbation, which usually has a smaller error. On the other hand, Kurihara (1968) and Mihailovic and Janjic (1986) used interpolation of the geopotential to calculate the g�opotential difference at constant pressure levels. Similarly, Danard et al. (1993) proposed a numerical scheme to calculate the pressure gradient in a sigma-coordinate in a barotropic atmosphere. Instead of reviewing the details of the existing methods, this paper focuses on tpe comparison . s of the pressure gradient calculated based upon the universal reference state and the local reference state presented here.
If the atmosphere is not uniform horizontally or a universal referenced state cannot be chosen properly, however, (which is often the case in a real atmosphere), a large error can still be introduced. In order to overcome this problem, the author proposes c�lculating pressure perturbation based upon the local reference state. This pressure _perturbation is then applied to calculate the pressure gradient equation in a sigma-coordinate. This method may be considered a combination of the Kurihara approach by which the pressure difference at a constant height between two vertical columns is calculated and the approach of subtraction from a reference pressure (e.g. Phillips, 1974 and others).
The results obtained from an idealized atmosphere and a real atmosphere simulation show that the error can be reduced significantly. Meanwhile, the interpretations and calcula tions are very straightforward, although more computation time is required when compared , with the methods using the universal reference state. This method has also been applied successfully in the Purdue Mesoscale Model Chem, 1993, 1994) to study the flow over complex terrain over Taiwan Chem, 1993, 1994) as well as_ the Rocky Mountains (Chem, 1994) .

THE BASIC EQUATIONS AND NUMERICAL METHOD
.

The hydrostatic equation for a dry air is
. .
In the Arakawa C -grids shown in Figure 1, the temperature is calculated at the levels of Tk l / 2 and Tk+ 1 ; 2 ; the pressure (p) and geopotential (</>) is defined at the levels of P k -1 , P k, and P k+ l · We can integrate (2.1) and obtain the geopotential: where . .
(2.5) which is different from (2.2a). Based on experience, the· author fi nds that the geopotential calculated from (2.3) is more accurate than from (2.2a) because (2.3) uses a linear piece-wise function while (2.2a) makes use of a step function. .
The schematic diagram of Arakawa C-grid: the ground is defined at the km-level.

The Horizontal Pressure Gradient
The x-component momentum equation is The advantage of using perturbation pressure p1 instead of p has briefly been discussed in the introduction. It is important to defi ne a reference state such that p' can be calculated appropriately. In an idealized case with a constant lapse ((3 k = (3 = constant), the reference pressure Pr can be defi ned as (2. 7), • where z 0 is the height at surface. The pressure perturbation p' becomes The pressure gradient in a sig�a-coordinate can then be easily obtained in accordance with the right.hand side of (2.6).
This method fails if the temperature is not uniform horizontally or if {3 is a function of height. Hence, a different method is proposed for calculating the reference pressure p k and pressure perturbation p'. · For any arbitrary terrain shown in Figure 1, the pressure gradient at u ( i-1 / 2, k+ 1/2) can be calculated as follows: The neighboring column (i.e., (i -1) th or (i) th column) with a lower terrain will be used as the local reference state. In this case, the (i) th column is chosen as the reference column. The pressure perturbation p' ( i -1, k ) is defined as (2.9) The value of p ( xi-l , z k ) is the pressure at the ( i -1, k ) grid. The pressure at the same height, but at the i-th column, p ( xi, z k ) can be calculated: Then p' ( x i-1 , Z k ) can be calculated from (2.9). It is also noted that p' (xi, Zm) = 0 in the entire column, because the (i) th column is chosen as the reference state. At this time, (2.6) can be calculated.
It is important that the calculation of the pressure gradient for u at the ( i -1/2) th column be completed before moving to the next column.
When we calculate the pressure gradient for u at the ( i + 1/2) th column, the. ( i + 1) th column is used as the local reference state because the terrain at (i + 1) is lower than that at ( i), as shown in Figure 1. Consequently, the calculation can be performed on a column by-column basis along the x-direction, then along the y-direction. The method proposed here does not require a universal reference state for an entire domain, as proposed by others.

NUMERICAL RESULTS FOR IDEALIZED ATMOSPHERE
We defi ne a domain of (30 �x, 12 km), consisting of 30 x 31 or 30 x 11 grids, where the horizontal increment �x = 1 km. Following Juan and Ogura (1990), we impose a horizontally-uniform atmosphere with a constant lapse rate, {3 = -6.3 x 10-3 K m -l. The sea surface pressure is 1013.25 mb, and the sea surface temperature is 288 K. A triangle mountain spaning over 6 �x with a peak of 3 km is also imposed at the center of the domain. Then, a th column and the loGation of the same height at the ( i) th column.
sigma coordinate, a = (p(x,z)-pt)/(psfc(x)-pt) is introduced, where psfc is the pressure at the surface, and pt is the pressure at 12 km. It is divided into 30 (or 10) layers with a unifor1n �a = 1/30 (or 1/10). With the total absence of wind, the perfect solution of -( 1/ p )(op' I ox) z should be zero everywhere in this case. The procedure in a numerical simulation is that we integrate the prognostic equations in order to calculate the values of T ( i, k ± 1 /2) and the surface pressure, psfc(i), for each time step. In this way, the pressure at the ( i, k) grid can be obtained: p(i , k ) p t+ a ( k ) * ( p sfc (i) -p t ).
The geopotential at p( i, k ) is then calculated from the hydrostatic equation (2.3), which is shown in Figure 3. Finally, the pressure gradient force, the Coriolis force, and the new velocity fields are calculated. .

Case A With 30 Layers Vertically
First, we introduce a (universal) reference �tate according to (2.7) with z0 = 0, T 0 = 288 K, p0 = 1013.25 mb and (3 = -6.3x io-3 K m-1 . Then, we calculate p' from (2.8) and -(1/p)(op ' /ox)z from (2.6). In order to avoid errors which can be generated by the computer itself, we use long-double-precision (Implicit Real* 16) at the IBM RISC-6000 workstation. The results, shown in the (x, f ) coordinate (Figure 4), reveal that the error of -(llp)(op'lox)z is almost zero ( 3.2x 10-2 m s-2 ) in the entire domain because we are able to find a perfect universal reference state. This method has also been successfully applied to study the formation of lee-vortices in Taiwan in an idealized atmosphere (Sun et al., 1991).

Case B With 30 Layers Vertically
The procedure in Case A is repeated, except that the local reference state defined in (2.10) is used, and (2.9) is employed to calculate p'. The error of 3.0x10-2 9 m s-2 in the results ( Figure 5) is comparable to 3.2x 10-2 9 m s-2 in Case A. This indicates that both methods are perfect for this particular case in that no error has been generated by either method in the entire domain.

Case C With 30 Layers Vertically
The same procedure as in Case A is duplicated, except that a different universal reference state with f3 = -6.0x 10-3 K m-1 in (2.7) is used. The error of -(1/p)(op' /ox)z reaches 2.5x 10-4 m s-2 ( Figure 6), which is.much larger than that in any other previous cases. This error generates a positive pressure gradient force over the eastern slope of the mountain and a negative one over the western slope. When /3= -3.0x 10-3 K m-1 is used in (2.7) as the reference state, the error of -(1/ p )(op' /ox)z reaches 2.25 x 10-3 m s-2 (not shown here).
Hence, with an inappropriate reference chosen to calculate the pressure perturbation, the error can be huge. Errors always exist as long as imperfections are present in the reference. It is worth noting that it is almost impossible to fi nd a perfect universal reference in a real atmosphere, but a perfect local reference state can always be determined according to the method proposed in this paper. In short, therefore, the error is very small when the local reference state is used.

Case D With 10 Layers Vertically
It is well know that the error of the pressure gradient force may also depend upon the vertical resolution. Figure 7 show the error reaches 1.5x 10-3 m s-2 using /3= -6.0x 10-3 K m -l as a universal reference. On the other hand,as shown in Figure 8 using a local reference, the error remains about 1.6x 10-2 9 m s-2 , which is much smaller than the results of Juang and Ogura. Furthermore, the method proposed here does not require that the pressure under the ground be calculated, as proposed by Juang and Ogura.

NUMERICAL RESULTS IN A REAL ATMOSPHERE
A simplified version of the Purdue Mesoscale Model (PMM), without radiation, con densation, and vertical mixing, is applied to simulate the cold front passing through Taiwan during the Taiwan Area Mesoscale Experiment in 1987. The ECMWF analysis of OOOOZ June 14 is used as the initial condition. The domain includes 106 x 96 grids with a horizontal space interval of 15 km (Figure 9). In the vertical direction, there are 25 layers extending from --.
the surface up to 100 mb. The detailed description of the PMM can be found in Sun et al., (1991) and Chem (1993, 1994). Figures lOa and b show the vertical cross-section of the temperature fi eld along the line AA' in Figure 9 after three hours' integration. It is noted that during the integration the horizontal wind is very light over Taiwan ( 2.5 mis) . The temperature perturbations over mountains are mainly caused by persistent vertical motion. The results (Figures lOa-b) show that over mountains in Taiwan and China, shortwaves are much stronger using a universal reference state (with f3 = -5.5 x 10-3 K m -l � an average value of the entire domain) than those using a local reference. The vertical motions (not shown here) are also stronger over Taiwan using a universal reference state, as expected. The value of -(llp)(8p'l8x)z obtained using the local Contours are 0.25x 10-2 9 from -3.0x 10-2 9 to 3.0x 10-2 9 . ..
.. Contours are from -1.60x 10-29 to 1.6x 10-29 m s-2 with interval of 0.2x 10-29 m s-2 (10 layers).  I  I  I  I  I  I  I  I  I  I  I  I   o: t   I  I  I  I   I  I  I  I  I  I  I  I  I  I  The PMM with full model physics has also been applied to integrate 30 days over the USA to study the 1988 drought (Sun et al.,t · 995). The simulations are in good agreement with the observations, which may indicate that the error does not grow with time in this proposed method. However, we believe that more cases should be tested.

CONCLUSION
The error of the pressure gradient over mountains can still be quite large even when a universal reference state is applied. This study shows that this error can be reduced drastically by using the local reference state to calculate the pressure perturbation. Although more computing time is required than when a universal reference state is used, the method proposed here is very straightforward and can defi nitely produce more accurate results. Furthermore, this method can also be applied in the sigma-z coordinate or other sigma-coordinates in nonhydrostatic systems. • 590 TAO, Vol.6, No.4, December 1995