Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
NCHelperEuler.cpp
Go to the documentation of this file.
1 #include "moab/MOABConfig.h"
2 #ifdef WIN32 /* windows */
3 #define _USE_MATH_DEFINES // For M_PI
4 #endif
5 
6 #include "NCHelperEuler.hpp"
7 #include "moab/FileOptions.hpp"
8 
9 #include <cmath>
10 #include <sstream>
11 
12 #ifdef WIN32
13 #ifdef size_t
14 #undef size_t
15 #endif
16 #endif
17 
18 namespace moab
19 {
20 
21 bool NCHelperEuler::can_read_file( ReadNC* readNC, int fileId )
22 {
23  std::vector< std::string >& dimNames = readNC->dimNames;
24 
25  // If dimension names "lon" AND "lat' exist then it could be either the Eulerian Spectral grid
26  // or the FV grid
27  if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "lon" ) ) != dimNames.end() ) &&
28  ( std::find( dimNames.begin(), dimNames.end(), std::string( "lat" ) ) != dimNames.end() ) )
29  {
30  // If dimension names "lon" AND "lat" AND "slon" AND "slat" exist then it should be the FV
31  // grid
32  if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "slon" ) ) != dimNames.end() ) &&
33  ( std::find( dimNames.begin(), dimNames.end(), std::string( "slat" ) ) != dimNames.end() ) )
34  return false;
35 
36  // Make sure it is CAM grid
37  std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
38  if( attIt == readNC->globalAtts.end() ) return false;
39  unsigned int sz = attIt->second.attLen;
40  std::string att_data;
41  att_data.resize( sz + 1 );
42  att_data[sz] = '\000';
43  int success =
44  NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
45  if( success ) return false;
46  if( att_data.find( "CAM" ) == std::string::npos ) return false;
47 
48  return true;
49  }
50 
51  return false;
52 }
53 
55 {
56  Interface*& mbImpl = _readNC->mbImpl;
57  std::vector< std::string >& dimNames = _readNC->dimNames;
58  std::vector< int >& dimLens = _readNC->dimLens;
59  std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
60  DebugOutput& dbgOut = _readNC->dbgOut;
61  bool& isParallel = _readNC->isParallel;
62  int& partMethod = _readNC->partMethod;
63  ScdParData& parData = _readNC->parData;
64 
65  // Look for names of center i/j dimensions
66  // First i
67  std::vector< std::string >::iterator vit;
68  unsigned int idx;
69  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lon" ) ) != dimNames.end() )
70  idx = vit - dimNames.begin();
71  else
72  {
73  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' dimension" );
74  }
75  iCDim = idx;
76  gCDims[0] = 0;
77  gCDims[3] = dimLens[idx] - 1;
78 
79  // Check i periodicity and set globallyPeriodic[0]
80  std::vector< double > til_vals( 2 );
81  ErrorCode rval = read_coordinate( "lon", gCDims[3] - 1, gCDims[3], til_vals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
82  if( std::fabs( 2 * til_vals[1] - til_vals[0] - 360 ) < 0.001 ) globallyPeriodic[0] = 1;
83 
84  // Now we can set gDims for i
85  gDims[0] = 0;
86  gDims[3] =
87  gCDims[3] + ( globallyPeriodic[0] ? 0 : 1 ); // Only if not periodic is vertex param max > elem param max
88 
89  // Then j
90  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lat" ) ) != dimNames.end() )
91  idx = vit - dimNames.begin();
92  else
93  {
94  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' dimension" );
95  }
96  jCDim = idx;
97  gCDims[1] = 0;
98  gCDims[4] = dimLens[idx] - 1;
99 
100  // For Eul models, will always be non-periodic in j
101  gDims[1] = 0;
102  gDims[4] = gCDims[4] + 1;
103 
104  // Try a truly 2D mesh
105  gDims[2] = -1;
106  gDims[5] = -1;
107 
108  // Look for time dimension
109  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
110  idx = vit - dimNames.begin();
111  else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
112  idx = vit - dimNames.begin();
113  else
114  {
115  MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
116  }
117  tDim = idx;
118  nTimeSteps = dimLens[idx];
119 
120  // Look for level dimension
121  if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
122  idx = vit - dimNames.begin();
123  else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
124  idx = vit - dimNames.begin();
125  else
126  {
127  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
128  }
129  levDim = idx;
130  nLevels = dimLens[idx];
131 
132  // Parse options to get subset
133  int rank = 0, procs = 1;
134 #ifdef MOAB_HAVE_MPI
135  if( isParallel )
136  {
137  ParallelComm*& myPcomm = _readNC->myPcomm;
138  rank = myPcomm->proc_config().proc_rank();
139  procs = myPcomm->proc_config().proc_size();
140  }
141 #endif
142  if( procs > 1 )
143  {
144  for( int i = 0; i < 6; i++ )
145  parData.gDims[i] = gDims[i];
146  for( int i = 0; i < 3; i++ )
147  parData.gPeriodic[i] = globallyPeriodic[i];
148  parData.partMethod = partMethod;
149  int pdims[3];
150 
151  rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval );
152  for( int i = 0; i < 3; i++ )
153  parData.pDims[i] = pdims[i];
154 
155  dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0] + 1, lDims[4] - lDims[1] + 1,
156  gDims[3] - gDims[0] + 1, gDims[4] - gDims[1] + 1 );
157  if( 0 == rank )
158  dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
159  8 * ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
160  }
161  else
162  {
163  for( int i = 0; i < 6; i++ )
164  lDims[i] = gDims[i];
166  }
167 
168  _opts.get_int_option( "IMIN", lDims[0] );
169  _opts.get_int_option( "IMAX", lDims[3] );
170  _opts.get_int_option( "JMIN", lDims[1] );
171  _opts.get_int_option( "JMAX", lDims[4] );
172 
173  // Now get actual coordinate values for vertices and cell centers
174  lCDims[0] = lDims[0];
175  if( locallyPeriodic[0] )
176  // If locally periodic, doesn't matter what global periodicity is, # vertex coords = # elem
177  // coords
178  lCDims[3] = lDims[3];
179  else
180  lCDims[3] = lDims[3] - 1;
181 
182  // For Eul models, will always be non-periodic in j
183  lCDims[1] = lDims[1];
184  lCDims[4] = lDims[4] - 1;
185 
186  // Resize vectors to store values later
187  if( -1 != lDims[0] ) ilVals.resize( lDims[3] - lDims[0] + 1 );
188  if( -1 != lCDims[0] ) ilCVals.resize( lCDims[3] - lCDims[0] + 1 );
189  if( -1 != lDims[1] ) jlVals.resize( lDims[4] - lDims[1] + 1 );
190  if( -1 != lCDims[1] ) jlCVals.resize( lCDims[4] - lCDims[1] + 1 );
191  if( nTimeSteps > 0 ) tVals.resize( nTimeSteps );
192 
193  // Now read coord values
194  std::map< std::string, ReadNC::VarData >::iterator vmit;
195  if( -1 != lCDims[0] )
196  {
197  if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
198  {
199  rval = read_coordinate( "lon", lCDims[0], lCDims[3], ilCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
200  }
201  else
202  {
203  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
204  }
205  }
206 
207  if( -1 != lCDims[1] )
208  {
209  if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
210  {
211  rval = read_coordinate( "lat", lCDims[1], lCDims[4], jlCVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" );
212  }
213  else
214  {
215  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
216  }
217  }
218 
219  if( -1 != lDims[0] )
220  {
221  if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
222  {
223  double dif = ( ilCVals[1] - ilCVals[0] ) / 2;
224  std::size_t i;
225  for( i = 0; i != ilCVals.size(); i++ )
226  ilVals[i] = ilCVals[i] - dif;
227  // The last one is needed only if not periodic
228  if( !locallyPeriodic[0] ) ilVals[i] = ilCVals[i - 1] + dif;
229  }
230  else
231  {
232  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
233  }
234  }
235 
236  if( -1 != lDims[1] )
237  {
238  if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
239  {
240  if( !isParallel || ( ( gDims[4] - gDims[1] ) == ( lDims[4] - lDims[1] ) ) )
241  {
242  std::string gwName( "gw" );
243  std::vector< double > gwVals( lDims[4] - lDims[1] - 1 );
244  rval = read_coordinate( gwName.c_str(), lDims[1], lDims[4] - 2, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
245  // Copy the correct piece
246  jlVals[0] = -( M_PI / 2 ) * 180 / M_PI;
247  std::size_t i = 0;
248  double gwSum = -1;
249  for( i = 1; i != gwVals.size() + 1; i++ )
250  {
251  gwSum += gwVals[i - 1];
252  jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
253  }
254  jlVals[i] = 90.0; // Using value of i after loop exits.
255  }
256  else
257  {
258  std::string gwName( "gw" );
259  double gwSum = 0;
260 
261  // If this is the first row
262  if( lDims[1] == gDims[1] )
263  {
264  std::vector< double > gwVals( lDims[4] );
265  rval = read_coordinate( gwName.c_str(), 0, lDims[4] - 1, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
266  // Copy the correct piece
267  jlVals[0] = -( M_PI / 2 ) * 180 / M_PI;
268  gwSum = -1;
269  for( std::size_t i = 1; i != jlVals.size(); i++ )
270  {
271  gwSum += gwVals[i - 1];
272  jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
273  }
274  }
275  // Or if it's the last row
276  else if( lDims[4] == gDims[4] )
277  {
278  std::vector< double > gwVals( lDims[4] - 1 );
279  rval = read_coordinate( gwName.c_str(), 0, lDims[4] - 2, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
280  // Copy the correct piece
281  gwSum = -1;
282  for( int j = 0; j != lDims[1] - 1; j++ )
283  gwSum += gwVals[j];
284  std::size_t i = 0;
285  for( ; i != jlVals.size() - 1; i++ )
286  {
287  gwSum += gwVals[lDims[1] - 1 + i];
288  jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
289  }
290  jlVals[i] = 90.0; // Using value of i after loop exits.
291  }
292  // It's in the middle
293  else
294  {
295  int start = lDims[1] - 1;
296  int end = lDims[4] - 1;
297  std::vector< double > gwVals( end );
298  rval = read_coordinate( gwName.c_str(), 0, end - 1, gwVals );MB_CHK_SET_ERR( rval, "Trouble reading 'gw' variable" );
299  gwSum = -1;
300  for( int j = 0; j != start - 1; j++ )
301  gwSum += gwVals[j];
302  std::size_t i = 0;
303  for( ; i != jlVals.size(); i++ )
304  {
305  gwSum += gwVals[start - 1 + i];
306  jlVals[i] = std::asin( gwSum ) * 180 / M_PI;
307  }
308  }
309  }
310  }
311  else
312  {
313  MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
314  }
315  }
316 
317  // Store time coordinate values in tVals
318  if( nTimeSteps > 0 )
319  {
320  if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
321  {
322  rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
323  }
324  else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
325  {
326  rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
327  }
328  else
329  {
330  // If expected time variable is not available, set dummy time coordinate values to tVals
331  for( int t = 0; t < nTimeSteps; t++ )
332  tVals.push_back( (double)t );
333  }
334  }
335 
336  dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
337  dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
338  ( lDims[3] - lDims[0] + 1 ) * ( lDims[4] - lDims[1] + 1 ) );
339 
340  // For each variable, determine the entity location type and number of levels
341  std::map< std::string, ReadNC::VarData >::iterator mit;
342  for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
343  {
344  ReadNC::VarData& vd = ( *mit ).second;
345 
346  // Default entLoc is ENTLOCSET
347  if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
348  {
349  if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
350  ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
352  }
353 
354  // Default numLev is 0
355  if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
356  }
357 
358  // For Eul models, slon and slat are "virtual" dimensions (not defined in the file header)
359  std::vector< std::string > ijdimNames( 4 );
360  ijdimNames[0] = "__slon";
361  ijdimNames[1] = "__slat";
362  ijdimNames[2] = "__lon";
363  ijdimNames[3] = "__lat";
364 
365  std::string tag_name;
366  Tag tagh;
367 
368  // __<dim_name>_LOC_MINMAX (for virtual slon, virtual slat, lon and lat)
369  for( unsigned int i = 0; i != ijdimNames.size(); i++ )
370  {
371  std::vector< int > val( 2, 0 );
372  if( ijdimNames[i] == "__slon" )
373  {
374  val[0] = lDims[0];
375  val[1] = lDims[3];
376  }
377  else if( ijdimNames[i] == "__slat" )
378  {
379  val[0] = lDims[1];
380  val[1] = lDims[4];
381  }
382  else if( ijdimNames[i] == "__lon" )
383  {
384  val[0] = lCDims[0];
385  val[1] = lCDims[3];
386  }
387  else if( ijdimNames[i] == "__lat" )
388  {
389  val[0] = lCDims[1];
390  val[1] = lCDims[4];
391  }
392  std::stringstream ss_tag_name;
393  ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
394  tag_name = ss_tag_name.str();
395  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 );
396  rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
397  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
398  }
399 
400  // __<dim_name>_LOC_VALS (for virtual slon, virtual slat, lon and lat)
401  // Assume all have the same data type as lon (expected type is float or double)
402  switch( varInfo["lon"].varDataType )
403  {
404  case NC_FLOAT:
405  case NC_DOUBLE:
406  break;
407  default:
408  MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
409  }
410 
411  for( unsigned int i = 0; i != ijdimNames.size(); i++ )
412  {
413  void* val = NULL;
414  int val_len = 0;
415  if( ijdimNames[i] == "__slon" )
416  {
417  val = &ilVals[0];
418  val_len = ilVals.size();
419  }
420  else if( ijdimNames[i] == "__slat" )
421  {
422  val = &jlVals[0];
423  val_len = jlVals.size();
424  }
425  else if( ijdimNames[i] == "__lon" )
426  {
427  val = &ilCVals[0];
428  val_len = ilCVals.size();
429  }
430  else if( ijdimNames[i] == "__lat" )
431  {
432  val = &jlCVals[0];
433  val_len = jlCVals.size();
434  }
435 
436  std::stringstream ss_tag_name;
437  ss_tag_name << ijdimNames[i] << "_LOC_VALS";
438  tag_name = ss_tag_name.str();
439  rval = mbImpl->tag_get_handle( tag_name.c_str(), 0, MB_TYPE_DOUBLE, tagh,
440  MB_TAG_CREAT | MB_TAG_SPARSE | MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
441  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 );
442  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
443  }
444 
445  // __<dim_name>_GLOBAL_MINMAX (for virtual slon, virtual slat, lon and lat)
446  for( unsigned int i = 0; i != ijdimNames.size(); i++ )
447  {
448  std::vector< int > val( 2, 0 );
449  if( ijdimNames[i] == "__slon" )
450  {
451  val[0] = gDims[0];
452  val[1] = gDims[3];
453  }
454  else if( ijdimNames[i] == "__slat" )
455  {
456  val[0] = gDims[1];
457  val[1] = gDims[4];
458  }
459  else if( ijdimNames[i] == "__lon" )
460  {
461  val[0] = gCDims[0];
462  val[1] = gCDims[3];
463  }
464  else if( ijdimNames[i] == "__lat" )
465  {
466  val[0] = gCDims[1];
467  val[1] = gCDims[4];
468  }
469  std::stringstream ss_tag_name;
470  ss_tag_name << ijdimNames[i] << "_GLOBAL_MINMAX";
471  tag_name = ss_tag_name.str();
472  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 );
473  rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
474  if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
475  }
476 
477  // Hack: create dummy variables, if needed, for dimensions with no corresponding coordinate
478  // variables
479  rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
480 
481  return MB_SUCCESS;
482 }
483 
484 } // namespace moab