URI:
       tFile.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
       ---
       tFile.cc (22490B)
       ---
            1 // Copyright (C) 2012, 2013, 2014, 2015, 2016, 2017, 2018, 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 <cassert>
           20 #include <cstdio>
           21 #include <memory>
           22 using std::shared_ptr;
           23 
           24 #include <petscvec.h>
           25 
           26 #include "File.hh"
           27 #include "pism/util/IceGrid.hh"
           28 #include "pism/util/pism_utilities.hh"
           29 #include "pism/util/VariableMetadata.hh"
           30 #include "pism/util/ConfigInterface.hh"
           31 #include "pism/util/Time.hh"
           32 #include "NC3File.hh"
           33 
           34 #include "pism/pism_config.hh"
           35 
           36 #if (Pism_USE_PARALLEL_NETCDF4==1)
           37 #include "NC4_Par.hh"
           38 #endif
           39 
           40 #if (Pism_USE_PNETCDF==1)
           41 #include "PNCFile.hh"
           42 #endif
           43 
           44 #if (Pism_USE_PIO==1)
           45 #include "ParallelIO.hh"
           46 #endif
           47 
           48 #include "pism/util/error_handling.hh"
           49 #include "pism/util/io/io_helpers.hh"
           50 
           51 namespace pism {
           52 
           53 struct File::Impl {
           54   MPI_Comm com;
           55   IO_Backend backend;
           56   io::NCFile::Ptr nc;
           57 };
           58 
           59 IO_Backend string_to_backend(const std::string &backend) {
           60   if (backend == "netcdf3") {
           61     return PISM_NETCDF3;
           62   }
           63   if (backend == "netcdf4_parallel") {
           64     return PISM_NETCDF4_PARALLEL;
           65   }
           66   if (backend == "pnetcdf") {
           67     return PISM_PNETCDF;
           68   }
           69   if (backend == "pio_pnetcdf") {
           70     return PISM_PIO_PNETCDF;
           71   }
           72   if (backend == "pio_netcdf") {
           73     return PISM_PIO_NETCDF;
           74   }
           75   if (backend == "pio_netcdf4c") {
           76     return PISM_PIO_NETCDF4C;
           77   }
           78   if (backend == "pio_netcdf4p") {
           79     return PISM_PIO_NETCDF4P;
           80   }
           81   throw RuntimeError::formatted(PISM_ERROR_LOCATION,
           82                                 "unknown or unsupported I/O backend: %s", backend.c_str());
           83 }
           84 
           85 // Chooses the best available I/O backend for reading from 'filename'.
           86 static IO_Backend choose_backend(MPI_Comm com, const std::string &filename) {
           87 
           88   std::string format;
           89   {
           90     // This is the rank-0-only purely-serial mode of accessing NetCDF files, but it
           91     // supports all the kinds of NetCDF, so this is fine.
           92     io::NC3File file(com);
           93 
           94     file.open(filename, PISM_READONLY);
           95     format = file.get_format();
           96     file.close();
           97   }
           98 
           99   if (format == "netcdf4") {
          100 #if (Pism_USE_PARALLEL_NETCDF4==1)
          101     return PISM_NETCDF4_PARALLEL;
          102 #endif
          103   } else {
          104 #if (Pism_USE_PNETCDF==1)
          105     return PISM_PNETCDF;
          106 #endif
          107   }
          108 
          109   // this choice is appropriate for both NetCDF-3 and NetCDF-4
          110   return PISM_NETCDF3;
          111 }
          112 
          113 static io::NCFile::Ptr create_backend(MPI_Comm com, IO_Backend backend, int iosysid) {
          114   int size = 1;
          115   MPI_Comm_size(com, &size);
          116 
          117   if (backend == PISM_NETCDF3) {
          118     return io::NCFile::Ptr(new io::NC3File(com));
          119   }
          120 #if (Pism_USE_PARALLEL_NETCDF4==1)
          121   if (backend == PISM_NETCDF4_PARALLEL) {
          122     return io::NCFile::Ptr(new io::NC4_Par(com));
          123   }
          124 #endif
          125 #if (Pism_USE_PNETCDF==1)
          126   if (backend == PISM_PNETCDF) {
          127     return io::NCFile::Ptr(new io::PNCFile(com));
          128   }
          129 #endif
          130 #if (Pism_USE_PIO==1)
          131   if (backend == PISM_PIO_PNETCDF or
          132       backend == PISM_PIO_NETCDF4P or
          133       backend == PISM_PIO_NETCDF4C or
          134       backend == PISM_PIO_NETCDF) {
          135     if (iosysid == -1) {
          136       throw RuntimeError::formatted(PISM_ERROR_LOCATION,
          137                                     "To use ParallelIO you have to pass iosysid to File");
          138     }
          139     return io::NCFile::Ptr(new io::ParallelIO(com, iosysid, backend));
          140   }
          141 #else
          142   (void) iosysid;               // silence a compiler warning
          143 #endif
          144   throw RuntimeError::formatted(PISM_ERROR_LOCATION,
          145                                 "unknown or unsupported I/O backend: %d", backend);
          146 }
          147 
          148 File::File(MPI_Comm com, const std::string &filename, IO_Backend backend, IO_Mode mode,
          149            int iosysid)
          150   : m_impl(new Impl) {
          151 
          152   if (filename.empty()) {
          153     throw RuntimeError::formatted(PISM_ERROR_LOCATION,
          154                                   "cannot open file: provided file name is empty");
          155   }
          156 
          157   if (backend == PISM_GUESS) {
          158     m_impl->backend = choose_backend(com, filename);
          159   } else {
          160     m_impl->backend = backend;
          161   }
          162 
          163   m_impl->com = com;
          164   m_impl->nc  = create_backend(m_impl->com, m_impl->backend, iosysid);
          165 
          166   this->open(filename, mode);
          167 }
          168 
          169 File::~File() {
          170   if (m_impl->nc and not filename().empty()) {
          171     try {
          172       // a file is still open, so we try to close it
          173       this->close();
          174     } catch (...) {
          175       // don't ever throw from here
          176       handle_fatal_errors(MPI_COMM_SELF);
          177     }
          178   }
          179   delete m_impl;
          180 }
          181 
          182 MPI_Comm File::com() const {
          183   return m_impl->com;
          184 }
          185 
          186 IO_Backend File::backend() const {
          187   return m_impl->backend;
          188 }
          189 
          190 void File::open(const std::string &filename, IO_Mode mode) {
          191   try {
          192 
          193     // opening for reading
          194     if (mode == PISM_READONLY) {
          195 
          196       m_impl->nc->open(filename, mode);
          197 
          198     } else if (mode == PISM_READWRITE_CLOBBER or mode == PISM_READWRITE_MOVE) {
          199 
          200       if (mode == PISM_READWRITE_MOVE) {
          201         io::move_if_exists(m_impl->com, filename);
          202       } else {
          203         io::remove_if_exists(m_impl->com, filename);
          204       }
          205 
          206       m_impl->nc->create(filename);
          207 
          208       int old_fill;
          209       m_impl->nc->set_fill(PISM_NOFILL, old_fill);
          210     } else if (mode == PISM_READWRITE) {                      // mode == PISM_READWRITE
          211 
          212       m_impl->nc->open(filename, mode);
          213 
          214       int old_fill;
          215       m_impl->nc->set_fill(PISM_NOFILL, old_fill);
          216     } else {
          217       throw RuntimeError::formatted(PISM_ERROR_LOCATION, "invalid mode: %d", mode);
          218     }
          219   } catch (RuntimeError &e) {
          220     e.add_context("opening or creating \"" + filename + "\"");
          221     throw;
          222   }
          223 }
          224 
          225 void File::remove_attribute(const std::string &variable_name, const std::string &att_name) const {
          226   try {
          227     m_impl->nc->del_att(variable_name, att_name);
          228   } catch (RuntimeError &e) {
          229     e.add_context("deleting the attribute %s:%s", variable_name.c_str(), att_name.c_str());
          230     throw;
          231   }
          232 }
          233 
          234 void File::close() {
          235   try {
          236     m_impl->nc->close();
          237   } catch (RuntimeError &e) {
          238     e.add_context("closing \"" + filename() + "\"");
          239     throw;
          240   }
          241 }
          242 
          243 void File::sync() const {
          244   try {
          245     m_impl->nc->sync();
          246   } catch (RuntimeError &e) {
          247     e.add_context("synchronizing \"" + filename() + "\"");
          248     throw;
          249   }
          250 }
          251 
          252 void File::redef() const {
          253   try {
          254     m_impl->nc->redef();
          255   } catch (RuntimeError &e) {
          256     e.add_context("switching to define mode; file \"" + filename() + "\"");
          257     throw;
          258   }
          259 }
          260 
          261 
          262 void File::enddef() const {
          263   try {
          264     m_impl->nc->enddef();
          265   } catch (RuntimeError &e) {
          266     e.add_context("switching to data mode; file \"" + filename() + "\"");
          267     throw;
          268   }
          269 }
          270 
          271 std::string File::filename() const {
          272   return m_impl->nc->filename();
          273 }
          274 
          275 
          276 //! \brief Get the number of records. Uses the length of an unlimited dimension.
          277 unsigned int File::nrecords() const {
          278   try {
          279     std::string dim;
          280     m_impl->nc->inq_unlimdim(dim);
          281 
          282     if (dim.empty()) {
          283       return 1;                 // one record
          284     } else {
          285       return this->dimension_length(dim);
          286     }
          287   } catch (RuntimeError &e) {
          288     e.add_context("getting the number of records in file \"" + filename() + "\"");
          289     throw;
          290   }
          291   return 0;                     // LCOV_EXCL_LINE
          292 }
          293 
          294 //! \brief Get the number of records of a certain variable. Uses the length of
          295 //! an associated "time" dimension.
          296 unsigned int File::nrecords(const std::string &name, const std::string &std_name,
          297                             units::System::Ptr unit_system) const {
          298   try {
          299     auto var = find_variable(name, std_name);
          300 
          301     if (not var.exists) {
          302       return 0;
          303     }
          304 
          305     auto dims = dimensions(var.name);
          306 
          307     for (auto d : dims) {
          308       AxisType dimtype = dimension_type(d, unit_system);
          309 
          310       if (dimtype == T_AXIS) {
          311         return this->dimension_length(d);
          312       }
          313     }
          314 
          315     return 1;                   // one record
          316   } catch (RuntimeError &e) {
          317     e.add_context("getting the number of records of variable '%s' ('%s') in '%s'",
          318                   name.c_str(), std_name.c_str(), filename().c_str());
          319     throw;
          320   }
          321   return 0;                     // LCOV_EXCL_LINE
          322 }
          323 
          324 
          325 //! \brief Find a variable using its standard name and/or short name.
          326 /*!
          327  * Sets "result" to the short name found.
          328  */
          329 VariableLookupData File::find_variable(const std::string &short_name, const std::string &std_name) const {
          330   VariableLookupData result;
          331   try {
          332     result.exists = false;
          333 
          334     if (not std_name.empty()) {
          335 
          336       int n_variables = nvariables();
          337 
          338       for (int j = 0; j < n_variables; ++j) {
          339         std::string
          340           name      = variable_name(j),
          341           attribute = read_text_attribute(name, "standard_name");
          342 
          343         if (attribute.empty()) {
          344           continue;
          345         }
          346 
          347         if (attribute == std_name) {
          348           if (not result.exists) {
          349             result.exists = true;
          350             result.found_using_standard_name = true;
          351             result.name = name;
          352           } else {
          353             throw RuntimeError::formatted(PISM_ERROR_LOCATION, "inconsistency in '%s': variables '%s' and '%s'\n"
          354                                           "have the same standard_name (%s)",
          355                                           filename().c_str(), result.name.c_str(),
          356                                           name.c_str(), attribute.c_str());
          357           }
          358         }
          359 
          360       } // end of the for loop
          361     } // end of if (not std_name.empty())
          362 
          363     if (not result.exists) {
          364       m_impl->nc->inq_varid(short_name, result.exists);
          365       if (result.exists) {
          366         result.name = short_name;
          367       } else {
          368         result.name = "";
          369       }
          370 
          371       result.found_using_standard_name = false;
          372     }
          373 
          374   } catch (RuntimeError &e) {
          375     e.add_context("searching for variable '%s' ('%s') in '%s'", short_name.c_str(), std_name.c_str(), filename().c_str());
          376     throw;
          377   }
          378 
          379   return result;
          380 }
          381 
          382 //! \brief Checks if a variable exists.
          383 bool File::find_variable(const std::string &name) const {
          384   try {
          385     bool exists = false;
          386     m_impl->nc->inq_varid(name, exists);
          387     return exists;
          388   } catch (RuntimeError &e) {
          389     e.add_context("searching for variable '%s' in '%s'", name.c_str(), filename().c_str());
          390     throw;
          391   }
          392 }
          393 
          394 std::vector<std::string> File::dimensions(const std::string &variable_name) const {
          395   try {
          396     std::vector<std::string> result;
          397     m_impl->nc->inq_vardimid(variable_name, result);
          398     return result;
          399   } catch (RuntimeError &e) {
          400     e.add_context("getting dimensions of variable '%s' in '%s'", variable_name.c_str(),
          401                   filename().c_str());
          402     throw;
          403   }
          404 }
          405 
          406 
          407 //! \brief Checks if a dimension exists.
          408 bool File::find_dimension(const std::string &name) const {
          409   try {
          410     bool exists = false;
          411     m_impl->nc->inq_dimid(name, exists);
          412     return exists;
          413   } catch (RuntimeError &e) {
          414     e.add_context("searching for dimension '%s' in '%s'", name.c_str(), filename().c_str());
          415     throw;
          416   }
          417 }
          418 
          419 //! \brief Get the length of a dimension.
          420 /*!
          421  * Sets result to 0 if a dimension does not exist.
          422  */
          423 unsigned int File::dimension_length(const std::string &name) const {
          424   try {
          425     if (find_dimension(name)) {
          426       unsigned int result = 0;
          427       m_impl->nc->inq_dimlen(name, result);
          428       return result;
          429     } else {
          430       return 0;
          431     }
          432   } catch (RuntimeError &e) {
          433     e.add_context("getting the length of dimension '%s' in '%s'", name.c_str(), filename().c_str());
          434     throw;
          435   }
          436 }
          437 
          438 //! \brief Get the "type" of a dimension.
          439 /*!
          440  * The "type" is one of X_AXIS, Y_AXIS, Z_AXIS, T_AXIS.
          441  */
          442 AxisType File::dimension_type(const std::string &name,
          443                               units::System::Ptr unit_system) const {
          444   try {
          445     if (not find_variable(name)) {
          446       throw RuntimeError(PISM_ERROR_LOCATION, "coordinate variable " + name + " is missing");
          447     }
          448 
          449     std::string
          450       axis          = read_text_attribute(name, "axis"),
          451       standard_name = read_text_attribute(name, "standard_name"),
          452       units         = read_text_attribute(name, "units");
          453 
          454     // check if it has units compatible with "seconds":
          455 
          456     units::Unit seconds(unit_system, "seconds");
          457     if (units::are_convertible(units::Unit(unit_system, units), seconds)) {
          458       return T_AXIS;
          459     }
          460 
          461     // check the standard_name attribute:
          462     if (standard_name == "time") {
          463       return T_AXIS;
          464     } else if (standard_name == "projection_x_coordinate") {
          465       return X_AXIS;
          466     } else if (standard_name == "projection_y_coordinate") {
          467       return Y_AXIS;
          468     }
          469 
          470     // check the axis attribute:
          471     if (axis == "T" or axis == "t") {
          472       return T_AXIS;
          473     } else if (axis == "X" or axis == "x") {
          474       return X_AXIS;
          475     } else if (axis == "Y" or axis == "y") {
          476       return Y_AXIS;
          477     } else if (axis == "Z" or axis == "z") {
          478       return Z_AXIS;
          479     }
          480 
          481     // check the variable name:
          482     if (name == "x" or name == "X" or
          483         name.find("x") == 0 or name.find("X") == 0) {
          484       return X_AXIS;
          485     } else if (name == "y" or name == "Y" or
          486                name.find("y") == 0 or name.find("Y") == 0) {
          487       return Y_AXIS;
          488     } else if (name == "z" or name == "Z" or
          489                name.find("z") == 0 or name.find("Z") == 0) {
          490       return Z_AXIS;
          491     } else if (name == "t" or name == "T" or name == "time" or
          492                name.find("t") == 0 or name.find("T") == 0) {
          493       return T_AXIS;
          494     }
          495 
          496     // we have no clue:
          497     return UNKNOWN_AXIS;
          498   } catch (RuntimeError &e) {
          499     e.add_context("getting the type of dimension '%s' in '%s'",
          500                   name.c_str(), filename().c_str());
          501     throw;
          502   }
          503   return UNKNOWN_AXIS;          // LCOV_EXCL_LINE
          504 }
          505 
          506 void File::define_dimension(const std::string &name, size_t length) const {
          507   try {
          508     m_impl->nc->def_dim(name, length);
          509   } catch (RuntimeError &e) {
          510     e.add_context("defining dimension '%s' in '%s'", name.c_str(), filename().c_str());
          511     throw;
          512   }
          513 }
          514 
          515 //! \brief Define a variable.
          516 void File::define_variable(const std::string &name, IO_Type nctype, const std::vector<std::string> &dims) const {
          517   try {
          518     m_impl->nc->def_var(name, nctype, dims);
          519 
          520     // FIXME: I need to write and tune chunk_dimensions that would be called below before we use
          521     // this.
          522     //
          523     /*
          524     // if it's not a spatial variable, we're done
          525     if (dims.size() < 2) {
          526       return;
          527     }
          528 
          529     std::vector<size_t> dim_lengths;
          530     for (unsigned int k = 0; k < dims.size(); ++k) {
          531       dim_lengths.push_back(this->dimension_length(dims[k]));
          532     }
          533 
          534     std::vector<size_t> chunk_dims = chunk_dimensions(nctype, dim_lengths);
          535 
          536     m_impl->nc->def_var_chunking(name, chunk_dims);
          537     */
          538 
          539   } catch (RuntimeError &e) {
          540     e.add_context("defining variable '%s' in '%s'", name.c_str(), filename().c_str());
          541     throw;
          542   }
          543 }
          544 
          545 //! \brief Get dimension data (a coordinate variable).
          546 std::vector<double>  File::read_dimension(const std::string &name) const {
          547   try {
          548     if (not find_variable(name)) {
          549       throw RuntimeError(PISM_ERROR_LOCATION, "coordinate variable not found");
          550     }
          551 
          552     unsigned int length = dimension_length(name);
          553 
          554     std::vector<double> result(length);
          555 
          556     read_variable(name, {0}, {length}, result.data());
          557 
          558     return result;
          559   } catch (RuntimeError &e) {
          560     e.add_context("reading dimension '%s' from '%s'", name.c_str(), filename().c_str());
          561     throw;
          562   }
          563 }
          564 
          565 //! \brief Append to the history global attribute.
          566 /*!
          567  * Use write_attribute("PISM_GLOBAL", "history", ...) to overwrite "history".
          568  */
          569 void File::append_history(const std::string &history) const {
          570   try {
          571     std::string old_history = read_text_attribute("PISM_GLOBAL", "history");
          572     redef();
          573     write_attribute("PISM_GLOBAL", "history", history + old_history);
          574   } catch (RuntimeError &e) {
          575     e.add_context("appending to the history attribute in \"" + filename() + "\"");
          576     throw;
          577   }
          578 }
          579 
          580 //! \brief Write a multiple-valued double attribute.
          581 void File::write_attribute(const std::string &var_name, const std::string &att_name, IO_Type nctype,
          582                            const std::vector<double> &values) const {
          583   try {
          584     redef();
          585     m_impl->nc->put_att_double(var_name, att_name, nctype, values);
          586   } catch (RuntimeError &e) {
          587     e.add_context("writing double attribute '%s:%s' in '%s'",
          588                   var_name.c_str(), att_name.c_str(), filename().c_str());
          589     throw;
          590   }
          591 }
          592 
          593 //! \brief Write a text attribute.
          594 void File::write_attribute(const std::string &var_name, const std::string &att_name,
          595                            const std::string &value) const {
          596   try {
          597     redef();
          598     // ensure that the string is null-terminated
          599     m_impl->nc->put_att_text(var_name, att_name, value + "\0");
          600   } catch (RuntimeError &e) {
          601     e.add_context("writing text attribute '%s:%s' in '%s'",
          602                   var_name.c_str(), att_name.c_str(), filename().c_str());
          603     throw;
          604   }
          605 }
          606 
          607 //! \brief Get a double attribute.
          608 std::vector<double> File::read_double_attribute(const std::string &var_name, const std::string &att_name) const {
          609   try {
          610     auto att_type = attribute_type(var_name, att_name);
          611 
          612     // Give an understandable error message if a string attribute was found when
          613     // a number (or a list of numbers) was expected. (We've seen datasets with
          614     // "valid_min" stored as a string...)
          615     if (att_type == PISM_CHAR) {
          616       std::string tmp = read_text_attribute(var_name, att_name);
          617 
          618       throw RuntimeError::formatted(PISM_ERROR_LOCATION,
          619                                     "attribute %s is a string '%s'; expected a number or a list of numbers",
          620                                     att_name.c_str(), tmp.c_str());
          621     } else {
          622       // In this case att_type might be PISM_NAT (if an attribute does not
          623       // exist), but read_double_attribute can handle that.
          624       std::vector<double> result;
          625       m_impl->nc->get_att_double(var_name, att_name, result);
          626       return result;
          627     }
          628   } catch (RuntimeError &e) {
          629     e.add_context("reading double attribute '%s:%s' from '%s'",
          630                   var_name.c_str(), att_name.c_str(), filename().c_str());
          631     throw;
          632   }
          633 }
          634 
          635 //! \brief Get a text attribute.
          636 std::string File::read_text_attribute(const std::string &var_name, const std::string &att_name) const {
          637   try {
          638     auto att_type = attribute_type(var_name, att_name);
          639     if (att_type != PISM_NAT and att_type != PISM_CHAR) {
          640       // attribute exists and is not a string
          641       throw RuntimeError::formatted(PISM_ERROR_LOCATION,
          642                                     "attribute %s is not a string", att_name.c_str());
          643     }
          644 
          645     std::string result;
          646     m_impl->nc->get_att_text(var_name, att_name, result);
          647     return result;
          648   } catch (RuntimeError &e) {
          649     e.add_context("reading text attribute '%s:%s' from %s", var_name.c_str(), att_name.c_str(), filename().c_str());
          650     throw;
          651   }
          652 }
          653 
          654 unsigned int File::nattributes(const std::string &var_name) const {
          655   try {
          656     int result = 0;
          657     m_impl->nc->inq_varnatts(var_name, result);
          658     return result;
          659   } catch (RuntimeError &e) {
          660     e.add_context("getting the number of attributes of variable '%s' in '%s'", var_name.c_str(), filename().c_str());
          661     throw;
          662   }
          663 }
          664 
          665 
          666 std::string File::attribute_name(const std::string &var_name, unsigned int n) const {
          667   try {
          668     std::string result;
          669     m_impl->nc->inq_attname(var_name, n, result);
          670     return result;
          671   } catch (RuntimeError &e) {
          672     e.add_context("getting the name of an attribute of variable '%s' in '%s'", var_name.c_str(), filename().c_str());
          673     throw;
          674   }
          675 }
          676 
          677 
          678 IO_Type File::attribute_type(const std::string &var_name, const std::string &att_name) const {
          679   try {
          680     IO_Type result;
          681     m_impl->nc->inq_atttype(var_name, att_name, result);
          682     return result;
          683   } catch (RuntimeError &e) {
          684     e.add_context("getting the type of an attribute of variable '%s' in '%s'", var_name.c_str(), filename().c_str());
          685     throw;
          686   }
          687 }
          688 
          689 
          690 void File::read_variable(const std::string &variable_name,
          691                            const std::vector<unsigned int> &start,
          692                            const std::vector<unsigned int> &count,
          693                           double *ip) const {
          694   try {
          695     m_impl->nc->get_vara_double(variable_name, start, count, ip);
          696   } catch (RuntimeError &e) {
          697     e.add_context("reading variable '%s' from '%s'", variable_name.c_str(), filename().c_str());
          698     throw;
          699   }
          700 }
          701 
          702 
          703 void File::write_variable(const std::string &variable_name,
          704                           const std::vector<unsigned int> &start,
          705                           const std::vector<unsigned int> &count,
          706                           const double *op) const {
          707   try {
          708     m_impl->nc->put_vara_double(variable_name, start, count, op);
          709   } catch (RuntimeError &e) {
          710     e.add_context("writing variable '%s' to '%s'", variable_name.c_str(), filename().c_str());
          711     throw;
          712   }
          713 }
          714 
          715 
          716 void File::write_distributed_array(const std::string &variable_name,
          717                                    const IceGrid &grid,
          718                                    unsigned int z_count,
          719                                    const double *input) const {
          720   try {
          721     unsigned int t_length = nrecords();
          722     assert(t_length > 0);
          723 
          724     m_impl->nc->write_darray(variable_name, grid, z_count, t_length - 1, input);
          725   } catch (RuntimeError &e) {
          726     e.add_context("writing distributed array '%s' to '%s'",
          727                   variable_name.c_str(), filename().c_str());
          728     throw;
          729   }
          730 }
          731 
          732 
          733 void File::read_variable_transposed(const std::string &variable_name,
          734                                     const std::vector<unsigned int> &start,
          735                                     const std::vector<unsigned int> &count,
          736                                     const std::vector<unsigned int> &imap, double *ip) const {
          737   try {
          738     m_impl->nc->get_varm_double(variable_name, start, count, imap, ip);
          739   } catch (RuntimeError &e) {
          740     e.add_context("reading variable '%s' from '%s'", variable_name.c_str(), filename().c_str());
          741     throw;
          742   }
          743 }
          744 
          745 unsigned int File::nvariables() const {
          746   int n_vars = 0;
          747 
          748   try {
          749     m_impl->nc->inq_nvars(n_vars);
          750   } catch (RuntimeError &e) {
          751     e.add_context("getting the number of variables in '%s'", filename().c_str());
          752     throw;
          753   }
          754 
          755   return n_vars;
          756 }
          757 
          758 std::string File::variable_name(unsigned int id) const {
          759   std::string result;
          760   try {
          761     m_impl->nc->inq_varname(id, result);
          762   } catch (RuntimeError &e) {
          763     e.add_context("getting the name of %d-th variable in '%s'", id, filename().c_str());
          764     throw;
          765   }
          766 
          767   return result;
          768 }
          769 
          770 
          771 } // end of namespace pism