src/io/NCWriteHelper.cpp

/*
 * NCWriteHelper.cpp
 *
 *  Created on: Mar 28, 2014
 *      Author: iulian
 */

#include "NCWriteHelper.hpp"
#include "NCWriteEuler.hpp"
#include "NCWriteFV.hpp"
#include "NCWriteHOMME.hpp"
#include "NCWriteMPAS.hpp"
#include "NCWriteGCRM.hpp"

#include "moab/WriteUtilIface.hpp"
#include "MBTagConventions.hpp"

#include <sstream>

#ifdef WIN32
#ifdef size_t
#undef size_t
#endif
#endif

namespace moab {

//! Get appropriate helper instance for WriteNC class; based on some info in the file set
NCWriteHelper* NCWriteHelper::get_nc_helper(WriteNC* writeNC, int fileId, const FileOptions& opts, EntityHandle fileSet)
{
  std::string& grid_type = writeNC->grid_type;
  if (grid_type == "CAM_EUL")
    return new (std::nothrow) NCWriteEuler(writeNC, fileId, opts, fileSet);
  else if (grid_type == "CAM_FV")
    return new (std::nothrow) NCWriteFV(writeNC, fileId, opts, fileSet);
  else if (grid_type == "CAM_SE")
    return new (std::nothrow) NCWriteHOMME(writeNC, fileId, opts, fileSet);
  else if (grid_type == "MPAS")
    return new (std::nothrow) NCWriteMPAS(writeNC, fileId, opts, fileSet);
  else if (grid_type == "GCRM")
    return new (std::nothrow) NCWriteGCRM(writeNC, fileId, opts, fileSet);

  // Unknown NetCDF grid
  return NULL;
}

ErrorCode NCWriteHelper::collect_variable_data(std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
{
  Interface*& mbImpl = _writeNC->mbImpl;
  std::vector<std::string>& dimNames = _writeNC->dimNames;
  std::vector<int>& dimLens = _writeNC->dimLens;
  std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
  std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
  std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
  DebugOutput& dbgOut = _writeNC->dbgOut;

  ErrorCode rval;

  usedCoordinates.clear();

  if (tstep_nums.empty() && nTimeSteps > 0) {
    // No timesteps input, get them all
    for (int i = 0; i < nTimeSteps; i++)
      tstep_nums.push_back(i);
  }

  for (size_t i = 0; i < var_names.size(); i++) {
    std::string varname = var_names[i];
    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(varname);
    if (vit == varInfo.end())
      MB_SET_ERR(MB_FAILURE, "Can't find variable " << varname);

    WriteNC::VarData& currentVarData = vit->second;

    dbgOut.tprintf(2, "    for variable %s varDims.size %d \n", varname.c_str(), (int)currentVarData.varDims.size());
    for (size_t j = 0; j < currentVarData.varDims.size(); j++) {
      std::string dimName = dimNames[currentVarData.varDims[j]];
      vit = varInfo.find(dimName);
      if (vit == varInfo.end())
        MB_SET_ERR(MB_FAILURE, "Can't find coordinate variable " << dimName);

      usedCoordinates.insert(dimName); // Collect those used, we will need to write them to the file
      dbgOut.tprintf(2, "    for variable %s need dimension %s with length %d\n", varname.c_str(), dimName.c_str(), dimLens[currentVarData.varDims[j]]);
    }

    // Process coordinate variables later
    if (usedCoordinates.find(varname) != usedCoordinates.end())
      continue;

    // Default has_tsteps is false
    if (std::find(currentVarData.varDims.begin(), currentVarData.varDims.end(), tDim) != currentVarData.varDims.end())
      currentVarData.has_tsteps = true;

    // Default numLev is 0
    if ((std::find(currentVarData.varDims.begin(), currentVarData.varDims.end(), levDim) != currentVarData.varDims.end()))
      currentVarData.numLev = nLevels;

    // Process set variables
    if (WriteNC::ENTLOCSET == currentVarData.entLoc) {
      if (currentVarData.has_tsteps) {
        // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
        // TBD
        MB_SET_ERR(MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet");
      }
      else {
        // Get the tag with varname
        Tag tag = 0;
        rval = mbImpl->tag_get_handle(varname.c_str(), tag);MB_CHK_SET_ERR(rval, "Can't find tag " << varname);
        currentVarData.varTags.push_back(tag); // Really, only one for these
        const void* data;
        int size;
        rval = mbImpl->tag_get_by_ptr(tag, &_fileSet, 1, &data, &size);MB_CHK_SET_ERR(rval, "Can't get data of tag " << varname);

        // Find the type of tag, and use it
        DataType type;
        rval = mbImpl->tag_get_data_type(tag, type);MB_CHK_SET_ERR(rval, "Can't get data type of tag " << varname);

        currentVarData.varDataType = NC_DOUBLE;
        if (MB_TYPE_INTEGER == type)
          currentVarData.varDataType = NC_INT;

        assert(0 == currentVarData.memoryHogs.size()); // Nothing so far
        currentVarData.memoryHogs.push_back((void*)data);

        if (currentVarData.varDims.empty()) {
          // Scalar variable
          currentVarData.writeStarts.push_back(0);
          currentVarData.writeCounts.push_back(1);
        }
        else {
          for (size_t j = 0; j < currentVarData.varDims.size(); j++) {
            currentVarData.writeStarts.push_back(0);
            currentVarData.writeCounts.push_back(dimLens[currentVarData.varDims[j]]);
          }
        }

        // Get variable size
        currentVarData.sz = 1;
        for (std::size_t idx = 0; idx != currentVarData.writeCounts.size(); idx++)
          currentVarData.sz *= currentVarData.writeCounts[idx];
      }
    } // if (WriteNC::ENTLOCSET == currentVarData.entLoc)
    // Process non-set variables
    else {
      Tag indexedTag = 0;

      if (currentVarData.has_tsteps) {
        for (unsigned int t = 0; t < tstep_nums.size(); t++) {
          std::stringstream ssTagNameWithIndex;
          ssTagNameWithIndex << varname << tstep_nums[t];
          rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);MB_CHK_SET_ERR(rval, "Can't find tag " << ssTagNameWithIndex.str());
          dbgOut.tprintf(2, "    found indexed tag %d with name %s\n", tstep_nums[t], ssTagNameWithIndex.str().c_str());
          currentVarData.varTags.push_back(indexedTag);
        }
      }
      else {
        // This should be a user-created non-set variable without timesteps
        // Treat it like having one, 0th, timestep
        std::stringstream ssTagNameWithIndex;
        ssTagNameWithIndex << varname << 0;
        rval = mbImpl->tag_get_handle(ssTagNameWithIndex.str().c_str(), indexedTag);MB_CHK_SET_ERR(rval, "Can't find tag " << ssTagNameWithIndex.str() << " for a user-created variable");
        dbgOut.tprintf(2, "    found indexed tag 0 with name %s\n", ssTagNameWithIndex.str().c_str());
        currentVarData.varTags.push_back(indexedTag);
      }

      // The type of the tag is fixed though
      DataType type;
      rval = mbImpl->tag_get_data_type(indexedTag, type);MB_CHK_SET_ERR(rval, "Can't get data type of tag " << varname);

      currentVarData.varDataType = NC_DOUBLE;
      if (MB_TYPE_INTEGER == type)
        currentVarData.varDataType = NC_INT;
    }
  } // for (size_t i = 0; i < var_names.size(); i++)

  // Process coordinate variables here
  // Check that for used coordinates we have found the tags
  for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
      setIt != usedCoordinates.end(); ++setIt) {
    const std::string& coordName = *setIt;

    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(coordName);
    if (vit == varInfo.end())
      MB_SET_ERR(MB_FAILURE, "Can't find coordinate variable " << coordName);

    WriteNC::VarData& varCoordData = vit->second;
    Tag coordTag = 0;
    rval = mbImpl->tag_get_handle(coordName.c_str(), coordTag);MB_CHK_SET_ERR(rval, "Can't find tag " << coordName);
    varCoordData.varTags.push_back(coordTag); // Really, only one for these

    const void* data;
    int sizeCoordinate;
    rval = mbImpl->tag_get_by_ptr(coordTag, &_fileSet, 1, &data, &sizeCoordinate);MB_CHK_SET_ERR(rval, "Can't get coordinate values of " << coordName);
    dbgOut.tprintf(2, "    found coordinate tag with name %s and length %d\n", coordName.c_str(),
        sizeCoordinate);

    // Find the type of tag, and use it
    DataType type;
    rval = mbImpl->tag_get_data_type(coordTag, type);MB_CHK_SET_ERR(rval, "Can't get data type of tag " << coordName);
    varCoordData.varDataType = NC_DOUBLE;
    if (MB_TYPE_INTEGER == type)
      varCoordData.varDataType = NC_INT;

    // Get dimension length (the only dimension of this coordinate variable, with the same name)
    assert(1 == varCoordData.varDims.size());
    int coordDimLen = dimLens[varCoordData.varDims[0]];

    if (dummyVarNames.find(coordName) != dummyVarNames.end()) {
      // For a dummy coordinate variable, the tag size is always 1
      // The number of coordinates should be set to dimension length, instead of 1
      assert(1 == sizeCoordinate);
      sizeCoordinate = coordDimLen;

      // No variable data to write
      data = NULL;
    }
    else {
      // The number of coordinates should be exactly the same as dimension length
      // However, if timesteps spread across files and time tag has been updated,
      // sizeCoordinate will be larger
      if (varCoordData.varDims[0] != tDim)
        assert(sizeCoordinate == coordDimLen);
    }

    // For time, the actual output size and values are determined by tstep_nums
    if (varCoordData.varDims[0] == tDim) {
      // Does not apply to dummy time tag (e.g. 'Time' tag of MPAS), when timesteps
      // spread across files
      if (NULL != data)
        assert(tstep_nums.size() > 0 && tstep_nums.size() <= (size_t)sizeCoordinate);

      sizeCoordinate = tstep_nums.size();

      if (NULL != data) {
        assert(NC_DOUBLE == varCoordData.varDataType);
        timeStepVals.resize(sizeCoordinate);
        for (unsigned int t = 0; t < tstep_nums.size(); t++)
          timeStepVals[t] = ((double*)data)[tstep_nums[t]];

        data = &timeStepVals[0];
      }
    }

    // This is the length
    varCoordData.sz = sizeCoordinate;
    varCoordData.writeStarts.resize(1);
    varCoordData.writeStarts[0] = 0;
    varCoordData.writeCounts.resize(1);
    varCoordData.writeCounts[0] = sizeCoordinate;

    assert(0 == varCoordData.memoryHogs.size()); // Nothing so far
    varCoordData.memoryHogs.push_back((void*)data);
  } // for (std::set<std::string>::iterator setIt ...

  return MB_SUCCESS;
}

ErrorCode NCWriteHelper::init_file(std::vector<std::string>& var_names, std::vector<std::string>& desired_names, bool append)
{
  std::vector<std::string>& dimNames = _writeNC->dimNames;
  std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
  std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
  std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;
  std::map<std::string, WriteNC::AttData>& globalAtts = _writeNC->globalAtts;
  DebugOutput& dbgOut = _writeNC->dbgOut;

  int tDim_in_dimNames = tDim;
  int levDim_in_dimNames = levDim;

  // If append mode, make sure we are in define mode; a simple open will not allow creation of new variables
  if (append) {
    int errcode = NCFUNC(redef)(_fileId);
    if (errcode != NC_NOERR)
      MB_SET_ERR(MB_FAILURE, "Can't open file in redefine mode");
  }

  // First initialize all coordinates, then fill VarData for actual variables (and dimensions)
  // Check that for used coordinates we have found the tags
  for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
      setIt != usedCoordinates.end(); ++setIt) {
    const std::string& coordName = *setIt;

    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(coordName);
    if (vit == varInfo.end())
      MB_SET_ERR(MB_FAILURE, "Can't find coordinate variable " << coordName);

    WriteNC::VarData& varCoordData = vit->second;
    varCoordData.varDims.resize(1);

    // If not append, create it for sure
    // If append, we might already have it, including the tag / variable with the same name
    /*
     * int ncmpi_inq_dimid(int ncid, const char *name, int *idp);
     */
    if (append) {
      int dimId;
      if (NCFUNC(inq_dimid)(_fileId, coordName.c_str(), &dimId) == NC_NOERR) { // If not found, create it later
        varCoordData.varDims[0] = dimId;
        dbgOut.tprintf(2, "    file already has coordName %s dim id is %d \n", coordName.c_str(), (int)varCoordData.varDims[0]);

        // Update tDim and levDim to actual dimension id
        if (coordName == dimNames[tDim_in_dimNames])
          tDim = varCoordData.varDims[0];
        else if (coordName == dimNames[levDim_in_dimNames])
          levDim = varCoordData.varDims[0];

        // Skip dummy coordinate variables (e.g. ncol)
        if (dummyVarNames.find(coordName) != dummyVarNames.end())
          continue;

        // Check that the coordinate is a variable too
        // Inquire for a variable with the same name
        int varId;
        if (NCFUNC(inq_varid)(_fileId, coordName.c_str(), &varId) != NC_NOERR)
          MB_SET_ERR(MB_FAILURE, "We do not have a variable with the same name " << coordName);
        // We should also check that this variable has one dimension, and it is dimId
        varCoordData.varId = varId;
        dbgOut.tprintf(2, "    file already has coordinate %s and varId is %d \n", coordName.c_str(), varId);

        continue; // Maybe more checks are needed here
      }
    }

    /* int nc_def_dim (int ncid, const char *name, size_t len, int *dimidp);
       * example:  status = nc_def_dim(fileId, "lat", 18L, &latid);
    */

    // Actually define a dimension
    if (NCFUNC(def_dim)(_fileId, coordName.c_str(), (size_t)varCoordData.sz,
        &varCoordData.varDims[0]) != NC_NOERR)
      MB_SET_ERR(MB_FAILURE, "Failed to generate dimension " << coordName);
    dbgOut.tprintf(2, "    for coordName %s dim id is %d \n", coordName.c_str(), (int)varCoordData.varDims[0]);

    // Update tDim and levDim to actual dimension id
    if (coordName == dimNames[tDim_in_dimNames])
      tDim = varCoordData.varDims[0];
    else if (coordName == dimNames[levDim_in_dimNames])
      levDim = varCoordData.varDims[0];

    // Create a variable with the same name, and its only dimension the one we just defined
    /*
     * int nc_def_var (int ncid, const char *name, nc_type xtype,
                       int ndims, const int dimids[], int *varidp);
       example: http://www.unidata.ucar.edu/software/netcdf/docs/netcdf-c/nc_005fdef_005fvar.html#nc_005fdef_005fvar
     */

    // Skip dummy coordinate variables (e.g. ncol)
    if (dummyVarNames.find(coordName) != dummyVarNames.end())
      continue;

    // Define a coordinate variable
    if (NCFUNC(def_var)(_fileId, coordName.c_str(), varCoordData.varDataType,
        1, &(varCoordData.varDims[0]), &varCoordData.varId) != NC_NOERR)
      MB_SET_ERR(MB_FAILURE, "Failed to create coordinate variable " << coordName);

    dbgOut.tprintf(2, "    for coordName %s variable id is %d \n", coordName.c_str(), varCoordData.varId);
  }

