# Data Layout for GCRM Output

11/19/08 5:15 PM

This document describes the horizontal and vertical grid used by the
GCRM, how variables are placed on the grid, the information required to
support visualization, and finally a data model for persisted data. The
data model is described in terms of its mapping to the netcdf data
model, though others are possible. The primary purpose is to document
the grid and reach consensus on the data model for the grid. Therefore
an emphasis is placed on the grid definition. A secondary purpose is to
define document the mapping of internal GCRM data structures to actual
output variables. This document is not intended to describe all the
actual data variables though examples of the more complex 2d and 3d
variables are provided.

## GCRM Horizontal Grid

The GCRM uses the geodesic grid. A detailed description of the grid is
beyond the scope of this document but can be found at the CSU web site -
<http://kiwi.atmos.colostate.edu/BUGS/geodesic/>.

The grid will be written using a morton ordering scheme (see
<http://en.wikipedia.org/wiki/Z-order_(curve)>). For our purpose, the
benefits of this scheme are: it supports large contiguous writes and the
data layout is not dependent upon the number of processors used to
create it. An example data set that demonstrates this layout can be
found at <http://climate.pnl.gov/data/>.

The number of grid points in a grid is determined by the equation

N=10x2\*\*(2r) + 2

where r is the recursion number. Some characteristics of the grid based
on the r value is shown in Figure 1. Target r value is r=11 or higher.

![](media/image1.png){width="5.997222222222222in"
height="2.4243055555555557in"}

Figure 1. Table of grid characteristics for various values of r.

Though there are 12 pentagons, all relevant arrays are dimensioned to 6
with one of the values repeated in the case of pentagons. Repeated
values may occur anywhere within the 6 elements. For example a neighbor
list may contain values 2,3,5,5,1900,1901.

## GCRM Vertical Dimension

A vertical column is partitioned using vertical coordinate surfaces.
 The positions of these surfaces are referred to as *interfaces*.    A
model *layer* is sandwiched between two interfaces (see Figure 2).
 Variables will be defined at both *interfaces* and *layers*.  For
example, eta (the horizontal component of vorticity) is defined at
interfaces, and zeta (the vertical component of vorticity) is defined in
the layers.  There will be 1 more interface than layer in a model
column.   For example, if we have a model with 4 layers, the values of
the interfaces might be:

> vertcoor_ifc = {0 km, 10 km, 20 km, 30 km, 40 km}
>
>  

and the values of the layers would then be:

> vertcoor_lyr = {5 km, 15 km, 25 km, 35 km} 

  \
Vertical *interfaces* and *layers* my either increase or decrease with
increasing height. For example, pressure is used as a vertical
coordinate, and it decreases with increasing height. The anticipated
number of *interfaces* is 96 for 4 km model runs.

In the initial versions of the GCRM, the grid metrics (areas, lengths
etc.)  will be the same for all levels. However, in future versions, we
will need to account for increasing metrics (stretching) for increasing
model layers \-- especially for a very deep atmosphere. This document
considers only the case of uniform
metrics.![](media/image2.png){width="5.501388888888889in"
height="3.6756944444444444in"}

Figure 2. Distribution of discrete variables on the grid.

## Variable Location On the Grid

GCRM data can be cell centered, edge centered, or associated with
corners as shown in Figure 2. The list below shows some examples of data
and their grid location, in the initial GCRM model.

> Prognostic variables:

-   Horizontal component of vorticity η ( \* 3 for edges)

-   Potential temperature θ

-   5 species of water

> Diagnostic Variables:

-   Vertical velocity *w*

-   Horizontal velocity V ( \* 3 for edges)

-   Vertical component of vorticity ζ ( \* 2 for corners)

Edges are defined as the midpoint (on the sphere) of the distance
between two cell centers. Corners are defined as the voronoi center of
three cell centers. Additionally, variables are defined at both
*interfaces* and *layers* as described above. The minimum amount of
information necessary to construct all needed locations are:
cell_neighbors, grid_center_lon, grid_center_lat, interfaces, and layers
plus the functions:

midpoint(cell1, cell2)

voronoi_corner(cell1,cell2,cell3)

The number of unique corners and edges on the grid can be determined
from the number of cells as shown in figure 2.

![](media/image3.png){width="5.501388888888889in"
height="3.4298611111111112in"}

### *Cell Centered Data*

Cell centered data is straight-forward to describe. Internally to the
model, 2d and 2d cell centered data will be defined as follows:

> real centervar(i,j,j,b)
>
> real centervar(i,j,k,layer,b)
>
> real centervar(i,j,k+1,b)

where i,j are indices of the horizontal grid within a block, k is the
vertical component, and b is the block. Externally, data will be
persisted with the following structure:

> centervar(cells)
>
> centervar(cells,layers)
>
> centervar(cells,interfaces)

### *Corner Data*

Internally to the model, corner data is associated with the unique
corners of a cell. As previously described, there are two unique corners
per cell. A model data structure may thus be defined as:

> real cornervar(2,i,j,k,b)

where the first index is the unique corners per cell and i,j,k, and b
are as previously described. Some data, such as wind may have multiple
components. For example:

> real cornervar(2,2,i,j,k,b)

stores the north/south component of wind in the 1^st^ value of the first
index whereas the east/west component of wind is stored in the 2^nd^
value of the first index. The other indices are as described above. When
data is persisted, each component of the first index will be written as
two variables in the output data files.

More generally, the form of variables is:

> real cornervar(index2,index1,i,j,k,b)

Index2 may be 2 in the case of the wind example above or three in the
case of vectors. The current plan is to separate all indices at index2
into separate variables but the api will be flexible enough to support
different choices for persistence. Index1 will be 2 for corners or 3 for
edges (discussed below).

The corner data can be persisted one of three ways:

1.  All corners per cells: In this approach, corner data is associated
    with all 6 corners of a given cell. This approach simplifies post
    processing and analysis at the cost of storing 6 values per cell
    instead of two. This data duplication creates unnecessary file size
    bloating and carries a negative impact on bandwidth. Howevver, it
    should be noted that unlike the other options, no additional grid
    description data is required. In c index order, a variable would be
    defined as follows:

> float cornervar(cells, layers, 6)
>
> No additional grid data is required.

2.  Unique corners per cell: This approach follows the design of the
    gcrm model. It minimizes the data munging required to write the
    data. Aditional grid data and a small algorithm are required to map
    the corners to their cells. In c index order, with the I,j values
    transformed into cell indices, we would have the following data
    arrays:

> float cornervar(cell, layer, 2)
>
> int grid_corner_map(cell, 6)
>
> The algorithm required to convert the value in grid_corner_map to
> separate cell and unique corner indices is as follows:
>
> unique_corner = corner % 2,
>
> cell = (corner-unique_corner)/2

3.  Arrays indexed by all unique corners: In this approache, data is
    associated with a "unique corner" index rather than being associated
    with individual cells. A mapping between cell corners and corner
    arry index must be provided. This strategy impacts subsetting in
    that the when requesting a region, multiple different dimensions and
    their coordinates must both be stored and considered in the
    subsetting. In c notation, the following arrays would be required:

> float cornervar(corner, layer)
>
> float grid_corner_lat(corner)
>
> float grid_corner_lon(corner)
>
> int grid_corner_map(cell, 6)
>
> For each cell, the corner map contains the index of the 6 corners into
> the corner dimension. The corner dimension lat/lon are required for
> completeness and easy subsetting.

Option 3 will be adopted because it is conceptually clear and generally
applicable to arbitrary grids.

### Edge Data

Storage of edge data has the same three options as corner data. Using
the preferred option of storing all edge data as a flat array, the
necessary variables are:

> float edgevar(edge, layer)
>
> float grid_edge_lat(edge)
>
> float grid_edge_lon(edge)
>
> int grid_edge_map(cell,6)

## Visualization

Corner and edge data create special difficulties for visualization.
Figure 3 shows considerations for tesselation when visualizating cell,
corner, and edge centered data. For corner variables, the algorithm is
to iterate over each of the cells and retrieve the list of corners. With
this information, the tesselation can be done in a straight forward
manner.

![](media/image4.png){width="5.997222222222222in"
height="3.2472222222222222in"}

Figure 3. Tesselations for cell, corner, and edge centered data.

![](media/image5.png){width="5.983333333333333in"
height="6.111805555555556in"}

Figure4: Tesselation scheme for edge-centered data.

The tesselation scheme for edge-centered data is slightly more
complicated that cell and corner-centered data. Figure4 shows the
desired tesselation. We need to create triangles for the inscribed
hexagon (in blue) and the yellow "star" triangles.

Creating the inscribed hexagon is easy: For a given cell, we find the
position of the center. Then we find the positions of the 6
edge-midpoints. This is sufficient information to create the blue
insribed hexagon. We also need to interpolate the edge data to the
center of the cell.

Creating the yellow "star" triangles is the tricky part. Refering to
Figure4 and taking the light yellow triangle as an example, we know 2
vertices of the triangle (e1 and e2). The question is how do we
determine e?

The algorithm for finding e can be described as follows (ref to Figure
4):

Assumptions:

iteration is performed over cells rather than edges

cell center and corner positions provided

cell neighbor list provided

edge positions provided/generated

mapping from edge to cells may be available

mapping of cell to unique corner and edge indices provided/generated

counter --clockwise ordering for cell neighbors (e.g. h = {h0,h1,h2...})

counter --clockwise ordering for edges (e.g. h = {e0, e1, e2...})

Algorithm for finding e, and hence the triangle (e1,e2, e). We are given
Cell h and edges e1, e2:

1\) Find the relevant neighboring cells h1 & h2\
         a) First we find all cell_neighbors of h\
            cell_neighbors(h) = {h0, h1, h2, h3, h4, h5}\
             - Go through all cell neighbors\
             - For each cell neighbor, go through all edges till you
