A. Equations of motion

The basic physical laws that govern atmospheric motions in the COLA GCM are the conservation laws for mass (dry air and moisture conserved separately), energy, and angular momentum expressed in the spherical shell geometry and rotating frame of reference of the Earth’s atmosphere. The equations that are employed include the equations of motion, continuity equations for dry air and water vapor, and the first law of thermodynamics. The equation for the vertical component of velocity (parallel to the gravitational force vector) is replaced by a diagnostic relationship on the assumption that, on the spatial and temporal scales of interest, motions are in hydrostatic balance. The complete set of equations which describe such a fluid system are termed the primitive equations. The model is global in extent and conforms to a spherical shell geometry whose vertical extent is sufficiently small that the distance from the center of the Earth may be assumed to be approximately constant (equal to the radius of the Earth). The distance above the surface of the Earth is included as an independent variable.

Because the atmosphere is assumed to be in hydrostatic balance, it is convenient to express the vertical variation in the primitive equations in terms of a pressure coordinate. In practice, as recommended by Phillips (1957), the vertical coordinate is pressure, p, normalized by its value at the surface, p_s:

so that the vertical coordinate, s , varies from a value of 1 at the Earth's surface to a value of 0 where the atmospheric pressure effectively vanishes. We write the conservation of angular momentum as follows:

where V is the horizontal wind vector, R is the dry air gas constant, T is the temperature, F is the geopotential, f is the Coriolis parameter, k is the unit vector in the vertical direction, and F is the vector sum of all frictional forces. This equation indicates that the acceleration of horizontal atmospheric motion results from the sum of the pressure gradient force, the Coriolis force due to the rotating frame of reference, and the sum of all frictional forces. By applying vector identities and the curl and divergence operators to this equation, we obtain the scalar vorticity and divergence equations:

where h is the absolute vorticity (h =z +f, where z is the relative vorticity), a is the Earth's radius, f is the latitude, l is the longitude, T₀ is the global mean temperature and A, B and E are defined as follows:

where U = u cos f , V = v cos f , u and v are the components of the horizontal wind vector in the longitudinal and latitudinal directions, respectively, T' is the deviation of temperature from its global mean T₀, Fl and Ff are the components of the friction force in the longitudinal and latitudinal directions, respectively, and E is the kinetic energy.

The first law of thermodynamics is written

where q is the potential temperature (q = T(p/p₀)k where p₀ is a standard pressure and k = R/c_p), Q is the heating rate, and c_p is the specific heat of air at constant pressure. By expanding the total derivative and rearranging terms, we can write the first law as:

a form that is more suitable for vertical differencing (see below).

The mass conservation law is expressed mathematically as the continuity equation. This equation may be vertically integrated to obtain an equation for the evolution of the surface pressure:

Note that the upper and lower boundary conditions of no flow normal to the ground or through the "top" of the atmosphere have been applied. Likewise, a conservation equation may be expressed for water vapor in the atmosphere given a set of sources and sinks (represented by S):

where q is the specific humidity.

The equations (3), (4), (9), (10) and (11) are the continuous form of the equations used to predict vorticity, divergence, temperature, surface pressure and atmospheric water vapor in the GCM. The geopotential that is needed in the divergence equation is obtained from vertical integration of the hydrostatic equation:

and, similarly, the vertical velocity is diagnosed by vertical integration of the divergence from the surface to the level in question:

While many texts provide derivations of these equations in a spherical shell geometry, the most relevant reference for the equations used in the COLA model is Sela (1980) or, alternatively, MRF88 (chapter 2). In addition to the above prognostic variables, the model also diagnoses clouds in four categories at all levels (see section IIC), sea ice temperature, and the land surface variables used in the biosphere model, namely, deep soil temperature, ground surface temperature, canopy temperature, three stores of soil moisture, ground water storage (ponding) and canopy water storage (either liquid or solid, i.e., snow) are also predicted (see section IIIB).

B. Model discretization

The equations described above are solved numerically. In order to place the equations in a form suitable for numerical solution, the four independent dimensions of the problem must be expressed in discrete coordinates. The time dimension is divided into evenly spaced time steps and differencing is carried out by means of a semi-implicit scheme described below. Variations in the horizontal dimensions are represented as coefficients of projection onto an orthonormal set of basis functions that are particularly appropriate for fluid motions in a spherical geometry, namely, the spherical harmonics.

