-
Notifications
You must be signed in to change notification settings - Fork 5
Example 1: Sample Glint2 Setup: ModelE
- Check out a version of Glint2 that has been tested with these examples:
git clone https://github.com/citibob/glint2.git git checkout examples
- Create an appropriate configuration file in configme/ directory for your system. Make sure that '--with-pism' is NOT used, but '--with-galahad' is. Other than that, follow standard build instructions.
Tutorial code is in the file runme.py. It is organized in a series of steps. Each step can be run individually, and output inspected afterwards.
The first step is to generate grid files for the GCM and ice grids. Grid definitions are stored in NetCDF format, and are generated by small C++ programs. Grid definition generation programs are provided for grids used at GISS; others may be easily added. For now, we will use the ones provided.
This step should be run by the command:
cd examples/example1 python runme.py 1
Commands given below are for explanatory purposes, and are included in runme.py step 1.
The program searise_g.cpp may be used to generate the grid used in the SeaRISE Assessment. To generate the standard 5km SeaRISE grid, it may be run as follows:
searise_g 5
The generated file describes all aspects of the grid in question, including the polygonal outline of every grid cell. This is not very interesting for a simple Cartesian grid; however, Glint2 is designed to be compatible with any type of grid, including unstructured meshes. The resulting grid file, with comments included, may be inspected via:
ncdump -h searise_g5.nc
netcdf searise_g5 {
dimensions:
one = 1 ;
grid.vertices.num_realized = 169724 ;
grid.cells.num_realized = 168861 ;
grid.cells.num_realized_plus1 = 168862 ;
grid.cells.num_vertex_refs = 675444 ;
two = 2 ;
three = 3 ;
grid.x_boundaries.length = 302 ;
grid.y_boundaries.length = 562 ;
variables:
int grid.info(one) ;
grid.info:version = 1 ;
grid.info:name = "searise" ;
grid.info:type = "XY" ;
grid.info:type.comment = "\n",
"The overall type of grid, controlling the C++ class\n",
"used to represent the grid. See Grid::Type in\n",
"slib/glint2/Grid.hpp" ;
grid.info:coordinates = "XY" ;
grid.info:coordinates.comment = "\n",
"The coordinate system used to represent grid vertices\n",
"(See Grid::Coordinates in slib/glint2/Grid.hpp. May be\n",
"either XY or LONLAT (longitude comes before latitude). \n",
"Note that this is different from grid.info.type. A\n",
"GENERIC grid, for example, could be expressed in either\n",
"XY or LONLAT coordinates." ;
grid.info:parameterization = "L0" ;
grid.info:parameterization.comment = "\n",
"Indicates how values are interpolated between grid\n",
"points (See Grid::Parameterization in \n",
"slib/glint2/Grid.hpp). Most finite difference models\n",
"will use L0, while finite element models would use L1\n",
"or something else." ;
grid.info:projection = "+proj=stere +lon_0=-39 +lat_0=90 +lat_ts=71.0 +ellps=WGS84" ;
grid.info:projection.comment = "\n",
"If grid.info.coordinates = XY, this indicates the\n",
"projection used to convert local XY coordinates to\n",
"LONLAT coordinates on the surface of the earth. See\n",
"http://trac.osgeo.org/proj/ Proj.4 for format of these\n",
"strings." ;
grid.info:cells.num_full = "168861" ;
grid.info:cells.num_full.comment = "\n",
"The total theoretical number of grid cells (polygons)\n",
"in this grid. Depending on grid.info:parameterization,\n",
"either cells or vertices will correspond to the\n",
"dimensionality of the grid\'s vector space." ;
grid.info:vertices.num_full = 169724 ;
grid.info:vertices.num_full.comment = "\n",
"The total theoretical of vertices (of polygons) on\n",
"this grid." ;
grid.info:grid.vertices.num_realized.comment = "\n",
"The number of \'realized\' cells in this grid. Only the\n",
"outlines of realized cells are computed and stored. \n",
"not all cells need to be realized. For example, a grid\n",
"file representing a GCM grid, in preparation for use\n",
"with ice models, would only need to realize GCM grid\n",
"cells that are close to the relevant ice sheets. In\n",
"this case, all grid cells are realized." ;
int grid.vertices.index(grid.vertices.num_realized) ;
grid.vertices.index:comment = "\n",
"For grids that index on cells (eg, L0): a dense,\n",
"zero-based 1D index used to identify each realized\n",
"cell. This will be used for vectors representing\n",
"fields on the grid." ;
double grid.vertices.xy(grid.vertices.num_realized, two) ;
int grid.cells.index(grid.cells.num_realized) ;
grid.cells.index:comment = "\n",
"For grids that index on vertices (eg, L1): a dense,\n",
"zero-based 1D index used to identify each realized\n",
"vertex. This will be used for vectors representing\n",
"fields on the grid." ;
int grid.cells.ijk(grid.cells.num_realized, three) ;
grid.cells.ijk:comment = "\n",
"OPTIONAL: Up to 3 dimensions can be used to assign a\n",
"\'real-world\' index to each grid cell. If\n",
"grid.info:type = EXCHANGE, then i and j correspond to\n",
"grid.vertices.index of the two overlapping source\n",
"cells." ;
double grid.cells.area(grid.cells.num_realized) ;
int grid.cells.vertex_refs(grid.cells.num_vertex_refs) ;
int grid.cells.vertex_refs_start(grid.cells.num_realized_plus1) ;
double grid.x_boundaries(grid.x_boundaries.length) ;
double grid.y_boundaries(grid.y_boundaries.length) ;
}
The program greenland_2x2_5.cpp may be used to generate areas of the (low-res) GCM grid that are close to the Greenland ice sheet. To generate and inspect, run:
greenland_2x2_5 ncdump -h greenland_2x2_5.nc
An exchange grid is the result of overlapping the polygons of two source grids. These overlap polygons are used in a variety of ways in the regridding algorithms, and can sometimes take a little bit of time to compute. Therefore, the exchange grid is to be pre-computed and viewed:
overlap greenland_2x2_5.nc searise_g5.nc ncdump -h greenland_2x2_5-searise_g5.nc
It's nice to see that the grids are correct. For this, do:
python runme.py 2
Then open the generated file grids.pdf. This plots three grids: the atmosphere grid, the ice grid (at 20km) and the exchange grid. In the exchange grid, notice that only regions where both grids overlap are included.
python runme.py 3
This retrieves the SeaRISE data file that provides high-resolution ice elevations, required for later work.
python runme.py 4 open elev2.png
This plots the elevations read out of the SeaRISE file. The code demonstrates how to plot fields on the ice grid.
python runme.py 5
This creates the file greenland_2x2_5-searise_g5-40-DISMAL.nc. All settings for a particular Glint2 configuration, as well as any ice models it might call through to, are stored in this file. Settings include:
* Details on the ice, GCM and exchange grids. Note that there may be more than one ice grid (eg, one for Greenland and one for Antarctica). * Ice sheet configuration. For each ice sheet, it lists the name of the ice sheet, the ice model used to run it, any parameters to be passed through to the ice model, etc (ice model specific). * Coupling parameters that affect how the GCM interacts with Glint2 and the ice models (GCM-specific).
The Glint2 configuration file is used for all subsequent action: both running with a GCM, and plotting results from that GCM. (Note: In a fully two-way coupled system, the elevation grid changes periodically. Software details of how this will be recorded and later matched with data for plotting purposes have not yet been worked out.)
The GCM may now be run to generate some elevation-classified data files. We won't do that in this tutorial, but instead use results from a sample run of GISS ModelE. See the files:
130516-Regrid_examples_e4f40-hc40-g5_JUL1956.aije4f40-hc40.nc 130516-Regrid_examples_e4f40-hc40-g5_JUL1956.ijhce4f40-hc40.nc
python runme.py 7
With a full Glint2 configuration file and GCM output in hand, we can now test the library's regridding capabilities. At its heart, Glint2 implements five regrid operators, including three forward and two reverse operators. The forward operators are returned as sparse matrices, the user uses them to regrid through matrix multiplication. They are accessed in the Python code via the lines:
RM = mm.hp_to_atm()
RXp = mm.iceinterp_to_atm('greenland', src='ICE')
XM = mm.hp_to_iceinterp('greenland', dest='ICE')
Note that the matrix names correspond to the Glint2 paper.
The two reverse transformations are implemented with quadratic optimization, and therefore cannot be represented by a simple matrix. They are accessed through Glint2 API calls, for example:
impm323 = mm.iceinterp_to_hp(
{'greenland' : impm32},
src='ICE',
qp_algorithm='SINGLE_QP')
impm313 = mm.atm_to_hp(impm31)
NOTES:
# Variables are named with numeric suffixes indicating the grid they are on. 1=GCM grid, 2=ice grid, 3=elevation grid # We know how to plot fields on the GCM or ice grid. But how does one plot a field on the elevation grid? We do so by regridding to the ice grid using "elevation class interpolation," and then plotting the resulting field on the ice grid. In elevation class interpolation, ice grid cells are assigned the value of the nearest elevation point (in elevation space). This should allow the user to see regions of the same color, which represent one elevation point in one GCM grid cell.
CHALLENGE:
# Try regridding the 'prec1' variable from the GCM to ice grid, and plot it.
Now that we've created grids and a Glint2 configuration file, and used them to regrid and plot, we have covered the basics of Glint2. There are more options, explained in the Glint2 paper.