IFS Documentation front page

I Observations
II Assimilation
III Dynamics
IV Physics
V Ensemble
VI Technical
VII Waves

2.2 Discretization

2.2.1 Vertical discretization

To represent the vertical variation of the dependent variables , , and , the atmosphere is divided into layers. These layers are defined by the pressures at the interfaces between them (the `half-levels'), and these pressures are given by

(2.11)

for . The and are constants whose values effectively define the vertical coordinate and is the surface pressure field.

The values of the and for all are stored in the GRIB header of all fields archived on model levels to allow the reconstruction of the `full-level' pressure associated with each model level (middle of layer) from ()by using (2.11) and the surface pressure field.

The prognostic variables are represented by their values at `full-level' pressures . Values for are not explicitly required by the model's vertical finite-difference scheme, which is described below.

The discrete analogue of the surface pressure tendency equation (2.10) is

(2.12)

where

.

(2.13)

From (2.11) we obtain

(2.14)

where is the divergence at level ,

(2.15)

and

.

(2.16)

The discrete analogue of (2.9) is

(2.17)

and from (2.11) we obtain

(2.18)

where is given by (2.14).

Vertical advection of a variable X is now given by

(2.19)

The discrete analogue of the hydrostatic equation (2.6) is

(2.20)

which gives

(2.21)

where is the geopotential at the surface. Full-level values of the geopotential, as required in the momentum equations (2.1) and (2.2), are given by

(2.22)

where and, for ,

,

(2.23)

The remaining part of the pressure gradient terms in (2.1) and (2.2) is given by

(2.24)

with given by (2.23) for all .

Finally, the energy conversion term in the thermodynamic equation (2.3) is discretized as

(2.25)

where , , is defined by (2.23) for , and

(2.26)

The reasons behind the various choices made in this vertical discretization scheme are discussed by Simmons and Burridge (1981); basically the scheme is designed to conserve angular momentum and energy, for frictionless adiabatic flow.

2.2.2 Finite-element vertical discretization

In CY24R3 the vertical discretization has been changed in the operational model from the finite-difference discretization in Lorenz staggering described in the previous subsection to a finite-element discretization using cubic B-splines as basis functions.

For the finite-element (FE) discretization, all variables (even pressure) are kept at the same levels (full levels), i.e. the values of pressure at full levels and not at half levels are required. Also, the values of the derivatives and at full levels are now needed, from which the vertical derivative of pressure can be computed according to . In the semi-Lagrangian version of the evolution equations these are the only vertical derivatives required. They are constant in time and linked to the definition of the vertical coordinate. It is therefore convenient to change the definition of the vertical coordinate and supply and at full levels (instead of and at half levels) such that

(2.27)

and

(2.28)

The two conditions of Eq.(2.28) above ensure that the integral of pressure from the top of the atmosphere to the surface yields exactly the surface pressure . These conditions have to be fulfilled to a good approximation with the numerical integration scheme used. Pressure at any full level can then be obtained by integrating Eq. (2.27) from the top of the atmosphere to the level in question.

The only operation in the vertical which has to be evaluated is the vertical integration. An integral operator based on a finite-element representation will be derived next.

Most of the integrals we have to evaluate are integrals from the top of the atmosphere to the individual model levels and to the surface. We therefore derive an operator in finite-element representation for this type of integral, i.e. an operator which returns the integral from the top of the atmosphere to each of the model levels and to the surface . The vertical integral in the hydrostatic equation (i.e. from the surface upwards) can be constructed by taking the difference of the integral from the top of the atmosphere to the model level in question minus the integral from the top to the surface.

Let and be two sets of linearly independent functions of compact support which can be used as basis functions to expand any function of the vertical coordinate given in the domain [0,1].

The vertical integral

can then be approximated as

(2.29)

where are the coefficients of the expansion of as a linear combination of the basis functions and are the coefficients of the expansion of as a linear combination of the basis functions .

We can then apply the Galerkin procedure to Eq. (2.29) by multiplying both sides of this equation by each function from a complete set of "test functions" and integrating over the vertical domain:

In matrix form this can be expressed as .

Incorporating into the above expression also the transformations from physical space to finite-element space and back, i.e. and , we obtain . Here and denote vectors in physical space composed mainly of the values of and , respectively, at the model levels: , , . The set of values F also includes the value at the surface of the model. Details of the basis functions chosen to implement the scheme as well as how to compute the projection matrices and are given in Untch and Hortal (2002) .

Matrix is the integration operator in finite-element formulation which, applied to a function given at full model levels, yields the integrals of this function from the top of the atmosphere to each individual full model level and to the surface. All the finite sums on the vertical levels in subsection 2.2.1 are replaced by vertical integrals computed by applying the matrix integration operator I. Moreover the quantities are no longer needed as the integration operator gives directly the value of the integral at the model levels (the half levels do not have any meaning in the FE discretization).

2.2.3 Time discretization

To introduce a discretization in time, together with a semi-implicit correction, we define the operators

,

where represents the value of a variable at time , the value at time , and the value at . In preparation for the semi-Lagrangian treatment to be developed in section 3, we also introduce the three-dimensional advection operator

(2.30)

Introducing the semi-implicit correction terms, Eqs. (2.1)-(2.4) become:

(2.31)

(2.32)

(2.33)

(2.34)

where is a parameter of the semi-implicit scheme; the classical scheme (Robert 1969) is recovered with . The semi-implicit correction terms are linearized versions of the pressure gradient terms in (2.1)-(2.2) and the energy conversion term in (2.3). Thus is a reference temperature (here chosen to be independent of vertical level), while and are matrices such that

,

(2.35)

.

(2.36)

where the half-level pressures appearing in (2.35) and (2.35) are reference values obtained from (2.11) by choosing a reference value of , and the coefficients are based on these reference values. The reference values adopted for the semi-implicit scheme are and .

The integrated surface pressure tendency equation (2.14) becomes

(2.37)

where

(2.38)

2.2.4 Horizontal grid

A novel feature of the model is the optional use of a reduced Gaussian grid, as described by Hortal and Simmons (1991). Thus, the number of points on each latitude row is chosen so that the local east-west grid length remains approximately constant, with the restriction that the number should be suitable for the FFT (). After some experimentation, the `fully reduced grid' option of Hortal and Simmons was implemented; all possible wavenumbers (up to the model's truncation limit) are used in the Legendre transforms. A small amount of noise in the immediate vicinity of the poles was removed by increasing the number of grid points in the three most northerly and southerly rows of the grid (from 6, 12 and 18 points in the original design of the T213 grid to 12, 16 and 20 points respectively). Courtier and Naughton (1994) have very recently reconsidered the design of reduced Gaussian grids.

