Add sea ice recentering

parm/soca/berror/soca_ensrecenter.yaml

@@ -14,6 +14,9 @@ increment variables:
 - sea_water_salinity
 - eastward_sea_water_velocity
 - northward_sea_water_velocity
+- sea_ice_area_fraction
+- sea_ice_thickness
+- sea_ice_snow_thickness
 set increment variables to zero:
 - eastward_sea_water_velocity
 - northward_sea_water_velocity
@@ -26,6 +29,57 @@ vertical geometry:
 
 add recentering increment: true
 recentering around deterministic: true
+sea ice recenter: true
+sea ice analysis:
+  pattern: '%mem%'
+  date: '{{ MARINE_WINDOW_END_ISO }}'
+  ocn_filename: "ocean.%mem%.nc"
+  ice_filename: "ice.%mem%.nc"
+  read_from_file: 1
+  basename: {{ ENSPERT_RELPATH }}/ens/
+  state variables:
+  - sea_water_potential_temperature
+  - sea_water_salinity
+  - eastward_sea_water_velocity
+  - northward_sea_water_velocity
+  - sea_water_cell_thickness
+  - sea_water_depth
+  - ocean_mixed_layer_thickness
+  - sea_ice_area_fraction
+  - sea_ice_thickness
+  - sea_ice_snow_thickness
+
+sea ice variable change:
+  variable change name: Soca2Cice
+  arctic:
+    seaice edge: 0.4
+    shuffle: true
+    rescale prior:
+      rescale: true
+      min hice: 0.5
+      min hsno: 0.1
+  antarctic:
+    seaice edge: 0.4
+    shuffle: true
+    rescale prior:
+      rescale: true
+      min hice: 0.5
+      min hsno: 0.1
+  cice background state:
+    pattern: '%mem%'
+    restart: ens/cice_model.res.%mem%.nc
+    ncat: 5
+    ice_lev: 7
+    sno_lev: 1
+  cice output:
+    pattern: '%mem%'
+    restart: ens/cice_model.res.output.%mem%.nc
+  output variables:
+  - sea_water_potential_temperature
+  - sea_water_salinity
+  - sea_ice_area_fraction
+  - sea_ice_thickness
+  - sea_ice_snow_thickness
 
 soca increments:  # Could also be states, but they are read as increments
   number of increments: {{ NMEM_ENS }}
@@ -34,8 +88,7 @@ soca increments:  # Could also be states, but they are read as increments
     date: '{{ MARINE_WINDOW_END_ISO }}'
     basename: ./ens/
     ocn_filename: 'ocean.%mem%.nc'
-    # TODO(AFE) fix ice recentering in cycled da
-    # ice_filename: 'ice.%mem%.nc'
+    ice_filename: 'ice.%mem%.nc'
     read_from_file: 3
 
 trajectory:
@@ -48,11 +101,13 @@ trajectory:
   - sea_water_cell_thickness
   - sea_water_depth
   - ocean_mixed_layer_thickness
+  - sea_ice_area_fraction
+  - sea_ice_thickness
+  - sea_ice_snow_thickness
   date: '{{ MARINE_WINDOW_END_ISO }}'
   basename: ./bkg/
   ocn_filename: ocean.bkg.f009.nc
-# TODO(AFE) fix ice recentering in cycled da
-  #ice_filename: cice.res.nc
+  ice_filename: ice.bkg.f009.nc
   read_from_file: 1
 
 output increment:

ush/soca/marine_recenter.py

@@ -65,6 +65,7 @@ def __init__(self, config: Dict) -> None:
 
         _window_begin = add_to_datetime(self.task_config.current_cycle, -to_timedelta(f"{self.task_config.assim_freq}H") / 2)
         _window_end = add_to_datetime(self.task_config.current_cycle, to_timedelta(f"{self.task_config.assim_freq}H") / 2)
