Apply land fraction as land mask for IMS snow depth preprocessing

parm/snow/jcb-base.yaml.j2

@@ -84,7 +84,7 @@ snow_obsdatain_path: "{{DATA}}/obs"
 snow_obsdatain_prefix: "{{OPREFIX}}"
 snow_obsdatain_suffix: ".tm00.bufr_d"
 snow_obsdataout_path: "{{DATA}}/diags"
-snow_prepobs_path: "{{DATA}}/prep"
+snow_prepobs_path: "{{ snow_prepobs_path }}"
 snow_obsdataout_prefix: diag_
 snow_obsdataout_suffix: "_{{ current_cycle | to_YMDH }}.nc"
 ims_scf_obs_suffix: "{{ ims_scf_obs_suffix }}"

utils/land/gdas_fv3jedi_calc_scf_to_ioda.cc

@@ -125,6 +125,10 @@ void gdasapp::CalcSCFtoIODA::calc_fcst_snow_density(fv3jedi::State & bkgState,
       bkg_density(jnode, 0) =
         std::max(80.0, std::min(120.0, tmp_density)) / 1000.0;
     }
+    // where (density < 0.0001) density = 0.08
+    if (bkg_density(jnode, 0) < 0.0001) {
+      bkg_density(jnode, 0) = 0.08;
+    }
   }
   // put the new density values back in the state
   bkgState.fromFieldSet(xBfs);
@@ -146,10 +150,12 @@ void gdasapp::CalcSCFtoIODA::calc_fcst_snow_cover_fraction(fv3jedi::State & bkgS
   auto bkg_snow_den = atlas::array::make_view<double, 2>(xBfs["snowDensity"]);
   auto bkg_snd = atlas::array::make_view<double, 2>(xBfs["totalSnowDepth"]);
   auto bkg_scf = atlas::array::make_view<double, 2>(xBfs["surface_snow_area_fraction"]);
+  auto bkg_landfrac = atlas::array::make_view<double, 2>(xBfs["fraction_of_land"]);
   // now compute snow cover fraction
   for (atlas::idx_t jnode = 0; jnode < xBfs["totalSnowDepth"].shape(0); ++jnode) {
     int vetfcs = static_cast<int>(bkg_vtype(jnode, 0));
-    if (vetfcs > 0) {
+    // if landcover fraction > 50%, but not land ice
+    if (bkg_landfrac(jnode, 0) > 0.5 && vetfcs != vtype_landice) {
       if (bkg_snd(jnode, 0) > 0.0f) {
         // snow is present, compute snow cover fraction
         float snowh = bkg_snd(jnode, 0)*0.001f;  // convert mm to m
@@ -660,6 +666,7 @@ void gdasapp::CalcSCFtoIODA::IMSscf::calcIMSsd(fv3jedi::State &state,
   // Get the vegetation type field from the state
   auto bkg_vtype = atlas::array::make_view<double, 2>(xBfs["vtype"]);
   auto bkg_snow_den = atlas::array::make_view<double, 2>(xBfs["snowDensity"]);
+  auto bkg_landfrac = atlas::array::make_view<double, 2>(xBfs["fraction_of_land"]);
   const auto bkg_idx =
       atlas::array::make_view<atlas::gidx_t, 1>(geom.functionSpace().global_index());
   std::vector<int> indices = geom.get_indices();
@@ -670,10 +677,10 @@ void gdasapp::CalcSCFtoIODA::IMSscf::calcIMSsd(fv3jedi::State &state,
       size_t local_j = fv3_j - (indices[2] - 1);
       size_t local_i = fv3_i - (indices[0] - 1);
       atlas::idx_t jnode = ((local_j)*(npx) + (local_i + 1)) - 1;
+      // if we have IMS data at this point
       if (abs(this->scfIMS[fv3_j][fv3_i] - nodata_float) > nodata_tol) {
-        // if we have IMS data at this point
-        if (bkg_vtype(jnode, 0) > 0) {
-          // if the model has land at this point
+        // if the model landcover fraction > 50%, but not land ice
+        if (bkg_landfrac(jnode, 0) > 0.5 && bkg_vtype(jnode, 0) != vtype_landice) {
           if (this->scfIMS[fv3_j][fv3_i] < 0.5f) {
             // if the IMS SCF is less than 0.5, set snow depth to 0
             this->sndIMS[fv3_j][fv3_i] = 0.0f;

utils/land/gdas_fv3jedi_calc_scf_to_ioda.h

@@ -42,6 +42,7 @@ namespace gdasapp {
                                                  // use this value when calculating SD
                                                  // to represent "full" coverage
       static constexpr float sndIMS_max = 300.0f;  // maximum snow depth derived from IMS
+      static constexpr int vtype_landice = 15;
       void netcdf_err(int error, const std::string &msg);
 
       // Add members and methods as needed
@@ -74,6 +75,7 @@ namespace gdasapp {
     float oberr_snd = 80.0f;
     static constexpr float nodata_float = -999.0f;
     static constexpr float nodata_tol = 0.1f;  // Tolerance for nodata checks
+    static constexpr int vtype_landice = 15;
 
    public:
     void run();