tNC3File.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tNC3File.cc (27337B)
---
1 // Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2019 PISM Authors
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include "NC3File.hh"
20
21 // The following is a stupid kludge necessary to make NetCDF 4.x work in
22 // serial mode in an MPI program:
23 #ifndef MPI_INCLUDED
24 #define MPI_INCLUDED 1
25 #endif
26 #include <netcdf.h>
27 #include <cstring> // memset
28 #include <cstdio> // stderr, fprintf
29
30 #include "pism/util/pism_utilities.hh" // join
31 #include "pism/util/error_handling.hh"
32
33 #include "pism_type_conversion.hh" // This has to be included *after* netcdf.h.
34
35 namespace pism {
36 namespace io {
37
38 //! \brief Prints an error message; for debugging.
39 static void check(const ErrorLocation &where, int return_code) {
40 if (return_code != NC_NOERR) {
41 throw RuntimeError(where, nc_strerror(return_code));
42 }
43 }
44
45 //! call MPI_Abort() if a NetCDF call failed
46 static void check_and_abort(MPI_Comm com, const ErrorLocation &where, int return_code) {
47 if (return_code != NC_NOERR) {
48 fprintf(stderr, "%s:%d: %s\n", where.filename, where.line_number, nc_strerror(return_code));
49 MPI_Abort(com, -1);
50 }
51 }
52
53 NC3File::NC3File(MPI_Comm c)
54 : NCFile(c), m_rank(0) {
55 MPI_Comm_rank(m_com, &m_rank);
56 }
57
58 NC3File::~NC3File() {
59 if (m_file_id >= 0) {
60 if (m_rank == 0) {
61 nc_close(m_file_id);
62 fprintf(stderr, "NC3File::~NC3File: NetCDF file %s is still open\n",
63 m_filename.c_str());
64 }
65 m_file_id = -1;
66 }
67 }
68
69 // open/create/close
70 void NC3File::open_impl(const std::string &fname, IO_Mode mode) {
71 int stat = NC_NOERR;
72
73 int open_mode = mode == PISM_READONLY ? NC_NOWRITE : NC_WRITE;
74
75 if (m_rank == 0) {
76 stat = nc_open(fname.c_str(), open_mode, &m_file_id);
77 }
78
79 MPI_Barrier(m_com);
80 MPI_Bcast(&m_file_id, 1, MPI_INT, 0, m_com);
81 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
82
83 check(PISM_ERROR_LOCATION, stat);
84 }
85
86 //! \brief Create a NetCDF file.
87 void NC3File::create_impl(const std::string &fname) {
88 int stat = NC_NOERR;
89
90 if (m_rank == 0) {
91 stat = nc_create(fname.c_str(), NC_CLOBBER | NC_64BIT_OFFSET, &m_file_id);
92 }
93
94 MPI_Barrier(m_com);
95 MPI_Bcast(&m_file_id, 1, MPI_INT, 0, m_com);
96 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
97
98 check(PISM_ERROR_LOCATION, stat);
99 }
100
101 //! \brief Close a NetCDF file.
102 void NC3File::close_impl() {
103 int stat = NC_NOERR;
104
105 if (m_rank == 0) {
106 stat = nc_close(m_file_id);
107 }
108
109 m_file_id = -1;
110
111 MPI_Barrier(m_com);
112 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
113
114 check(PISM_ERROR_LOCATION, stat);
115 }
116
117
118 void NC3File::sync_impl() const {
119 int stat = NC_NOERR;
120
121 if (m_rank == 0) {
122 stat = nc_sync(m_file_id);
123 }
124
125 MPI_Barrier(m_com);
126 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
127
128 check(PISM_ERROR_LOCATION, stat);
129 }
130
131
132 //! \brief Exit define mode.
133 void NC3File::enddef_impl() const {
134 int stat = NC_NOERR;
135
136 int header_size = 200 * 1024;
137
138 if (m_rank == 0) {
139 stat = nc__enddef(m_file_id, header_size, 4, 0, 4);
140 }
141
142 MPI_Barrier(m_com);
143 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
144
145 check(PISM_ERROR_LOCATION, stat);
146 }
147
148 //! \brief Enter define mode.
149 void NC3File::redef_impl() const {
150 int stat = NC_NOERR;
151
152 if (m_rank == 0) {
153 stat = nc_redef(m_file_id);
154 }
155
156 MPI_Barrier(m_com);
157 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
158
159 check(PISM_ERROR_LOCATION, stat);
160 }
161
162
163 //! \brief Define a dimension.
164 void NC3File::def_dim_impl(const std::string &name, size_t length) const {
165 int stat = NC_NOERR;
166
167 if (m_rank == 0) {
168 int dimid;
169 stat = nc_def_dim(m_file_id, name.c_str(), length, &dimid);
170 }
171
172 MPI_Barrier(m_com);
173 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
174
175 check(PISM_ERROR_LOCATION, stat);
176 }
177
178 void NC3File::inq_dimid_impl(const std::string &dimension_name, bool &exists) const {
179 int stat, flag = -1;
180
181 if (m_rank == 0) {
182 stat = nc_inq_dimid(m_file_id, dimension_name.c_str(), &flag);
183
184 flag = (stat == NC_NOERR) ? 1 : 0;
185 }
186 MPI_Barrier(m_com);
187 MPI_Bcast(&flag, 1, MPI_INT, 0, m_com);
188
189 exists = (flag == 1);
190 }
191
192
193 //! \brief Get a dimension length.
194 void NC3File::inq_dimlen_impl(const std::string &dimension_name, unsigned int &result) const {
195 int stat = NC_NOERR;
196
197 if (m_rank == 0) {
198 int dimid;
199 size_t length;
200
201 stat = nc_inq_dimid(m_file_id, dimension_name.c_str(), &dimid);
202
203 if (stat == NC_NOERR) {
204 stat = nc_inq_dimlen(m_file_id, dimid, &length);
205 result = static_cast<unsigned int>(length);
206 }
207 }
208
209 MPI_Barrier(m_com);
210 MPI_Bcast(&result, 1, MPI_UNSIGNED, 0, m_com);
211 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
212
213 check(PISM_ERROR_LOCATION, stat);
214 }
215
216 //! \brief Get an unlimited dimension.
217 void NC3File::inq_unlimdim_impl(std::string &result) const {
218 int stat = NC_NOERR;
219 std::vector<char> dimname(NC_MAX_NAME + 1, 0);
220
221 if (m_rank == 0) {
222 int dimid;
223 stat = nc_inq_unlimdim(m_file_id, &dimid);
224
225 // nc_inq_unlimdim() sets dimid to -1 if there is no unlimited dimension
226 if (stat == NC_NOERR and dimid != -1) {
227 stat = nc_inq_dimname(m_file_id, dimid, dimname.data());
228 } else {
229 // leave dimname empty
230 stat = NC_NOERR;
231 }
232 }
233
234 MPI_Barrier(m_com);
235
236 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
237 MPI_Bcast(dimname.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
238
239 check(PISM_ERROR_LOCATION, stat);
240
241 result = dimname.data();
242 }
243
244 //! \brief Define a variable.
245 void NC3File::def_var_impl(const std::string &name,
246 IO_Type nctype,
247 const std::vector<std::string> &dims) const {
248 int stat = NC_NOERR;
249
250 if (m_rank == 0) {
251 std::vector<int> dimids;
252 int varid;
253
254 for (auto d : dims) {
255 int dimid;
256 stat = nc_inq_dimid(m_file_id, d.c_str(), &dimid);
257 if (stat != NC_NOERR) {
258 break;
259 }
260 dimids.push_back(dimid);
261 }
262
263 if (stat == NC_NOERR) {
264 stat = nc_def_var(m_file_id, name.c_str(), pism_type_to_nc_type(nctype),
265 static_cast<int>(dims.size()), &dimids[0], &varid);
266 }
267 }
268
269 MPI_Barrier(m_com);
270 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
271
272 check(PISM_ERROR_LOCATION, stat);
273 }
274
275 void NC3File::get_varm_double_impl(const std::string &variable_name,
276 const std::vector<unsigned int> &start,
277 const std::vector<unsigned int> &count,
278 const std::vector<unsigned int> &imap, double *op) const {
279 return this->get_var_double(variable_name,
280 start, count, imap, op, true);
281 }
282
283 void NC3File::get_vara_double_impl(const std::string &variable_name,
284 const std::vector<unsigned int> &start,
285 const std::vector<unsigned int> &count,
286 double *op) const {
287 std::vector<unsigned int> dummy;
288 return this->get_var_double(variable_name,
289 start, count, dummy, op, false);
290 }
291
292 //! \brief Get variable data.
293 void NC3File::get_var_double(const std::string &variable_name,
294 const std::vector<unsigned int> &start_input,
295 const std::vector<unsigned int> &count_input,
296 const std::vector<unsigned int> &imap_input, double *ip,
297 bool transposed) const {
298 std::vector<unsigned int> start = start_input;
299 std::vector<unsigned int> count = count_input;
300 std::vector<unsigned int> imap = imap_input;
301 const int start_tag = 1,
302 count_tag = 2,
303 data_tag = 3,
304 imap_tag = 4,
305 chunk_size_tag = 5;
306 int stat = NC_NOERR, com_size, ndims = static_cast<int>(start.size());
307 std::vector<double> processor_0_buffer;
308 MPI_Status mpi_stat;
309 unsigned int local_chunk_size = 1,
310 processor_0_chunk_size = 0;
311
312 if (not transposed) {
313 imap.resize(ndims);
314 }
315
316 // get the size of the communicator
317 MPI_Comm_size(m_com, &com_size);
318
319 // compute the size of a local chunk
320 for (int k = 0; k < ndims; ++k) {
321 local_chunk_size *= count[k];
322 }
323
324 // compute the maximum and send it to processor 0; this is the size of the
325 // buffer processor 0 will need
326 MPI_Reduce(&local_chunk_size, &processor_0_chunk_size, 1, MPI_UNSIGNED, MPI_MAX, 0, m_com);
327
328 // now we need to send start, count and imap data to processor 0 and receive data
329 if (m_rank == 0) {
330 // Note: this could be optimized: if processor_0_chunk_size <=
331 // max(local_chunk_size) we can avoid allocating this buffer. The inner for
332 // loop will have to be re-ordered, though.
333 processor_0_buffer.resize(processor_0_chunk_size);
334
335 // MPI calls below require C datatypes (so that we don't have to worry
336 // about sizes of size_t and ptrdiff_t), so we make local copies of start,
337 // count, and imap to use in the nc_get_varm_double() call.
338 std::vector<size_t> nc_start(ndims), nc_count(ndims);
339 std::vector<ptrdiff_t> nc_imap(ndims), nc_stride(ndims);
340 int varid;
341
342 stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
343 check_and_abort(m_com, PISM_ERROR_LOCATION, stat);
344
345 for (int r = 0; r < com_size; ++r) {
346
347 if (r != 0) {
348 // Note: start, count, imap, and local_chunk_size on processor zero are
349 // used *before* they get overwritten by these calls
350 MPI_Recv(&start[0], ndims, MPI_UNSIGNED, r, start_tag, m_com, &mpi_stat);
351 MPI_Recv(&count[0], ndims, MPI_UNSIGNED, r, count_tag, m_com, &mpi_stat);
352 MPI_Recv(&imap[0], ndims, MPI_UNSIGNED, r, imap_tag, m_com, &mpi_stat);
353 MPI_Recv(&local_chunk_size, 1, MPI_UNSIGNED, r, chunk_size_tag, m_com, &mpi_stat);
354 }
355
356 // This for loop uses start, count and imap passed in as arguments when r
357 // == 0. For r > 0 they are overwritten by MPI_Recv calls above.
358 for (int k = 0; k < ndims; ++k) {
359 nc_start[k] = start[k];
360 nc_count[k] = count[k];
361 nc_imap[k] = imap[k];
362 nc_stride[k] = 1; // fill with ones; this way it works even with
363 // NetCDF versions with a bug affecting the
364 // stride == NULL case.
365 }
366
367 if (transposed) {
368 stat = nc_get_varm_double(m_file_id, varid, &nc_start[0], &nc_count[0], &nc_stride[0], &nc_imap[0],
369 &processor_0_buffer[0]);
370 } else {
371 stat = nc_get_vara_double(m_file_id, varid, &nc_start[0], &nc_count[0],
372 &processor_0_buffer[0]);
373 }
374 check_and_abort(m_com, PISM_ERROR_LOCATION, stat);
375
376 if (r != 0) {
377 MPI_Send(&processor_0_buffer[0], local_chunk_size, MPI_DOUBLE, r, data_tag, m_com);
378 } else {
379 for (unsigned int k = 0; k < local_chunk_size; ++k) {
380 ip[k] = processor_0_buffer[k];
381 }
382 }
383 } // end of the for loop
384
385 } else {
386 MPI_Send(&start[0], ndims, MPI_UNSIGNED, 0, start_tag, m_com);
387 MPI_Send(&count[0], ndims, MPI_UNSIGNED, 0, count_tag, m_com);
388 MPI_Send(&imap[0], ndims, MPI_UNSIGNED, 0, imap_tag, m_com);
389 MPI_Send(&local_chunk_size, 1, MPI_UNSIGNED, 0, chunk_size_tag, m_com);
390
391 MPI_Recv(ip, local_chunk_size, MPI_DOUBLE, 0, data_tag, m_com, &mpi_stat);
392 }
393 }
394
395 void NC3File::put_vara_double_impl(const std::string &variable_name,
396 const std::vector<unsigned int> &start_input,
397 const std::vector<unsigned int> &count_input,
398 const double *op) const {
399 // make copies of start and count so that we can use them in MPI_Recv() calls below
400 std::vector<unsigned int> start = start_input;
401 std::vector<unsigned int> count = count_input;
402 const int start_tag = 1,
403 count_tag = 2,
404 data_tag = 3,
405 chunk_size_tag = 4;
406 int stat = NC_NOERR, com_size = 0, ndims = static_cast<int>(start.size());
407 std::vector<double> processor_0_buffer;
408 MPI_Status mpi_stat;
409 unsigned int local_chunk_size = 1,
410 processor_0_chunk_size = 0;
411
412 // get the size of the communicator
413 MPI_Comm_size(m_com, &com_size);
414
415 // compute the size of a local chunk
416 for (int k = 0; k < ndims; ++k) {
417 local_chunk_size *= count[k];
418 }
419
420 // compute the maximum and send it to processor 0; this is the size of the
421 // buffer processor 0 will need
422 MPI_Reduce(&local_chunk_size, &processor_0_chunk_size, 1, MPI_UNSIGNED, MPI_MAX, 0, m_com);
423
424 // now we need to send start and count data to processor 0 and receive data
425 if (m_rank == 0) {
426 processor_0_buffer.resize(processor_0_chunk_size);
427
428 // MPI calls below require C datatypes (so that we don't have to worry about sizes of
429 // size_t and ptrdiff_t), so we make local copies of start and count to use in the
430 // nc_get_vara_double() call.
431 std::vector<size_t> nc_start(ndims), nc_count(ndims);
432 std::vector<ptrdiff_t> nc_stride(ndims);
433 int varid;
434
435 stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
436 check_and_abort(m_com, PISM_ERROR_LOCATION, stat);
437
438 for (int r = 0; r < com_size; ++r) {
439
440 if (r != 0) {
441 // Note: start, count, and local_chunk_size on processor zero are used *before*
442 // they get overwritten by these calls
443 MPI_Recv(&start[0], ndims, MPI_UNSIGNED, r, start_tag, m_com, &mpi_stat);
444 MPI_Recv(&count[0], ndims, MPI_UNSIGNED, r, count_tag, m_com, &mpi_stat);
445 MPI_Recv(&local_chunk_size, 1, MPI_UNSIGNED, r, chunk_size_tag, m_com, &mpi_stat);
446
447 MPI_Recv(&processor_0_buffer[0], local_chunk_size, MPI_DOUBLE, r, data_tag, m_com, &mpi_stat);
448 } else {
449 for (unsigned int k = 0; k < local_chunk_size; ++k) {
450 processor_0_buffer[k] = op[k];
451 }
452 }
453
454 // This for loop uses start and count passed in as arguments when r == 0. For r > 0
455 // they are overwritten by MPI_Recv calls above.
456 for (int k = 0; k < ndims; ++k) {
457 nc_start[k] = start[k];
458 nc_count[k] = count[k];
459 nc_stride[k] = 1; // fill with ones; this way it works even with
460 // NetCDF versions with a bug affecting the
461 // stride == NULL case.
462 }
463
464 stat = nc_put_vara_double(m_file_id, varid, &nc_start[0], &nc_count[0],
465 &processor_0_buffer[0]);
466 check_and_abort(m_com, PISM_ERROR_LOCATION, stat);
467 } // end of the for loop
468 } else {
469 MPI_Send(&start[0], ndims, MPI_UNSIGNED, 0, start_tag, m_com);
470 MPI_Send(&count[0], ndims, MPI_UNSIGNED, 0, count_tag, m_com);
471 MPI_Send(&local_chunk_size, 1, MPI_UNSIGNED, 0, chunk_size_tag, m_com);
472
473 MPI_Send(const_cast<double*>(op), local_chunk_size, MPI_DOUBLE, 0, data_tag, m_com);
474 }
475 }
476
477 //! \brief Get the number of variables.
478 void NC3File::inq_nvars_impl(int &result) const {
479 int stat = NC_NOERR;
480
481 if (m_rank == 0) {
482 stat = nc_inq_nvars(m_file_id, &result);
483 }
484 MPI_Barrier(m_com);
485
486 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
487 check(PISM_ERROR_LOCATION, stat);
488
489 MPI_Bcast(&result, 1, MPI_INT, 0, m_com);
490 }
491
492 //! \brief Get dimensions a variable depends on.
493 void NC3File::inq_vardimid_impl(const std::string &variable_name,
494 std::vector<std::string> &result) const {
495 int stat, ndims, varid = -1;
496 std::vector<int> dimids;
497
498 if (m_rank == 0) {
499 stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
500
501 if (stat == NC_NOERR) {
502 stat = nc_inq_varndims(m_file_id, varid, &ndims);
503 }
504 }
505
506 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
507 check(PISM_ERROR_LOCATION, stat);
508
509 MPI_Bcast(&ndims, 1, MPI_INT, 0, m_com);
510
511 if (ndims == 0) {
512 result.clear();
513 return;
514 }
515
516 result.resize(ndims);
517 dimids.resize(ndims);
518
519 if (m_rank == 0) {
520 stat = nc_inq_vardimid(m_file_id, varid, &dimids[0]);
521 }
522
523 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
524 check(PISM_ERROR_LOCATION, stat);
525
526 MPI_Barrier(m_com);
527
528 for (int k = 0; k < ndims; ++k) {
529 std::vector<char> name(NC_MAX_NAME + 1, 0);
530
531 if (m_rank == 0) {
532 stat = nc_inq_dimname(m_file_id, dimids[k], name.data());
533 }
534
535 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
536 check(PISM_ERROR_LOCATION, stat);
537
538 MPI_Barrier(m_com);
539 MPI_Bcast(name.data(), name.size(), MPI_CHAR, 0, m_com);
540
541 result[k] = name.data();
542 }
543 }
544
545 //! \brief Get the number of attributes of a variable.
546 /*!
547 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
548 */
549 void NC3File::inq_varnatts_impl(const std::string &variable_name, int &result) const {
550 int stat = NC_NOERR;
551
552 if (m_rank == 0) {
553 int varid = get_varid(variable_name);
554
555 if (varid >= NC_GLOBAL) {
556 stat = nc_inq_varnatts(m_file_id, varid, &result);
557 } else {
558 stat = varid; // LCOV_EXCL_LINE
559 }
560 }
561 MPI_Barrier(m_com);
562
563 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
564 check(PISM_ERROR_LOCATION, stat);
565
566 MPI_Bcast(&result, 1, MPI_INT, 0, m_com);
567 }
568
569 //! \brief Finds a variable and sets the "exists" flag.
570 void NC3File::inq_varid_impl(const std::string &variable_name, bool &exists) const {
571 int stat, flag = -1;
572
573 if (m_rank == 0) {
574 stat = nc_inq_varid(m_file_id, variable_name.c_str(), &flag);
575
576 flag = (stat == NC_NOERR) ? 1 : 0;
577 }
578 MPI_Barrier(m_com);
579 MPI_Bcast(&flag, 1, MPI_INT, 0, m_com);
580
581 exists = (flag == 1);
582 }
583
584 void NC3File::inq_varname_impl(unsigned int j, std::string &result) const {
585 int stat = NC_NOERR;
586 std::vector<char> varname(NC_MAX_NAME + 1, 0);
587
588 if (m_rank == 0) {
589 stat = nc_inq_varname(m_file_id, j, varname.data());
590 }
591
592 MPI_Barrier(m_com);
593
594 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
595 MPI_Bcast(varname.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
596
597 check(PISM_ERROR_LOCATION, stat);
598
599 result = varname.data();
600 }
601
602 //! \brief Gets a double attribute.
603 /*!
604 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
605 */
606 void NC3File::get_att_double_impl(const std::string &variable_name,
607 const std::string &att_name,
608 std::vector<double> &result) const {
609 int stat = NC_NOERR, len = 0;
610
611 int varid = get_varid(variable_name);
612
613 // Read and broadcast the attribute length:
614 if (m_rank == 0) {
615 size_t attlen = 0;
616
617 if (varid >= NC_GLOBAL) {
618 stat = nc_inq_attlen(m_file_id, varid, att_name.c_str(), &attlen);
619 } else {
620 stat = varid; // LCOV_EXCL_LINE
621 }
622
623 if (stat == NC_NOERR) {
624 len = static_cast<int>(attlen);
625 } else {
626 len = 0;
627 }
628 }
629 MPI_Bcast(&len, 1, MPI_INT, 0, m_com);
630
631 if (len == 0) {
632 result.clear();
633 return;
634 }
635
636 result.resize(len);
637
638 // Now read data and broadcast stat to see if we succeeded:
639 if (m_rank == 0) {
640 stat = nc_get_att_double(m_file_id, varid, att_name.c_str(), &result[0]);
641 }
642 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
643
644 check(PISM_ERROR_LOCATION, stat);
645
646 // Broadcast data
647 MPI_Bcast(&result[0], len, MPI_DOUBLE, 0, m_com);
648 }
649
650 // Get a text (character array) attribute on rank 0.
651 static int get_att_text(int ncid, int varid, const std::string &att_name,
652 std::string &result) {
653 int stat = NC_NOERR;
654
655 size_t attlen = 0;
656 stat = nc_inq_attlen(ncid, varid, att_name.c_str(), &attlen);
657 if (stat != NC_NOERR) {
658 result = "";
659 return 0;
660 }
661
662 std::vector<char> buffer(attlen + 1, 0);
663 stat = nc_get_att_text(ncid, varid, att_name.c_str(), &buffer[0]);
664 if (stat == NC_NOERR) {
665 result = &buffer[0];
666 } else {
667 result = "";
668 }
669
670 return 0;
671 }
672
673 // Get a string attribute on rank 0. In "string array" attributes array elements are concatenated
674 // using "," as the separator.
675 static int get_att_string(int ncid, int varid, const std::string &att_name,
676 std::string &result) {
677 int stat = NC_NOERR;
678
679 size_t attlen = 0;
680 stat = nc_inq_attlen(ncid, varid, att_name.c_str(), &attlen);
681 if (stat != NC_NOERR) {
682 result = "";
683 return 0;
684 }
685
686 std::vector<char*> buffer(attlen + 1, 0);
687 stat = nc_get_att_string(ncid, varid, att_name.c_str(), &buffer[0]);
688 if (stat == NC_NOERR) {
689 std::vector<std::string> strings(attlen);
690 for (size_t k = 0; k < attlen; ++k) {
691 strings[k] = buffer[k];
692 }
693 result = join(strings, ",");
694 } else {
695 result = "";
696 }
697 stat = nc_free_string(attlen, &buffer[0]);
698
699 return stat;
700 }
701
702
703 //! \brief Gets a text attribute.
704 /*!
705 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
706 */
707 void NC3File::get_att_text_impl(const std::string &variable_name,
708 const std::string &att_name, std::string &result) const {
709 int stat = NC_NOERR;
710
711 // Read and broadcast the attribute length:
712 if (m_rank == 0) {
713
714 int varid = get_varid(variable_name);
715
716 if (varid >= NC_GLOBAL) {
717 nc_type nctype = NC_NAT;
718 stat = nc_inq_atttype(m_file_id, varid, att_name.c_str(), &nctype);
719
720 if (stat == NC_NOERR) {
721 switch (nctype) {
722 case NC_CHAR:
723 stat = pism::io::get_att_text(m_file_id, varid, att_name, result);
724 break;
725 case NC_STRING:
726 stat = pism::io::get_att_string(m_file_id, varid, att_name, result);
727 break;
728 default:
729 result = "";
730 stat = NC_NOERR;
731 }
732 } else if (stat == NC_ENOTATT) {
733 result = "";
734 stat = NC_NOERR;
735 }
736 } else {
737 stat = varid; // LCOV_EXCL_LINE
738 }
739 }
740 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
741 check(PISM_ERROR_LOCATION, stat);
742
743 int len = result.size();
744 MPI_Bcast(&len, 1, MPI_INT, 0, m_com);
745
746 result.resize(len);
747 MPI_Bcast(&result[0], len, MPI_CHAR, 0, m_com);
748 }
749
750
751 //! \brief Writes a double attribute.
752 /*!
753 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
754 */
755 void NC3File::put_att_double_impl(const std::string &variable_name, const std::string &att_name,
756 IO_Type nctype, const std::vector<double> &data) const {
757 int stat = NC_NOERR;
758
759 if (m_rank == 0) {
760 int varid = get_varid(variable_name);
761
762 if (varid >= NC_GLOBAL) {
763 stat = nc_put_att_double(m_file_id, varid, att_name.c_str(),
764 pism_type_to_nc_type(nctype), data.size(), &data[0]);
765 } else {
766 stat = varid; // LCOV_EXCL_LINE
767 }
768 }
769
770 MPI_Barrier(m_com);
771 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
772
773 check(PISM_ERROR_LOCATION, stat);
774 }
775
776
777
778 //! \brief Writes a text attribute.
779 /*!
780 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
781 */
782 void NC3File::put_att_text_impl(const std::string &variable_name, const std::string &att_name,
783 const std::string &value) const {
784 int stat = NC_NOERR;
785
786 if (m_rank == 0) {
787 int varid = get_varid(variable_name);
788
789 if (varid >= NC_GLOBAL) {
790 stat = nc_put_att_text(m_file_id, varid, att_name.c_str(), value.size(), value.c_str());
791 } else {
792 stat = varid; // LCOV_EXCL_LINE
793 }
794 }
795
796 MPI_Barrier(m_com);
797 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
798
799 check(PISM_ERROR_LOCATION, stat);
800 }
801
802 //! \brief Gets the name of a numbered attribute.
803 /*!
804 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
805 */
806 void NC3File::inq_attname_impl(const std::string &variable_name, unsigned int n, std::string &result) const {
807 int stat = NC_NOERR;
808 std::vector<char> name(NC_MAX_NAME + 1, 0);
809
810 if (m_rank == 0) {
811 int varid = get_varid(variable_name);
812
813 if (varid >= NC_GLOBAL) {
814 stat = nc_inq_attname(m_file_id, varid, n, name.data()); check(PISM_ERROR_LOCATION, stat);
815 } else {
816 stat = varid; // LCOV_EXCL_LINE
817 }
818 }
819 MPI_Barrier(m_com);
820 MPI_Bcast(name.data(), NC_MAX_NAME, MPI_CHAR, 0, m_com);
821 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
822
823 check(PISM_ERROR_LOCATION, stat);
824
825 result = name.data();
826 }
827
828 //! \brief Gets the type of an attribute.
829 /*!
830 * Use "PISM_GLOBAL" as the "variable_name" to get the number of global attributes.
831 */
832 void NC3File::inq_atttype_impl(const std::string &variable_name, const std::string &att_name, IO_Type &result) const {
833 int stat, tmp;
834
835 if (m_rank == 0) {
836 int varid = get_varid(variable_name);
837
838 if (varid >= NC_GLOBAL) {
839 // In NetCDF 3.6.x nc_type is an enum; in 4.x it is 'typedef int'.
840 nc_type nctype = NC_NAT;
841 stat = nc_inq_atttype(m_file_id, varid, att_name.c_str(), &nctype);
842 if (stat == NC_ENOTATT) {
843 tmp = NC_NAT;
844 stat = NC_NOERR;
845 } else {
846 tmp = static_cast<int>(nctype);
847 }
848 } else {
849 stat = varid; // LCOV_EXCL_LINE
850 }
851 }
852 MPI_Barrier(m_com);
853 MPI_Bcast(&tmp, 1, MPI_INT, 0, m_com);
854
855 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
856 check(PISM_ERROR_LOCATION, stat);
857
858 result = nc_type_to_pism_type(tmp);
859 }
860
861
862 //! \brief Sets the fill mode.
863 void NC3File::set_fill_impl(int fillmode, int &old_modep) const {
864 int stat = NC_NOERR;
865
866 if (m_rank == 0) {
867 stat = nc_set_fill(m_file_id, fillmode, &old_modep);
868 }
869
870 MPI_Barrier(m_com);
871 MPI_Bcast(&old_modep, 1, MPI_INT, 0, m_com);
872 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
873
874 check(PISM_ERROR_LOCATION, stat);
875 }
876
877 std::string NC3File::get_format() const {
878 int format;
879
880 if (m_rank == 0) {
881 int stat = nc_inq_format(m_file_id, &format); check(PISM_ERROR_LOCATION, stat);
882 }
883 MPI_Barrier(m_com);
884 MPI_Bcast(&format, 1, MPI_INT, 0, m_com);
885
886 switch(format) {
887 case NC_FORMAT_CLASSIC:
888 case NC_FORMAT_64BIT_OFFSET:
889 return "netcdf3";
890 case NC_FORMAT_64BIT_DATA:
891 return "cdf5";
892 case NC_FORMAT_NETCDF4:
893 case NC_FORMAT_NETCDF4_CLASSIC:
894 default:
895 return "netcdf4";
896 }
897 }
898
899 void NC3File::del_att_impl(const std::string &variable_name, const std::string &att_name) const {
900 int stat = NC_NOERR;
901
902 if (m_rank == 0) {
903 int varid = get_varid(variable_name);
904
905 if (varid >= NC_GLOBAL) {
906 stat = nc_del_att(m_file_id, varid, att_name.c_str());
907 }
908 }
909
910 MPI_Barrier(m_com);
911 MPI_Bcast(&stat, 1, MPI_INT, 0, m_com);
912
913 check(PISM_ERROR_LOCATION, stat);
914 }
915
916 /*!
917 * return the varid corresponding to a variable.
918 *
919 * If the value returned is NC_GLOBAL or greater, it is a varid, otherwise it is an error
920 * code.
921 */
922 int NC3File::get_varid(const std::string &variable_name) const {
923 if (variable_name == "PISM_GLOBAL") {
924 return NC_GLOBAL;
925 }
926
927 if (m_rank == 0) {
928 int varid = -2;
929 int stat = nc_inq_varid(m_file_id, variable_name.c_str(), &varid);
930
931 if (stat == NC_NOERR) {
932 return varid;
933 } else {
934 return stat;
935 }
936 } else {
937 return -2; // this value will not be used
938 }
939 }
940
941 } // end of namespace io
942 } // end of namespace pism