1. Basic Model Equations
The model is based on the usual expressions of conservation of mass, momentum, energy and moisture. In order to take advantage of the spectral technique in the horizontal, we replace the momentum equations by the vorticity and divergence equations (Bourke, 1972), thus eliminating the difficulties associated with the spectral representation of vector quantities on a sphere (Machenhauer, 1972).

The vertical coordinate sigma, , is described by Phillips (1959) and is depicted in Fig.2.1. Differential operators in this coordinate are implemented by finite differences where interface values are assumed to be averages of their bracketing layers, which ensures quadratic conservation (Arakawa, 1972). After the vertical discretizations are introduced into the governing equations, the derivation of the horizontal spectral form of the equations can be carried out using the appropriate equations in flux form.

Before deriving the detailed spectral equations, it is useful to present the basic equations after discretization in s. A list of symbols appears in Appendix 2B.

  1. Momentum Equation
The continuous momentum equation is

where the I operator is the horizontal gradient in the s system and represents dissipative processes in the model. We expand the total time derivative in sigma coordinates, and make use of vector analysis identities to write

The vertical derivative is written in finite differences using

and approximating the derivative by




This equation will be the starting point for the derivation of the divergence and vorticity equations.

  1. Thermodynamic Equation
The thermodynamic equation

where H is the heating rate per unit mass, is written in flux form as in Eqn (2.4). Applying the vertical finite differencing approximation as in Eq. (2.4) and assuming

we have



  1. Surface Pressure Equation
This equation is obtained from the continuity equation by integration over the full sigma domain:


we have

At this point it is natural to introduce the vertical velocity equation. This equation is obtained from the continuity equation discretized in the vertical,

and substituting the time derivative from Eq. 2.15:

This recursive relation starts from(surface)=0 and satisfies (top)=0.

  1. Water Substance Equation
We assume that the mixing ratio q is conserved subject to sources and sinks represented by S. Then

Transforming Eq. (2.18) into flux form as in Eq. (2.3) and performing the vertical derivative approximation yields

This completes the presentation of the basic equations discretized in the vertical. The differentiated momentum equation for divergence and vorticity are presented next.

  1. Divergence Equation
We operate on Eq. (2.1.6) with and find


is the horizontal divergence in layer k,

is the kinetic energy per unit mass in layer k,

is the absolute vorticity,


The purpose of linearizing the temperature about a sigma dependent basic state profile is to facilitate a semi-implicit time integration. Eq. (2.20) can also be written as

where the term Xk represents the nonlinear terms

and the Laplacian in Eq. (2.28) contains linear terms only.

The geopotential Fk appearing in Eq. (2.21) is diagnosed from the temperature field using the Arakawa (1972) scheme for the hydrostatic equation (Appendix 2A):

for k = 2,...K, and

where is the surface geopotential. The system (2.30-2.31) is solved in the model for the geopotentials by matrix inversion.

  1. Vorticity Equation
We operate on Eq. (2.6) with and find

where Ak, Bk are given by Eqs. (2.25) and (2.26). It should be noted that Eq. (2.32) does not depend on either the full temperature or its departure from the basic state. This feature is used in the code implementation.

Summarizing, the model prognostic equations are:

Temperature Eq. (2.1.9)
Surface pressure Eq. (2.15)
Moisture Eq. (2.19)
Divergence Eq. (2.20)
Vorticity Eq. (2.32)

The diagnostic equations are:

Hydrostatic Eqs. (2.30) and (2.31)

Vertical velocity Eq. (2.17)

  1. Overview of the Spectral Technique
    1. Definition of Expansion in Spherical Harmonics
In the following, all prognostic variables as well as the tendencies will be assumed to possess a spherical harmonic representation of the form

where the spherical harmonics are given by

( computed in subroutine PLN.)

The upper limit of the second sum is left as a general function of l. The model can be integrated with any resolution, and the specific resolution is determined by a control array, set at the beginning of the integration.