The spherical harmonics are the eigenfunctions of the spherical form of Laplace's equation that governs tidal motions. Because the atmosphere is well stratified and motions are assumed hydrostatic, the surface spherical harmonics are used instead of the full three dimensional spherical harmonics. This affords analytic evaluation of derivatives in the horizontal directions. The application of this technique to atmospheric motion was suggested by Platzman (1960) and made practical for reasonable resolutions by means of the spectral transform method (Orszag, 1970). The spherical harmonics may be written as the product of Fourier components and associated Legendre functions. The Fourier components (e^iml ) represent variations in the longitudinal direction, being a natural choice given the periodic nature of the spherical shell geometry in this direction. Variations in the latitudinal direction are represented by the associated Legendre polynomials (P_mⁿ) expressed as functions of the sine of latitude. Because the associated Legendre polynomials satisfy the geometrically determined boundary conditions at the poles and because they satisfy Laplace's equation governing tidal motions on the sphere they are a natural choice for representing latitudinal variations. The spherical harmonic of order m and degree n is the product of (P_mⁿ) and (e^iml ). Any piecewise continuous scalar field on the surface of a sphere may be represented as a series of these orthogonal spherical harmonics. In practice, the series is truncated at a given two-dimensional wave number, with higher order terms in the series being ignored.

The spherical harmonics which are retained in the truncated series may be represented in a two dimensional (m,n) wave number space, where m is the zonal wave number and n-m is the number of zeroes of P_mⁿ, and the truncation may have various shapes. The two most commonly used truncations are triangular in which the sum of the zonal and latitudinal limits of truncation is a constant and rhomboidal in which the two truncation limits are the same. The triangular truncation has the advantage that all basis functions having the same scale are either included or dropped. The rhomboidal truncation was the first one adopted for numerical treatment of the primitive equations because it was easier to code for optimal performance on fast vector supercomputers (Bourke, 1972).

The spectral method of solution, in which prognostic variables are expressed in series of orthogonal basis functions, has the drawback that nonlinear terms in the equations involve quadratic and triple products (and therefore integrals of quadratic and triple products) which require a great deal of computation which exceeds available computer resources at high model resolution. Orszag (1970) suggested a method for making this computation feasible. By representing the quadratic terms on a grid to which the prognostic variables have already been transformed, the contributions of these terms to the time tendencies could be computed using the back transform. By using suitable numerical algorithms to gain computational efficiency (fast Fourier transform) and retain the exactitude of the spectral treatment (e.g. by using Gaussian quadrature in the integral with respect to latitude), it was shown that the spectral transform method is computationally superior to the direct spectral method (Orszag, 1970; Bourke, 1972). The grid on which the most computationally efficient transform may be effected is called the Gaussian grid to recall the form of quadrature that is employed. A more complete and rigorous development of spherical harmonics, the spectral method and the spectral transform method may be found in Haltiner and Williams (1980).

The elevation of the Earth’s surface is represented by interpolating the Navy 10’ by 10’ observed topography to the model Gaussian grid, iteratively filtering and transforming the result to its spectral representation (Fennessy et al., 1994).

Vertical derivatives are evaluated by finite differencing in s . The vertical differencing is written in an energy conserving form (Arakawa, 1972). Terms involving derivatives with respect to s are expressed in flux form, e.g.,

The first term may then be approximated as:

where the subscript k refers to the k^th vertical (discrete) level. The hat is a vertical averaging operator:

As a result, terms involving vertical derivatives in equations (3), (4), (9), (10) and (11) may be expressed in finite difference form by analogy with the above. The thermodynamic equation is rearranged for energy conservation as follows:

The arrangement of variables on the discrete vertical levels is shown in Fig. 1. A full derivation of the vertical finite difference equations is given in MRF88 (chapter 2).

