diff --git a/modules/util/util.py b/modules/util/util.py
index 166744d8199eda8a81ed7a6e134206be4e48f710..733bdfcdb1f54f4825142ef7c35313d67b1c867c 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):
+                         has_time=False, use_nan=False):
     outfile_name = output_dir + 'icing_prediction_'+clvrx_str_time+'.nc'
     rootgrp = Dataset(outfile_name, 'w', format='NETCDF4')
 
@@ -828,18 +828,22 @@ def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
 
     for flvl in flt_lvls:
         probs = probs_dct[flvl]
+        probs = np.where(probs < 0.5, np.nan, probs)
         prob_s.append(probs)
 
         icing_prob_ds = rootgrp.createVariable('icing_probability_level_'+flt_level_ranges_str[flvl], '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)
         if has_time:
             probs = probs.reshape((1, y.shape[0], x.shape[0]))
         icing_prob_ds[:,] = probs
 
     prob_s = np.stack(prob_s, axis=-1)
-    max_prob = np.max(prob_s, axis=2)
+    max_prob = np.nanmax(prob_s, axis=2)
     if has_time:
         max_prob = max_prob.reshape(1, y.shape[0], x.shape[0])
 
@@ -849,7 +853,7 @@ def write_icing_file_nc4(clvrx_str_time, output_dir, preds_dct, probs_dct,
     icing_prob_ds.setncattr('missing', -1.0)
     icing_prob_ds[:,] = max_prob
 
-    max_lvl = np.argmax(prob_s, axis=2)
+    max_lvl = np.nanargmax(prob_s, axis=2)
     if has_time:
         max_lvl = max_lvl.reshape((1, y.shape[0], x.shape[0]))