diff --git a/DESCRIPTION b/DESCRIPTION
index 2530c9b9d96f98d66017aedf71cb46ec1ebd9eb6..0c5151978a47f47d02a0e1dfe2ea051c28b6d17a 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.2.1
+Version: 0.3.0
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut"), comment = c(ORCID = "0000-0001-8568-3071")),
@@ -29,7 +29,7 @@ Suggests:
     testthat,
     markdown,
     rmarkdown
-License: Apache License 2.0
+License: GPL-3
 URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
 BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/-/issues
 Encoding: UTF-8
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/NEWS.md b/NEWS.md
index 5b5d43ded490ccb02b55d0f2548994e426d51a68..dd95442f5ba4fc20e1e46608107618ec7440c725 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,10 @@
+# 0.3.0 (Release date: 2023-02-28)
+- SelBox() and ShiftLon() to accept non-numerical data input  
+- SelBox() uses the latitude and longitude dimension name instead of index  
+- WeightedMean() uses multiApply::Apply inside  
+
 # 0.2.1 (Release date: 2022-12-01)
-- Fix the mistake that function "WeightedCells" was not included in the last submission.
+- Fix the mistake that function "WeightedCells" was not included in the last submission.  
 
 # 0.2.0 (Release date: 2022-11-04)
 - New functions: ShiftLon, WeightedCells  
diff --git a/R/SelBox.R b/R/SelBox.R
index 5b7cd5982583f278142a8321c250539115676027..d743aef4728aef3cd1e2ddd72b1e888f2dc87da2 100644
--- a/R/SelBox.R
+++ b/R/SelBox.R
@@ -1,102 +1,154 @@
-#'Select apatial region from multidimensional arrays
+#'Select spatial region from multidimensional arrays
 #'
-#'@description This function subsets an spatial region from spatial data giving a vector with the maximum and minimum of latitudes and longitudes of the selected region.
+#'@description Subset a spatial region from spatial data giving a vector with
+#'  the maximum and minimum of latitudes and longitudes of the selected region.
 #'
 #'@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}'. 
-#'@param lat Numeric vector of latitude locations of the cell centers of the grid of \code{data}'.
-#'@param region A vector of length four indicating the minimum longitude, the maximum longitude, the minimum latitude and the maximum latitude.
-#'@param londim An integer number indicating the position of the longitude dimension in the \code{data} object. If NULL (by deafault), the function search for a dimension call 'lon' in the \code{data} input.
-#'@param latdim An integer number indicating the position of the latitude dimension in the \code{data} object.  If NULL (by deafault), the function search for a dimension call 'lat' in the \code{data} input.
+#'@param lon Numeric vector of longitude locations of the cell centers of the
+#'  grid of \code{data}'. 
+#'@param lat Numeric vector of latitude locations of the cell centers of the
+#'  grid of \code{data}'.
+#'@param region A vector of length four indicating the minimum longitude, the
+#'  maximum longitude, the minimum latitude and the maximum latitude.
+#'@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 mask A matrix with the same spatial dimensions of \code{data}.
 #'
 #'@return A list of length 4:
 #'\itemize{
-#'  \item\code{$data}{An array with the same dimensions as the input \code{data} array, but with spatial dimension reduced to the selected \code{region}}
-#'  \item\code{$lat}{A vector with the new corresponding latitudes for the selected \code{region}}
-#'  \item\code{$lon}{A vector with the new corresponding longitudes for the selected \code{region}}
-#'  \item\code{$mask}{If parameter \code{mask} is supplied, an array  with reduced length of the dimensions to the selected \code{region}. Otherwise, a NULL element is returned.}}
+#'  \item\code{$data}{
+#'    An array with the same dimensions as the input \code{data} array, but with
+#'    spatial dimension reduced to the selected \code{region}}
+#'  \item\code{$lat}{A vector with the new corresponding latitudes for the 
+#'    selected \code{region}}
+#'  \item\code{$lon}{A vector with the new corresponding longitudes for the 
+#'    selected \code{region}}
+#'  \item\code{$mask}{If parameter \code{mask} is supplied, an array with 
+#'    reduced length of the dimensions to the selected \code{region}. Otherwise, 
+#'    a NULL element is returned.}}
 #'
 #'@examples 