  // Now look at requested variables, and update from the index in dimNames to the actual dimension id
  for (size_t i = 0; i < var_names.size(); i++) {
    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(var_names[i]);
    if (vit == varInfo.end())
      MB_SET_ERR(MB_FAILURE, "Can't find requested variable " << var_names[i]);

    // Skip coordinate variables
    if (usedCoordinates.find(var_names[i]) != usedCoordinates.end())
      continue;

    WriteNC::VarData& variableData = vit->second;

    // The index is for dimNames; we need to find out the actual dimension id (from above)
    int numDims = (int)variableData.varDims.size();
    for (int j = 0; j < numDims; j++) {
      std::string dimName = dimNames[variableData.varDims[j]];
      std::map<std::string, WriteNC::VarData>::iterator vit2 = varInfo.find(dimName);
      if (vit2 == varInfo.end())
        MB_SET_ERR(MB_FAILURE, "Can't find requested coordinate variable " << dimName);

      WriteNC::VarData& coordData = vit2->second;
      // Index in dimNames to actual dimension id
      variableData.varDims[j] = coordData.varDims[0]; // This one, being a coordinate, is the only one
      dbgOut.tprintf(2, "          dimension with index %d name %s has ID %d \n",
          j, dimName.c_str(), variableData.varDims[j]);
    }

    // Define the variable now:
    int errCode = NCFUNC(def_var)(_fileId, desired_names[i].c_str(), variableData.varDataType,
        (int)variableData.varDims.size(), &(variableData.varDims[0]),
        &variableData.varId);
    if (errCode != NC_NOERR)
      MB_SET_ERR(MB_FAILURE, "Failed to create requested variable " << desired_names[i]);

    dbgOut.tprintf(2, "    for variable %s with desired name %s variable id is %d \n", var_names[i].c_str(),
        desired_names[i].c_str(), variableData.varId);
    // Now define the variable, with all dimensions
  }

