1. Introduction

The general concept used to consider optimization of the computational system of atmospheric prediction equations is to find, if possible, eigensolutions of the most sophisticated linearized form of the nonlinear equations. If these solutions henceforth to be denoted as 'structures' can be used in a spectral expansion, then the numerical integrations will be carried out using the best representation available. If a finite grid of points is needed, then a grid is sought on which the structures to the linear system are also eigensolutions. Since these structures represent variability in the space domain, they are scale dependent. It is not always practical or even possible to use structures for the most complex system, thus we search for structures which satisfy at least some linear form of the equations. The structures can also be selected from observations provided they have significant amplitude in the largest scales or they might be selected from EOF analysis of data.

As examples of this process, the earliest models were two dimensional and expanded in Fourier series on an equally spaced grid; the structures satisfied the barotropic vorticity equation (BVE) in linear form since the finite difference grid yields the same spectral expansion functions. Expansion of the flow field on these structures showed rapid decay with scale and when applied in an integration, the short term forecasts were reasonable. Subsequently the models were moved onto the sphere and the appropriate structure functions were recognized as Legendre polynomials which satisfied the BVE. In this format scaling analysis suggested a two dimensional index for the representation and triangular (T) truncation as optimal. Projection of data onto these structures showed rapid decrease of amplitude with decreasing scale and truncation experiments identified sensitivity of closure to scale.

With the evolution of baroclinic models, which implies the coupling of horizontal and vertical structures and their scales, most models use either finite difference or spectral (Legendre polynomial) representation. All models use a form of finite differencing in the vertical with no apparent scale consideration; the choice of number and distribution of vertical points has remained at the discretion of the modeler. Some studies of vertical structures have been made and have suggested coupling clues.

Our approach to studying the optimization problem is to first select an appropriate model and convert it to computational form either in physical space via a finite difference grid, a spectral expansion, or a mix of the two. Although we will transfer our results to the nonlinear system, we do the analysis on the linearized version of the model. Truncation is a critical feature of the analysis because the nonlinear interactions are scale dependent and energy amplitude shows a significant decrease with scale. For any of the representations, truncation occurs at the short scale end of the spectrum. The characteristic structures which satisfy the exact or computational linear system should be independent and thus their properties provide the optimization information sought.

2. Characteristic structures

As noted above, the characteristic structures are selected to represent the dependent variables of the linear system and in that framework can be written most generally in the following form:

(1),

where represents the variable and F_n,m,k is a three dimensional characteristic structure. Note that the indices (n,m,k) can be considered scaling indices which have unique properties defined by the particular structures chosen. The total variance in can be determined from a truncated subset of the defined structures. Since the indices are scale dependent, truncation can be made scale selective and knowing the scale dependence of the amplitudes of the structures, nonlinear exchanges can be minimized by optimal truncation.

The most general linear system representing the atmospheric prediction equations would yield only a highly confusing set of coupled structure functions in three dimensions. Thus we identify those structures which will represent simpler versions of the system. As an example if one considers the BVE on the sphere, only horizontal structures ensue and they are seen to be products of Fourier series in longitude and Legendre polynomials in latitude; i.e., the scales in the separate dimensions uncouple. The scaling properties of this system have been studied in some detail ( see Baer, 1972 for example). For the shallow water equations (SWE) the structures also uncouple, but the structures in latitude become Hough functions, and a vertical structure equation evolves which must be solved for the vertical structures. The scaling here is somewhat more complicated and depends on the vertical structures selected.

A simple view of three dimensional scaling may be had by considering quasi-geostrophic flow on the sphere. In this case, if s=s(n,m,k), and the differential operator L is applied to Equation (1) , L(F_s) = -a²s²F_s and the eigenvalues are defined as s²(n,k) = n(n+1) + a²f_o²/gh_k . Note that the three dimensional index (s) depends only on the latitudinal index (n) and the vertical index (k), where h_k is the kth equivalent depth determined from the vertical structure equation and associated with a vertical structure function as yet to be determined.

a. Horizontal structures