The expansion coefficients are functions of time and sigma and the are Legendre functions of the first kind as defined in the list of symbols and given by

where x=sin and the normalization of the is such that

The can be generated from the recursion relation


(Subroutine EPSILON)

When horizontal gradients require north-south derivatives, they can be obtained from

  1. Grid to Spectral Transform
(Subroutine GOZRIM, or "let's differentiate" in Hebrew)

It is now assumed that grid values of a given field are known, and it is required to compute the field's expansion coefficients. In practice, the field could represent a history variable or a tendency, in which case the field will be computed with quadratic accuracy.

Let the field be represented by Eq. (2.33). Multiplication of this equation by and a surface integration over the sphere yields:

The numerical evaluation of the integrals in Eq. (2.40) is carried out in two steps:

Step 1.

We define the field's Fourier coefficients at a given latitude as

This integral can be evaluated using a discrete Fourier transform, provided F is a trigonometric polynomial.

In general, if f(x) is a trigonometric polynomial of degree not exceeding M-1

is an exact evaluation of the integral (Krylov, 1962). Since the model's variables are assumed to have a spherical harmonic expansion, they are represented in the zonal direction by a trigonometric polynomial of degree J. Quadratic terms will contain pow ers of at most 2J, and the integrand in Eq. (2.41) will be of degree not exceeding 3J. For exact integration in the zonal direction, we therefore require at least 3J+1 points around a latitude circle.

Step 2.

The integration in the meridional direction is performed by Gaussian Quadrature (Krylov, p.108).

The Gaussian weights are

where the Xk are the zeros of for some N which will be determined later. Then

is an exact integration of the function y, provided y=y(x) is a polynomial of degree not exceeding 2N-1, and the function y is known at the points, see Fig. 2.2 and Table 2.1.

Using this method

For a rhomboidal truncation, N(5J+1) / 2 , while for a triangular one, N(3J+1) / 2. For an arbitrary resolution, the maximum power of must be evaluated, and the condition imposed by the Gaussian Quadrature requirement must be applied.

  1. Spectral to Grid Transform
The computation of grid values from the spectral expansion given by Eq. (2.33) is carried out in two steps. First, for real F, we take advantage of the relation

to derive

Here, the limit N(l) is equal to J for triangular truncation, and to J+l for rhomboidal truncation. It is assumed that grid values are required at a given latitude (corresponding Xk). The Legendre polynomials can be computed at this latitude, and the partial sums

are evaluated. Eq. (2.47) can, therefore, be written as

If equally spaced values of F are required at the latitude , a discrete Fourier Transform can be applied to Eq. (2.49).

with inverse

  1. Velocity and Divergence-Vorticity Relations
While the kinematic history variables in the model are divergence and vorticity, the wind velocity is required, too. In order to compute the spectral coefficients of the scaled velocity, we use the decomposition of a vector field into its rotational an d divergent components. We introduce the stream function and velocity potential and write

Clearly, all the results in terms of and can be expressed in terms of and D since


We now assume that and are expressed by

Substitution of the above series into the U equation (2.53) and (2.49) yields

since can be expressed as a linear combination of and, , (2.39) U can be expressed as a spherical harmonic series. We find, after some algebra,


It is noted that for each zonal wave number, the meridional derivative introduces an additional term. Thus, the series of U and V are

It is pointed out that the above method is not the only way to compute the velocity from the divergence and vorticity. It is possible, at each Gaussian latitude, to use the Fourier coefficients of and (or and D) and to compute the Fourier coefficients of U and V directly from Eq. (2.58).

If we let




are the Fourier coefficients of U at the given latitude. A similar expression can be found for the V component. The above method is useful when computer memory is limited and the usual trade off between memory and computations applies. In this case, the computation of the Fourier coefficients of U or V requires two summations per latitude, while the method using the spherical harmonic expansion of U or V requires only one summation.

  1. Spectral Form of the Gradient
At this point, it is convenient to introduce the computation of the gradient of a function F:

