fathom/moab-docs-develop/NCHelperEuler_8cpp_source.html

00001 #include "moab/MOABConfig.h"
00002 #ifdef WIN32               /* windows */
00003 #define _USE_MATH_DEFINES  // For M_PI
00004 #endif
00005
00006 #include "NCHelperEuler.hpp"
00007 #include "moab/FileOptions.hpp"
00008
00009 #include
00010 #include
00011
00012 #ifdef WIN32
00013 #ifdef size_t
00014 #undef size_t
00015 #endif
00016 #endif
00017
00018 namespace moab
00019 {
00020
00021 bool NCHelperEuler::can_read_file( ReadNC* readNC, int fileId )
00022 {
00023     std::vector< std::string >& dimNames = readNC->dimNames;
00024
00025     // If dimension names "lon" AND "lat' exist then it could be either the Eulerian Spectral grid
00026     // or the FV grid
00027     if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "lon" ) ) != dimNames.end() ) &&
00028         ( std::find( dimNames.begin(), dimNames.end(), std::string( "lat" ) ) != dimNames.end() ) )
00029     {
00030         // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV
00031         // grid
00032         if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "slon" ) ) != dimNames.end() ) &&
00033             ( std::find( dimNames.begin(), dimNames.end(), std::string( "slat" ) ) != dimNames.end() ) )
00034             return false;
00035
00036         // Make sure it is CAM grid
00037         std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
00038         if( attIt == readNC->globalAtts.end() ) return false;
00039         unsigned int sz = attIt->second.attLen;
00040         std::string att_data;
00041         att_data.resize( sz + 1 );
00042         att_data[sz] = '\000';
00043         int success =
00044             NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
00045         if( success ) return false;
00046         if( att_data.find( "CAM" ) == std::string::npos ) return false;
00047
00048         return true;
00049     }
00050
00051     return false;
00052 }
00053
00054 ErrorCode NCHelperEuler::init_mesh_vals()
00055 {
00056     Interface*& mbImpl                                = _readNC->mbImpl;
00057     std::vector< std::string >& dimNames              = _readNC->dimNames;
00058     std::vector< int >& dimLens                       = _readNC->dimLens;
00059     std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00060     DebugOutput& dbgOut                               = _readNC->dbgOut;
00061     bool& isParallel                                  = _readNC->isParallel;
00062     int& partMethod                                   = _readNC->partMethod;
00063     ScdParData& parData                               = _readNC->parData;
00064
00065     // Look for names of center i/j dimensions
00066     // First i
00067     std::vector< std::string >::iterator vit;
00068     unsigned int idx;
00069     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lon" ) ) != dimNames.end() )
00070         idx = vit - dimNames.begin();
00071     else
00072     {
00073         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' dimension" );
00074     }
00075     iCDim     = idx;
00076     gCDims[0] = 0;
00077     gCDims[3] = dimLens[idx] - 1;
00078
00079     // Check i periodicity and set globallyPeriodic[0]
00080     std::vector< double > til_vals( 2 );
00081     ErrorCode rval = read_coordinate( "lon", gCDims[3] - 1, gCDims[3], til_vals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
00082     if( std::fabs( 2 * til_vals[1] - til_vals[0] - 360 ) < 0.001 ) globallyPeriodic[0] = 1;
00083
00084     // Now we can set gDims for i
00085     gDims[0] = 0;
00086     gDims[3] =
00087         gCDims[3] + ( globallyPeriodic[0] ? 0 : 1 );  // Only if not periodic is vertex param max > elem param max
00088
00089     // Then j
00090     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lat" ) ) != dimNames.end() )
00091         idx = vit - dimNames.begin();
00092     else
00093     {
00094         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' dimension" );
00095     }
00096     jCDim     = idx;
00097     gCDims[1] = 0;
00098     gCDims[4] = dimLens[idx] - 1;
00099
00100     // For Eul models, will always be non-periodic in j
00101     gDims[1] = 0;
00102     gDims[4] = gCDims[4] + 1;
00103
00104     // Try a truly 2D mesh
00105     gDims[2] = -1;
00106     gDims[5] = -1;
00107
00108     // Look for time dimension
00109     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
00110         idx = vit - dimNames.begin();
00111     else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
00112         idx = vit - dimNames.begin();
00113     else
00114     {
00115         MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
00116     }
00117     tDim       = idx;
00118     nTimeSteps = dimLens[idx];
00119
00120     // Look for level dimension
00121     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
00122         idx = vit - dimNames.begin();
00123     else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
00124         idx = vit - dimNames.begin();
00125     else
00126     {
00127         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
00128     }
00129     levDim  = idx;
00130     nLevels = dimLens[idx];
00131
00132     // Parse options to get subset
00133     int rank = 0, procs = 1;
00134 #ifdef MOAB_HAVE_MPI
00135     if( isParallel )
00136     {
00137         ParallelComm*& myPcomm = _readNC->myPcomm;
00138         rank                   = myPcomm->proc_config().proc_rank();
00139         procs                  = myPcomm->proc_config().proc_size();
00140     }
00141 #endif
00142     if( procs > 1 )
00143     {
00144         for( int i = 0; i < 6; i++ )
00145             parData.gDims[i] = gDims[i];
00146         for( int i = 0; i < 3; i++ )
00147             parData.gPeriodic[i] = globallyPeriodic[i];
00148         parData.partMethod = partMethod;
00149         int pdims[3];
00150
00151         rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval );
00152         for( int i = 0; i < 3; i++ )
00153             parData.pDims[i] = pdims[i];
00154
00155         dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
00156                         gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1 );
00157         if( 0 == rank )
00158             dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
00159                             8 * ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
00160     }
00161     else
00162     {
00163         for( int i = 0; i < 6; i++ )
00164             lDims[i] = gDims[i];
00165         locallyPeriodic[0] = globallyPeriodic[0];
00166     }
00167
00168     _opts.get_int_option( "IMIN", lDims[0] );
00169     _opts.get_int_option( "IMAX", lDims[3] );
00170     _opts.get_int_option( "JMIN", lDims[1] );
00171     _opts.get_int_option( "JMAX", lDims[4] );
00172
00173     // Now get actual coordinate values for vertices and cell centers
00174     lCDims[0] = lDims[0];
00175     if( locallyPeriodic[0] )
00176         // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem
00177         // coords
00178         lCDims[3] = lDims[3];
00179     else
00180         lCDims[3] = lDims[3] - 1;
00181
00182     // For Eul models, will always be non-periodic in j
00183     lCDims[1] = lDims[1];
00184     lCDims[4] = lDims[4] - 1;
00185
00186     // Resize vectors to store values later
00187     if( -1 != lDims[0] ) ilVals.resize( lDims[3] - lDims[0] + 1 );
00188     if( -1 != lCDims[0] ) ilCVals.resize( lCDims[3] - lCDims[0] + 1 );
00189     if( -1 != lDims[1] ) jlVals.resize( lDims[4] - lDims[1] + 1 );
00190     if( -1 != lCDims[1] ) jlCVals.resize( lCDims[4] - lCDims[1] + 1 );
00191     if( nTimeSteps > 0 ) tVals.resize( nTimeSteps );
00192
00193     // Now read coord values
00194     std::map< std::string, ReadNC::VarData >::iterator vmit;
00195     if( -1 != lCDims[0] )
00196     {
00197         if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00198         {
00199             rval = read_coordinate( "lon", lCDims[0], lCDims[3], ilCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
00200         }
00201         else
00202         {
00203             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
00204         }
00205     }
00206
00207     if( -1 != lCDims[1] )
00208     {
00209         if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00210         {
00211             rval = read_coordinate( "lat", lCDims[1], lCDims[4], jlCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" );
00212         }
00213         else
00214         {
00215             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
00216         }
00217     }
00218
00219     if( -1 != lDims[0] )
00220     {
00221         if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00222         {
00223             double dif = ( ilCVals[1] - ilCVals[0] ) / 2;
00224             std::size_t i;
00225             for( i = 0; i != ilCVals.size(); i++ )
00226                 ilVals[i] = ilCVals[i] - dif;
00227             // The last one is needed only if not periodic
00228             if( !locallyPeriodic[0] ) ilVals[i] = ilCVals[i - 1] + dif;
00229         }
00230         else
00231         {
00232             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
00233         }
00234     }
00235
00236     if( -1 != lDims[1] )
00237     {
00238         if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00239         {
00240             if( !isParallel || ( ( gDims[4] - gDims[1] ) == ( lDims[4] - lDims[1] ) ) )
00241             {
00242                 std::string gwName( "gw" );
00243                 std::vector< double > gwVals( lDims[4] - lDims[1] - 1 );
00244                 rval = read_coordinate( gwName.c_str(), lDims[1], lDims[4] - 2, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
00245                 // Copy the correct piece
00246                 jlVals[0]     = -( M_PI / 2 ) * 180 / M_PI;
00247                 std::size_t i = 0;
00248                 double gwSum  = -1;
00249                 for( i = 1; i != gwVals.size() + 1; i++ )
00250                 {
00251                     gwSum += gwVals[i - 1];
00252                     jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
00253                 }
00254                 jlVals[i] = 90.0;  // Using value of i after loop exits.
00255             }
00256             else
00257             {
00258                 std::string gwName( "gw" );
00259                 double gwSum = 0;
00260
00261                 // If this is the first row
00262                 if( lDims[1] == gDims[1] )
00263                 {
00264                     std::vector< double > gwVals( lDims[4] );
00265                     rval = read_coordinate( gwName.c_str(), 0, lDims[4] - 1, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
00266                     // Copy the correct piece
00267                     jlVals[0] = -( M_PI / 2 ) * 180 / M_PI;
00268                     gwSum     = -1;
00269                     for( std::size_t i = 1; i != jlVals.size(); i++ )
00270                     {
00271                         gwSum += gwVals[i - 1];
00272                         jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
00273                     }
00274                 }
00275                 // Or if it's the last row
00276                 else if( lDims[4] == gDims[4] )
00277                 {
00278                     std::vector< double > gwVals( lDims[4] - 1 );
00279                     rval = read_coordinate( gwName.c_str(), 0, lDims[4] - 2, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
00280                     // Copy the correct piece
00281                     gwSum = -1;
00282                     for( int j = 0; j != lDims[1] - 1; j++ )
00283                         gwSum += gwVals[j];
00284                     std::size_t i = 0;
00285                     for( ; i != jlVals.size() - 1; i++ )
00286                     {
00287                         gwSum += gwVals[lDims[1] - 1 + i];
00288                         jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
00289                     }
00290                     jlVals[i] = 90.0;  // Using value of i after loop exits.
00291                 }
00292                 // It's in the middle
00293                 else
00294                 {
00295                     int start = lDims[1] - 1;
00296                     int end   = lDims[4] - 1;
00297                     std::vector< double > gwVals( end );
00298                     rval = read_coordinate( gwName.c_str(), 0, end - 1, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
00299                     gwSum = -1;
00300                     for( int j = 0; j != start - 1; j++ )
00301                         gwSum += gwVals[j];
00302                     std::size_t i = 0;
00303                     for( ; i != jlVals.size(); i++ )
00304                     {
00305                         gwSum += gwVals[start - 1 + i];
00306                         jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
00307                     }
00308                 }
00309             }
00310         }
00311         else
00312         {
00313             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
00314         }
00315     }
00316
00317     // Store time coordinate values in tVals
00318     if( nTimeSteps > 0 )
00319     {
00320         if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00321         {
00322             rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
00323         }
00324         else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00325         {
00326             rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
00327         }
00328         else
00329         {
00330             // If expected time variable is not available, set dummy time coordinate values to tVals
00331             for( int t = 0; t < nTimeSteps; t++ )
00332                 tVals.push_back( (double)t );
00333         }
00334     }
00335
00336     dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
00337     dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
00338                     ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
00339
00340     // For each variable, determine the entity location type and number of levels
00341     std::map< std::string, ReadNC::VarData >::iterator mit;
00342     for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
00343     {
00344         ReadNC::VarData& vd = ( *mit ).second;
00345
00346         // Default entLoc is ENTLOCSET
00347         if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
00348         {
00349             if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
00350                 ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
00351                 vd.entLoc = ReadNC::ENTLOCFACE;
00352         }
00353
00354         // Default numLev is 0
00355         if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
00356     }
00357
00358     // For Eul models, slon and slat are "virtual" dimensions (not defined in the file header)
00359     std::vector< std::string > ijdimNames( 4 );
00360     ijdimNames[0] = "__slon";
00361     ijdimNames[1] = "__slat";
00362     ijdimNames[2] = "__lon";
00363     ijdimNames[3] = "__lat";
00364
00365     std::string tag_name;
00366     Tag tagh;
00367
00368     // ___LOC_MINMAX (for virtual slon, virtual slat, lon and lat)
00369     for( unsigned int i = 0; i != ijdimNames.size(); i++ )
00370     {
00371         std::vector< int > val( 2, 0 );
00372         if( ijdimNames[i] == "__slon" )
00373         {
00374             val[0] = lDims[0];
00375             val[1] = lDims[3];
00376         }
00377         else if( ijdimNames[i] == "__slat" )
00378         {
00379             val[0] = lDims[1];
00380             val[1] = lDims[4];
00381         }
00382         else if( ijdimNames[i] == "__lon" )
00383         {
00384             val[0] = lCDims[0];
00385             val[1] = lCDims[3];
00386         }
00387         else if( ijdimNames[i] == "__lat" )
00388         {
00389             val[0] = lCDims[1];
00390             val[1] = lCDims[4];
00391         }
00392         std::stringstream ss_tag_name;
00393         ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
00394         tag_name = ss_tag_name.str();
00395         rval     = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00396         rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00397         if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
00398     }
00399
00400     // ___LOC_VALS (for virtual slon, virtual slat, lon and lat)
00401     // Assume all have the same data type as lon (expected type is float or double)
00402     switch( varInfo["lon"].varDataType )
00403     {
00404         case NC_FLOAT:
00405         case NC_DOUBLE:
00406             break;
00407         default:
00408             MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
00409     }
00410
00411     for( unsigned int i = 0; i != ijdimNames.size(); i++ )
00412     {
00413         void* val   = NULL;
00414         int val_len = 0;
00415         if( ijdimNames[i] == "__slon" )
00416         {
00417             val     = &ilVals[0];
00418             val_len = ilVals.size();
00419         }
00420         else if( ijdimNames[i] == "__slat" )
00421         {
00422             val     = &jlVals[0];
00423             val_len = jlVals.size();
00424         }
00425         else if( ijdimNames[i] == "__lon" )
00426         {
00427             val     = &ilCVals[0];
00428             val_len = ilCVals.size();
00429         }
00430         else if( ijdimNames[i] == "__lat" )
00431         {
00432             val     = &jlCVals[0];
00433             val_len = jlCVals.size();
00434         }
00435
00436         std::stringstream ss_tag_name;
00437         ss_tag_name << ijdimNames[i] << "_LOC_VALS";
00438         tag_name = ss_tag_name.str();
00439         rval     = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_DOUBLE, tagh,
00440                                            MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00441         rval = mbImpl->tag_set_by_ptr( tagh, &_fileSet, 1, &val, &val_len );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00442         if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
00443     }
00444
00445     // ___GLOBAL_MINMAX (for virtual slon, virtual slat, lon and lat)
00446     for( unsigned int i = 0; i != ijdimNames.size(); i++ )
00447     {
00448         std::vector< int > val( 2, 0 );
00449         if( ijdimNames[i] == "__slon" )
00450         {
00451             val[0] = gDims[0];
00452             val[1] = gDims[3];
00453         }
00454         else if( ijdimNames[i] == "__slat" )
00455         {
00456             val[0] = gDims[1];
00457             val[1] = gDims[4];
00458         }
00459         else if( ijdimNames[i] == "__lon" )
00460         {
00461             val[0] = gCDims[0];
00462             val[1] = gCDims[3];
00463         }
00464         else if( ijdimNames[i] == "__lat" )
00465         {
00466             val[0] = gCDims[1];
00467             val[1] = gCDims[4];
00468         }
00469         std::stringstream ss_tag_name;
00470         ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
00471         tag_name = ss_tag_name.str();
00472         rval     = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00473         rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00474         if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
00475     }
00476
00477     // Hack: create dummy variables, if needed, for dimensions with no corresponding coordinate
00478     // variables
00479     rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
00480
00481     return MB_SUCCESS;
00482 }
00483
00484 }  // namespace moab