diff --git a/modules/deeplearning/cloud_opd_fcn_abi.py b/modules/deeplearning/cloud_opd_fcn_abi.py
index 3c57441b965305ee8e74699f5bd2d51e44bfec34..02baf1909b03639c245ee1d786a1d7a09b881aed 100644
--- a/modules/deeplearning/cloud_opd_fcn_abi.py
+++ b/modules/deeplearning/cloud_opd_fcn_abi.py
@@ -177,22 +177,29 @@ def get_min_max_std(grd_k):
 
 
 def get_cldy_frac_opd(cld_prob, opd):
-    cld_prob = np.where(np.isnan(cld_prob), 0, cld_prob)
+    cld_prob = np.where(np.isnan(cld_prob), 0.0, cld_prob)
+    opd = np.where(np.isnan(opd), 0.0, opd)
     cld = np.where(cld_prob < 0.5, 0, 1)
-    opd[cld == 0] = 0.0
 
-    s = cld[:, 0::4, 0::4] + cld[:, 1::4, 0::4] + cld[:, 2::4, 0::4] + cld[:, 3::4, 0::4] + \
-        cld[:, 0::4, 1::4] + cld[:, 1::4, 1::4] + cld[:, 2::4, 1::4] + cld[:, 3::4, 1::4] + \
-        cld[:, 0::4, 2::4] + cld[:, 1::4, 2::4] + cld[:, 2::4, 2::4] + cld[:, 3::4, 2::4] + \
-        cld[:, 0::4, 3::4] + cld[:, 1::4, 3::4] + cld[:, 2::4, 3::4] + cld[:, 3::4, 3::4]
+    cnt_cld = cld[:, 0::4, 0::4] + cld[:, 1::4, 0::4] + cld[:, 2::4, 0::4] + cld[:, 3::4, 0::4] + \
+              cld[:, 0::4, 1::4] + cld[:, 1::4, 1::4] + cld[:, 2::4, 1::4] + cld[:, 3::4, 1::4] + \
+              cld[:, 0::4, 2::4] + cld[:, 1::4, 2::4] + cld[:, 2::4, 2::4] + cld[:, 3::4, 2::4] + \
+              cld[:, 0::4, 3::4] + cld[:, 1::4, 3::4] + cld[:, 2::4, 3::4] + cld[:, 3::4, 3::4]
 
-    cldy_opd = np.sum([opd[:, 0::4, 0::4], opd[:, 1::4, 0::4], opd[:, 2::4, 0::4], opd[:, 3::4, 0::4],
-                       opd[:, 0::4, 1::4], opd[:, 1::4, 1::4], opd[:, 2::4, 1::4], opd[:, 3::4, 1::4],
-                       opd[:, 0::4, 2::4], opd[:, 1::4, 2::4], opd[:, 2::4, 2::4], opd[:, 3::4, 2::4],
-                       opd[:, 0::4, 3::4], opd[:, 1::4, 3::4], opd[:, 2::4, 3::4], opd[:, 3::4, 3::4]], axis=0)
+    opd_sum = np.sum([opd[:, 0::4, 0::4], opd[:, 1::4, 0::4], opd[:, 2::4, 0::4], opd[:, 3::4, 0::4],
+                      opd[:, 0::4, 1::4], opd[:, 1::4, 1::4], opd[:, 2::4, 1::4], opd[:, 3::4, 1::4],
+                      opd[:, 0::4, 2::4], opd[:, 1::4, 2::4], opd[:, 2::4, 2::4], opd[:, 3::4, 2::4],
+                      opd[:, 0::4, 3::4], opd[:, 1::4, 3::4], opd[:, 2::4, 3::4], opd[:, 3::4, 3::4]], axis=0)
 
-    s = np.where(s == 0, 16, s)
-    cldy_opd /= s
+    opd[cld == 0] = 0.0
+    cld_opd_sum = np.sum([opd[:, 0::4, 0::4], opd[:, 1::4, 0::4], opd[:, 2::4, 0::4], opd[:, 3::4, 0::4],
+                          opd[:, 0::4, 1::4], opd[:, 1::4, 1::4], opd[:, 2::4, 1::4], opd[:, 3::4, 1::4],
+                          opd[:, 0::4, 2::4], opd[:, 1::4, 2::4], opd[:, 2::4, 2::4], opd[:, 3::4, 2::4],
+                          opd[:, 0::4, 3::4], opd[:, 1::4, 3::4], opd[:, 2::4, 3::4], opd[:, 3::4, 3::4]], axis=0)
+
+    cldy_opd = np.zeros(cnt_cld.shape, dtype=opd.dtype)
+    cldy_opd[cnt_cld == 0] = opd_sum[cnt_cld == 0] / 16
+    cldy_opd[cnt_cld != 0] = cld_opd_sum[cnt_cld != 0] / cnt_cld[cnt_cld != 0]
 
     return cldy_opd