![]() |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
00001 /*
00002 * addncdata.cpp
00003 * this tool will take an existing h5m file and add data from an nc type file
00004 * will support mainly showing the data associated with unstructured meshes (Homme, MPAS) with Visit
00005 *
00006 * example of usage:
00007 * ./mbaddnc -i wholeFineATM.h5m -n surfdata_ne11np4_simyr1850_c160614.nc -o
00008 * whole_LONGXY_surfdata.h5m -v LONGXY
00009 *
00010 * Basically, will output a new h5m file (whole_LONGXY_surfdata.h5m), which has an extra tag,
00011 * corresponding to the variable LONGXY from the file surfdata_ne11np4_simyr1850_c160614.nc;
00012 * matching is based on the global ids between what we think is the order on the original file
00013 * (wholeFineATM.h5m) and the order of surfdata_ne11np4_simyr1850_c160614.nc
00014 *
00015 * file wholeFineATM.h5m is obtained from a coupled run in e3sm, with the ne 11, np 4,
00016 *
00017 * add an option to output the nc data to the original coarse ATM file, the one that also has the
00018 * GLOBAL_DOFS tag with the global DOFs of the GLL points
00019 */
00020
00021 #include "moab/ProgOptions.hpp"
00022 #include "moab/Core.hpp"
00023
00024 #include "netcdf.h"
00025 #include
00026 #include
00027
00028 using namespace moab;
00029
00030 // copy from ReadNCDF.cpp some useful macros for reading from a netcdf file (exodus?)
00031 // ncFile is an integer initialized when opening the nc file in read mode
00032
00033 int ncFile;
00034
00035 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00036
00037 #define GET_DIM( ncdim, name, val ) \
00038 { \
00039 int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) ); \
00040 if( NC_NOERR == gdfail ) \
00041 { \
00042 size_t tmp_val; \
00043 gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val ); \
00044 if( NC_NOERR != gdfail ) \
00045 { \
00046 MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
00047 } \
00048 else \
00049 ( val ) = tmp_val; \
00050 } \
00051 else \
00052 ( val ) = 0; \
00053 }
00054
00055 #define GET_DIMB( ncdim, name, varname, id, val ) \
00056 INS_ID( name, varname, id ); \
00057 GET_DIM( ncdim, name, val );
00058
00059 #define GET_VAR( name, id, dims ) \
00060 { \
00061 ( id ) = -1; \
00062 int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
00063 if( NC_NOERR == gvfail ) \
00064 { \
00065 int ndims; \
00066 gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
00067 if( NC_NOERR == gvfail ) \
00068 { \
00069 ( dims ).resize( ndims ); \
00070 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
00071 if( NC_NOERR != gvfail ) \
00072 { \
00073 MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable dimension IDs" ); \
00074 } \
00075 } \
00076 } \
00077 }
00078
00079 #define GET_1D_INT_VAR( name, id, vals ) \
00080 { \
00081 GET_VAR( name, id, vals ); \
00082 if( -1 != ( id ) ) \
00083 { \
00084 size_t ntmp; \
00085 int ivfail = nc_inq_dimlen( ncFile, ( vals )[0], &ntmp ); \
00086 if( NC_NOERR != ivfail ) \
00087 { \
00088 MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
00089 } \
00090 ( vals ).resize( ntmp ); \
00091 size_t ntmp1 = 0; \
00092 ivfail = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
00093 if( NC_NOERR != ivfail ) \
00094 { \
00095 MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
00096 } \
00097 } \
00098 }
00099
00100 #define GET_1D_DBL_VAR( name, id, vals ) \
00101 { \
00102 std::vector< int > dum_dims; \
00103 GET_VAR( name, id, dum_dims ); \
00104 if( -1 != ( id ) ) \
00105 { \
00106 size_t ntmp; \
00107 int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
00108 if( NC_NOERR != dvfail ) \
00109 { \
00110 MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
00111 } \
00112 ( vals ).resize( ntmp ); \
00113 size_t ntmp1 = 0; \
00114 dvfail = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
00115 if( NC_NOERR != dvfail ) \
00116 { \
00117 MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
00118 } \
00119 } \
00120 }
00121
00122 #define GET_1D_FLT_VAR( name, id, vals ) \
00123 { \
00124 std::vector< int > dum_dims; \
00125 GET_VAR( name, id, dum_dims ); \
00126 if( -1 != ( id ) ) \
00127 { \
00128 size_t ntmp; \
00129 int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp ); \
00130 if( NC_NOERR != dvfail ) \
00131 { \
00132 MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); \
00133 } \
00134 ( vals ).resize( ntmp ); \
00135 size_t ntmp1 = 0; \
00136 dvfail = nc_get_vara_float( ncFile, id, &ntmp1, &ntmp, &( vals )[0] ); \
00137 if( NC_NOERR != dvfail ) \
00138 { \
00139 MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << ( name ) ); \
00140 } \
00141 } \
00142 }
00143
00144 int main( int argc, char* argv[] )
00145 {
00146 ErrorCode rval;
00147 ProgOptions opts;
00148
00149 std::string inputfile, outfile( "out.h5m" ), netcdfFile, variable_name, sefile_name;
00150
00151 opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
00152 opts.addOpt< std::string >( "netcdfFile,n", "netcdf file aligned with the mesh input file", &netcdfFile );
00153 opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );
00154
00155 opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name );
00156
00157 opts.addOpt< std::string >( "sefile,s", "spectral elements file (coarse SE mesh)", &sefile_name );
00158 opts.parseCommandLine( argc, argv );
00159
00160 Core* mb = new Core();
00161
00162 rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );
00163
00164 std::cout << " opened " << inputfile << " with initial h5m data.\n";
00165 // open the netcdf file, and see if it has that variable we are looking for
00166
00167 Range nodes;
00168 rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );
00169
00170 Range edges;
00171 rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );
00172
00173 Range cells;
00174 rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );
00175
00176 std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n";
00177
00178 // construct maps between global id and handles
00179 std::map< int, EntityHandle > vGidHandle;
00180 std::map< int, EntityHandle > eGidHandle;
00181 std::map< int, EntityHandle > cGidHandle;
00182 std::vector< int > gids;
00183 Tag gid;
00184 rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
00185 gids.resize( nodes.size() );
00186 rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" );
00187 int i = 0;
00188 for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
00189 {
00190 vGidHandle[gids[i++]] = *vit;
00191 }
00192
00193 gids.resize( edges.size() );
00194 rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" );
00195 i = 0;
00196 for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
00197 {
00198 eGidHandle[gids[i++]] = *vit;
00199 }
00200
00201 gids.resize( cells.size() );
00202 rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" );
00203 i = 0;
00204 for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
00205 {
00206 cGidHandle[gids[i++]] = *vit;
00207 }
00208
00209 // Open netcdf/exodus file
00210 int fail = nc_open( netcdfFile.c_str(), 0, &ncFile );
00211 if( NC_NOWRITE != fail )
00212 {
00213 MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf II file " << netcdfFile );
00214 }
00215
00216 std::cout << " opened " << netcdfFile << " with new data \n";
00217 std::vector< int > dims;
00218 int nc_var;
00219
00220 size_t recs;
00221 char recname[NC_MAX_NAME + 1];
00222
00223 std::cout << " looking for variable " << variable_name << "\n";
00224 GET_VAR( variable_name.c_str(), nc_var, dims );
00225 std::cout << " it has " << dims.size() << " dimensions\n";
00226
00227 int dimIndex = -1; // index of the dimension of interest
00228 bool vertex_data = false;
00229 bool cell_data = false;
00230 for( size_t j = 0; j < dims.size(); j++ )
00231 {
00232 fail = nc_inq_dim( ncFile, dims[j], recname, &recs );
00233 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
00234 std::string name_dim( recname );
00235 std::cout << " dimension index " << j << " in file: " << dims[j] << " name: " << name_dim << " recs:" << recs
00236 << "\n";
00237 if( recs == nodes.size() )
00238 {
00239 dimIndex = j;
00240 vertex_data = true;
00241 }
00242 if( recs == cells.size() )
00243 {
00244 dimIndex = j;
00245 cell_data = true;
00246 }
00247 }
00248
00249 int otherDim = 1 - dimIndex; // used only if 2 dimensions ; could be 0 or 1;
00250 size_t size_tag = 1; // good for one dimension
00251 std::vector< int > evals; // size of size_tag
00252 std::vector< double > dvals; // size of size_tag
00253 nc_type dataType;
00254
00255 // read the variable, and set it to the tag
00256 fail = nc_inq_vartype( ncFile, nc_var, &dataType );
00257 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable type" );
00258 DataType mbtype = MB_TYPE_DOUBLE;
00259 bool float_var = false;
00260 if( NC_INT == dataType )
00261 mbtype = MB_TYPE_INTEGER;
00262 else if( NC_DOUBLE != dataType && NC_FLOAT != dataType )
00263 MB_CHK_SET_ERR( MB_FAILURE, "unknown type" );
00264
00265 if( NC_FLOAT == dataType ) float_var = true;
00266
00267 int time_id = -1;
00268 fail = nc_inq_varid( ncFile, "time", &time_id );
00269 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get time variable id" );
00270 std::vector< float > times;
00271 if( NC_NOERR == fail )
00272 {
00273 int ii;
00274 GET_1D_FLT_VAR( "time", ii, times );
00275 }
00276 Tag newTag;
00277
00278 if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data ) )
00279 {
00280
00281 if( dims.size() == 2 )
00282 {
00283 fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag );
00284 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
00285 }
00286
00287 int def_val = 0;
00288 rval = mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag, MB_TAG_CREAT | MB_TAG_DENSE,
00289 &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" );
00290
00291 if( NC_INT == dataType )
00292 {
00293 std::vector< int > vals;
00294 if( 1 == dims.size() )
00295 {
00296 GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals );
00297 if( vertex_data )
00298 {
00299 for( size_t k = 0; k < vals.size(); k++ )
00300 {
00301 EntityHandle vh = vGidHandle[k + 1];
00302 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
00303 }
00304 }
00305 else // cell_data
00306 {
00307 for( size_t k = 0; k < vals.size(); k++ )
00308 {
00309 EntityHandle ch = cGidHandle[k + 1]; // cell handle
00310 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
00311 }
00312 }
00313 }
00314 else // dims.size() == 2
00315 {
00316 // Single var for all coords
00317 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
00318 count[dimIndex] = recs;
00319 count[otherDim] = size_tag;
00320 vals.resize( recs * size_tag );
00321 fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] );
00322 evals.resize( size_tag );
00323
00324 if( vertex_data )
00325 {
00326 for( size_t k = 0; k < recs; k++ )
00327 {
00328 EntityHandle vh = vGidHandle[k + 1];
00329 size_t start_in_vals = k * size_tag, stride = 1;
00330 for( size_t j = 0; j < size_tag; j++ )
00331 evals[j] = vals[start_in_vals + j * stride];
00332 rval = mb->tag_set_data( newTag, &vh, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
00333 }
00334 }
00335 else // cell_data
00336 {
00337 for( size_t k = 0; k < recs; k++ )
00338 {
00339 EntityHandle ch = cGidHandle[k + 1]; // cell handle
00340 size_t start_in_vals = k * size_tag, stride = 1;
00341 for( size_t j = 0; j < size_tag; j++ )
00342 evals[j] = vals[start_in_vals + j * stride];
00343 rval = mb->tag_set_data( newTag, &ch, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
00344 }
00345 }
00346 }
00347 }
00348 else
00349 {
00350 std::vector< double > vals;
00351 if( 1 == dims.size() )
00352 {
00353 GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals );
00354 if( vertex_data )
00355 {
00356 for( size_t k = 0; k < vals.size(); k++ )
00357 {
00358 EntityHandle vh = vGidHandle[k + 1];
00359 rval = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
00360 }
00361 }
00362 else // cell_data
00363 {
00364 for( size_t k = 0; k < vals.size(); k++ )
00365 {
00366 EntityHandle ch = cGidHandle[k + 1]; // cell handle
00367 rval = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
00368 }
00369 }
00370 }
00371 else // dims.size() == 2
00372 {
00373 // Single var for all coords
00374 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
00375 count[dimIndex] = recs;
00376 count[otherDim] = size_tag;
00377 vals.resize( recs * size_tag );
00378 fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] );
00379 dvals.resize( size_tag );
00380
00381 if( vertex_data )
00382 {
00383 for( size_t k = 0; k < recs; k++ )
00384 {
00385 EntityHandle vh = vGidHandle[k + 1];
00386 size_t start_in_vals = k * size_tag, stride = 1;
00387 for( size_t j = 0; j < size_tag; j++ )
00388 dvals[j] = vals[start_in_vals + j * stride];
00389 rval = mb->tag_set_data( newTag, &vh, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
00390 }
00391 }
00392 else // cell_data
00393 {
00394 for( size_t k = 0; k < recs; k++ )
00395 {
00396 EntityHandle ch = cGidHandle[k + 1]; // cell handle
00397 size_t start_in_vals = k * size_tag, stride = 1;
00398 for( size_t j = 0; j < size_tag; j++ )
00399 dvals[j] = vals[start_in_vals + j * stride];
00400 rval = mb->tag_set_data( newTag, &ch, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
00401 }
00402 }
00403 }
00404 }
00405 }
00406 else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 &&
00407 mbtype == MB_TYPE_DOUBLE ) // the last one is the vertex
00408 {
00409 // the case when the last dim is ncol (for homme type mesh))
00410 size_t dim0, dim1; // dim 2 is ncol..
00411 char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1];
00412 fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 );
00413 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
00414 fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 );
00415 if( NC_NOERR != fail ) MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension" );
00416 std::string name0( recname0 );
00417 std::string name1( recname1 );
00418 std::string timestr( "time" );
00419 size_t start[3] = { 0, 0, 0 };
00420 size_t count[3] = { 1, 0, 0 };
00421 count[1] = dim1;
00422 count[2] = nodes.size();
00423 // create a few tags, with name inserted
00424 if( name0.compare( "time" ) == 0 )
00425 {
00426
00427 std::vector< double > dvalues( dim1 * count[2] );
00428 std::vector< float > fvalues( dim1 * count[2] );
00429 std::vector< double > transp( dim1 * count[2] );
00430 // get the variable time
00431 for( size_t k = 0; k < dim0; k++ )
00432 {
00433 // create a tag for each time, and
00434 std::stringstream tag_name;
00435 tag_name << variable_name << "_t" << times[k];
00436 std::vector< double > defvals( dim1, 0. );
00437 rval = mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag,
00438 MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );
00439 start[0] = k;
00440
00441 if( float_var )
00442 {
00443 fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] );
00444 // now arrange them in a tag, transpose data
00445 for( size_t ii = 0; ii < dim1; ii++ )
00446 for( size_t j = 0; j < count[2]; j++ )
00447 {
00448 transp[j * dim1 + ii] = fvalues[ii * count[2] + j];
00449 }
00450 }
00451 else // double
00452 {
00453 fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] );
00454 // now arrange them in a tag, transpose data
00455 for( size_t ii = 0; ii < dim1; ii++ )
00456 for( size_t j = 0; j < count[2]; j++ )
00457 {
00458 transp[j * dim1 + ii] = dvalues[ii * count[2] + j];
00459 }
00460 }
00461 for( size_t ii = 0; ii < nodes.size(); ii++ )
00462 {
00463 EntityHandle vh = vGidHandle[ii + 1];
00464 rval = mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] );MB_CHK_SET_ERR( rval, "can't set tag on nodes" );
00465 }
00466 }
00467 }
00468 }
00469 rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
00470 std::cout << " wrote file " << outfile << "\n";
00471
00472 // now, if s option, load the coarse mesh and put data on each element, according to a matrix
00473 if( !sefile_name.empty() )
00474 {
00475 // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated
00476 Core* mb2 = new Core();
00477 rval = mb2->load_file( sefile_name.c_str() );MB_CHK_SET_ERR( rval, "can't load spectral element file" );
00478 std::cout << " loaded spectral file " << sefile_name << "\n";
00479 // look for GLOBAL_DOFS tag
00480 Tag gdofeTag;
00481 rval = mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag );MB_CHK_SET_ERR( rval, "file does not have GLOBAL_DOFS file" );
00482 int sizeTag;
00483 rval = mb2->tag_get_length( gdofeTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag" );
00484 int np = (int)sqrt( 1.0 * sizeTag );
00485 std::cout << " size of tag: " << sizeTag << " np = " << np << "\n";
00486 std::vector< int > gdofs;
00487 Range cells2;
00488 rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells on spectral mesh" );
00489 gdofs.resize( cells2.size() * sizeTag );
00490 rval = mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] );MB_CHK_SET_ERR( rval, "can't get global dofs tag" );
00491 // create a new tag for element data arranged as DFIELD
00492
00493 std::vector< double > dfield;
00494 dfield.resize( sizeTag, 0.0 );
00495 Tag newTag2;
00496 rval = mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2, MB_TAG_CREAT | MB_TAG_DENSE,
00497 &dfield[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );
00498
00499 int i1 = 0; // index in the gdofs array, per element
00500
00501 // get the tag values from the other moab core, for newTag
00502 int dataTagLen;
00503 rval = mb->tag_get_length( newTag, dataTagLen );MB_CHK_SET_ERR( rval, "can't get size of newTag" );
00504 //
00505 std::vector< double > oldData;
00506 oldData.resize( dataTagLen * nodes.size() ); //
00507
00508 // get the "old" values
00509 rval = mb->tag_get_data( newTag, nodes, &oldData[0] );MB_CHK_SET_ERR( rval, "can't get old values" );
00510 for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ )
00511 {
00512 EntityHandle cel = *it;
00513 // gdofs per element are gdofs[i:i+np*np];
00514 for( int k = 0; k < sizeTag; k++ )
00515 {
00516 int gdof = gdofs[i1 + k];
00517 EntityHandle node = vGidHandle[gdof];
00518 int index = nodes.index( node );
00519 // get last layer of data
00520 double val = oldData[index * dataTagLen + dataTagLen - 1]; //
00521 dfield[k] = val;
00522 }
00523 i1 = i1 + sizeTag;
00524 rval = mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] );MB_CHK_SET_ERR( rval, "can't set new tag" );
00525 }
00526
00527 // write the appended file with the new field:
00528 rval = mb2->write_file( "atm2.h5m" );MB_CHK_SET_ERR( rval, "can't write new spectral file" );
00529 std::cout << " wrote file atm2.h5m \n";
00530 }
00531 return 0;
00532 }
00533