diff --git a/DESCRIPTION b/DESCRIPTION
index cbcbdb22375b4d81c1b5604a71730bb9ea718172..bd980e3dceb2dd4949ff38d3324011944cd262d9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -27,7 +27,7 @@ License: Apache License 2.0
 URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
 BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/-/issues
 Encoding: UTF-8
-RoxygenNote: 5.0.0
+RoxygenNote: 7.2.0
 Suggests:
     knitr,
     testthat,
diff --git a/R/Subset.R b/R/Subset.R
index bad89e7b4a200e2e1591155141a177f629dd6b1a..c077c2a96330f6fd9d236ac6c037493fb5a4a42c 100644
--- a/R/Subset.R
+++ b/R/Subset.R
@@ -1,22 +1,44 @@
 #'Subset a Data Array
 #'
-#'This function allows to subset (i.e. slice, take a chunk of) an array, in a similar way as done in the function \code{take()} in the package plyr. There are two main inprovements:\cr\cr The input array can have dimension names, either in \code{names(dim(x))} or in the attribute 'dimensions', and the dimensions to subset along can be specified via the parameter \code{along} either with integer indices or either by their name.\cr\cr There are additional ways to adjust which dimensions are dropped in the resulting array: either to drop all, to drop none, to drop only the ones that have been sliced or to drop only the ones that have not been sliced.\cr\cr If an array is provided without dimension names, dimension names taken from the parameter \code{dim_names} will be added to the array.
-#'This function computes the threshold based on a quantile value for each day of the year of the daily data input.
+#'This function allows to subset (i.e. slice, take a chunk of) an array, in a 
+#'similar way as done in the function \code{take()} in the package plyr. There
+#'are two main snprovements:\cr\cr First, the input array can have dimension 
+#'names, either in \code{names(dim(x))} or in the attribute 'dimensions'. If 
+#'both exist, the attribute 'dimensions' is prioritized. The dimensions to 
+#'subset along can be specified via the parameter \code{along} either with
+#'integer indices or either by their name.\cr\cr Second, there are additional
+#'ways to adjust which dimensions are dropped in the resulting array: either to
+#'drop all, to drop none, to drop only the ones that have been sliced or to drop
+#'only the ones that have not been sliced.\cr\cr 
 #'
-#'@param x A multidimensional array to be sliced. It can have dimension names either in \code{names(dim(x))} or either in the attribute 'dimensions'.
-#'@param along Vector with references to the dimensions to take the subset from: either integers or dimension names.
-#'@param indices List of indices to take from each dimension specified in 'along'. If a single dimension is specified in 'along' the indices can be directly provided as a single integer or as a vector.
-#'@param drop Whether to drop all the dimensions of length 1 in the resulting array, none, only those that are specified in 'along', or only those that are not specified in 'along'. The possible values are, respectively: 'all' or TRUE, 'none' or FALSE, 'selected', and 'non-selected'. 
+#'@param x A named multidimensional array to be sliced. It can have dimension 
+#'  names either in \code{names(dim(x))} or in the attribute 'dimensions'.
+#'@param along A vector with references to the dimensions to take the subset 
+#'  from: either integers or dimension names.
+#'@param indices A list of indices to take from each dimension specified in 
+#'  'along'. If a single dimension is specified in 'along', it can be directly
+#'  provided as an integer or a vector.
+#'@param drop Whether to drop all the dimensions of length 1 in the resulting 
+#'  array, none, only those that are specified in 'along', or only those that 
+#'  are not specified in 'along'. The possible values are: 'all' or TRUE, 'none'
+#'  or FALSE, 'selected', and 'non-selected'. The default value is FALSE.
 #'
-#'@return An array with similar dimensions as the \code{x} input, but with trimmed or dropped dimensions.
+#'@return An array with similar dimensions as the \code{x} input, but with 
+#'  trimmed or dropped dimensions.
 #'
 #'@examples
-#'##Example synthetic data:
+#'#Example synthetic data:
+#'# Dimension has name already
 #'data <- 1:(2 * 3 * 372 * 1)
 #'dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
 #'data_subset <- Subset(data, c('time', 'model'), 
 #'                      list(1:10, TRUE), drop = 'selected')
 #'dim(data_subset)
