wiki:ReaderLogic
Last modified 5 years ago Last modified on 06/06/12 23:22:09

Proposed re-write of MOAB NetCDF reader to read many types of CESM and other grids

NOTE: This is pseudocode mixed with C++.

Keep load_file the same up to the call to read_header:

    // Read the header (num dimensions, dimensions, num variables, global attribs)
  rval = read_header();
  ERRORR(rval, " ");

After that, add code to see if CF convention is being followed.

    Inquire for a global attribute called "conventions" or "Conventions"
        If NOT ( "conventions" global attribute exists AND contains the string "CF" )    
                fatal error: "File not following known convention"
        endif

Next change is to substitue the call to init_ijkt_vals with new code.

Setting up the grid is different for almost every grid in CESM. So need one init function per grid. Replace this:

    // Get bounds on ijk space
  rval = init_ijkt_vals(opts, scdi);

  if (MB_SUCCESS != rval) {
      // try reading ucd format (HOMME) data
    rval = init_ucd_mesh(opts);
    if (MB_SUCCESS == rval) {
      ucdMesh = true;
    }
  }
  ERRORR(rval, "Trouble initializing ijk values.");

with

if global attribute "source" has the string CAM then
   if dimension names "lon" AND "lat" AND "slon" AND "slat" exist then   !  Its the FV grid
      rval = init_FVCDscd_vals(opts, scdi)
   else if global attribute "np" exists then  ! its the HOMME grid
      rval = init_HOMMEucd_vals
      if (MB_SUCCESS == rval) {
        ucdMesh = true;
      }
   else if dimension names  "lon" and "lat' exist then ! its the Eulerian Spectral grid
      rval = init_EulSpcscd_vals(opts, scdi)
   else
      error "unknown CAM grid"
else
      error "unknown grid"   (will fill this in later for POP, CICE and CLM)
endif

ERRORR(rval, "Trouble initializing mesh values.");

init_FVCDscd_vals is based on the current version of init_ijkt_vals but with less guessing because we know what the grid is.

ErrorCode ReadNC::init_FVCDscd_vals(const FileOptions &opts, ScdInterface *scdi)
{
    // construct grid for Finite Volume CD-grid

    // get  i/j/k dimensions
  gDims[0] = gDims[3] = -1;
  std::vector<std::string>::iterator vit;
  unsigned int idx;
  if ((vit = std::find(dimNames.begin(), dimNames.end(), "slon")) != dimNames.end())
  else ERRORR(MB_FAILURE, "Couldn't find slon variable.");
  iDim = idx;
  gDims[3] = dimVals[idx]-1;
  gDims[0] = 0;
  iName = dimNames[idx];

   if ((vit = std::find(dimNames.begin(), dimNames.end(), "slat")) != dimNames.end())
  else ERRORR(MB_FAILURE, "Couldn't find slat variable.");
  jDim = idx;
  gDims[4] = dimVals[idx]-1+2;    // add 2 for the pole points
  gDims[1] = 0;
  jName = dimNames[idx];

  gDims[2] = gDims[5] = -1;     //  ilev actually has the number of interfaces which is what we want.
  if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end()) {
    idx = vit-dimNames.begin();
    gDims[5] = dimVals[idx]-1, gDims[2] = 0, kName = std::string("ilev");
    kDim = idx;
  }

// we know time variable is called "time" in a CESM CF file.
  if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
    idx = vit-dimNames.begin();
  else ERRORR(MB_FAILURE, "Couldn't find time variable.");
  tDim = idx;
  tMax = dimVals[idx]-1;
  tMin = 0;
  tName = dimNames[idx];

NEW: To determine what type of variable it is, this logic can be used in both of the structured atmosphere readers.

  set all variables to be on a set by default.
  if variable's dimensions contain at least ("lat" AND "lon")
    it should be on a quad
  if variable's dimensions contain at least ("slat" AND "lon")
    it should be on a north/south edge
  if variable's dimensions contain at least  ("lat" AND "slon")
     it should be on an east/west edge

code from "initialize parameter bounds" to the end of the function is the same EXCEPT the reading of the j coordinate variableswhich should be modified as follows:

  if (lDims[1] != -1) {
    if ((vmit = varInfo.find(jName)) != varInfo.end() && (*vmit).second.varDims.size() == 1) {

      if(lDims[1] == gDims[1]) then   // If this is the first row
         // need to read one less then available and read it into a dummy var
         create dummy var of size (lDims[4]-lDims[1])
         rval = read_coordinate(jName.c_str(), lDims[1] , lDims[4]-1  dummyVar);
          ERRORR(rval, "Trouble reading y variable.");
         // copy the correct piece
         jlVals[0]=-90.0;
         for(i=1; i< lDims[4]+1; i++)
              jlVals[i]=dummyVar[i-1];
      

      else if (lDims[4] == gDims[4] ) then     // or if its the last row
         create dummy var of size (lDims[4]-lDims[1])
         rval = read_coordinate(jName.c_str(), lDims[1]-1 , lDims[4]  dummyVar);
          ERRORR(rval, "Trouble reading y variable.");
         // copy the correct piece
         for(i=0; i< lDims[4]; i++)
              jlVals[i]=dummyVar[i];

         jlVals[i] = 90.0   // usin value of i after loop exits.

      else  // its in the middle
        rval = read_coordinate(jName.c_str(), lDims[1]-1 , lDims[4]-1  jlVals);
        ERRORR(rval, "Trouble reading y variable.");
      endif
    }
    else {
      ERRORR(MB_FAILURE, "Couldn't find y coordinate.");
    }
  }

The call to check_vert_hexes, as its currently written, should only be done if ucdMesh is false. Now that the vertices are created correctly (as above) the hexes should be correct.

?? The call to create_vert_hexes should be modified to also create quads (cell centers) ??

?? where do we set the coordinates of the quads ??

The call to read_variables, this line should be modified.

 rval = mbImpl->get_entities_by_dimension(file_set, 0, verts);

The "0" should be a "1" or "2" to get the quads or edges (??? not sure which or if we need both ???)

read_variable_allocate, called by read_variables, needs to be modified. This line gets a pointer to the internal representation of a tag on vertices. The pointer is later used in read_variables.

rval = mbImpl->tag_iterate(vdatas[i].varTags[t], verts.begin(), verts.end(), count, data);

The "verts.begin(),verts.end()" needs to be changed to "quads.begin(),quads.end()" to get a pointer to the quads.

Need to be aware of the order of the quads. It has to match the order in the NetCDF file IF we read directly in to the quad pointer. If the order is wrong, read in to a temporary buffer and then load the MOAB data correctly.

I think the code called "Create partition set" can stay the same.

In create_tags, the block of code that checks tags names "lon, "lat" should be pulled out and put in to a grid-specific version.

NEW: create_tags can be moved to earlier in load_file