From 8ca585bf2f54e8f1b3b7a79e8c5c45919f36b117 Mon Sep 17 00:00:00 2001 From: AndrewEichmann-NOAA <58948505+AndrewEichmann-NOAA@users.noreply.github.com> Date: Mon, 7 Aug 2023 20:28:30 -0400 Subject: [PATCH 1/2] removed duplicate variable specification (#551) * removed duplicate variable specification * pep8 compliance --- scripts/exgdas_global_marine_analysis_vrfy.py | 55 +++++++++---------- ush/soca/soca_vrfy.py | 17 +++--- 2 files changed, 34 insertions(+), 38 deletions(-) diff --git a/scripts/exgdas_global_marine_analysis_vrfy.py b/scripts/exgdas_global_marine_analysis_vrfy.py index 88a33f49b..fa72a148b 100755 --- a/scripts/exgdas_global_marine_analysis_vrfy.py +++ b/scripts/exgdas_global_marine_analysis_vrfy.py @@ -48,11 +48,11 @@ config = plotConfig(grid_file=grid_file, data_file=data_file, lats=np.arange(-60, 60, 10), - variables_zonal=['Temp', 'Salt'], - variables_horiz=['Temp', 'Salt', 'ave_ssh'], - allbounds={'Temp': [-0.5, 0.5], - 'Salt': [-0.1, 0.1], - 'ave_ssh': [-0.1, 0.1]}, + variables_zonal={'Temp': [-0.5, 0.5], + 'Salt': [-0.1, 0.1]}, + variables_horiz={'Temp': [-0.5, 0.5], + 'Salt': [-0.1, 0.1], + 'ave_ssh': [-0.1, 0.1]}, colormap='RdBu', comout=os.path.join(comout, 'vrfy', 'incr')) ocnIncPlotter = statePlotter(config) @@ -66,10 +66,9 @@ config = plotConfig(grid_file=grid_file, data_file=data_file, lats=np.arange(-60, 60, 10), - variables_horiz=['aicen', 'hicen', 'hsnon'], - allbounds={'aicen': [-0.2, 0.2], - 'hicen': [-0.5, 0.5], - 'hsnon': [-0.1, 0.1]}, + variables_horiz={'aicen': [-0.2, 0.2], + 'hicen': [-0.5, 0.5], + 'hsnon': [-0.1, 0.1]}, colormap='RdBu', projs=['North', 'South'], comout=os.path.join(comout, 'vrfy', 'incr')) @@ -83,10 +82,9 @@ data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.iceana.nc') config = plotConfig(grid_file=grid_file, data_file=data_file, - variables_horiz=['aicen', 'hicen', 'hsnon'], - allbounds={'aicen': [0.0, 1.0], - 'hicen': [0.0, 4.0], - 'hsnon': [0.0, 0.5]}, + variables_horiz={'aicen': [0.0, 1.0], + 'hicen': [0.0, 4.0], + 'hsnon': [0.0, 0.5]}, colormap='jet', projs=['North', 'South', 'Global'], comout=os.path.join(comout, 'vrfy', 'ana')) @@ -100,10 +98,9 @@ data_file = os.path.join(com_ice_history, f'{RUN}.t{gcyc}z.icef006.nc') config = plotConfig(grid_file=grid_file, data_file=data_file, - variables_horiz=['aice_h', 'hs_h', 'hi_h'], - allbounds={'aice_h': [0.0, 1.0], - 'hs_h': [0.0, 4.0], - 'hi_h': [0.0, 0.5]}, + variables_horiz={'aice_h': [0.0, 1.0], + 'hs_h': [0.0, 4.0], + 'hi_h': [0.0, 0.5]}, colormap='jet', projs=['North', 'South', 'Global'], comout=os.path.join(comout, 'vrfy', 'bkg')) @@ -117,10 +114,9 @@ data_file = os.path.join(comout, f'{RUN}.t'+cyc+'z.ocnana.nc') config = plotConfig(grid_file=grid_file, data_file=data_file, - variables_horiz=['ave_ssh', 'Temp', 'Salt'], - allbounds={'ave_ssh': [-1.8, 1.3], - 'Temp': [-1.8, 34.0], - 'Salt': [30, 38]}, + variables_horiz={'ave_ssh': [-1.8, 1.3], + 'Temp': [-1.8, 34.0], + 'Salt': [30, 38]}, colormap='jet', comout=os.path.join(comout, 'vrfy', 'ana')) ocnAnaPlotter = statePlotter(config) @@ -133,10 +129,9 @@ data_file = os.path.join(com_ocean_history, f'{RUN}.t{gcyc}z.ocnf006.nc') config = plotConfig(grid_file=grid_file, data_file=data_file, - variables_horiz=['ave_ssh', 'Temp', 'Salt'], - allbounds={'ave_ssh': [-1.8, 1.3], - 'Temp': [-1.8, 34.0], - 'Salt': [30, 38]}, + variables_horiz={'ave_ssh': [-1.8, 1.3], + 'Temp': [-1.8, 34.0], + 'Salt': [30, 38]}, colormap='jet', comout=os.path.join(comout, 'vrfy', 'bkg')) ocnBkgPlotter = statePlotter(config) @@ -150,11 +145,11 @@ config = plotConfig(grid_file=grid_file, data_file=data_file, lats=np.arange(-60, 60, 10), - variables_zonal=['Temp', 'Salt'], - variables_horiz=['Temp', 'Salt', 'ave_ssh'], - allbounds={'Temp': [0, 2], - 'Salt': [0, 0.2], - 'ave_ssh': [0, 0.1]}, + variables_zonal={'Temp': [0, 2], + 'Salt': [0, 0.2]}, + variables_horiz={'Temp': [0, 2], + 'Salt': [0, 0.2], + 'ave_ssh': [0, 0.1]}, colormap='jet', comout=os.path.join(comout, 'vrfy', 'bkgerr')) bkgErrPlotter = statePlotter(config) diff --git a/ush/soca/soca_vrfy.py b/ush/soca/soca_vrfy.py index 248505711..2c8d602a0 100755 --- a/ush/soca/soca_vrfy.py +++ b/ush/soca/soca_vrfy.py @@ -19,14 +19,13 @@ def plotConfig(grid_file=[], data_file=[], variable=[], levels=[], - allbounds=[], bounds=[], colormap=[], max_depth=np.nan, max_depths=[700.0, 5000.0], comout=[], - variables_horiz=[], - variables_zonal=[], + variables_horiz={}, + variables_zonal={}, lat=np.nan, lats=np.arange(-60, 60, 10), proj='set me', @@ -41,7 +40,6 @@ def plotConfig(grid_file=[], config['fields file'] = data_file config['levels'] = [1] config['colormap'] = colormap - config['all bounds'] = allbounds config['bounds'] = bounds config['lats'] = lats # all the lats to plot config['lat'] = lat # the lat being currently plotted @@ -140,15 +138,18 @@ def plot(self): for max_depth in self.config['max depths']: self.config['max depth'] = max_depth - for variable in self.config['zonal variables']: - bounds = self.config['all bounds'][variable] + variableBounds = self.config['zonal variables'] + for variable in variableBounds.keys(): + bounds = variableBounds[variable] self.config.update({'variable': variable, 'bounds': bounds}) plotZonalSlice(self.config) ####################################### # Horizontal slices for proj in self.config['projs']: - for variable in self.config['horiz variables']: - bounds = self.config['all bounds'][variable] + + variableBounds = self.config['horiz variables'] + for variable in variableBounds.keys(): + bounds = variableBounds[variable] self.config.update({'variable': variable, 'bounds': bounds, 'proj': proj}) plotHorizontalSlice(self.config) From 7b671463d4714ae92799f27a6101fff9fbe5c87c Mon Sep 17 00:00:00 2001 From: Guillaume Vernieres Date: Mon, 7 Aug 2023 20:31:12 -0400 Subject: [PATCH 2/2] Processing of the offline ensemble (#549) * wip * wip * wip * post proc ensemble of incr * wip * fixed tests, apply new tool to B * code tidy * ... --- parm/soca/berror/saber_blocks.yaml | 4 +- parm/soca/berror/soca_apply_steric.yaml | 66 ++-- parm/soca/berror/soca_clim_ens_perts.yaml | 2 +- parm/soca/berror/soca_postproc_stddev.yaml | 31 ++ parm/soca/variational/3dvarfgat.yaml | 1 - parm/soca/variational/socaincr2mom6.yaml | 3 +- scripts/exgdas_global_marine_analysis_bmat.sh | 15 +- .../exgdas_global_marine_analysis_chkpt.sh | 2 +- scripts/exgdas_global_marine_analysis_prep.py | 20 +- test/soca/gw/CMakeLists.txt | 2 +- test/soca/testinput/incr_handler.yaml | 54 +++ test/soca/testinput/socaincr2mom6.yaml | 5 +- utils/CMakeLists.txt | 12 +- utils/gdas_incr_handler.cc | 8 + utils/gdas_incr_handler.h | 72 ++++ utils/gdas_postprocincr.h | 308 ++++++++++++++++++ utils/gdas_socaincr2mom6.cc | 8 - utils/gdas_socaincr2mom6.h | 112 ------- 18 files changed, 545 insertions(+), 180 deletions(-) create mode 100644 parm/soca/berror/soca_postproc_stddev.yaml create mode 100644 test/soca/testinput/incr_handler.yaml create mode 100644 utils/gdas_incr_handler.cc create mode 100644 utils/gdas_incr_handler.h create mode 100644 utils/gdas_postprocincr.h delete mode 100644 utils/gdas_socaincr2mom6.cc delete mode 100644 utils/gdas_socaincr2mom6.h diff --git a/parm/soca/berror/saber_blocks.yaml b/parm/soca/berror/saber_blocks.yaml index d039355da..53fc75567 100644 --- a/parm/soca/berror/saber_blocks.yaml +++ b/parm/soca/berror/saber_blocks.yaml @@ -56,8 +56,8 @@ components: read_from_file: 1 date: '{{ATM_WINDOW_MIDDLE}}' basename: ./static_ens/ - ocn_filename: 'ocn.filtered.%mem%.incr.{{ATM_WINDOW_BEGIN}}.nc' - ice_filename: 'ice.filtered.%mem%.incr.{{ATM_WINDOW_BEGIN}}.nc' + ocn_filename: 'ocn.pert.steric.%mem%.{{ATM_WINDOW_BEGIN}}.nc' + ice_filename: 'ice.pert.ens.%mem%.{{ATM_WINDOW_BEGIN}}.PT0S.nc' state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon] pattern: '%mem%' nmembers: ${CLIM_ENS_SIZE} diff --git a/parm/soca/berror/soca_apply_steric.yaml b/parm/soca/berror/soca_apply_steric.yaml index 62e9ded2b..ba8e894bd 100644 --- a/parm/soca/berror/soca_apply_steric.yaml +++ b/parm/soca/berror/soca_apply_steric.yaml @@ -1,43 +1,49 @@ -input geometry: +geometry: geom_grid_file: soca_gridspec.nc mom6_input_nml: mom_input.nml fields metadata: fields_metadata.yaml -output geometry: - geom_grid_file: soca_gridspec.nc - mom6_input_nml: mom_input.nml - fields metadata: fields_metadata.yaml +date: '{{ATM_WINDOW_BEGIN}}' + +layers variable: [hocn] + +increment variables: [tocn, socn, uocn, vocn, ssh, hocn] + +set increment variables to zero: [ssh] + +vertical geometry: + read_from_file: 1 + basename: ./INPUT/ + ocn_filename: MOM.res.nc + date: '{{ATM_WINDOW_BEGIN}}' + +soca increments: + number of increments: ${CLIM_ENS_SIZE} + pattern: '%mem%' + template: + date: '{{ATM_WINDOW_BEGIN}}' + basename: ./static_ens/ + ocn_filename: 'ocn.pert.ens.%mem%.{{ATM_WINDOW_BEGIN}}.PT0S.nc' + read_from_file: 1 linear variable change: - input variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon] - output variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon] - do inverse: false linear variable changes: - linear variable change name: BkgErrFILT ocean_depth_min: 500 # zero where ocean is shallower than 500m - rescale_bkgerr: 0.3 # rescale perturbation + rescale_bkgerr: 1.0 # rescale perturbation efold_z: 1500.0 # Apply exponential decay - - linear variable change name: BalanceSOCA # linear steric height from (T,S) perturbation - -increments: -- date: '{{ATM_WINDOW_BEGIN}}' - input variables: [tocn, socn, ssh, uocn, vocn, hocn, cicen, hicen, hsnon] - input: - read_from_file: 1 - basename: ./static_ens/ - ocn_filename: 'ocn.bal.ens.MEMNUM.{{ATM_WINDOW_BEGIN}}.PT0S.nc' - ice_filename: 'ice.bal.ens.MEMNUM.{{ATM_WINDOW_BEGIN}}.PT0S.nc' - date: '{{ATM_WINDOW_BEGIN}}' - state variables: [ssh, tocn, socn, uocn, vocn, cicen, hicen, hsnon] + - linear variable change name: BalanceSOCA trajectory: - read_from_file: 1 + state variables: [tocn, socn, uocn, vocn, ssh, hocn, layer_depth, mld] + date: '{{ATM_WINDOW_BEGIN}}' basename: ./INPUT/ ocn_filename: MOM.res.nc - ice_filename: cice.res.nc - date: '{{ATM_WINDOW_BEGIN}}' - state variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh, hocn, mld, layer_depth] - output: - datadir: ./static_ens - exp: filtered.MEMNUM - type: incr - date: '{{ATM_WINDOW_BEGIN}}' + read_from_file: 1 + +output increment: + datadir: ./static_ens/ + date: '{{ATM_WINDOW_BEGIN}}' + exp: tmp + type: incr + output file: 'ocn.pert.steric.%mem%.{{ATM_WINDOW_BEGIN}}.nc' + pattern: '%mem%' diff --git a/parm/soca/berror/soca_clim_ens_perts.yaml b/parm/soca/berror/soca_clim_ens_perts.yaml index 90d0abb36..b51decdb2 100644 --- a/parm/soca/berror/soca_clim_ens_perts.yaml +++ b/parm/soca/berror/soca_clim_ens_perts.yaml @@ -30,6 +30,6 @@ ensemble: recentered output: datadir: ./static_ens - exp: bal + exp: pert type: ens date: '{{ATM_WINDOW_BEGIN}}' diff --git a/parm/soca/berror/soca_postproc_stddev.yaml b/parm/soca/berror/soca_postproc_stddev.yaml new file mode 100644 index 000000000..6720e1b5a --- /dev/null +++ b/parm/soca/berror/soca_postproc_stddev.yaml @@ -0,0 +1,31 @@ +geometry: + geom_grid_file: soca_gridspec.nc + mom6_input_nml: mom_input.nml + fields metadata: fields_metadata.yaml + +date: '{{ATM_WINDOW_BEGIN}}' + +layers variable: [hocn] + +increment variables: [tocn, socn, uocn, vocn, ssh, hocn] + +set increment variables to zero: [ssh] + +vertical geometry: + read_from_file: 1 + basename: ./INPUT/ + ocn_filename: MOM.res.nc + date: '{{ATM_WINDOW_BEGIN}}' + +soca increment: + date: '{{ATM_WINDOW_BEGIN}}' + basename: ./static_ens/ + ocn_filename: 'ocn.orig_ens_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' + read_from_file: 1 + +output increment: + datadir: ./ + date: '{{ATM_WINDOW_BEGIN}}' + exp: filtered + type: incr + output file: 'ocn.orig_ens_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc' diff --git a/parm/soca/variational/3dvarfgat.yaml b/parm/soca/variational/3dvarfgat.yaml index de59a03d8..2e4371586 100644 --- a/parm/soca/variational/3dvarfgat.yaml +++ b/parm/soca/variational/3dvarfgat.yaml @@ -34,7 +34,6 @@ variational: fields metadata: ./fields_metadata.yaml ninner: !ENV ${SOCA_NINNER} gradient norm reduction: 1e-10 - test: on diagnostics: departures: ombg diff --git a/parm/soca/variational/socaincr2mom6.yaml b/parm/soca/variational/socaincr2mom6.yaml index 75101244e..bbada8283 100644 --- a/parm/soca/variational/socaincr2mom6.yaml +++ b/parm/soca/variational/socaincr2mom6.yaml @@ -20,8 +20,9 @@ soca increment: ocn_filename: 'ocn.3dvarfgat_pseudo.incr.{{ATM_WINDOW_MIDDLE}}.nc' read_from_file: 1 -mom6 iau increment: +output increment: datadir: ./ date: '{{ATM_WINDOW_BEGIN}}' exp: mom6_iau type: incr + output file: inc.nc diff --git a/scripts/exgdas_global_marine_analysis_bmat.sh b/scripts/exgdas_global_marine_analysis_bmat.sh index 28dab42cc..107a610e3 100755 --- a/scripts/exgdas_global_marine_analysis_bmat.sh +++ b/scripts/exgdas_global_marine_analysis_bmat.sh @@ -114,11 +114,11 @@ fi ################################################################################ # Compute the ens std. dev, ignore ssh variance # TODO (G): Implement what's below into one single oops application -# 0 - Compute moments of original ensemble, used at the diag of the first +# 0 - Compute moments of original ensemble, used at the diag of the first # component of the hybrid B # 1 - Ensemble perturbations: # o Vertically Interpolate to the deterministic layers -# o Recenter around 0 to create an ensemble of perturbations +# o Recenter around 0 to create an ensemble of perurbations # 2 - Filter members and apply the linear steric height balance to each members # 3 - Copy h from deterministic to unbalanced perturbations # 4 - Compute moments of converted ensemble perturbations @@ -131,6 +131,15 @@ if [ $err -gt 0 ]; then exit $err fi +# Zero out std. dev of ssh to use the balance as a strong constraint +# TODO: Apply the inverse of the balance +clean_yaml soca_clim_ens_moments.yaml +$APRUN_OCNANAL $JEDI_BIN/gdas_incr_handler.x soca_postproc_stddev.yaml +export err=$?; err_chk +if [ $err -gt 0 ]; then + exit $err +fi + # Compute ensemble perturbations, vertically remap to cycle's vertical geometry clean_yaml soca_clim_ens_perts.yaml $APRUN_OCNANAL $JEDI_BIN/soca_ensrecenter.x soca_clim_ens_perts.yaml @@ -141,7 +150,7 @@ fi # Vertical filtering of the 3D perturbations and recompute the steric height perturbation clean_yaml soca_apply_steric.yaml -$APRUN_OCNANAL $JEDI_BIN/soca_convertincrement.x soca_apply_steric.yaml +$APRUN_OCNANAL $JEDI_BIN/gdas_incr_handler.x soca_apply_steric.yaml export err=$?; err_chk if [ $err -gt 0 ]; then exit $err diff --git a/scripts/exgdas_global_marine_analysis_chkpt.sh b/scripts/exgdas_global_marine_analysis_chkpt.sh index be506379d..8b585693d 100755 --- a/scripts/exgdas_global_marine_analysis_chkpt.sh +++ b/scripts/exgdas_global_marine_analysis_chkpt.sh @@ -55,7 +55,7 @@ EOF --out "${mom6_iau_incr}" \ --nsst_yaml "nsst.yaml" else - $APRUN_OCNANAL ${JEDI_BIN}/gdas_socaincr2mom6.x socaincr2mom6.yaml + $APRUN_OCNANAL ${JEDI_BIN}/gdas_incr_handler.x socaincr2mom6.yaml fi export err=$? if [ $err -gt 0 ]; then diff --git a/scripts/exgdas_global_marine_analysis_prep.py b/scripts/exgdas_global_marine_analysis_prep.py index a6bef66bc..73f67b4e9 100755 --- a/scripts/exgdas_global_marine_analysis_prep.py +++ b/scripts/exgdas_global_marine_analysis_prep.py @@ -354,12 +354,18 @@ def find_clim_ens(input_date): # generate YAMLS file for diag of clim. ens. B berror_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror') -logging.info(f"---------------- generate soca_clim_moments.yaml") +logging.info(f"---------------- generate soca_clim_ens_moments.yaml") berr_yaml = os.path.join(anl_dir, 'soca_clim_ens_moments.yaml') berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_clim_ens_moments.yaml') config = YAMLFile(path=berr_yaml_template) config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) -config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) +config.save(berr_yaml) + +logging.info(f"---------------- generate soca_postproc_stddev.yaml") +berr_yaml = os.path.join(anl_dir, 'soca_postproc_stddev.yaml') +berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_postproc_stddev.yaml') +config = YAMLFile(path=berr_yaml_template) +config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) config.save(berr_yaml) logging.info(f"---------------- generate soca_clim_ens_perts.yaml") @@ -376,16 +382,6 @@ def find_clim_ens(input_date): config = YAMLFile(path=berr_yaml_template) config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get) config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get) -member_template = config['increments'][0] -ocn_fname_tmpl = member_template['input']['ocn_filename'] -ice_fname_tmpl = member_template['input']['ice_filename'] -config['increments'] = [] -for n in range(1, clim_ens_size+1): - member = copy.deepcopy(member_template) - member['input']['ocn_filename'] = member['input']['ocn_filename'].replace("MEMNUM", str(n)) - member['input']['ice_filename'] = member['input']['ice_filename'].replace("MEMNUM", str(n)) - member['output']['exp'] = member['output']['exp'].replace("MEMNUM", str(n)) - config['increments'].append(member) config.save(berr_yaml) logging.info(f"---------------- generate soca_ensweights.yaml") diff --git a/test/soca/gw/CMakeLists.txt b/test/soca/gw/CMakeLists.txt index e1b20f602..acb7d3d54 100644 --- a/test/soca/gw/CMakeLists.txt +++ b/test/soca/gw/CMakeLists.txt @@ -78,7 +78,7 @@ foreach(jjob ${jjob_list}) endforeach() # Test gdas/oops applications -set(ctest_list "socahybridweights" "socaincr2mom6") +set(ctest_list "socahybridweights" "incr_handler") foreach(ctest ${ctest_list}) set(TEST ${ctest}) set(EXEC ${PROJECT_BINARY_DIR}/bin/gdas_${ctest}.x) diff --git a/test/soca/testinput/incr_handler.yaml b/test/soca/testinput/incr_handler.yaml new file mode 100644 index 000000000..fb81932f4 --- /dev/null +++ b/test/soca/testinput/incr_handler.yaml @@ -0,0 +1,54 @@ +geometry: + mom6_input_nml: mom_input.nml + fields metadata: ./fields_metadata.yaml + +date: 2018-04-15T09:00:00Z + +layers variable: [hocn] + +increment variables: [tocn, socn, uocn, vocn, ssh, hocn] + +set increment variables to zero: [uocn, vocn, ssh] + +vertical geometry: + date: 2018-04-15T09:00:00Z + basename: ./INPUT/ + ocn_filename: MOM.res.nc + read_from_file: 1 + +soca increment: + date: 2018-04-15T09:00:00Z + basename: ./Data/ + ocn_filename: 'ocn.3dvarfgat_pseudo.incr.2018-04-15T12:00:00Z.nc' + read_from_file: 1 + +#TODO: add one more ctest to check the snippet below +#soca increments: +# number of increments: 1 +# pattern: incr +# template: +# date: 2018-04-15T09:00:00Z +# basename: ./static_ens/ +# ocn_filename: 'ocn.1.nc' +# read_from_file: 1 + +linear variable change: + linear variable changes: + - linear variable change name: BkgErrFILT + ocean_depth_min: 500 # zero where ocean is shallower than 500m + rescale_bkgerr: 1.0 # rescale perturbation + efold_z: 1500.0 # Apply exponential decay + - linear variable change name: BalanceSOCA + trajectory: + state variables: [tocn, socn, uocn, vocn, ssh, hocn, layer_depth, mld] + date: 2018-04-15T09:00:00Z + basename: ./INPUT/ + ocn_filename: MOM.res.nc + read_from_file: 1 + +output increment: + datadir: ./ + date: 2018-04-15T09:00:00Z + exp: mom6_iau + type: incr + output file: inc.nc #inc.%mem%.nc diff --git a/test/soca/testinput/socaincr2mom6.yaml b/test/soca/testinput/socaincr2mom6.yaml index a5b5cf600..3a78c8164 100644 --- a/test/soca/testinput/socaincr2mom6.yaml +++ b/test/soca/testinput/socaincr2mom6.yaml @@ -6,7 +6,7 @@ date: 2018-04-15T09:00:00Z layers variable: [hocn] -increment variables: [tocn, socn, uocn, vocn, ssh] +increment variables: [tocn, socn, uocn, vocn, ssh, hocn] vertical geometry: date: 2018-04-15T09:00:00Z @@ -20,8 +20,9 @@ soca increment: ocn_filename: 'ocn.3dvarfgat_pseudo.incr.2018-04-15T12:00:00Z.nc' read_from_file: 1 -mom6 iau increment: +output increment: datadir: ./ date: 2018-04-15T09:00:00Z exp: mom6_iau type: incr + output file: inc.nc diff --git a/utils/CMakeLists.txt b/utils/CMakeLists.txt index fa229981b..86be649b2 100644 --- a/utils/CMakeLists.txt +++ b/utils/CMakeLists.txt @@ -5,14 +5,14 @@ find_package(oops REQUIRED) find_package(atlas REQUIRED) find_package(soca REQUIRED) -ecbuild_add_executable( TARGET gdas_socaincr2mom6.x - SOURCES gdas_socaincr2mom6.cc ) - -target_compile_features( gdas_socaincr2mom6.x PUBLIC cxx_std_17) -target_link_libraries( gdas_socaincr2mom6.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) +# Increment post processing +ecbuild_add_executable( TARGET gdas_incr_handler.x + SOURCES gdas_incr_handler.cc gdas_postprocincr.h) +target_compile_features( gdas_incr_handler.x PUBLIC cxx_std_17) +target_link_libraries( gdas_incr_handler.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) +# Hybrid-Weight ecbuild_add_executable( TARGET gdas_socahybridweights.x SOURCES gdas_socahybridweights.cc ) - target_compile_features( gdas_socahybridweights.x PUBLIC cxx_std_17) target_link_libraries( gdas_socahybridweights.x PUBLIC NetCDF::NetCDF_CXX oops atlas soca) diff --git a/utils/gdas_incr_handler.cc b/utils/gdas_incr_handler.cc new file mode 100644 index 000000000..65ce9ac43 --- /dev/null +++ b/utils/gdas_incr_handler.cc @@ -0,0 +1,8 @@ +#include "gdas_incr_handler.h" +#include "oops/runs/Run.h" + +int main(int argc, char ** argv) { + oops::Run run(argc, argv); + gdasapp::SocaIncrHandler incrhandler; + return run.execute(incrhandler); +} diff --git a/utils/gdas_incr_handler.h b/utils/gdas_incr_handler.h new file mode 100644 index 000000000..1d3d7cf19 --- /dev/null +++ b/utils/gdas_incr_handler.h @@ -0,0 +1,72 @@ +#include +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.h" + +#include "oops/base/PostProcessor.h" +#include "oops/mpi/mpi.h" +#include "oops/runs/Application.h" +#include "oops/util/DateTime.h" +#include "oops/util/Duration.h" +#include "oops/util/Logger.h" + +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" +#include "soca/LinearVariableChange/LinearVariableChange.h" +#include "soca/State/State.h" + +# include "gdas_postprocincr.h" + +namespace gdasapp { + + class SocaIncrHandler : public oops::Application { + public: + explicit SocaIncrHandler(const eckit::mpi::Comm & comm = oops::mpi::world()) + : Application(comm) {} + static const std::string classname() {return "gdasapp::SocaIncrHandler";} + + int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { + + /// Setup the soca geometry + const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); + oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl; + const soca::Geometry geom(geomConfig, this->getComm()); + + + // Initialize the post processing + PostProcIncr postProcIncr(fullConfig, geom, this->getComm()); + + oops::Log::info() << "soca increments: " << std::endl << postProcIncr.inputIncrConfig_ << std::endl; + + // Process list of increments + int result = 0; + for (size_t i = 1; i < postProcIncr.ensSize_+1; ++i) { + oops::Log::info() << postProcIncr.inputIncrConfig_ << std::endl; + + // At the very minimum, we run this script to append the layers state, so do that! + soca::Increment incr = postProcIncr.appendLayer(i); + + // Zero out specified fields + incr = postProcIncr.setToZero(incr); + + // Apply linear change of variables + incr = postProcIncr.applyLinVarChange(incr); + + // Save final increment + result = postProcIncr.save(incr, i); + } + return result; + } + // ----------------------------------------------------------------------------- + private: + util::DateTime dt_; + + // ----------------------------------------------------------------------------- + std::string appname() const { + return "gdasapp::SocaIncrHandler"; + } + }; +} // namespace gdasapp diff --git a/utils/gdas_postprocincr.h b/utils/gdas_postprocincr.h new file mode 100644 index 000000000..613c9cd3c --- /dev/null +++ b/utils/gdas_postprocincr.h @@ -0,0 +1,308 @@ +#ifndef GDAS_POSTPROCINCR_H +#define GDAS_POSTPROCINCR_H + +#include +#include + +#include "eckit/config/LocalConfiguration.h" + +#include "atlas/field.h" + +#include "oops/base/PostProcessor.h" +#include "oops/mpi/mpi.h" +#include "oops/util/ConfigFunctions.h" +#include "oops/util/DateTime.h" +#include "oops/util/Logger.h" + +#include "soca/Geometry/Geometry.h" +#include "soca/Increment/Increment.h" +#include "soca/LinearVariableChange/LinearVariableChange.h" +#include "soca/State/State.h" + +namespace gdasapp { + +class PostProcIncr { +public: + // Constructor + PostProcIncr(const eckit::Configuration & fullConfig, const soca::Geometry& geom, + const eckit::mpi::Comm & comm) + : dt_(getDate(fullConfig)), + layerVar_(getLayerVar(fullConfig)), + geom_(geom), + Layers_(getLayerThickness(fullConfig, geom)), + xTraj_(getTraj(fullConfig, geom)), + comm_(comm), + ensSize_(1), + pattern_(){ + + oops::Log::info() << "Date: " << std::endl << dt_ << std::endl; + + // Increment variables + oops::Variables socaIncrVar(fullConfig, "increment variables"); + ASSERT(socaIncrVar.size() >= 1); + socaIncrVar_ = socaIncrVar; + + // Input increments configuration + if ( fullConfig.has("soca increments.template") ) { + fullConfig.get("soca increments.template", inputIncrConfig_); + fullConfig.get("soca increments.number of increments", ensSize_); + fullConfig.get("soca increments.pattern", pattern_); + } else { + fullConfig.get("soca increment", inputIncrConfig_); + } + + // Output incrememnt configuration + eckit::LocalConfiguration outputIncrConfig(fullConfig, "output increment"); + outputIncrConfig_ = outputIncrConfig; + + // Variables that should be set to 0 + setToZero_ = false; + if ( fullConfig.has("set increment variables to zero") ) { + oops::Variables socaZeroIncrVar(fullConfig, "set increment variables to zero"); + socaZeroIncrVar_ = socaZeroIncrVar; + setToZero_ = true; + } + + // Linear variable changes + doLVC_ = false; + if ( fullConfig.has("linear variable change") ) { + eckit::LocalConfiguration lvcConfig(fullConfig, "linear variable change"); + lvcConfig_ = lvcConfig; + doLVC_ = true; + } + } + + // Append layer thicknesses to increment + // TODO: There's got to be a better way to append a variable. + soca::Increment appendLayer(const int n){ + oops::Log::info() << "==========================================" << std::endl; + oops::Log::info() << "====== Append Layers" << std::endl; + + // initialize the soca increment + soca::Increment socaIncr(geom_, socaIncrVar_, dt_); + eckit::LocalConfiguration memberConfig; //inputIncrConfig_); + memberConfig = inputIncrConfig_; + + // replace templated string if necessary + if (not pattern_.empty()) { + util::seekAndReplace(memberConfig, pattern_, std::to_string(n)); + oops::Log::info() << "oooooooooooooooooooooooooooooooooooo" << memberConfig << std::endl; + } + + // read the soca increment + socaIncr.read(memberConfig); + oops::Log::info() << "-------------------- input increment: " << std::endl; + oops::Log::info() << socaIncr << std::endl; + + // concatenate variables + oops::Variables outputIncrVar(socaIncrVar_); + outputIncrVar += layerVar_; + oops::Log::info() << "-------------------- outputIncrVar: " << std::endl; + oops::Log::info() << outputIncrVar << std::endl; + + // append layer variable to the soca increment + atlas::FieldSet socaIncrFs; + socaIncr.toFieldSet(socaIncrFs); + socaIncr.updateFields(outputIncrVar); + + // pad layer increment with zeros + soca::Increment layerThick(Layers_); + atlas::FieldSet socaLayerThickFs; + oops::Log::info() << "-------------------- thickness field: " << std::endl; + oops::Log::info() << layerThick << std::endl; + layerThick.toFieldSet(socaLayerThickFs); + layerThick.updateFields(outputIncrVar); + + // append layer thinckness to increment + socaIncr += layerThick; + oops::Log::info() << "-------------------- output increment: " << std::endl; + oops::Log::info() << socaIncr << std::endl; + + return socaIncr; + } + + // Set specified variables to 0 + soca::Increment setToZero(soca::Increment socaIncr) { + oops::Log::info() << "==========================================" << std::endl; + if (not this->setToZero_) { + oops::Log::info() << "====== no variables to set to 0.0" << std::endl; + return socaIncr; + } + 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; + auto view = atlas::array::make_view(field); + view.assign(0.0); + } + } + socaIncr.fromFieldSet(socaIncrFs); + oops::Log::info() << "-------------------- increment with zero'ed out fields: " << std::endl; + oops::Log::info() << socaIncr << std::endl; + return socaIncr; + } + + // Apply linear variable changes + soca::Increment applyLinVarChange(soca::Increment socaIncr) { + oops::Log::info() << "==========================================" << std::endl; + if (not this->doLVC_) { + return socaIncr; + } + oops::Log::info() << "====== applying specified change of variables" << std::endl; + soca::LinearVariableChangeParameters params; + params.deserialize(lvcConfig_); + oops::Log::info() << params << std::endl; + soca::LinearVariableChange lvc(this->geom_, params); + oops::Log::info() << "traj: " << xTraj_ << std::endl; + lvc.changeVarTraj(xTraj_, socaIncrVar_); + lvc.changeVarTL(socaIncr, socaIncrVar_); + oops::Log::info() << "incr after var change: " << std::endl << socaIncr << std::endl; + + return socaIncr; + } + + // Save increment + int save(soca::Increment socaIncr, int ensMem = 1) { + oops::Log::info() << "==========================================" << std::endl; + oops::Log::info() << "-------------------- save increment: " << std::endl; + socaIncr.write(outputIncrConfig_); + + // wait for everybody to be done + comm_.barrier(); + + // Change soca standard output name to something specified in the config + int result = 0; + if ( comm_.rank() == 0 ) { + // get the output directory + std::string dataDir; + outputIncrConfig_.get("datadir", dataDir); + // get the output file name + std::string outputFileName; + outputIncrConfig_.get("output file", outputFileName); + outputFileName = dataDir + "/" + outputFileName; + if (outputIncrConfig_.has("pattern")) { + std::string pattern; + outputIncrConfig_.get("pattern", pattern); + outputFileName = this->swapPattern(outputFileName, pattern, std::to_string(ensMem)); + oops::Log::info() << "-------------------- pattern: " << pattern << std::endl; + } + const char* charPtrOut = outputFileName.c_str(); + + std::string incrFname = this->socaFname(); + const char* charPtr = incrFname.c_str(); + + oops::Log::info() << "rename: " << incrFname << " to " << outputFileName << std::endl; + result = std::rename(charPtr, charPtrOut); + } + return result; + } + + // ----------------------------------------------------------------------------- + + // Initializers + // ----------------------------------------------------------------------------- + // Date from config + util::DateTime getDate(const eckit::Configuration& fullConfig) const { + std::string strdt; + fullConfig.get("date", strdt); + return util::DateTime(strdt); + } + // get the layer variable + oops::Variables getLayerVar(const eckit::Configuration& fullConfig) const { + oops::Variables layerVar(fullConfig, "layers variable"); + ASSERT(layerVar.size() == 1); + return layerVar; + } + // Read the layer thickness from the relevant background + soca::Increment getLayerThickness(const eckit::Configuration& fullConfig, + const soca::Geometry& geom) const { + soca::Increment layerThick(geom, getLayerVar(fullConfig), getDate(fullConfig)); + const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry"); + layerThick.read(vertGeomConfig); + oops::Log::info() << "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, "linear variable change.trajectory"); + soca::State traj(geom, trajConfig); + oops::Log::info() << "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; + } + } + + // ----------------------------------------------------------------------------- + + // Utility functions + // ----------------------------------------------------------------------------- + // Recreate the soca filename from the configuration + // TODO: Change this in soca? + // TODO: Hard-coded for ocean, implement for seaice as well + std::string socaFname() { + std::string datadir; + outputIncrConfig_.get("datadir", datadir); + std::filesystem::path pathToResolve = datadir; + std::string exp; + outputIncrConfig_.get("exp", exp); + std::string outputType; + outputIncrConfig_.get("type", outputType); + std::string incrFname = std::filesystem::canonical(pathToResolve); + incrFname += "/ocn." + exp + "." + outputType + "." + dt_.toString() + ".nc"; + + return incrFname; + } + + // Function to replace all occurrences of a pattern in a string with a replacement + std::string swapPattern(const std::string& input, + const std::string& pattern, + const std::string& replacement) { + std::string result = input; + size_t startPos = 0; + + while ((startPos = result.find(pattern, startPos)) != std::string::npos) { + result.replace(startPos, pattern.length(), replacement); + startPos += replacement.length(); + } + + return result; +} + + +public: + util::DateTime dt_; // valid date of increment + oops::Variables layerVar_; // layer variable + const soca::Increment Layers_; // layer thicknesses + const soca::Geometry & geom_; + const eckit::mpi::Comm & comm_; + // std::vector inputIncrConfig_; + eckit::LocalConfiguration inputIncrConfig_; + eckit::LocalConfiguration outputIncrConfig_; + eckit::LocalConfiguration zeroIncrConfig_; + eckit::LocalConfiguration lvcConfig_; + oops::Variables socaIncrVar_; + bool setToZero_; + bool doLVC_; + const soca::State xTraj_; + oops::Variables socaZeroIncrVar_; + int ensSize_; + std::string pattern_; +}; +} // namespace gdasapp + +#endif // GDAS_POSTPROCINCR_H diff --git a/utils/gdas_socaincr2mom6.cc b/utils/gdas_socaincr2mom6.cc deleted file mode 100644 index 0b1230f56..000000000 --- a/utils/gdas_socaincr2mom6.cc +++ /dev/null @@ -1,8 +0,0 @@ -#include "gdas_socaincr2mom6.h" -#include "oops/runs/Run.h" - -int main(int argc, char ** argv) { - oops::Run run(argc, argv); - gdasapp::SocaIncr2Mom6 socaincr2mom6; - return run.execute(socaincr2mom6); -} diff --git a/utils/gdas_socaincr2mom6.h b/utils/gdas_socaincr2mom6.h deleted file mode 100644 index c644b4c01..000000000 --- a/utils/gdas_socaincr2mom6.h +++ /dev/null @@ -1,112 +0,0 @@ -#include -#include -#include - -#include "eckit/config/LocalConfiguration.h" - -#include "atlas/field.h" - -#include "oops/base/PostProcessor.h" -#include "oops/mpi/mpi.h" -#include "oops/runs/Application.h" -#include "oops/util/DateTime.h" -#include "oops/util/Duration.h" -#include "oops/util/Logger.h" - -#include "soca/Geometry/Geometry.h" -#include "soca/Increment/Increment.h" - -namespace gdasapp { - - class SocaIncr2Mom6 : public oops::Application { - public: - explicit SocaIncr2Mom6(const eckit::mpi::Comm & comm = oops::mpi::world()) - : Application(comm) {} - static const std::string classname() {return "gdasapp::SocaIncr2Mom6";} - - int execute(const eckit::Configuration & fullConfig, bool /*validate*/) const { - - /// Setup the soca geometry - const eckit::LocalConfiguration geomConfig(fullConfig, "geometry"); - oops::Log::info() << "geometry: " << std::endl << geomConfig << std::endl; - const soca::Geometry geom(geomConfig, this->getComm()); - - /// Setup the vertical geometry from the background (layer thicknesses) - // get date - std::string strdt; - fullConfig.get("date", strdt); - util::DateTime dt = util::DateTime(strdt); - - // get layer thickness variable name - oops::Variables layerVar(fullConfig, "layers variable"); - ASSERT(layerVar.size() == 1); - - // read layer thicknesses from the relevant background - soca::Increment layerThick(geom, layerVar, dt); - const eckit::LocalConfiguration vertGeomConfig(fullConfig, "vertical geometry"); - layerThick.read(vertGeomConfig); - oops::Log::info() << "layerThick: " << std::endl << layerThick << std::endl; - - // Setup the soca increment - // get the increment variables - oops::Variables socaIncrVar(fullConfig, "increment variables"); - ASSERT(socaIncrVar.size() >= 1); - - // read the soca increment - soca::Increment socaIncr(geom, socaIncrVar, dt); - const eckit::LocalConfiguration socaIncrConfig(fullConfig, "soca increment"); - socaIncr.read(socaIncrConfig); - oops::Log::info() << "socaIncr: " << std::endl << socaIncr << std::endl; - - /// Create the MOM6 IAU increment - // concatenate variables - oops::Variables mom6IauVar(socaIncrVar); - mom6IauVar += layerVar; - oops::Log::info() << "mom6IauVar: " << std::endl << mom6IauVar << std::endl; - - // append layer variable to soca increment - atlas::FieldSet socaIncrFs; - socaIncr.toFieldSet(socaIncrFs); - socaIncr.updateFields(mom6IauVar); - oops::Log::info() << "MOM6 incr: " << std::endl << socaIncr << std::endl; - - // pad layer increment with zeros - atlas::FieldSet socaLayerThickFs; - layerThick.toFieldSet(socaLayerThickFs); - layerThick.updateFields(mom6IauVar); - oops::Log::info() << "h: " << std::endl << layerThick << std::endl; - - // append layer thinckness to increment - socaIncr += layerThick; - oops::Log::info() << "MOM6 IAU increment: " << std::endl << socaIncr << std::endl; - - // Save MOM6 IAU Increment - const eckit::LocalConfiguration mom6IauConfig(fullConfig, "mom6 iau increment"); - socaIncr.write(mom6IauConfig); - - // TODO: the "checkpoint" script expects the ocean increment output to - // be in "inc.nc". Remove what's below, eventually - std::string datadir; - mom6IauConfig.get("datadir", datadir); - std::filesystem::path pathToResolve = datadir; - std::string exp; - mom6IauConfig.get("exp", exp); - std::string outputType; - mom6IauConfig.get("type", outputType); - std::string incrFname = std::filesystem::canonical(pathToResolve); - incrFname += "/ocn." + exp + "." + outputType + "." + dt.toString() + ".nc"; - const char* charPtr = incrFname.c_str(); - oops::Log::info() << "rename: " << incrFname << " to " << "inc.nc" << std::endl; - int result = std::rename(charPtr, "inc.nc"); - - return result; - } - // ----------------------------------------------------------------------------- - private: - std::string appname() const { - return "gdasapp::SocaIncr2Mom6"; - } - // ----------------------------------------------------------------------------- - }; - -} // namespace gdasapp