diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 4ac49859b09e1ab09243200bb5989c0d0f7a3714..12b13ae6365d65421ae353b9de4b50874c9b3881 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -878,55 +878,30 @@ def analyze2(filename, filename_l1b):
 
     iint = f['icing_intensity'][:]
     icing_alt = f['flight_altitude'][:]
-    cld_top_hgt = f['cld_height_acha'][:, :, :]
-    cld_top_tmp = f['cld_temp_acha'][:, :, :]
-    sc_cld_frac = f['supercooled_cloud_fraction'][:, :, :]
-    cld_mask = f['cloud_mask'][:, :, :]
-    cld_reff = f['cld_reff_dcomp'][:, :, :]
-    cld_emiss = f['cld_emiss_acha'][:, :, :]
-    # cld_phase = f['cloud_phase'][:, :, :]
-    # cld_opd = f['cld_opd_dcomp'][:, :, :]
-    # solzen = f['solar_zenith_angle'][:, :, :]
-
-    bt_10_4 = f_l1b['temp_10_4um_nom'][:, :, :]
-
-    cth = np.reshape(cld_top_hgt, [iint.size, -1])
-    ctt = np.reshape(cld_top_tmp, [iint.size, -1])
-    scf = np.reshape(sc_cld_frac, [iint.size, -1])
-    bt = np.reshape(bt_10_4, [iint.size, -1])
-    cmsk = np.where(cld_mask > 1, 1, 0)
-    cmsk = np.reshape(cmsk, [iint.size, -1])
-    cld_frac = np.sum(cmsk, axis=1) / 256
-
-    cth = np.nanmedian(cth, axis=1)
-    ctt = np.nanmedian(ctt, axis=1)
-    scf = np.nanmedian(scf, axis=1)
-    bt = np.nanmedian(bt, axis=1)
-    # cth = np.nanmax(cth, axis=1)
-    # ctt = np.nanmax(ctt, axis=1)
-    # scf = np.nanmax(scf, axis=1)
-    # cth = np.nanmin(cth, axis=1)
-    # ctt = np.nanmin(ctt, axis=1)
-    # scf = np.nanmin(scf, axis=1)
-
-    keep_0 = icing_alt < 5000.0
-    keep_1 = cth < 6000.0
-    keep_2 = ctt < 274.5
-    keep_3 = np.logical_and(icing_alt > cth - 1000.0, icing_alt < cth + 1000.0)
-    # keep = keep_0 & keep_1 & keep_2 & keep_3
-    # keep = keep_0 & keep_1 & keep_2
-    # keep = keep_0 & keep_2
-    keep = keep_0
-
-    iint = iint[keep]
-    icing_alt = icing_alt[keep]
-    cth = cth[keep]
-    ctt = ctt[keep]
-    scf = scf[keep]
-    bt = bt[keep]
-    cfrac = cld_frac[keep]
-
-    icing_alt += np.random.normal(scale=100, size=icing_alt.size)
+    num_obs = iint.size
+
+    cld_top_hgt = f['cld_height_acha'][:, :, :].reshape([num_obs, -1])
+    cld_top_tmp = f['cld_temp_acha'][:, :, :].reshape([num_obs, -1])
+    sc_cld_frac = f['supercooled_cloud_fraction'][:, :, :].reshape([num_obs, -1])
+    cld_mask = f['cloud_mask'][:, :, :].reshape([num_obs, -1])
+    cld_opd = f['cld_opd_dcomp'][:, :, :].reshape([num_obs, -1])
+    # cld_opd = f['cld_opd_acha'][:, :, :].reshape([num_obs, -1])
+    cld_reff = f['cld_reff_dcomp'][:, :, :].reshape([num_obs, -1])
+    # cld_reff = f['cld_reff_acha'][:, :, :].reshape([num_obs, -1])
+    cld_emiss = f['cld_emiss_acha'][:, :, :].reshape([num_obs, -1])
+    # cld_phase = f['cloud_phase'][:, :, :].reshape([num_obs, -1])
+    # solzen = f['solar_zenith_angle'][:, :, :].reshape([num_obs, -1])
+
+    cld_mask = np.where(cld_mask > 1, 1, 0)
+    cld_frac = np.sum(cld_mask, axis=1) / 256
+
+    bt_10_4 = f_l1b['temp_10_4um_nom'][:, :, :].reshape([num_obs, -1])
+
+    keep_mask = f['FOV_mask'][:, :, :].reshape([num_obs, -1])
+    keep_mask = keep_mask.astype(np.bool)
+
+
+    # icing_alt += np.random.normal(scale=100, size=icing_alt.size)
 
     iint = np.where(iint == -1, 0, iint)
     iint = np.where(iint != 0, 1, iint)
@@ -938,8 +913,10 @@ def analyze2(filename, filename_l1b):
     f.close()
     f_l1b.close()
 
-    return icing_alt[no_ice], icing_alt[ice], cth[no_ice], cth[ice], ctt[no_ice], ctt[ice], scf[no_ice], scf[ice]
-    # return icing_alt[no_ice], icing_alt[ice], cth[no_ice], cth[ice], ctt[no_ice], ctt[ice], bt[no_ice], bt[ice]
+    return (icing_alt[no_ice], icing_alt[ice],
+            cld_top_hgt[no_ice][keep_mask], cld_top_hgt[ice][keep_mask],
+            cld_top_tmp[no_ice][keep_mask], cld_top_tmp[ice][keep_mask],
+            sc_cld_frac[no_ice][keep_mask], sc_cld_frac[ice][keep_mask])
 
 
 # --------------------------------------------