find e1. Return the cell that contains e1 (Cell h1 in our case).\
\
b) Based on counter-clockwise ordering of cell neighbors, return h2, the
cell that immediately follows h1 in the ordering.

2\) Obtain the set of counter-clockwise edges E {e2, e, \...} in h2.
Return the edge that immediately follows e2 in E. This corresponds to
the desired edge "e" in our case.\
\
We need to take care to create the "star" triangles only once, and not
three times. This can be done by inserting the triangle only when the
current cell index is the minimum of the 3 cells.

Note that if a mapping between an edge and its two cells is also
available, this will speed up the algorithm. This would be defined as
follows:

> int grid_edge_cells(edges,2)

It is also necessary to get the complete description of the cell bounds
in order to draw the polygons. This information is provided by an cell
dimensioned array that provides coordinates for each cell corner.

The vertical dimension requires being able to look up an elevation value
for each level. Currently no scaling or offsetting is planned. A
spherical projection is done on the height and will need to take into
account whether the z direction is moving up or down. This information
must be provided as part of the metadata.

## Subsetting

Current subsetting tools are geared toward lat/lon grids. Because of the
monotonic nature of the coordinates, dimensions can be readily
subsettted by index or value. For the geodesic grid, new software is
necessary. For cell centered data, standard names can be used to find
the lat/lon variables and use those to subset by coordinate value.
However, for this to work for corner and edge data, the data must be
indexed using the same scheme or subsetting must be performed as a
multi-step process. For example, in option 3 in the corner section
above, the cell centered data is subsetted based on cell centers, the
corner data is subsetted based on corner coordinates, and edge data is
subsetted based on edge coordinates. This allows existing tools that
work based on generic dimension processing to work on this data with the
caveat that at the edges of the grid region, cells will be missing some
edge/corner data. However, to ensure that all corner and edge variables
associated with a cell that is within the selection range, the corners
and edges should be computed from the selected cells. More
significantly, all index arrays will need to have the indices updated.

