using diffusion for correlation

parm/soca/berror/saber_blocks.yaml

@@ -1,25 +1,35 @@
 covariance model: hybrid
 components:
 - covariance:
+    # This will setup B = K D C_v C_h C_h^{T} C_v^{T} D K^{T}
     covariance model: SABER
     saber central block:
-      saber block name: ID
-    # This will setup B = K D C_v C_h C_h^{T} C_v^{T} D K^{T}
+      saber block name: EXPLICIT_DIFFUSION
+      active variables: [tocn, socn, ssh]
+      geometry:
+        mom6_input_nml: mom_input.nml
+        fields metadata: ./fields_metadata.yaml
+      group mapping:
+      - name: ocean
+        variables: [tocn, socn, socn, ssh]
+      - name: ice
+        variables: [cicen]
+      read:
+        groups:
+        - name: ocean
+          horizontal:
+            filename: hz_ocean.nc
+          vertical:
+            filename: vt_ocean.nc
+        - name: ice
+          horizontal:
+            filename: hz_ice.nc
+
     linear variable change:
       input variables: [cicen, hicen, hsnon, socn, tocn, uocn, vocn, ssh]
       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

parm/soca/berror/soca_parameters_diffusion_hz.yaml

@@ -0,0 +1,37 @@
+geometry: &geom
+  mom6_input_nml: mom_input.nml
+  fields metadata: ./fields_metadata.yaml
+
+background:
+  read_from_file: 1
+  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]
+
+background error:
+  covariance model: SABER
+  saber central block:
+    saber block name: EXPLICIT_DIFFUSION
+    geometry: *geom
+    calibration:
+      normalization:
+        method: randomization
+        iterations: 10000
+
+      groups:
+      - name: ocean
+        horizontal:
+          from file:
+            filename: ocn.cor_rh.incr.0001-01-01T00:00:00Z.nc
+            variable name: ave_ssh
+        write:
+          filename: hz_ocean.nc
+
+      - name: ice
+        horizontal:
+          as gaussian: true
+          fixed value: 200000.0
+        write:
+          filename: hz_ice.nc

parm/soca/berror/soca_parameters_diffusion_vt.yaml

@@ -0,0 +1,30 @@
+geometry: &geom
+  mom6_input_nml: mom_input.nml
+  fields metadata: ./fields_metadata.yaml
+
+background:
+  read_from_file: 1
+  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]
+
+background error:
+  covariance model: SABER
+  saber central block:
+    saber block name: EXPLICIT_DIFFUSION
+    geometry: *geom
+    calibration:
+      normalization:
+        # NOTE, not actually used here, since the normalization spec is only used for hz
+        method: randomization  #< other option is "brute force"
+        iterations: 1000       #< in the real world you'll want to use 1e4 or so
+
+      groups:
+      - name: vt_ocean
+        vertical:
+          as gaussian: true
+          fixed value: 15.0
+        write:
+          filename: vt_ocean.nc

parm/soca/berror/soca_setcorscales.yaml

@@ -4,37 +4,13 @@ resolution:
 
 date: 0001-01-01T00:00:00Z
 
-corr variables: [socn, tocn, ssh, uocn, vocn, cicen, hicen, hsnon]
+corr variables: [ssh]
 
 scales:
   vert layers: 5 # in units of layer
-  socn:
-    rossby mult: 1.00
-    min grid mult: 2.0
-  tocn:
-    rossby mult: 1.00
-    min grid mult: 2.0
-  uocn:
-    rossby mult: 1.00
-    min grid mult: 2.0
-  vocn:
-    rossby mult: 1.00
-    min grid mult: 2.0
   ssh:
     rossby mult: 1.00
     min grid mult: 2.0
-  cicen:
-    rossby mult: 0.0
-    min grid mult: 2.0
-    min: 50.0
-  hicen:
-    rossby mult: 0.0
-    min grid mult: 4.0
-    min: 150000.0
-  hsnon:
-    rossby mult: 0.0
-    min grid mult: 4.0
-    min: 150000.0
 
 rh output:
   datadir: ./

scripts/exgdas_global_marine_analysis_bmat.sh

@@ -109,6 +109,30 @@ if [ $err -gt 0  ]; then
     exit $err
 fi
 
+################################################################################
+# Set localization scales for the hybrid en. var.
+$APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml
+export err=$?; err_chk
+if [ $err -gt 0  ]; then
+    exit $err
+fi
+
+################################################################################
+# Initialize diffusion blocks
+clean_yaml soca_parameters_diffusion_hz.yaml
+$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_hz.yaml
+export err=$?; err_chk
+if [ $err -gt 0  ]; then
+    exit $err
+fi
+
+clean_yaml soca_parameters_diffusion_vt.yaml
+$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x soca_parameters_diffusion_vt.yaml
+export err=$?; err_chk
+if [ $err -gt 0  ]; then
+    exit $err
+fi
+
 ################################################################################
 # Compute convolution coefs for L
 clean_yaml soca_bump_loc.yaml

scripts/exgdas_global_marine_analysis_prep.py

@@ -373,6 +373,33 @@ def find_clim_ens(input_date):
 locscales_yaml_dst = os.path.join(stage_cfg['stage_dir'], 'soca_setlocscales.yaml')
 FileHandler({'copy': [[locscales_yaml_src, locscales_yaml_dst]]}).sync()
 
+################################################################################
+# copy yaml for correlation length scales
+
+logging.info(f"---------------- generate soca_setcorscales.yaml")
+corscales_yaml_src = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_setcorscales.yaml')
+corscales_yaml_dst = os.path.join(stage_cfg['stage_dir'], 'soca_setcorscales.yaml')
+FileHandler({'copy': [[corscales_yaml_src, corscales_yaml_dst]]}).sync()
+
+################################################################################
+# copy yaml for diffusion initialization
+
+logging.info(f"---------------- generate soca_parameters_diffusion_hz.yaml")
+diffu_hz_yaml = os.path.join(anl_dir, 'soca_parameters_diffusion_hz.yaml')
+diffu_hz_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror')
+diffu_hz_yaml_template = os.path.join(berror_yaml_dir, 'soca_parameters_diffusion_hz.yaml')
+config = YAMLFile(path=diffu_hz_yaml_template)
+config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
+config.save(diffu_hz_yaml)
+
+logging.info(f"---------------- generate soca_parameters_diffusion_vt.yaml")
+diffu_vt_yaml = os.path.join(anl_dir, 'soca_parameters_diffusion_vt.yaml')
+diffu_vt_yaml_dir = os.path.join(gdas_home, 'parm', 'soca', 'berror')
+diffu_vt_yaml_template = os.path.join(berror_yaml_dir, 'soca_parameters_diffusion_vt.yaml')
+config = YAMLFile(path=diffu_vt_yaml_template)
+config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
+config.save(diffu_vt_yaml)
+
 ################################################################################
 # generate yaml for bump/nicas (used for correlation and/or localization)