diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 9be1317175d3905b19ccaced01554ff18b913130..0160478b5059a28c6944c6e0c493bb5f4b3d7989 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -970,6 +970,7 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
     f = h5py.File(filename, 'r')
     icing_alt = f['icing_altitude'][:]
     cld_top_hgt = f['cld_height_acha'][:, y_a:y_b, x_a:x_b]
+    cld_geo_dz = f['cld_geo_thick'][:, y_a:y_b, x_a:x_b]
     cld_phase = f['cloud_phase'][:, y_a:y_b, x_a:x_b]
 
     if day_night == 'DAY':
@@ -985,9 +986,9 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
     bt_11um = f_l1b['temp_11_0um_nom'][:, y_a:y_b, x_a:x_b]
 
     if icing:
-        mask, idxs, num_tested = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, day_night=day_night)
+        mask, idxs, num_tested = apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, day_night=day_night)
     else:
-        mask, idxs, num_tested = apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, day_night=day_night)
+        mask, idxs, num_tested = apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, sol_zen, sat_zen, day_night=day_night)
 
     keep_idxs = []
 
@@ -1037,7 +1038,7 @@ def run_qc(filename, filename_l1b, day_night='ANY', pass_thresh_frac=0.20, icing
     return len(icing_alt), len(mask), len(keep_idxs)
 
 
-def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
+def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
 
     if day_night == 'DAY':
         opd_thick_threshold = 20
@@ -1050,6 +1051,7 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
     num_obs = len(icing_alt)
     cld_mask = cld_mask.reshape((num_obs, -1))
     cld_top_hgt = cld_top_hgt.reshape((num_obs, -1))
+    cld_geo_dz = cld_geo_dz.reshape((num_obs, -1))
     cld_phase = cld_phase.reshape((num_obs, -1))
     cld_opd = cld_opd.reshape((num_obs, -1))
     bt_11um = bt_11um.reshape((num_obs, -1))
@@ -1072,7 +1074,8 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
         keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3)  # cloudy
         keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
         keep_2 = np.invert(np.isnan(bt_11um[i,]))
-        keep_3 = np.invert(np.isnan(cld_opd[i,]))
+        keep_3 = np.invert(np.isnan(cld_geo_dz[i,]))
+        # keep_3 = np.invert(np.isnan(cld_opd[i,]))
         keep = keep_0 & keep_1 & keep_2 & keep_3
         num_keep = np.sum(keep)
         if num_keep == 0:
@@ -1080,8 +1083,10 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
 
         # Test 1
         keep = np.where(keep, (cld_top_hgt[i,] + closeness) > icing_alt[i], False)
-        # if (np.sum(keep) / num_keep) < 0.80:
-        #     continue
+        if (np.sum(keep) / num_keep) < 0.80:
+            continue
+        num_keep = np.sum(keep)
+        keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,]) < icing_alt[i], False)
 
         # # Test2
         # keep = np.where(keep,
@@ -1105,7 +1110,7 @@ def apply_qc_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask,
     return mask, idxs, num_tested
 
 
-def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
+def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_geo_dz, cld_phase, cld_opd, cld_mask, bt_11um, solzen, satzen, day_night='ANY'):
 
     if day_night == 'DAY':
         opd_thick_threshold = 20
@@ -1118,6 +1123,7 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mas
     num_obs = len(icing_alt)
     cld_mask = cld_mask.reshape((num_obs, -1))
     cld_top_hgt = cld_top_hgt.reshape((num_obs, -1))
+    cld_geo_dz = cld_geo_dz.reshape((num_obs, -1))
     cld_phase = cld_phase.reshape((num_obs, -1))
     cld_opd = cld_opd.reshape((num_obs, -1))
     bt_11um = bt_11um.reshape((num_obs, -1))
@@ -1140,17 +1146,19 @@ def apply_qc_no_icing_pireps(icing_alt, cld_top_hgt, cld_phase, cld_opd, cld_mas
         keep_0 = np.logical_or(cld_mask[i,] == 2, cld_mask[i,] == 3)  # cloudy
         keep_1 = np.invert(np.isnan(cld_top_hgt[i,]))
         keep_2 = np.invert(np.isnan(bt_11um[i,]))
+        keep_3 = np.invert(np.isnan(cld_geo_dz[i,]))
         # keep_3 = np.invert(np.isnan(cld_opd[i,]))
-        # keep = keep_0 & keep_1 & keep_2 & keep_3
-        keep = keep_0 & keep_1 & keep_2
+        keep = keep_0 & keep_1 & keep_2 & keep_3
         num_keep = np.sum(keep)
         if num_keep == 0:
             continue
 
         # Exp
         keep = np.where(keep, (cld_top_hgt[i,] + closeness) > icing_alt[i], False)
-        #if (np.sum(keep) / num_keep) < 0.80:
-        #    continue
+        if (np.sum(keep) / num_keep) < 0.80:
+           continue
+        num_keep = np.sum(keep)
+        keep = np.where(keep, (cld_top_hgt[i,] - cld_geo_dz[i,]) < icing_alt[i], False)
 
         # keep = np.where(keep, np.logical_and(bt_11um[i,] > 228.0, bt_11um[i,] < 273.0), False)