Make aerosol diag B able to handle multiple resolutions + add in more config choices

parm/aero/berror/aero_diagb.yaml.j2

@@ -1,3 +1,5 @@
+{% set offset_td = "+6H" | to_timedelta %}
+{% set background_time = current_cycle | add_to_datetime(offset_td) %}
 geometry:
   fms initialization:
     namelist filename: ./fv3jedi/fmsmpp.nml
@@ -6,13 +8,13 @@ geometry:
   layout:
   - {{ layout_x }}
   - {{ layout_y }}
-  npx: {{ npx_ges }}
-  npy: {{ npy_ges }}
-  npz: {{ npz_ges }}
+  npx: {{ npx_anl }}
+  npy: {{ npy_anl }}
+  npz: {{ npz_anl }}
   field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
-date: '{{ current_cycle | to_isotime }}' 
+date: '{{ background_time | to_isotime }}'
 background:
-  datetime: '{{ current_cycle | to_isotime }}' 
+  datetime: '{{ background_time | to_isotime }}'
   filetype: fms restart
   state variables: [mass_fraction_of_sulfate_in_air,
                     mass_fraction_of_hydrophobic_black_carbon_in_air,
@@ -25,22 +27,37 @@ background:
                     mass_fraction_of_sea_salt002_in_air, mass_fraction_of_sea_salt003_in_air,
                     mass_fraction_of_sea_salt004_in_air]
   datapath: ./bkg
-  filename_core: '{{ current_cycle | to_fv3time }}.fv_core.res.nc'
-  filename_trcr: '{{ current_cycle | to_fv3time }}.fv_tracer.res.nc'
-  filename_cplr: '{{ current_cycle | to_fv3time }}.coupler.res'
+  filename_core: '{{ background_time | to_fv3time }}.anlres.fv_core.res.nc'
+  filename_trcr: '{{ background_time | to_fv3time }}.anlres.fv_tracer.res.nc'
+  filename_cplr: '{{ background_time | to_fv3time }}.anlres.coupler.res'
 
 background error:
   filetype: fms restart
   datapath: ./stddev
   filename_trcr: stddev.fv_tracer.res.nc
   filename_cplr: stddev.coupler.res
 
-climate background error stddev:
-  filetype: fms restart
-  skip coupler file: true
-  datapath: ./clm_stddev
-  filename_trcr: stddev.fv_tracer.res.nc
-  filename_cplr: stddev.coupler.res
+climate background error:
+  geometry:
+    fms initialization:
+      namelist filename: ./fv3jedi/fmsmpp.nml
+      field table filename: ./fv3jedi/field_table
+    akbk: ./fv3jedi/akbk.nc4
+    layout:
+    - {{ layout_x }}
+    - {{ layout_y }}
+    npx: {{ npx_clim_b }}
+    npy: {{ npy_clim_b }}
+    npz: {{ npz_ges }}
+    field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
+  climate background error stddev:
+    filetype: fms restart
+    skip coupler file: true
+    datapath: ./clm_stddev
+    filename_trcr: stddev.fv_tracer.res.nc
+    filename_cplr: stddev.coupler.res
+  diagb weight: {{ aero_diagb_weight }}
+  staticb rescaling factor: {{ aero_staticb_rescaling_factor }}
 
 variables:
   name:
@@ -59,11 +76,9 @@ variables:
   - mass_fraction_of_sea_salt003_in_air
   - mass_fraction_of_sea_salt004_in_air
 
-rescale: 2.0    # rescales the filtered std. dev. by "rescale"
-diagb weight: 0.5
-staticb rescaling factor: 2.0
-number of halo points: 4
-number of neighbors: 16
+rescale: {{ aero_diagb_rescale }}    # rescales the filtered std. dev. by "rescale"
+number of halo points: {{ aero_diagb_n_halo }}
+number of neighbors: {{ aero_diagb_n_neighbors }}
 simple smoothing:
-  horizontal iterations: 0 
-  vertical iterations: 0
+  horizontal iterations: {{ aero_diagb_smooth_horiz_iter }}
+  vertical iterations: {{ aero_diagb_smooth_vert_iter }}

parm/aero/berror/aero_diffusionparm.yaml.j2

@@ -1,3 +1,5 @@
+{% set offset_td = "+6H" | to_timedelta %}
+{% set background_time = current_cycle | add_to_datetime(offset_td) %}
 geometry:
   fms initialization:
     namelist filename: ./fv3jedi/fmsmpp.nml
@@ -6,18 +8,18 @@ geometry:
   layout:
   - {{ layout_x }}
   - {{ layout_y }}
-  npx: {{ npx_ges }}
-  npy: {{ npy_ges }}
-  npz: {{ npz_ges }}
+  npx: {{ npx_anl }}
+  npy: {{ npy_anl }}
+  npz: {{ npz_anl }}
   field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
-date: '{{ current_cycle | to_isotime }}' 
+date: '{{ background_time | to_isotime }}' 
 background:
-  datetime: '{{ current_cycle | to_isotime }}' 
+  datetime: '{{ background_time | to_isotime }}' 
   filetype: fms restart
   datapath: ./bkg