Since is easily computed (2.63), it is convenient to compute

Introducing a spherical harmonic series for F, we find

where, again, the meridional derivative introduces for each l an extra term in the expansion.

  1. The Spectral Form of the Model's Equations
    1. The Vorticity Equation
Equation (2.25) is transformed into spectral form by assuming that all prognostic variables have spherical harmonic expansions that can be evaluated on the Gaussian Grid, while and possess Fourier series representations at each latitude of the Gaussian Grid, see Fig. 2.2 .

The tendency of a particular vorticity spectral component is obtained from Eq. (2.25) by multiplication with and integration over the sphere:


The zonal integration is performed first, and integration by parts with respect to yields

Eq. (2.71) is integrated by a Gaussian Quadrature, and the final form may be written as


Wj are the Gaussian weights (2.43), and N is the number of Gaussian latitudes. Eq. (2.73) is used to forecast the vorticity.

  1. The Divergence Equation
Eq. (2.20) is transformed into spectral form in three steps:

a) The terms

are treated in the same fashion as the corresponding terms in the vorticity equation. The spectral form of this term is, therefore, (subroutine MSUM)

b) The Laplacian of the kinetic energy is accumulated separately, thus, if we denote


(subroutine FL22), where the kinetic energy was first computed on the Gaussian Grid and then decomposed as

c) The linear terms.


The nonlinear part of Eq. (2.28), the divergence equation as written for the semi-implicit computation, may, therefore be written as

The linear part will be further manipulated in the semi-implicit computation.

  1. The Thermodynamic Equation
Before considering the spectral form of the thermodynamic equation (2.9), we linearize it for a semi-implicit time integration.

We consider Eq. (2.9) and let

In the above equation, we have defined, for coding convenience

and the vertical integral of the divergence, the horizontal advection of log , and its vertical integral are given by:

We now group the terms that are linear in D. This is achieved by noting that

The linearizable terms are the last two terms in Eq. (2.82).

If we define the matrix B with its row Bk as

where represents the vector of components Dj. Note that Bk is independent from horizontal location. It is computed once in subroutine BMCM. We now group the nonlinear terms in

We may write

The spectral form of the last equation is obtained by computing the nonlinear terms Yk on the grid and subsequently expressing them in a spherical harmonic series

The final form of the thermodynamic equation is

This form will be used in the section describing the time integration scheme.

  1. The Surface Pressure Equation
The nonlinear component of (2.15)

is computed on the Gaussian grid and coefficients are obtained, then

The surface pressure equation in spectral from is, therefore,



  1. The Moisture Equation
Equation (2.19) is written in the form

The R.H.S. of Eq. (2.96) is computed using the transform method and when a full scan of the Gaussian grid is completed, the moisture tendency is available. Note that, as with he thermodynamic equation, we have separated the horizontal advection of moisture into two terms. Then the moisture tendency has two contributions:

The tendency is computed from two contributions.



The reason for separating the full moisture tendency into two parts is the ease of implementing the divergence operator in spectral form using the transform method. It is pointed out that this operator must be implemented anyway in order to compute the vorticity tendency.

  1. The Time Integration Scheme
In order to employ a longer time step, a semi-implicit time integration is applied to the divergence, temperature and surface pressure equations. The vorticity and moisture equations are integrated explicitly.

The coupled system of equations can be written as

In the above, x, y and z are the nonlinear contributions to the divergence, vorticity and the natural log of the surface pressure. Horizontal diffusion is introduced as a numerical smoothing effect in the form of double Laplacian operator in Eqs. (2.90) and (2.91). Kd is a divergence diffusion coefficient and is set to


plays that role in the thermodynamic equation and is set to .

The matrix B is defined by grouping together all the linear terms depending on the divergence, and the column is the basic state temperature, chosen to be isothermal and equal to 300oK in all layers.

It is now assumed that all pertinent variables possess a spherical harmonic expansion, and the double Laplacian of temperature is approximated on pressure surfaces by

