Add method to write increment fields (#89)

* Convert state field write to a generic nemo field writer. Add code to write increments.
* Detect fieldtype from the atlas field when writing.
* Implement test for increment writing and also nemo field set write.
* Setup increment configuration parameters.
* move StateIOUtils to utilities/IOUtils (#105)

---------

Co-authored-by: Toby Searle <[email protected]>

src/orca-jedi/CMakeLists.txt

@@ -21,8 +21,6 @@ model/ModelBiasIncrement.h
 model/Model.h
 state/State.h
 state/State.cc
-state/StateIOUtils.h
-state/StateIOUtils.cc
 nemo_io/NemoFieldReader.h
 nemo_io/NemoFieldReader.cc
 nemo_io/NemoFieldWriter.h
@@ -35,6 +33,8 @@ nemo_io/WriteServer.cc
 utilities/OrcaModelTraits.h
 utilities/Types.h
 utilities/Types.cc
+utilities/IOUtils.h
+utilities/IOUtils.cc
 variablechanges/VariableChangeParameters.h
 variablechanges/VariableChange.h
 )

src/orca-jedi/increment/Increment.cc

@@ -28,10 +28,11 @@
 
 #include "ufo/GeoVaLs.h"
 
+#include "orca-jedi/utilities/IOUtils.h"
 #include "orca-jedi/geometry/Geometry.h"
 #include "orca-jedi/state/State.h"
-#include "orca-jedi/state/StateIOUtils.h"
 #include "orca-jedi/increment/Increment.h"
+#include "orca-jedi/increment/IncrementParameters.h"
 
 #include "atlas/mesh.h"
 #include "atlas-orca/grid/OrcaGrid.h"
@@ -397,11 +398,16 @@ void Increment::random() {
 }
 
 /// \brief Apply Dirac delta functions to configuration specified points.
-void Increment::dirac(const eckit::Configuration & conf) {
+void Increment::dirac(const eckit::Configuration & config) {
+  dirac(oops::validateAndDeserialize<OrcaDiracParameters>(config));
+}
+
+/// \brief Apply Dirac delta functions to params specified points.
+void Increment::dirac(const OrcaDiracParameters & params) {
 // Adding a delta function at points specified by ixdir, iydir, izdir
-  const std::vector<int> & ixdir = conf.getIntVector("ixdir");
-  const std::vector<int> & iydir = conf.getIntVector("iydir");
-  const std::vector<int> & izdir = conf.getIntVector("izdir");
+  const std::vector<int> & ixdir = params.ixdir;
+  const std::vector<int> & iydir = params.iydir;
+  const std::vector<int> & izdir = params.izdir;
 
   ASSERT(ixdir.size() == iydir.size() && ixdir.size() == izdir.size());
   int ndir = ixdir.size();
@@ -526,10 +532,25 @@ void Increment::read(const eckit::Configuration & conf) {
   throw eckit::NotImplemented(err_message, Here());
 }
 
-void Increment::write(const eckit::Configuration & conf) const {
-  std::string err_message =
-      "orcamodel::Increment::write not implemented";
-  throw eckit::NotImplemented(err_message, Here());
+/// \brief write out increments fields to a file using params specified filename.
+void Increment::write(const OrcaIncrementParameters & params) const {
+  oops::Log::debug() << "orcamodel::increment::write" << std::endl;
+
+  std::string output_filename = params.nemoFieldFile.value();
+  if (output_filename == "")
+    throw eckit::BadValue(std::string("orcamodel::writeIncrementFieldsToFile:: ")
+        + "file name not specified", Here());
+
+  auto nemo_field_path = eckit::PathName(output_filename);
+  oops::Log::debug() << "Increment::write to filename "
+                     << nemo_field_path << std::endl;
+
+  writeFieldsToFile(nemo_field_path, *geom_, time_, incrementFields_);
+}
+
+/// \brief write out increments fields to a file using config specified filename.
+void Increment::write(const eckit::Configuration & config) const {
+  write(oops::validateAndDeserialize<OrcaIncrementParameters>(config));
 }
 
 /// \brief Print some basic information about the self increment object.

src/orca-jedi/increment/Increment.h

@@ -21,6 +21,7 @@
 
 #include "orca-jedi/geometry/Geometry.h"
 #include "orca-jedi/state/State.h"
+#include "orca-jedi/increment/IncrementParameters.h"
 
 #include "eckit/exception/Exceptions.h"
 
@@ -78,6 +79,7 @@ class Increment : public util::Printable,
   double dot_product_with(const Increment &) const;
   void schur_product_with(const Increment &);
   void random();
+  void dirac(const OrcaDiracParameters &);
   void dirac(const eckit::Configuration &);
 
 /// ATLAS
@@ -96,7 +98,9 @@ class Increment : public util::Printable,
   };
 
   void read(const eckit::Configuration &);
+  void write(const OrcaIncrementParameters &) const;
   void write(const eckit::Configuration &) const;
+
   void print(std::ostream & os) const override;
   struct stats stats(const std::string & field_name) const;
   double norm() const;

src/orca-jedi/increment/IncrementParameters.h

@@ -0,0 +1,35 @@
+/*
+ * (C) British Crown Copyright 2024 Met Office
+ */
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "eckit/exception/Exceptions.h"
+#include "oops/util/DateTime.h"
+#include "oops/base/Variables.h"
+#include "oops/base/ParameterTraitsVariables.h"
+#include "oops/util/parameters/Parameter.h"
+#include "oops/util/parameters/RequiredParameter.h"
+
+namespace orcamodel {
+
+class OrcaIncrementParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(OrcaIncrementParameters, oops::Parameters)
+
+ public:
+  oops::RequiredParameter<std::string> nemoFieldFile{"output path", this};
+  oops::RequiredParameter<util::DateTime> date{"date", this};
+};
+
+class OrcaDiracParameters : public oops::Parameters {
+  OOPS_CONCRETE_PARAMETERS(OrcaDiracParameters, oops::Parameters)
+
+ public:
+  oops::RequiredParameter<std::vector<int>> ixdir{"x indices", this};
+  oops::RequiredParameter<std::vector<int>> iydir{"y indices", this};
+  oops::RequiredParameter<std::vector<int>> izdir{"z indices", this};
+};
+}  //  namespace orcamodel

src/orca-jedi/state/State.cc

@@ -39,7 +39,7 @@
 #include "orca-jedi/variablechanges/VariableChange.h"
 #include "orca-jedi/increment/Increment.h"
 #include "orca-jedi/state/State.h"
-#include "orca-jedi/state/StateIOUtils.h"
+#include "orca-jedi/utilities/IOUtils.h"
 #include "orca-jedi/utilities/Types.h"
 
 
@@ -78,9 +78,12 @@ State::State(const Geometry & geom,
   if (params_.analyticInit.value().value_or(false)) {
     this->analytic_init(*geom_);
   } else {
-    readFieldsFromFile(params_, *geom_, validTime(), "background",
-       stateFields_);
-    readFieldsFromFile(params_, *geom_, validTime(), "background variance",
+    std::string nemo_file_name;
+    nemo_file_name = params.nemoFieldFile.value();
+    readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background",
+         stateFields_);
+    nemo_file_name = params.varianceFieldFile.value().value_or("");
+    readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background variance",
        stateFields_);
   }
   geom_->log_status();
@@ -185,7 +188,8 @@ void State::read(const OrcaStateParameters & params) {
     throw eckit::UserError(msg.str(), Here());
   }
 
-  readFieldsFromFile(params, *geom_, validTime(), "background",
+  std::string nemo_file_name = params.nemoFieldFile.value();
+  readFieldsFromFile(nemo_file_name, *geom_, validTime(), "background",
       stateFields_);
   oops::Log::trace() << "State(ORCA)::read done" << std::endl;
 }
@@ -225,7 +229,8 @@ void State::setupStateFields() {
 
 void State::write(const OrcaStateParameters & params) const {
   oops::Log::trace() << "State(ORCA)::write starting" << std::endl;
-  writeFieldsToFile(params, *geom_, validTime(), stateFields_);
+  writeFieldsToFile(params.outputNemoFieldFile.value().value_or(""), *geom_, validTime(),
+      stateFields_);
 }
 
 void State::write(const eckit::Configuration & config) const {

src/orca-jedi/state/StateIOUtils.cc → src/orca-jedi/utilities/IOUtils.cc

@@ -2,7 +2,7 @@
  * (C) British Crown Copyright 2024 Met Office
  */
 
-#include "orca-jedi/state/StateIOUtils.h"
+#include "orca-jedi/utilities/IOUtils.h"
 
 #include <sstream>
 #include <vector>
@@ -15,31 +15,28 @@
 #include "atlas/field/MissingValue.h"
 
 #include "orca-jedi/geometry/Geometry.h"
-#include "orca-jedi/state/State.h"
 #include "orca-jedi/nemo_io/WriteServer.h"
+#include "orca-jedi/utilities/Types.h"
 
 #define NEMO_FILL_TOL 1e-6
 
 namespace orcamodel {
 
+/// \brief Read atlas fields to a file.
+/// \param nemo_file_name The name of the netCDF file to read.
+/// \param geom The orcamodel::Geometry of the model data.
+/// \param valid_date The datetime of data to write to the file.
+/// \param variable_type The type of variable (e.g "background" for background data).
+/// \param fs The atlas FieldSet collection of fields to write.
 void readFieldsFromFile(
-  const OrcaStateParameters & params,
+  const std::string & nemo_file_name,
   const Geometry & geom,
   const util::DateTime & valid_date,
   const std::string & variable_type,
   atlas::FieldSet & fs) {
     oops::Log::trace() << "orcamodel::readFieldsFromFile:: start for valid_date"
                        << " " << valid_date << std::endl;
 
-    // Open Nemo Feedback file
-    std::string nemo_file_name;
-    if (variable_type == "background") {
-      nemo_file_name = params.nemoFieldFile.value();
-    }
-    if (variable_type == "background variance") {
-      nemo_file_name = params.varianceFieldFile.value().value_or("");
-    }
-
     auto nemo_field_path = eckit::PathName(nemo_file_name);
     oops::Log::debug() << "orcamodel::readFieldsFromFile:: nemo_field_path "
                        << nemo_field_path << std::endl;
@@ -77,7 +74,7 @@ void readFieldsFromFile(
         };
         ApplyForFieldType(populate,
                           geom.fieldPrecision(fieldName),
-                          std::string("State(ORCA)::readFieldsFromFile '")
+                          std::string("orcamodel::readFieldsFromFile '")
                             + nemoName + "' field type not recognised");
         // Add a halo exchange following read to fill out halo points
         geom.functionSpace().haloExchange(field);
@@ -125,21 +122,29 @@ template void populateField<float>(
   ReadServer & nemo_reader,
   atlas::Field & field);
 
+/// \brief write out atlas fields to a file.
+/// \param output_file_name The name of the netCDF file.
+/// \param geom The orcamodel::Geometry of the model data.
+/// \param valid_date The datetime of data to write to the file.
+/// \param fs The atlas FieldSet collection of fields to write.
 void writeFieldsToFile(
-  const OrcaStateParameters & params,
+  const std::string & output_file_name,
   const Geometry & geom,
   const util::DateTime & valid_date,
   const atlas::FieldSet & fs) {
-    oops::Log::trace() << "orcamodel::writeFieldsToFile:: start for valid_date"
+    oops::Log::trace() << "orcamodel::writeStateFieldsToFile:: start for valid_date"
                        << " " << valid_date << std::endl;
 
-    std::string output_filename =
-      params.outputNemoFieldFile.value().value_or("");
-    if (output_filename == "")
-      throw eckit::BadValue(std::string("orcamodel::writeFieldsToFile:: ")
+    if (output_file_name == "")
+      throw eckit::BadValue(std::string("orcamodel::writeStateFieldsToFile:: ")
           + "file name not specified", Here());
 
+    auto nemo_field_path = eckit::PathName(output_file_name);
+    oops::Log::debug() << "orcamodel::writeStateFieldsToFile:: "
+                       << nemo_field_path << std::endl;
+
     std::map<std::string, std::string> varCoordTypeMap;
+
     {
       const oops::Variables vars = geom.variables();
       const std::vector<std::string> coordSpaces =
@@ -148,7 +153,6 @@ void writeFieldsToFile(
         varCoordTypeMap[vars[i].name()] = coordSpaces[i];
     }
 
-    auto nemo_field_path = eckit::PathName(output_filename);
     oops::Log::debug() << "orcamodel::writeFieldsToFile:: nemo_field_path "
                        << nemo_field_path << std::endl;
     std::vector<util::DateTime> datetimes = {valid_date};
@@ -170,10 +174,13 @@ void writeFieldsToFile(
           writer.write_vol_var<T>(nemoName, 0, field_mv, field_view);
         }
       };
+
+      oops::Log::debug() << "field " << fieldName
+                         << " data type " << field.datatype().str() << std::endl;
       ApplyForFieldType(write,
-                        geom.fieldPrecision(fieldName),
-                        std::string("State(ORCA)::writeFieldsToFile '")
-                          + nemoName + "' field type not recognised.");
+                        field.datatype(),
+                        std::string("orcamodel::writeFieldsToFile '")
+                          + nemoName + "' field data type not recognised.");
     }
 }
 

src/orca-jedi/state/StateIOUtils.h → src/orca-jedi/utilities/IOUtils.h

@@ -10,21 +10,18 @@
 #include "atlas/field.h"
 
 #include "orca-jedi/geometry/Geometry.h"
-#include "orca-jedi/state/State.h"
-#include "orca-jedi/state/StateParameters.h"
-#include "orca-jedi/utilities/Types.h"
 #include "orca-jedi/nemo_io/ReadServer.h"
 
 namespace orcamodel {
 
 void readFieldsFromFile(
-  const OrcaStateParameters & params,
+  const std::string & nemo_file_name,
   const Geometry & geom,
   const util::DateTime & valid_date,
   const std::string & variable_type,
   atlas::FieldSet & fs);
 void writeFieldsToFile(
-  const OrcaStateParameters & params,
+  const std::string & output_file_name,
   const Geometry & geom,
   const util::DateTime & valid_date,
   const atlas::FieldSet & fs);

src/orca-jedi/utilities/Types.h

@@ -10,6 +10,8 @@
 
 #include "eckit/exception/Exceptions.h"
 
+#include "atlas/array/DataType.h"
+
 namespace orcamodel {
 
 //// \brief Enum type for obs variable data types
@@ -31,6 +33,19 @@ void ApplyForFieldType(const Functor& functor, FieldDType field_type,
   }
 }
 
+/// \brief Apply a function for a given Atlas DataType
+template<typename Functor>
+void ApplyForFieldType(const Functor& functor, atlas::DataType datatype,
+     const std::string& error_message) {
+  if (datatype.str() == "real32") {
+    functor(float{});
+  } else if (datatype.str() == "real64") {
+    functor(double{});
+  } else {
+    throw eckit::BadParameter(error_message);
+  }
+}
+
 // create a mapping between C++ types and NetCDF type objects
 template <typename T>
 struct NetCDFTypeMap {