An Introduction to Finite Differencing

This chapter on finite differencing appears oddly placed in the early part of a text on spectral modeling. Finite differences are still traditionally used for vertical differencing and for time differencing. Therefore, we feel that an introduction to finite-differencing methods is quite useful. Furthermore, the student reading this chapter has the opportunity to compare these methods with the spectral method which will be developed in later chapters. One may use Taylor’s expansion of a given function about a single point to approximate the derivative(s) at that point. Derivatives in the equation involving a function are replaced by finite difference approximations. The values of the function are known at discrete points in both space and time. The resulting equation is then solved algebraically with appropriate restrictions. Suppose u is a function of x possessing derivatives of all orders in the interval (x — n∆x, x + n∆x). Then we can obtain the values of u at points x ± n∆ x, where n is any integer, in terms of the value of the function and its derivatives at point x, that is, u(x) and its higher derivatives.

Mathematical Aspects of Spectral Models

In this chapter we provide an introduction to the topic of spherical harmonics as basis functions for a global spectral model. The spherical harmonics are made up of trigonometric functions along the zonal direction and associated Legendre functions in the meridional direction. A number of properties of these functions need to be understood for the formulation of a spectral model. This chapter describes some useful properties that will be used to illustrate the procedure for the representation of data sets over a sphere with spherical harmonics as basis functions. The calculations of Fourier and Legendre transforms and their inverse transforms are an important part of global spectral modeling, and these are covered in some detail in this chapter. Finally, this chapter addresses the formulation of two simple spectral models. One of these is a single-level barotropic model, and the other is a shallow-water model.

Introduction

The spectral modeling approach to numerical weather prediction is being practiced in many parts of the world. Historically, the spectral approach came into atmospheric sciences from studies of geomagnetism, where it was introduced by Elsasser in the late nineteenth century. The first attempt of spectral representation of data sets and its use via simple vorticity conserving models came in the 1940s. The works of Neamtan (1946) from the University of Manitoba and Haurwitz (1940) and Craig (1945) from New York University were pioneering during this era. These were still close to linear problems, lacking any formalism for addressing the nonlinear advective dynamics. It was in the late 1950s when we saw the emergence of formal proposals for the solution of the nonlinear barotropic vorticity equation. Pioneering work from the University of Chicago (Platzman 1960; Baer 1964) explored what is sometimes called the interaction coefficients approach for the nonlinear problem. Around the same time at MIT we saw the elucidation of what are now called the low-order systems. This pioneering work of Lorenz (1960b) and Saltzman (1959) brought to us the first exposure to simple nonlinear systems and the concept of chaos. These simple three-component systems demonstrated some of the essentials of nonlinear dynamics and the growth of errors arising from initial state uncertainties. The interaction coefficients approach to the solution of the weather forecast problem led to unmanageably large memory requirements that were not easily amenable to the then available, or even to the present, memory of computers. It was during the mid-1950s when the Cooley-Tukey algorithm (1965) emerged and provided a breakthrough via the fast Fourier transform. This was exploited and demonstrated to provide accurate representations of the quadratic terms for fast computation of the nonlinear advective dynamics by contributions from Eliasen et al. (1970) and Machenhauer (1974) from the University of Copenhagen and from Orszag (1970) at MIT. Thereafter, we saw a rapid development of global spectral models in many parts of the world, especially Australia, Canada, England, Japan, and the United States.

Spectral Energetics

In this chapter we present spectral energetics. This is a useful tool for the interpretation of model output. It can be used to interpret both short term weather evolution and climate time scales. These same procedures can be used with both real atmospheric data and model output. A comparison of energetics histories can be very useful for the assessment of model performance. In the first section of this chapter we derive the equations for atmospheric energetics. In the following section, a method for the representation of these equations in the Fourier domain is introduced, and the equations are derived in the one-dimensional (zonal) wavenumber domain. In the last section, we view the problem in two-dimensional wavenumber domain and derive expressions for barotropic energy exchanges and baroclinic energy conversions in this framework. Some sample results of the energetics in two-dimensional wavenumber domain are also presented in this section. In this section we consider a system of basic equations in spherical coordinates and derive the relevant energy equations for the zonally averaged flow and the eddy flow. These derivations are essentially based on the work of Saltzman (1957).

