diff --git a/NAMESPACE b/NAMESPACE
index e3226a5a05d98a2802001d21c5e6708f9d79b700..47ee4dc92d747aebb9591c68cd42617cd21a0648 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -20,5 +20,4 @@ export(WeightedMean)
 import(PCICt)
 import(climdex.pcic)
 import(multiApply)
-importFrom(plyr,aaply)
 importFrom(stats,quantile)
diff --git a/R/WeightedMean.R b/R/WeightedMean.R
index 393f56f680ffba9c459d69884043ed55f79a036b..ffcd27383c04b4aee189435564b9e92f810e7cbe 100644
--- a/R/WeightedMean.R
+++ b/R/WeightedMean.R
@@ -1,101 +1,173 @@
 #'Calculate spatial area-weighted average of multidimensional arrays
 #'
+#'This function computes a spatial area-weighted average of n-dimensional arrays 
+#'being possible to select a region and to add a mask to be applied when 
+#'computing the average.
 #'
-#'This function computes a spatial area-weighted average of n-dimensional arrays being possible to select a region and to add a mask to be applied when computing the average.
+#'@param data A numeric array with named dimensions, representing the data to be
+#'  applied the weights. It should have at least the latitude dimension and it 
+#'  can have more other dimensions.
+#'@param lon A numeric vector of longitude locations of the cell centers of the 
+#'  grid of \code{data}. This vector must be of the same length as the longitude 
+#'  dimension in the parameter \code{data} (in degrees).
+#'@param lat A numeric vector of latitude locations of the cell centers of the 
+#'  grid of \code{data}. This vector must be of the same length as the latitude 
+#'  dimension in the parameter \code{data} (in degrees). 
+#'@param region A vector of length four indicating the minimum longitude, the 
+#'  maximum longitude, the minimum latitude and the maximum latitude of the 
+#'  region to be averaged.
+#'@param mask A matrix with the same spatial dimensions of \code{data}. It can 
+#'  contain either a) TRUE where the value at that position is to be accounted 
+#'  for and FALSE where not, or b) numeric values, where those greater or equal 
+#'  to 0.5 are to be accounted for, and those smaller are not. Attention: if the 
+#'  longitude and latitude dimensions of the data and mask coincide in length, 
+#'  the user must ensure the dimensions of the mask are in the same order as the 
+#'  dimensions in the array provided in the parameter \code{data}.
+#'@param londim A character string indicating the name of the longitudinal
+#'  dimension. The default value is 'lon'.
+#'@param latdim A character string indicating the name of the latitudinal
+#'  dimension. The default value is 'lat'.
+#'@param na.rm A logical value indicating whether missing values should be 
+#'  stripped before the computation proceeds, by default it is set to TRUE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
 #'
-#'@param data An array with minimum two dimensions of latitude and longitude.
-#'@param lon Numeric vector of longitude locations of the cell centers of the grid of \code{data}. This vector must be the same length as the longitude dimension in the parameter \code{data}.
-#'@param lat Numeric vector of latitude locations of the cell centers of the grid of \code{data}. This vector must be the same length as the latitude dimension in the parameter \code{data}. 
-#'@param region A vector of length four indicating the minimum longitude, the maximum longitude, the minimum latitude and the maximum latitude of the region to be averaged.
-#'@param mask A matrix with the same spatial dimensions of \code{data}. It can contain either a) TRUE where the value at that position is to be accounted for and FALSE where not, or b) numeric values, where those greater or equal to 0.5 are to be accounted for, and those smaller are not. Attention: if the longitude and latitude dimensions of the data and mask coincide in length, the user must ensure the dimensions of the mask are in the same order as the dimensions in the array provided in the parameter \code{data}.
-#'@param londim An integer number indicating the position of the longitude dimension in the \code{data} object.
-#'@param latdim An integer number indicating the position of the latitude dimension in the \code{data} object.
+#'@return An array, matrix or vector containig the area-weighted average with 
+#'the same dimensions as \code{data}, except for the spatial longitude and 
+#'latitude dimensions, which disappear. 
 #'
-#'@return An array, matrix or vector containig the area-weighted average with the same dimensions as \code{data}, except for the spatial longitude and latitude dimensions, which disappear.
-#'
-#'@importFrom plyr aaply
 #'@examples