## NetCDF Mapping

Given all the information from the previous sections, we can define a
binding to our netcdf. We rely on CF conventions
([[http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.3/]{.underline}](http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.3/))
where appropriate, and extend them as necessary.

In addition, the standard practice of providing dimension variables will
be followed and can be leveraged by applications. For example, given the
dimension layers, there will also be a variable named layers the
provides the values for each layer.

At this time, it is expected that each 3d variable will be stored in its
own time series file with the full grid definition. The reasons for this
choice at this time include:

-   Time series data for individual variables is likely to make data
    access the most convenient for the majority of operations.

-   Size limits imposed by netcdf3 currently restrict the size of
    variables due to integer offsets

-   Many processing tools do not yet support external grids

Long term, our goal, which is supported by our IO API design, is to
allow the modeler to aggregate variables together into files in any way
they choose which includes control over whether or not a grid variables
are written to each file or to one file. This will require conversion to
netcdf4 or the upgrade of netcdf3 to handle large arrays. The best guess
for the latter is early to mid 2009.

The largest grid size that can be supported with current netcdf
libraries is r=11 which represents a resolution of about 3.5 km.

Initially, variables will be written as floats. Support for double
precision will be added later if necessary. Similarly, indices will be
written as ints initially. Though longs will ultimately be required,
PNetCDF does not currently support them nor does netcdf generally
support dimensions that are longs.

