Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NCHelper.cpp
Go to the documentation of this file.
1 #include "NCHelper.hpp"
2 #include "NCHelperEuler.hpp"
3 #include "NCHelperFV.hpp"
4 #include "NCHelperDomain.hpp"
5 #include "NCHelperScrip.hpp"
6 #include "NCHelperHOMME.hpp"
7 #include "NCHelperMPAS.hpp"
8 #include "NCHelperGCRM.hpp"
9 #include "NCHelperESMF.hpp"
10 
11 #include <sstream>
12 
13 #include "MBTagConventions.hpp"
14 
15 #ifdef WIN32
16 #ifdef size_t
17 #undef size_t
18 #endif
19 #endif
20 
21 namespace moab
22 {
23 
24 ReadNC::NCFormatType NCHelper::get_nc_format( ReadNC* readNC, int fileId )
25 {
26  // Check if CF convention is being followed
27  bool is_CF = false;
28 
29  std::map< std::string, ReadNC::AttData >& globalAtts = readNC->globalAtts;
30  std::map< std::string, ReadNC::AttData >::iterator attIt = globalAtts.find( "conventions" );
31  if( attIt == globalAtts.end() ) attIt = globalAtts.find( "Conventions" );
32 
33  if( attIt != globalAtts.end() )
34  {
35  unsigned int sz = attIt->second.attLen;
36  std::string att_data;
37  att_data.resize( sz + 1 );
38  att_data[sz] = '\000';
39  int success =
40  NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
41  if( 0 == success && att_data.find( "CF" ) != std::string::npos ) is_CF = true;
42  }
43 
44  // Depending on the format, update FileOptions as well
45  if( NCHelperMPAS::can_read_file( readNC ) )
46  return ReadNC::NC_FORMAT_MPAS;
47  else if( NCHelperScrip::can_read_file( readNC, fileId ) )
48  return ReadNC::NC_FORMAT_SCRIP;
49  else if( NCHelperESMF::can_read_file( readNC ) )
50  return ReadNC::NC_FORMAT_ESMF;
51  else if( NCHelperDomain::can_read_file( readNC, fileId ) ) // && is_CF )
52  return ReadNC::NC_FORMAT_DOMAIN;
53  else if( NCHelperHOMME::can_read_file( readNC, fileId ) )
54  return ReadNC::NC_FORMAT_HOMME;
55  else if( NCHelperEuler::can_read_file( readNC, fileId ) && is_CF )
56  return ReadNC::NC_FORMAT_EULER;
57  else if( NCHelperGCRM::can_read_file( readNC ) )
58  return ReadNC::NC_FORMAT_GCRM;
59  else if( NCHelperFV::can_read_file( readNC, fileId ) && is_CF )
60  return ReadNC::NC_FORMAT_FV;
61  else // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
62  return ReadNC::NC_FORMAT_UNKNOWN_TYPE;
63 }
64 
65 //! Get appropriate file read options depending on the format of the NC file
66 std::string NCHelper::get_default_ncformat_options( ReadNC::NCFormatType format )
67 {
68  switch( format )
69  {
70  case moab::ReadNC::NC_FORMAT_GCRM: // GCRM format reader
71  case moab::ReadNC::NC_FORMAT_MPAS: // MPAS format reader
72  return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;"
73  "PARALLEL_RESOLVE_SHARED_ENTS;NO_EDGES;NO_MIXED_ELEMENTS;VARIABLE=;";
74  case moab::ReadNC::NC_FORMAT_SCRIP: // SCRIP format reader
75  return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;";
76  case moab::ReadNC::NC_FORMAT_ESMF: // ESMF unstructured format reader
77  return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;VARIABLE=;";
78  case moab::ReadNC::NC_FORMAT_DOMAIN: // Climate Domain reader
79  return "PARALLEL=READ_PART;PARTITION_METHOD=SQIJ;VARIABLE=;";
80  case moab::ReadNC::NC_FORMAT_HOMME: // HOMME format reader
81  return "PARALLEL=READ_PART;PARTITION_METHOD=TRIVIAL;PARALLEL_RESOLVE_SHARED_ENTS;";
82  case moab::ReadNC::NC_FORMAT_EULER: // Euler format reader
83  case moab::ReadNC::NC_FORMAT_FV: // FV climate format reader
84  return "PARALLEL=READ_PART;PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
85  "PARTITION_METHOD=SQIJ;VARIABLE=;";
86  default:
87  return "PARALLEL=READ_PART;PARTITION_METHOD=RCBZOLTAN;PARALLEL_RESOLVE_SHARED_ENTS;";
88  }
89 }
90 
91 NCHelper* NCHelper::get_nc_helper( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
92 {
93  ReadNC::NCFormatType nctype = NCHelper::get_nc_format( readNC, fileId );
94 
95  switch( nctype )
96  {
97  case ReadNC::NC_FORMAT_MPAS: // MPAS format reader
98  return new( std::nothrow ) NCHelperMPAS( readNC, fileId, opts, fileSet );
99  case ReadNC::NC_FORMAT_SCRIP: // SCRIP format reader
100  // SCRIP can be CF or non CF, if it comes from MPAS :)
101  return new( std::nothrow ) NCHelperScrip( readNC, fileId, opts, fileSet );
102  case ReadNC::NC_FORMAT_ESMF: // ESMF unstructured format reader
103  return new( std::nothrow ) NCHelperESMF( readNC, fileId, opts, fileSet );
104  case ReadNC::NC_FORMAT_DOMAIN: // Climate Domain reader
105  return new( std::nothrow ) NCHelperDomain( readNC, fileId, opts, fileSet );
106  case ReadNC::NC_FORMAT_HOMME: // HOMME format reader
107  // For a HOMME connectivity file, there might be no CF convention
108  return new( std::nothrow ) NCHelperHOMME( readNC, fileId, opts, fileSet );
109  case ReadNC::NC_FORMAT_GCRM: // GCRM format reader
110  return new( std::nothrow ) NCHelperGCRM( readNC, fileId, opts, fileSet );
111  case ReadNC::NC_FORMAT_EULER: // Euler format reader
112  return new( std::nothrow ) NCHelperEuler( readNC, fileId, opts, fileSet );
113  case ReadNC::NC_FORMAT_FV: // FV climate format reader
114  return new( std::nothrow ) NCHelperFV( readNC, fileId, opts, fileSet );
115  default: // Unknown NetCDF grid (will fill this in later for POP, CICE and CLM)
116  return nullptr;
117  }
118 }
119 
120 ErrorCode NCHelper::create_conventional_tags( const std::vector< int >& tstep_nums )
121 {
122  Interface*& mbImpl = _readNC->mbImpl;
123  std::vector< std::string >& dimNames = _readNC->dimNames;
124  std::vector< int >& dimLens = _readNC->dimLens;
125  std::map< std::string, ReadNC::AttData >& globalAtts = _readNC->globalAtts;
126  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
127  DebugOutput& dbgOut = _readNC->dbgOut;
128  int& partMethod = _readNC->partMethod;
129  ScdInterface* scdi = _readNC->scdi;
130 
131  ErrorCode rval;
132  std::string tag_name;
133 
134  // <__NUM_DIMS>
135  Tag numDimsTag = 0;
136  tag_name = "__NUM_DIMS";
137  int numDims = dimNames.size();
138  rval = mbImpl->tag_get_handle( tag_name.c_str(), 1, MB_TYPE_INTEGER, numDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
139  rval = mbImpl->tag_set_data( numDimsTag, &_fileSet, 1, &numDims );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
140  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
141 
142  // <__NUM_VARS>
143  Tag numVarsTag = 0;
144  tag_name = "__NUM_VARS";
145  int numVars = varInfo.size();
146  rval = mbImpl->tag_get_handle( tag_name.c_str(), 1, MB_TYPE_INTEGER, numVarsTag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
147  rval = mbImpl->tag_set_data( numVarsTag, &_fileSet, 1, &numVars );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
148  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
149 
150  // <__DIM_NAMES>
151  Tag dimNamesTag = 0;
152  tag_name = "__DIM_NAMES";
153  std::string dimnames;
154  unsigned int dimNamesSz = dimNames.size();
155  for( unsigned int i = 0; i != dimNamesSz; i++ )
156  {
157  dimnames.append( dimNames[i] );
158  dimnames.push_back( '\0' );
159  }
160  int dimnamesSz = dimnames.size();
161  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, dimNamesTag,
162  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
163  const void* ptr = dimnames.c_str();
164  rval = mbImpl->tag_set_by_ptr( dimNamesTag, &_fileSet, 1, &ptr, &dimnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
165  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
166 
167  // <__DIM_LENS>
168  Tag dimLensTag = 0;
169  tag_name = "__DIM_LENS";
170  int dimLensSz = dimLens.size();
171  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_INTEGER, dimLensTag,
172  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
173  ptr = &( dimLens[0] );
174  rval = mbImpl->tag_set_by_ptr( dimLensTag, &_fileSet, 1, &ptr, &dimLensSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
175  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
176 
177  // <__VAR_NAMES>
178  Tag varNamesTag = 0;
179  tag_name = "__VAR_NAMES";
180  std::string varnames;
181  std::map< std::string, ReadNC::VarData >::iterator mapIter;
182  for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
183  {
184  varnames.append( mapIter->first );
185  varnames.push_back( '\0' );
186  }
187  int varnamesSz = varnames.size();
188  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varNamesTag,
189  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
190  ptr = varnames.c_str();
191  rval = mbImpl->tag_set_by_ptr( varNamesTag, &_fileSet, 1, &ptr, &varnamesSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
192  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
193 
194  // __<dim_name>_LOC_MINMAX (for time)
195  for( unsigned int i = 0; i != dimNamesSz; i++ )
196  {
197  if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
198  {
199  // some files might have Time dimension as 0, it will not appear in any
200  // variables, so skip it
201  if( nTimeSteps == 0 ) continue;
202  std::stringstream ss_tag_name;
203  ss_tag_name << "__" << dimNames[i] << "_LOC_MINMAX";
204  tag_name = ss_tag_name.str();
205  Tag tagh = 0;
206  std::vector< int > val( 2, 0 );
207  val[0] = 0;
208  val[1] = nTimeSteps - 1;
209  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 );
210  rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
211  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
212  }
213  }
214 
215  // __<dim_name>_LOC_VALS (for time)
216  for( unsigned int i = 0; i != dimNamesSz; i++ )
217  {
218  if( dimNames[i] == "time" || dimNames[i] == "Time" || dimNames[i] == "t" )
219  {
220  std::vector< int > val;
221  if( !tstep_nums.empty() )
222  val = tstep_nums;
223  else
224  {
225  // some files might have Time dimension as 0, it will not appear in any
226  // variables, so skip it
227  if( tVals.empty() ) continue;
228  val.resize( tVals.size() );
229  for( unsigned int j = 0; j != tVals.size(); j++ )
230  val[j] = j;
231  }
232  Tag tagh = 0;
233  std::stringstream ss_tag_name;
234  ss_tag_name << "__" << dimNames[i] << "_LOC_VALS";
235  tag_name = ss_tag_name.str();
236  rval = mbImpl->tag_get_handle( tag_name.c_str(), val.size(), MB_TYPE_INTEGER, tagh,
237  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
238  rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
239  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
240  }
241  }
242 
243  // __<var_name>_DIMS
244  for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
245  {
246  Tag varNamesDimsTag = 0;
247  std::stringstream ss_tag_name;
248  ss_tag_name << "__" << mapIter->first << "_DIMS";
249  tag_name = ss_tag_name.str();
250  unsigned int varDimSz = varInfo[mapIter->first].varDims.size();
251  if( varDimSz == 0 ) continue;
252  std::vector< Tag > varDimTags( varDimSz );
253  for( unsigned int i = 0; i != varDimSz; i++ )
254  {
255  Tag tmptag = 0;
256  std::string tmptagname = dimNames[varInfo[mapIter->first].varDims[i]];
257  rval = mbImpl->tag_get_handle( tmptagname.c_str(), 0, MB_TYPE_OPAQUE, tmptag, MB_TAG_ANY );MB_CHK_SET_ERR( rval, "Trouble getting tag " << tmptagname );
258  varDimTags[i] = tmptag;
259  }
260  // rval = mbImpl->tag_get_handle(tag_name.c_str(), varDimSz, MB_TYPE_HANDLE,
261  // varNamesDimsTag, MB_TAG_SPARSE | MB_TAG_CREAT); We should not use MB_TYPE_HANDLE for Tag
262  // here. Tag is a pointer, which is 4 bytes on 32 bit machines and 8 bytes on 64 bit
263  // machines. Normally, entity handle is 8 bytes on 64 bit machines, but it can also be
264  // configured to 4 bytes.
265  rval = mbImpl->tag_get_handle( tag_name.c_str(), varDimSz * sizeof( Tag ), MB_TYPE_OPAQUE, varNamesDimsTag,
266  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
267  rval = mbImpl->tag_set_data( varNamesDimsTag, &_fileSet, 1, &( varDimTags[0] ) );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
268  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
269  }
270 
271  // <PARTITION_METHOD>
272  Tag part_tag = scdi->part_method_tag();
273  if( !part_tag ) MB_SET_ERR( MB_FAILURE, "Trouble getting PARTITION_METHOD tag" );
274  rval = mbImpl->tag_set_data( part_tag, &_fileSet, 1, &partMethod );MB_CHK_SET_ERR( rval, "Trouble setting data to PARTITION_METHOD tag" );
275  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
276 
277  // <__GLOBAL_ATTRIBS>
278  tag_name = "__GLOBAL_ATTRIBS";
279  Tag globalAttTag = 0;
280  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, globalAttTag,
281  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
282  std::string gattVal;
283  std::vector< int > gattLen;
284  rval = create_attrib_string( globalAtts, gattVal, gattLen );MB_CHK_SET_ERR( rval, "Trouble creating global attribute string" );
285  const void* gattptr = gattVal.c_str();
286  int globalAttSz = gattVal.size();
287  rval = mbImpl->tag_set_by_ptr( globalAttTag, &_fileSet, 1, &gattptr, &globalAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
288  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
289 
290  // <__GLOBAL_ATTRIBS_LEN>
291  tag_name = "__GLOBAL_ATTRIBS_LEN";
292  Tag globalAttLenTag = 0;
293  if( gattLen.size() == 0 ) gattLen.push_back( 0 );
294  rval = mbImpl->tag_get_handle( tag_name.c_str(), gattLen.size(), MB_TYPE_INTEGER, globalAttLenTag,
295  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
296  rval = mbImpl->tag_set_data( globalAttLenTag, &_fileSet, 1, &gattLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
297  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
298 
299  // __<var_name>_ATTRIBS and __<var_name>_ATTRIBS_LEN
300  for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
301  {
302  std::stringstream ssTagName;
303  ssTagName << "__" << mapIter->first << "_ATTRIBS";
304  tag_name = ssTagName.str();
305  Tag varAttTag = 0;
306  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, varAttTag,
307  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
308 
309  std::string varAttVal;
310  std::vector< int > varAttLen;
311  if( mapIter->second.numAtts < 1 )
312  {
313  if( dummyVarNames.find( mapIter->first ) != dummyVarNames.end() )
314  {
315  // This variable is a dummy coordinate variable
316  varAttVal = "DUMMY_VAR";
317  }
318  else
319  {
320  // This variable has no attributes
321  varAttVal = "NO_ATTRIBS";
322  }
323  }
324  else
325  {
326  rval = create_attrib_string( mapIter->second.varAtts, varAttVal, varAttLen );MB_CHK_SET_ERR( rval, "Trouble creating attribute string for variable " << mapIter->first );
327  }
328  const void* varAttPtr = varAttVal.c_str();
329  int varAttSz = varAttVal.size();
330  if( 0 == varAttSz ) varAttSz = 1;
331  rval = mbImpl->tag_set_by_ptr( varAttTag, &_fileSet, 1, &varAttPtr, &varAttSz );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
332  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
333 
334  ssTagName << "_LEN";
335  tag_name = ssTagName.str();
336  Tag varAttLenTag = 0;
337  if( 0 == varAttLen.size() ) varAttLen.push_back( 0 );
338  rval = mbImpl->tag_get_handle( tag_name.c_str(), varAttLen.size(), MB_TYPE_INTEGER, varAttLenTag,
339  MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
340  rval = mbImpl->tag_set_data( varAttLenTag, &_fileSet, 1, &varAttLen[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
341  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
342  }
343 
344  // <__VAR_NAMES_LOCATIONS>
345  tag_name = "__VAR_NAMES_LOCATIONS";
346  Tag varNamesLocsTag = 0;
347  std::vector< int > varNamesLocs( varInfo.size() );
348  rval = mbImpl->tag_get_handle( tag_name.c_str(), varNamesLocs.size(), MB_TYPE_INTEGER, varNamesLocsTag,
349  MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
350  for( mapIter = varInfo.begin(); mapIter != varInfo.end(); ++mapIter )
351  {
352  varNamesLocs[std::distance( varInfo.begin(), mapIter )] = mapIter->second.entLoc;
353  }
354  rval = mbImpl->tag_set_data( varNamesLocsTag, &_fileSet, 1, &varNamesLocs[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
355  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
356 
357  // <__MESH_TYPE>
358  Tag meshTypeTag = 0;
359  tag_name = "__MESH_TYPE";
360  std::string meshTypeName = get_mesh_type_name();
361 
362  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_OPAQUE, meshTypeTag,
363  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
364  ptr = meshTypeName.c_str();
365  int leng = meshTypeName.size();
366  rval = mbImpl->tag_set_by_ptr( meshTypeTag, &_fileSet, 1, &ptr, &leng );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
367  dbgOut.tprintf( 2, "Conventional tag %s created\n", tag_name.c_str() );
368 
369  return MB_SUCCESS;
370 }
371 
372 ErrorCode NCHelper::update_time_tag_vals()
373 {
374  Interface*& mbImpl = _readNC->mbImpl;
375  std::vector< std::string >& dimNames = _readNC->dimNames;
376 
377  ErrorCode rval;
378 
379  // The time tag might be a dummy one (e.g. 'Time' for MPAS)
380  std::string time_tag_name = dimNames[tDim];
381  if( dummyVarNames.find( time_tag_name ) != dummyVarNames.end() ) return MB_SUCCESS;
382 
383  Tag time_tag = 0;
384  const void* data = NULL;
385  int time_tag_size = 0;
386  rval = mbImpl->tag_get_handle( time_tag_name.c_str(), 0, MB_TYPE_DOUBLE, time_tag, MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble getting tag " << time_tag_name );
387  rval = mbImpl->tag_get_by_ptr( time_tag, &_fileSet, 1, &data, &time_tag_size );MB_CHK_SET_ERR( rval, "Trouble getting data of tag " << time_tag_name );
388  const double* time_tag_vals = static_cast< const double* >( data );
389 
390  // Merge tVals (read from current file) to existing time tag
391  // Assume that time_tag_vals and tVals are both sorted
392  std::vector< double > merged_time_vals;
393  merged_time_vals.reserve( time_tag_size + nTimeSteps );
394  int i = 0;
395  int j = 0;
396 
397  // Merge time values from time_tag_vals and tVals
398  while( i < time_tag_size && j < nTimeSteps )
399  {
400  if( time_tag_vals[i] < tVals[j] )
401  merged_time_vals.push_back( time_tag_vals[i++] );
402  else
403  merged_time_vals.push_back( tVals[j++] );
404  }
405 
406  // Append remaining time values of time_tag_vals (if any)
407  while( i < time_tag_size )
408  merged_time_vals.push_back( time_tag_vals[i++] );
409 
410  // Append remaining time values of tVals (if any)
411  while( j < nTimeSteps )
412  merged_time_vals.push_back( tVals[j++] );
413 
414  data = &merged_time_vals[0];
415  time_tag_size = merged_time_vals.size();
416  rval = mbImpl->tag_set_by_ptr( time_tag, &_fileSet, 1, &data, &time_tag_size );MB_CHK_SET_ERR( rval, "Trouble setting data to tag " << time_tag_name );
417 
418  return MB_SUCCESS;
419 }
420 
421 ErrorCode NCHelper::read_variables_setup( std::vector< std::string >& var_names,
422  std::vector< int >& tstep_nums,
423  std::vector< ReadNC::VarData >& vdatas,
424  std::vector< ReadNC::VarData >& vsetdatas )
425 {
426  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
427  std::vector< std::string >& dimNames = _readNC->dimNames;
428 
429  std::map< std::string, ReadNC::VarData >::iterator mit;
430 
431  // If empty read them all (except ignored variables)
432  if( var_names.empty() )
433  {
434  for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
435  {
436  ReadNC::VarData vd = ( *mit ).second;
437 
438  // If read all variables at once, skip ignored ones
439  if( ignoredVarNames.find( vd.varName ) != ignoredVarNames.end() ) continue;
440 
441  // Coordinate variables (include dummy ones) were read to the file set by default
442  if( std::find( dimNames.begin(), dimNames.end(), vd.varName ) != dimNames.end() ) continue;
443 
444  if( vd.entLoc == ReadNC::ENTLOCSET )
445  vsetdatas.push_back( vd );
446  else
447  vdatas.push_back( vd );
448  }
449  }
450  else
451  {
452  // Read specified variables (might include ignored ones)
453  for( unsigned int i = 0; i < var_names.size(); i++ )
454  {
455  mit = varInfo.find( var_names[i] );
456  if( mit != varInfo.end() )
457  {
458  ReadNC::VarData vd = ( *mit ).second;
459 
460  // Upon creation of dummy coordinate variables, tag values have already been set
461  if( dummyVarNames.find( vd.varName ) != dummyVarNames.end() ) continue;
462 
463  if( vd.entLoc == ReadNC::ENTLOCSET )
464  vsetdatas.push_back( vd );
465  else
466  vdatas.push_back( vd );
467  }
468  else
469  {
470  MB_SET_ERR( MB_FAILURE, "Couldn't find specified variable " << var_names[i] );
471  }
472  }
473  }
474 
475  if( tstep_nums.empty() && nTimeSteps > 0 )
476  {
477  // No timesteps input, get them all
478  for( int i = 0; i < nTimeSteps; i++ )
479  tstep_nums.push_back( i );
480  }
481 
482  if( !tstep_nums.empty() )
483  {
484  for( unsigned int i = 0; i < vdatas.size(); i++ )
485  {
486  vdatas[i].varTags.resize( tstep_nums.size(), 0 );
487  vdatas[i].varDatas.resize( tstep_nums.size() );
488  // NC reader assumes that non-set variables always have timesteps
489  assert( std::find( vdatas[i].varDims.begin(), vdatas[i].varDims.end(), tDim ) != vdatas[i].varDims.end() );
490  vdatas[i].has_tsteps = true;
491  }
492 
493  for( unsigned int i = 0; i < vsetdatas.size(); i++ )
494  {
495  if( ( std::find( vsetdatas[i].varDims.begin(), vsetdatas[i].varDims.end(), tDim ) !=
496  vsetdatas[i].varDims.end() ) &&
497  ( vsetdatas[i].varName != dimNames[tDim] ) )
498  {
499  // Set variables with timesteps: e.g. xtime(Time) or xtime(Time, StrLen)
500  vsetdatas[i].varTags.resize( tstep_nums.size(), 0 );
501  vsetdatas[i].varDatas.resize( tstep_nums.size() );
502  vsetdatas[i].has_tsteps = true;
503  }
504  else
505  {
506  // Set variables without timesteps: no time dimension, or time itself
507  vsetdatas[i].varTags.resize( 1, 0 );
508  vsetdatas[i].varDatas.resize( 1 );
509  vsetdatas[i].has_tsteps = false;
510  }
511  }
512  }
513 
514  return MB_SUCCESS;
515 }
516 
517 ErrorCode NCHelper::read_variables_to_set( std::vector< ReadNC::VarData >& vdatas, std::vector< int >& tstep_nums )
518 {
519  Interface*& mbImpl = _readNC->mbImpl;
520  DebugOutput& dbgOut = _readNC->dbgOut;
521 
522  ErrorCode rval = read_variables_to_set_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read set variables" );
523 
524  // Finally, read into that space
525  int success;
526  for( unsigned int i = 0; i < vdatas.size(); i++ )
527  {
528  // Note, for set variables without timesteps, loop one time and then break
529  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
530  {
531  void* data = vdatas[i].varDatas[t];
532 
533  // Set variables with timesteps, e.g. xtime(Time) or xtime(Time, StrLen)
534  if( vdatas[i].has_tsteps )
535  {
536  // Set readStart for each timestep along time dimension
537  vdatas[i].readStarts[0] = tstep_nums[t];
538  }
539 
540  switch( vdatas[i].varDataType )
541  {
542  case NC_BYTE:
543  case NC_CHAR:
544  success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
545  &vdatas[i].readCounts[0], (char*)data );
546  if( success )
547  MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
548  break;
549  case NC_SHORT:
550  case NC_INT:
551  success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
552  &vdatas[i].readCounts[0], (int*)data );
553  if( success )
554  MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
555  break;
556  case NC_INT64:
557  success = NCFUNCAG( _vara_long )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
558  &vdatas[i].readCounts[0], (long*)data );
559  if( success )
560  MB_SET_ERR( MB_FAILURE, "Failed to read long data for variable " << vdatas[i].varName );
561  break;
562  case NC_FLOAT:
563  case NC_DOUBLE:
564  success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
565  &vdatas[i].readCounts[0], (double*)data );
566  if( success )
567  MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
568  break;
569  default:
570  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
571  }
572 
573  dbgOut.tprintf( 2, "Setting data for variable %s, time step %d\n", vdatas[i].varName.c_str(),
574  tstep_nums[t] );
575  rval = mbImpl->tag_set_by_ptr( vdatas[i].varTags[t], &_fileSet, 1, &data, &vdatas[i].sz );MB_CHK_SET_ERR( rval, "Trouble setting tag data for variable " << vdatas[i].varName );
576 
577  // Memory pointed by pointer data can be deleted, as tag_set_by_ptr() has already copied
578  // the tag values
579  switch( vdatas[i].varDataType )
580  {
581  case NC_BYTE:
582  case NC_CHAR:
583  delete[](char*)data;
584  break;
585  case NC_SHORT:
586  case NC_INT:
587  delete[](int*)data;
588  break;
589  case NC_INT64:
590  delete[](long*)data;
591  break;
592  case NC_FLOAT:
593  case NC_DOUBLE:
594  delete[](double*)data;
595  break;
596  default:
597  break;
598  }
599  vdatas[i].varDatas[t] = NULL;
600 
601  // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
602  // StrLen)
603  if( !vdatas[i].has_tsteps ) break;
604  }
605  }
606 
607  // Debug output, if requested
608  if( 1 == dbgOut.get_verbosity() )
609  {
610  dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
611  for( unsigned int i = 1; i < vdatas.size(); i++ )
612  dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
613  dbgOut.tprintf( 1, "\n" );
614  }
615 
616  return rval;
617 }
618 
619 ErrorCode NCHelper::read_coordinate( const char* var_name, int lmin, int lmax, std::vector< double >& cvals )
620 {
621  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
622  std::map< std::string, ReadNC::VarData >::iterator vmit = varInfo.find( var_name );
623  if( varInfo.end() == vmit ) MB_SET_ERR( MB_FAILURE, "Couldn't find variable " << var_name );
624 
625  assert( lmin >= 0 && lmax >= lmin );
626  NCDF_SIZE tstart = lmin;
627  NCDF_SIZE tcount = lmax - lmin + 1;
628  NCDF_DIFF dum_stride = 1;
629  int success;
630 
631  // Check size
632  if( (std::size_t)tcount != cvals.size() ) cvals.resize( tcount );
633 
634  // Check to make sure it's a float or double
635  switch( ( *vmit ).second.varDataType )
636  {
637  case NC_FLOAT:
638  case NC_DOUBLE:
639  // Read float as double
640  success =
641  NCFUNCAG( _vars_double )( _fileId, ( *vmit ).second.varId, &tstart, &tcount, &dum_stride, &cvals[0] );
642  if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << var_name );
643  break;
644  default:
645  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_name );
646  }
647 
648  return MB_SUCCESS;
649 }
650 
651 ErrorCode NCHelper::get_tag_to_set( ReadNC::VarData& var_data, int tstep_num, Tag& tagh )
652 {
653  Interface*& mbImpl = _readNC->mbImpl;
654  DebugOutput& dbgOut = _readNC->dbgOut;
655  int& tStepBase = _readNC->tStepBase;
656 
657  if( tStepBase > 0 ) tstep_num += tStepBase;
658 
659  std::ostringstream tag_name;
660  if( var_data.has_tsteps )
661  tag_name << var_data.varName << tstep_num;
662  else
663  tag_name << var_data.varName;
664 
665  ErrorCode rval = MB_SUCCESS;
666  tagh = 0;
667  switch( var_data.varDataType )
668  {
669  case NC_BYTE:
670  case NC_CHAR:
671  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_OPAQUE, tagh,
672  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
673  break;
674  case NC_SHORT:
675  case NC_INT:
676  case NC_INT64: // this is a big stretch; we should introduce LONG tag
677  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_INTEGER, tagh,
678  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
679  break;
680  case NC_FLOAT:
681  case NC_DOUBLE:
682  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), 0, MB_TYPE_DOUBLE, tagh,
683  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
684  break;
685  default:
686  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
687  }
688 
689  dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
690 
691  return rval;
692 }
693 
694 ErrorCode NCHelper::get_tag_to_nonset( ReadNC::VarData& var_data, int tstep_num, Tag& tagh, int num_lev )
695 {
696  Interface*& mbImpl = _readNC->mbImpl;
697  DebugOutput& dbgOut = _readNC->dbgOut;
698  int& tStepBase = _readNC->tStepBase;
699 
700  if( tStepBase > 0 ) tstep_num += tStepBase;
701 
702  std::ostringstream tag_name;
703  tag_name << var_data.varName << tstep_num;
704 
705  ErrorCode rval = MB_SUCCESS;
706  tagh = 0;
707  switch( var_data.varDataType )
708  {
709  case NC_BYTE:
710  case NC_CHAR:
711  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_OPAQUE, tagh,
712  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
713  break;
714  case NC_SHORT:
715  case NC_INT:
716  case NC_INT64:
717  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_INTEGER, tagh,
718  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
719  break;
720  case NC_FLOAT:
721  case NC_DOUBLE:
722  rval = mbImpl->tag_get_handle( tag_name.str().c_str(), num_lev, MB_TYPE_DOUBLE, tagh,
723  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating tag " << tag_name.str() );
724  break;
725  default:
726  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << var_data.varName );
727  }
728 
729  dbgOut.tprintf( 2, "Tag %s created\n", tag_name.str().c_str() );
730 
731  return rval;
732 }
733 
734 ErrorCode NCHelper::create_attrib_string( const std::map< std::string, ReadNC::AttData >& attMap,
735  std::string& attVal,
736  std::vector< int >& attLen )
737 {
738  int success;
739  std::stringstream ssAtt;
740  unsigned int sz = 0;
741  std::map< std::string, ReadNC::AttData >::const_iterator attIt = attMap.begin();
742  for( ; attIt != attMap.end(); ++attIt )
743  {
744  ssAtt << attIt->second.attName;
745  ssAtt << '\0';
746  void* attData = NULL;
747  switch( attIt->second.attDataType )
748  {
749  case NC_BYTE:
750  case NC_CHAR:
751  sz = attIt->second.attLen;
752  attData = (char*)malloc( sz );
753  success = NCFUNC( get_att_text )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
754  (char*)attData );
755  if( success )
756  MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for attribute " << attIt->second.attName );
757  ssAtt << "char;";
758  break;
759  case NC_SHORT:
760  sz = attIt->second.attLen * sizeof( short );
761  attData = (short*)malloc( sz );
762  success = NCFUNC( get_att_short )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
763  (short*)attData );
764  if( success )
765  MB_SET_ERR( MB_FAILURE, "Failed to read short data for attribute " << attIt->second.attName );
766  ssAtt << "short;";
767  break;
768  case NC_INT:
769  sz = attIt->second.attLen * sizeof( int );
770  attData = (int*)malloc( sz );
771  success = NCFUNC( get_att_int )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
772  (int*)attData );
773  if( success )
774  MB_SET_ERR( MB_FAILURE, "Failed to read int data for attribute " << attIt->second.attName );
775  ssAtt << "int;";
776  break;
777  case NC_INT64: // be careful here
778  sz = attIt->second.attLen * sizeof( long );
779  attData = (long*)malloc( sz );
780  success = NCFUNC( get_att_long )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
781  (long*)attData );
782  if( success )
783  MB_SET_ERR( MB_FAILURE, "Failed to read int data for attribute " << attIt->second.attName );
784  ssAtt << "long;";
785  break;
786  case NC_FLOAT:
787  sz = attIt->second.attLen * sizeof( float );
788  attData = (float*)malloc( sz );
789  success = NCFUNC( get_att_float )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
790  (float*)attData );
791  if( success )
792  MB_SET_ERR( MB_FAILURE, "Failed to read float data for attribute " << attIt->second.attName );
793  ssAtt << "float;";
794  break;
795  case NC_DOUBLE:
796  sz = attIt->second.attLen * sizeof( double );
797  attData = (double*)malloc( sz );
798  success = NCFUNC( get_att_double )( _fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
799  (double*)attData );
800  if( success )
801  MB_SET_ERR( MB_FAILURE, "Failed to read double data for attribute " << attIt->second.attName );
802  ssAtt << "double;";
803  break;
804  default:
805  MB_SET_ERR( MB_FAILURE, "Unexpected data type for attribute " << attIt->second.attName );
806  }
807  char* tmpc = (char*)attData;
808  for( unsigned int counter = 0; counter != sz; ++counter )
809  ssAtt << tmpc[counter];
810  free( attData );
811  ssAtt << ';';
812  attLen.push_back( ssAtt.str().size() - 1 );
813  }
814  attVal = ssAtt.str();
815 
816  return MB_SUCCESS;
817 }
818 
819 ErrorCode NCHelper::create_dummy_variables()
820 {
821  Interface*& mbImpl = _readNC->mbImpl;
822  std::vector< std::string >& dimNames = _readNC->dimNames;
823  std::vector< int >& dimLens = _readNC->dimLens;
824  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
825  DebugOutput& dbgOut = _readNC->dbgOut;
826 
827  // Hack: look at all dimensions, and see if we have one that does not appear in the list of
828  // varInfo names Right now, candidates are from unstructured meshes, such as ncol (HOMME) and
829  // nCells (MPAS) For each of them, create a dummy coordinate variable with a sparse tag to store
830  // the dimension length
831  for( unsigned int i = 0; i < dimNames.size(); i++ )
832  {
833  // If there is a variable with this dimension name, skip
834  if( varInfo.find( dimNames[i] ) != varInfo.end() ) continue;
835 
836  // Create a dummy coordinate variable
837  int sizeTotalVar = varInfo.size();
838  std::string var_name( dimNames[i] );
839  ReadNC::VarData& data = varInfo[var_name];
840  data.varName = var_name;
841  data.varId = sizeTotalVar;
842  data.varTags.resize( 1, 0 );
843  data.varDataType = NC_INT;
844  data.varDims.resize( 1 );
845  data.varDims[0] = (int)i;
846  data.numAtts = 0;
847  data.entLoc = ReadNC::ENTLOCSET;
848  dummyVarNames.insert( var_name );
849  dbgOut.tprintf( 2, "Dummy coordinate variable created for dimension %s\n", var_name.c_str() );
850 
851  // Create a corresponding sparse tag
852  Tag tagh;
853  ErrorCode rval = mbImpl->tag_get_handle( var_name.c_str(), 0, MB_TYPE_INTEGER, tagh,
854  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating tag for dummy coordinate variable " << var_name );
855 
856  // Tag value is the dimension length
857  const void* ptr = &dimLens[i];
858  // Tag size is 1
859  int size = 1;
860  rval = mbImpl->tag_set_by_ptr( tagh, &_fileSet, 1, &ptr, &size );MB_CHK_SET_ERR( rval, "Trouble setting tag data for dummy coordinate variable " << var_name );
861 
862  dbgOut.tprintf( 2, "Sparse tag created for dimension %s\n", var_name.c_str() );
863  }
864 
865  return MB_SUCCESS;
866 }
867 
868 ErrorCode NCHelper::read_variables_to_set_allocate( std::vector< ReadNC::VarData >& vdatas,
869  std::vector< int >& tstep_nums )
870 {
871  std::vector< int >& dimLens = _readNC->dimLens;
872  DebugOutput& dbgOut = _readNC->dbgOut;
873 
874  ErrorCode rval = MB_SUCCESS;
875 
876  for( unsigned int i = 0; i < vdatas.size(); i++ )
877  {
878  // Set up readStarts and readCounts
879  if( vdatas[i].has_tsteps )
880  {
881  // First: time
882  vdatas[i].readStarts.push_back( 0 ); // This value is timestep dependent, will be set later
883  vdatas[i].readCounts.push_back( 1 );
884 
885  // Next: other dimensions
886  for( unsigned int idx = 1; idx != vdatas[i].varDims.size(); idx++ )
887  {
888  vdatas[i].readStarts.push_back( 0 );
889  vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
890  }
891  }
892  else
893  {
894  if( vdatas[i].varDims.empty() )
895  {
896  // Scalar variable
897  vdatas[i].readStarts.push_back( 0 );
898  vdatas[i].readCounts.push_back( 1 );
899  }
900  else
901  {
902  for( unsigned int idx = 0; idx != vdatas[i].varDims.size(); idx++ )
903  {
904  vdatas[i].readStarts.push_back( 0 );
905  vdatas[i].readCounts.push_back( dimLens[vdatas[i].varDims[idx]] );
906  }
907  }
908  }
909 
910  // Get variable size
911  vdatas[i].sz = 1;
912  for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
913  vdatas[i].sz *= vdatas[i].readCounts[idx];
914 
915  // Note, for set variables without timesteps, loop one time and then break
916  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
917  {
918  dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
919 
920  if( tstep_nums[t] >= dimLens[tDim] )
921  {
922  MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
923  }
924 
925  // Get the tag to read into
926  if( !vdatas[i].varTags[t] )
927  {
928  rval = get_tag_to_set( vdatas[i], tstep_nums[t], vdatas[i].varTags[t] );MB_CHK_SET_ERR( rval, "Trouble getting tag to set variable " << vdatas[i].varName );
929  }
930 
931  switch( vdatas[i].varDataType )
932  {
933  case NC_BYTE:
934  case NC_CHAR:
935  vdatas[i].varDatas[t] = new char[vdatas[i].sz];
936  break;
937  case NC_SHORT:
938  case NC_INT:
939  vdatas[i].varDatas[t] = new int[vdatas[i].sz];
940  break;
941  case NC_INT64:
942  vdatas[i].varDatas[t] = new long[vdatas[i].sz];
943  break;
944  case NC_FLOAT:
945  case NC_DOUBLE:
946  vdatas[i].varDatas[t] = new double[vdatas[i].sz];
947  break;
948  default:
949  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
950  }
951 
952  // Loop continues only for set variables with timesteps, e.g. xtime(Time) or xtime(Time,
953  // StrLen)
954  if( !vdatas[i].has_tsteps ) break;
955  }
956  }
957 
958  return rval;
959 }
960 
961 ErrorCode ScdNCHelper::check_existing_mesh()
962 {
963  Interface*& mbImpl = _readNC->mbImpl;
964 
965  // Get the number of vertices
966  int num_verts;
967  ErrorCode rval = mbImpl->get_number_entities_by_dimension( _fileSet, 0, num_verts );MB_CHK_SET_ERR( rval, "Trouble getting number of vertices" );
968 
969  /*
970  // Check against parameters
971  // When ghosting is used, this check might fail (to be updated later)
972  if (num_verts > 0) {
973  int expected_verts = (lDims[3] - lDims[0] + 1) * (lDims[4] - lDims[1] + 1) * (-1 == lDims[2] ?
974  1 : lDims[5] - lDims[2] + 1); if (num_verts != expected_verts) { MB_SET_ERR(MB_FAILURE, "Number
975  of vertices doesn't match");
976  }
977  }
978  */
979 
980  // Check the number of elements too
981  int num_elems;
982  rval = mbImpl->get_number_entities_by_dimension( _fileSet, ( -1 == lCDims[2] ? 2 : 3 ), num_elems );MB_CHK_SET_ERR( rval, "Trouble getting number of elements" );
983 
984  /*
985  // Check against parameters
986  // When ghosting is used, this check might fail (to be updated later)
987  if (num_elems > 0) {
988  int expected_elems = (lCDims[3] - lCDims[0] + 1) * (lCDims[4] - lCDims[1] + 1) * (-1 ==
989  lCDims[2] ? 1 : (lCDims[5] - lCDims[2] + 1)); if (num_elems != expected_elems) {
990  MB_SET_ERR(MB_FAILURE, "Number of elements doesn't match");
991  }
992  }
993  */
994 
995  return MB_SUCCESS;
996 }
997 
998 ErrorCode ScdNCHelper::create_mesh( Range& faces )
999 {
1000  Interface*& mbImpl = _readNC->mbImpl;
1001  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
1002  const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
1003  DebugOutput& dbgOut = _readNC->dbgOut;
1004  ScdInterface* scdi = _readNC->scdi;
1005  ScdParData& parData = _readNC->parData;
1006 
1007  Range tmp_range;
1008  ScdBox* scd_box;
1009 
1010  ErrorCode rval =
1011  scdi->construct_box( HomCoord( lDims[0], lDims[1], lDims[2], 1 ), HomCoord( lDims[3], lDims[4], lDims[5], 1 ),
1012  NULL, 0, scd_box, locallyPeriodic, &parData, true );MB_CHK_SET_ERR( rval, "Trouble creating scd vertex sequence" );
1013 
1014  // Add verts to tmp_range first, so we can duplicate global ids in vertex ids
1015  tmp_range.insert( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
1016 
1017  if( mpFileIdTag )
1018  {
1019  int count;
1020  void* data;
1021  rval = mbImpl->tag_iterate( *mpFileIdTag, tmp_range.begin(), tmp_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file ID tag on local vertices" );
1022  assert( count == scd_box->num_vertices() );
1023  int* fid_data = (int*)data;
1024  rval = mbImpl->tag_iterate( mGlobalIdTag, tmp_range.begin(), tmp_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global ID tag on local vertices" );
1025  assert( count == scd_box->num_vertices() );
1026  int* gid_data = (int*)data;
1027  for( int i = 0; i < count; i++ )
1028  fid_data[i] = gid_data[i];
1029  }
1030 
1031  // Then add box set and elements to the range, then to the file set
1032  tmp_range.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
1033  tmp_range.insert( scd_box->box_set() );
1034  rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Couldn't add new vertices to current file set" );
1035 
1036  dbgOut.tprintf( 1, "scdbox %d quads, %d vertices\n", scd_box->num_elements(), scd_box->num_vertices() );
1037 
1038  // Set the vertex coordinates
1039  double *xc, *yc, *zc;
1040  rval = scd_box->get_coordinate_arrays( xc, yc, zc );MB_CHK_SET_ERR( rval, "Couldn't get vertex coordinate arrays" );
1041 
1042  int i, j, k, il, jl, kl;
1043  int dil = lDims[3] - lDims[0] + 1;
1044  int djl = lDims[4] - lDims[1] + 1;
1045  assert( dil == (int)ilVals.size() && djl == (int)jlVals.size() &&
1046  ( -1 == lDims[2] || lDims[5] - lDims[2] + 1 == (int)levVals.size() ) );
1047 
1048  for( kl = lDims[2]; kl <= lDims[5]; kl++ )
1049  {
1050  k = kl - lDims[2];
1051  for( jl = lDims[1]; jl <= lDims[4]; jl++ )
1052  {
1053  j = jl - lDims[1];
1054  for( il = lDims[0]; il <= lDims[3]; il++ )
1055  {
1056  i = il - lDims[0];
1057  unsigned int pos = i + j * dil + k * dil * djl;
1058  xc[pos] = ilVals[i];
1059  yc[pos] = jlVals[j];
1060  zc[pos] = ( -1 == lDims[2] ? 0.0 : levVals[k] );
1061  }
1062  }
1063  }
1064 
1065 #ifndef NDEBUG
1066  int num_verts =
1067  ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) * ( -1 == lDims[2] ? 1 : lDims[5] - lDims[2] + 1 );
1068  std::vector< int > gids( num_verts );
1069  Range verts( scd_box->start_vertex(), scd_box->start_vertex() + scd_box->num_vertices() - 1 );
1070  rval = mbImpl->tag_get_data( mGlobalIdTag, verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
1071  int vmin = *( std::min_element( gids.begin(), gids.end() ) ),
1072  vmax = *( std::max_element( gids.begin(), gids.end() ) );
1073  dbgOut.tprintf( 1, "Vertex gids %d-%d\n", vmin, vmax );
1074 #endif
1075 
1076  // Add elements to the range passed in
1077  faces.insert( scd_box->start_element(), scd_box->start_element() + scd_box->num_elements() - 1 );
1078 
1079  if( 2 <= dbgOut.get_verbosity() )
1080  {
1081  assert( scd_box->boundary_complete() );
1082  EntityHandle dum_ent = scd_box->start_element();
1083  rval = mbImpl->list_entities( &dum_ent, 1 );MB_CHK_SET_ERR( rval, "Trouble listing first hex" );
1084 
1085  std::vector< EntityHandle > connect;
1086  rval = mbImpl->get_connectivity( &dum_ent, 1, connect );MB_CHK_SET_ERR( rval, "Trouble getting connectivity" );
1087 
1088  rval = mbImpl->list_entities( &connect[0], connect.size() );MB_CHK_SET_ERR( rval, "Trouble listing element connectivity" );
1089  }
1090 
1091  Range edges;
1092  mbImpl->get_adjacencies( faces, 1, true, edges, Interface::UNION );
1093 
1094  // Create COORDS tag for quads
1095  rval = create_quad_coordinate_tag();MB_CHK_SET_ERR( rval, "Trouble creating COORDS tag for quads" );
1096 
1097  return MB_SUCCESS;
1098 }
1099 
1100 ErrorCode ScdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
1101 {
1102  std::vector< ReadNC::VarData > vdatas;
1103  std::vector< ReadNC::VarData > vsetdatas;
1104 
1105  ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
1106 
1107  if( !vsetdatas.empty() )
1108  {
1109  rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
1110  }
1111 
1112  if( !vdatas.empty() )
1113  {
1114  rval = read_scd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1115  }
1116 
1117  return MB_SUCCESS;
1118 }
1119 
1120 ErrorCode ScdNCHelper::read_scd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
1121  std::vector< int >& tstep_nums )
1122 {
1123  Interface*& mbImpl = _readNC->mbImpl;
1124  std::vector< int >& dimLens = _readNC->dimLens;
1125  DebugOutput& dbgOut = _readNC->dbgOut;
1126 
1127  Range* range = NULL;
1128 
1129  // Get vertices
1130  Range verts;
1131  ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
1132  assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
1133 
1134  Range edges;
1135  rval = mbImpl->get_entities_by_dimension( _fileSet, 1, edges );MB_CHK_SET_ERR( rval, "Trouble getting edges in current file set" );
1136 
1137  // Get faces
1138  Range faces;
1139  rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_SET_ERR( rval, "Trouble getting faces in current file set" );
1140  assert( "Should only have a single face subrange, since they were read in one shot" && faces.psize() == 1 );
1141 
1142 #ifdef MOAB_HAVE_MPI
1143  moab::Range faces_owned;
1144  bool& isParallel = _readNC->isParallel;
1145  if( isParallel )
1146  {
1147  ParallelComm*& myPcomm = _readNC->myPcomm;
1148  rval = myPcomm->filter_pstatus( faces, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &faces_owned );MB_CHK_SET_ERR( rval, "Trouble getting owned faces in current file set" );
1149  }
1150  else
1151  faces_owned = faces; // Not running in parallel, but still with MPI
1152 #endif
1153 
1154  for( unsigned int i = 0; i < vdatas.size(); i++ )
1155  {
1156  // Support non-set variables with 4 dimensions like (time, lev, lat, lon)
1157  assert( 4 == vdatas[i].varDims.size() );
1158 
1159  // For a non-set variable, time should be the first dimension
1160  assert( tDim == vdatas[i].varDims[0] );
1161 
1162  // Set up readStarts and readCounts
1163  vdatas[i].readStarts.resize( 4 );
1164  vdatas[i].readCounts.resize( 4 );
1165 
1166  // First: time
1167  vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
1168  vdatas[i].readCounts[0] = 1;
1169 
1170  // Next: lev
1171  vdatas[i].readStarts[1] = 0;
1172  vdatas[i].readCounts[1] = vdatas[i].numLev;
1173 
1174  // Finally: lat (or slat) and lon (or slon)
1175  switch( vdatas[i].entLoc )
1176  {
1177  case ReadNC::ENTLOCVERT:
1178  // Vertices
1179  vdatas[i].readStarts[2] = lDims[1];
1180  vdatas[i].readCounts[2] = lDims[4] - lDims[1] + 1;
1181  vdatas[i].readStarts[3] = lDims[0];
1182  vdatas[i].readCounts[3] = lDims[3] - lDims[0] + 1;
1183  range = &verts;
1184  break;
1185  case ReadNC::ENTLOCNSEDGE:
1186  case ReadNC::ENTLOCEWEDGE:
1187  case ReadNC::ENTLOCEDGE:
1188  // Not implemented yet, set a global error
1189  MB_SET_GLB_ERR( MB_NOT_IMPLEMENTED, "Reading edge data is not implemented yet" );
1190  case ReadNC::ENTLOCFACE:
1191  // Faces
1192  vdatas[i].readStarts[2] = lCDims[1];
1193  vdatas[i].readCounts[2] = lCDims[4] - lCDims[1] + 1;
1194  vdatas[i].readStarts[3] = lCDims[0];
1195  vdatas[i].readCounts[3] = lCDims[3] - lCDims[0] + 1;
1196 #ifdef MOAB_HAVE_MPI
1197  range = &faces_owned;
1198 #else
1199  range = &faces;
1200 #endif
1201  break;
1202  default:
1203  MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
1204  }
1205 
1206  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
1207  {
1208  dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
1209 
1210  if( tstep_nums[t] >= dimLens[tDim] )
1211  {
1212  MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] );
1213  }
1214 
1215  // Get the tag to read into
1216  if( !vdatas[i].varTags[t] )
1217  {
1218  rval = get_tag_to_nonset( vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev );MB_CHK_SET_ERR( rval, "Trouble getting tag to non-set variable " << vdatas[i].varName );
1219  }
1220 
1221  // Get ptr to tag space
1222  void* data;
1223  int count;
1224  rval = mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate tag for non-set variable " << vdatas[i].varName );
1225  assert( (unsigned)count == range->size() );
1226  vdatas[i].varDatas[t] = data;
1227  }
1228 
1229  // Get variable size
1230  vdatas[i].sz = 1;
1231  for( std::size_t idx = 0; idx != vdatas[i].readCounts.size(); idx++ )
1232  vdatas[i].sz *= vdatas[i].readCounts[idx];
1233  }
1234 
1235  return rval;
1236 }
1237 
1238 ErrorCode ScdNCHelper::read_scd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
1239  std::vector< int >& tstep_nums )
1240 {
1241  DebugOutput& dbgOut = _readNC->dbgOut;
1242 
1243  ErrorCode rval = read_scd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
1244 
1245  // Finally, read into that space
1246  int success;
1247  for( unsigned int i = 0; i < vdatas.size(); i++ )
1248  {
1249  std::size_t sz = vdatas[i].sz;
1250 
1251  // A typical supported variable: float T(time, lev, lat, lon)
1252  // For tag values, need transpose (lev, lat, lon) to (lat, lon, lev)
1253  size_t ni = vdatas[i].readCounts[3]; // lon or slon
1254  size_t nj = vdatas[i].readCounts[2]; // lat or slat
1255  size_t nk = vdatas[i].readCounts[1]; // lev
1256 
1257  for( unsigned int t = 0; t < tstep_nums.size(); t++ )
1258  {
1259  // Tag data for this timestep
1260  void* data = vdatas[i].varDatas[t];
1261 
1262  // Set readStart for each timestep along time dimension
1263  vdatas[i].readStarts[0] = tstep_nums[t];
1264 
1265  switch( vdatas[i].varDataType )
1266  {
1267  case NC_BYTE:
1268  case NC_CHAR: {
1269  std::vector< char > tmpchardata( sz );
1270  success = NCFUNCAG( _vara_text )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1271  &vdatas[i].readCounts[0], &tmpchardata[0] );
1272  if( success )
1273  MB_SET_ERR( MB_FAILURE, "Failed to read byte/char data for variable " << vdatas[i].varName );
1274  if( vdatas[i].numLev > 1 )
1275  // Transpose (lev, lat, lon) to (lat, lon, lev)
1276  kji_to_jik( ni, nj, nk, data, &tmpchardata[0] );
1277  else
1278  {
1279  for( std::size_t idx = 0; idx != tmpchardata.size(); idx++ )
1280  ( (char*)data )[idx] = tmpchardata[idx];
1281  }
1282  break;
1283  }
1284  case NC_SHORT:
1285  case NC_INT: {
1286  std::vector< int > tmpintdata( sz );
1287  success = NCFUNCAG( _vara_int )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1288  &vdatas[i].readCounts[0], &tmpintdata[0] );
1289  if( success )
1290  MB_SET_ERR( MB_FAILURE, "Failed to read short/int data for variable " << vdatas[i].varName );
1291  if( vdatas[i].numLev > 1 )
1292  // Transpose (lev, lat, lon) to (lat, lon, lev)
1293  kji_to_jik( ni, nj, nk, data, &tmpintdata[0] );
1294  else
1295  {
1296  for( std::size_t idx = 0; idx != tmpintdata.size(); idx++ )
1297  ( (int*)data )[idx] = tmpintdata[idx];
1298  }
1299  break;
1300  }
1301  case NC_INT64: {
1302  std::vector< long > tmpintdata( sz );
1303  success = NCFUNCAG( _vara_long )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1304  &vdatas[i].readCounts[0], &tmpintdata[0] );
1305  if( success )
1306  MB_SET_ERR( MB_FAILURE, "Failed to read long data for variable " << vdatas[i].varName );
1307  if( vdatas[i].numLev > 1 )
1308  // Transpose (lev, lat, lon) to (lat, lon, lev)
1309  kji_to_jik( ni, nj, nk, data, &tmpintdata[0] );
1310  else
1311  {
1312  for( std::size_t idx = 0; idx != tmpintdata.size(); idx++ )
1313  ( (long*)data )[idx] = tmpintdata[idx];
1314  }
1315  break;
1316  }
1317  case NC_FLOAT:
1318  case NC_DOUBLE: {
1319  std::vector< double > tmpdoubledata( sz );
1320  success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &vdatas[i].readStarts[0],
1321  &vdatas[i].readCounts[0], &tmpdoubledata[0] );
1322  if( success )
1323  MB_SET_ERR( MB_FAILURE, "Failed to read float/double data for variable " << vdatas[i].varName );
1324  if( vdatas[i].numLev > 1 )
1325  // Transpose (lev, lat, lon) to (lat, lon, lev)
1326  kji_to_jik( ni, nj, nk, data, &tmpdoubledata[0] );
1327  else
1328  {
1329  for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
1330  ( (double*)data )[idx] = tmpdoubledata[idx];
1331  }
1332  break;
1333  }
1334  default:
1335  MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
1336  }
1337  }
1338  }
1339 
1340  // Debug output, if requested
1341  if( 1 == dbgOut.get_verbosity() )
1342  {
1343  dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
1344  for( unsigned int i = 1; i < vdatas.size(); i++ )
1345  dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
1346  dbgOut.tprintf( 1, "\n" );
1347  }
1348 
1349  return rval;
1350 }
1351 
1352 ErrorCode ScdNCHelper::create_quad_coordinate_tag()
1353 {
1354  Interface*& mbImpl = _readNC->mbImpl;
1355 
1356  Range ents;
1357  ErrorCode rval = mbImpl->get_entities_by_type( _fileSet, moab::MBQUAD, ents );MB_CHK_SET_ERR( rval, "Trouble getting quads" );
1358 
1359  std::size_t numOwnedEnts = 0;
1360 #ifdef MOAB_HAVE_MPI
1361  Range ents_owned;
1362  bool& isParallel = _readNC->isParallel;
1363  if( isParallel )
1364  {
1365  ParallelComm*& myPcomm = _readNC->myPcomm;
1366  rval = myPcomm->filter_pstatus( ents, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &ents_owned );MB_CHK_SET_ERR( rval, "Trouble getting owned quads" );
1367  numOwnedEnts = ents_owned.size();
1368  }
1369  else
1370  {
1371  numOwnedEnts = ents.size();
1372  ents_owned = ents;
1373  }
1374 #else
1375  numOwnedEnts = ents.size();
1376 #endif
1377 
1378  if( numOwnedEnts == 0 ) return MB_SUCCESS;
1379 
1380  assert( numOwnedEnts == ilCVals.size() * jlCVals.size() );
1381  std::vector< double > coords( numOwnedEnts * 3 );
1382  std::size_t pos = 0;
1383  for( std::size_t j = 0; j != jlCVals.size(); ++j )
1384  {
1385  for( std::size_t i = 0; i != ilCVals.size(); ++i )
1386  {
1387  pos = j * ilCVals.size() * 3 + i * 3;
1388  coords[pos] = ilCVals[i];
1389  coords[pos + 1] = jlCVals[j];
1390  coords[pos + 2] = 0.0;
1391  }
1392  }
1393  std::string tag_name = "COORDS";
1394  Tag tagh = 0;
1395  rval = mbImpl->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, tagh, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating COORDS tag" );
1396 
1397  void* data;
1398  int count;
1399 #ifdef MOAB_HAVE_MPI
1400  rval = mbImpl->tag_iterate( tagh, ents_owned.begin(), ents_owned.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate COORDS tag on quads" );
1401 #else
1402  rval = mbImpl->tag_iterate( tagh, ents.begin(), ents.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate COORDS tag on quads" );
1403 #endif
1404  assert( count == (int)numOwnedEnts );
1405  double* quad_data = (double*)data;
1406  std::copy( coords.begin(), coords.end(), quad_data );
1407 
1408  return MB_SUCCESS;
1409 }
1410 
1411 ErrorCode UcdNCHelper::read_variables( std::vector< std::string >& var_names, std::vector< int >& tstep_nums )
1412 {
1413  std::vector< ReadNC::VarData > vdatas;
1414  std::vector< ReadNC::VarData > vsetdatas;
1415 
1416  ErrorCode rval = read_variables_setup( var_names, tstep_nums, vdatas, vsetdatas );MB_CHK_SET_ERR( rval, "Trouble setting up to read variables" );
1417 
1418  if( !vsetdatas.empty() )
1419  {
1420  rval = read_variables_to_set( vsetdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to set" );
1421  }
1422 
1423  if( !vdatas.empty() )
1424  {
1425 #ifdef MOAB_HAVE_PNETCDF
1426  // With pnetcdf support, we will use async read
1427  rval = read_ucd_variables_to_nonset_async( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1428 #else
1429  // Without pnetcdf support, we will use old read
1430  rval = read_ucd_variables_to_nonset( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble reading variables to verts/edges/faces" );
1431 #endif
1432  }
1433 
1434  return MB_SUCCESS;
1435 }
1436 
1437 static double tolerance = 1.e-12;
1438 // Used to put points in an STL tree-based container
1439 bool NCHelper::Node3D::operator<( const NCHelper::Node3D& other ) const
1440 {
1441  if( coords[0] <= other.coords[0] - tolerance )
1442  return true;
1443  else if( coords[0] >= other.coords[0] + tolerance )
1444  return false;
1445  if( coords[1] <= other.coords[1] - tolerance )
1446  return true;
1447  else if( coords[1] >= other.coords[1] + tolerance )
1448  return false;
1449  if( coords[2] <= other.coords[2] - tolerance )
1450  return true;
1451  else if( coords[0] >= other.coords[0] + tolerance )
1452  return false;
1453 
1454  return false;
1455 }
1456 /*bool NCHelper::Node3D::operator == ( const NCHelper::Node3D& other ) const
1457  return ( !(*this < other) && ! ( other < *this) ) ;*/
1458 } // namespace moab