-#'##Example synthetic data 1:
+#'# Example 1:
 #'data <- 1:(2 * 3 * 4 * 5)
 #'dim(data) <- c(lon = 2, lat = 3, time = 4, model = 5)
 #'lat <- c(1, 10, 20)
 #'lon <- c(1, 10)
-#'
-#'a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL, 
-#'                  mask = NULL, londim = 1, latdim = 2)
-#'str(a)
+#'a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL)
 #'
 #'mask <- c(0, 1, 0, 1, 0, 1)
 #'dim(mask) <- c(lon = 2, lat = 3)
-#'a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL, 
-#'                  mask = mask, londim = 1, latdim = 2)
-#'str(a)
+#'a <- WeightedMean(data = data, lon = lon, lat = lat, mask = mask)
 #'
 #'region <- c(1, 10, 1, 10)
 #'a <- WeightedMean(data = data, lon = lon, lat = lat, region = region, 
-#'                  mask = mask, londim = 1, latdim = 2)
-#'str(a)
+#'                  mask = mask)
 #'
-#'##Example synthetic data:
+#'# Example 2:
 #'data <- 1:(2 * 3 * 4)
-#'dim(data) <- c(lon = 2, lat = 3, time=4)
+#'dim(data) <- c(lon = 2, lat = 3, time = 4)
 #'lat <- c(1, 10, 20)
 #'lon <- c(1, 10)
-#'
-#'a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL,
-#'                  mask = NULL, londim = 1, latdim = 2)
-#'str(a)
+#'a <- WeightedMean(data = data, lon = lon, lat = lat)
+#' 
+#'@import multiApply
 #'@export
-WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, londim = NULL, latdim = NULL){
-  if (is.null(data) | is.null(lon) | is.null(lat)) {
-    stop("Parameter 'data', 'lon' and 'lat' cannot be NULL.")
+WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, 
+                         londim = 'lon', latdim = 'lat', na.rm = TRUE, 
+                         ncores = NULL) {
+  # Check inputs
+  # data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
   }
-  if (!is.numeric(data) | !is.numeric(lon)| !is.numeric(lat)) {
-    stop("Parameter 'data', 'lon' and 'lat' must be a numeric.")
+  if (!is.array(data)) {
+    stop("Parameter 'data' must be a numeric array.")
   }
-  if (length(dim(data)) < 2) {
-    stop("Parameter 'data' needs to have dimensions lon and lat.")
+  dim_names <- names(dim(data))
+  if (is.null(dim_names)) {
+    stop("Parameter 'data' must have dimension names.")
   }
-  if (!is.null(dim(lat)) | !is.null(dim(lon))) {
-    stop("Parameter 'lon' and lat' need to be a vector.")
+  # lon, lat
+  if (is.null(lon) | is.null(lat)) {
+    stop("Parameters 'lon' and 'lat' cannot be NULL.")
   }
-  dim_names <- names(dim(data))
-  dims <- 1 : length(dim(data))
-  if (is.null(londim)) {
-    if (!is.null(dim_names)) {
-      londim <- which(dim_names == 'lon')
-    }
-    if(is.null(londim)) {
-      londim <- which(dim(data) == length(lon))
-    }
-    if (length(londim) == 0)  {
-      stop("No longitudinal dimension provided in parameter 'londim' nor as attribute of parameter 'data'.")
-    }
+  if (!is.numeric(lon) | !is.numeric(lat)) {
+    stop("Parameters 'lon' and 'lat' must be numeric.")
   }
-  if (is.null(latdim)) {
-    if (!is.null(dim_names)) {
-      latdim <- which(dim_names == 'lat')
+  if (!is.null(dim(lon)) | !is.null(dim(lat))) {
+    if (length(dim(lon)) == 1 & length(dim(lat)) == 1) {
+      lon <- as.vector(lon)
+      lat <- as.vector(lat)
+    } else {
+      stop("Parameters 'lon' and 'lat' need to be a vector.")
     }
-    if (is.null(latdim)) {
-      latdim <- which(dim(data) == length(lat))
-    }  
-    if (length(latdim) == 0) {
-        stop("No latitudinal dimension provided in parameter 'latdim' nor as attribute of parameter 'data'.")
+  }
+  # londim
+  if (is.numeric(londim)) {
+    warning("Numeric 'londim' is deprecated, use dimension names instead. The ", 
+            "corresponding dimension name will be assigned.")
+    londim <- dim_names[londim]
+  }
+  if (!is.character(londim)) {
+    stop("Parameter 'londim' must be a character string.")
+  }
+  if (length(londim) > 1) {
+    warning("Parameter 'londim' must be of length 1. Only the first value ", 
+            "will be used.")
+    londim <- londim[1]
+  }
+  if (!londim %in% names(dim(data))) {
+    stop("Parameter 'londim' is not found in 'data'.")
+  }
+  if (dim(data)[londim] != length(lon)) {
+    stop(paste0("The longitudinal dimension of parameter 'data' must be of the ", 
+                "same length as parameter 'lon'."))
+  }
+  # latdim
+  if (is.numeric(latdim)) {
+    warning("Numeric 'latdim' is deprecated, use dimension names instead. The ", 
+            "corresponding dimension name will be assigned.")
+    latdim <- dim_names[latdim]
+  }
+  if (!is.character(latdim)) {
+    stop("Parameter 'latdim' must be a character string.")
+  }
+  if (length(latdim) > 1) {
+    warning("Parameter 'latdim' must be of length 1. Only the first value ", 
+            "will be used.")
+    latdim <- latdim[1]
+  }
+  if (!latdim %in% names(dim(data))) {
+    stop("Parameter 'latdim' is not found in 'data'.")
+  }
+  if (dim(data)[latdim] != length(lat)) {
+    stop(paste0("The latitudinal dimension of parameter 'data' must be of the ", 
+                "same length as parameter 'lat'."))
+  }
+  lon_pos <- which(dim_names == londim)
+  lat_pos <- which(dim_names == latdim)
+  # region
+  if (!is.null(region)) {
+    if (length(region) != 4) {
+      stop(paste0("The region argument has to be of length four indicating ", 
+                  "the minimum longitude, the maximum longitude, the minimum ", 
+                  "latitude and the maximum latitude of the region to be averaged."))
     }
   }
-  if (londim == latdim) {
-    stop("Parameter 'londim' and 'latdim' cannot be equal.")
+  # mask
+  if (!is.null(mask)) {
+    if (!all(dim(data)[which(names(dim(data)) %in% c(londim, latdim))] %in% 
+             dim(mask)[which(names(dim(mask)) %in% c(londim, latdim))])) {
+      stop("Parameter 'mask' must have the same spatial dimensions of data.")
+    }
   }
-  if (dim(data)[londim] != length(lon)){
-    stop("The longitudinal dimension of parameter 'data' must be the same length of parameter 'lon'.")
+  # na.rm
+  if (!is.logical(na.rm) | length(na.rm) > 1) {
+    stop("Parameter 'na.rm' must be one logical value.")
   }
-  if (dim(data)[latdim] != length(lat)){
-    stop("The latitudinal dimension of parameter 'data' must be the same length of parameter 'lat'.")
+  # ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | length(ncores) > 1) {
+      stop("Parameter 'ncores' must be either NULL or a positive integer.")
+    } else if (ncores %% 1 != 0 | ncores <= 0) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
   }
+
   nlon <- length(lon)
   nlat <- length(lat)
   if (!is.null(region)) {
-    aux <- SelBox(data, lon = lon, lat = lat, region = region, londim = londim, latdim = latdim, mask = mask)
+    aux <- SelBox(data, lon = lon, lat = lat, region = region, londim = lon_pos, 
+                  latdim = lat_pos, mask = mask)
     data <- aux$data
     lon <- aux$lon
     lat <- aux$lat
@@ -103,24 +175,8 @@ WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, londim = NU
       mask <- aux$mask
     }
   }
-  if (length(dim(data)) > 2) {
-    wtmean <- aaply(data, .margins = dims[c(-londim, -latdim)], 
-                          .fun = .WeightedMean, lon, lat, mask, .drop = FALSE)
-    dim(wtmean) <- dim(data)[-c(latdim,londim)]
-  } else {
-    wtmean <- .WeightedMean(data, lon = lon, lat = lat, mask = mask)
-  }
-  if(length(dim(data)) > 3) {
-    if (is.null(dim_names)) {
-      dim_names <- paste0("dim", 1:length(dim(data)))
-    }
-    names(dim(wtmean)) = dim_names[-c(londim, latdim)]
-  } else {
-    attributes(wtmean) <- NULL
-  }
-  wtmean
-}
-.WeightedMean <- function(data, lon, lat, mask) {
+
+  # Compute the weights
   cosphi <- t(array(cos(lat * pi / 180), dim = c(length(lat), length(lon))))
   nblat <- length(lat)
   nblon <- length(lon)
@@ -131,23 +187,32 @@ WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, londim = NU
   dlat <- c(dlat, dlat[1])
   dlat <- t(array(dlat, dim = c(nblat, nblon)))
   weight <- (dlon * dlat * cosphi) 
-  if ((is.null(mask) == FALSE) && all(dim(data) != dim(mask))) {
-    stop("The provided parameter 'mask' must have the same size as the parameter 'data'.")
-  }
-  if ((nblat == dim(data)[1]) && (nblon == dim(data)[2])) {
-    data <- t(data)
-    if (!is.null(mask)) {
-      mask <- t(mask)
-    }
-  } 
-  if ((nblon != dim(data)[1]) || (nblat != dim(data)[2])) {
-    stop("The parameter 'data' needs to have at least two dimensions of 'lon' and 'lat'.")
+
+  if (is.null(mask)) {
+    res <- Apply(data = list(data), 
+                 target_dims = c(londim, latdim), 
+                 fun = .WeightedMean,
+                 mask = NULL, 
+                 weight = weight, 
+                 na.rm = na.rm, 
+                 ncores = ncores)$output1
+  } else {
+    res <- Apply(data = list(data, mask), 
+                 target_dims = c(londim, latdim), 
+                 fun = .WeightedMean,
+                 weight = weight, 
+                 na.rm = na.rm, 
+                 ncores = ncores)$output1
   }
-  if(!is.null(mask)) {
+  return(res)
+}
+
+.WeightedMean <- function(data, mask = NULL, weight, na.rm = TRUE) {
+  if (!is.null(mask)) {
     data[mask < 0.5] <- NA
   }
   weight[is.na(data)] <- NA
-  coeff <- sum(weight, na.rm = TRUE)
-  mean <- sum(weight * data, na.rm = TRUE) / coeff
-  output <- mean
+  coeff <- sum(weight, na.rm = na.rm)
+  mean <- sum(weight * data, na.rm = na.rm) / coeff
+  return(mean)
 }
diff --git a/man/WeightedMean.Rd b/man/WeightedMean.Rd
index ac3411cba7fc3fadea3f557220d7a6607aedad12..ad59157e2e39007c727ba7ccdcb806dda3ce989b 100644
--- a/man/WeightedMean.Rd
+++ b/man/WeightedMean.Rd
@@ -10,60 +10,80 @@ WeightedMean(
   lat,
   region = NULL,
   mask = NULL,
-  londim = NULL,
-  latdim = NULL
+  londim = "lon",
+  latdim = "lat",
+  na.rm = TRUE,
+  ncores = NULL
 )
 }
 \arguments{
-\item{data}{An array with minimum two dimensions of latitude and longitude.}
+\item{data}{A numeric array with named dimensions, representing the data to be
+applied the weights. It should have at least the latitude dimension and it 
+can have more other dimensions.}
 
-\item{lon}{Numeric vector of longitude locations of the cell centers of the grid of \code{data}. This vector must be the same length as the longitude dimension in the parameter \code{data}.}
+\item{lon}{A numeric vector of longitude locations of the cell centers of the 
+grid of \code{data}. This vector must be of the same length as the longitude 
+dimension in the parameter \code{data} (in degrees).}
 
-\item{lat}{Numeric vector of latitude locations of the cell centers of the grid of \code{data}. This vector must be the same length as the latitude dimension in the parameter \code{data}.}
+\item{lat}{A numeric vector of latitude locations of the cell centers of the 
+grid of \code{data}. This vector must be of the same length as the latitude 
+dimension in the parameter \code{data} (in degrees).}
 
-\item{region}{A vector of length four indicating the minimum longitude, the maximum longitude, the minimum latitude and the maximum latitude of the region to be averaged.}
+\item{region}{A vector of length four indicating the minimum longitude, the 
+maximum longitude, the minimum latitude and the maximum latitude of the 
+region to be averaged.}
 
-\item{mask}{A matrix with the same spatial dimensions of \code{data}. It can contain either a) TRUE where the value at that position is to be accounted for and FALSE where not, or b) numeric values, where those greater or equal to 0.5 are to be accounted for, and those smaller are not. Attention: if the longitude and latitude dimensions of the data and mask coincide in length, the user must ensure the dimensions of the mask are in the same order as the dimensions in the array provided in the parameter \code{data}.}
+\item{mask}{A matrix with the same spatial dimensions of \code{data}. It can 
+contain either a) TRUE where the value at that position is to be accounted 
+for and FALSE where not, or b) numeric values, where those greater or equal 
+to 0.5 are to be accounted for, and those smaller are not. Attention: if the 
+longitude and latitude dimensions of the data and mask coincide in length, 
+the user must ensure the dimensions of the mask are in the same order as the 
+dimensions in the array provided in the parameter \code{data}.}
 
-\item{londim}{An integer number indicating the position of the longitude dimension in the \code{data} object.}
+\item{londim}{A character string indicating the name of the longitudinal
+dimension. The default value is 'lon'.}
 
-\item{latdim}{An integer number indicating the position of the latitude dimension in the \code{data} object.}
+\item{latdim}{A character string indicating the name of the latitudinal
+dimension. The default value is 'lat'.}
+
+\item{na.rm}{A logical value indicating whether missing values should be 
+stripped before the computation proceeds, by default it is set to TRUE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
 }
 \value{
-An array, matrix or vector containig the area-weighted average with the same dimensions as \code{data}, except for the spatial longitude and latitude dimensions, which disappear.
+An array, matrix or vector containig the area-weighted average with 
+the same dimensions as \code{data}, except for the spatial longitude and 
+latitude dimensions, which disappear.
 }
 \description{
-This function computes a spatial area-weighted average of n-dimensional arrays being possible to select a region and to add a mask to be applied when computing the average.
+This function computes a spatial area-weighted average of n-dimensional arrays 
+being possible to select a region and to add a mask to be applied when 
+computing the average.
 }
 \examples{
-##Example synthetic data 1:
+# Example 1:
 data <- 1:(2 * 3 * 4 * 5)
 dim(data) <- c(lon = 2, lat = 3, time = 4, model = 5)
 lat <- c(1, 10, 20)
 lon <- c(1, 10)
-
-a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL, 
-                 mask = NULL, londim = 1, latdim = 2)
-str(a)
+a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL)
 
 mask <- c(0, 1, 0, 1, 0, 1)
 dim(mask) <- c(lon = 2, lat = 3)
-a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL, 
-                 mask = mask, londim = 1, latdim = 2)
-str(a)
+a <- WeightedMean(data = data, lon = lon, lat = lat, mask = mask)
 
 region <- c(1, 10, 1, 10)
 a <- WeightedMean(data = data, lon = lon, lat = lat, region = region, 
-                 mask = mask, londim = 1, latdim = 2)
-str(a)
+                 mask = mask)
 
-##Example synthetic data:
+# Example 2:
 data <- 1:(2 * 3 * 4)
-dim(data) <- c(lon = 2, lat = 3, time=4)
+dim(data) <- c(lon = 2, lat = 3, time = 4)
 lat <- c(1, 10, 20)
 lon <- c(1, 10)
+a <- WeightedMean(data = data, lon = lon, lat = lat)
 
-a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL,
-                 mask = NULL, londim = 1, latdim = 2)
-str(a)
 }
