Add u,v to the MOM6 IAU increment if it is not part of the DA control…

… variables (#1417)

Mostly as the title says, this branch will be used for the next batch of
experiments (cp4).
Only a draft because I need to update the #'s for `jcb-algorithm` and
`jcb-gdas`.

# Automated CI tests to run in Global Workflow
<!-- Which Global Workflow CI tests are required to adequately test this
PR? -->
- [ ] atm_jjob <!-- JEDI atm single cycle DA !-->
- [ ] C96C48_ufs_hybatmDA <!-- JEDI atm cycled DA !-->
- [ ] C96C48_hybatmaerosnowDA  <!-- JEDI aero/snow cycled DA !-->
- [x] C48mx500_3DVarAOWCDA <!-- JEDI low-res marine 3DVar cycled DA !-->
- [x] C48mx500_hybAOWCDA <!-- JEDI marine hybrid envar cycled DA !-->
- [ ] C96C48_hybatmDA <!-- GSI atm cycled DA !-->

parm/soca/marine-jcb-base.yaml

@@ -24,7 +24,15 @@ minimizer: RPCG
 final_diagnostics_departures: oman
 final_prints_frequency: PT3H
 number_of_outer_loops: 1
-analysis_variables: [sea_ice_area_fraction, sea_ice_thickness, sea_ice_snow_thickness, sea_water_salinity, sea_water_potential_temperature, eastward_sea_water_velocity, northward_sea_water_velocity, sea_surface_height_above_geoid]
+analysis_variables:
+- sea_ice_area_fraction
+- sea_ice_thickness
+- sea_ice_snow_thickness
+- sea_water_salinity
+- sea_water_potential_temperature
+#- eastward_sea_water_velocity
+#- northward_sea_water_velocity
+- sea_surface_height_above_geoid
 
 # Model things
 # ------------
@@ -49,7 +57,7 @@ marine_forecast_timestep: PT3H
 # Background error model
 background_error_file: '{{berror_model}}'
 marine_number_ensemble_members: '{{NMEM_ENS}}'
-marine_stddev_time: '{{MARINE_WINDOW_MIDDLE}}'
+marine_stddev_time: '{{ MARINE_WINDOW_MIDDLE  | to_isotime }}'
 
 # Observations
 observations: all_observations

utils/soca/gdas_ens_handler.h

@@ -151,14 +151,14 @@ namespace gdasapp {
           oops::Log::info() << "recentered incr " << i << ":" << incr << std::endl;
 
           // Append the vertical geometry (for MOM6 IAU)
-          soca::Increment mom6_incr = postProcIncr.appendLayer(incr);
-          oops::Log::info() << "incr " << i << ":" << mom6_incr << std::endl;
+          postProcIncr.appendLayer(incr);
+          oops::Log::info() << "incr " << i << ":" << incr << std::endl;
 
           // Set variables to zero if specified in the configuration
           postProcIncr.setToZero(incr);
 
           // Save the increments used to initialize the ensemble forecast
-          result = postProcIncr.save(mom6_incr, i+1);
+          result = postProcIncr.save(incr, i+1);
         }
         return result;
       }

utils/soca/gdas_incr_handler.h

@@ -35,6 +35,10 @@ namespace gdasapp {
       oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl;
       const soca::Geometry geom(geomConfig, this->getComm());
 
+      // Check that we are using at least 2 mpi tasks
+      if (this->getComm().size() < 2) {
+        throw eckit::BadValue("This application requires at least 2 MPI tasks", Here());
+      }
 
       // Initialize the post processing
       PostProcIncr postProcIncr(fullConfig, geom, this->getComm());
@@ -50,18 +54,20 @@ namespace gdasapp {
         // Read increment from file
         soca::Increment incr = postProcIncr.read(i);
 
-        // At the very minimum, we run this script to append the layers state, so do that!
-        soca::Increment incrWithLayer = postProcIncr.appendLayer(incr);
-        oops::Log::debug() << "========= after appending layer:" << std::endl;
-        oops::Log::debug() << incrWithLayer << std::endl;
+        // Append variables to the increment
+        oops::Variables extraVars(postProcIncr.socaZeroIncrVar_);
+        extraVars += postProcIncr.layerVar_;
+        soca::Increment incr_mom6 = postProcIncr.appendVar(incr, extraVars);
 
         // Zero out specified fields
-        postProcIncr.setToZero(incrWithLayer);
+        postProcIncr.setToZero(incr_mom6);
+        oops::Log::debug() << "========= after appending variables:" << std::endl;
+        oops::Log::debug() << incr_mom6 << std::endl;
 
         // Save final increment
-        result = postProcIncr.save(incrWithLayer, i);
+        result = postProcIncr.save(incr_mom6, i);
         oops::Log::debug() << "========= after appending layer and after saving:" << std::endl;
-        oops::Log::debug() << incrWithLayer << std::endl;
+        oops::Log::debug() << incr_mom6 << std::endl;
       }
       return result;
     }

