streamlined and more diag for static B

CMakeLists.txt

@@ -71,7 +71,7 @@ if(BUILD_GDASBUNDLE)
   endif()
 
 # Core JEDI repositories
-  ecbuild_bundle( PROJECT oops  GIT "https://github.com/jcsda/oops.git"  BRANCH feature/fsdivide-skipabort)
+  ecbuild_bundle( PROJECT oops  GIT "https://github.com/jcsda/oops.git"  BRANCH feature/fsdivide-skipabort-noci)
   ecbuild_bundle( PROJECT vader GIT "https://github.com/jcsda/vader.git" BRANCH develop )
   ecbuild_bundle( PROJECT saber GIT "https://github.com/jcsda/saber.git" BRANCH develop )
   ecbuild_bundle( PROJECT ioda  GIT "https://github.com/jcsda/ioda.git" BRANCH develop  )

parm/soca/berror/saber_blocks.yaml

@@ -10,16 +10,6 @@ components:
       output variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh]
       linear variable changes:
 
-      - linear variable change name: HorizFiltSOCA
-        niter: 5
-        filter variables: [cicen, socn, tocn, ssh]
-
-      - linear variable change name: VertConvSOCA
-        Lz_min: 10.0
-        Lz_mld: 0
-        Lz_mld_max: 1.0
-        scale_layer_thick: 5
-
       - linear variable change name: BkgErrFILT
         ocean_depth_min: 500 # [m]
         rescale_bkgerr: 1.0
@@ -28,8 +18,8 @@ components:
       - linear variable change name: BkgErrSOCA
         read_from_file: 3
         basename: ./static_ens/
-        ocn_filename: 'ocn.orig_ens_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc'
-        ice_filename: 'ice.orig_ens_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc'
+        ocn_filename: 'ocn.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc'
+        ice_filename: 'ice.bkgerr_stddev.incr.{{ATM_WINDOW_BEGIN}}.nc'
         date: '{{ATM_WINDOW_MIDDLE}}'
         t_min: 0.1
         t_max: 10.0
@@ -44,8 +34,7 @@ components:
         standard deviation: true
 
       - linear variable change name: BalanceSOCA
-#        ksshts:
-#          nlayers: 10
+
   weight:
     value: 1.0
 - covariance:
@@ -55,8 +44,8 @@ components:
         read_from_file: 1
         date: '{{ATM_WINDOW_MIDDLE}}'
         basename: ./static_ens/
-        ocn_filename: 'ocn.pert.steric.%mem%.{{ATM_WINDOW_BEGIN}}.nc'
-        ice_filename: 'ice.pert.ens.%mem%.{{ATM_WINDOW_BEGIN}}.PT0S.nc'
+        ocn_filename: 'ocn.pert.steric.%mem%.nc'
+        ice_filename: 'ice.%mem%.nc'
         state variables: [tocn, socn, ssh, uocn, vocn, cicen, hicen, hsnon]
       pattern: '%mem%'
       nmembers: ${CLIM_ENS_SIZE}

parm/soca/berror/soca_ensb.yaml