Spectral truncation of the prognostic mass variable, log of surface pressure, results in non-conservation of atmospheric mass. In the COLA GCM, the surface pressure is adjusted uniformly over the globe at periodic simulated intervals (typically one month) to ensure that the total mass of the atmosphere remains constant. The spectral truncation of specific humidity results in negative values of this variable, under certain circumstances. A scheme (ECMWF, 1988) in which moisture is borrowed from lower atmospheric layers is employed to restore negative values of specific humidity to zero. Both conservation restoration schemes are described in greater detail in Schneider and Kinter (1994).

The time domain of the model is expressed in discrete steps of a fixed length. Given a set of initial conditions for the model prognostic variables, their time evolution is computed by "marching" forward using the time derivatives (as expressed in equations 1 - 5). Since the equations describe fluid flow which supports wave motions, the size of the time step is constrained by the speed of the fastest such wave due to the Courant-Friedrichs-Levy (CFL) computational instability associated with aliasing that wave (e.g. Haltiner and Williams, 1980). Due to the presence of fast gravity waves in the atmosphere, this constraint can result in an extremely small time step that is prohibitive in terms of computational speed, so some means of overcoming this limitation needs to be employed. In many GCMs, the method of choice is to treat the terms in the equations associated with gravity waves separately from the other terms. In particular, the pressure gradient term in the momentum equation (or its analog in the vorticity and divergence equations) and the divergence term in the surface pressure tendency equation are stepped forward in time implicitly while the other terms are treated explicitly (e.g., leap frog scheme). An implicit time differencing scheme is one in which the value of a variable at a future time appears on the right hand side of the equation, so that a system of simultaneous equations (whose order is the number of spatial degrees of freedom) must be solved at each time step. Its advantage is that it is unconditionally stable (no computational instability), so that the time step can be arbitrarily large. Its disadvantages are the huge computational expense of inverting a matrix at each time step and the phase error introduced by the differencing scheme which grows with the size of the time step (Haltiner and Williams, 1980). Thus, a balance must be struck between computational efficiency and numerical accuracy. A scheme in which part of the equation is treated implicitly is called semi-implicit. The COLA model uses such a semi-implicit scheme for the divergence, continuity and hydrostatic equations. A full derivation of the semi-implicit divergence, temperature, and moisture equations is given in MRF88 (chapter 2).

A drawback of the explicit leap frog scheme used to march the vorticity and moisture equations is the tendency for the solution to diverge at alternative time steps. To suppress the time splitting, a Robert (1969) time filter is employed.

C. Horizontal diffusion

Horizontal diffusion is necessary to control small scale "noise" which would otherwise arise in the model. The origins of the noise include (i) the effects of a finite spectral truncation, which interrupt the downscale cascade of energy, (ii) small scale gravity waves caused by the sub-grid scale physical processes, and (iii) other, purely computational effects. Since the larger, well resolved scales should not be affected, a scale selective biharmonic type diffusion is utilized. Experience at NCEP has shown that for the variables of moisture and temperature, this diffusion must be along pressure surfaces and not along the model s surfaces, for otherwise excessive precipitation in regions of steep terrain results (MRF88). For vorticity and divergence, diffusion along s surfaces seems to be adequate.

The relationship between a horizontal derivative on constant pressure surfaces of any variable Y and the corresponding derivative on constant s surfaces is:

where p is the log of surface pressure, x is one of the horizontal directions, and Y is an arbitrary scalar field. The subscripts p and s refer to derivatives along constant pressure and s surfaces, respectively. Applying (19) twice in both horizontal directions, we obtain:

Experiments have shown that an adequate (and computationally efficient) approximation to the full equation (20) is given by retaining only the first and third terms on the right hand side of (20), and by replacing Y with its horizontal average (denoted by an overbar) in the third term. Expressing Y in spectral form, the (m,n) component of the diffusion along pressure surfaces may be written:

where a is the radius of the Earth. Applying (21) twice, and again assuming constant vertical derivatives, we obtain:

 The above formulation (22) is used for temperature and moisture. For vorticity and divergence, the vertical derivative term embodying the corrections to ensure diffusion along a constant pressure surface is omitted.

The diffusion coefficient is set to a value that causes the highest resolved total wavenumber component to damp to 1/e of its initial value in T hours in the absence of all other processes. For divergence, the value of T is set to 21 hours, and for temperature, vorticity and moisture, the corresponding time is 28 hours.