utils/soca/gdas_postprocincr.h

@@ -92,42 +92,49 @@ class PostProcIncr {
   }
 
   // -----------------------------------------------------------------------------
-  // Append layer thicknesses to increment
-  // TODO(guillaume): There's got to be a better way to append a variable.
-  soca::Increment appendLayer(const soca::Increment& socaIncr) {
+  // Append variable to increment
+  soca::Increment appendVar(const soca::Increment& socaIncr, const oops::Variables varToAppend) {
     oops::Log::info() << "==========================================" << std::endl;
-    oops::Log::info() << "======  Append Layers" << std::endl;
+    oops::Log::info() << "======  Append " << varToAppend << std::endl;
 
     // make a copy of the input increment
     soca::Increment socaIncrOut(socaIncr);
 
     // concatenate variables
     oops::Variables outputIncrVar(socaIncrVar_);
-    outputIncrVar += layerVar_;
+    outputIncrVar += varToAppend;
     oops::Log::debug() << "-------------------- outputIncrVar: " << std::endl;
     oops::Log::debug() << outputIncrVar << std::endl;
 
-    // append layer variable to the soca increment
+    // append variable to the soca increment
     atlas::FieldSet socaIncrFs;
     socaIncrOut.toFieldSet(socaIncrFs);
     socaIncrOut.updateFields(outputIncrVar);
 
     // pad layer increment with zeros
-    soca::Increment layerThick(Layers_);
-    atlas::FieldSet socaLayerThickFs;
-    oops::Log::debug() << "-------------------- thickness field: " << std::endl;
-    oops::Log::debug() << layerThick << std::endl;
-    layerThick.toFieldSet(socaLayerThickFs);
-    layerThick.updateFields(outputIncrVar);
-
-    // append layer thinckness to increment
-    socaIncrOut += layerThick;
+    soca::Increment incrToAppend(Layers_);  // Assumes that Layers_ contains varToAppend
+    atlas::FieldSet incrToAppendFs;
+    oops::Log::debug() << "-------------------- incrToAppend fields: " << std::endl;
+    oops::Log::debug() << incrToAppend << std::endl;
+    incrToAppend.toFieldSet(incrToAppendFs);
+    incrToAppend.updateFields(outputIncrVar);
+
+    // append variables to increment
+    socaIncrOut += incrToAppend;
     oops::Log::debug() << "-------------------- output increment: " << std::endl;
     oops::Log::debug() << socaIncrOut << std::endl;
 
     return socaIncrOut;
   }
 
+  // -----------------------------------------------------------------------------
+  // Append layer thicknesses to increment
+  soca::Increment appendLayer(soca::Increment& socaIncr) {
+    // Append layer thicknesses to the increment
+    soca::Increment socaIncrOut = appendVar(socaIncr, layerVar_);
+    return socaIncrOut;
+  }
+
   // -----------------------------------------------------------------------------
   // Set specified variables to 0
 
@@ -138,18 +145,18 @@ class PostProcIncr {
       return;
     }
     oops::Log::info() << "======      Set specified increment variables to 0.0" << std::endl;
-    std::cout << socaZeroIncrVar_ << std::endl;
 
     atlas::FieldSet socaIncrFs;
     socaIncr.toFieldSet(socaIncrFs);
 
+
     for (auto & field : socaIncrFs) {
       // only works if rank is 2
       ASSERT(field.rank() == 2);
 
       // Set variable to zero
       if (socaZeroIncrVar_.has(field.name())) {
-        std::cout << "setting " << field.name() << " to 0" << std::endl;
+        oops::Log::info() << "setting " << field.name() << " to 0" << std::endl;
         auto view = atlas::array::make_view<double, 2>(field);
         view.assign(0.0);
       }
@@ -170,7 +177,7 @@ class PostProcIncr {
     soca::LinearVariableChange lvc(this->geom_, lvcConfig);
     lvc.changeVarTraj(xTraj, socaIncrVar_);
     lvc.changeVarTL(socaIncr, socaIncrVar_);
-    oops::Log::info() << "$%^#& in var change:" << socaIncr << std::endl;
+    oops::Log::info() << " in var change:" << socaIncr << std::endl;
   }
 
   // -----------------------------------------------------------------------------
@@ -179,6 +186,7 @@ class PostProcIncr {
   int save(soca::Increment& socaIncr, int ensMem = 1) {
     oops::Log::info() << "==========================================" << std::endl;
     oops::Log::info() << "-------------------- save increment: " << std::endl;
+    oops::Log::info() << socaIncr << std::endl;
     socaIncr.write(outputIncrConfig_);
 
     // wait for everybody to be done