@@ -0,0 +1,86 @@
+geometry:
+  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, cicen, hicen, hsnon]
+
+set increment variables to zero: [uocn, vocn, ssh]
+
+vertical geometry:
+  date: '{{ATM_WINDOW_BEGIN}}'
+  basename: ./INPUT/
+  ocn_filename: MOM.res.nc
+  read_from_file: 3
+
+soca increments:
+  number of increments: ${CLIM_ENS_SIZE}
+  pattern: '%mem%'
+  template:
+    date: '{{ATM_WINDOW_BEGIN}}'
+    basename: ./static_ens/
+    ocn_filename: 'ocn.%mem%.nc'
+    ice_filename: 'ice.%mem%.nc'
+    read_from_file: 3
+
+steric height:
+  linear variable changes:
+  - linear variable change name: BalanceSOCA  # Only the steric balance is applied
+
+ssh output:
+  unbalanced:
+    datadir: ./static_ens
+    date: '{{ATM_WINDOW_BEGIN}}'
+    exp: ssh_unbal_stddev
+    type: incr
+
+  steric:
+    datadir: ./static_ens
+    date: '{{ATM_WINDOW_BEGIN}}'
+    exp: ssh_steric_stddev
+    type: incr
+
+  total:
+    datadir: ./static_ens
+    date: '{{ATM_WINDOW_BEGIN}}'
+    exp: ssh_total_stddev
+    type: incr
+
+  explained variance:
+    datadir: ./static_ens
+    date: '{{ATM_WINDOW_BEGIN}}'
+    exp: steric_explained_variance
+    type: incr
+
+background error output:
+  datadir: ./static_ens
+  date: '{{ATM_WINDOW_BEGIN}}'
+  exp: bkgerr_stddev
+  type: incr
+
+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, cicen, hicen, hsnon]
+  date: '{{ATM_WINDOW_BEGIN}}'
+  basename: ./INPUT/
+  ocn_filename: MOM.res.nc
+  ice_filename: cice.res.nc
+  read_from_file: 1
+
+output increment:
+  datadir: ./static_ens
+  date: '{{ATM_WINDOW_BEGIN}}'
+  exp: trash
+  type: incr
+  output file: 'ocn.pert.steric.%mem%.nc'
+  pattern: '%mem%'

scripts/exgdas_global_marine_analysis_bmat.sh

@@ -123,42 +123,14 @@ fi
 # 3 - Copy h from deterministic to unbalanced perturbations
 # 4 - Compute moments of converted ensemble perturbations
 
-# Compute ensemble moments
+# Process climatological ensemble
 clean_yaml soca_clim_ens_moments.yaml
-$APRUN_OCNANAL $JEDI_BIN/soca_ensmeanandvariance.x soca_clim_ens_moments.yaml
+$APRUN_OCNANAL $JEDI_BIN/gdas_ens_handler.x soca_ensb.yaml
 export err=$?; err_chk
 if [ $err -gt 0  ]; then
     exit $err
 fi
 
-# Zero out std. dev of ssh to use the balance as a strong constraint
-# and append h
-# 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
-export err=$?; err_chk
-if [ $err -gt 0  ]; then
-    exit $err
-fi
-
-# Vertical filtering of the 3D perturbations and recompute the steric height
-# referenced to the deterministic background (f003)
-clean_yaml 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
-fi
-
-
 ################################################################################
 # Correlation and Localization operators
 shopt -s nullglob
@@ -226,14 +198,6 @@ if [ $err -gt 0  ]; then
     exit $err
 fi
 
-################################################################################
-# Create ensemble of perturbations for the cycle
-
-# Use ensemble recenter with "zerocenter"
-
-# Apply inverse balance to perturbations, use deterministic background as trajectory
-
-
 ################################################################################
 set +x
 if [ $VERBOSE = "YES" ]; then

scripts/exgdas_global_marine_analysis_post.py

@@ -60,11 +60,16 @@ def list_all_files(dir_in, dir_out, wc='*', fh_list=[]):
                        os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocninc.nc')])
 
 # Copy of the diagonal of the background error for the cycle
-post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.orig_ens_stddev.incr.{bdate}.nc'),
+post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.bkgerr_stddev.incr.{bdate}.nc'),
                        os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.bkgerr_stddev.nc')])
-post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ice.orig_ens_stddev.incr.{bdate}.nc'),
+post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ice.bkgerr_stddev.incr.{bdate}.nc'),
                        os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ice.bkgerr_stddev.nc')])
 