The following netcdf definition is based on horizontal resolution of
r=11. Variables in **bold** are new additions from original sample
files. The exact mapping for corner and cell data is still up for
debate. Two options are under consideration. A cell-centric view is show
in blue. A unique corner/edge view is shown in purple.

The following table summarizes the extra grid information and the
additional size requirements for the alternative grid descriptions. The
approximate extra size for the unique corner/edge view is due to the
duplication of grid corner lat/lon.

  ------------------------ ------------------------ ---------------------
  **[Variable]{.mark}**    **[Formula]{.mark}**     **[Size
                                                    (floats)]{.mark}**

  grid_edge_lat            3 \* (cells-2)           .5 GB

  grid_edge_lon            3 \* (cells-2)           .5 GB

  cell_corners             6 \* cells               1 GB

  cell_edges               6 \* cells               1 GB

  grid_edge_lat2           3 \* (cells-2)           .5 GB

  grid_edge_lon2           3 \* (cells-2)           .5 GB

  grid_corner_lat2         2 \* (cells-2)           .3 GB

  grid_corner_lon2         2 \* (cells-2)           .3 GB

  corner_cells             3\* 2 \* (cells-2)       1.0 GB

  edge_cells               2\* 3 \* (cells-2)       1.0 GB

  *Total*                                           3.0 GB /3.6 GB
  ------------------------ ------------------------ ---------------------

Dimensions:

> **interfaces**: 96\
> **layers**: 95 \# interfaces -1\
> cells: 41943042
>
> **corners**: 83886080 \# total unique corners in the grid; 2 \* (cells
> -- 2)
>
> **edges**: 125829120 \# total unique edges in the grid; 3 \* (cells -
> 2)
>
> cellcorners: 6 \# max corners per cell
>
> celledges: 6 \# max edges per cell
>
> cellneighbors: 6 \# max neighbors per cell

Variables:

float **interfaces**(interfaces)\
interfaces:units = "Pa";

interfaces:**positive** = "down";

interfaces::**axis** = "Z";

float **layers**(layers)\
layers:units = "Pa";