  // Define global attributes (exactly copied from the original file for the time being)
  // Should we modify some of them (e.g. revision_Id) later?
  std::map<std::string, WriteNC::AttData>::iterator attIt;
  for (attIt = globalAtts.begin(); attIt != globalAtts.end(); ++attIt) {
    const std::string& attName = attIt->first;
    WriteNC::AttData& attData = attIt->second;
    NCDF_SIZE& attLen = attData.attLen;
    nc_type& attDataType = attData.attDataType;
    const std::string& attValue = attData.attValue;

    switch (attDataType) {
      case NC_BYTE:
      case NC_CHAR:
        if (NC_NOERR != NCFUNC(put_att_text)(_fileId, NC_GLOBAL, attName.c_str(), attLen, attValue.c_str()))
          MB_SET_ERR(MB_FAILURE, "Failed to define text type attribute");
        break;
      case NC_DOUBLE:
        if (NC_NOERR != NCFUNC(put_att_double)(_fileId, NC_GLOBAL, attName.c_str(), NC_DOUBLE, 1, (double*)attValue.c_str()))
          MB_SET_ERR(MB_FAILURE, "Failed to define double type attribute");
        break;
      case NC_FLOAT:
        if (NC_NOERR != NCFUNC(put_att_float)(_fileId, NC_GLOBAL, attName.c_str(), NC_FLOAT, 1, (float*)attValue.c_str()))
          MB_SET_ERR(MB_FAILURE, "Failed to define float type attribute");
        break;
      case NC_INT:
        if (NC_NOERR != NCFUNC(put_att_int)(_fileId, NC_GLOBAL, attName.c_str(), NC_INT, 1, (int*)attValue.c_str()))
          MB_SET_ERR(MB_FAILURE, "Failed to define int type attribute");
        break;
      case NC_SHORT:
        if (NC_NOERR != NCFUNC(put_att_short)(_fileId, NC_GLOBAL, attName.c_str(), NC_SHORT, 1, (short*)attValue.c_str()))
          MB_SET_ERR(MB_FAILURE, "Failed to define short type attribute");
        break;
      default:
        MB_SET_ERR(MB_FAILURE, "Unknown attribute data type");
    }
  }