-#'## Example with synthetic data:
+#'# Example with synthetic data:
 #'data <- 1:(20 * 3 * 2 * 4)
 #'dim(data) <- c(lon = 20, lat = 3, time = 2, model = 4)
 #'lon <- seq(2, 40, 2)
 #'lat <- c(1, 5, 10)
 #'
 #'a <- SelBox(data = data, lon = lon, lat = lat, region = c(2, 20, 1, 5), 
-#'            londim = 1, latdim = 2, mask = NULL)
-#'str(a)
+#'            londim = "lon", latdim = "lat", mask = NULL)
 #'@export
-SelBox <- function(data, lon, lat, region, londim = NULL, latdim = NULL, mask = NULL) {
-  if (is.null(data) | is.null(lon) | is.null(lat) | is.null(region)){
-    stop("Parameters 'data', 'lon', 'lat' or 'region' cannot be NULL.")
-  }
-  if (!is.numeric(data) | !is.numeric(lon) | !is.numeric(lat) | !is.numeric(region)){
-    stop("Parameters 'data', 'lon', 'lat' or 'region' must be numeric.")
+SelBox <- function(data, lon, lat, region, londim = 'lon', latdim = 'lat', 
+                   mask = NULL) {
+  # Check inputs
+  # data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
   }
   if (!is.array(data) && !is.matrix(data)) {
-    stop("Parameter 'data' must be an array or matrix.")
-  }
-  if (!is.null(dim(lat)) | !is.null(dim(lon))) {
-    stop("Parameter 'lon' and lat' need to be a vector.")
-  }
-  if (length(region) != 4){
-    stop("The region argument has to be a vector of length four indicating the minimum longitude, the maximum longitude, the minimum latitude and the maximum latitude.")
-  }
-  dims <- 1:length(dim(data))
-  if (is.null(londim)) {
-    if ("lon" %in% names(dim(data))) {
-      londim <- which(names(dim(data)) == "lon")
-    } else if (length(lon) %in% dim(data)) {
-      londim <- which(dim(data) == length(lon))
-      if (length(londim) > 1) {
-        stop("More than one dimension of the parameter 'data' has the same length as 'lon' parameter.")
-      }
-    } else {
-      stop("Non of the dimensions of the parameter 'data' are of the same length as 'lon'.")
-    }
-  } 
-  if (is.null(latdim)) {
-    if ("lat" %in% names(dim(data))) {
-      latdim <- which(names(dim(data)) == "lat")
-    } else if (length(lat) %in% dim(data)) {
-      latdim <- which(dim(data) == length(lat))
-      if (length(latdim) > 1) {
-        stop("More than one dimension of the parameter 'data' has the same length as 'lat' parameter.")
-      }
+    stop("Parameter 'data' must be a numeric array or matrix.")
+  }
+  dim_names <- names(dim(data))
+  if (is.null(dim_names)) {
+    stop("Parameter 'data' must have dimension names.")
+  }
+
+  # lon, lat
+  if (is.null(lon) | is.null(lat)) {
+    stop("Parameters 'lon' and 'lat' cannot be NULL.")
+  }
+  if (!is.numeric(lon) | !is.numeric(lat)) {
+    stop("Parameters 'lon' and 'lat' must be numeric.")
+  }
+  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("Non of the dimensions of the parameter 'data' are of the same length as 'lat' parameter.")
+      stop("Parameters 'lon' and 'lat' need to be a vector.")
     }
   }
-  if (londim == latdim) {
-    stop("Parameter 'londim' and 'latdim' cannot be equal.")
+  # region
+  if (is.null(region)) {
+    stop("Parameter 'region' cannot be NULL.")
+  }
+  if (!is.numeric(region)) {
+    stop("Parameter 'region' must be numeric.")
+  }
+  if (length(region) != 4) {
+    stop("The region argument has to be a vector of length four indicating the ",
+         "minimum longitude, the maximum longitude, the minimum latitude and ", 
+         "the maximum latitude.")
+  }
+  # 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 (dim(data)[londim] != length(lon)){
-    stop("The longitudinal dimension of parameter 'data' must be of the same length of parameter  'lon'.")
+  if (length(londim) > 1) {
+    warning("Parameter 'londim' must be of length 1. Only the first value ", 
+            "will be used.")
+    londim <- londim[1]
   }
-  if (dim(data)[latdim] != length(lat)){
-    stop("The latitudinal dimension of parameter 'data' must be of the same length of parameter  'lat'.")
+  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'."))
+  }
+  # 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 (region[3] <= region[4]) {
-    LatIdx <- which( lat >= region[3] & lat <= region[4])
+    LatIdx <- which(lat >= region[3] & lat <= region[4])
   } else {
     LatIdx <- which(lat <= region[3] | lat >= region[4])
   }
-  #if (region[1] <= region[2]) {
-  #  LonIdx <- which(lon >= region[1] & lon <= region[2])
-  #} else {
-  #  LonIdx <- which(lon >= region[1] | lon <= region[2])
-  #}
+  # if (region[1] <= region[2]) {
+  #   LonIdx <- which(lon >= region[1] & lon <= region[2])
+  # } else {
+  #   LonIdx <- which(lon >= region[1] | lon <= region[2])
+  # }
   LonIdx <- Lon2Index(lon, lonmin = region[1], lonmax = region[2])
 
-  data <- Subset(data, along = londim, indices = LonIdx, drop = "none")
-  data <- Subset(data, along = latdim, indices = LatIdx, drop = "none")
+  data <- Subset(data, along = c(londim, latdim), indices = list(LonIdx, LatIdx), drop = "none")
+
   if (!is.null(mask)) {
-    mask <- Subset(mask, along = latdim - length(dim(data)) + 2, indices = LatIdx, drop = "none")
-    mask <- Subset(mask, along = londim - length(dim(data)) + 2, indices = LonIdx, drop = "none")
+    mask <- Subset(mask, along = c(londim, latdim), indices = list(LonIdx, LatIdx), drop = "none")
   } else {
     mask <- NULL
   }
-  list(data = data, lon = lon[LonIdx], lat = lat[LatIdx], mask = mask)
+  return(list(data = data, lon = lon[LonIdx], lat = lat[LatIdx], mask = mask))
 }
