Add code changes for the implementation of the low-resolution B

utils/soca/gdas_ens_handler.h

@@ -75,14 +75,22 @@ namespace gdasapp {
     // -----------------------------------------------------------------------------
 
     int execute(const eckit::Configuration & fullConfig) const {
-      // Setup the soca geometry
+      // Setup the native geometry of the ens. members,
+      // currently assumed to be the same as the deterministic
       const eckit::LocalConfiguration geomConfig(fullConfig, "geometry");
       oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl;
       const soca::Geometry geom(geomConfig, this->getComm());
 
+      // Setup the ensemble B (and output) soca geometry
+      oops::Log::info() << "====================== ens B geometry" << std::endl;
+      const std::string outputGeometryKey = fullConfig.has("output geometry")
+                        ? "output geometry"  // keep things backward compatible for now
+                        : "geometry";        // and default to the input geometry
+      const eckit::LocalConfiguration geomOutConfig(fullConfig, outputGeometryKey);
+      const soca::Geometry geomOut(geomOutConfig, this->getComm());
 
       // Initialize the post processing
-      PostProcIncr postProcIncr(fullConfig, geom, this->getComm());
+      PostProcIncr postProcIncr(fullConfig, geom, this->getComm(), geomOut);
 
       oops::Log::info() << "soca increments: " << std::endl
                         << postProcIncr.inputIncrConfig_ << std::endl;
@@ -93,7 +101,9 @@ namespace gdasapp {
         oops::Log::info() << postProcIncr.inputIncrConfig_ << std::endl;
 
         // Assume z* output is the same for the trajectory and the state
-        ensMembers.push_back(postProcIncr.read(i));
+        soca::Increment fullResIncr = postProcIncr.read(i);
+        soca::Increment lowResIncr(geomOut, fullResIncr);  // interp to low resolution geometry
+        ensMembers.push_back(lowResIncr);
       }
 
       // Check if we need to recenter the increment around the deterministic
@@ -103,9 +113,9 @@ namespace gdasapp {
       }
 
       // Compute ensemble moments
-      soca::Increment ensMean(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
-      soca::Increment ensStd(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
-      soca::Increment ensVariance(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
+      soca::Increment ensMean(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_);
+      soca::Increment ensStd(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_);
+      soca::Increment ensVariance(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_);
       gdasapp_ens_utils::ensMoments(ensMembers, ensMean, ensStd, ensVariance);
       oops::Log::info() << "mean: " << ensMean << std::endl;
       oops::Log::info() << "std: " << ensStd << std::endl;
@@ -123,14 +133,15 @@ namespace gdasapp {
 
       // Initialize the trajectories used in the linear variable changes
       const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory");
-      soca::State determTraj(geom, trajConfig);     // deterministic trajectory
+      soca::State determTrajNative(geom, trajConfig);     // deterministic trajectory at full res
+      soca::State determTraj(geomOut, determTrajNative);  // interpolated deterministic trajectory
       soca::State ensMeanTraj(determTraj);          // trajectory defined as the ens. mean
       ensMeanTraj.zero();
       ensMeanTraj += ensMean;
 
       // Compute the recentering increment as the difference between
       // the ensemble mean and the deterministic
-      soca::Increment recenteringIncr(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
+      soca::Increment recenteringIncr(geomOut, postProcIncr.socaIncrVar_, postProcIncr.dt_);
       recenteringIncr.diff(determTraj, ensMeanTraj);
       postProcIncr.setToZero(recenteringIncr);
 
@@ -192,7 +203,7 @@ namespace gdasapp {
 
         // Save total ssh
         oops::Log::info() << "ssh ensemble member "  << i << std::endl;
-        soca::Increment ssh_tmp(geom, socaSshVar, postProcIncr.dt_);
+        soca::Increment ssh_tmp(geomOut, socaSshVar, postProcIncr.dt_);
         ssh_tmp = ensMembers[i];
         sshTotal.push_back(ssh_tmp);
 
@@ -257,9 +268,9 @@ namespace gdasapp {
       }
 
       // Compute ensemble moments for total ssh
-      soca::Increment sshMean(geom, socaSshVar, postProcIncr.dt_);
-      soca::Increment sshStd(geom, socaSshVar, postProcIncr.dt_);
-      soca::Increment sshTotalVariance(geom, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshMean(geomOut, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshStd(geomOut, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshTotalVariance(geomOut, socaSshVar, postProcIncr.dt_);
       gdasapp_ens_utils::ensMoments(sshTotal, sshMean, sshStd, sshTotalVariance);
       oops::Log::info() << "mean ssh total: " << sshMean << std::endl;
       oops::Log::info() << "std ssh total: " << sshStd << std::endl;
@@ -269,7 +280,7 @@ namespace gdasapp {
       // Compute ensemble moments for steric ssh
       sshMean.zero();
       sshStd.zero();
-      soca::Increment sshStericVariance(geom, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshStericVariance(geomOut, socaSshVar, postProcIncr.dt_);
       gdasapp_ens_utils::ensMoments(sshSteric, sshMean, sshStd, sshStericVariance);
       oops::Log::info() << "mean steric ssh: " << sshMean << std::endl;
       oops::Log::info() << "std steric ssh: " << sshStd << std::endl;
@@ -278,8 +289,8 @@ namespace gdasapp {
 
       // Compute ensemble moments for non-steric ssh
       sshMean.zero();
-      soca::Increment sshNonStericVariance(geom, socaSshVar, postProcIncr.dt_);
-      soca::Increment sshNonStericStd(geom, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshNonStericVariance(geomOut, socaSshVar, postProcIncr.dt_);
+      soca::Increment sshNonStericStd(geomOut, socaSshVar, postProcIncr.dt_);
       sshNonStericStd.zero();
       gdasapp_ens_utils::ensMoments(sshNonSteric, sshMean, sshNonStericStd, sshNonStericVariance);
       oops::Log::info() << "mean non-steric ssh: " << sshMean << std::endl;
@@ -308,7 +319,7 @@ namespace gdasapp {
       ensStd.write(bkgErrOutputConfig);
 
       // Explained variance by steric height R=1-SS(non-steric ssh)/SS(total ssh)
-      soca::Increment varianceRatio(geom, socaSshVar, postProcIncr.dt_);
+      soca::Increment varianceRatio(geomOut, socaSshVar, postProcIncr.dt_);
       varianceRatio = sshNonStericVariance;
       atlas::FieldSet varianceRatioFs;
       varianceRatio.toFieldSet(varianceRatioFs);
@@ -319,7 +330,7 @@ namespace gdasapp {
       util::divideFieldSets(varianceRatioFs, sshTotalVarianceFs, sshTotalVarianceFs);
       varianceRatio.fromFieldSet(varianceRatioFs);
 
-      soca::Increment stericExplainedVariance(geom, socaSshVar, postProcIncr.dt_);
+      soca::Increment stericExplainedVariance(geomOut, socaSshVar, postProcIncr.dt_);
       stericExplainedVariance.ones();
       stericExplainedVariance -= varianceRatio;
 

utils/soca/gdas_postprocincr.h

@@ -22,17 +22,38 @@
 
 namespace gdasapp {
 
+// -----------------------------------------------------------------------------
+/*! \class PostProcIncr
+    \brief This class handles the processing of increments in the GDAS application.
+
+    The PostProcIncr class is responsible for managing the configuration and processing
+    of increments, including reading configuration parameters, initializing variables,
+    applying specified variable changes and handling input and output.
+*/
+// -----------------------------------------------------------------------------
 class PostProcIncr {
  public:
   // -----------------------------------------------------------------------------
-  // Constructor
-
+  // Constructors
+
+  /**
+   * @brief Constructor for the PostProcIncr class.
+   *
+   * This constructor initializes the PostProcIncr object using the provided configuration, geometry, and MPI communicator.
+   * It sets up various parameters and configurations required for post-processing increments.
+   *
+   * @param fullConfig The full configuration object containing various settings.
+   * @param geom The native geometry of the increments output.
+   * @param comm The MPI communicator.
+   * @param geomProc The geometry to perform the processing on.
+   */
   PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry& geom,
-               const eckit::mpi::Comm & comm)
+               const eckit::mpi::Comm & comm, const soca::Geometry& geomProc)
     : dt_(getDate(fullConfig)),
       layerVar_(getLayerVar(fullConfig)),
       geom_(geom),
-      Layers_(getLayerThickness(fullConfig, geom)),
+      geomProc_(geomProc),
+      Layers_(getLayerThickness(fullConfig, geom, geomProc)),
       comm_(comm),
       ensSize_(1),
       pattern_() {
@@ -66,6 +87,20 @@ class PostProcIncr {
     }
   }
 
+  /**
+   * @brief Constructor for the PostProcIncr class when the compute/processing geometry
+   * is the same as the native geometry.
+   *
+   * This constructor delegates the initialization to the main constructor
+   *
+   * @param fullConfig The full configuration object.
+   * @param geom The geometry object.
+   * @param comm The MPI communicator.
+   */
+  PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry & geom,
+               const eckit::mpi::Comm & comm)
+      : PostProcIncr(fullConfig, geom, comm, geom) {}
+
   // -----------------------------------------------------------------------------
   // Read ensemble member n
 
@@ -75,7 +110,7 @@ class PostProcIncr {
 
     // initialize the soca increment
     soca::Increment socaIncr(geom_, socaIncrVar_, dt_);
-    eckit::LocalConfiguration memberConfig;  // inputIncrConfig_);
+    eckit::LocalConfiguration memberConfig;
     memberConfig = inputIncrConfig_;
 
     // replace templated string if necessary
@@ -88,8 +123,10 @@ class PostProcIncr {
     oops::Log::debug() << "-------------------- input increment: " << std::endl;
     oops::Log::debug() << socaIncr << std::endl;
 
-    return socaIncr;
-  }
+
+    soca::Increment socaIncrOut(geomProc_, socaIncr);
+    return socaIncrOut;
+    }
 
   // -----------------------------------------------------------------------------
   // Append variable to increment
@@ -174,7 +211,7 @@ class PostProcIncr {
                          const soca::State& xTraj) {
     oops::Log::info() << "==========================================" << std::endl;
     oops::Log::info() << "======      applying specified change of variables" << std::endl;
-    soca::LinearVariableChange lvc(this->geom_, lvcConfig);
+    soca::LinearVariableChange lvc(this->geomProc_, lvcConfig);
     lvc.changeVarTraj(xTraj, socaIncrVar_);
     lvc.changeVarTL(socaIncr, socaIncrVar_);
     oops::Log::info() << " in var change:" << socaIncr << std::endl;
@@ -236,27 +273,15 @@ class PostProcIncr {
   }
   // Read the layer thickness from the relevant background
   soca::Increment getLayerThickness(const eckit::Configuration& fullConfig,
-                                    const soca::Geometry& geom) const {
+                                    const soca::Geometry& geom,
+                                    const soca::Geometry& geomProc) const {
     soca::Increment layerThick(geom, getLayerVar(fullConfig), getDate(fullConfig));
     const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry");
     layerThick.read(vertGeomConfig);
-    oops::Log::debug() << "layerThick: " << std::endl << layerThick << std::endl;
-    return layerThick;
-  }
-  // Initialize the trajectory
-  soca::State getTraj(const eckit::Configuration& fullConfig,
-                      const soca::Geometry& geom) const {
-    if ( fullConfig.has("linear variable change") ) {
-      const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory");
-      soca::State traj(geom, trajConfig);
-      oops::Log::debug() << "traj:" << traj << std::endl;
-      return traj;
-    } else {
-      oops::Variables tmpVar(fullConfig, "layers variable");
-      util::DateTime tmpDate(getDate(fullConfig));
-      soca::State traj(geom, tmpVar, tmpDate);
-      return traj;
-    }
+    soca::Increment layerThickOut(geomProc, layerThick);
+
+    oops::Log::debug() << "layerThickOut: " << std::endl << layerThickOut << std::endl;
+    return layerThickOut;
   }
 
   // -----------------------------------------------------------------------------
@@ -300,7 +325,8 @@ class PostProcIncr {
   util::DateTime dt_;                  // valid date of increment
   oops::Variables layerVar_;           // layer variable
   const soca::Increment Layers_;       // layer thicknesses
-  const soca::Geometry & geom_;
+  const soca::Geometry & geom_;        // Native geometry
+  const soca::Geometry & geomProc_;    // Geometry to perform processing on
   const eckit::mpi::Comm & comm_;
   //  std::vector<eckit::LocalConfiguration> inputIncrConfig_;
   eckit::LocalConfiguration inputIncrConfig_;

utils/soca/gdas_soca_diagb.h

@@ -204,6 +204,14 @@ namespace gdasapp {
       oops::Log::info() << "Background:" << std::endl;
       oops::Log::info() << xb << std::endl;
 
+      // Setup the output soca geometry
+      oops::Log::info() << "====================== output geometry" << std::endl;
+      const std::string outputGeometryKey = fullConfig.has("output geometry")
+                        ? "output geometry"  // keep things backward compatible for now
+                        : "geometry";        // and default to the input geometry
+      const eckit::LocalConfiguration geomOutConfig(fullConfig, outputGeometryKey);
+      const soca::Geometry geomOut(geomOutConfig, this->getComm());
+
       /// Create the mesh connectivity (Copy/paste of Francois's stuff)
       // --------------------------------------------------------------
       // Build edges, then connections between nodes and edges
@@ -485,7 +493,8 @@ namespace gdasapp {
 
       // Save the background error
       const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error");
-      bkgErr.write(bkgErrorConfig);
+      soca::Increment bkgErrOut(geomOut, bkgErr);
+      bkgErrOut.write(bkgErrorConfig);
 
       return 0;
     }