  // Take it out of define mode
  if (NC_NOERR != NCFUNC(enddef)(_fileId))
    MB_SET_ERR(MB_FAILURE, "Failed to close define mode");

  return MB_SUCCESS;
}

ErrorCode NCWriteHelper::write_values(std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
{
  std::set<std::string>& usedCoordinates = _writeNC->usedCoordinates;
  std::set<std::string>& dummyVarNames = _writeNC->dummyVarNames;
  std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;

  std::vector<WriteNC::VarData> vdatas;
  std::vector<WriteNC::VarData> vsetdatas;

  // For set variables, include coordinates used by requested var_names
  for (std::set<std::string>::iterator setIt = usedCoordinates.begin();
      setIt != usedCoordinates.end(); ++setIt) {
    const std::string& coordName = *setIt;

    // Skip dummy coordinate variables (if any)
    if (dummyVarNames.find(coordName) != dummyVarNames.end())
      continue;

    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(coordName);
    if (vit == varInfo.end()) {
      MB_SET_ERR(MB_FAILURE, "Can't find coordinate variable " << coordName);
    }

     vsetdatas.push_back(vit->second);
  }

  // Collect non-set and set variables from requested var_names
  for (unsigned int i = 0; i < var_names.size(); i++) {
    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(var_names[i]);
    if (vit == varInfo.end()) {
      MB_SET_ERR(MB_FAILURE, "Can't find requested variable " << var_names[i]);
    }

    WriteNC::VarData& variableData = vit->second;
    if (WriteNC::ENTLOCSET == variableData.entLoc) {
      // Used coordinates has all ready been included
      if (usedCoordinates.find(var_names[i]) != usedCoordinates.end())
        continue;

      vsetdatas.push_back(variableData);
    }
    else
      vdatas.push_back(variableData);
  }

  // Assume that the data ranges do not overlap across processors
  // While overlapped writing might still work, we should better not take that risk
  write_nonset_variables(vdatas, tstep_nums);

  // Use independent I/O mode put, since this write is only for the root processor
  write_set_variables(vsetdatas, tstep_nums);

  return MB_SUCCESS;
}

ErrorCode NCWriteHelper::write_set_variables(std::vector<WriteNC::VarData>& vsetdatas, std::vector<int>& /* tstep_nums */)
{
  int success;

// CAUTION: if the NetCDF ID is from a previous call to ncmpi_create rather than ncmpi_open,
// all processors need to call ncmpi_begin_indep_data(). If only the root processor does so,
// ncmpi_begin_indep_data() call will be blocked forever :(
 #ifdef MOAB_HAVE_PNETCDF
   // Enter independent I/O mode
   success = NCFUNC(begin_indep_data)(_fileId);
   if (success)
     MB_SET_ERR(MB_FAILURE, "Failed to begin independent I/O mode");
 #endif

   int rank = 0;
 #ifdef MOAB_HAVE_MPI
   bool& isParallel = _writeNC->isParallel;
   if (isParallel) {
     ParallelComm*& myPcomm = _writeNC->myPcomm;
     rank = myPcomm->proc_config().proc_rank();
   }
 #endif
   if (0 == rank) {
     for (unsigned int i = 0; i < vsetdatas.size(); i++) {
       WriteNC::VarData& variableData = vsetdatas[i];

       // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
       if (variableData.has_tsteps) {
         MB_SET_ERR(MB_NOT_IMPLEMENTED, "Writing set variables with timesteps is not implemented yet");
       }

       switch (variableData.varDataType) {
         case NC_DOUBLE:
           // Independent I/O mode put
           success = NCFUNCP(_vara_double)(_fileId, variableData.varId, &variableData.writeStarts[0],
                     &variableData.writeCounts[0], (double*)(variableData.memoryHogs[0]));
           if (success)
             MB_SET_ERR(MB_FAILURE, "Failed to write double data for variable " << variableData.varName);
           break;
         case NC_INT:
           // Independent I/O mode put
           success = NCFUNCP(_vara_int)(_fileId, variableData.varId, &variableData.writeStarts[0],
                     &variableData.writeCounts[0], (int*)(variableData.memoryHogs[0]));
           if (success)
             MB_SET_ERR(MB_FAILURE, "Failed to write int data for variable " << variableData.varName);
           break;
         default:
           MB_SET_ERR(MB_NOT_IMPLEMENTED, "Writing non-double or non-int data is not implemented yet");
       }
     }
   }

 #ifdef MOAB_HAVE_PNETCDF
   // End independent I/O mode
   success = NCFUNC(end_indep_data)(_fileId);
   if (success)
     MB_SET_ERR(MB_FAILURE, "Failed to end independent I/O mode");
 #endif

   return MB_SUCCESS;
}

ErrorCode ScdNCWriteHelper::collect_mesh_info()
{
  Interface*& mbImpl = _writeNC->mbImpl;
  std::vector<std::string>& dimNames = _writeNC->dimNames;
  std::vector<int>& dimLens = _writeNC->dimLens;

  ErrorCode rval;

  // Look for time dimension
  std::vector<std::string>::iterator vecIt;
  if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
    tDim = vecIt - dimNames.begin();
  else if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
    tDim = vecIt - dimNames.begin();
  else {
    MB_SET_ERR(MB_FAILURE, "Couldn't find 'time' or 't' dimension");
  }
  nTimeSteps = dimLens[tDim];

  // Get number of levels
  if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
    levDim = vecIt - dimNames.begin();
  else if ((vecIt = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
    levDim = vecIt - dimNames.begin();
  else {
    MB_SET_ERR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension");
  }
  nLevels = dimLens[levDim];

  // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
  Tag convTag = 0;
  rval = mbImpl->tag_get_handle("__slon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag,
                                MB_TAG_ANY);MB_CHK_SET_ERR(rval, "Trouble getting conventional tag __slon_LOC_MINMAX");
  int val[2];
  rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);MB_CHK_SET_ERR(rval, "Trouble getting data of conventional tag __slon_LOC_MINMAX");
  lDims[0] = val[0];
  lDims[3] = val[1];

  rval = mbImpl->tag_get_handle("__slat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag,
                                MB_TAG_ANY);MB_CHK_SET_ERR(rval, "Trouble getting conventional tag __slat_LOC_MINMAX");
  rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);MB_CHK_SET_ERR(rval, "Trouble getting data of conventional tag __slat_LOC_MINMAX");
  lDims[1] = val[0];
  lDims[4] = val[1];

  rval = mbImpl->tag_get_handle("__lon_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag,
                                MB_TAG_ANY);MB_CHK_SET_ERR(rval, "Trouble getting conventional tag __lon_LOC_MINMAX");
  rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);MB_CHK_SET_ERR(rval, "Trouble getting data of conventional tag __lon_LOC_MINMAX");
  lCDims[0] = val[0];
  lCDims[3] = val[1];

  rval = mbImpl->tag_get_handle("__lat_LOC_MINMAX", 0, MB_TYPE_INTEGER, convTag,
                                MB_TAG_ANY);MB_CHK_SET_ERR(rval, "Trouble getting conventional tag __lat_LOC_MINMAX");
  rval = mbImpl->tag_get_data(convTag, &_fileSet, 1, val);MB_CHK_SET_ERR(rval, "Trouble getting data of conventional tag __lat_LOC_MINMAX");
  lCDims[1] = val[0];
  lCDims[4] = val[1];

  // Get local faces
  rval = mbImpl->get_entities_by_dimension(_fileSet, 2, localCellsOwned);MB_CHK_SET_ERR(rval, "Trouble getting local faces in current file set");
  assert(!localCellsOwned.empty());

#ifdef MOAB_HAVE_MPI
  bool& isParallel = _writeNC->isParallel;
  if (isParallel) {
    ParallelComm*& myPcomm = _writeNC->myPcomm;
    int procs = myPcomm->proc_config().proc_size();
    if (procs > 1) {
      rval = myPcomm->filter_pstatus(localCellsOwned, PSTATUS_NOT_OWNED, PSTATUS_NOT);MB_CHK_SET_ERR(rval, "Trouble getting owned faces in current file set");
    }
  }
#endif

  return MB_SUCCESS;
}