diff --git a/tests/testthat/test-WeightedCells.R b/tests/testthat/test-WeightedCells.R
index dcf171d6f29eecfdcb2819aec595d6b279eafc64..a5e9dd156e75f4833cf1f6707e293eb7e5a6ebf2 100644
--- a/tests/testthat/test-WeightedCells.R
+++ b/tests/testthat/test-WeightedCells.R
@@ -1,4 +1,4 @@
-context("s2dv::WeightedCells tests")
+context("ClimProjDiags::WeightedCells tests")
 
 ##############################################
 
@@ -66,7 +66,6 @@ test_that("1. Input checks", {
     WeightedCells(data1, lat1, ncores = 1.5),
     "Parameter 'ncores' must be either NULL or a positive integer."
   )
-
 })
 
 ##############################################
diff --git a/tests/testthat/test-WeightedMean.R b/tests/testthat/test-WeightedMean.R
index 2c5a23176bfe912072b9d12fbbdbafeb55615c5d..24144db92269ada883dc6408e34c0b55710f503a 100644
--- a/tests/testthat/test-WeightedMean.R
+++ b/tests/testthat/test-WeightedMean.R
@@ -1,4 +1,7 @@
-context("WeightedMean tests")
+context("ClimProjDiags::WeightedMean tests")
+
+##############################################
+
 set.seed(1)
 dat1 <- array(rnorm(10000), dim = c(lat = 50, lon = 100))
 lat1 <- seq(-90, 90, length.out = 50)
