diff --git a/R/Clim.R b/R/Clim.R
index d4773d5f1783e9b0ce633e910a19201a07a73e60..ce9ba2ed17b283d9a137df32e8b06275e7691f2e 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -190,27 +190,23 @@ Clim <- function(exp, obs, time_dim = 'sdate', dat_dim = c('dataset', 'member'),
   #----------------------------------
   # Per-pair: Remove all sdate if not complete along dat_dim
   if (!is.null(dat_dim)) {
-    pos <- which(names(dim(obs)) %in% dat_dim)
-
-    na_array <- array(0, dim = dim(exp)[-pos])
-    dat_dim_comb <- plyr::alply(expand.grid(lapply(dim(exp)[dat_dim], seq, 1)), 1)
-    for (i_dat in 1:length(dat_dim_comb)) {
-      na_array <- na_array + is.na(Subset(exp, dat_dim, dat_dim_comb[[i_dat]], drop = 'selected'))
-    }
-    dat_dim_comb <- plyr::alply(expand.grid(lapply(dim(obs)[dat_dim], seq, 1)), 1)
-    for (i_dat in 1:length(dat_dim_comb)) {
-      na_array <- na_array + is.na(Subset(obs, dat_dim, dat_dim_comb[[i_dat]], drop = 'selected'))
+    pos <- rep(0, length(dat_dim))
+    for (i in 1:length(dat_dim)) {  #[dat, sdate]
+      ## dat_dim: [dataset, member]
+      pos[i] <- which(names(dim(obs)) == dat_dim[i])
     }
+    outrows_exp <- MeanDims(exp, pos, na.rm = FALSE) + 
+                   MeanDims(obs, pos, na.rm = FALSE)
+    outrows_obs <- outrows_exp
 
-    na_array_2 <- na_array
     for (i in 1:length(pos)) {
-      na_array <- InsertDim(na_array, pos[i], dim(exp)[pos[i]])
-      na_array_2 <- InsertDim(na_array_2, pos[i], dim(obs)[pos[i]])
+      outrows_exp <- InsertDim(outrows_exp, pos[i], dim(exp)[pos[i]])
+      outrows_obs <- InsertDim(outrows_obs, pos[i], dim(obs)[pos[i]])
     }
-    exp[which(na_array != 0)] <- NA
-    obs[which(na_array_2 != 0)] <- NA
-
+    exp[which(is.na(outrows_exp))] <- NA
+    obs[which(is.na(outrows_obs))] <- NA
   }
+
   #-----------------------------------
 
   if (method == 'clim') {
diff --git a/R/MeanDims.R b/R/MeanDims.R
index 9e5dd49fe47ed21008ed0c95da4a2f3163d98acd..56a304dbefb98d359c0458ea44ee1ee450a016b4 100644
--- a/R/MeanDims.R
+++ b/R/MeanDims.R
@@ -10,9 +10,10 @@
 #'  not (FALSE).
 #'@param drop A logical value indicating whether to keep the averaged 
 #'  dimension (FALSE) or drop it (TRUE). The default value is TRUE.
-#'@return An array with the same dimension as parameter 'data' except the 
-#'  'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if 'drop' is
-#'  FALSE, 'dims' will be preserved and the length will be 1.
+#'@return A numeric or an array with the same dimension as parameter 'data'
+#'  except the 'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if
+#' 'drop' is FALSE, 'dims' will be preserved and the length will be 1. If all
+#'  the dimensions are averaged out, a numeric is returned.
 #'
 #'@examples
 #'a <- array(rnorm(24), dim = c(dat = 2, member = 3, time = 4))
@@ -69,7 +70,7 @@ MeanDims <- function(data, dims, na.rm = FALSE, drop = TRUE) {
 
   if (length(dims) == length(dim_data)) {
     if (drop) {
-      data <- as.array(mean(data, na.rm = na.rm))
+      data <- mean(data, na.rm = na.rm)
     } else {
       data <- array(mean(data, na.rm = na.rm), 
                     dim = rep(1, length(dim_data)))
@@ -79,8 +80,8 @@ MeanDims <- function(data, dims, na.rm = FALSE, drop = TRUE) {
     if (is.character(dims)) {
       dims <- which(names(dim_data) %in% dims)
     }
-    pos <- (1:length(dim_data))[-dims]
-    data <- apply(data, pos, mean, na.rm = na.rm)
+    data <- aperm(data, c(dims, (1:length(dim_data))[-dims]))
+    data <- colMeans(data, dims = length(dims), na.rm = na.rm)
 
     # If data is vector
     if (is.null(dim(data))) {
diff --git a/man/MeanDims.Rd b/man/MeanDims.Rd
index 721866b9fbc58225741f89a5866445d600ad1dfb..94369da0bdd255126766d347ba940a5c1d02a963 100644
--- a/man/MeanDims.Rd
+++ b/man/MeanDims.Rd
@@ -19,9 +19,10 @@ not (FALSE).}
 dimension (FALSE) or drop it (TRUE). The default value is TRUE.}
 }
 \value{
-An array with the same dimension as parameter 'data' except the 
- 'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if 'drop' is
- FALSE, 'dims' will be preserved and the length will be 1.
+A numeric or an array with the same dimension as parameter 'data'
+ except the 'dims' dimensions. If 'drop' is TRUE, 'dims' will be removed; if
+'drop' is FALSE, 'dims' will be preserved and the length will be 1. If all
+ the dimensions are averaged out, a numeric is returned.
 }
 \description{
 This function returns the mean of an array along a set of dimensions and 
diff --git a/tests/testthat/test-MeanDims.R b/tests/testthat/test-MeanDims.R
index a431b995422d6cff9050d5f7a2248da0f15451bd..c043c784b8f302ae8c5ab53f5a00a206601ba4fc 100644
--- a/tests/testthat/test-MeanDims.R
+++ b/tests/testthat/test-MeanDims.R
@@ -78,6 +78,14 @@ test_that("2. Output checks: dat1", {
     dim(MeanDims(dat1, dims = c('sdate', 'ftime'))),
     c(dat = 1)
   )
+  expect_equal(
+    dim(MeanDims(dat1, dims = c('sdate', 'ftime', 'dat'))),
+    NULL
+  )
+  expect_equal(
+    dim(MeanDims(dat1, dims = 1:3)),
+    NULL
+  )
   expect_equal(
   MeanDims(dat1, dims = c('sdate'))[1:2],
   c(3, 8)
@@ -96,7 +104,7 @@ test_that("2. Output checks: dat1", {
   )
   expect_equal(
   dim(MeanDims(dat1, dims = 1:3, drop = T)),
-  1
+  NULL
   )
   expect_equal(
   as.vector(drop(MeanDims(dat1, dims = 1:3, drop = F))),
@@ -153,6 +161,10 @@ test_that("5. Output checks: dat4", {
     dim(MeanDims(dat4, dims = 1, drop = F)),
     1
   )
+  expect_equal(
+    dim(MeanDims(dat4, dims = 1, drop = T)),
+    NULL
+  )
 })
 ##############################################