The coupled system, in spectral form, may now be written as

where for brevity the spectral indexing was omitted. In the above, we defined matrices

and the vectors

where the are the sigma thicknesses of the model's layers, and the matrix A is the hydrostatic metric

It is noted that the term has now been absorbed in the nonlinear term X

We now proceed to integrate in time the coupled system, using a centered time differencing scheme. The linear terms are treated implicitly. If we designate

Then, we have

and in terms of variables

Similarly, for the thermodynamic equation

while the surface pressure equation yields




Substitution of into the divergence equation yields

collecting terms in and substituting for from the surface pressure equation

The right hand side of the last equation is known, and Dm turns out to be quite regular, we can therefore solve for . is computed from

Once is known, , Q+ and , more easily computed.

The time integration of the vorticity and moisture equations is performed on the equations


where W is the full vorticity tendency, R is the part of the tendency, resulting from advection, surface fluxes and vertical diffusion while Sq represent moisture phase changes. The numerical implementation of the moisture equation uses a splitting method where the preliminary moisture prediction is computed from the contribution of advection, surface fluxes, horizontal and vertical diffusion. These preliminary values are then compacted, in conjunction with the temperature field, to account for phase changes due to convective and large scale precipitation. Shallow convection is also modeled during the second segment of the time integration.

  1. Normal Modes Initialization
A normal mode of the model's linear version is a solution of , D, T, and lnp* having the same temporal behavior. The number of such modes is a function of the model's resolution, vertical and horizontal. In order to determine all possible model oscillations, a linear model version must be first determined and cast in spectral form. For separability of variables, we consider perturbations about a resting atmosphere. The basic state vertical temperature profile is specified by a U.S. Standard Atmosphere.

We begin with the linearized vorticity equation:


equation (2.113) can be written as

The spectral form of equation (2.115) is obtained by substituting the spherical harmonic expansions of the dependent variables, multiplying by the conjugate harmonic of the desired component and integrating over the sphere. After some manipulation, the spectral linear vorticity equation in each model layer reads:

We now turn to the linearized divergence equation

The geopotential F in equation (2.117) does not include its orographic component. The effects of mountains are introduced during the initialization through the nonlinear part of the divergence tendency.

In order to simplify the vertical separation of variables we define the composite variable

Using the notation in equation (2.114), the linear divergence equation may be written

Inspecting equations (2.115) and (2.116) and recognizing their analogy with equation (2.119), the spectral form of the linearized divergence equation may be written for each model layer:

To close the system of equations, the linear thermodynamic and continuity equations are needed. These are


Equations (2.121) and (2.122) can now be used to express the composite variable W in terms of T and D. Consider the tendency of W. Using the definition (2.118) we may write

Substituting from equations (2.121) and (2.122) and using the hydrostatic relation , equation (2.123) takes the form

The definition of the composite variable W thus allows us to deal with only one vertically coupled equation. Let

Since is linear in , we may write.

Equation (2.125) closes the system of linearized model equations.

In order to effect a separation of variables it is natural to expand all pertinent vertical column vectors in the basis vectors defined by the eigenvectors of the vertical coupling matrix G. Let = (F1, ..., Fk) be an arbitrary vector and let , with j=1, ..., K, be the eigenvectors of G.

We write the eigenvector expansion

and need to specify the coefficients aj. Since G is not generally symmetric, the are not always orthonormal; they are, however, orthogonal to the eigenvectors , of , if the eigenvalues of G are distinct. Therefore,

Equations (2.126) and (2.127) will be used to separate the linear solutions of equations (2.116), (2.120), and (2.125) as well as to express tendencies in the initialization procedure. To complete the normal modes computation, let -(gh) be the eigenvalues of G . Equation (2.125) may then be written as

and the separation of equation (2.125) has been accomplished.

