diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 41335c3badf217394270b5ed8f68d4a484fb0328..bd98f4bb175a5ed92cf40c48cc1ac843f52727b7 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -10,7 +10,6 @@ ci:
   services:
     - docker:19.03.1-dind
   script:
-    - docker --version
     - docker login --username $CI_REGISTRY_USER --password $CI_REGISTRY_PASSWORD $CI_REGISTRY
     - cd ci
     - docker build --tag $CI_REGISTRY_IMAGE .
diff --git a/mvcm/__init__.py b/mvcm/__init__.py
index cf1769a6f4541a3edf53ad893e7e76eea2dc6e1e..256c8a2600f357e56f089341d24549be53c85be4 100644
--- a/mvcm/__init__.py
+++ b/mvcm/__init__.py
@@ -1,3 +1,3 @@
 """Set up version."""
 
-__version__ = "0.2.3"
+__version__ = "0.2.4"
diff --git a/mvcm/main.py b/mvcm/main.py
index 2e1a0e8142e2f7ba2eae4c9c9e970667cb6ced90..7049ec8d37815efc6f4f5ac15bf03eb76db84ff1 100644
--- a/mvcm/main.py
+++ b/mvcm/main.py
@@ -12,6 +12,7 @@ import time
 
 import numpy as np
 import xarray as xr
+from netCDF4 import Dataset
 from pkg_resources import get_distribution  # type: ignore
 from rich.logging import RichHandler
 
