Dr. Gavin A. Schmidt

Using tracer moments within the GISS GCM

Version 1.2 Dec 1999
Gavin Schmidt and Gary Russell


The GISS GCM is unique in using prognostic first and second order moments as well as the mean for all the tracers (including water vapour and potential enthalpy) within each grid box. These were introduced to improve the advection of tracers, and have not as yet truly been used to their maximum potential. In theory, the quadratic (linear) moments provide sub-grid scale horizontal resolution three (twice) times more detailed than the coarser resolution of the U and V fields. That is to say that using the linear moments is (almost) equivalent to four boxes within each grid cell, while using the quadratic moments gives 9 boxes per grid cell. (Strictly speaking, the moments are not exactly equivalent to the higher resolution, information is lost going from 9 horizontal boxes to the 6 horizontal moments. However, this is a minor point). This takes 4x5 resolution down to 1.3x1.7 resolution for the tracers. In the vertical it can also double or triple the resolution.

This paper is intended to be a simple introduction to using the moments more effectively in the GCM (i.e. not just in the dynamics). We will use a generic concentration tracer TR(I,J,L) to describe how the moments should be changed based on common processes in the GCM. We should note that although water vapour Q and potential enthalpy T are saved as concentrations, the tracers (such as Radon, O-18 and the like) are generally saved as masses. This is simply a unit issue and does not effect any of the following advice. In the GCM, for quantities in mass units, the moments are generally denoted with an M in the variable name (such as QXM).

Figure 1
Figure 1. The quadratic and linear profiles within a gird box.

For each tracer there is the mean value TR0; the three first order (linear) moments, TRX, TRY, and TRZ; and the 6 second order (quadratic) moments TRXX, TRXY, TRYY, TRYZ, TRZZ, and TRZX (some older code uses the variable name TRXZ instead). They are scaled so that they are generally the same order of magnitude as the mean. A good way to visualize them in the linear case (for a one dimensional tracer) is to think of the value of the tracer at three points in the box (Figure 1). The left hand most value is TR-TRX, the mid point is TR and the right hand most value, TR+TRX. As a function of the distance x from the mid point (scaled so that the box is of length 2x with midpoint x=0), the tracer has the value

TR(x)= TR0 + x TRX, -1 < x < 1

Including the second order moments is also straightforward,

TR(x)= TR0 + x TRX + (3/2 x2-1/2)TRXX, -1 < x < 1

Note that the mid point value is now TR0 - TRXX/2. These equations can of course be generalized in the three directions (making sure to include the cross terms for the quadratic case).

TR(x,y,z)= TR0 + x TRX + y TRY + z TRZ + xy TRXY + yz TRYZ + zx TRZX + (3/2 x2-1/2) TRXX + (3/2 y2-1/2) TRYY + (3/2 z2-1/2) TRZZ
with -1 < x,y,z < 1. For mass tracers, each equation should be divided by the air mass to get concentration.

There are two stages to using the moments more effectively. Firstly, by making sure that all operations on the mean fields are carried over consistently to the moments, and secondly, actually using the sub-grid scale resolution to perform more detailed physics within each grid box. The following sections outline the details of how this should be done.

Rules for manipulating moments consistently