+#'# Use attributes 'dimensions'
+#'data <- array(1:(2 * 3 * 372 * 1), dim = c(2, 3, 372, 1))
+#'attributes(data)[['dimensions']] <- c('lat', 'lon', 'time', 'model')
+#'data_subset <- Subset(data, c('lon', 'lat'), list(1, 1), drop = TRUE)
+#'dim(data_subset)
 #'
 #'@export
 Subset <- function(x, along, indices, drop = FALSE) {
@@ -29,7 +51,10 @@ Subset <- function(x, along, indices, drop = FALSE) {
   dim_names <- attr(x, 'dimensions')
   if (!is.character(dim_names)) {
     dim_names <- names(dim(x))
+  } else {
+    names(dim(x)) <- dim_names
   }
+
   if (!is.character(dim_names)) {
     if (any(sapply(along, is.character))) {
       stop("The input array 'x' doesn't have labels for the dimensions but the parameter 'along' contains dimension names.")
@@ -54,7 +79,14 @@ Subset <- function(x, along, indices, drop = FALSE) {
 
   # Check indices
   if (!is.list(indices)) {
-    indices <- list(indices)
+    if (length(along) == 1) {
+      indices <- list(indices)
+    } else {
+      stop("Parameter 'indices' should be a list.")
+    }
+  }
+  if (length(indices) != length(along)) {
+    stop("Parameter 'along' and 'indices' should have the same length.")
   }
 
   # Check parameter drop
@@ -90,25 +122,42 @@ Subset <- function(x, along, indices, drop = FALSE) {
     if (length(dim_names_to_remove) > 0) {
       dim_names <- dim_names[-dim_names_to_remove]
     }
+    # If there is one dim left, subset won't have dimension (but it should have one). Add it back
+    if (is.null(dim(subset))) {
+      if (!identical(dim_names, character(0))) {
+        # If there is one dim left, subset won't have dimension (but it should 
+        # have one). Add it back.
+        dim(subset) <- dim(x)[dim_names]
+      } else {  # a number left
+        dim(subset) <- 1
+      }
+    }
   }
 
   # Amend the final dimensions and put dimnames and attributes
   metadata <- attributes(x)
   metadata[['dim']] <- dim(subset)
   if (length(dims_to_drop) > 0) {
-    metadata[['dim']] <- metadata[['dim']][-dims_to_drop]
+    if (length(dims_to_drop) == length(metadata[['dim']])) { # a number left
+      metadata[['dim']] <- 1
+    } else {
+      metadata[['dim']] <- metadata[['dim']][-dims_to_drop]
+    }
     if (is.character(dim_names)) {
-      names(metadata[['dim']]) <- dim_names[-dims_to_drop]
+      if (!identical(dim_names[-dims_to_drop], character(0))) {
+        names(metadata[['dim']]) <- dim_names[-dims_to_drop]
+      }
       if ('dimensions' %in% names(attributes(x))) {
         metadata[['dimensions']] <- dim_names[-dims_to_drop]
       }
     }
-  } else if (is.character(dim_names)) {
+  } else if (is.character(dim_names) & !identical(dim_names, character(0))) {
     names(metadata[['dim']]) <- dim_names
     if ('dimensions' %in% names(attributes(x))) {
       metadata[['dimensions']] <- dim_names
     }
   }
+
   attributes(subset) <- metadata
   subset
 }
diff --git a/man/AnoAgree.Rd b/man/AnoAgree.Rd
index c929352d71c27790644570ed71b937a22a7072d9..63c1450039013f94bee9c4e66d29c0923fcdab84 100644
--- a/man/AnoAgree.Rd
+++ b/man/AnoAgree.Rd
@@ -34,4 +34,3 @@ a <- rnorm(6)
 agree <- AnoAgree(ano = a, membersdim = 1, na.rm = TRUE, ncores = NULL)
 print(agree)
 }
-
diff --git a/man/ArrayToList.Rd b/man/ArrayToList.Rd
index 9a7d181f1926cea32625c95a6861e0deaaaff7f0..9951ae13d42a263ce58fdfd20d954c8d23f84ed2 100644
--- a/man/ArrayToList.Rd
+++ b/man/ArrayToList.Rd
@@ -38,4 +38,3 @@ str(datalist)
 \seealso{
 \link[s2dv]{PlotLayout}
 }
-
diff --git a/man/Climdex.Rd b/man/Climdex.Rd
index 669570ab98f8060630d9df2ec90911deec345a1b..372712c240998579f23b70b8fbfd998656ce00f7 100644
--- a/man/Climdex.Rd
+++ b/man/Climdex.Rd
@@ -4,8 +4,16 @@
 \alias{Climdex}
 \title{Wrapper for applying the climdex routine  ETCCDI climate change indices to n-dimensional arrays.}
 \usage{
-Climdex(data, metric, threshold = NULL, base.range = NULL, dates = NULL,
-  timedim = NULL, calendar = NULL, ncores = NULL)
+Climdex(
+  data,
+  metric,
+  threshold = NULL,
+  base.range = NULL,
+  dates = NULL,
+  timedim = NULL,
+  calendar = NULL,
+  ncores = NULL
+)
 }
 \arguments{
 \item{data}{A numeric n-dimensional array containing daily maximum or minimum temperature, wind speed or precipitation amount.}
@@ -68,4 +76,3 @@ David Bronaugh for the Pacific Climate Impacts Consortium (2015).
  climdex.pcic: PCIC Implementation of Climdex Routines. R package
  version 1.1-6. http://CRAN.R-project.org/package=climdex.pcic
 }
-
diff --git a/man/CombineIndices.Rd b/man/CombineIndices.Rd
index 95abc5ed4f5564527b93bd9646221b7c4b1a8d3d..6d032cddcba684a8375572aa6c259da2c9d5e982 100644
--- a/man/CombineIndices.Rd
+++ b/man/CombineIndices.Rd
@@ -33,4 +33,3 @@ dim(b) <- c(lon = 2, lat = 3, mod = 4)
 comb_ind <- CombineIndices(indices = list(a, b), weights = c(2, 1), operation = "add")
 print(comb_ind)
 }
-
diff --git a/man/DTRIndicator.Rd b/man/DTRIndicator.Rd
index 33e7f507882ba4583c5796645b05991a1d93af22..e8be18bb79aa379757334633414c7ea5438dec0f 100644
--- a/man/DTRIndicator.Rd
+++ b/man/DTRIndicator.Rd
@@ -4,8 +4,16 @@
 \alias{DTRIndicator}
 \title{Diurnal temperature range indicator (DTR) of multidimensional arrays}
 \usage{
-DTRIndicator(tmax, tmin, ref, by.seasons = TRUE, dates = NULL,
-  timedim = NULL, calendar = NULL, ncores = NULL)
+DTRIndicator(
+  tmax,
+  tmin,
+  ref,
+  by.seasons = TRUE,
+  dates = NULL,
+  timedim = NULL,
+  calendar = NULL,
+  ncores = NULL
+)
 }
 \arguments{
 \item{tmax}{A numeric multidimensional array containing daily maximum temperature.}
@@ -60,4 +68,3 @@ aa <- DTRIndicator(tmax, tmin, ref = a, by.seasons = FALSE, ncores = NULL)
 str(aa)
 dim(aa$indicator)
 }
-
diff --git a/man/DTRRef.Rd b/man/DTRRef.Rd
index 6e32e31a0085e87305b3e30c601dec7fe28fa6bc..daea8c2390d609adb7055b306d56ad7b870d846a 100644
--- a/man/DTRRef.Rd
+++ b/man/DTRRef.Rd
@@ -4,8 +4,16 @@
 \alias{DTRRef}
 \title{Diurnal temperature range of multidimensional arrays}
 \usage{
-DTRRef(tmax, tmin, by.seasons = TRUE, dates = NULL, timedim = NULL,
-  calendar = NULL, na.rm = TRUE, ncores = NULL)
+DTRRef(
+  tmax,
+  tmin,
+  by.seasons = TRUE,
+  dates = NULL,
+  timedim = NULL,
+  calendar = NULL,
+  na.rm = TRUE,
+  ncores = NULL
+)
 }
 \arguments{
 \item{tmax}{A numeric multidimensional array containing daily maximum temperature.}
@@ -68,4 +76,3 @@ dim(tmin) <- c(2, 3, 365)
 a <- DTRRef(tmax, tmin, by.seasons = FALSE, dates = time,  timedim = 3, ncores = NULL)
 str(a)
 }
-
diff --git a/man/DailyAno.Rd b/man/DailyAno.Rd
index d4033f79e505fbd3df46a6304be677212c811f1b..3bf59f6cd8941a724a7723c8ae6f16d72394e904 100644
--- a/man/DailyAno.Rd
+++ b/man/DailyAno.Rd
@@ -4,8 +4,7 @@
 \alias{DailyAno}
 \title{Daily anomalies}
 \usage{
-DailyAno(data, jdays = NULL, dates = NULL, calendar = NULL,
-  na.rm = TRUE)
+DailyAno(data, jdays = NULL, dates = NULL, calendar = NULL, na.rm = TRUE)
 }
 \arguments{
 \item{data}{A vector of daily data.}
@@ -31,4 +30,3 @@ jdays <- c(rep(1, 5), rep(2, 5))
 daily_anomaly <- DailyAno(data = data, jdays = jdays, na.rm = TRUE)
 print(daily_anomaly)
 }
-
diff --git a/man/Extremes.Rd b/man/Extremes.Rd
index 0939fa67275433b18b9f62df84150911ea72e352..ddc57e8cab54aafb7acd93620d5bc51e7b512699 100644
--- a/man/Extremes.Rd
+++ b/man/Extremes.Rd
@@ -4,9 +4,18 @@
 \alias{Extremes}
 \title{Sum of spell lengths exceeding daily threshold for n-dimensional arrays}
 \usage{
-Extremes(data, threshold, op = ">", min.length = 6,
-  spells.can.span.years = TRUE, max.missing.days = 5, dates = NULL,
-  timedim = NULL, calendar = NULL, ncores = NULL)
+Extremes(
+  data,
+  threshold,
+  op = ">",
+  min.length = 6,
+  spells.can.span.years = TRUE,
+  max.missing.days = 5,
+  dates = NULL,
+  timedim = NULL,
+  calendar = NULL,
+  ncores = NULL
+)
 }
 \arguments{
 \item{data}{A n-dimensional array containing daily data.}
@@ -58,4 +67,3 @@ a <- Extremes(data, threshold = threshold, op = ">", min.length = 6, spells.can.
              max.missing.days = 5, ncores = NULL)
 str(a)
 }
-
diff --git a/man/Lon2Index.Rd b/man/Lon2Index.Rd
index b2012a61f97126ee886ba942f3e5020961d29994..dbae6dcbc533716959433999b059eba7276c1799 100644
--- a/man/Lon2Index.Rd
+++ b/man/Lon2Index.Rd
@@ -33,4 +33,3 @@ pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
 lon[pos]
 
 }
-
diff --git a/man/SeasonSelect.Rd b/man/SeasonSelect.Rd
index a71fd242c9a71d0647fa10fad5fc1a440823af70..33066eed0e4ce2cd5b1a1abbe92930b43854975c 100644
--- a/man/SeasonSelect.Rd
+++ b/man/SeasonSelect.Rd
@@ -45,4 +45,3 @@ attr(data, 'Variables')$common[[2]]$dim[[3]]$vals <- time
 a <- SeasonSelect(data = data, season = 'JJA')
 str(a)
 }
