diff --git a/modules/Loading/testing_recipes/recipe_era5land.yml b/modules/Loading/testing_recipes/recipe_era5land.yml
index 55da61e5f970039391319604b3026bb075c7b037..c99906b68c37e65c19c98ab062947d6f2d078e37 100644
--- a/modules/Loading/testing_recipes/recipe_era5land.yml
+++ b/modules/Loading/testing_recipes/recipe_era5land.yml
@@ -33,7 +33,7 @@ Analysis:
     Skill:
       metric: RPS RPSS CRPS CRPSS FRPSS BSS10 BSS90 EnsCorr Corr
     Probabilities:
-      percentiles: [[1/3, 2/3], [1/10, 9/10]] #, [1/4, 2/4, 3/4]]
+      percentiles: [[1/4, 2/4, 3/4]]
     Indicators:
       index: no
   ncores: 1
diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R
index d6b959c8a035ce8b7d09f8074339859f951947dc..7ed2947fc2fb104b1d9d59459847e66728899def 100644
--- a/modules/Skill/Skill.R
+++ b/modules/Skill/Skill.R
@@ -181,11 +181,14 @@ compute_probabilities <- function(data, recipe) {
   } else {
     ncores <- recipe$Analysis$ncores
   }
-  if (is.null(recipe$Analysis$remove_NAs)) {
-    na.rm = T
-  } else {
-    na.rm = recipe$Analysis$remove_NAs
-  }
+  
+  ## TODO: Remove commented lines and include warning if quantile()
+  ## can not accept na.rm = FALSE
+  # if (is.null(recipe$Analysis$remove_NAs)) {
+  #   na.rm = T
+  # } else {
+  #   na.rm = recipe$Analysis$remove_NAs
+  # }
   
   named_probs <- list()
   named_quantiles <- list()
@@ -194,37 +197,35 @@ compute_probabilities <- function(data, recipe) {
 	 "thresholds are provided in the recipe.")
   } else {
     for (element in recipe$Analysis$Workflow$Probabilities$percentiles) {
+      # Parse thresholds in recipe
       thresholds <- sapply(element, function (x) eval(parse(text = x)))
+      # Compute quantiles and probability bins
       probs <- Compute_probs(data$data, thresholds,
 			     ncores = ncores,
-			     na.rm = na.rm)
+			     na.rm = TRUE)
+      # Separate quantiles into individual arrays and name them percentile_xx
       for (i in seq(1:dim(probs$quantiles)['bin'][[1]])) {
-	named_quantiles <- append(named_quantiles,
-				  list(ClimProjDiags::Subset(probs$quantiles,
-							     'bin', i)))
-	names(named_quantiles)[length(named_quantiles)] <- paste0("percentile_",
-					    as.integer(thresholds[i]*100))
+        named_quantiles <- append(named_quantiles,
+        			  list(probs$quantiles[i, , , , , , ,]))
+        names(named_quantiles)[length(named_quantiles)] <- paste0("percentile_",
+        				    as.integer(thresholds[i]*100))
       }
+      # Separate probability bins into individual arrays and name them:
+      # prob_b = prob. below, prob_a = prob. above, prob_xx_to_yy
       for (i in seq(1:dim(probs$probs)['bin'][[1]])) {
-	if (i == 1) {
-	  name_i <- paste0("prob_b", as.integer(thresholds[1]*100))
-	} else if (i == dim(probs$probs)['bin'][[1]]) {
-	  name_i <- paste0("prob_a", as.integer(thresholds[i-1]*100))
-	} else {
-	  name_i <- paste0("prob_", as.integer(thresholds[i-1]*100), "_to_",
-			   as.integer(thresholds[i]*100))
-	}
-	named_probs <- append(named_probs,
-			      list(ClimProjDiags::Subset(probs$probs,
-							 'bin', i)))
-	names(named_probs)[length(named_probs)] <- name_i
+        if (i == 1) {
+          name_i <- paste0("prob_b", as.integer(thresholds[1]*100))
+        } else if (i == dim(probs$probs)['bin'][[1]]) {
+          name_i <- paste0("prob_a", as.integer(thresholds[i-1]*100))
+        } else {
+          name_i <- paste0("prob_", as.integer(thresholds[i-1]*100), "_to_",
+        		   as.integer(thresholds[i]*100))
+        }
+        named_probs <- append(named_probs, list(probs$probs[i, , , , , , , ,]))
+        names(named_probs)[length(named_probs)] <- name_i
       }
-      # remove(probs)
     }
   }
-  named_quantiles <- lapply(named_quantiles, function(x) {.drop_dims(x)})
-  named_probs <- lapply(named_probs, function(x) {.drop_dims(x)})
-
   print("##### PERCENTILES AND PROBABILITY CATEGORIES COMPUTED #####")
   return(list(probs=named_probs, percentiles=named_quantiles))
 }
diff --git a/modules/Skill/s2s.probs.R b/modules/Skill/s2s.probs.R
index a8cab57c16618a163a8526331ebcabdaa78b31d4..889330453d15bb3e90e7f14de5eb5b456697643d 100644
--- a/modules/Skill/s2s.probs.R
+++ b/modules/Skill/s2s.probs.R
@@ -23,7 +23,7 @@ Compute_probs <- function(data, thresholds,
       if (any(!is.na(x))){
         colMeans(convert2prob(as.vector(x), threshold = as.vector(t)))
       } else {
-        c(NA, NA, NA)
+	rep(NA, dim(t)[['bin']] + 1) # vector with as many NAs as probability bins.
       }
     }
   }