@@ -19,53 +22,227 @@ dat4 <- array(rnorm(10000), dim = c(lat = 180, lon = 360, sdate = 2, time = 1))
 lat4 <- seq(-89.5, 89.5, length.out = 180)
 lon4 <- seq(-179.5, 179.5, length.out = 360)
 
+set.seed(2)
+mask <- array(rnorm(180*360), dim = c(lat = 180, lon = 360))
+region <- c(1, 55, 2, 67)
+
+dat5 <- dat4
+dat5[1,1,1,1] <- NA
+
+##############################################
+
+test_that("1. Input checks", {
+  # data, lon, lat
+  expect_error(
+    WeightedMean(c(), lon = 1:10, lat = 1:10),
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    WeightedMean(1:10, lon = 1:10, lat = 1:10),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    WeightedMean(array(1:10), lon = 1:10, lat = 1:10),
+    "Parameter 'data' must have dimension names."
+  )
+  # lon, lat
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10)), lon = NULL, lat = 1:10),
+    "Parameters 'lon' and 'lat' cannot be NULL."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10)), lon = 'a', lat = 1:10),
+    "Parameters 'lon' and 'lat' must be numeric."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10)), 
+                 lon = 1:10, lat = array(1:10, dim = c(lon = 10, lat = 1))),
+    "Parameters 'lon' and 'lat' need to be a vector."
+  )
+  # londim
+  expect_warning(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, londim = 1, latdim = 2),
+    "Numeric 'londim' is deprecated, use dimension names instead. The ", 
+    "corresponding dimension name will be assigned."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, londim = list(1,2), latdim = 2),
+    "Parameter 'londim' must be a character string."
+  )
+  expect_warning(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, londim = c('lon', 'lat'), latdim = 'lat'),
+    "Parameter 'londim' must be of length 1. Only the first value ", 
+    "will be used."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lons = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, londim = 'lon', latdim = 'lat'),
+    "Parameter 'londim' is not found in 'data'."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:11, lat = 1, londim = 'lon', latdim = 'lat'),
+    paste0("The longitudinal dimension of parameter 'data' must be of the ", 
+           "same length as parameter 'lon'.")
+  )
+  # latdim
+  expect_warning(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, latdim = 2),
+    "Numeric 'latdim' is deprecated, use dimension names instead. The ", 
+    "corresponding dimension name will be assigned."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, latdim = list(1,2)),
+    "Parameter 'latdim' must be a character string."
+  )
+  expect_warning(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1, latdim = c('lat', 'lon')),
+    "Parameter 'latdim' must be of length 1. Only the first value ", 
+    "will be used."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10, lats = 1)), 
+                 lon = 1:10, lat = 1, londim = 'lon', latdim = 'lat'),
+    "Parameter 'latdim' is not found in 'data'."
+  )
+  expect_error(
+    WeightedMean(array(1:10, dim = c(lon = 10, lat = 1)), 
+                 lon = 1:10, lat = 1:2, londim = 'lon', latdim = 'lat'),
+    paste0("The latitudinal dimension of parameter 'data' must be of the ", 
+           "same length as parameter 'lat'.")
+  )
+  # region
+  expect_error(
+    WeightedMean(dat1, lat = lat1, lon = lon1, region = 1:3),
+    paste0("The region argument has to be of length four indicating ", 
+           "the minimum longitude, the maximum longitude, the minimum ", 
+           "latitude and the maximum latitude of the region to be averaged.")
+  )
+  # mask
+  expect_error(
+    WeightedMean(dat1, lat = lat1, lon = lon1, mask = 1:3),
+    "Parameter 'mask' must have the same spatial dimensions of data."
+  )
+  # na.rm
+  expect_error(
+    WeightedMean(dat1, lat = lat1, lon = lon1, na.rm = 1:3),
+    "Parameter 'na.rm' must be one logical value."
+  )
+  # ncores
+  expect_error(
+    WeightedMean(dat1, lat = lat1, lon = lon1, ncores = 'a'),
+    "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+  expect_error(
+    WeightedMean(dat1, lat = lat1, lon = lon1, ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+})
+
+##############################################
+test_that("2. Output test: dat1", {
+  expect_equal(
+    dim(WeightedMean(dat1, lat = lat1, lon = lon1)),
+    NULL
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat1, lat = lat1, lon = lon1)),
+    -0.009785971,
+    tolerance = 0.0001
+  )
+})
+
+##############################################
 
