diff --git a/R/Ano.R b/R/Ano.R
index e0a69db232230da28d8dde61fc1d92d07f123de5..c4c70a330bf00074632ec18f1b6ea243e621ced8 100644
--- a/R/Ano.R
+++ b/R/Ano.R
@@ -21,10 +21,10 @@
 #'clim <- Clim(sampleData$mod, sampleData$obs)
 #'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 #'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-#'\donttest{
+#'\dontrun{
 #'PlotAno(ano_exp, ano_obs, startDates, 
 #'        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
-#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+#'        legends = 'ERSST', biglab = FALSE)
 #'}
 #'@import multiApply
 #'@export
diff --git a/R/Ano_CrossValid.R b/R/Ano_CrossValid.R
index 22e710a1a3ded8b428eb7493a1dd0c0aec489d18..99205020b7412afdf90513ad18dedf0b9302320d 100644
--- a/R/Ano_CrossValid.R
+++ b/R/Ano_CrossValid.R
@@ -47,7 +47,7 @@
 #'\dontrun{
 #'PlotAno(anomalies$exp, anomalies$obs, startDates, 
 #'        toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'), 
-#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano_crossvalid.eps')
+#'        legends = 'ERSST', biglab = FALSE)
 #'}
 #'@import multiApply
 #'@importFrom ClimProjDiags Subset
diff --git a/R/Clim.R b/R/Clim.R
index 21f97b67ebcd586efb284f6654523b51f53eadd7..d01cf5d91ec219dd7b9124abccabb392d1f35215 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -56,11 +56,11 @@
 #'example(Load)
 #'clim <- Clim(sampleData$mod, sampleData$obs)
 #'clim2 <- Clim(sampleData$mod, sampleData$obs, method = 'kharin', memb = FALSE)
-#'\donttest{
+#'\dontrun{
 #'PlotClim(clim$clim_exp, clim$clim_obs, 
 #'         toptitle = paste('sea surface temperature climatologies'), 
 #'         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-#'         listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
+#'         listobs = c('ERSST'), biglab = FALSE)
 #'}
 #'@importFrom abind adrop
 #'@importFrom ClimProjDiags Subset
diff --git a/R/MeanDims.R b/R/MeanDims.R
index 8d5d77b045876324d3c81c27311ee1809c0b4504..9e5dd49fe47ed21008ed0c95da4a2f3163d98acd 100644
--- a/R/MeanDims.R
+++ b/R/MeanDims.R
@@ -8,17 +8,21 @@
 #'  dimensions to average.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
 #'  not (FALSE).
-#'@return An array with the same dimension as parameter 'data' except the 'dims' 
-#'  dimensions. 
-#'  removed.
+#'@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.
 #'
 #'@examples
-#'a <- array(rnorm(24), dim = c(2, 3, 4))
-#'MeanDims(a, 2)
-#'MeanDims(a, c(2, 3))
+#'a <- array(rnorm(24), dim = c(dat = 2, member = 3, time = 4))
+#'ens_mean <- MeanDims(a, 'member')
+#'dim(ens_mean)
+#'ens_time_mean <- MeanDims(a, c(2, 3), drop = FALSE)
+#'dim(ens_time_mean)
 #'@import multiApply
 #'@export