As a rule, horizontal structures used in models over the sphere are the associated Legendre polynomials. They have simple and well defined properties and they are scaled using a two dimensional index which is the ordinal number of the functions. For more detail, see Baer (1972). These functions do not solve a very complex differential equation, but show realistic distribution of energy with scale when applied to atmospheric data, as comprehensively summarized by Wiin-Nielsen and Chen (1993). We shall discuss this issue somewhat further subsequently.

b. Vertical structures

The vertical structure problem has been addressed by a number of researchers as for example Lindzen and Fox-Rabinowitz (1989) and Baer and Ji (1989, BJ). The method of BJ to which we will refer here found the solutions to the vertical structure equation (see also Kasahara and Puri, 1981) using very high resolution. An optimal grid with K points was then established on which the first K structures to the high resolution solution were also the solution to the difference equations of the problem on the optimal grid. These eigenstructures together with their equivalent depths were used with the Legendre polynomials to define the three dimensional characteristic structures, F_m,n,k.

3. Optimal structures, gridding and validation

Our assumption is that the best predictions would come about if the model equations were expanded in the "best" characteristic structures and the set of structures was truncated with a three dimensional index such as s defined above. For model applications, it is traditional to expand spectrally in the horizontal surface using the associated Legendre polynomials as the horizontal structure functions, but the vertical domain is generally gridded. Our approach here is to use the vertical eigenstructures with the optimal grid defined above. Thus we have a mixed spectral-gridded system, but we have structures in all dimensions and we have a three dimensional scaling index. Combinations of n and k giving curves of constant s are seen on Fig. 1 and represent the basis for 3-D scale truncation. For any desired truncation in s, simply select the n and k indices from the intercepts of that s curve. Although this gives somewhat too large a truncation, it approaches the optimum. If the expansion were entirely spectral, one could stay on the s-curve for all n and k.

It remains to be established whether this optimization concept actually works. Thus a number of experiments were undertaken to test the feasibility of the technique.

FIG. 1: Three dimensional index s as a function of the horizontal (n) and vertical (k) indices as determined from quasi-geostrophic theory ( see text).

a. Tests with an 'optimal' vertical grid

The first test was made with the NCAR/CCM0B and was designed simply to check if moving the vertical levels to their optimum positions would impact on the prediction. No connection was made between the horizontal and vertical truncation. The model was run in its standard form with R15, L9 configuration with only the level positions altered. Comparison of twin experiments using the standard levels and test (optimal) levels showed distinct improvement in the forecast of the higher vertical modes out to 4 days with the optimally gridded model. These results were reported in BJ.

This success led to more detailed testing of the optimum grid, detailed by Baer and Zhu (1992). In that study, the NCAR/CCM1 was used with T31, L12. Again the integrations for standard and optimal levels were compared but in this case a number of different initial conditions were used so that statistics of the results wiuld be more robust; additionally comparison with observations was possible. The results showed clearly that the component of the temperature in the external vertical mode and the component of the vorticity in the internal vertical modes were predicted more accurately using the optimal grid. This implies that the optimal grid is representing the vertical structure of the system with more fidelity.

b. Tests of three-dimensional optimal scaling