-test_that("dat1", {
-expect_equal(
-dim(WeightedMean(dat1, lat = lat1, lon = lon1)),
-NULL
-)
-expect_equal(
-as.vector(WeightedMean(dat1, lat = lat1, lon = lon1)),
--0.009785971,
-tolerance = 0.0001
-)
+test_that("3. Output test: dat2", {
+  expect_equal(
+    dim(WeightedMean(dat2, lat = lat2, lon = lon2)),
+    c(sdate = 2)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat2, lat = lat2, lon = lon2)),
+    c(-0.005799676, -0.007599831),
+    tolerance = 0.0001
+  )
 })
 
-test_that("dat2", {
-expect_equal(
-dim(WeightedMean(dat2, lat = lat2, lon = lon2)),
-NULL
-)
-expect_equal(
-as.vector(WeightedMean(dat2, lat = lat2, lon = lon2)),
-c(-0.005799676, -0.007599831),
-tolerance = 0.0001
-)
+test_that("4. Output test: dat3", {
+  expect_equal(
+    dim(WeightedMean(dat3, lat = lat3, lon = lon3)),
+    c(sdate = 2)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat3, lat = lat3, lon = lon3)),
+    c(-0.0253997, 0.0132251),
+    tolerance = 0.0001
+  )
+})
+
+##############################################
+
+test_that("5. Output test: dat4", {
+  expect_equal(
+    dim(WeightedMean(dat4, lat = lat4, lon = lon4)),
+    c(sdate = 2, time = 1)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat4, lat = lat4, lon = lon4)),
+    c(-0.005799676, -0.007599831),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    dim(WeightedMean(dat4, lat = lat4, lon = lon4, mask = mask, region = region)),
+    c(sdate = 2, time = 1)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat4, lat = lat4, lon = lon4, 
+                           mask = mask, region = region)),
+    c( 0.03610427, 0.02197599),
+    tolerance = 0.0001
+  )
 })
 