2.2.5 Time-stepping procedure

The time-stepping procedure for the Eulerian - version of the model follows closely that outlined by Temperton (1991) for the shallow-water equations. At the start of a time-step, the model state at time is defined by the values of , , , and on the Gaussian grid. To compute the semi-implicit corrections, the values of divergence , and are also held on the grid, where and

.

(2.39)

The model state at time is defined by the spectral coefficients of , , , and . Legendre transforms followed by Fourier transforms are then used to compute , , , , , , and at time on the model grid. Additional Fourier transforms are used to compute the corresponding values of , . , and . The meridional gradients of and are obtained using the relationships

.

All the information is then available to evaluate the terms at time on the left-hand sides of (2.31)-(2.34) and (2.37), and thus to compute `provisional' tendencies of the model variables. These tendencies (together with values of the variables at are supplied to the physical parametrization routines, which increment the tendencies with their respective contributions. The semi-implicit correction terms evaluated at time-levels () and are then added to the tendencies. Ignoring the horizontal diffusion terms (which are handled later in spectral space), and grouping together the terms which have been computed on the grid, (2.31)-(2.34) and (2.37) can be written in the form

(2.40)

(2.41)

(2.42)

(2.43)

(2.44)

The right-hand sides - are transformed to spectral space via Fourier transforms followed by Gaussian integration. The curl and divergence of (2.40) and (2.41) are then computed in spectral space, leading to

(2.45)

.

(2.46)

Eqs. (2.42), (2.44) and (2.46) can then be combined with the aid of (2.39) to obtain an equation of the form

(2.47)

for each zonal wavenumber and total wavenumber , where the matrix

(2.48)

couples all the values of in a vertical column. Once has been found, the calculation of and can be completed, while and have already been obtained from (2.43) and (2.45).

Finally, a `fractional step' approach is used to implement the horizontal diffusion of vorticity, divergence, temperature and specific humidity. A simple linear diffusion of order is applied along the hybrid coordinate surfaces:

(2.49)

where , or . It is applied in spectral space to the values such that if is the spectral coefficient of prior to diffusion, then the diffused value is given by

(2.50)

A modified form of (2.50) is also used for the temperature , to approximate diffusion on surfaces of constant pressure rather than on the sloping hybrid coordinate surfaces (Simmons, 1987). The operational version of the model uses fourth-order horizontal diffusion

2.2.6 Time filtering

To avoid decoupling of the solutions at odd and even time steps, a Robert filter (Asselin 1972) is applied at each timestep. The time-filtering is defined by

(2.51)

where the subscript denotes a filtered value, and , and represent values at , and , respectively.

Because of the scanning structure of the model, it is convenient to apply the time-filtering in grid-point space, and to split (2.51) into two parts:

(2.52)

(2.53)

The `partially filtered' values computed by (2.52) are stored on a grid-point work file and passed from one time-step to the next. Thus, the information available after the transforms to grid-point space consists of partially filtered values at time together with unfiltered values at time . The filtering of the fields can then be completed via (2.53), which after shifting by one timestep becomes:

.

(2.54)

The computations described in Section 2.2.5 are performed using these fully filtered values at time and the unfiltered values at time . Once (2.54) has been implemented, values of are also available to implement (2.52) for the partially filtered values to be passed on to the next timestep.

2.2.7 Remarks

Ritchie (1988) noted that for a spectral model of the shallow-water equations, the - form and the - form gave identical results (apart from round-off error). In extending this work to a multi-level model, Ritchie (1991) found that this equivalence was not maintained. This was in fact a result of some analytic manipulations in the vertical, used to eliminate between the variables in solving the equations of the semi-implicit scheme, which were not exactly matched by the finite-element vertical discretization of Ritchie's model.

In the case of the model described here, the corresponding elimination between the variables is purely algebraic, and the equivalence between the - form and the - form is maintained apart from one small exception due to the use of the hybrid vertical coordinate. In the - model, the gradients of the geopotential are computed in grid-point space (from the spectrally computed gradients of , and ), while in the - model itself is computed and transformed separately into spectral space, where its Laplacian is added into the divergence equation. Since is not a quadratic function of the model variables there is some aliasing, which is different for the two versions of the model. In practice the differences between the - model and the - model were found to be very small, and in the case of a pure sigma-coordinate the two models would be algebraically equivalent.

The - model is nevertheless considerably more economical than its - counterpart in terms of the number of Legendre transforms required. In addition to the transform of referred to above, four Legendre transforms are saved in the treatment of the wind fields using the procedures described by Temperton (1991) for the shallow-water equations. The number of multi-level Legendre transforms is thereby reduced from 17 to 12 per time-step.

2.2.8 as spectral variable

In preparation for a further reduction in the number of Legendre transforms required by the semi-Lagrangian version of the model, the modified Eulerian version includes an option to keep the virtual temperature , rather than the temperature , as the spectral variable. In the time-stepping procedure, Legendre transforms followed by Fourier transforms are used to compute , and at time on the model grid; the corresponding values of , and are then computed using the corresponding values of , and . The thermodynamic equation (2.3) is then stepped forward in time exactly as before. After the physical parametrization routines, the `provisional' value of is combined with to compute a provisional value of . The semi-implicit correction terms evaluated at time-levels and are then added to the provisional value of , just before the transform back to spectral space.

There are corresponding slight changes in the semi-implicit correction terms. The linearized hydrostatic matrix in (2.31)-(2.32) and (2.39) now operates on rather than on . From the point of view of the semi-implicit scheme, (2.33) has implicitly been replaced by an equation of the form

(2.55)

although as explained above it is not necessary to formulate or compute the missing terms explicitly. Hence, (2.42) is replaced by

(2.56)

and the solution of the semi-implicit equations in spectral space proceeds just as before.

This change of spectral variable results in only insignificant changes to a 10-day model forecast, but permits useful economies in the semi-Lagrangian version to be described in the next chapter.