Add code changes for the implementation of the low-resolution B (#1441)

As the title says.

## Dependencies
- NOAA-EMC/jcb-gdas#66
- NOAA-EMC/global-workflow#3238
- Update of the static files in glopara:
NOAA-EMC/global-workflow#3239

## Work done in this PR
Interp to low-resolution:
- [ ] variance partitioning
- [ ] preparation of the perturbations for the envar
- [ ] hybrid weights
- [ ] ... 

## Testing
Visual check showing the variance explained by the steric height (30
members)

![steic-explained-variance](https://github.com/user-attachments/assets/274c3fdb-94c3-42fc-ad6e-90867b39a778)
left is native resolution, right is on the 1/2 deg grid.

**"Science" note:** The explained variance above is almost the opposite
to what I would have expected, and what we saw in the previous offline
ensemble, or even what we see in tendencies. Check if there's a bug in
the explained variance calculation.

# 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 !-->

Author: guillaumevernieres

parm/jcb-gdas

@@ -1,22 +1,16 @@
  &MOM_input_nml
         output_directory = './',
-        input_filename = 'r'
+        input_filename = 'r',
+        {% if not simple_geom %}
         restart_input_dir = 'INPUT/',
         restart_output_dir = 'RESTART/',
-        parameter_filename = 'MOM_input' /
+        {% endif %}
+        parameter_filename = '{{ mom_input }}'
+ /
 
  &diag_manager_nml
  /
 
- &ocean_solo_nml
-            months = 0
-            days   = 1
-            date_init = {{ date_init }},
-            hours = 0
-            minutes = 0
-            seconds = 0
-            calendar = 'NOLEAP' /
-
  &fms_io_nml
       max_files_w=100
       checksum_required=.false.

parm/soca/fms/input.nml.j2

@@ -1,9 +1,11 @@
 mkdir:
 - "{{ DATA }}/INPUT"
+- "{{ DATA }}/anl_geom"
 ######################################
 # fix files to copy
 ######################################
 copy:
+# Deterministic resource files
 - ["{{ SOCA_INPUT_FIX_DIR }}/rossrad.nc", "{{ DATA }}/rossrad.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/field_table", "{{ DATA }}/field_table"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/diag_table", "{{ DATA }}/diag_table"]
@@ -23,3 +25,6 @@ copy:
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/geothermal_davies2013_v1.nc", "{{ DATA }}/INPUT/geothermal_davies2013_v1.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/interpolate_zgrid_46L.nc", "{{ DATA }}/INPUT/interpolate_zgrid_46L.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/C384_mosaic.nc", "{{ DATA }}/INPUT/C384_mosaic.nc"]
+# Analysis resource files
+- ["{{ SOCA_ANL_GEOM }}/soca_gridspec.nc", "{{ DATA }}/anl_geom/soca_gridspec.nc"]
+- ["{{ SOCA_ANL_GEOM }}/MOM_input", "{{ DATA }}/anl_geom/MOM_input"]

parm/soca/soca_fix_stage_025.yaml.j2

@@ -1,9 +1,11 @@
 mkdir:
 - "{{ DATA }}/INPUT"
+- "{{ DATA }}/anl_geom"
 ######################################
 # fix files to copy
 ######################################
 copy:
+# Deterministic resource files
 - ["{{ SOCA_INPUT_FIX_DIR }}/rossrad.nc", "{{ DATA }}/rossrad.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/field_table", "{{ DATA }}/field_table"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/diag_table", "{{ DATA }}/diag_table"]
@@ -23,3 +25,6 @@ copy:
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/topo_edits_011818.nc", "{{ DATA }}/INPUT/topo_edits_011818.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/MOM_channels_SPEAR", "{{ DATA }}/INPUT/MOM_channels_SPEAR"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/seawifs_1998-2006_smoothed_2X.nc", "{{ DATA }}/INPUT/seawifs_1998-2006_smoothed_2X.nc"]
+# Analysis resource files
+- ["{{ SOCA_ANL_GEOM }}/soca_gridspec.nc", "{{ DATA }}/anl_geom/soca_gridspec.nc"]
+- ["{{ SOCA_ANL_GEOM }}/MOM_input", "{{ DATA }}/anl_geom/MOM_input"]

parm/soca/soca_fix_stage_100.yaml.j2

@@ -1,9 +1,11 @@
 mkdir:
 - "{{ DATA }}/INPUT"
+- "{{ DATA }}/anl_geom"
 ######################################
 # fix files to copy
 ######################################
 copy:
+# Deterministic resource files
 - ["{{ SOCA_INPUT_FIX_DIR }}/rossrad.nc", "{{ DATA }}/rossrad.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/field_table", "{{ DATA }}/field_table"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/diag_table", "{{ DATA }}/diag_table"]
@@ -15,3 +17,6 @@ copy:
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/ocean_hgrid.nc", "{{ DATA }}/INPUT/ocean_hgrid.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/ocean_mosaic.nc", "{{ DATA }}/INPUT/ocean_mosaic.nc"]
 - ["{{ SOCA_INPUT_FIX_DIR }}/INPUT/ocean_topog.nc", "{{ DATA }}/INPUT/ocean_topog.nc"]
+# Analysis resource files
+- ["{{ SOCA_ANL_GEOM }}/soca_gridspec.nc", "{{ DATA }}/anl_geom/soca_gridspec.nc"]
+- ["{{ SOCA_ANL_GEOM }}/MOM_input", "{{ DATA }}/anl_geom/MOM_input"]

parm/soca/soca_fix_stage_500.yaml.j2

@@ -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_ens_handler.h

@@ -2,6 +2,7 @@
 
 #include <iostream>
 #include <string>
+#include <vector>
 
 #include "eckit/config/LocalConfiguration.h"
 
@@ -35,6 +36,12 @@ namespace gdasapp {
       oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl;
       const soca::Geometry geom(geomConfig, this->getComm());
 
+      // Domains to save
+      std::vector<std::string> domains = {"ocn", "ice"};
+      if (fullConfig.has("domains")) {
+        fullConfig.get("domains", domains);
+      }
+
       // 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());
@@ -65,7 +72,7 @@ namespace gdasapp {
         oops::Log::debug() << incr_mom6 << std::endl;
 
         // Save final increment
-        result = postProcIncr.save(incr_mom6, i);
+        result = postProcIncr.save(incr_mom6, i, domains);
         oops::Log::debug() << "========= after appending layer and after saving:" << std::endl;
         oops::Log::debug() << incr_mom6 << std::endl;
       }