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

@@ -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}

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%'

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}}'

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'

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
 

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

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

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

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")

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)

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

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

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)

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);
+}