ErrorCode ScdNCWriteHelper::collect_variable_data(std::vector<std::string>& var_names, std::vector<int>& tstep_nums)
{
  NCWriteHelper::collect_variable_data(var_names, tstep_nums);

  std::map<std::string, WriteNC::VarData>& varInfo = _writeNC->varInfo;

  for (size_t i = 0; i < var_names.size(); i++) {
    std::string varname = var_names[i];
    std::map<std::string, WriteNC::VarData>::iterator vit = varInfo.find(varname);
    if (vit == varInfo.end())
      MB_SET_ERR(MB_FAILURE, "Can't find variable " << varname);

    WriteNC::VarData& currentVarData = vit->second;
#ifndef NDEBUG
    std::vector<int>& varDims = currentVarData.varDims;
#endif

    // Skip set variables, which were already processed in NCWriteHelper::collect_variable_data()
    if (WriteNC::ENTLOCSET == currentVarData.entLoc)
      continue;

    // Set up writeStarts and writeCounts (maximum number of dimensions is 4)
    currentVarData.writeStarts.resize(4);
    currentVarData.writeCounts.resize(4);
    unsigned int dim_idx = 0;

    // First: time
    if (currentVarData.has_tsteps) {
      // Non-set variables with timesteps
      // 4 dimensions like (time, lev, lat, lon)
      // 3 dimensions like (time, lat, lon)
      assert(4 == varDims.size() || 3 == varDims.size());

      // Time should be the first dimension
      assert(tDim == varDims[0]);

      currentVarData.writeStarts[dim_idx] = 0; // This value is timestep dependent, will be set later
      currentVarData.writeCounts[dim_idx] = 1;
      dim_idx++;
    }
    else {
      // Non-set variables without timesteps
      // 3 dimensions like (lev, lat, lon)
      // 2 dimensions like (lat, lon)
      assert(3 == varDims.size() || 2 == varDims.size());
    }

    // Next: lev
    if (currentVarData.numLev > 0) {
      // Non-set variables with levels
      // 4 dimensions like (time, lev, lat, lon)
      // 3 dimensions like (lev, lat, lon)
      assert(4 == varDims.size() || 3 == varDims.size());

      currentVarData.writeStarts[dim_idx] = 0;
      currentVarData.writeCounts[dim_idx] = currentVarData.numLev;
      dim_idx++;
    }
    else {
      // Non-set variables without levels
      // 3 dimensions like (time, lat, lon)
      // 2 dimensions like (lat, lon)
      assert(3 == varDims.size() || 2 == varDims.size());
    }

    // Finally: lat and lon
    switch (currentVarData.entLoc) {
      case WriteNC::ENTLOCFACE:
        // Faces
        currentVarData.writeStarts[dim_idx] = lCDims[1];
        currentVarData.writeCounts[dim_idx] = lCDims[4] - lCDims[1] + 1;
        currentVarData.writeStarts[dim_idx + 1] = lCDims[0];
        currentVarData.writeCounts[dim_idx + 1] = lCDims[3] - lCDims[0] + 1;
        break;
      default:
        MB_SET_ERR(MB_FAILURE, "Unexpected entity location type for variable " << varname);
    }
    dim_idx += 2;

    // Get variable size
    currentVarData.sz = 1;
    for (std::size_t idx = 0; idx < dim_idx; idx++)
      currentVarData.sz *= currentVarData.writeCounts[idx];
  } // for (size_t i = 0; i < var_names.size(); i++)

  return MB_SUCCESS;
}