-
diff --git a/man/SelBox.Rd b/man/SelBox.Rd
index 07e92b20e28a41e872059f311378f8e1f90b7573..38e3547bf6ae6d6e699e78821aaefb1311896535 100644
--- a/man/SelBox.Rd
+++ b/man/SelBox.Rd
@@ -43,4 +43,3 @@ a <- SelBox(data = data, lon = lon, lat = lat, region = c(2, 20, 1, 5),
            londim = 1, latdim = 2, mask = NULL)
 str(a)
 }
-
diff --git a/man/Subset.Rd b/man/Subset.Rd
index 634cfe2a60d07b8d31eeafa4885fca01afd71b42..f50c5da0bb4098fb407a0335904a70147fa47c70 100644
--- a/man/Subset.Rd
+++ b/man/Subset.Rd
@@ -7,28 +7,49 @@
 Subset(x, along, indices, drop = FALSE)
 }
 \arguments{
-\item{x}{A multidimensional array to be sliced. It can have dimension names either in \code{names(dim(x))} or either in the attribute 'dimensions'.}
+\item{x}{A named multidimensional array to be sliced. It can have dimension 
+names either in \code{names(dim(x))} or in the attribute 'dimensions'.}
 
-\item{along}{Vector with references to the dimensions to take the subset from: either integers or dimension names.}
+\item{along}{A vector with references to the dimensions to take the subset 
+from: either integers or dimension names.}
 
-\item{indices}{List of indices to take from each dimension specified in 'along'. If a single dimension is specified in 'along' the indices can be directly provided as a single integer or as a vector.}
+\item{indices}{A list of indices to take from each dimension specified in 
+'along'. If a single dimension is specified in 'along', it can be directly
+provided as an integer or a vector.}
 
-\item{drop}{Whether to drop all the dimensions of length 1 in the resulting array, none, only those that are specified in 'along', or only those that are not specified in 'along'. The possible values are, respectively: 'all' or TRUE, 'none' or FALSE, 'selected', and 'non-selected'.}
+\item{drop}{Whether to drop all the dimensions of length 1 in the resulting 
+array, none, only those that are specified in 'along', or only those that 
+are not specified in 'along'. The possible values are: 'all' or TRUE, 'none'
+or FALSE, 'selected', and 'non-selected'. The default value is FALSE.}
 }
 \value{
-An array with similar dimensions as the \code{x} input, but with trimmed or dropped dimensions.
+An array with similar dimensions as the \code{x} input, but with 
+ trimmed or dropped dimensions.
 }
 \description{
-This function allows to subset (i.e. slice, take a chunk of) an array, in a similar way as done in the function \code{take()} in the package plyr. There are two main inprovements:\cr\cr The input array can have dimension names, either in \code{names(dim(x))} or in the attribute 'dimensions', and the dimensions to subset along can be specified via the parameter \code{along} either with integer indices or either by their name.\cr\cr There are additional ways to adjust which dimensions are dropped in the resulting array: either to drop all, to drop none, to drop only the ones that have been sliced or to drop only the ones that have not been sliced.\cr\cr If an array is provided without dimension names, dimension names taken from the parameter \code{dim_names} will be added to the array.
-This function computes the threshold based on a quantile value for each day of the year of the daily data input.
+This function allows to subset (i.e. slice, take a chunk of) an array, in a 
+similar way as done in the function \code{take()} in the package plyr. There
+are two main snprovements:\cr\cr First, the input array can have dimension 
+names, either in \code{names(dim(x))} or in the attribute 'dimensions'. If 
+both exist, the attribute 'dimensions' is prioritized. The dimensions to 
+subset along can be specified via the parameter \code{along} either with
+integer indices or either by their name.\cr\cr Second, there are additional
+ways to adjust which dimensions are dropped in the resulting array: either to
+drop all, to drop none, to drop only the ones that have been sliced or to drop
+only the ones that have not been sliced.\cr\cr
 }
 \examples{
-##Example synthetic data:
+#Example synthetic data:
+# Dimension has name already
 data <- 1:(2 * 3 * 372 * 1)
 dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
 data_subset <- Subset(data, c('time', 'model'), 
                      list(1:10, TRUE), drop = 'selected')
 dim(data_subset)
+# Use attributes 'dimensions'
+data <- array(1:(2 * 3 * 372 * 1), dim = c(2, 3, 372, 1))
+attributes(data)[['dimensions']] <- c('lat', 'lon', 'time', 'model')
+data_subset <- Subset(data, c('lon', 'lat'), list(1, 1), drop = TRUE)
+dim(data_subset)
 
 }