-MeanDims <- function(data, dims, na.rm = FALSE) {
+MeanDims <- function(data, dims, na.rm = FALSE, drop = TRUE) {
 
   # Check inputs 
   ## data
@@ -54,23 +58,42 @@ MeanDims <- function(data, dims, na.rm = FALSE) {
   if (!is.logical(na.rm) | length(na.rm) > 1) {
     stop("Parameter 'na.rm' must be one logical value.")
   }
+  ## drop
+  if (!is.logical(drop) | length(drop) > 1) {
+    stop("Parameter 'drop' must be one logical value.")
+  }
+
   ###############################
   # Calculate MeanDims
   dim_data <- dim(data)
 
   if (length(dims) == length(dim_data)) {
-    data <- mean(data, na.rm = na.rm)
+    if (drop) {
+      data <- as.array(mean(data, na.rm = na.rm))
+    } else {
+      data <- array(mean(data, na.rm = na.rm), 
+                    dim = rep(1, length(dim_data)))
+      names(dim(data)) <- names(dim_data)
+    }
   } else {
     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)
-  
+
+    # If data is vector
     if (is.null(dim(data))) {
       data <- array(data, dim = dim_data[-dims])
     }
+    if (!drop) {
+      restore_dim <- as.array(rep(1, length(dims)))
+      names(restore_dim) <- names(dim_data)[dims]
+      data <- array(data, dim = c(dim(data), restore_dim))
+      data <- Reorder(data, names(dim_data))
+    }
   }
+
   return(data)
 }
 
diff --git a/R/RatioSDRMS.R b/R/RatioSDRMS.R
index b74fbb494798eba739d07396053275e32b504486..85313030dfd0bedaeb7b78c44423dbfb559c45dc 100644
--- a/R/RatioSDRMS.R
+++ b/R/RatioSDRMS.R
@@ -44,11 +44,10 @@
 #'rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
 #'rsdrms_plot[, , 2, ] <- rsdrms$ratio
 #'rsdrms_plot[, , 4, ] <- rsdrms$p.val
-#'\donttest{
+#'\dontrun{
 #'PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
 #'            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
-#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE, 
-#'            fileout = 'tos_rsdrms.eps')
+#'            listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
 #'}
 #'
 #'@import multiApply
diff --git a/R/Smoothing.R b/R/Smoothing.R
index d5fd2a5eac95d11d7217b05d2c5b795f8792d71d..1b31e6594858262fbea5353392e04a7435e2aec8 100644
--- a/R/Smoothing.R
+++ b/R/Smoothing.R
@@ -25,10 +25,9 @@
 #'smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
 #'smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
 #'smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
-#'  \donttest{
+#'  \dontrun{
 #'PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
-#'        toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
-#'        fileout = "tos_smoothed_ano.png")
+#'        toptitle = "Smoothed Mediterranean mean SST", ytitle = "K")
 #'  }
 #'@import plyr multiApply
 #'@export
diff --git a/R/Spread.R b/R/Spread.R
index 4b3bc6b1167534dc9004a91d99da686851ed4611..9f1624a3f34934b4de61e1523044c20ab7db28af 100644
--- a/R/Spread.R
+++ b/R/Spread.R
@@ -54,27 +54,27 @@
 #'                                                   name = 'member')
 #'spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
 #'
-#'\donttest{
+#'\dontrun{
 #'PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Inter-Quartile Range between ensemble members",
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_iqr.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Maximum minus minimum of the members", 
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_maxmin.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Standard deviation of the members", 
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_sd.png')
+#'            hlines = c(0))
 #'PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
 #'            toptitle = "Median Absolute Deviation of the members",
 #'            ytitle = "K", monini = 11, limits = NULL, 
 #'            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-#'            hlines = c(0), fileout = 'tos_mad.png')
+#'            hlines = c(0))
 #'}
 #'
 #'@import multiApply
diff --git a/man/Ano.Rd b/man/Ano.Rd
index d15ffd14bdbfbb52a3975d8146267edb20eca0b9..e30467e0f3070aacae112aa9086322f3320f2903 100644
--- a/man/Ano.Rd
+++ b/man/Ano.Rd
@@ -32,9 +32,9 @@ example(Load)
 clim <- Clim(sampleData$mod, sampleData$obs)
 ano_exp <- Ano(sampleData$mod, clim$clim_exp)
 ano_obs <- Ano(sampleData$obs, clim$clim_obs)
-\donttest{
+\dontrun{
 PlotAno(ano_exp, ano_obs, startDates, 
        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
-       legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+       legends = 'ERSST', biglab = FALSE)
 }
 }
diff --git a/man/Ano_CrossValid.Rd b/man/Ano_CrossValid.Rd
index 1e91528335c2e68b812594311f1446a8d7eabeb3..d2234a1623901efa767aba83bbd4c5eb532940a7 100644
--- a/man/Ano_CrossValid.Rd
+++ b/man/Ano_CrossValid.Rd
@@ -69,6 +69,6 @@ anomalies <- Ano_CrossValid(sampleData$mod, sampleData$obs)
 \dontrun{
 PlotAno(anomalies$exp, anomalies$obs, startDates, 
        toptitle = paste('anomalies'), ytitle = c('K', 'K', 'K'), 
-       legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano_crossvalid.eps')
+       legends = 'ERSST', biglab = FALSE)
 }
 }
