[CF-metadata] WRF staggered grids and vertical coordinates

From: John Caron <caron>
Date: Wed, 29 Oct 2003 17:01:56 -0700

Further thoughts on WRF staggered and vertical coordinates

Most quantities in the file are on a normal grid, say
    float N(t, z, y, x);

This is a projection, and we do have the 2D lat, lon values:
   float lat(t, y, x);
   float lon(t, y, x);

However there are some quantities on a staggered x coordinate (we'll
call x_stag), or on a staggered y coordinate (we'll call y_stag).
Staggrered means that they are offset by half a grid cell, ie these are
the "edges" of the grids:

    float U(z, y, x_stag);
    float V(z, y_stag, x);

The usual thing to do would be to also have lat, lon arrays for these 2
cases:

   float lat2(t, y, x_stag);
   float lon2(t, y, x_stag);

   float lat3(t, y_stag, x);
   float lon3(t, y_stag, x);

(there are none using y-stag, x_stag coordinates, but in principle i
suppose there could be).


Now to add in the vertical coordinates:

The z coordinate is "eta values on half (mass) levels". Here is how the
WRF experts describe how to transform to pressure:

 standard_name = "terrain-following_eta_coordinate"

 Definition: (but these are not coordinate definitions)

   press(n,k,j,i) = P(n,k,j,i) + PB(n,k,j,i)

 where press is the pressure at model gridpoint (k,j,i) and time (n),
 P is the perturbation pressure and PB the base state pressure
 at gridpoint (k,j,i) and time (n). These pressures are defined at
 non-staggered eta levels, which are defined by ZNU.


So we have these fields in the file to satisfy the above definition:
   float P(t, z, y, x);
   float PB(t, z, y, x);

However this will not work for the staggered grids, so we would need:
   float P2(t, z, y, x_stag);
   float PB3(t, z, y, x_stag);

   float P3(t, z, y_stag, x);
   float PB3(t, z, y_stag, x);

The z coordinate also comes in a staggered form, we'll call z_stag, so
we also need

   float P4(t, z_stag, y, x);
   float PB4(t, z_stag, y, x);

To summarize, there are staggered coordinates for x, y, and z. In
principle one could create 8 different grids, in practice it appears
that only 4 are used.
This still seems like a mess, and I doubt we really want all these
variations explicitly added, in the form of 6 extra 4D variables, and 4
extra 3D variables.

In actual practice, the advice is to interpolate back to the normal grid:

 U and V, together with temp and moisture variables are on non-staggered
 vertical grid (we sometimes call it half-eta levels), and pressure too.
 Only vertical motion, W, and geopotential height fields (PH and PHB) are
 output on staggered vertical grid (we sometimes call these levels full-eta
 levels). When we do post-processing using other software, we typically
 average W and PH/PHB to half-eta levels to make all variables colocate
 vertically. This is a simple arithmetic average.

 Then there is the horizontal staggering, where U is not colocated with
 V and all other variables. It may be easier to consider averaging U and V
 to the non-staggered horizontal grid, since when one computes meteorological
 wind fields (whether it is wind speed, or vectors), the wind components
 need to be colocated. On the other hand, if we'd want to do diagnostics, we may
 want to use the U and V on their original staggering. An example maybe the
 calculation of divergence or vorticity. So we may have the need to provide
 height and pressure on U and V grid. To do that, again we can simply average
 the nearby grid values to those at U and V.


So one way out is that variables on staggered coordinates can use the
normal grid coordinates as their coordinate system, by either
interpolating the variable to the normal grid coordinates, or i suppose,
interpolate the normal grid coordinates to the staggered ones. This is
what we might mean by "stagger", that such interpolation is reasonable.
If its not, then you have to specifiy the coordinate systems
independently and in all their gorious detail.

So concretely, if the above proposal was implemented, we are allowed to
use non-staggered lat, lon and Z coordinates for staggered variables.
The client program would be in charge of the interpolation.

   float x(x);
   float x_stag(x_stag);
      x_stag:stagger="x";
   float y(y);
   float y_stag(y_stag);
      y_stag:stagger="y";
   float lat(t, y, x);
   float lon(t, y, x);

   float ZNU(t, z);
      ZNU:long_name = "eta values on half (mass) levels";
      ZNU:positive = "down";
      ZNU:standard_name = "terrain-following_eta_coordinate";
      ZNU:formula_terms = "pressure_perturbation: P pressure_base: PB";

   float ZNW(t, z_stag);
      ZNW:long_name = "eta values on full (w) levels";
      ZNW:positive = "down";
      ZNW:standard_name = "terrain-following_eta_coordinate";
      ZNW:formula_terms = "pressure_perturbation: P pressure_base: PB";
  // P, PB not on z_stag

   float N(t, z, y, x);
      N:coordinates = "lat lon ZNU"; // normal case

    float U(t, z, y, x_stag);
      U:coordinates = "lat lon ZNU"; // lat, lon not on y, x_stag

   float V(t, z, y_stag, x);
      V:coordinates = "lat lon ZNU"; // lat, lon not on y_stag, x

   float W(t, z_stag, y, x);
      V:coordinates = "lat lon ZNW";
Received on Wed Oct 29 2003 - 17:01:56 GMT