diff --git a/modules/icing/pirep_goes.py b/modules/icing/pirep_goes.py
index 39e5a34ecd4a74cf0bd448e7b663c611d209b5aa..8c08cd616ef089c0cfb59ba9dc6cd30ca9a515c5 100644
--- a/modules/icing/pirep_goes.py
+++ b/modules/icing/pirep_goes.py
@@ -2200,6 +2200,7 @@ def run_icing_predict(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=h
     for fidx, fname in enumerate(clvrx_files):
         h5f = h5py.File(fname, 'r')
         dto = clvrx_ds.get_datetime(fname)
+        ts = dto.timestamp()
         clvrx_str_time = dto.strftime('%Y-%m-%d_%H:%M')
 
         data_dct, ll, cc = make_for_full_domain_predict(h5f, name_list=train_params, satellite=satellite, domain=domain)
@@ -2255,7 +2256,8 @@ def run_icing_predict(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=h
 
             preds_day_dct, probs_day_dct = run_evaluate_static(day_grd_dct, num_day_tiles, day_night='DAY',
                                                                ckpt_dir_s_path=day_model_path, prob_thresh=prob_thresh,
-                                                               use_flight_altitude=use_flight_altitude)
+                                                               use_flight_altitude=use_flight_altitude,
+                                                               flight_levels=flight_levels)
             day_idxs = np.array(day_idxs)
             for flvl in flight_levels:
                 day_preds = preds_day_dct[flvl]
@@ -2277,7 +2279,8 @@ def run_icing_predict(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=h
 
             preds_nght_dct, probs_nght_dct = run_evaluate_static(nght_grd_dct, num_nght_tiles, day_night='NIGHT',
                                                                  ckpt_dir_s_path=night_model_path, prob_thresh=prob_thresh,
-                                                                 use_flight_altitude=use_flight_altitude)
+                                                                 use_flight_altitude=use_flight_altitude,
+                                                                 flight_levels=flight_levels)
             nght_idxs = np.array(nght_idxs)
             for flvl in flight_levels:
                 nght_preds = preds_nght_dct[flvl]
@@ -2298,3 +2301,162 @@ def run_icing_predict(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=h
 
         print('Done: ', clvrx_str_time)
         h5f.close()
+
+
+def run_icing_predict_image(clvrx_dir='/Users/tomrink/data/clavrx/RadC/', output_dir=homedir,
+                            day_model_path=model_path_day, night_model_path=model_path_night,
+                            prob_thresh=0.5, satellite='GOES16', domain='CONUS', day_night='AUTO',
+                            l1b_andor_l2='BOTH', use_flight_altitude=True,
+                            extent=[-105, -70, 15, 50],
+                            pirep_file='/Users/tomrink/data/pirep/pireps_202109200000_202109232359.csv',
+                            obs_lons=None, obs_lats=None, obs_times=None, obs_alt=None, flight_level=None):
+
+    flight_levels = [0, 1, 2, 3, 4]
+
+    if pirep_file is not None:
+        ice_dict, no_ice_dict, neg_ice_dict = setup(pirep_file)
+
+    alt_lo, alt_hi = 0.0, 15000.0
+    if flight_level is not None:
+        alt_lo, alt_hi = flt_level_ranges[flight_level]
+
+    day_train_params = get_training_parameters(day_night='DAY', l1b_andor_l2=l1b_andor_l2)
+    nght_train_params = get_training_parameters(day_night='NIGHT', l1b_andor_l2=l1b_andor_l2)
+
+    if day_night == 'AUTO':
+        train_params = list(set(day_train_params + nght_train_params))
+    elif day_night == 'DAY':
+        train_params = day_train_params
+    elif day_night == 'NIGHT':
+        train_params = nght_train_params
+
+    if satellite == 'H08':
+        clvrx_ds = CLAVRx_H08(clvrx_dir)
+    else:
+        clvrx_ds = CLAVRx(clvrx_dir)
+    clvrx_files = clvrx_ds.flist
+
+    for fidx, fname in enumerate(clvrx_files):
+        h5f = h5py.File(fname, 'r')
+        dto = clvrx_ds.get_datetime(fname)
+        ts = dto.timestamp()
+        clvrx_str_time = dto.strftime('%Y-%m-%d_%H:%M')
+
+        data_dct, ll, cc = make_for_full_domain_predict(h5f, name_list=train_params, satellite=satellite, domain=domain)
+
+        if fidx == 0:
+            num_elems = len(cc)
+            num_lines = len(ll)
+            nav = get_navigation(satellite, domain)
+            lons_2d, lats_2d, x_rad, y_rad = get_lon_lat_2d_mesh(nav, ll, cc, offset=8)
+
+        ancil_data_dct, _, _ = make_for_full_domain_predict(h5f, name_list=
+                            ['solar_zenith_angle', 'sensor_zenith_angle', 'cld_height_acha', 'cld_geo_thick'],
+                            satellite=satellite, domain=domain)
+
+        satzen = ancil_data_dct['sensor_zenith_angle']
+        solzen = ancil_data_dct['solar_zenith_angle']
+        day_idxs = []
+        nght_idxs = []
+        for j in range(num_lines):
+            for i in range(num_elems):
+                k = i + j*num_elems
+                if not check_oblique(satzen[k]):
+                    continue
+                if is_day(solzen[k]):
+                    day_idxs.append(k)
+                else:
+                    nght_idxs.append(k)
+
+        num_tiles = num_lines * num_elems
+        num_day_tiles = len(day_idxs)
+        num_nght_tiles = len(nght_idxs)
+
+        # initialize output arrays
+        probs_2d_dct = {flvl: None for flvl in flight_levels}
+        preds_2d_dct = {flvl: None for flvl in flight_levels}
+        for flvl in flight_levels:
+            fd_preds = np.zeros(num_lines * num_elems, dtype=np.int8)
+            fd_preds[:] = -1
+            fd_probs = np.zeros(num_lines * num_elems, dtype=np.float32)
+            fd_probs[:] = -1.0
+            preds_2d_dct[flvl] = fd_preds
+            probs_2d_dct[flvl] = fd_probs
+
+        if (day_night == 'AUTO' or day_night == 'DAY') and num_day_tiles > 0:
+
+            day_data_dct = {name: [] for name in day_train_params}
+            for name in day_train_params:
+                for k in day_idxs:
+                    day_data_dct[name].append(data_dct[name][k])
+            day_grd_dct = {name: None for name in day_train_params}
+            for ds_name in day_train_params:
+                day_grd_dct[ds_name] = np.stack(day_data_dct[ds_name])
+
+            preds_day_dct, probs_day_dct = run_evaluate_static(day_grd_dct, num_day_tiles, day_night='DAY',
+                                                               ckpt_dir_s_path=day_model_path, prob_thresh=prob_thresh,
+                                                               use_flight_altitude=use_flight_altitude,
+                                                               flight_levels=flight_levels)
+            day_idxs = np.array(day_idxs)
+            for flvl in flight_levels:
+                day_preds = preds_day_dct[flvl]
+                day_probs = probs_day_dct[flvl]
+                fd_preds = preds_2d_dct[flvl]
+                fd_probs = probs_2d_dct[flvl]
+                fd_preds[day_idxs] = day_preds[:]
+                fd_probs[day_idxs] = day_probs[:]
+
+        if (day_night == 'AUTO' or day_night == 'NIGHT') and num_nght_tiles > 0:
+
+            nght_data_dct = {name: [] for name in nght_train_params}
+            for name in nght_train_params:
+                for k in nght_idxs:
+                    nght_data_dct[name].append(data_dct[name][k])
+            nght_grd_dct = {name: None for name in nght_train_params}
+            for ds_name in nght_train_params:
+                nght_grd_dct[ds_name] = np.stack(nght_data_dct[ds_name])
+
+            preds_nght_dct, probs_nght_dct = run_evaluate_static(nght_grd_dct, num_nght_tiles, day_night='NIGHT',
+                                                                 ckpt_dir_s_path=night_model_path, prob_thresh=prob_thresh,
+                                                                 use_flight_altitude=use_flight_altitude,
+                                                                 flight_levels=flight_levels)
+            nght_idxs = np.array(nght_idxs)
+            for flvl in flight_levels:
+                nght_preds = preds_nght_dct[flvl]
+                nght_probs = probs_nght_dct[flvl]
+                fd_preds = preds_2d_dct[flvl]
+                fd_probs = probs_2d_dct[flvl]
+                fd_preds[nght_idxs] = nght_preds[:]
+                fd_probs[nght_idxs] = nght_probs[:]
+
+        for flvl in flight_levels:
+            fd_preds = preds_2d_dct[flvl]
+            fd_probs = probs_2d_dct[flvl]
+            preds_2d_dct[flvl] = fd_preds.reshape((num_lines, num_elems))
+            probs_2d_dct[flvl] = fd_probs.reshape((num_lines, num_elems))
+
+        # write_icing_file_nc4(clvrx_str_time, output_dir, preds_2d_dct, probs_2d_dct,
+        #                      x_rad, y_rad, lons_2d, lats_2d, cc, ll, satellite=satellite, domain=domain)
+
+        dto, _ = get_time_tuple_utc(ts)
+        dto_0 = dto - datetime.timedelta(minutes=30)
+        dto_1 = dto + datetime.timedelta(minutes=30)
+        ts_0 = dto_0.timestamp()
+        ts_1 = dto_1.timestamp()
+
+        if pirep_file is not None:
+            _, keep_lons, keep_lats, _ = time_filter_3(ice_dict, ts_0, ts_1, alt_lo, alt_hi)
+        elif obs_times is not None:
+            keep = np.logical_and(obs_times >= ts_0, obs_times < ts_1)
+            keep = np.where(keep, np.logical_and(obs_alt >= alt_lo, obs_alt < alt_hi), False)
+            keep_lons = obs_lons[keep]
+            keep_lats = obs_lats[keep]
+        else:
+            keep_lons = None
+            keep_lats = None
+
+        make_icing_image(None, probs_2d_dct[0], None, None, clvrx_str_time, satellite, domain,
+                         ice_lons_vld=keep_lons, ice_lats_vld=keep_lats, extent=extent)
+
+        print('Done: ', clvrx_str_time)
+        h5f.close()
\ No newline at end of file