diff --git a/man/Clim.Rd b/man/Clim.Rd
index 78559bdbefc9c5253d510f2c7eee1fbc743324aa..50ec0d9c90d7e2dc1b0486de820e657ee84aa325 100644
--- a/man/Clim.Rd
+++ b/man/Clim.Rd
@@ -84,10 +84,10 @@ the 'exp' and 'obs' are excluded when computing the climatologies.
 example(Load)
 clim <- Clim(sampleData$mod, sampleData$obs)
 clim2 <- Clim(sampleData$mod, sampleData$obs, method = 'kharin', memb = FALSE)
-\donttest{
+\dontrun{
 PlotClim(clim$clim_exp, clim$clim_obs, 
         toptitle = paste('sea surface temperature climatologies'), 
         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-        listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
+        listobs = c('ERSST'), biglab = FALSE)
 }
 }
diff --git a/man/MeanDims.Rd b/man/MeanDims.Rd
index f70b78b391f1205a30614735a8af875a64c9b608..721866b9fbc58225741f89a5866445d600ad1dfb 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 = FALSE, drop = TRUE)
 }
 \arguments{
 \item{data}{An array to be averaged.}
@@ -14,18 +14,23 @@ dimensions to average.}
 
 \item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
 not (FALSE).}
+
+\item{drop}{A logical value indicating whether to keep the averaged 
+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. 
- removed.
+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.
 }
 \description{
 This function returns the mean of an array along a set of dimensions and 
 preserves the dimension names if it has.
 }
 \examples{
-a <- array(rnorm(24), dim = c(2, 3, 4))
-MeanDims(a, 2)
-MeanDims(a, c(2, 3))
+a <- array(rnorm(24), dim = c(dat = 2, member = 3, time = 4))
+ens_mean <- MeanDims(a, 'member')
+dim(ens_mean)
+ens_time_mean <- MeanDims(a, c(2, 3), drop = FALSE)
+dim(ens_time_mean)
 }
diff --git a/man/RatioSDRMS.Rd b/man/RatioSDRMS.Rd
index 7dbd68283599edcfa599768a1946b00ce89fdea1..f1f6f3ddea92b777910d2fdc0176e9fd02a5f67a 100644
--- a/man/RatioSDRMS.Rd
+++ b/man/RatioSDRMS.Rd
@@ -67,11 +67,10 @@ rsdrms <- RatioSDRMS(sampleData$mod, sampleData$obs)
 rsdrms_plot <- array(dim = c(dim(rsdrms$ratio)[1:2], 4, dim(rsdrms$ratio)[3]))
 rsdrms_plot[, , 2, ] <- rsdrms$ratio
 rsdrms_plot[, , 4, ] <- rsdrms$p.val
-\donttest{
+\dontrun{
 PlotVsLTime(rsdrms_plot, toptitle = "Ratio ensemble spread / RMSE", ytitle = "", 
            monini = 11, limits = c(-1, 1.3), listexp = c('CMIP5 IC3'), 
-           listobs = c('ERSST'), biglab = FALSE, siglev = TRUE, 
-           fileout = 'tos_rsdrms.eps')
+           listobs = c('ERSST'), biglab = FALSE, siglev = TRUE) 
 }
 
 }
diff --git a/man/Smoothing.Rd b/man/Smoothing.Rd
index 8d4a55871654d6159691f896bbed85434fe94da4..5769aa2b3763834651398bbecbd1268236264967 100644
--- a/man/Smoothing.Rd
+++ b/man/Smoothing.Rd
@@ -37,9 +37,8 @@ smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
 smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
 smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
 smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
- \donttest{
+ \dontrun{
 PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
-       toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
-       fileout = "tos_smoothed_ano.png")
+       toptitle = "Smoothed Mediterranean mean SST", ytitle = "K")
  }
 }