// Write CAM-EUL and CAM-FV non-set variables on non-shared quads (e.g. T)
// We assume that there are no variables on vertices and we do not support
// variables on edges (e.g. US in CAM-FV) for the time being
ErrorCode ScdNCWriteHelper::write_nonset_variables(std::vector<WriteNC::VarData>& vdatas, std::vector<int>& tstep_nums)
{
  Interface*& mbImpl = _writeNC->mbImpl;

  int success;

  // For each indexed variable tag, write a time step data
  for (unsigned int i = 0; i < vdatas.size(); i++) {
    WriteNC::VarData& variableData = vdatas[i];

    // Assume this variable is on faces for the time being
    switch (variableData.entLoc) {
      case WriteNC::ENTLOCFACE:
        // Faces
        break;
      default:
        MB_SET_ERR(MB_FAILURE, "Unexpected entity location type for variable " << variableData.varName);
    }

    unsigned int num_timesteps;
    unsigned int lat_idx = 0;
    unsigned int lon_idx = 1;
    if (variableData.has_tsteps) {
      // Non-set variables with timesteps
      // 4 dimensions like (time, lev, lat, lon)
      // 3 dimensions like (time, lat, lon)
      num_timesteps = tstep_nums.size();
      lat_idx++;
      lon_idx++;
    }
    else {
      // Non-set variables without timesteps
      // 3 dimensions like (lev, lat, lon)
      // 2 dimensions like (lat, lon)
      num_timesteps = 1;
    }

    unsigned int num_lev;
    if (variableData.numLev > 0) {
      // Non-set variables with levels
      // 4 dimensions like (time, lev, lat, lon)
      // 3 dimensions like (lev, lat, lon)
      num_lev = variableData.numLev;
      lat_idx++;
      lon_idx++;
    }
    else {
      // Non-set variables without levels
      // 3 dimensions like (time, lat, lon)
      // 2 dimensions like (lat, lon)
      num_lev = 1;
    }

    size_t ni = variableData.writeCounts[lon_idx]; // lon
    size_t nj = variableData.writeCounts[lat_idx]; // lat

    // At each timestep, we need to transpose tag format (lat, lon, lev) back
    // to NC format (lev, lat, lon) for writing
    for (unsigned int t = 0; t < num_timesteps; t++) {
      // We will write one time step, and count will be one; start will be different
      // Use tag_iterate to get tag data (assume that localCellsOwned is contiguous)
      // We should also transpose for level so that means deep copy for transpose
      if (tDim == variableData.varDims[0])
        variableData.writeStarts[0] = t; // This is start for time
      int count;
      void* dataptr;
      ErrorCode rval = mbImpl->tag_iterate(variableData.varTags[t], localCellsOwned.begin(), localCellsOwned.end(),
                                           count, dataptr);MB_CHK_SET_ERR(rval, "Failed to iterate tag on owned faces");
      assert(count == (int)localCellsOwned.size());

      // Now transpose and write tag data
      // Use collective I/O mode put (synchronous write) for the time being, we can try
      // nonblocking put (request aggregation) later
      switch (variableData.varDataType) {
        case NC_DOUBLE: {
          std::vector<double> tmpdoubledata(ni*nj*num_lev);
          if (num_lev > 1)
            // Transpose (lat, lon, lev) back to (lev, lat, lon)
            jik_to_kji(ni, nj, num_lev, &tmpdoubledata[0], (double*)(dataptr));
          success = NCFUNCAP(_vara_double)(_fileId, variableData.varId,
                    &variableData.writeStarts[0], &variableData.writeCounts[0],
                    &tmpdoubledata[0]);
          if (success)
            MB_SET_ERR(MB_FAILURE, "Failed to write double data for variable " << variableData.varName);
          break;
        }
        default:
          MB_SET_ERR(MB_NOT_IMPLEMENTED, "Writing non-double data is not implemented yet");
      }
    }
  }

  return MB_SUCCESS;
}

} /* namespace moab */