NCWriteHelper.cpp

Go to the documentation of this file.

/*
 * 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 

#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::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];

    // ___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 */