diff --git a/modules/util/util.py b/modules/util/util.py
index 9bb6db2b034a08daf8f94db021f8355309ebd572..11c74becad8654ac1b531725724b82ea9bcdac0c 100644
--- a/modules/util/util.py
+++ b/modules/util/util.py
@@ -788,7 +788,7 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon
 
 def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
                          x, y, lons, lats, elems, lines, satellite='GOES16', domain='CONUS',
-                         has_time=False, use_nan=False):
+                         has_time=False, use_nan=False, prob_thresh=0.5):
     outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc'
     rootgrp = Dataset(outfile_name, 'w', format='NETCDF4')
 
@@ -839,20 +839,27 @@ def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
             icing_prob_ds.setncattr('missing', np.nan)
         if has_time:
             probs = probs.reshape((1, y.shape[0], x.shape[0]))
+        if use_nan:
+            probs = np.where(probs < prob_thresh, np.nan, probs)
         icing_prob_ds[:,] = probs
 
     prob_s = np.stack(prob_s, axis=-1)
-    max_prob = np.nanmax(prob_s, axis=2)
+    max_prob = np.max(prob_s, axis=2)
+    if use_nan:
+        max_prob = np.where(max_prob < prob_thresh, np.nan, max_prob)
     if has_time:
         max_prob = max_prob.reshape(1, y.shape[0], x.shape[0])
 
     icing_prob_ds = rootgrp.createVariable('max_icing_probability_column', 'f4', var_dim_list)
     icing_prob_ds.setncattr('coordinates', geo_coords)
     icing_prob_ds.setncattr('grid_mapping', 'Projection')
-    icing_prob_ds.setncattr('missing', -1.0)
+    if not use_nan:
+        icing_prob_ds.setncattr('missing', -1.0)
+    else:
+        icing_prob_ds.setncattr('missing', np.nan)
     icing_prob_ds[:,] = max_prob
 
-    max_lvl = np.nanargmax(prob_s, axis=2)
+    max_lvl = np.argmax(prob_s, axis=2)
     if has_time:
         max_lvl = max_lvl.reshape((1, y.shape[0], x.shape[0]))