From afca3d03748ff7d4940725da48b1a03bd84f0dd3 Mon Sep 17 00:00:00 2001 From: Anna Shlyaeva Date: Thu, 6 Feb 2025 12:59:28 +0000 Subject: [PATCH 1/2] Add sea ice recentering --- parm/soca/berror/soca_ensrecenter.yaml | 63 ++++++++++++++++++++++++-- ush/soca/marine_recenter.py | 63 ++++++++++++++++++++++++-- utils/soca/gdas_ens_handler.h | 31 ++++++++++++- 3 files changed, 146 insertions(+), 11 deletions(-) diff --git a/parm/soca/berror/soca_ensrecenter.yaml b/parm/soca/berror/soca_ensrecenter.yaml index b7d8e2032..ffe81f655 100644 --- a/parm/soca/berror/soca_ensrecenter.yaml +++ b/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: diff --git a/ush/soca/marine_recenter.py b/ush/soca/marine_recenter.py index 4c4a0c285..9f7ad38a2 100644 --- a/ush/soca/marine_recenter.py +++ b/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() diff --git a/utils/soca/gdas_ens_handler.h b/utils/soca/gdas_ens_handler.h index c60723db5..1f8d0b4b8 100644 --- a/utils/soca/gdas_ens_handler.h +++ b/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" @@ -19,6 +20,7 @@ #include "soca/Geometry/Geometry.h" #include "soca/Increment/Increment.h" #include "soca/LinearVariableChange/LinearVariableChange.h" +#include "soca/VariableChange/VariableChange.h" #include "soca/State/State.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; } From 18b2109aa66d722d6a39aa0e20a74ffed5efb93e Mon Sep 17 00:00:00 2001 From: Anna Shlyaeva Date: Thu, 6 Feb 2025 21:18:35 +0000 Subject: [PATCH 2/2] coding norms (!) --- ush/soca/marine_recenter.py | 6 +++--- utils/soca/gdas_ens_handler.h | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/ush/soca/marine_recenter.py b/ush/soca/marine_recenter.py index 9f7ad38a2..dd7ff3afd 100644 --- a/ush/soca/marine_recenter.py +++ b/ush/soca/marine_recenter.py @@ -186,10 +186,10 @@ def initialize(self): 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")) + "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")) + "cice_model.res.output."+str(mem)+".nc")) ens_cice_list.append([fname_in, fname_out]) FileHandler({'copy': ens_cice_list}).sync() @@ -298,6 +298,6 @@ def finalize(self): 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')]) + 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() diff --git a/utils/soca/gdas_ens_handler.h b/utils/soca/gdas_ens_handler.h index 1f8d0b4b8..d3f22d15f 100644 --- a/utils/soca/gdas_ens_handler.h +++ b/utils/soca/gdas_ens_handler.h @@ -20,8 +20,8 @@ #include "soca/Geometry/Geometry.h" #include "soca/Increment/Increment.h" #include "soca/LinearVariableChange/LinearVariableChange.h" -#include "soca/VariableChange/VariableChange.h" #include "soca/State/State.h" +#include "soca/VariableChange/VariableChange.h" #include "gdas_postprocincr.h"