Initialization Procedures

In this chapter we describe two of the most commonly used initialization procedures. These are the dynamic normal mode initialization and the physical initialization methods. Historically, initialization for primitive equation models started from a hierarchy of static initialization methods. These include balancing the mass and the wind fields using a linear or nonlinear balance equation (Charney 1955; Phillips 1960), variational techniques for such adjustments satisfying the constraints of the model equations (Sasaki 1958), and dynamic initialization involving forward and backward integration of the model over a number of cycles to suppress high frequency gravity oscillations before the start of the integration (Miyakoda and Moyer 1968; Nitta and Hovermale 1969; Temperton 1976). A description of these classical methods can be found in textbooks such as Haltiner and Williams (1980). Basically, these methods invoke a balanced relationship between the mass and motion fields. However, it was soon realized that significant departures from the balance laws do occur over the tropics and the upperlevel jet stream region. It was also noted that such departures can be functions of the heat sources and sinks and dynamic instabilities of the atmosphere. The procedure called nonlinear normal mode initialization with physics overcomes some of these difficulties. Physical initialization is a powerful method that permits the incorporation of realistic rainfall distribution in the model’s initial state. This is an elegant and successful initialization procedure based on selective damping of the normal modes of the atmosphere, where the high-frequency gravity modes are suppressed while the slow-moving Rossby modes are left untouched. Williamson (1976) used the normal modes of a shallow water model for initialization by setting the initial amplitudes of the high frequency gravity modes equal to zero. Machenhauer (1977) and Baer (1977) developed the procedure for nonlinear normal mode initialization (NMI), which takes into account the nonlinearities in the model equations. Kitade (1983) incorporated the effect of physical processes in this initialization procedure. We describe here the normal mode initialization procedure. Essentially following Kasahara and Puri (1981), we first derive the equations for vertical and horizontal modes of the linearized form of the model equations.

Physical Processes

In this chapter we present some of the physical processes that are used in numerical weather prediction modeling. Grid-point models, based on finite differences, and spectral models both generally treat the physical processes in the same manner. The vertical columns above the horizontal grid points (the transform grid for the spectral models) are the ones along which estimates of the effects of the physical processes are made. In this chapter we present a treatment of the planetary boundary layer, including a discussion on the surface similarity theory. Also covered is the cumulus parameterization problem in terms of the Kuo scheme and the Arakawa- Schubert sheme. Large-scale condensation and radiative transfer in clear and cloudy skies are the final topics reviewed. There are at least three types of fluxes that one deals with, namely momentum, sensible heat, and moisture. Furthermore, one needs to examine separately the land and ocean regions. In this section we present the socalled bulk aerodynamic methods as well as the similarity analysis approach for the estimation of the surface fluxes. The radiation code in a numerical weather prediction model is usually coupled to the calculation of the surface energy balance. This will be covered later in Section 8.5.6. This surface energy balance is usually carried out over land areas, where one balances the net radiation against the surface fluxes of heat and moisture for the determination of soil temperature. Over oceans, the sea-surface temperatures are prescribed where the surface energy balance is implicit. Thus it is quite apparent that what one does in the parameterization of the planetary boundary layer has to be integrated with the radiative parameterization in a consistent manner.

Multilevel Global Spectral Model

Since the 1970s, the spectral method has become an increasingly popular technique for global numerical weather prediction. Global numerical models formulated using the spectral technique are used worldwide for both research and operational purposes. The success of the spectral technique can be attributed to the spectral transform technique developed independently by Eliasen et al. (1970) and Orszag (1970), and later refined by Bourke (1972). Prior to the introduction of the transform technique, the nonlinear terms were computed using a very tedious process called the interaction coefficients method. This method required large amounts of computer resources as well as enormous bookkeeping. The transform technique facilitates the computation of the nonlinear terms, as discussed later in Section 7.4. Furthermore, the Galerkin method discussed in Chapter 4 is widely used in most spectral models and provides us with alias-free computation of the nonlinear terms. The transform technique enables the current spectral models to be competitive in terms of computational overhead with respect to their grid-point counterparts. The transform technique is also well-suited for incorporating the terms dealing with physics in the prediction scheme. There are a number of advantages to using the spectral technique over the conventional grid-point method. However, we will not get into this discussion here. It should be noted that the model truncation limit specifies the scale of the shortest wavelength that can be resolved by the model. In the following section, we discuss the two most widely used truncations in a spectral model.

