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.
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
with
and 
Then
This equation will be the starting point for the derivation of the divergence and vorticity equations.
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
where
and
Since
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.
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.
and find
where
is the horizontal divergence in layer k,
is the kinetic energy per unit mass in layer k,
is the absolute vorticity,
and
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.
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)
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
where
(Subroutine EPSILON)
When horizontal gradients require north-south derivatives, they can be obtained from
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.
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
and velocity potential 
and write
Clearly, all the results in terms of and can be expressed in terms of 
and D since
and
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,
Similarly,
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
where
then
Thus,
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.
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.
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:
where
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
where
Wj are the Gaussian weights (2.43), and N is the number of Gaussian latitudes. Eq. (2.73) is used to forecast the vorticity.
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
then
(subroutine FL22), where the kinetic energy was first computed on the Gaussian Grid and then decomposed as
c) The linear terms.
then
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.
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.
is computed on the Gaussian grid and coefficients 
are obtained, then
The surface pressure equation in spectral from is, therefore,
where
and
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.
and
Then
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.
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
while
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
Let
Then
and
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
and
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.
, 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:
Letting
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
and
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:
Let
Then
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.
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
where
Defining
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.
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
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 .