+# Copy of the ssh diagnostics
+for string in ['ssh_steric_stddev', 'ssh_unbal_stddev', 'ssh_total_stddev', 'steric_explained_variance']:
+    post_file_list.append([os.path.join(anl_dir, 'static_ens', f'ocn.{string}.incr.{bdate}.nc'),
+                           os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.{string}.nc')])
+
 # Copy the ice and ocean increments
 post_file_list.append([os.path.join(anl_dir, 'Data', f'ocn.3dvarfgat_pseudo.incr.{mdate}.nc'),
                        os.path.join(com_ocean_analysis, f'{RUN}.t{cyc}z.ocn.incr.nc')])

scripts/exgdas_global_marine_analysis_prep.py

@@ -334,56 +334,17 @@ def find_clim_ens(input_date):
 gridgen_yaml_dst = os.path.abspath(os.path.join(stage_cfg['stage_dir'], 'gridgen.yaml'))
 FileHandler({'copy': [[gridgen_yaml_src, gridgen_yaml_dst]]}).sync()
 
-
 ################################################################################
-# generate YAML file for parametric diag of B
-
-logging.info(f"---------------- generate parametric_stddev_b.yaml")
-berr_yaml = os.path.join(anl_dir, 'parametric_stddev_b.yaml')
-berr_yaml_template = os.path.join(gdas_home,
-                                  'parm',
-                                  'soca',
-                                  'berror',
-                                  'parametric_stddev_b.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)
-
-################################################################################
-# generate YAMLS file for diag of clim. ens. B
+# generate the YAML file for the post processing of the clim. ens. B
 berror_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror')
 
-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')
+logging.info(f"---------------- generate soca_ensb.yaml")
+berr_yaml = os.path.join(anl_dir, 'soca_ensb.yaml')
+berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_ensb.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_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")
-berr_yaml = os.path.join(anl_dir, 'soca_clim_ens_perts.yaml')
-berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_clim_ens_perts.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_apply_steric.yaml")
-berr_yaml = os.path.join(anl_dir, 'soca_apply_steric.yaml')
-berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_apply_steric.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_ensweights.yaml")
 berr_yaml = os.path.join(anl_dir, 'soca_ensweights.yaml')
 berr_yaml_template = os.path.join(berror_yaml_dir, 'soca_ensweights.yaml')

utils/soca/gdas_ens_handler.h

@@ -44,7 +44,7 @@ namespace gdasapp_ens_utils {
       }
       ensMean *= 1.0 / nm;
 
-      // Standard deviation
+      // Variance
       const double rk = 1.0/(nm - 1.0);
       ensVar.zero();
       for (const soca::Increment& incr : ensMembers) {
@@ -210,7 +210,7 @@ namespace gdasapp {
       oops::Log::info() << "filtered mean: " << ensMean << std::endl;
       oops::Log::info() << "filtered std: " << ensStd << std::endl;
       eckit::LocalConfiguration bkgErrOutputConfig(fullConfig, "background error output");
-      sshStd.write(nonStericSshOutputConfig);
+      ensStd.write(bkgErrOutputConfig);
 
       // Explained variance by steric height R=1-SS(non-steric ssh)/SS(total ssh)
       soca::Increment varianceRatio(geom, socaSshVar, postProcIncr.dt_);
@@ -230,45 +230,11 @@ namespace gdasapp {
 
       eckit::LocalConfiguration expvarSshOutputConfig(fullConfig, "ssh output.explained variance");
       stericExplainedVariance.write(expvarSshOutputConfig);
-
-
       oops::Log::info() << "explained variance: " << stericExplainedVariance << std::endl;
 
       return result;
     }
 
-    /*    // -----------------------------------------------------------------------------
-    // Compute moments
-    // TODO(G): use EnsembleIncrememt after the parameter class is removed?
-    void ensMoments(const std::vector<soca::Increment> ensMembers,
-                    soca::Increment& ensMean, soca::Increment& ensStd) const {
-      double nm = static_cast<double>(ensMembers.size());
-
-      // Mean
-      ensMean.zero();
-      for (const soca::Increment& incr : ensMembers) {
-        ensMean += incr;
-      }
-      ensMean *= 1.0 / nm;
-
-      // Standard deviation
-      const double rk = 1.0/(nm - 1.0);
-      ensStd.zero();
-      for (const soca::Increment& incr : ensMembers) {
-        oops::Log::info() << "members: " << incr << std::endl;
-        soca::Increment pert(incr);
-        pert -= ensMean;
-        pert.schur_product_with(pert);
-        ensStd += pert;
-      }
-      ensStd.axpy(rk, ensStd);
-      atlas::FieldSet ensStdFs;
-      ensStd.toFieldSet(ensStdFs);
-      util::sqrtFieldSet(ensStdFs);
-      ensStd.fromFieldSet(ensStdFs);
-    }
-    */
-    // -----------------------------------------------------------------------------
    private:
     util::DateTime dt_;