We now note that equations (2.116), (2.120), and (2.128) are invariant under a substitution of all column vectors by their eigenvector expansions. The resulting equations hold for the expansion coefficients. If we assume t hat a given mode behaves in time according to eiwt, we may continue to deal with equations (2.116), (2.120), and (2.128) and consider their solutions multiplied by e, the model's normal modes.

The matrices arising from equations (2.116), (2.120), and (2.128) are not symmetric and, therefore, do not possess orthonormal eigenvectors. The scaling

leads to symmetric matrices relating the primed variables. Applying this scaling to equations (2.116), (2.120), and (2.128), we find

Equations (2.5.19), (2.5.20), and (2.5.21) may be written in matrix form:



The horizontal eigenvalue problem may now be written as

Note that matrix E depends on the eigenvalues of G and, therefore, equation (2.152) must be solved for each eigenvalue of matrix G. The symbol -(gh) is chosen to reflect the analogy of the shallow water problem with a nonrotating earth version of the linear model. For W=0, equations (2.5.7) and (2.5.14) combine to produce the wave equations with speed . It is pointed out that for positive eigenvalues of G the semi-implicit computation would be unstable.

At this point, the tools required to implement a nonlinear normal modes Machenhauer initialization scheme have been assembled. The procedure consists of the following steps: Initial data from the Optimal Interpolation analysis are transferred to the model, and the full tendencies of the model's variables are calculated. These tendencies are transformed into variables analogous to those used in the normal modes computation, namely the composite variable W, and scaled, primed variables. At this stage, the transformed tendencies are expanded in terms of normal modes, precomputed in permanent files, and the changes to gravity modes with periods less than 2 days are computed. These changes, , where Gm represents a given gravity mode with period are transformed back to the model's variables and are used to adjust the input data. This process does not ensure the vanishing of gravity modes' tendencies in a subsequent tendency computation due to the nonlinear nature of the problem. It was experimentally determined that two iterations using four vertical modes produce acceptable changes in the initial conditions and very smooth surface pressure time integrations.

In order to capture the effects of diabatic heating, the model is first integrated for seven time steps, and the temperature change during that time is translated into a heating rate. This heating rate is added to the tendency in the thermodynamic equation during the initialization procedure. A typical behavior of the surface pressure at a point is depicted in Fig. 2.3.

Table 2.1

Co-latitudes and gaussian weights for rhomboidal 40 and triangular 80 spectral truncation.


The Arakawa vertical finite differencing scheme is derived by writing the energy equation in the continuous and discretized forms. The requirement that both forms be similar leads to the finite differencing form of the hydrostatic equation.

We begin with the continuous equations. The momentum equations in flux form in the sigma coordinate system are

The vertical finite differencing considerations are independent of the choice of horizontal coordinates. Using the continuity equation

the kinetic energy equation is derived:

In equation (2.A.4) is the kinetic energy per unit mass, p is the pressure divided by 100 cb, and its substantial derivative.

To complete the energy equation, we add the theromodynamic equation,

to the kinetic energy equation (2.A.4) to obtain

In the last equations, . The vanishing of the right-hand side in equation (2.A.6) is a consequence of the hydrostatic assumption.

We now turn to the discretized equations using, at this point, values at interfaces as well as layers. For layer k, we express (2.A.1) and (2.A.2) as



and multiplying the thermodynamic equation by , we have

The total energy equation is now obtained by combining equations (2.A.7), (2.A.8), and (2.A.11) in a manner similar to the continuous case; it reads

Equation (2.A.12) would be analogous to the continuous energy equation (2.A.6) if the right-hand side vanished for all layers, and for all possible sk. This would occur if

and for k = 1, ..., K

It can be shown that following Arakawa and writing , will be conserved in time, with similar results for other variables so related. With this specification for the interface potential temperatures, equations (2.A.13) and (2.A.14) for k = 2, ..., K, reduce to

and are easily eliminated from equations (2.A.16) leaving a relation between layer temperatures and layer geopotentials. To close the system, equation (2.A.15) is manipulated using the Brown (1974) and Phillips (1974) specification of layer pressures:

to obtain, for k = 1, ..., K,

When equation (2.A.18) is summed over all layers and identified as the fixed geopotential at the ground, we have

thus completing the finite difference version of the hydrostatic equation. The transformation of the equations from flux form to the form used in the model's formulation is straightforward.


a = radius of the earth

Ek =

Cp = specific heat at constant pressure for moist air

Dk = divergence

E = 1/2 ( u2 + v2 ) , kinetic energy per unit mass

f = Coriolis parameter

g = acceleration of gravity

= vertical unit vector

k = index of vertical layer

K = total number of layers

l = zonal wave number

n = meridional index

= pressure

= Legendre function of the first kind

= surface pressure

R = gas constant

q = specific humidity

t = time

T = temperature

u = zonal velocity

U =

v = meridional velocity

V =

w = Gaussian quadrature weight

z = height

= latitude


l = longitude

c = relative vorticity

f = geopotential height

= geopotential height of the earth's surface

x = velocity potential

q = potential temperature

y = stream function

W = angular velocity of the earth

= value defined at interface between layers

= ; finite difference vertical integral

Dj = layer sigma thickness


Arakawa, A., and Y. Mintz, 1974: The UCLA Atmospheric Circulation Model, Department of Meteorology, University of California.

Ballish, A. B., 1980: Initialization, theory and application to the NMC spectral model. Doctoral Dissertation, University of Maryland.

Bourke, W., 1974: a multi-level spectral model. Formulation and hemispheric integrations. Mon. Wea. Rev., 102, 687-701.

Brown, J. A., 1974: On vertical differencing in the sigma system. NMC. Office Note 92, U. S. Department of Commerce, NOAA, NWS, Washington, DC.

Eliasen, E., B. Machenhauer, and E. Rasmussen, 1970: On a numerical method for integration of the hydrodynamic equations with a spectral representation of the horizontal fields. Institute for Theoretical Meteorology, University of Copenhagen, Report No. 2, Copenhagen, Denmark.

Krylov, V. I., 1962: Approximate calculation of integrals, ACM Monograph Series, The MacMillan Company, New York, NY.

Kuo, H. L., 1965: On formation and intensification of tropical cyclones through latent heat release by cumulus convection. J. Atmos. Sci., 22, 40-63.

Machenhauer, B. and E. Rasmussen, 1972: On the integration of the spectral hydrodynamical equations by a transform method. Institute for Theoretical Meteorology, University of Copenhagen, Report No. 3, Copenhagen, Denmark.

Machenhauer, B., 1977: On the dynamics of gravity oscillations in a shallow water model, with applications to normal mode initialization. Beitrage zur Physik der Atmosphere, 50, pp. 253-271.

Orszag, S. A., 1970: Transform methods for calculation of vector coupled sums: Application to the spectral form of the vorticity equations. J. Atmos. Sci., 27, 890-895.

Phillips, N. A., 1959: Numerical integration of the primitive equation on the hemisphere. Mon. Wea. Rev., 87, 333-345. _______________, 1974: Application of Arakawa's energy conserving layer model to operational numerical weather prediction. NMC Office Note 104, August 1975, U. S. Department of Commerce, NOAA, NWS.

Robert, A. J., 1965: The behavior of planetary waves in an atmospheric model based on spherical harmonics. Arctic Meteor. Res. Group, McGill University Publ. Meteor., No. 77, pp. 59-62.

_____________, 1969: Integration of a spectral model of the atmosphere by the implicit method. Proc. WMO/IUGG Symposium on Numerical Weather Prediction, Tokyo, 26 November - 4 December, 1968. Japan Meteorological Agency, Tokyo.

Fig. 2.1 Vertical structure of the 18 layer model . Layer pressure , p , in millibars .

Fig. 2.2 (a) The Gaussian grid. (b) The spectral truncation.

Fig. 2.3 Surface pressure ( mb ) during 0-12 hour forecast at a Gaussian grid point with ( smoothest line ) and without initialization .