What Is a Spectral Model?

If one takes a closed system of the basic meteorological equations and introduces within this system a finite expansion of the dependent variables using functions such as double Fourier or Fourier-Legendre functions in space, then the use of the orthogonality properties of these spatial functions enables one to obtain a set of coupled nonlinear ordinary differential equations for the coefficients of these functions. These coefficients are functions of time and the vertical coordinate, since the horizontal spatial dependence has been removed by taking a Fourier or a Fourier-Legendre transform of the equations. The coupled nonlinear ordinary differential equations for the coefficients are usually solved by simple time-differencing and vertical finite-differencing schemes. The mapping of the solution requires the multiplication of the coefficients with the spatial functions summed over a set of chosen finite spatial basis functions. This is what defines spectral modeling. Meteorological application of the spectral method was initiated by Silberman (1954), who studied the nondivergent barotropic vorticity equation in the spherical coordinate system using the spectral technique. In its earlier days, the spectral method was particularly suitable for low-resolution simple models. The equations of these simple models involved nonlinear terms evaluated at each time step. Evaluation of the nonlinear terms was performed using the interaction coefficient method and thus required large memory allocations, which was an undesirable proposition. However, with the introduction of the transform method, developed independently by Eliasen et al. (1970) and Orszag (1970), the method for evaluation of these nonlinear terms changed completely. This transform method also made it feasible to include nonadiabatic effects in the model equations. For the past couple of decades, the spectral method has be come an increasingly popular technique for studies of general circulation and numerical weather prediction at the operational and research centers. This method forms the basis for spectral modeling, and it is easy to understand if the reader has some background in linear algebra. We have a set of linearly independent functions θi(x), which are called the basis functions. The dependent variables of the problem are represented by a finite sum of these basis functions.

Time-Differencing Schemes

Atmospheric models generally require the solutions of partial differential equations. In spectral models, the governing partial differential equations reduce to a set of coupled ordinary nonlinear differential equations where the dependent variables contain derivatives with respect to time as well. To march forward in time in numerical weather prediction, one needs to use a time-differencing scheme. Although much sophistication has emerged for the spatial derivatives (i.e., second- and fourth-order differencing), the time derivative has remained constructed mostly around the first- and second-order accurate schemes. Higher-order schemes in time require the specification of more than a single initial state, which has been considered to be rather cumbersome. Therefore, following the current state of the art, we focus on the first- and second-order accurate schemes. However, higher-order schemes, especially for long-term integrations such as climate modeling, deserve examination. We start with the differential equation dF/dt = G, where F = F(t) and G = G(t). If the exact solution of the above equation can be expressed by trigonometric functions, then our problem would be to choose an appropriate time step in order to obtain a solution which behaves properly; that is, it remains bounded with time. This is illustrated in Fig. 3.1. We next show that: (1) if an improper time step is chosen, then the approximate finite difference solution may become unbounded, and (2) if a proper time step is chosen, then the finite difference solution will behave quite similar to the exact solution. The stability or instability of a numerical scheme will be discussed for a single Fourier wave. This would also be valid for a somewhat more general case, since the total solution is a linear combination of sine and cosine functions, which are all bounded. We next define an amplification factor |λ|, the magnitude of which would determine whether a scheme is stable or not.

Lower-Order Spectral Model

This system was first developed by Lorenz (1960b). It is an elegant system that provides an introduction to the concepts of spectral modeling, based on the use of double Fourier series representations of the basic equations in a doubly periodic domain. Here we examine the barotropic vorticity equation.