diff --git a/man/Spread.Rd b/man/Spread.Rd
index 26e289e919178e5b945a2477e0f266992c7fcf08..e26bc14551fcbdbdfa5eaad2a9628119c9b65210 100644
--- a/man/Spread.Rd
+++ b/man/Spread.Rd
@@ -74,27 +74,27 @@ smooth_ano_exp_m_sub <- smooth_ano_exp - InsertDim(MeanDims(smooth_ano_exp, 'mem
                                                   name = 'member')
 spread <- Spread(smooth_ano_exp_m_sub, compute_dim = c('member', 'sdate'))
 
-\donttest{
+\dontrun{
 PlotVsLTime(Reorder(spread$iqr, c('dataset', 'stats', 'ftime')), 
            toptitle = "Inter-Quartile Range between ensemble members",
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_iqr.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$maxmin, c('dataset', 'stats', 'ftime')), 
            toptitle = "Maximum minus minimum of the members", 
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_maxmin.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$sd, c('dataset', 'stats', 'ftime')), 
            toptitle = "Standard deviation of the members", 
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_sd.png')
+           hlines = c(0))
 PlotVsLTime(Reorder(spread$mad, c('dataset', 'stats', 'ftime')), 
            toptitle = "Median Absolute Deviation of the members",
            ytitle = "K", monini = 11, limits = NULL, 
            listexp = c('CMIP5 IC3'), listobs = c('ERSST'), biglab = FALSE, 
-           hlines = c(0), fileout = 'tos_mad.png')
+           hlines = c(0))
 }
 
 }
diff --git a/tests/testthat/test-MeanDims.R b/tests/testthat/test-MeanDims.R
index 9c7c5666a6246759ea7abed3f478bd85f8749f2a..a431b995422d6cff9050d5f7a2248da0f15451bd 100644
--- a/tests/testthat/test-MeanDims.R
+++ b/tests/testthat/test-MeanDims.R
@@ -57,7 +57,10 @@ test_that("1. Input checks", {
   MeanDims(dat1, dims = 'ftime', na.rm = na.omit),
   "Parameter 'na.rm' must be one logical value."
   )
-
+  expect_error(
+  MeanDims(dat1, dims = 'ftime', drop = 'selected'),
+  "Parameter 'drop' must be one logical value."
+  )
 })
 
 ##############################################
@@ -79,6 +82,26 @@ test_that("2. Output checks: dat1", {
   MeanDims(dat1, dims = c('sdate'))[1:2],
   c(3, 8)
   )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:2, drop = F)),
+  c(dat = 1, sdate = 1, ftime = 4)
+  )
+  expect_equal(
+  as.vector(drop(MeanDims(dat1, dims = 1:2, drop = F))),
+  as.vector(MeanDims(dat1, dims = 1:2, drop = T))
+  )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:3, drop = F)),
+  c(dat = 1, sdate = 1, ftime = 1)
+  )
+  expect_equal(
+  dim(MeanDims(dat1, dims = 1:3, drop = T)),
+  1
+  )
+  expect_equal(
+  as.vector(drop(MeanDims(dat1, dims = 1:3, drop = F))),
+  as.vector(MeanDims(dat1, dims = 1:3, drop = T))
+  )
 
 })
 
@@ -126,6 +149,10 @@ test_that("5. Output checks: dat4", {
     length(MeanDims(dat4, dims = 1)),
     1
   )
+  expect_equal(
+    dim(MeanDims(dat4, dims = 1, drop = F)),
+    1
+  )
 })
 ##############################################
 
diff --git a/tests/testthat/test-RatioRMS.R b/tests/testthat/test-RatioRMS.R
index b70d6fb5404998d7d737579a57012dd378383317..11df46ccf04d14e305059163e5c5d670906e57bc 100644
--- a/tests/testthat/test-RatioRMS.R
+++ b/tests/testthat/test-RatioRMS.R
@@ -114,12 +114,12 @@ names(RatioRMS(exp2_1, exp2_2, obs2)),
 c('ratiorms', 'p.val')
 )
 expect_equal(
-RatioRMS(exp2_1, exp2_2, obs2)$p.val,
+as.vector(RatioRMS(exp2_1, exp2_2, obs2)$p.val),
 0.7418331,
 tolerance = 0.0001
 )
 expect_equal(
-RatioRMS(exp2_1, exp2_2, obs2)$ratiorms,
+as.vector(RatioRMS(exp2_1, exp2_2, obs2)$ratiorms),
 0.8931399,
 tolerance = 0.0001
 )