-
diff --git a/man/Threshold.Rd b/man/Threshold.Rd
index cec12e68bc83ea93767f034c6aa417b5076e81b0..a0fa10aec9ff6aab2e2993a9cbfb2be4338f309b 100644
--- a/man/Threshold.Rd
+++ b/man/Threshold.Rd
@@ -4,8 +4,15 @@
 \alias{Threshold}
 \title{Daily thresholds based on quantiles for n-dimensional arrays}
 \usage{
-Threshold(data, dates = NULL, calendar = NULL, base.range = NULL,
-  qtiles = 0.9, ncores = NULL, na.rm = FALSE)
+Threshold(
+  data,
+  dates = NULL,
+  calendar = NULL,
+  base.range = NULL,
+  qtiles = 0.9,
+  ncores = NULL,
+  na.rm = FALSE
+)
 }
 \arguments{
 \item{data}{A numeric n-dimensional array containing daily data.}
@@ -42,4 +49,3 @@ attr(data, 'Variables')$dat1$time <- time
 a <- Threshold(data, dates = NULL, base.range = NULL, qtiles = 0.9, ncores = NULL)
 str(a)
 }
-
diff --git a/man/WaveDuration.Rd b/man/WaveDuration.Rd
index b0682523c96b1d604c09c6c397b8040523c7e4e7..a43c6b9252de6976d55a2db1096a5b971c375eb5 100644
--- a/man/WaveDuration.Rd
+++ b/man/WaveDuration.Rd
@@ -4,8 +4,16 @@
 \alias{WaveDuration}
 \title{Heat and cold waves duration for n-dimensional arrays}
 \usage{
-WaveDuration(data, threshold, op = ">", spell.length = 6,
-  by.seasons = TRUE, dates = NULL, calendar = NULL, ncores = NULL)
+WaveDuration(
+  data,
+  threshold,
+  op = ">",
+  spell.length = 6,
+  by.seasons = TRUE,
+  dates = NULL,
+  calendar = NULL,
+  ncores = NULL
+)
 }
 \arguments{
 \item{data}{A numeric n-dimensional array containing daily maximum or minimum temperature}
@@ -48,4 +56,3 @@ threshold <- rep(40, 31)
 a <- WaveDuration(data, threshold, op = ">", spell.length = 6, by.seasons = TRUE, ncores = NULL)
 str(a)
 }