In this section we list the most common operations to the tracers and how the gradients should be manipulated consistently. The overriding principle is that operators should be applied to the whole of the profile (equation 3). Since it is linear, the operators can be applied to each of the moments independently. The secondary principle is that the profile should approximate the least squares quadratic fit to the tracer distribution, hence any change that is not uniform within a grid box can be incorporated by adjusting the profile accordingly.

    Complete mixing of two boxes:
    All gradients should also be mixed, except those in the direction of the mixing which should be set to zero. For instance, during dry convective mixing set TRZ=TRZZ=TRYZ=TRZX=0 in each box and each other moment to its mass-weighted average over both boxes.

    Removing a fraction of tracer uniformly over a grid box:
    All moments should be multiplied by (1-a), where a is the fraction of the tracer that is removed. For instance, this should be done for radioactive decay, or condensation of water vapour.

    Adding tracer uniformly over a grid box:
    This is somewhat more subjective, but for stability the gradients should not be changed at all. Otherwise, gradients can become unreasonably large and destabilize the model. The only time this possibly should cause a change for the gradients is when the process that is increasing the tracer is a function of the tracer concentration itself (such as an auto-catalytic reaction?). This asymmetry between adding and removing tracer uniformly is slightly diffusive. In the Michael Prather's chemistry models, the changes in moments of each species depend on the moments of the other species involved in its production.

    Defining a vertically moving plume:
    All horizontal moments should be preserved for the plume. All vertical moments (including TRYZ and TRZX) and should be set to zero within the plume.

    Removing a vertically moving plume:
    If the plume originates uniformly from the box, all moments should be reduced by (1-a), where a is the fraction of the tracer that is removed. If the plume originates somewhere specific then the vertical moments should be re-calculated according to the least square fit (see /HOSTS/Isis/u/cmglr/DOC/PLUMES.060777 for exact details on how this is done for the linear scheme).

    Adding a vertically moving plume:
    If it is inserted uniformly, all moments (in mass units) within the plume should be added to the moments (again, in mass units) of the receiving grid box. Again, if the plume is inserted in a specific point the vertical moments should be re-calculated according to the least square fit (see /HOSTS/Isis/u/cmglr/DOC/PLUMES.060777).

    Diffusing tracer within a grid box:
    For diffusion in one dimension (for instance vertically), the same operator that applies to the mean should also apply to the moments in the unrelated directions (i.e QX, QY, QXX, QYY, QXY). The boundary conditions on these gradients should be no-flux at the surface (unless of course some sub-grid scale processes have been calculated using moments as well). For the moments involving Z, (i.e. QZ, QZZ, QZX, and QYZ) the equations for modifying them come from looking at weighted integrals of the exact diffusion equation within the box and depend on the values of the diffusion and the fluxes at the top and bottom of the box. See Appendix A for the relevant equations.

    Using a separately calculated sub-gridscale profile:
    The PBL code and the radiation schemes (?) use vertical profiles that are calculated separately from the moments. These should be used to update the vertical moments so that the rest of the code can use them. All that is needed is the least squares fit of equation 2 in the vertical, to the internally calculated profile. Possibly, the vertical profile should be used as to adjust the initial values for the internal profiles (since other processes may have changed the vertical moments since the last time the internal profiles were calculated). See Appendix B for how the equations for this.

    Adding a point source (usually at the surface):
    When tracer sources occur at the surface without any associated mass flux, a delta function is the best approximation of the change. A least squares fit to the new profile would then give the new gradients. However, (depending on the amount of tracer present in the box already relative to the addition), such a scheme would violate the positive definite criteria. Hence the changes to the gradients in this case must be modified. We choose to modify them based on two criteria: firstly, that at no place does the implied tracer concentration decrease, and secondly, that the changes are monotonic to the opposite edge. With these constraints, for a point source of amount DTR at the bottom surface (z=-1):
    TR0 = TR0 + DTR ; TRZ = TRZ - 3/2 DTR ; TRZZ = TRZZ + 1/2 DTR
There are a couple of further points that need to be addressed. If a tracer is positive definite (like Q), then the sub-gridscale profile also should be positive definite. However, this is not absolutely guaranteed by the dynamics. Applying the following limits (taken from Prather, 1986) will ensure this is the case, however, it does not ensure that the profile is smooth across grid-boxes. If these limits need to be applied, it is possible that the profile has already become very noisy and the numerics already compromised.

TRX =min(1.5 TR0,max(-1.5 TR0, TRX))
TRXX=min(2 TR0-| TRX|/3,max(| TRX|- TR0, TRXX))
TRXY=min( TR0,max(- TR0, TRXY))
TRZX=min( TR0,max(- TR0, TRXZ))

In order to check the fully resolved profile, the program ra:/u/exec/qusplt will take as input lines of the form TR0 TRX TRXX and plot the resulting profile. It will work with either piped input, or read from a file. The output is displayed on the screen and saved in x.ps (using aplotX). See file ra:/u/exec/qusplt.doc for more info.

Incorporating non-uniform sub-gridscale physics

There are two ways to incorporate processes that are not uniform over the grid box. Firstly, we can breakdown the grid box into component sub-boxes, or we can directly manipulate the tracer moments. It is more straightforward to split the box into sub-boxes, but for column physics and the like, this leads to a large increase in the time required for the routines. Changing the gradients directly (as is done in DYNAM) is more efficient, but also more complicated.

A note of caution: Although the tracer concentration profile can be extended to the edges of any grid box, the errors in using the profile increase with x, and the linear or quadratic profile is not guaranteed to be continuous across the edge to the next box. Hence it is more conservative to use values of x less than unit magnitude. The C-grid model code generally uses x=+/-2/sqrt(12)=0.57735 since those are the points at which the linear least squares fit crosses the quadratic profile. These points also have nice properties for evaluating integrals of the tracer over the box. However, for other cases it may be more straightforward to evenly split a box in two (x=+/-1/2) or three (x=-2/3,0,2/3).

