APPENDIX A
 
VERTICAL w- IMPLICIT SOLVER

 

The pressure and vertical velocity are coupled in the vertical through a Crank-Nicholson scheme. This approach removes the vertical sound constraint on the small time step. The method follows that of Durran and Klemp (1983) and the ARPS (1995). The discretized equations for  and w are (from Chapter 3):

 

                                            

                                                                          (A.1)

                                            

                                                                           (A.2)

 

The time weighting term  is used to reduce the small to large time step instability and is applied to w in the  equation and to  in the w equation. A value of 0.6 is sufficient following the results of Durran and Klemp (1983). All the large time step forcing terms and the known small time step terms in (A.1) and (A.2) are grouped into wforce and pforce.

                                             (A.3)

                                            

                                                                 (A.4)

In (A.4) the total non-dimensional pressure at the current small time step multiplies the divergence term and the heating term. These terms are explicit and remove some of the gain obtained from the implicit application. The effect on the time step criterion is small since perturbation pressures are rarely greater than 3kPa or 1/30th of the base state pressure. Equations (A.1) and (A.2) are rewritten as:

 

                                                                 (A.5)

                                                                 (A.6)

 

Eliminating the pressure produces an equation for w at the future time step in terms of known terms.

                                                                

A(k), B(k), C(k), and D(k) are given by:

 

                                                                

                                                                

                                                                

                                                                

 

where,                                                     

                                                                

                                                                

                                                                

                                                                

The tridiagonal system is solved with appropriate boundary conditions. The solver is applied from 3 to nz-2. The top boundary condition for the rigid lid case is w=0. For the linear radiation condition between pressure and w the formulation from Durran and Klemp (1983) is used. The details of the upper radiation condition implementation are given in Appendix B. The performance of the implicit solver is roughly 2.0 times slower than the explicit version for a given small step. The implicit method is computationally effective when the ratio of the horizontal to vertical grid spacing is:

                                                                

 

 

APPENDIX B
 
LINEAR HYDROSTATIC w- TOP BOUNDARY CONDITION

 

Vertical energy propagation from hydrostatically forced gravity waves can be significant. A damping layer in combination with a rigid lid is commonly used in mesoscale numerical models to prevent reflection from a rigid lid. The rigid lid-sponge combination is effective (Klemp and Lilly, 1978) in preventing significant reflection. But for strong mountain waves the thickness of a properly designed damping layer can be as much as 1/2 of the model vertical domain. Another method applied by Klemp and Durran (1983) as well as others applies a linearized hydrostatic analytical relation between the vertical velocity and pressure variables at the top boundary. The advantage to this method is the number of grid points in the vertical is significantly reduced without degradation of the numerical solution. Fewer vertical grid points correspond to smaller memory and CPU cycle requirements. Since the relation between vertical velocity and pressure is linearized, the application is somewhat limited. Tests of the condition for a finite amplitude mountain wave reveals solution sensitivity to the upper radiation condition. The condition is more applicable to cases in which wave breaking takes place below the upper boundary.

The formulation implemented in ARPI3D follows Klemp and Durran (1983). For a linearized hydrostatic Boussinesq set of 3-dimensional equations, to which the derivation is not reproduced here, the Fourier transformed analytical relation between w and  is

                                                                  (B.1)

The relation (B.1) is applied at k=nz-2. The pressure and velocity points are staggered in the vertical, so (B.1) is approximate. An average of w at nz-1 and nz-2 proved to be unstable and the above relation is adopted. The outline of the solution method is:

    1. Apply the reduction step of the tridiagonal solver from Appendix A to w.
2. Starting with the pressure equation, obtain a relation for w(nz-2) in terms of

w(nz-1) and (nz-2). Substitute the result into the reduced step (1). Step

(1) becomes an equation relating w(nz-1) and  (nz-2).

3. Transform (2) into Fourier space.

    1. Apply (B.1) to (3) obtaining an equation for w(nz-1) in Fourier space it

      2. terms of known quantities.

    2. Solve (4) for w(nz-1), perform a reverse transform and apply the result to
the back-substitution phase of the tridiagonal solver to recover w at all

other levels.

 

 The above method is applied to both the vertically explicit and implicit w-solving techniques in ARPI3D. Step (2) uses the discretized pressure equation:

                                    

                                                                                                   (B.2)

The variable pforce includes all the known terms in the pressure equation including the big time step advection and heating terms and the small time step updated horizontal velocity divergence and a weighted portion (1-) of the vertical velocity divergence from the previous small time step. The reduction phase of the tridiagonal solver provides a relation between w(nz-2) and w(nz-1).

 

                                                                                     (B.3)

 

The coefficients c and d have been modified from the reduction steps. Substituting (B.3) into (B.2) gives:

 

                                                            (B.4)

 

Where A1 and A2 are the coefficients of w(nz-1) and w(nz-2) in (B.2), respectively. The last step is the transformation of (B.4) into Fourier space and the substitution of (B.1) into (B.4) to obtain an equation for w(nz-1) in terms of known quantities. The vertical velocity at k=nz-1 is then transformed back to real space and the remaining vertical velocities at all other levels are obtained using the back substitution phase of the tridiagonal solver. Pressure is then computed using the updated three-dimensional velocity field.

 

 

APPENDIX C
STREAMLINE METHOD AND TESTS

Two-dimensional streamline-trajectory computations used in Chapter 4 to validate the numerical model formulation are described and tested here. The trajectory model applies the predictor-corrector method to the three-dimensional velocity field. No generality is lost in presenting only the x-direction, as this analysis can be equally applied in the y and z directions. In the x-direction, a parcel trajectory is defined as:  (C.1)

where (A.1) can be integrated to obtain the new parcel location given the current parcel velocity and position. A simple approach to computing the new parcel location would be to estimate the velocity at the present location and apply a time increment obtaining the distance traveled and the new position. Assuming that the velocity field changes slowly in space the deviation from the true parcel movement would be small. However, for the case of strong spatial velocity gradients large errors are possible from simply estimating the initial velocity as the average velocity over the parcel displacement. To improve the accuracy of the scheme, the predicted parcel movement is used to obtain an updated or corrected velocity field. The initial parcel position  is used to obtain a linearly interpolated velocity . The time step is applied to the interpolated velocity field and an updated intermediate position  is obtained. The relationship describing the parcel displacement is

 

(C.2)
 where  is the current interpolated velocity field. The corrected velocity field is determined using:
(C.3)
 With the updated velocity, Equation (C.2) is reapplied n times producing a final parcel location  for each time step. For the n=1 case, there is no velocity correction and only Equation (C.2) is used. This process can be repeated n times but little is gained after the 2nd application. For a slowly varying velocity field, n=2 produces a convergent streamline pattern. This method is similar to the Adams-Moulton presented by Hornbeck, (1975). Tests of this trajectory method were conducted for a radially varying velocity field. The error is a function of the time step (dt) and the number of iterations (n). Figure C.1 displays the streamlines using a radially varying two-dimensional wind field for two different number of corrector steps. The radial velocity fields are defined by:

                                                                                 and 

where L is the width and H is the height of the test domain, xcenter and zcenter are the central locations, and  and  are the velocity magnitudes.

 

 

 

(a)

 

 

 

 

 

 

(b)

 

 

 

 

 

 

 

Figure C.1. Streamline test pattern for =20 m/s, =20m/s, dt=25 seconds, dx=400m, dz=250 m , for (a) n=1, and (b) n=2. The CFL criteria is approximately 0.2.