-
diff --git a/man/WeightedMean.Rd b/man/WeightedMean.Rd
index 23f31c6c415b83ecd829fba0cc27c3c82d750504..ac3411cba7fc3fadea3f557220d7a6607aedad12 100644
--- a/man/WeightedMean.Rd
+++ b/man/WeightedMean.Rd
@@ -4,8 +4,15 @@
 \alias{WeightedMean}
 \title{Calculate spatial area-weighted average of multidimensional arrays}
 \usage{
-WeightedMean(data, lon, lat, region = NULL, mask = NULL, londim = NULL,
-  latdim = NULL)
+WeightedMean(
+  data,
+  lon,
+  lat,
+  region = NULL,
+  mask = NULL,
+  londim = NULL,
+  latdim = NULL
+)
 }
 \arguments{
 \item{data}{An array with minimum two dimensions of latitude and longitude.}
@@ -60,4 +67,3 @@ a <- WeightedMean(data = data, lon = lon, lat = lat, region = NULL,
                  mask = NULL, londim = 1, latdim = 2)
 str(a)
 }
-
diff --git a/tests/testthat/test-Subset.R b/tests/testthat/test-Subset.R
new file mode 100644
index 0000000000000000000000000000000000000000..42f609a14bb366240d77ec1b50c5480aaeeafe55
--- /dev/null
+++ b/tests/testthat/test-Subset.R
@@ -0,0 +1,286 @@
+context("Subset tests")
+
+dat1 <- array(1:20, dim = c(dat = 1, lat = 2, lon = 10))
+
+dat2 <- dat1
+attributes(dat2)[['units']] <- 'K'
+
+dat3 <- array(1:20, dim = c(1, 2, 10))
+attributes(dat3)[['dimensions']] <- c('dat', 'lat', 'lon')
+
+dat4 <- dat1
+attributes(dat4)[['dimensions']] <- c('dataset', 'latitude', 'longitude')
+
+#==============================================
+
+test_that("1. Sanity checks", {
+
+expect_error(
+Subset(x = 1:4),
+"Input array 'x' must be a numeric array."
+)
+expect_error(
+Subset(x = array(1:20, dim = c(1, 2, 10)), along = 'a'),
+"The input array 'x' doesn't have labels for the dimensions but the parameter 'along' contains dimension names."
+)
+expect_error(
+Subset(x = dat1, along = T),
+"All provided dimension indices in 'along' must be integers or character strings."
+)
+expect_error(
+Subset(x = dat1, along = c('dat', 'dat')),
+"The parameter 'along' must not contain repeated dimension names."
+)
+expect_error(
+Subset(x = dat1, along = 'dataset'),
+"Could not match all dimension names in 'indices' with dimension names in input array 'x'."
+)
+expect_error(
+Subset(x = dat1, 'dat', 1, drop = 'yes'),
+"Parameter 'drop' must be one of TRUE, FALSE, 'all', 'selected', 'non-selected', 'none'."
+)
+
+})
+
+test_that("2. dat1", {
+# drop: dat
+expect_equal(
+dim(Subset(dat1, 'dat', 1, drop = FALSE)),
+dim(dat1)
+)
+expect_equal(
+as.vector(Subset(dat1, 'dat', 1, drop = FALSE)),
+as.vector(dat1)
+)
+expect_equal(
+attributes(Subset(dat1, 'dat', 1, drop = FALSE)),
+attributes(dat1)
+)
+expect_equal(
+dim(Subset(dat1, 'dat', 1, drop = TRUE)),
+dim(dat1)[-1]
+)
+expect_equal(
+as.vector(Subset(dat1, 'dat', 1, drop = TRUE)),
+as.vector(dat1)
+)
+expect_equal(
+as.vector(Subset(dat1, 'dat', 1, drop = TRUE)),
+as.vector(Subset(dat1, 1, 1, drop = TRUE))
+)
+expect_equal(
+dim(Subset(dat1, 'dat', 1, drop = 'selected')),
+dim(dat1)[-1]
+)
+expect_equal(
+as.vector(Subset(dat1, 'dat', 1, drop = 'selected')),
+as.vector(dat1)
+)
+# drop: lat
+expect_equal(
+dim(Subset(dat1, 'lat', 1, drop = FALSE)),
+c(dat = 1, lat = 1, lon = 10)
+)
+expect_equal(
+as.vector(Subset(dat1, 'lat', 1, drop = FALSE)),
+as.vector(dat1[, 1, ])
+)
+expect_equal(
+dim(Subset(dat1, 'lat', 1, drop = TRUE)),
+dim(dat1)[3]
+)
+expect_equal(
+as.vector(Subset(dat1, 'lat', 1, drop = TRUE)),
+as.vector(dat1[, 1, ])
+)
+expect_equal(
+Subset(dat1, 'lat', 1, drop = TRUE),
+Subset(dat1, 2, 1, drop = TRUE)
+)
+expect_equal(
+dim(Subset(dat1, 'lat', 1, drop = 'selected')),
+dim(dat1)[c(1, 3)]
+)
+expect_equal(
+as.vector(Subset(dat1, 'lat', 1, drop = 'selected')),
+as.vector(dat1[, 1, ])
+)
+# drop: lat, lon
+expect_equal(
+dim(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = FALSE)),
+c(dat = 1, lat = 1, lon = 1)
+)
+expect_equal(
+as.vector(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = FALSE)),
+3
+)
+expect_equal(
+dim(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = TRUE)),
+1
+)
+expect_equal(
+as.vector(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = TRUE)),
+3
+)
+expect_equal(
+Subset(dat1, c('lat', 'lon'), list(1, 2), drop = TRUE),
+Subset(dat1, c(2, 3), list(1, 2), drop = TRUE)
+)
+expect_equal(
+dim(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = 'selected')),
+c(dat = 1)
+)
+expect_equal(
+as.vector(Subset(dat1, c('lat', 'lon'), list(1, 2), drop = 'selected')),
+3
+)
+# drop: dat, lat, lon
+expect_equal(
+dim(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE)),
+c(dat = 1, lat = 1, lon = 1)
+)
+expect_equal(
+as.vector(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE)),
+3
+)
+expect_equal(
+Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE),
+Subset(dat1, c(1:3), list(1, 1, 2), drop = FALSE)
+)
+expect_equal(
+dim(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE)),
+1
+)
+expect_equal(
+as.vector(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE)),
+3
+)
+expect_equal(
+Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE),
+Subset(dat1, c(1:3), list(1, 1, 2), drop = TRUE)
+)
+expect_equal(
+class(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+'array'
+)
+expect_equal(
+dim(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+1
+)
+expect_equal(
+as.vector(Subset(dat1, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+3
+)
+
+})
+
+test_that("3. dat2", {
+# Same as dat1 but with attributes
+
+# drop: lat
+expect_equal(
+dim(Subset(dat2, 'lat', 1, drop = TRUE)),
+dim(dat1)[3]
+)
+expect_equal(
+as.vector(Subset(dat2, 'lat', 1, drop = TRUE)),
+as.vector(dat1[, 1, ])
+)
+expect_equal(
+attributes(Subset(dat2, 'lat', 1, drop = TRUE)),
+list(dim = c(lon = 10), units = 'K')
+)
+
+# drop: dat, lat, lon
+expect_equal(
+dim(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE)),
+c(dat = 1, lat = 1, lon = 1)
+)
+expect_equal(
+as.vector(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE)),
+3
+)
+expect_equal(
+attributes(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = FALSE)),
+list(dim = c(dat = 1, lat = 1, lon = 1), units = 'K')
+)
+expect_equal(
+dim(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE)),
+1
+)
+expect_equal(
+as.vector(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE)),
+3
+)
+expect_equal(
+attributes(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = TRUE)),
+list(dim = 1, units = 'K')
+)
+expect_equal(
+class(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+'array'
+)
+expect_equal(
+dim(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+1
+)
+expect_equal(
+as.vector(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+3
+)
+expect_equal(
+attributes(Subset(dat2, c('dat', 'lat', 'lon'), list(1, 1, 2), drop = 'selected')),
+list(dim = 1, units = 'K')
+)
+
+})
+
+test_that("4. dat3", {
+# drop: dat
+expect_equal(
+Subset(dat1, 'dat', 1, drop = FALSE),
+Subset(dat3, 'dat', 1, drop = FALSE),
+check.attributes = F
+)
+expect_equal(
+attributes(Subset(dat3, 'dat', 1, drop = FALSE)),
+list(dim = c(dat = 1, lat = 2, lon = 10), dimensions = c('dat', 'lat', 'lon'))
+)
+
+# drop: lat, lon
+expect_equal(
+Subset(dat1, c('lat', 'lon'), list(1, 2), drop = TRUE),
+Subset(dat3, c('lat', 'lon'), list(1, 2), drop = TRUE),
+check.attributes = F
+)
+expect_equal(
+attributes(Subset(dat3, c('lat', 'lon'), list(1, 2), drop = TRUE)),
+list(dim = 1, dimensions = c('dat', 'lat', 'lon'))
+)
+
+})
+
+test_that("5. dat4", {
+
+# drop: lat
+expect_equal(
+Subset(dat4, 'latitude', 1, drop = FALSE),
+Subset(dat1, 'lat', 1, drop = FALSE),
+check.attributes = F
+)
+expect_equal(
+attributes(Subset(dat4, 'latitude', 1, drop = FALSE)),
+list(dim = c(dataset = 1, latitude = 1, longitude = 10), dimensions = c('dataset', 'latitude', 'longitude'))
+)
+# drop: lat, lon
+expect_equal(
+Subset(dat4, c('latitude', 'longitude'), list(1, 2), drop = TRUE),
+Subset(dat1, c('lat', 'lon'), list(1, 2), drop = TRUE),
+check.attributes = F
+)
+expect_equal(
+attributes(Subset(dat4, c('latitude', 'longitude'), list(1, 2), drop = TRUE)),
+list(dim = 1, dimensions = c('dataset', 'latitude', 'longitude'))
+)
+
+})