diff --git a/mvcm/main.py b/mvcm/main.py
index 51c284a11d06dad48554131ea8877f5963c6e405..7b78f3d1b1e2eb062d94b2cf432bf1c129146fa6 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -297,18 +297,16 @@ def main(
         cmin_g1, bits["01"] = my_scene.test_11um(m15_name, cmin_g1, bits["01"])
         # this test needs to be implemented properly,
         # for now is commented since it's only used for night granules
-        # cmin_g1, bits["02"] = my_scene.surface_temperature_test(
-        #     "", viirs_data, cmin_g1, bits["02"]
-        # )
+        # cmin_g1, bits["02"] = my_scene.surface_temperature_test("", viirs_data, cmin_g1, bits["02"])
         cmin_g1, bits["03"] = my_scene.sst_test("M15", "M16", cmin_g1, bits["03"])
 
         # Group 2
         cmin_g2, bits["04"] = my_scene.bt_diff_86_11um("M14-M15", cmin_g2, bits["04"])
         cmin_g2, bits["05"] = my_scene.test_11_12um_diff("M15-M16", cmin_g2, bits["05"])
         cmin_g2, bits["06"] = my_scene.variability_11um_test(m15_name, cmin_g2, bits["06"])
-        # cmin_g2, bits["07"] = my_scene.bt_difference_11_4um_test_ocean(
-        #     "M15-M12", cmin_g2, bits["07"]
-        # )
+        cmin_g2, bits["07"] = my_scene.bt_difference_11_4um_test_ocean(
+            "M15-M12", cmin_g2, bits["07"]
+        )
         # cmin_g2, bits["08"] = my_scene.bt_difference_11_4um_test_land(
         #     "M15-M16", "M15-M12", cmin_g2, bits["08"]
         # )
@@ -513,7 +511,7 @@ def main(
         "integer_cloud_mask": {"dims": ("x", "y"), "data": integer_cloud_mask},
     }
 
-    debug_data = {
+    old_debug_data = {
         "ndvi": {"dims": ("x", "y"), "data": viirs_data.ndvi.values},
         "sunglint_bits": {"dims": ("x", "y"), "data": sunglint_bits},
         "eco": {"dims": ("x", "y"), "data": viirs_data.ecosystem.values},
@@ -563,6 +561,8 @@ def main(
         "M15-M12": {"dims": ("x", "y"), "data": viirs_data["M15-M12"].values},
         "M15-M16": {"dims": ("x", "y"), "data": viirs_data["M15-M16"].values},
         "Ocean_Day": {"dims": ("x", "y"), "data": viirs_data.Ocean_Day.values},
+        "Land_Day": {"dims": ("x", "y"), "data": viirs_data.Land_Day.values},
+        "Day_Snow": {"dims": ("x", "y"), "data": viirs_data.Day_Snow.values},
         "Polar_Day_Ocean": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Ocean.values},
         "Polar_Day_Land": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Land.values},
         "Polar_Day_Desert": {"dims": ("x", "y"), "data": viirs_data.Polar_Day_Desert.values},
@@ -575,13 +575,32 @@ def main(
     }
 
     if debug is True:
-        output_dict.update(debug_data)
+        # debug_output.update(old_debug_data)
         sf = {k: {"dims": ("x", "y"), "data": scene_flags[k]} for k in scene_flags}
-        output_dict = output_dict | sf
+        old_debug_data = old_debug_data | sf
+
+    debug_output = xr.Dataset.from_dict(old_debug_data)
+    debug_output.to_netcdf("/ships19/hercules/pveglio/debug_data.nc", mode="w")
+
+    all_bits = np.zeros((bit["01"].shape[0], bit["01"].shape[1], 22))
+    for i in range(len(_bitlist)):
+        all_bits[:, :, i] = bit[_bitlist[i]]
+
+    all_bands = np.zeros((viirs_data.M01.shape[0], viirs_data.M01.shape[1], 16))
+    for j in range(len(_mod_bands)):
+        all_bands[:, :, j] = viirs_data[_mod_bands[j]].values
+
+    debug_data = {
+        "bits": {"dims": ("x", "y", "testbits"), "data": all_bits},
+        "bands": {"dims": ("x", "y", "nbands"), "data": all_bands},
+    }
+
+    if debug is True:
+        output_dict.update(debug_data)
 
     output_data = xr.Dataset.from_dict(output_dict)
 
-    save_output(output_data, attrs, out_file, compression=5)
+    save_output(output_data, attrs, out_file, compression=5, debug=debug)
 
     # comp = {"zlib": True, "complevel": 5}
     # if type(ds_out) is xr.Dataset:
diff --git a/mvcm/spectral_tests.py b/mvcm/spectral_tests.py
index 999f1c69ac3b4266d4d770f4d0263011207769cb..14f03549dfceaa56718a212c2b40c6347e8576ae 100644
--- a/mvcm/spectral_tests.py
+++ b/mvcm/spectral_tests.py
@@ -582,7 +582,16 @@ class CloudTests:
                 )
                 bits["test"][idx] = 1
                 kwargs["confidence"][self.scene_idx] = conf.conf_test_dble(rad[self.scene_idx], thr)
-
+            elif self.scene_name in ["Day_Snow"]:
+                f_rad = np.array(rad[self.scene_idx], dtype=np.float32)
+                locut = np.array(thr[0, :], dtype=np.float32)
+                midpt = np.array(thr[1, :], dtype=np.float32)
+                hicut = np.array(thr[2, :], dtype=np.float32)
+                power = np.array(thr[3, :], dtype=np.float32)
+                bits["test"][self.scene_idx] = 1
+                kwargs["confidence"][self.scene_idx] = anc.py_conf_test(
+                    f_rad, locut, hicut, power, midpt
+                )
             else:
                 kwargs["confidence"][self.scene_idx] = conf.conf_test_new(rad[self.scene_idx], thr)
 
diff --git a/mvcm/write_output.py b/mvcm/write_output.py
index 98aebe93c9277d0354e71d8aa09ffbf84d3796fe..b74ddbf4904fa9ad31949f2fe577505c192efe68 100644
--- a/mvcm/write_output.py
+++ b/mvcm/write_output.py
@@ -23,7 +23,7 @@ _attributes_list = [
 logger = logging.getLogger(__name__)
 
 
-def save_output(data: xr.Dataset, attrs: dict, fname: str, compression: int) -> None:
+def save_output(data: xr.Dataset, attrs: dict, fname: str, compression: int, debug: bool) -> None:
     """Save the output data in a netcdf file.
 
     Parameters
@@ -170,6 +170,16 @@ def save_output(data: xr.Dataset, attrs: dict, fname: str, compression: int) ->
         },
     }
 
+    if debug is True:
+        geophysical_data["Bits"] = {
+            "dims": ("number_of_lines", "number_of_pixels", "bits"),
+            "data": data.bits.values,
+        }
+        geophysical_data["Bands"] = {
+            "dims": ("number_of_lines", "number_of_pixels", "nbands"),
+            "data": data.bands.values,
+        }
+
     logger.debug("geophysical_data dictionary created successfully")
 
     attributes = {