The first section of this documnet provides in-depth background to spherical harmonics and the spectral transform technique by discussing technical aspects and terminology essential to a clear understanding and informed usage of fields archived as spectral coefficients by ECMWF. In spectral models, the spectral space representation of fields is tightly coupled to the corresponding physical latitude-longitude transform grid (a Gaussian grid), and vice versa. The second section discusses the ERA-Interim transform grid, which in this case is a reduced N128 Gaussian grid.
over a truncated wavenumber space where is the zonal (east-west) wavenumber, and can be regarded as a "total" wavenumber whereby gives the number of zeroes between the poles (not including the poles) of the spherical harmonic function . Hence, can be interpreted as a type of meridional (north-south) wavenumber. In the spherical harmonic expansion, the coordinate variables are represented by where is longitude, ( being latitude) and time. (We have omitted a vertical coordinate for simplicity of presentation.) In addition, the polynomials defined by
are associated Legendre polynomials of order and degree , and
are complex-valued spectral coefficients where is the complex conjugate operator. Finally, for simplicity of presentation we set , resulting in a single constant spectral truncation parameter which corresponds to what is commonly referred to as triangular truncation. (The usual convention is to use rather than to designate triangular truncation.) In the () wavenumber space prescribes a triangular region of spherical harmonic modes indicated by the filled circles in the diagram below. Modes outside of this triangle are set to (open circles).
|Spherical Harmonic Wavenumber Space|
Other types of truncation may also be used, such as rhomboidal, but triangular truncation is the most common and is used in ERA-Interim (ECMWF, 2002), and for example the NCAR CCM (Kiehl et al., 1998) and CSM (Boville and Gent, 1998). Triangular truncation is frequently referred to as isotropic in the sense that every position and direction on the sphere is treated identically, that is, spectral solutions obtained using triangular truncation are invariant with respect to a coordinate rotation. Washington and Parkinson (1986), and Hack (1992), discuss many aspects of spectral truncation in more detail.
Equation (1) represents the spectral synthesis of a scalar field from a truncated series of spectral coefficients and spherical harmonic functions . Conversely, in a spectral analysis stage, the spectral coefficients are obtained by a discretized version of
where the inner integral (highlighted in light blue)
is a forward Fourier transform applied in the zonal (east-west) direction. The forward Fourier transform is computed at each circle of latitude using a discrete fast Fourier transform (FFT). The outer integral
is evaluated in the meridional (north-south) direction using Gaussian quadrature
where denotes Gaussian grid points in the meridional direction, the corresponding Gaussian weight at point , and the number of Gaussian grid points in the meridional direction. The are given by the roots of the Legendre polynomial and the by
The Gaussian grid points are synonymous with Gaussian latitudes , the relation being , and the number of Gaussian grid points is synonymous with the number of Gaussian latitudes.
The spectral analysis stage represented by Equation (2) can for illustrative purposes be constructed sequentially from two arrays as shown in the following diagram. If we choose equally spaced longitudes and Gaussian latitudes, then to the columns of the array in the left panel we assign the FFTs, one for each circle of latitude. The length of the FFTs is (due to the Nyquist frequency limit, only half of the possible Fourier modes are retained, , and the Fourier "mode" is the mean value of at and time ). In the array in the right panel, we store the spherical harmonic spectral coefficients obtained by Gaussian quadrature of the Fourier coefficients , associated Legendre polynomials , and Gaussian weights . As an example, we show the how the spectral coefficient is computed from the sum of the products of (the in the left panel), , and . (In passing we note that represents the average value of at time and fixed level.)
|Array of Complex Fourier
|Array of Complex Spherical Harmonic
Observe that the dimensions of the spectral coefficient array in the right panel are , and that a triangular truncation has been applied. Also, keep in mind that there are no modes to compute in the region, as for the associated Legendre polynomials . Hence, values of the spectral coefficient array at () indices outside of the light-gray region need not be computed and are set to . In general, the triangular truncation is chosen such that , or where is the number of east-west grid points, to avoid aliasing potentially extraneous small spatial scales onto large spatial scales. (However, this is not to say that modes beyond cannot be or are not computed indeed, modes beyond may be retained depending on the post-processing conventions of the operational, reanalysis, or modeling center, with the implicit understanding that the appropriate truncation was applied during model integrations.)
Once spectral coefficients are obtained during a spectral analysis stage, the equivalent physical representation can be obtained by a spectral synthesis via Equation (1). To follow up on our example, we may choose to focus on only a single mode such that . The real part of the spherical harmonic function is shown in the following figure. Noting that we observe a pattern of alternating negative and positive regions with zonal wavenumber 2 (with nodes at ±45°E and ±135°E) and 3 meridional nodes at 0°N and ±35.3°N. Of course, in a complete spectral synthesis we are summing over many different spherical harmonic modes (multiplied by the corresponding spectral coefficient).
|Real part of Spherical Harmonic Function|
|blue represents negative values, red positive values|
Equations (1) and (2) (spectral synthesis and spectral analysis) constitute a transformation pair in which a scalar variable at time may be represented in physical space, Equation (1), or spectral space, Equation (2). Both are equivalent representations of and each possess context-dependent computational advantages and disadvantages. The transformation pair is bound together by the transform grid, constituting the physical, or grid-point space, defined by the Gaussian latitudes and equally spaced longitudes. Throughout the preceding discussion, we have chosen to emphasize the time-dependent nature of the scalar variable by adhering to notation that retains explicit reference to time, especially in the spectral coefficients . This allows us to conclude this section by discussing the broad outline of a typical forecast (integration) cycle in modern spectral GCMs.
GCMs have as their basis the primitive equations which describe atmospheric dynamics and thermodynamics employing six equations in six unknowns three momentum equations relating the east, north, and vertical components of velocity () to the gradient of pressure (Newton's second law of motion in a noninertial reference frame), a mass continuity equation relating local changes in density to divergence or convergence of atmospheric mass, a thermodynamic equation relating temperature to the storage and conversion of thermal and other forms of energy into work (first law of thermodynamics), and an equation of state relating , , and (the ideal gas law). Through various approximations and combinations, the primitive equations are reduced to a set of model equations that are of two types, these being prognostic equations of scalar variables (these predictive variables are hereafter referred to as prognostic variables), and diagnostic equations of variables that are most conveniently computed from the prognostic variables. In a typical GCM, the resulting prognostic variables are vorticity , divergence , temperature , and surface pressure , with specific humidity added for the inclusion of atmospheric moisture. The horizontal components of both the prognostic and diagnostic sets of model equations are then subjected to the spectral transform method, converting these equations to an equivalent representation in spectral space, i.e. in terms of spectral coefficients.
To examine the outlines of a model forecast cycle let us focus on a single prognostic variable, in this case vorticity . The prognostic equation for vorticity contains linear terms , horizontal derivatives , nonlinear terms (usually products of individual space- and time-dependent terms), and parameterizations of forcing terms . In the spectral method, the forecast cycle begins in grid-point space at time and involves the point-by-point multiplication of the individual parts of nonlinear terms, as well as calculation of physical parameterizations (see lower left panel in diagram below). (Linear operations and horizontal derivatives are to be carried out exclusively in spectral space.) In the second phase of a forecast cycle (upper left panel), nonlinear terms and physical parameterizations are transformed to spectral space. The spectral form of nonlinear terms and physical parameterizations , as well as the spectral form of linear operations and horizontal derivatives , are then used to compute the right-hand side of the prognostic equation for vorticity, i.e. . In the third phase of the forecast cycle (upper right panel), the spectral coefficients are integrated forward in time one time step to time . In addition, from the diagnostic set of equations, the spectral coefficients of the individual parts of nonlinear terms are computed as functions of ().
Finally, the forecast cycle is completed by transforming the individual parts of nonlinear terms from spectral space to grid-point space (lower right panel). In ERA-Interim, at time the spectral form of the prognostic variables and the grid-point form of parameterized and derived fields are archived together which results in some variables being represented by spherical harmonics, and some variables represented on the physical transform grid the ERA-Interim physical transform grid is the topic of the next section2.
2Our discussion of an integration cycle in a spectral model has purposely omitted reference to data assimilation which combines objective analysis and initialization, key elements in ECMWF operational and reanalysis programs. For a brief overview of data assimilation at ECMWF, see pages 462-468 in Holton (1992).
For Gaussian grids in general, one may intuit that as the pole is approached, the curvilinear distance between meridians of longitude decreases markedly (e.g. from 78.170 km just north and south of the Equator to 0.733 km near the poles for regular ), introducing a distortion in the area bounded by the four sides of a grid cell. Because of the isotropic nature of triangular truncation, it has been posited that a prescribed and sufficient grid resolution at the Equator should also serve as a sufficient grid resolution over the entire sphere, alleviating the poleward distortion of grid cells. Beginning with Machenhauer (1979), this premise led to the investigation of reduced Gaussian grids in which the number of longitudes per circle of latitude decreases as one approaches the poles. The operational aspects of using reduced Gaussian grids was more fully developed at ECMWF by Hortal and Simmons (1991), providing the basis for using a reduced Gaussian grid for ERA-40 and a reduced Gaussian grid for ERA-Interim. (Hortal and Simmons, 1991, found no significant loss of accuracy in using the reduced grid for short- and medium-range forecasts of vorticity and geopotential height compared to the full grid. In addition to no significant loss of accuracy, there are real reductions in computational time (20-25%) and storage requirements (30%) for reduced compared to regular at T106 spectral truncation. The reader is referred to Hortal and Simmons, 1991, for further details.)
We provide graphical and numerical representations of reduced and regular Gaussian grids via the following two links. The reduced Gaussian grid is used in ERA-Interim.
Boville, B. A., and P. R. Gent, 1998: The NCAR Climate System Model, version one. J. Climate, 11, 1115-1130.
Bourke, W., 1972: An efficient, one-level primitive-equation spectral model. Mon. Wea. Rev., 100, 683-689.
Eliasen, E., B. Machenhauer, and E. Rasmussen, 1970: On a numerical method for integration of the hydrodynamical equations with a spectral representation of the horizontal fields. Rep. No. 2, Institut for Teoretisk Meteorologi, Københavns Universitet, Copenhagen, Denmark, 35 pp.
European Centre for Medium-Range Weather Forecasts, 2002: The ERA-40 Archive. Reading, ECMWF, 40 pp. (On-line and A4 PostScript versions are available from the ECMWF ERA-40 Project Plan link.)
Hack, J. J., 1992: Climate system simulation: basic numerical and computational concepts. Climate System Modeling, K. E. Trenberth, Ed., Cambridge University Press, 283-318.
Hack, J. J., and R. Jakob, 1992: Description of a global shallow water model based on the spectral transform method. NCAR Tech. Note NCAR/TN-343+STR, 39 pp.
Hack, J. J., B. A. Boville, B. P. Briegleb, J. T. Kiehl, P. J. Rasch, and D. L. Williamson, 1993: Description of the NCAR Community Climate Model (CCM2). NCAR Tech. Note NCAR/TN-382+STR, 108 pp.
Haurwitz, B., 1940: The motion of atmospheric disturbances on the spherical earth. J. Mar. Res., 3, 254-267.
Holton, J. R., 1992: An introduction to Dynamic Meteorology. Second Edition. Academic Press, 511 pp.
Hortal, M., and A. J. Simmons, 1991: Use of reduced Gaussian grids in spectral models. Mon. Wea. Rev., 119, 1057-1074.
Kiehl, J. T., J. J. Hack, G. B. Bonan, B. A. Boville, B. P. Briegleb, D. L. Williamson, and P. J. Rasch, 1996: Description of the NCAR Community Climate Model (CCM3). NCAR Tech. Note NCAR/TN-420+STR, 152 pp.
Kiehl, J. T., J. J. Hack, G. B. Bonan, B. A. Boville, D. L. Williamson, and P. J. Rasch, 1998: The National Center for Atmospheric Research Community Climate Model: CCM3. J. Climate, 11, 1131-1149.
Machenhauer, B., 1979: The Spectral Method. Numerical Methods used in Atmospheric Models, Vol. II, GARP Publication Series 17, World Meteorological Organization, Geneva, 121-275.
Orszag, S. A., 1970: Transform method for calculation of vector coupled sums: Application to the spectral form of the vorticity equation. J. Atmos. Sci., 27, 890-895.
Sardeshmukh, P. D., and B. J. Hoskins, 1984: Spatial smoothing on the sphere. Mon. Wea. Rev., 112, 2524-2529.
Swarztrauber, P. N., 1993: The vector harmonic transform method for solving partial differential equations in spherical geometry. Mon. Wea. Rev., 121, 3415-3437.
Temperton, C., 1991: On scalar and vector transform methods for global spectral methods. Mon. Wea. Rev., 119, 1303-1307.
Washington, W. M., and C. L. Parkinson, 1986: An Introduction to Three-Dimensional Climate Modeling. University Science Books, 422 pp.