@@ -217,7 +218,11 @@ def main(
                     return
             logger.info(f"All bands found in file {file_names['IMG02']}. The code will run.")
 
-    # logger.info(f"Reading VNP02: {file_names['MOD02']}")
+    with Dataset(file_names["MOD03"]) as f:
+        # time_coverage_start = f.getncattr("time_coverage_start")
+        # time_coverage_end = f.getncattr("time_coverage_end")
+        # orbit_number = f.getncattr("OrbitNumber")
+        attrs = {attr: f.getncattr(attr) for attr in f.ncattrs()}
 
     viirs_data = rd.get_data(satellite, sensor, file_names, hires=use_hires)
 
@@ -483,33 +488,8 @@ def main(
     #                                                          cmin,
     #                                                          bits['r06'])[idx]
 
-    # cmin_3d[i, :, :] = cmin
-    # i += 1
-
-    # cmin_final = np.min(cmin_3d, axis=0)
-    # bit['15'] = np.max(temp_bit, axis=0)
     cmin_final = cmin
 
-    """
-    idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['water'] == 1) &
-                     (scene_flags['ice'] == 0))
-    cmin[idx] = restore.shallow_water(sh_water_bits, cmin)[idx]
-
-    idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) &
-                     (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) &
-                     (cmin <= 0.95))
-    cmin[idx] = restore.land(land_bits, cmin)[idx]
-
-    idx = np.nonzero((scene_flags['day'] == 1) & (scene_flags['land'] == 1) &
-                     (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) &
-                     (scene_flags['coast'] == 1))
-    cmin[idx] = restore.coast(coast_bits, cmin)[idx]
-
-    idx = np.nonzero((scene_flags['night'] == 1) & (scene_flags['land'] == 1) &
-                     (scene_flags['snow'] == 0) & (scene_flags['ice'] == 0) &
-                     (cmin <= 0.95))
-    cmin[idx] = restore.land_night(land_night_bits, cmin)[idx]
-    """
     # bits translation MVCM-python -> MVCM-C
     # 01  13  test_11um
     # 02  27  surface_temperature_test
@@ -523,29 +503,6 @@ def main(
     # 15  16  test_1_38um_high_clouds
     # 16  17  thin_cirrus_4_12um_BTD_test
 
-    #     if scene_name in ['Ocean_Day']:
-    #         # total_bit = np.full(viirs_data.M01.shape, 4)
-    #         total_bit = bit1 + bit2 + bit4 + bit11
-    #         sunglint_angle = thresholds['Sun_Glint']['bounds'][3]
-    #         scene_flags = scn.find_scene(viirs_data, sunglint_angle)
-    #
-    #         # idx = np.nonzero((scene_flags['water'] == 1) &
-    #                            (scene_flags['ice'] == 0) &
-    #         #                  (scene_flags['uniform'] == 1) &
-    #                            (cmin <= 0.99) & (cmin >= 0.05))
-    #         # cmin[idx] = restoral.spatial(viirs_data, thresholds['Sun_Glint'],
-    #                                        scene_flags, cmin)[idx]
-    #
-    #         idx = np.nonzero((scene_flags['water'] == 1) &
-    #                          (scene_flags['sunglint'] == 1) &
-    #                          (scene_flags['uniform'] == 1) & (cmin <= 0.95))
-    #         cmin[idx] = restoral.sunglint(viirs_data, thresholds['Sun_Glint'],
-    #                                       total_bit, cmin)[idx]
-    #
-    #     # MVCM = tst.ComputeTests(viirs_data, scene_name, thresholds)
-    #     # cmin = MVCM.clear_sky_restoral(cmin)
-    #
-    # # return cmin
     csc = cmin_final
     cloud_mask = np.zeros(csc.shape)
     cloud_mask[csc > 0.99] = 3
@@ -560,6 +517,7 @@ def main(
         "longitude": {"dims": ("x", "y"), "data": viirs_data.longitude.values},
         "confidence": {"dims": ("x", "y"), "data": cmin_final},
         "cloud_mask": {"dims": ("x", "y"), "data": cloud_mask},
+        "sensor_zenith": {"dims": ("x", "y"), "data": viirs_data.sensor_zenith.values},
     }
     d_debug = {
         "ndvi": {"dims": ("x", "y"), "data": viirs_data.ndvi.values},
@@ -632,54 +590,17 @@ def main(
 
     ds_out = xr.Dataset.from_dict(out_dict)
 
-    comp = dict(zlib=True, complevel=5)
+    comp = {"zlib": True, "complevel": 5}
     if type(ds_out) is xr.Dataset:
         for var in ds_out:
             ds_out[var].encoding.update(comp)
     if type(ds_out) is xr.DataArray:
         ds_out.encoding.update(comp)
 
-    # if use_hires is False:
-    #     ds_out.to_netcdf('test_w_restorals.nc')
-    # else:
-    #     ds_out.to_netcdf('test_hires.nc')
-    ds_out.to_netcdf(out_file)
-
-
-"""
-def test_main():
-    rad1 = [
-        [255, 260, 265, 248, 223],
-        [278, 285, 270, 268, 256],
-        [275, 273, 266, 254, 259],
-    ]
-    rad2 = [
-        [270, 273, 271, 268, 265],
-        [277, 286, 275, 277, 269],
-        [280, 281, 272, 270, 267],
-    ]
-    thresh_file = (
-        "/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml"
-    )
-
-    with open(thresh_file) as f:
-        text = f.read()
-
-    rad1 = np.array(rad1)
-    rad2 = np.array(rad2)
+    for attr in attrs:
+        ds_out.attrs[attr] = attrs[attr]
 
-    confidence = np.zeros((2, rad1.shape[0], rad1.shape[1]))
-
-    thresholds = yml.safe_load(text)
-
-    confidence[0, :, :] = tests.test_11um(rad1, thresholds["Daytime_Ocean"])
-    confidence[1, :, :] = tests.test_11_4diff(rad1, rad2, thresholds["Daytime_Ocean"])
-
-    print(f"Confidence[0,:,:]: \n {confidence[0, :, :]}")
-    print(f"Confidence[1,:,:]: \n {confidence[1, :, :]}")
-
-    return confidence
-"""
+    ds_out.to_netcdf(out_file)
 
 
 if __name__ == "__main__":
diff --git a/mvcm/preprocess_thresholds.py b/mvcm/preprocess_thresholds.py
index 20aa06786b26e2f040e33e3fa13a10aefdaee0cf..571f8ede7b0e326e941640cb94a87dbc5e65883b 100644
--- a/mvcm/preprocess_thresholds.py
+++ b/mvcm/preprocess_thresholds.py
@@ -16,7 +16,7 @@ _DTR = np.pi / 180
 logger = logging.getLogger(__name__)
 
 
-def prepare_11_12um_thresholds(thresholds: Dict, dim1: int) -> Dict:
+def prepare_11_12um_thresholds(thresholds: dict, dim1: int) -> dict:
     """Prepare 11-12um test thresholds."""
     coeff_values = np.empty((dim1, 2))
     coeff_values[:, 0] = np.full(dim1, thresholds["coeffs"][0])
@@ -108,9 +108,11 @@ def thresholds_11_12um(
         elif scene in ["Day_Snow"]:
             hicut = locut - (thr_dict["cmult"] * locut)
         else:
-            raise ValueError("Scene name not valid")
+            err_msg = "Scene name not valid"
+            raise ValueError(err_msg)
     else:
-        raise ValueError("Scene name not valid")
+        err_msg = "Scene name not valid"
+        raise ValueError(err_msg)
 
     thr_out = np.dstack((locut, midpt, hicut, np.ones(locut.shape), np.ones(locut.shape)))
     return np.squeeze(thr_out.T)