diff --git a/modules/Skill/compute_probs.R b/modules/Skill/compute_probs.R
index a662df1483569c2bc3bd38d71640473bd16484fc..1c17b35872597c8215337b436e036b0026edf7f8 100644
--- a/modules/Skill/compute_probs.R
+++ b/modules/Skill/compute_probs.R
@@ -4,28 +4,32 @@ compute_probs <- function(data, quantiles,
                           probs_dims=list('ensemble', 'bin'),
                           split_factor=1, na.rm=FALSE) {
 
-  if (na.rm == FALSE) {
-    c2p <- function(x, t) {
-      # If the array contains any NA values, return NA
-      if (any(is.na(x))) {
-	rep(NA, dim(t)[['bin']] + 1)
+  # Define na.rm behavior
+  if (na.rm) {
+    .c2p <- function(x, t) {
+      if (any(!is.na(x))) {
+	# If the array contains any non-NA values, call convert2prob
+        colMeans(convert2prob(as.vector(x), threshold = as.vector(t)))
       } else {
-       colMeans(convert2prob(as.vector(x), threshold = as.vector(t)))
+	# If the array only contains NAs, return NA vector
+	rep(NA, dim(t)[['bin']] + 1) # vector with as many NAs as prob bins.
       }
     }
   } else {
-    c2p <- function(x, t) {
-      if (any(!is.na(x))) { # If the array contains some non-NA values
+    .c2p <- function(x, t) {
+      if (any(is.na(x))) {
+	# If the array contains any NA values, return NA vector
+	rep(NA, dim(t)[['bin']] + 1)
+      } else {
+	# If there are no NAs, call convert2prob
         colMeans(convert2prob(as.vector(x), threshold = as.vector(t)))
-      } else { # If the array contains NAs only
-	rep(NA, dim(t)[['bin']] + 1) # vector with as many NAs as prob bins.
       }
     }
   }
   
   probs <- Apply(data = list(x = data, t = quantiles),
                  target_dims = probs_dims,
-                 c2p,
+                 .c2p,
                  output_dims = "bin",
                  split_factor = split_factor,
                  ncores = ncores)[[1]]
diff --git a/modules/Skill/compute_quants.R b/modules/Skill/compute_quants.R
index 60ad981fb377c928641a43996566cef1e8b4bbdf..8c89e87e4bf41d2a77a0ddb1a3196a8f9643cdbd 100644
--- a/modules/Skill/compute_quants.R
+++ b/modules/Skill/compute_quants.R
@@ -5,23 +5,26 @@ compute_quants <- function(data, thresholds,
                            probs_dims=list('ensemble', 'bin'),
                            split_factor=1, na.rm=FALSE) {
 
-  if (na.rm == FALSE) {
-    get_quantiles <- function(x, t) {
+  # Define na.rm behavior
+  if (na.rm) {
+    .get_quantiles <- function(x, t) {
+      quantile(as.vector(x), t, na.rm = TRUE)
+    }
+  } else {
+    .get_quantiles <- function(x, t) {
       if (any(is.na(x))) {
+	# If the array contains any NA values, return NA vector
 	rep(NA, length(t))
       } else {
+        # If there are no NAs, call quantile()
 	quantile(as.vector(x), t, na.rm = FALSE)
       }
     }
-  } else {
-    get_quantiles <- function(x, t) {
-      quantile(as.vector(x), t, na.rm = TRUE)
-    }
   }
                             
   quantiles <- Apply(data,
                      target_dims = quantile_dims,
-		     function(x, t) {get_quantiles(as.vector(x), thresholds)},
+		     function(x, t) {.get_quantiles(as.vector(x), thresholds)},
                      output_dims = "bin",
                      ncores = ncores, 
                      split_factor = split_factor)[[1]]