diff --git a/R/ShiftLon.R b/R/ShiftLon.R
index 88b034d929f246dc677ab80af4fac430d80e527b..1467108581c8bfe5adc5b1a2248fdd20db3bce0d 100644
--- a/R/ShiftLon.R
+++ b/R/ShiftLon.R
@@ -35,15 +35,20 @@
 #'  
 #'@import multiApply
 #'@export
-ShiftLon <- function(data, lon,  westB, lon_dim = 'lon', ncores = NULL) {
+ShiftLon <- function(data, lon, westB, lon_dim = 'lon', ncores = NULL) {
   
   # Check inputs
   ## data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
-  if (!is.array(data) | !is.numeric(data)) {
-    stop("Parameter 'data' must be a numeric array.")
+  if (is.vector(data)) {
+    data <- as.array(data)
+    names(dim(data)) <- lon_dim
+    warning("Parameter 'data' is a vector. Transfer it to an array and assign ", lon_dim, " as dimension name.")
+  }
+  if (!is.array(data)) {
+    stop("Parameter 'data' must be an array.")
   }
   if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
     stop("Parameter 'data' must have dimension names.")
@@ -106,13 +111,15 @@ ShiftLon <- function(data, lon,  westB, lon_dim = 'lon', ncores = NULL) {
   } else {
     new.lon <- c(lon[first:length(lon)],lon[1:(first-1)])
   }
-  ## Order to  monotonically increasing
+  ## Order to monotonically increasing
   if (!all(diff(new.lon) > 0)) {
     new.lon[(which(diff(new.lon) < 0) + 1):length(new.lon)] <- new.lon[(which(diff(new.lon) < 0) + 1):length(new.lon)] + 360
     
   }
   
   # Shifting the data
+  ori_dim <- dim(data)
+
   output <- Apply(data = data,
                   target_dims = lon_dim,
                   fun = .ShiftLon,
@@ -122,13 +129,16 @@ ShiftLon <- function(data, lon,  westB, lon_dim = 'lon', ncores = NULL) {
   
   # Shift new.lon back to start at westB
   new.lon <- new.lon + shft_westB_back
-  
+
+  # Change dimension order back
+  output <- aperm(output, match(names(ori_dim), names(dim(output))))
+
   return(list(data = output, lon = new.lon))
 }
 
 .ShiftLon <- function(data, lon, new.lon) {
   # data: [lon]
-  new.data <- NA * data
+  new.data <- data
   new.data[new.lon %in% lon] <- data[lon %in% new.lon, drop = F]
   new.data[!new.lon %in% lon] <- data[!lon %in% new.lon, drop = F]
   
diff --git a/R/WeightedMean.R b/R/WeightedMean.R
index 393f56f680ffba9c459d69884043ed55f79a036b..0932599ed672451ed7bd895abe3ed79aea257b22 100644
--- a/R/WeightedMean.R
+++ b/R/WeightedMean.R
@@ -1,101 +1,169 @@
 #'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'."))
+  }
+  # 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 = londim, 
+                  latdim = latdim, mask = mask)
     data <- aux$data
     lon <- aux$lon
     lat <- aux$lat
@@ -103,24 +171,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 +183,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/SelBox.Rd b/man/SelBox.Rd
index 38e3547bf6ae6d6e699e78821aaefb1311896535..1bc496beff0c876ff14323e80984321d290df3e0 100644
--- a/man/SelBox.Rd
+++ b/man/SelBox.Rd
@@ -2,44 +2,55 @@
 % Please edit documentation in R/SelBox.R
 \name{SelBox}
 \alias{SelBox}
-\title{Select apatial region from multidimensional arrays}
+\title{Select spatial region from multidimensional arrays}
 \usage{
-SelBox(data, lon, lat, region, londim = NULL, latdim = NULL, mask = NULL)
+SelBox(data, lon, lat, region, londim = "lon", latdim = "lat", mask = NULL)
 }
 \arguments{
 \item{data}{An array with minimum two dimensions of latitude and longitude.}
 
-\item{lon}{Numeric vector of longitude locations of the cell centers of the grid of \code{data}'.}
+\item{lon}{Numeric vector of longitude locations of the cell centers of the
+grid of \code{data}'.}
 
-\item{lat}{Numeric vector of latitude locations of the cell centers of the grid of \code{data}'.}
+\item{lat}{Numeric vector of latitude locations of the cell centers of the
+grid of \code{data}'.}
 
-\item{region}{A vector of length four indicating the minimum longitude, the maximum longitude, the minimum latitude and the maximum latitude.}
+\item{region}{A vector of length four indicating the minimum longitude, the
+maximum longitude, the minimum latitude and the maximum latitude.}
 
-\item{londim}{An integer number indicating the position of the longitude dimension in the \code{data} object. If NULL (by deafault), the function search for a dimension call 'lon' in the \code{data} input.}
+\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.  If NULL (by deafault), the function search for a dimension call 'lat' in the \code{data} input.}
+\item{latdim}{A character string indicating the name of the latitudinal
+dimension. The default value is 'lat'.}
 
 \item{mask}{A matrix with the same spatial dimensions of \code{data}.}
 }
 \value{
 A list of length 4:
 \itemize{
- \item\code{$data}{An array with the same dimensions as the input \code{data} array, but with spatial dimension reduced to the selected \code{region}}
- \item\code{$lat}{A vector with the new corresponding latitudes for the selected \code{region}}
- \item\code{$lon}{A vector with the new corresponding longitudes for the selected \code{region}}
- \item\code{$mask}{If parameter \code{mask} is supplied, an array  with reduced length of the dimensions to the selected \code{region}. Otherwise, a NULL element is returned.}}
+ \item\code{$data}{
+   An array with the same dimensions as the input \code{data} array, but with
+   spatial dimension reduced to the selected \code{region}}
+ \item\code{$lat}{A vector with the new corresponding latitudes for the 
+   selected \code{region}}
+ \item\code{$lon}{A vector with the new corresponding longitudes for the 
+   selected \code{region}}
+ \item\code{$mask}{If parameter \code{mask} is supplied, an array with 
+   reduced length of the dimensions to the selected \code{region}. Otherwise, 
+   a NULL element is returned.}}
 }
 \description{
-This function subsets an spatial region from spatial data giving a vector with the maximum and minimum of latitudes and longitudes of the selected region.
+Subset a spatial region from spatial data giving a vector with
+ the maximum and minimum of latitudes and longitudes of the selected region.
 }
 \examples{
-## Example with synthetic data:
+# Example with synthetic data:
 data <- 1:(20 * 3 * 2 * 4)
 dim(data) <- c(lon = 20, lat = 3, time = 2, model = 4)
 lon <- seq(2, 40, 2)
 lat <- c(1, 5, 10)
 
 a <- SelBox(data = data, lon = lon, lat = lat, region = c(2, 20, 1, 5), 
-           londim = 1, latdim = 2, mask = NULL)
-str(a)
+           londim = "lon", latdim = "lat", mask = NULL)
 }
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-Selbox.R b/tests/testthat/test-Selbox.R
new file mode 100644
index 0000000000000000000000000000000000000000..fce7f0690d0f866881f57c8aaefb1a167b84fc6b
--- /dev/null
+++ b/tests/testthat/test-Selbox.R
@@ -0,0 +1,201 @@
+context("ClimProjDiags::SelBox tests")
+
+##############################################
+
+# dat1
+dat1 <- 1:(2 * 3 * 2 * 4)
+dim(dat1) <- c(lon = 2, lat = 3, time = 2, model = 4)
+lon <- seq(2, 40, 20)
+lat <- c(1, 5, 10)
+region <-  c(1, 2, 1, 1)
+mask <- array(c(0, 1, 0, 1, 0, 1), dim = c(lon = 2, lat = 3))
+
+# dat2
+dat2 <- 1:(2 * 3 * 2 * 4)
+dim(dat2) <- c(lons = 2, lats = 3, time = 2, model = 4)
+mask2 <- array(c(0, 1, 0, 1, 0, 1), dim = c(lons = 2, lats = 3))
+
+# dat3
+dat3 <- 1:(20 * 3)
+dim(dat3) <- c(longitude = 20, latitude = 3)
+lon3 <- seq(2, 40, 2)
+lat3 <- c(1, 5, 10)
+region3 <- c(2, 20, 1, 5)
+region3_2 <- c(15, 5, 1, 2)
+mask3 <- array(c(0, 1, 0, 1, 0, 1), dim = c(latitude = 3, longitude = 20))
+
+##############################################
+
+test_that("1. Input checks", {
+  # data
+  expect_error(
+    SelBox(c(), lon = 1:10, lat = 1:10),
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    SelBox(1:10, lon = 1:10, lat = 1:10),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    SelBox(array(1:10), lon = 1:10, lat = 1:10),
+    "Parameter 'data' must have dimension names."
+  )
+  # lon, lat
+  expect_error(
+    SelBox(array(1:10, dim = c(lon = 10)), lon = NULL, lat = 1:10),
+    "Parameters 'lon' and 'lat' cannot be NULL."
+  )
+  expect_error(
+    SelBox(array(1:10, dim = c(lon = 10)), lon = 'a', lat = 1:10),
+    "Parameters 'lon' and 'lat' must be numeric."
+  )
+  expect_error(
+    SelBox(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."
+  )
+  # region
+  expect_error(
+    SelBox(dat1, lat = lat, lon = lon, region = NULL),
+    "Parameter 'region' cannot be NULL."
+  )
+  expect_error(
+    SelBox(dat1, lat = lat, lon = lon, region = 'a'),
+    "Parameter 'region' must be numeric."
+  )
+  expect_error(
+    SelBox(dat1, lat = lat, lon = lon, region = 1:3),
+    paste0("The region argument has to be a vector of length four indicating the ",
+           "minimum longitude, the maximum longitude, the minimum latitude and ", 
+           "the maximum latitude.")
+  )
+  # londim
+  expect_warning(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region,
+           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(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           lon = 1:10, lat = 1, londim = list(1,2), latdim = 2),
+    "Parameter 'londim' must be a character string."
+  )
+  expect_warning(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           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(
+    SelBox(array(1:10, dim = c(lons = 10, lat = 1)), region = region, 
+           lon = 1:10, lat = 1, londim = 'lon', latdim = 'lat'),
+    "Parameter 'londim' is not found in 'data'."
+  )
+  expect_error(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           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(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           lon = 1:10, lat = 1, latdim = 2),
+    "Numeric 'latdim' is deprecated, use dimension names instead. The ", 
+    "corresponding dimension name will be assigned."
+  )
+  expect_error(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           lon = 1:10, lat = 1, latdim = list(1,2)),
+    "Parameter 'latdim' must be a character string."
+  )
+  expect_warning(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           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(
+    SelBox(array(1:10, dim = c(lon = 10, lats = 1)), region = region, 
+           lon = 1:10, lat = 1, londim = 'lon', latdim = 'lat'),
+    "Parameter 'latdim' is not found in 'data'."
+  )
+  expect_error(
+    SelBox(array(1:10, dim = c(lon = 10, lat = 1)), region = region, 
+           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'.")
+  )
+  # mask
+  expect_error(
+    SelBox(dat1, lat = lat, lon = lon, region = region, mask = 1:3),
+    "Parameter 'mask' must have the same spatial dimensions of data."
+  )
+  expect_error(
+    SelBox(dat1, lat = lat, lon = lon, region = region, mask = mask2),
+    "Parameter 'mask' must have the same spatial dimensions of data."
+  )
+})
+
+##############################################
+
+test_that("2. Output checks", {
+  # dat2
+  expect_equal(
+    names(SelBox(dat2, lat = lat, lon = lon, region = region, 
+                 londim = 'lons', latdim = 'lats', mask = mask2)), 
+    c('data', 'lon', 'lat', 'mask')
+  )
+  expect_equal(
+    dim(SelBox(dat2, lat = lat, lon = lon, region = region, 
+               londim = 'lons', latdim = 'lats')$data), 
+    c(lons = 1, lats = 1, time = 2, model = 4)
+  )
+  expect_equal(
+    dim(SelBox(dat2, lat = lat, lon = lon, region = region, 
+               londim = 'lons', latdim = 'lats', mask = mask2)$mask), 
+    c(lons = 1, lats = 1)
+  )
+  expect_equal(
+    as.vector(SelBox(dat2, lat = lat, lon = lon, region = region, 
+              londim = 'lons', latdim = 'lats', mask = mask2)$data), 
+    c(1, 7, 13, 19, 25, 31, 37, 43)
+  )
+  expect_equal(
+    SelBox(dat2, lat = lat, lon = lon, region = region, 
+           londim = 'lons', latdim = 'lats', mask = mask2)$lat, 
+    1
+  )
+  # dat3
+  expect_equal(
+    dim(SelBox(dat3, lat = lat3, lon = lon3, region = region3, 
+               londim = 'longitude', latdim = 'latitude')$data), 
+    c(longitude = 10, latitude = 2)
+  )
+  expect_equal(
+    dim(SelBox(dat3, lat = lat3, lon = lon3, region = region3, 
+               londim = 'longitude', latdim = 'latitude', mask = mask3)$mask), 
+    c(latitude = 2, longitude = 10)
+  )
+  expect_equal(
+    SelBox(dat3, lat = lat3, lon = lon3, region = region3, 
+           londim = 'longitude', latdim = 'latitude', mask = mask3)$data, 
+    array(c(1:10, 21:30), dim = c(longitude = 10, latitude = 2))
+  )
+  expect_equal(
+    SelBox(dat3, lat = lat3, lon = lon3, region = region3, 
+           londim = 'longitude', latdim = 'latitude', mask = mask3)$lon, 
+    seq(2, 20, 2)
+  )
+  expect_equal(
+    SelBox(dat3, lat = lat3, lon = lon3, region = region3, 
+           londim = 'longitude', latdim = 'latitude', mask = mask3)$lat, 
+    c(1, 5)
+  )
+  expect_equal(
+    SelBox(dat3, lat = lat3, lon = lon3, region = region3_2,
+           londim = 'longitude', latdim = 'latitude', mask = mask3)$lon,
+    c(seq(16, 40, 2), 2, 4)
+  )
+})
diff --git a/tests/testthat/test-ShiftLon.R b/tests/testthat/test-ShiftLon.R
index efa358ebb90cb2ac44903dd28f12492fab7a4cd1..e414b5c9641fec7a05abd21cd980879459c77851 100644
--- a/tests/testthat/test-ShiftLon.R
+++ b/tests/testthat/test-ShiftLon.R
@@ -17,6 +17,11 @@ dat3 <- array(1:512, dim = c(lon = 512, lat = 216, time = 2))
 lon3 <- seq(0, 360, length.out = 513)[1:512]
 lat3 <- seq(-90, 90, length.out = 216)
 
+# dat4: vector, logical
+set.seed(1)
+dat4 <- as.logical(sample(c(0, 1), replace = TRUE, size = 360))
+lon4 <- 0:359
+
 ##############################################
 test_that("1. Input checks", {
   # data
@@ -24,14 +29,15 @@ test_that("1. Input checks", {
     ShiftLon(c()),
     "Parameter 'data' cannot be NULL."
   )
-  expect_error(
-    ShiftLon(1:10),
-    "Parameter 'data' must be a numeric array."
-  )
   expect_error(
    ShiftLon(array(1:10)),
    "Parameter 'data' must have dimension names."
   )
+  expect_warning(
+  ShiftLon(1:10, lon = 1:10, westB = 2),
+  "Parameter 'data' is a vector. Transfer it to an array and assign lon as dimension name."
+  )
+
   # lon_dim
   expect_error(
     ShiftLon(dat1, lon1, lon_dim = 1),
@@ -340,3 +346,15 @@ test_that("4. dat3", {
   )
   
 })
+
+
+##############################################
+test_that("5. dat4", {
+
+  westB <- -10
+expect_equal(
+suppressWarnings(as.vector(ShiftLon(dat4, lon4, westB)$data)),
+dat4[c(351:360, 1:350)]
+)
+
+})
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..3edfcaec471e3caadeddb9e5fd8e2cdb5e646daf 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
+  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.02119318, 0.01728634),
+    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.02119318, 0.01728634),
+    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
+  )
 })