> layers:**positive** = "down";
>
> layers::**axis** = "Z";
>
> float grid_center_lat(cells)
>
> grid_center_lat:long_name = \"Latitude of cell center\" ;
>
> grid_center_lat:units = \"radians\" ;
>
> grid_center_lat:standard_name = \"latitude\" ;
>
> grid_center_lat:bounds = \"grid_corner_lat\" ;
>
> float grid_corner_lat(cells,cellcorners)
>
> float grid_center_lon(cells)
>
> grid_center_lon:long_name = \"Longitude of cell center\" ;
>
> grid_center_lon:units = \"radians\" ;
>
> grid_center_lon:standard_name = \"longitude\" ;
>
> grid_center_lon:bounds = \"grid_corner_lon\" ;
>
> float grid_corner_lon(cells, cellcorners)
>
> float area(cells)
>
> grid_area:long_name = \"cell area\" ;
>
> grid_area:standard_name = \"area\" ;
>
> grid_area:units = \"square radians\" ;
>
> int cell_neighbors(cells, cellneighbors)
>
> cell_neighbors:long_name = \"Index of cell neighbors\" ;
>
> cell_neighbors:**traverse** = "counter-clockwise" ;
>
> float **grid_edge_lat***(*cells, celledges)
>
> grid_edge_lat:long_name = \"Mid-point latitude of edge \" ;
>
> grid_edge_lat:standard_name = "latitude" ;
>
> grid_edge_lat:units = \"radians\" ;
>
> float **grid_edge_lon**(cells, celledges)
>
> grid_edge_lon:long_name = \"Mid-point longitude of edge\" ;
>
> grid_edge_lon:standard_name = "latitude" ;
>
> grid_edge_lon:units = \"radians\" ;
>
> int **cell_corners**(cells, cellcorners)
>
> cell_corners:**traverse** = "counter-clockwise" ;
>
> cell_corners:long_name = "mapping of unique corner indices for each
> cell";
>
> int **cell_edges**(cells, celledges)
>
> cell_edges::**traverse** = "counterclockwise" ;
>
> cell_edges:long_name = "mapping of unique edge indices for each cell";
>
> float **grid_corner_lat2**(corners):
>
> grid_corner_lat2:long_name = "TBD";
>
> float **grid_corner_lon2**(corners):
>
> grid_corner_lon2:long_name = ""TBD";
>
> *float* **grid_edge_lat2(**edges):
>
> grid_edge_lat2:long_name = "TBD";
>
> float **grid_edge_lon2**(edges):
>
> grid_edge_lon2:long_name = ""TBD";
>
> int **corner_cells**(corners, 3):
>
> corner_cells:long_name = "list of cells per unique corner";
>
> int **edge_cells**(edges, 2):
>
> edge_cells:long_name = "list of cells per unique edge";
>
> int ***grid_global_index***(cells): \# To reconstruct full globe from
> regional files IFF global files not used
>
> grid_global_index:long_name = "cells index into global cell array"
>
> Some example variables on this grid
>
> float centervar1(time, cells, layers): \# a cell centered variable
> residing on a layer
>
> centervar1:coordinates = "grid_center_lat grid_center_lon";
>
> float centervar2(time, cells, interfaces): \# a cell centered variable
> residing on interface
>
> centervar2:coordinates = "grid_center_lat grid_center_lon";
>
> float cornervar(time, corners, layers)
>
> cornervar:coordinates = "grid_corner_lat grid_corner_lon";
>
> float edgevar(time, edges, layers)
>
> edgevar:coordinates = "grid_edge_lat grid_edge_lon";
>
> \# global attributes:\
>          :Conventions = \"CF-1.3\" ;\
>         :institution = \"\<specifies where the original data was
> produced\>\" ;\
>         :title = \"\<a succint description of what is in the
> dataset\>\" ;\
>         :history = \"\<audit trail\>\" ;\
>         :source = \"\<method of production of the original data\>\" ;\
>         :references = \"\<published or web-based references that
> describe the data or methods to produce it\>\" ;\
>        :grid = \"geodesic\" ;
>
> :**grid_cell_index_origin**: = \"0\" ; \# Identifies index of first
> cell (i.e. origin)