That selecting a vertical grid which is more consistent with relevant vertical characteristic structures shows prediction improvement is gratifying. However it does not give insight into the coupling of horizontal and vertical structures nor if there is a systematic way to truncate in all dimensions simultaneously. We suggested such an approach above using the quasi-geostrophic system, by implying that one might truncate using the s index. We tested this hypothesis using the optimal vertical levels and a spectral expansion in the horizontal applied to an atmospheric model based on an analysis scheme by Ji and Baer (1992) which defined an error function applicable to any model integration results. If the selection of (N,K) truncation for a chosen model was optimum (assuming triangular truncation), that error function would be a minimum when compared to a truncation of the model with (N',K), N'_N. Clearly this is relevant for any K and indicates how the truncation of n"N and k"K are coupled.

We tested this hypothesis using the NCAR/CCM0B with dynamics only. Six horizontal resolutions (T15, T21, T26, T31, T35, T42) were tested for L9 using optimally chosen levels. Ten integrations with different initial conditions for each truncation were run out to 8 days. The ensemble average of the ten integrations showed a distinct minimum at N~27, near the value one finds from Fig.1. Additional tests were subsequently made with the NCAR/CCM1 using 12 levels. Horizontal truncations of T21, T31, T42, T63, T95 were tested and the minimum appeared at N~42. Although slightly larger than the value suggested by Fig.1, it is quite reasonable.

Our final test was with the NMC/MRF Research Test Version (1990) model. If the technique has value to forecasting, we felt it should be tested on a true forecast model. The model had 18 levels including four in the boundary layer and had a top at 20 hpa. We reset the levels above the boundary to the optimal levels as defined herein. Three horizontal levels were tested (T40, T62, T80) and six 10-day integrations were performed for each resolution. Statistics of the integration results confirmed what we had seen from the CCM results. However, the errors tended to plateau after the T62 resolution; i.e., the T80 results were similar to the T62 results but both showed improvement over the T40 integrations. The implication of this observation is that T62 or thereabouts is perhaps the optimum resolution although using higher resolution, aside from wasting computer time, does not degrade the forecast. Interestingly, Fig.1 ( based on the equation for s) indicates an optimum value for L=18 to be N~58.

c. Additional considerations of horizontal structures

Consideration of the shallow water equations has already indicated that the Legendre polynomials are not the characteristic structures to the linearized horizontal part of those equations but that Hough functions are. Although using the Hough functions as expansion functions in a spectral model is not feasible, a grid in latitude on which a subset of the Hough functions are orthogonal in a finite difference sense (using a process similar to the one used in BJ for the vertical) might optimize a model's predictive skill. Using a finite difference form of the horizontal shallow water equations and applying both a fourth-order and a finite element scheme on a uniform grid, we demonstrated that in addition to the Hough modes, other modes evolved which did not satisfy the original differential equations. These modes--we denote them as false modes--must somehow be removed to optimize the integrations and can possibly be eliminated by a judicious choice of grid intervals similar in concept to the technique utilized in the vertical domain.

4. Conclusions and ongoing studies

We have suggested an approach to systematically truncating the atmospheric prediction equations based on scaling considerations which might optimize a numerical prediction since it implies a coupling of all the space dimensions, thereby minimizing aliasing errors. In addition, we have demonstrated by experiments with several models that the process has an effect and that effect seems to yield the expected improvements. Yet the development to date is somewhat too primitive since it is based on linearization with no mean flow. It is evident that the atmosphere has both a systematic mean flow and shear. If such currents are introduced into the linearized model, the vertical and horizontal equations do not separate and a linear system in two dimensions must be considered.

Such a system is currently under study. We have isolated characteristic structures in both latitude and height and are identifying their properties. We will select a subset of these structures which are as close to the solutions of the differential system as possible and use them to test a numerical model either in an expansion mode or on a suitably optimized grid. We anticipate that results from these experiments will yield the optimum truncation.

Baer, F., 1972: An alternate scale representation of atmospheric energy spectra. J. Atmos. Sci., 29, 649-664.

---, and M. Ji, 1989: Optimal vertical discretization for atmospheric models. Mon. Wea. Rev., 117, 391-406.

---, and Y. Zhu, 1992. Forecast accuracy with optimum vertical model truncation. Mon. Wea. Rev., 120, 2579-2591.

Ji, M., and F. Baer, 1992: Three-dimensional scaling and consistent truncation of global atmospheric models. Mon. Wea. Rev., 120, 131-148.

Kasahara, A., and K. Puri, 1981: Spectral representation of three-dimensional global data by expansion in normal mode functions. Mon. Wea. Rev., 109, 37-51.

Lindzen, R., and M. Fox-Robinovitz, 1989: Consistent vertical and horizontal resolution. Mon. Wea. Rev., 117, 2575-2583.

Wiin-Nielsen, A., and T-C. Chen, 1993: Fundamentals of Atmospheric Energetics.Oxford, 376pp.