If we decide to split the box into sub-boxes, define the horizontal sub-boxes using IQ=1,3 and JQ=1,3 and set T(IQ,JQ)=TR(2*(IQ-2)/3, 2*(JQ-2)/3)) using equation 3 in two dimensions. The physics can now be applied for each (IQ,JQ) pair. Recalculating the moments is straightforward although there may be round off errors to beware of (see /HOSTS/Isis/u/gavin/gissgcm/doc/sigfig.doc for more details). In the vertical a similar process can be used. For example, testing for stability within a grid box can be done twice (between z=-2/3 and z=0 and again between z=0 and z=2/3). This could allow for more low level convective cloud for instance.

Examples of manipulating the gradients directly (i.e. the Prather chemical schemes) depend heavily on the process involved and no examples will be given here.


Even if much of the physics is calculated at the original resolution, there is no excuse for not manipulating the gradients in a consistent manner. Eventually, everywhere the tracers are modified some thought must be given to the effect of the moments. Examples of how this works can be found in the condensation and moist convection routines ( CB263.S), in SURFCE ( SB354BM15), DRYCNV ( PB354AM15), and in FILTER ( MB354AM15).

It is possible that the noise introduced by the inconsistency in treating the gradients is beneficial to the model, increasing diffusion or triggering more convection. However, it does not seem a good idea to rely on these mechanisms. If more diffusion is wanted, then it should be explicitly entered in a controlled, and hopefully more physical, manner.

Appendix A: Modifying the moments during diffusion

We will consider vertical diffusion with a non-dimensional z (-1 < z < 1), and with K in units of per second. We define K'=(K1+K-1)/2, dK=(K1-K-1) and similarly for the fluxes F = K dQ/dz (however the gradient is calculated). Then if dt is the timestep, the equations for the moments are,

QZX+=QZX- - 3 dt K' QZX
QYZ+=QYZ- - 3 dt K' QYZ
QZ+=QZ- + dt (3F' - 3 K' QZ - 1.5 dK QZZ)
QZZ+=QZZ-+dt (2.5 dF - 2.5 dK QZ - 15 K' QZZ)

Note that if there are asymmetries in the fluxes or diffusion coefficient, the first and second moments are entangled. If we assume an implicit calculation, the new moments are

QZX+= QZX- / (1 + 3 dt K'), (and similarly for QYZ) and

QZ+ = [(1+15 dt K') ( QZ-+3 dt F') - 1.5 dt dK ( QZZ-+2.5 dt dF)]/ D
QZZ+= [(1+3 dt K') ( QZZ-+2.5 dt dF) - 2.5 dt dK ( QZ-+3 dt F')]/ D

where D=(1 + 3 dt K')(1 + 15 dt K')- 3.75 dt 2 dK2 , (D > 0 if dt >= 0 and K >= 0). This simplifies considerably in the case of constant diffusivity (dK = 0).

In dimensional terms, every K should be replaced with 4 K/dz2 , every F with 2F/ dz where dz is the height of the box in whatever unit is appropriate (either m (height coordinates) for K in m2/s, or kg/m2 (mass coordinates) for K in units kg2/s m4).
(Thanks to Max Kelley for helping correct this section).

Appendix B: Quadratic least squares fit of an internal profile

This should be useful for using changes in internal profiles (such as used in the PBL code) into changes to the vertical moments. For a uniformly spaced internal vertical profile tr(zi), define departures from the mean tr'(zi). If the z coordinate is scaled such that -1 < zi <1 over the grid box, the quadratic least squares fit is given by,

TRZ = (T1S22 - T2 S12)/(S11 S22 - S122)
TRZZ=(T1 S12 - T2 S11)/(S11 S22 - S122)
T1=sumi zi tr'(zi), T2=1/2 sumi (1 - 3 zi2) tr'(zi)
S11=sumi zi2, S12=1/2 sumi zi (1 - 3 zi2),
S22=1/4 sumi (1 - 3 zi2)

If the zi are not uniform then all the sums will include extra weighting terms.

Simple case for half-box profiles

If the internal vertical profile is valid from the mid point down (i.e. over only half the box), and an even weighting between the upper and lower half of the boxes is used, the equations are very simple:

TRZ+ = 3/2 ( tz + 1/2 TR0- + 1/3 TRZ- + 1/8 TRZZ-)
TRZZ+ = 5/2 ( tzz + 1/8 TRZ- + 1/5 TRZZ-)

where tz,tzz are the integrals of z*t(z) and 1/2 (3z2-1)*t(zz) over the interval [-1,0].