-test_that("dat3", {
-expect_equal(
-dim(WeightedMean(dat3, lat = lat3, lon = lon3)),
-NULL
-)
-expect_equal(
-as.vector(WeightedMean(dat3, lat = lat3, lon = lon3)),
-c(-0.0253997, 0.0132251),
-tolerance = 0.0001
-)
+##############################################
 
+test_that("6. Output test: dat4 (mask, region, ncores)", {
+  expect_equal(
+    dim(WeightedMean(dat4, lat = lat4, lon = lon4, mask = mask, region = region)),
+    c(sdate = 2, time = 1)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat4, lat = lat4, lon = lon4, mask = mask, 
+                           region = region, ncores = 1)),
+    c( 0.03610427, 0.02197599),
+    tolerance = 0.0001
+  )
 })
 
+##############################################
 
-test_that("dat4", {
-expect_equal(
-dim(WeightedMean(dat4, lat = lat4, lon = lon4)),
-c(sdate = 2, time = 1)
-)
-expect_equal(
-as.vector(WeightedMean(dat4, lat = lat4, lon = lon4)),
-c(-0.005799676, -0.007599831),
-tolerance = 0.0001
-)
+test_that("7. Output test: dat5", {
+  expect_equal(
+    dim(WeightedMean(dat5, lat = lat4, lon = lon4, na.rm = FALSE)),
+    c(sdate = 2, time = 1)
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat5, lat = lat4, lon = lon4, na.rm = FALSE)),
+    c(NA, -0.007599831),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat5, lat = lat4, lon = lon4, na.rm = TRUE)),
+    as.vector(WeightedMean(dat4, lat = lat4, lon = lon4, na.rm = TRUE)),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(WeightedMean(dat4, lat = lat4, lon = lon4, na.rm = FALSE)),
+    c(-0.005799676, -0.007599831),
+    tolerance = 0.0001
+  )
 })