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. 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
- 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
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 -
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 profilesIf 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].