diff --git a/R/MeanDims.R b/R/MeanDims.R
index 2da3144c55259551d7f7cb30fef36f26776d0e17..4b22d518cc5a7f1dab6320de16c7186115c8c6c8 100644
--- a/R/MeanDims.R
+++ b/R/MeanDims.R
@@ -3,26 +3,28 @@
 #'This function returns the mean of an array along a set of dimensions and 
 #'preserves the dimension names if it has.
 #'
-#'@details It is recommended to use \code{'apply(x, dim, mean)'} to improve the
-#'  efficiency when the dimension to be averaged is only one. 
-#'
 #'@param data An array to be averaged.
 #'@param dims A vector of numeric or charactor string, indicating along which 
 #'  dimensions to average.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
-#'  not (FALSE). The default value is FALSE.
-#'
+#'  not (FALSE).
+#'@param ncores A integer indicating the number of cores to use in parallel computation.
 #'@return An array with the same dimension as parameter 'data' except the 'dims' 
 #'  dimensions. 
 #'  removed.
 #'
+#'@keywords datagen
+#'@author History:\cr
+#'0.1  -  2011-04  (V. Guemas, \email{vguemas@@ic3.cat})  -  Original code\cr
+#'1.0  -  2013-09  (N. Manubens, \email{nicolau.manubens@@ic3.cat})  -  Formatting to R CRAN\cr
+#'1.1  -  2015-03  (N. Manubens, \email{nicolau.manubens@@ic3.cat})  -  Improved memory usage
+#'3.0  -  2020-01  (A. Ho, \email{an.ho@bsc.es})  -  Preserve dimension names
 #'@examples
-#'a <- array(rnorm(24), dim = c(a = 2, b= 3, c = 4))
-#'print(dim(MeanDims(a, 2)))
-#'print(dim(MeanDims(a, c(2, 3))))
-#'print(dim(MeanDims(a, c('a', 'b'))))
+#'a <- array(rnorm(24), dim = c(2, 3, 4))
+#'MeanDims(a, 2)
+#'MeanDims(a, c(2, 3))
 #'@export
-MeanDims <- function(data, dims, na.rm = FALSE) {
+MeanDims <- function(data, dims, na.rm = FALSE, ncores = NULL) {
 
   # Check inputs 
   ## data
@@ -63,28 +65,18 @@ MeanDims <- function(data, dims, na.rm = FALSE) {
 
   ###############################
   # Calculate MeanDims
-
-  ## Change character dims into indices
-  if (is.character(dims)) {
-    tmp <- rep(0, length(dims))
-    for (i in 1:length(dims)) {
-      tmp[i] <- which(names(dim(data)) == dims[i])
+  if (length(dims) == length(dim(data)) || length(dim(data)) == 1) {
+    res <- mean(data, na.rm = na.rm)
+  } else if (length(dims) == 1) {
+    if (is.character(dims)) {
+      dims <- which(names(dim(data)) == dims)
     }
-  dims <- tmp
-  }
-
-  if (length(dim(data)) == 1) {
-    res <- mean(data, na.rm = na.rm) 
+    pos <- (1:length(dim(data)))[-dims]
+    res <- apply(data, pos, mean, na.rm = na.rm)
   } else {
-
-    margins <- setdiff(c(1:length(dim(data))), dims)
-    res <- as.array(apply(data, margins, mean, na.rm = na.rm))
-    if (!is.null(names(dim(data))[margins]) & is.null(names(dim(res)))) {
-      names(dim(res)) <- names(dim(data))[margins]
-    }
+    res <- Apply(list(data), target_dims = list(dims), fun = mean,
+                 na.rm = na.rm, ncores = ncores)$output1
   }
-
   return(res)
-
 }
 
diff --git a/man/MeanDims.Rd b/man/MeanDims.Rd
index 2e6022f7b5bdaf41d9e37606a4b0838b04f236a8..f2000239679377d0cd131b9fe6fa8f1a696976f6 100644
--- a/man/MeanDims.Rd
+++ b/man/MeanDims.Rd
@@ -4,7 +4,7 @@
 \alias{MeanDims}
 \title{Average an array along multiple dimensions}
 \usage{
-MeanDims(data, dims, na.rm = FALSE)
+MeanDims(data, dims, na.rm = TRUE, ncores = NULL)
 }
 \arguments{
 \item{data}{An array to be averaged.}
@@ -13,7 +13,9 @@ MeanDims(data, dims, na.rm = FALSE)
 dimensions to average.}
 
 \item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
-not (FALSE). The default value is FALSE.}
+not (FALSE).}
+
+\item{ncores}{A integer indicating the number of cores to use in parallel computation.}
 }
 \value{
 An array with the same dimension as parameter 'data' except the 'dims' 
@@ -24,14 +26,17 @@ An array with the same dimension as parameter 'data' except the 'dims'
 This function returns the mean of an array along a set of dimensions and 
 preserves the dimension names if it has.
 }
-\details{
-It is recommended to use \code{'apply(x, dim, mean)'} to improve the
- efficiency when the dimension to be averaged is only one.
-}
 \examples{
-a <- array(rnorm(24), dim = c(a = 2, b= 3, c = 4))
-print(dim(MeanDims(a, 2)))
-print(dim(MeanDims(a, c(2, 3))))
-print(dim(MeanDims(a, c('a', 'b'))))
+a <- array(rnorm(24), dim = c(2, 3, 4))
+MeanDims(a, 2)
+MeanDims(a, c(2, 3))
+}
+\author{
+History:\cr
+0.1  -  2011-04  (V. Guemas, \email{vguemas@ic3.cat})  -  Original code\cr
+1.0  -  2013-09  (N. Manubens, \email{nicolau.manubens@ic3.cat})  -  Formatting to R CRAN\cr
+1.1  -  2015-03  (N. Manubens, \email{nicolau.manubens@ic3.cat})  -  Improved memory usage
+3.0  -  2020-01  (A. Ho, \email{an.ho@bsc.es})  -  Preserve dimension names
 }
+\keyword{datagen}