-  filename_core: '{{ current_cycle | to_fv3time }}.fv_core.res.nc'
-  filename_trcr: '{{ current_cycle | to_fv3time }}.fv_tracer.res.nc'
-  filename_cplr: '{{ current_cycle | to_fv3time }}.coupler.res'
+  filename_core: '{{ background_time | to_fv3time }}.anlres.fv_core.res.nc'
+  filename_trcr: '{{ background_time | to_fv3time }}.anlres.fv_tracer.res.nc'
+  filename_cplr: '{{ background_time | to_fv3time }}.anlres.coupler.res'
   state variables: [mass_fraction_of_sulfate_in_air,
                     mass_fraction_of_hydrophobic_black_carbon_in_air,
                     mass_fraction_of_hydrophilic_black_carbon_in_air,

parm/aero/berror/aero_interp.yaml.j2

@@ -0,0 +1,50 @@
+{% set offset_td = "+6H" | to_timedelta %}
+{% set background_time = current_cycle | add_to_datetime(offset_td) %}
+input geometry:
+  fms initialization:
+    namelist filename: ./fv3jedi/fmsmpp.nml
+    field table filename: ./fv3jedi/field_table
+  akbk: ./fv3jedi/akbk.nc4
+  layout:
+  - {{ layout_x }}
+  - {{ layout_y }}
+  npx: {{ npx_ges }}
+  npy: {{ npy_ges }}
+  npz: {{ npz_ges }}
+  field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
+output geometry:
+  fms initialization:
+    namelist filename: ./fv3jedi/fmsmpp.nml
+    field table filename: ./fv3jedi/field_table
+  akbk: ./fv3jedi/akbk.nc4
+  layout:
+  - {{ layout_x }}
+  - {{ layout_y }}
+  npx: {{ npx_anl }}
+  npy: {{ npy_anl }}
+  npz: {{ npz_anl }}
+  field metadata override: ./fv3jedi/fv3jedi_fieldmetadata_restart.yaml
+states:
+- input:
+    datetime: '{{ background_time | to_isotime }}'
+    filetype: fms restart
+    state variables: [mass_fraction_of_sulfate_in_air,
+                        mass_fraction_of_hydrophobic_black_carbon_in_air,
+                        mass_fraction_of_hydrophilic_black_carbon_in_air,
+                        mass_fraction_of_hydrophobic_organic_carbon_in_air,
+                        mass_fraction_of_hydrophilic_organic_carbon_in_air,
+                        mass_fraction_of_dust001_in_air, mass_fraction_of_dust002_in_air,
+                        mass_fraction_of_dust003_in_air, mass_fraction_of_dust004_in_air,
+                        mass_fraction_of_dust005_in_air, mass_fraction_of_sea_salt001_in_air,
+                        mass_fraction_of_sea_salt002_in_air, mass_fraction_of_sea_salt003_in_air,
+                        mass_fraction_of_sea_salt004_in_air]
+    datapath: ./bkg
+    filename_core: '{{ background_time | to_fv3time }}.fv_core.res.nc'
+    filename_trcr: '{{ background_time | to_fv3time }}.fv_tracer.res.nc'
+    filename_cplr: '{{ background_time | to_fv3time }}.coupler.res'
+  output:
+    filetype: fms restart
+    datapath: ./bkg
+    filename_core: 'anlres.fv_core.res.nc'
+    filename_trcr: 'anlres.fv_tracer.res.nc'
+    filename_cplr: 'anlres.coupler.res'

utils/chem/chem_diagb.h

@@ -132,6 +132,7 @@ namespace gdasapp {
 
       // Loop through variables
       for (auto & var : chemVars.variables()) {
+        nodeColumns.haloExchange(xbFs[var]);
         oops::Log::info() << "====================== std dev for " << var << std::endl;
         auto bkg = atlas::array::make_view<double, 2>(xbFs[var]);
         auto stdDevBkg = atlas::array::make_view<double, 2>(bkgErrFs[var]);
@@ -230,54 +231,66 @@ namespace gdasapp {
       }
 
       // Rescale
-      double rescale;
-      fullConfig.get("rescale", rescale);
-      util::multiplyFieldSet(bkgErrFs, rescale);
-
-      // Hybrid diagb_weight coefficient
-      // std = diagb_weight*diagb_std + (1-diagb_weight)*Climat_std
-      double diagb_weight;
-      fullConfig.get("diagb weight", diagb_weight);
-      // Initialize and read the climatological background error standard deviation field
-      oops::Log::info() << "====================== read climat bkg error std dev" << std::endl;
-      fv3jedi::Increment ClimBkgErrorStdDev(geom, chemVars, cycleDate);
-      ClimBkgErrorStdDev.zero();
-      const eckit::LocalConfiguration ClimBkgErrorStdDevConfig(fullConfig,
-      "climate background error stddev");
-      ClimBkgErrorStdDev.read(ClimBkgErrorStdDevConfig);
-      atlas::FieldSet ClimBkgErrorStdDevFs;
-      ClimBkgErrorStdDev.toFieldSet(ClimBkgErrorStdDevFs);
-
-     // Replace negative values with zeros
-     for (const auto &var : chemVars.variables()) {
-        auto ClimBkgErrView = atlas::array::make_view<double, 2>(ClimBkgErrorStdDevFs[var]);
-        for (atlas::idx_t jnode = 0; jnode < ClimBkgErrorStdDevFs[var].shape(0); ++jnode) {
-          for (atlas::idx_t level = 0; level < ClimBkgErrorStdDevFs[var].shape(1); ++level) {
-            ClimBkgErrView(jnode, level) = std::max(ClimBkgErrView(jnode, level), 0.0);
-          }
-        }
+      if (fullConfig.has("rescale")) {
+        double rescale;
+        fullConfig.get("rescale", rescale);
+        util::multiplyFieldSet(bkgErrFs, rescale);
       }
 
-      // Staticb rescale
-      double rescale_staticb;
-      fullConfig.get("staticb rescaling factor", rescale_staticb);
+      bkgErr.fromFieldSet(bkgErrFs);
+
+      // Hybrid B option
+      if (fullConfig.has("climate background error")) {
+        const eckit::LocalConfiguration ClimBConfig(fullConfig, "climate background error");
+        // Hybrid diagb_weight coefficient
+        // std = diagb_weight*diagb_std + (1-diagb_weight)*Climat_std
+        double diagb_weight;
+        ClimBConfig.get("diagb weight", diagb_weight);
+        // Initialize and read the climatological background error standard deviation field
+        oops::Log::info() << "====================== read climat bkg error std dev" << std::endl;
+        const eckit::LocalConfiguration climBGeomConfig(ClimBConfig, "geometry");
+        const fv3jedi::Geometry climBGeom(climBGeomConfig, this->getComm());
+        fv3jedi::Increment ClimBkgErrorStdDevOrig(climBGeom, chemVars, cycleDate);
+        ClimBkgErrorStdDevOrig.zero();
+        const eckit::LocalConfiguration ClimBkgErrorStdDevConfig(ClimBConfig,
+        "climate background error stddev");
+        ClimBkgErrorStdDevOrig.read(ClimBkgErrorStdDevConfig);
+        // interpolate to background resolution
+        fv3jedi::Increment ClimBkgErrorStdDev(geom, ClimBkgErrorStdDevOrig);
+        atlas::FieldSet ClimBkgErrorStdDevFs;
+        ClimBkgErrorStdDev.toFieldSet(ClimBkgErrorStdDevFs);
+
+        // Replace negative values with zeros
+        for (const auto &var : chemVars.variables()) {
+            auto ClimBkgErrView = atlas::array::make_view<double, 2>(ClimBkgErrorStdDevFs[var]);
+            for (atlas::idx_t jnode = 0; jnode < ClimBkgErrorStdDevFs[var].shape(0); ++jnode) {
+              for (atlas::idx_t level = 0; level < ClimBkgErrorStdDevFs[var].shape(1); ++level) {
+                ClimBkgErrView(jnode, level) = std::max(ClimBkgErrView(jnode, level), 0.0);
+              }
+            }
+          }
+
+        // Staticb rescale
+        double rescale_staticb;
+        ClimBConfig.get("staticb rescaling factor", rescale_staticb);
 
-      // Combine diagb and climatological background errors
-      fv3jedi::Increment stddev_hybrid(geom, chemVars, cycleDate);
-      stddev_hybrid.zero();
+        // Combine diagb and climatological background errors
+        fv3jedi::Increment stddev_hybrid(geom, chemVars, cycleDate);
+        stddev_hybrid.zero();
 
-      // Convert FieldSets to States for accumulation
-      fv3jedi::State bkgErrState(geom, chemVars, cycleDate);
-      bkgErrState.fromFieldSet(bkgErrFs);
+        // Convert FieldSets to States for accumulation
+        fv3jedi::State bkgErrState(geom, chemVars, cycleDate);
+        bkgErrState.fromFieldSet(bkgErrFs);
 
-      fv3jedi::State ClimBkgErrorStdDevState(geom, chemVars, cycleDate);
-      ClimBkgErrorStdDevState.fromFieldSet(ClimBkgErrorStdDevFs);
+        fv3jedi::State ClimBkgErrorStdDevState(geom, chemVars, cycleDate);
+        ClimBkgErrorStdDevState.fromFieldSet(ClimBkgErrorStdDevFs);
 
-      // Accumulate the fields with the given weights
-      stddev_hybrid.accumul(diagb_weight, bkgErrState);
-      stddev_hybrid.accumul((1.0 - diagb_weight)*rescale_staticb, ClimBkgErrorStdDevState);
-      // Use the hybrid increment
-      bkgErr = stddev_hybrid;
+        // Accumulate the fields with the given weights
+        stddev_hybrid.accumul(diagb_weight, bkgErrState);
+        stddev_hybrid.accumul((1.0 - diagb_weight)*rescale_staticb, ClimBkgErrorStdDevState);
+        // Use the hybrid increment
+        bkgErr = stddev_hybrid;
+      }
 
       // Save the background error
       const eckit::LocalConfiguration bkgErrorConfig(fullConfig, "background error");