+        _enspert_relpath = os.path.relpath(self.task_config.DATAens, self.task_config.DATA)
 
         local_dict = AttrDict({'window_begin': f"{window_begin.strftime('%Y-%m-%dT%H:%M:%SZ')}",
                                'PARMsoca': os.path.join(self.task_config.PARMgfs, 'gdas', 'soca'),
@@ -80,6 +81,8 @@ def __init__(self, config: Dict) -> None:
                                'stage_dir': DATA,
                                'soca_input_fix_dir': self.task_config.SOCA_INPUT_FIX_DIR,
                                'NMEM_ENS': self.task_config.NMEM_ENS,
+                               'GDUMP_ENS': self.task_config.GDUMP_ENS,
+                               'ENSPERT_RELPATH': _enspert_relpath,
                                'MARINE_WINDOW_LENGTH': f"PT{config['assim_freq']}H",
                                'recen_yaml_template': os.path.join(berror_yaml_dir, 'soca_ensrecenter.yaml'),
                                'recen_yaml_file': os.path.join(DATA, 'soca_ensrecenter.yaml'),
@@ -124,6 +127,9 @@ def initialize(self):
         # stage backgrounds
         bkg_list = parse_j2yaml(self.task_config.MARINE_DET_STAGE_BKG_YAML_TMPL, self.task_config)
         FileHandler(bkg_list).sync()
+        # stage ensemble backgrounds for soca2cice
+        ens_bkg_list = parse_j2yaml(self.task_config.MARINE_ENSDA_STAGE_BKG_YAML_TMPL, self.task_config)
+        FileHandler(ens_bkg_list).sync()
 
 #        ################################################################################
 #        # Copy initial condition
@@ -155,6 +161,37 @@ def initialize(self):
                 ens_member_list.append([fname_in, fname_out])
 
         FileHandler({'copy': ens_member_list}).sync()
+        # stage ensemble ice restarts
+        # make a copy of the CICE6 restart
+        logger.info("---------------- Stage ensemble CICE restarts")
+        # set the restart date, dependent on the cycling type
+        if self.task_config.DOIAU:
+            # forecast initialized at the begining of the DA window
+            fcst_begin = self.task_config.MARINE_WINDOW_BEGIN_ISO
+            rst_date = self.task_config.MARINE_WINDOW_BEGIN.strftime('%Y%m%d.%H%M%S')
+        else:
+            # forecast initialized at the middle of the DA window
+            fcst_begin = self.task_config.MARINE_WINDOW_MIDDLE_ISO
+            rst_date = self.task_config.MARINE_WINDOW_MIDDLE.strftime('%Y%m%d.%H%M%S')
+        ens_cice_list = []
+        for mem in range(1, nmem_ens+1):
+            mem_dir = os.path.join(self.task_config.ROTDIR,
+                                   f'enkf{RUN}.{gPDYstr}',
+                                   f'{gcyc}',
+                                   f'mem{str(mem).zfill(3)}',
+                                   'model',
+                                   'ice',
+                                   'restart')
+            mem_dir_real = os.path.realpath(mem_dir)
+            f00 = f'{rst_date}.cice_model.res.nc'
+            fname_in = os.path.abspath(os.path.join(mem_dir_real, f00))
+            fname_out = os.path.realpath(os.path.join(self.task_config.ens_dir,
+                                         "cice_model.res."+str(mem)+".nc"))
+            ens_cice_list.append([fname_in, fname_out])
+            fname_out = os.path.realpath(os.path.join(self.task_config.ens_dir,
+                                         "cice_model.res.output."+str(mem)+".nc"))
+            ens_cice_list.append([fname_in, fname_out])
+        FileHandler({'copy': ens_cice_list}).sync()
 
         ################################################################################
         # generate the YAML file for recenterer
@@ -228,10 +265,6 @@ def finalize(self):
         mem_dir_list = []
         copy_list = []
 
-        # Skip the analysis insertion into the CICE restart for now
-        # TODO (G): Add this back in when we have hardened the soca to cice
-        #           change of variable
-
         # Copy the recentering increment files to the member COMROOT directories
         for mem in range(1, nmem_ens+1):
             mem_dir = os.path.join(self.task_config.ROTDIR,
@@ -242,9 +275,29 @@ def finalize(self):
                                    'ocean')
             mem_dir_real = os.path.realpath(mem_dir)
             mem_dir_list.append(mem_dir_real)
-
             copy_list.append([f'ocn.recenter.incr.{str(mem)}.nc',
                               os.path.join(mem_dir_real, incr_file)])
 
         FileHandler({'mkdir': mem_dir_list}).sync()
         FileHandler({'copy': copy_list}).sync()
+
+        # Copy the CICE restart files to the member COMROOT directories
+        if os.getenv('DOIAU') == "YES":
+            cice_rst_date = self.task_config.MARINE_WINDOW_BEGIN.strftime('%Y%m%d.%H%M%S')
+        else:
+            cice_rst_date = cdate.strftime('%Y%m%d.%H%M%S')
+        mem_dir_list = []
+        copy_list = []
+        for mem in range(1, nmem_ens+1):
+            mem_dir = os.path.join(self.task_config.ROTDIR,
+                                   f'enkf{RUN}.{PDYstr}',
+                                   f'{cyc}',
+                                   f'mem{str(mem).zfill(3)}',
+                                   'analysis',
+                                   'ice')
+            mem_dir_real = os.path.realpath(mem_dir)
+            mem_dir_list.append(mem_dir_real)
+            copy_list.append([f'ens/cice_model.res.output.{str(mem)}.nc',
+                              os.path.join(mem_dir_real, f'{cice_rst_date}.cice_model_anl.res.nc')])
+        FileHandler({'mkdir': mem_dir_list}).sync()
+        FileHandler({'copy': copy_list}).sync()

utils/soca/gdas_ens_handler.h

@@ -11,6 +11,7 @@
 #include "oops/base/PostProcessor.h"
 #include "oops/mpi/mpi.h"
 #include "oops/runs/Application.h"
+#include "oops/util/ConfigFunctions.h"
 #include "oops/util/DateTime.h"
 #include "oops/util/Duration.h"
 #include "oops/util/FieldSetOperations.h"
@@ -20,6 +21,7 @@
 #include "soca/Increment/Increment.h"
 #include "soca/LinearVariableChange/LinearVariableChange.h"
 #include "soca/State/State.h"
+#include "soca/VariableChange/VariableChange.h"
 
 #include "gdas_postprocincr.h"
 
@@ -147,7 +149,7 @@ namespace gdasapp {
 
       // Check if we're only re-centering the ensemble fcst around the det.
       bool recenterOnly = fullConfig.getBool("recentering around deterministic", false);
-
+      bool seaiceRecenter = fullConfig.getBool("sea ice recenter", false);
       // Save increments and exit if all we're doing is re-centering
       // the ensemble fcst around the det.
       if (recenterOnly) {
@@ -169,7 +171,32 @@ namespace gdasapp {
           postProcIncr.setToZero(incr);
 
           // Save the increments used to initialize the ensemble forecast
-          result = postProcIncr.save(mom6_incr, i+1, {"ocn"});
+          result = postProcIncr.save(mom6_incr, i+1);
+
+          // recenter ice if needed
+          if (seaiceRecenter) {
+            // read state
+            eckit::LocalConfiguration ensmem_config(fullConfig, "sea ice analysis");
+            std::string pattern;
+            fullConfig.get("sea ice analysis.pattern", pattern);
+            // replace templated string if necessary
+            if (!pattern.empty()) {
+              util::seekAndReplace(ensmem_config, pattern, std::to_string(i+1));
+            }
+            // assuming that this option is only used when recentering increment is on the
+            // same geometry
+            soca::State ens_an(geomOut, ensmem_config);
+            oops::Log::info() << "recentering ice state " << i << ":" << ens_an << std::endl;
+            ens_an += recenteringIncr;
+            oops::Log::info() << "recentered ice state " << i << ":" << ens_an << std::endl;
+            // set up variable change
+            eckit::LocalConfiguration varchange_config(fullConfig, "sea ice variable change");
+            util::seekAndReplace(varchange_config, pattern, std::to_string(i+1));
+            soca::VariableChange varchange(varchange_config, geomOut);
+            oops::Variables varout(varchange_config, "output variables");
+            // output happens inside soca2cice?
+            varchange.changeVar(ens_an, varout);
+          }
         }
         return result;
       }