diff --git a/DESCRIPTION b/DESCRIPTION
index 4151ba48204d6c18077b6610edda77ed8ad2c0a8..835c96b1fb7f0fbb8d34997b911fc4e9942196e2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -27,7 +27,8 @@ Depends:
     R (>= 3.6.0)
 Imports:
     multiApply (>= 2.1.1),
-    stats
+    stats, 
+    ClimProjDiags
 Suggests:
     testthat,
     CSTools,
@@ -40,4 +41,4 @@ URL: https://earth.bsc.es/gitlab/es/csindicators/
 BugReports: https://earth.bsc.es/gitlab/es/csindicators/-/issues
 Encoding: UTF-8
 RoxygenNote: 7.2.0
-Config/testthat/edition: 3
\ No newline at end of file
+Config/testthat/edition: 3
diff --git a/NAMESPACE b/NAMESPACE
index d80accbbe3b6781acb27cb8c953a9c6b00a1824a..0a16d4d2dfb5888adf0007be1556694a7f87b899 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,6 +26,7 @@ export(TotalTimeExceedingThreshold)
 export(WindCapacityFactor)
 export(WindPowerDensity)
 import(multiApply)
+importFrom(ClimProjDiags,Subset)
 importFrom(stats,approxfun)
 importFrom(stats,ecdf)
 importFrom(stats,quantile)
diff --git a/R/AccumulationExceedingThreshold.R b/R/AccumulationExceedingThreshold.R
index e346b5310504a9e0a9b0fbbf289d332d4cbeb8cb..083b76cf07e3fdcc3d44246676b96b18e2b3fb14 100644
--- a/R/AccumulationExceedingThreshold.R
+++ b/R/AccumulationExceedingThreshold.R
@@ -49,15 +49,32 @@
 #'
 #'@return An 's2dv_cube' object containing the aggregated values in the element
 #'\code{data} with dimensions of the input parameter 'data' except the dimension
-#'where the indicator has been computed.
+#'where the indicator has been computed. The 'Dates' array is updated to 
+#'the dates corresponding to the beginning of the aggregated time period. A new 
+#'element called 'time_bounds' will be added into the 'attrs' element in the 
+#''s2dv_cube' object. It consists of a list containing two elements, the start 
+#'and end dates of the aggregated period with the same dimensions of 'Dates' 
+#'element.
+#' 
 #'@examples 
 #'exp <- NULL
-#'exp$data <- array(rnorm(216)*200, dim = c(dataset = 1, member = 2, sdate = 3, 
-#'                  ftime = 9, lat = 2, lon = 2))
+#'exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*100), 
+#'                  c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
 #'class(exp) <- 's2dv_cube'
-#'DOT <- CST_AccumulationExceedingThreshold(exp, threshold = 280)
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'exp$attrs$Dates <- Dates
+#'AT <- CST_AccumulationExceedingThreshold(data = exp, threshold = 100,
+#'                                         start = list(21, 4), 
+#'                                         end = list(21, 6))
 #'
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FALSE,
                                                start = NULL, end = NULL, time_dim = 'ftime',
@@ -78,7 +95,7 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
 
   if (length(op) == 1) {
     if (inherits(threshold, 's2dv_cube')) {
-        threshold <- threshold$data
+      threshold <- threshold$data
     }
   } else if (length(op) == 2) {
     if (inherits(threshold[[1]], 's2dv_cube')) {
@@ -89,17 +106,38 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
     }
   }
 
-  total <- AccumulationExceedingThreshold(data$data, dates = data$attrs$Dates,
+  Dates <- data$attrs$Dates
+  total <- AccumulationExceedingThreshold(data = data$data, dates = Dates,
                                           threshold = threshold, op = op, diff = diff,
                                           start = start, end = end, time_dim = time_dim,
                                           na.rm = na.rm, ncores = ncores)
   data$data <- total
-  if (!is.null(start) && !is.null(end)) {
-    data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
-                                            start = start, end = end,
-                                            time_dim = time_dim, 
-                                            ncores = ncores)
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
   }
+
   return(data)
 }
 #'Accumulation of a variable when Exceeding (not exceeding) a Threshold
@@ -133,9 +171,9 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
 #'@param diff A logical value indicating whether to accumulate the difference 
 #'  between data and threshold (TRUE) or not (FALSE by default). It can only be
 #'  TRUE if a unique threshold is used.
-#'@param dates A vector of dates or a multidimensional array with of dates with 
-#'  named dimensions matching the dimensions on parameter 'data'. By default it
-#'  is NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to define the initial date of the period
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is 
@@ -161,12 +199,6 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
 #'# Assuming data is already (tasmax + tasmin)/2 - 10
 #'data <- array(rnorm(5 * 3 * 214 * 2, mean = 25, sd = 3),
 #'                    c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
-#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
-#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
-#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
-#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
-#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
-#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
 #'GDD <- AccumulationExceedingThreshold(data, threshold = 0, start = list(1, 4),
 #'                                      end = list(31, 10))
 #'@import multiApply
@@ -228,8 +260,8 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
   }
   if (length(op) == 2) {
     if (length(op) != length(threshold)) {
-      stop(paste0("If 'op' is a  pair of logical operators parameter 'threshold' ",
-                  "also has to be a pair of values."))
+      stop("If 'op' is a  pair of logical operators parameter 'threshold' ",
+           "also has to be a pair of values.")
     }
     if (!is.numeric(threshold[[1]]) | !is.numeric(threshold[[2]])) {
       stop("Parameter 'threshold' must be numeric.")
@@ -240,7 +272,8 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
 
     if (!is.array(threshold[[1]]) && length(threshold[[1]]) > 1) {
       if (dim(data)[time_dim] != length(threshold[[1]])) {
-        stop("If parameter 'threshold' is a vector it must have the same length as data any time dimension.")
+        stop("If parameter 'threshold' is a vector it must have the same ", 
+             "length as data any time dimension.")
       } else {
         dim(threshold[[1]]) <- length(threshold[[1]])
         dim(threshold[[2]]) <- length(threshold[[2]])
@@ -265,8 +298,8 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
       if (any(names(dim(threshold[[1]])) %in% names(dim(data)))) {
         common_dims <- dim(threshold[[1]])[names(dim(threshold[[1]])) %in% names(dim(data))]
         if (!all(common_dims == dim(data)[names(common_dims)])) {
-          stop(paste0("Parameter 'data' and 'threshold' must have same length of ",
-                      "all common dimensions."))
+          stop("Parameter 'data' and 'threshold' must have same length of ",
+               "all common dimensions.")
         }
       }
     } else if (length(threshold[[1]]) == 1) {
@@ -276,7 +309,8 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
   } else {
     if (!is.array(threshold) && length(threshold) > 1) {
       if (dim(data)[time_dim] != length(threshold)) {
-        stop("If parameter 'threshold' is a vector it must have the same length as data time dimension.")
+        stop("If parameter 'threshold' is a vector it must have the same ",
+             "length as data time dimension.")
       } else {
         dim(threshold) <- length(threshold)
         names(dim(threshold)) <- time_dim
@@ -288,8 +322,8 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
       if (any(names(dim(threshold)) %in% names(dim(data)))) {
         common_dims <- dim(threshold)[names(dim(threshold)) %in% names(dim(data))]
         if (!all(common_dims == dim(data)[names(common_dims)])) {
-          stop(paste0("Parameter 'data' and 'threshold' must have same length of ",
-                      "all common dimensions."))
+          stop("Parameter 'data' and 'threshold' must have same length of ",
+               "all common dimensions.")
         }
       }
     } else if (length(threshold) == 1) {
@@ -313,27 +347,41 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
       if (length(op) == 1) {
         if (time_dim %in% names(dim(threshold))) {
           if (dim(threshold)[time_dim] == dim(data)[time_dim]) {
-            threshold <- SelectPeriodOnData(threshold, dates, start, end,
-                                            time_dim = time_dim, ncores = ncores)
+            threshold <- SelectPeriodOnData(data = threshold, dates = dates, 
+                                            start = start, end = end,
+                                            time_dim = time_dim, 
+                                            ncores = ncores)
           }
         }
       } else if (length(op) == 2) {
         if (time_dim %in% names(dim(threshold[[1]]))) {
           if (dim(threshold[[1]])[time_dim] == dim(data)[time_dim]) {
-            threshold[[1]] <- SelectPeriodOnData(threshold[[1]], dates, start, end,
-                                                 time_dim = time_dim, ncores = ncores)
-            threshold[[2]] <- SelectPeriodOnData(threshold[[2]], dates, start, end,
-                                                 time_dim = time_dim, ncores = ncores)
+            threshold[[1]] <- SelectPeriodOnData(data = threshold[[1]], 
+                                                 dates = dates, start = start, 
+                                                 end = end, time_dim = time_dim, 
+                                                 ncores = ncores)
+            threshold[[2]] <- SelectPeriodOnData(data = threshold[[2]], dates = dates, 
+                                                 start = start, end = end,
+                                                 time_dim = time_dim, 
+                                                 ncores = ncores)
           }
         }
       }
-      data <- SelectPeriodOnData(data, dates, start, end, 
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
+
     }
   }
   # diff
   if (length(op) == 2 & diff == TRUE) {
-    stop("Parameter 'diff' can't be TRUE if the parameter 'threshold' is a range of values.")
+    stop("Parameter 'diff' can't be TRUE if the parameter 'threshold' is a ",
+         "range of values.")
   } else  if (diff == TRUE) {
     if (length(threshold) != 1) {
       stop("Parameter 'diff' can't be TRUE if the parameter 'threshold' is not a scalar.")
@@ -345,8 +393,6 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
     threshold <- 0
   }
 
-  ###
-
   if (length(op) > 1) {
     thres1 <- threshold[[1]]
     thres2 <- threshold[[2]]
diff --git a/R/PeriodAccumulation.R b/R/PeriodAccumulation.R
index d181d8eeaa71100948ffa1594fce8d6ebadd376f..3b0d33d8ced0e1a730967cc4a5c36eef2be6ada0 100644
--- a/R/PeriodAccumulation.R
+++ b/R/PeriodAccumulation.R
@@ -30,8 +30,14 @@
 #'@param ncores An integer indicating the number of cores to use in parallel 
 #'  computation.
 #'
-#'@return A 's2dv_cube' object containing the indicator in the element
-#'\code{data}.
+#'@return An 's2dv_cube' object containing the indicator in the element
+#'\code{data} with dimensions of the input parameter 'data' except the dimension
+#'where the accumulation has been computed (specified with 'time_dim'). The 
+#''Dates' array is updated to the dates corresponding to the beginning of the 
+#'aggregated time period. A new element called 'time_bounds' will be added into 
+#'the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+#'containing two elements, the start and end dates of the aggregated period with 
+#'the same dimensions of 'Dates' element.
 #'
 #'@examples
 #'exp <- NULL
@@ -39,14 +45,16 @@
 #'                  ftime = 9, lat = 2, lon = 2))
 #'class(exp) <- 's2dv_cube'
 #'TP <- CST_PeriodAccumulation(exp)
-#'exp$data <- array(rnorm(5 * 3 * 214 * 2),
-#'                    c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
-#'exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
-#'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
-#'                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
-#'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
-#'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
-#'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'exp$data <- array(rnorm(5 * 3 * 214 * 2), 
+#'                  c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'exp$attrs$Dates <- Dates
 #'SprR <- CST_PeriodAccumulation(exp, start = list(21, 4), end = list(21, 6))
 #'dim(SprR$data)
 #'head(SprR$attrs$Dates)
@@ -55,6 +63,7 @@
 #'head(HarR$attrs$Dates)
 #'
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 CST_PeriodAccumulation <- function(data, start = NULL, end = NULL,
                                    time_dim = 'ftime', na.rm = FALSE,
@@ -73,14 +82,31 @@ CST_PeriodAccumulation <- function(data, start = NULL, end = NULL,
     }
   }
 
-  total <- PeriodAccumulation(data$data, dates = data$attrs$Dates, start, end,
+  Dates <- data$attrs$Dates
+  total <- PeriodAccumulation(data$data, dates = Dates, start, end,
                               time_dim = time_dim, na.rm = na.rm, ncores = ncores)
   data$data <- total
-  if (!is.null(start) && !is.null(end)) {
-     data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
-                                             start = start, end = end, 
-                                             time_dim = time_dim, 
-                                             ncores = ncores)
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end, 
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates  
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(Dates, time_dim, 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(Dates, time_dim, dim(Dates)[time_dim], drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
   }
   return(data)
 }
@@ -98,9 +124,9 @@ CST_PeriodAccumulation <- function(data, start = NULL, end = NULL,
 #'}
 #'
 #'@param data A multidimensional array with named dimensions.
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -128,15 +154,17 @@ CST_PeriodAccumulation <- function(data, start = NULL, end = NULL,
 #'TP <- PeriodAccumulation(exp, time_dim = 'ftime')
 #'data <- array(rnorm(5 * 3 * 214 * 2),
 #'              c(memb = 5, sdate = 3, time = 214, lon = 2)) 
-#'# ftime tested
 #'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
 #'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
 #'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
 #'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
 #'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
 #'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
-#'SprR <- PeriodAccumulation(data, dates = Dates, start = list(21, 4), end = list(21, 6))
-#'HarR <- PeriodAccumulation(data, dates = Dates, start = list(21, 8), end = list(21, 10))
+#'dim(Dates) <- c(sdate = 3, time = 214)
+#'SprR <- PeriodAccumulation(data, dates = Dates, start = list(21, 4), 
+#'                           end = list(21, 6))
+#'HarR <- PeriodAccumulation(data, dates = Dates, start = list(21, 8), 
+#'                           end = list(21, 10))
 #'
 #'@import multiApply
 #'@export
@@ -153,14 +181,23 @@ PeriodAccumulation <- function(data, dates = NULL, start = NULL, end = NULL,
     dim(data) <- length(data)
     names(dim(data)) <- time_dim
   }
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
       }
-      data <- SelectPeriodOnData(data, dates, start, end, 
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data, dates, start, end, 
+                                   time_dim = time_dim, ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }
   total <- Apply(list(data), target_dims = time_dim, fun = sum,
diff --git a/R/PeriodMean.R b/R/PeriodMean.R
index 85a12a74272726c3694425fe1728181409d956b1..db6a78f273965f19ac056e3b4f5d577e7b239963 100644
--- a/R/PeriodMean.R
+++ b/R/PeriodMean.R
@@ -31,19 +31,33 @@
 #'  computation.
 #'
 #'@return An 's2dv_cube' object containing the indicator in the element 
-#'  \code{data}.
+#'\code{data} with dimensions of the input parameter 'data' except the 
+#'dimension where the mean has been computed (specified with 'time_dim'). The 
+#''Dates' array is updated to the dates corresponding to the beginning of the 
+#'aggregated time period. A new element called 'time_bounds' will be added into 
+#'the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+#'containing two elements, the start and end dates of the aggregated period with 
+#'the same dimensions of 'Dates' element.
 #'
 #'@examples
 #'exp <- NULL
-#'exp$data <- array(rnorm(45), dim = c(member = 7, ftime = 8))
+#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, ftime = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, ftime = 3)
+#'exp$attrs$Dates <- Dates
 #'class(exp) <- 's2dv_cube'
-#'exp$attrs$Dates <- c(seq(as.Date("01-07-1993", "%d-%m-%Y", tz = 'UTC'),
-#'                         as.Date("01-08-1993","%d-%m-%Y", tz = 'UTC'), "day"),
-#'                     seq(as.Date("01-07-1994", "%d-%m-%Y", tz = 'UTC'),
-#'                         as.Date("01-08-1994","%d-%m-%Y", tz = 'UTC'), "day"))
-#'SA <- CST_PeriodMean(exp)
+#'
+#'SA <- CST_PeriodMean(exp, start = list(01, 12), end = list(01, 01))
 #'
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 CST_PeriodMean <- function(data, start = NULL, end = NULL,
                            time_dim = 'ftime', na.rm = FALSE,
@@ -62,14 +76,35 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
     }
   }
 
-  total <- PeriodMean(data = data$data, dates = data$attrs$Dates, start, end,
+  Dates <- data$attrs$Dates
+  total <- PeriodMean(data = data$data, dates = Dates, start = start, end = end,
                       time_dim = time_dim, na.rm = na.rm, ncores = ncores)
+  
   data$data <- total
-  if (!is.null(start) && !is.null(end)) {
-     data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
-                                             start = start, end = end,
-                                             time_dim = time_dim, 
-                                             ncores = ncores)
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
   }
   return(data)
 }
@@ -87,9 +122,9 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
 #'}
 #'
 #'@param data A multidimensional array with named dimensions.
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -100,7 +135,7 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
 #'  the period and the final month of the period. By default it is set to NULL 
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the dimension to 
-#'  compute the indicator. By default, it is set to 'ftime'. More than one 
+#'  compute the indicator. By default, it is set to 'time'. More than one 
 #'  dimension name matching the dimensions provided in the object 
 #'  \code{data$data} can be specified.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
@@ -112,8 +147,17 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
 #'indicator in the element \code{data}.
 #'
 #'@examples
-#'exp <- array(rnorm(56), dim = c(member = 7, ftime = 8))
-#'SA <- PeriodMean(exp, time_dim = 'ftime')
+#'data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'SA <- PeriodMean(data, dates = Dates, start = list(01, 12), end = list(01, 01))
 #'
 #'@import multiApply
 #'@export
@@ -130,21 +174,29 @@ PeriodMean <- function(data, dates = NULL, start = NULL, end = NULL,
     dim(data) <- length(data)
     names(data) <- time_dim
   }
-  if (is.null(dates)) {
-    warning("Parameter 'dates' is NULL and the Average of the ",
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
               "full data provided in 'data' is computed.")
-  } else {
-    if (!is.null(start) && !is.null(end)) {
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
-       stop("Parameter 'start' and 'end' must be lists indicating the ",
-             "day and the month of the period start and end.")
+      stop("Parameter 'start' and 'end' must be lists indicating the ",
+            "day and the month of the period start and end.")
+      }
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
       }
-      data <- SelectPeriodOnData(data, dates, start, end,
-                                 time_dim = time_dim, ncores = ncores)
     }
   }
   total <- Apply(list(data), target_dims = time_dim, fun = mean,
                  na.rm = na.rm, ncores = ncores)$output1
+  return(total)
 }
 
 
diff --git a/R/QThreshold.R b/R/QThreshold.R
index 49217dd20e8986ecb786d8f4232212743132a989..c5089df13556930fc7913da922a5a6c2fffa3466 100644
--- a/R/QThreshold.R
+++ b/R/QThreshold.R
@@ -57,15 +57,19 @@
 #'exp$data <- array(abs(rnorm(112)*26), dim = c(member = 7, sdate = 8, ftime = 2))
 #'class(exp) <- 's2dv_cube'
 #'exp_probs <- CST_QThreshold(exp, threshold)
-#'exp$data <- array(rnorm(5 * 3 * 214 * 2),
-#'                    c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+#'
+#'exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+#'                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 #'exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
-#'exp_probs <- CST_QThreshold(exp, threshold)
+#'dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
+#'class(exp) <- 's2dv_cube'
+#'exp_probs <- CST_QThreshold(exp, threshold, start = list(21, 4), 
+#'                            end = list(21, 6))
 #' 
 #'@import multiApply
 #'@export
@@ -92,6 +96,7 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
                       start, end, time_dim = time_dim, memb_dim = memb_dim,
                       sdate_dim = sdate_dim, ncores = ncores)
   data$data <- probs
+  data$dims <- dim(probs)
   if (!is.null(start) && !is.null(end)) {
     data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
                                             start = start, end = end, 
@@ -127,9 +132,9 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
 #'@param threshold A multidimensional array with named dimensions in the same 
 #'  units as parameter 'data' and with the common dimensions of the element 
 #'  'data' of the same length.
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -140,7 +145,7 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
 #'  the period and the final month of the period. By default it is set to NULL 
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the temporal 
-#'  dimension. By default, it is set to 'ftime'. More than one dimension name 
+#'  dimension. By default, it is set to 'time'. More than one dimension name 
 #'  matching the dimensions provided in the object \code{data$data} can be 
 #'  specified. This dimension is required to subset the data in a requested 
 #'  period.
@@ -157,14 +162,24 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
 #'@examples
 #'threshold = 25
 #'data <- array(rnorm(5 * 3 * 20 * 2, mean = 26), 
-#'              c(member = 5, sdate = 3, time = 20, lon = 2)) 
-#'thres_q <- QThreshold(data, threshold)
+#'              c(member = 5, sdate = 3, time = 214, lon = 2)) 
+#'
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, time = 214)
+#'
+#'thres_q <- QThreshold(data, threshold, dates = Dates, time_dim = 'time',
+#'                      start = list(21, 4), end = list(21, 6))
 #' 
 #'@import multiApply
 #'@export
 QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
-                       time_dim = 'time', memb_dim = 'member', sdate_dim = 'sdate',
-                       ncores = NULL) {
+                       time_dim = 'time', memb_dim = 'member', 
+                       sdate_dim = 'sdate', ncores = NULL) {
   # Initial checks
   ## data
   if (is.null(data)) {
@@ -202,8 +217,12 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
   if (is.null(memb_dim)) {
     memb_dim <- 99999
   }
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
@@ -220,10 +239,16 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
           }
         }
       }
-      data <- SelectPeriodOnData(data, dates, start, end, 
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }
+  
   if (length(threshold) == 1) {
     if (memb_dim %in% names(dim(data))) {
       probs <- Apply(list(data), target_dims = c(memb_dim, sdate_dim),
diff --git a/R/Threshold.R b/R/Threshold.R
index 3122c1289394f9e6039c2e462b1776f448a49519..73e3715a9c2b9ad558869406128c50368ca02726 100644
--- a/R/Threshold.R
+++ b/R/Threshold.R
@@ -8,7 +8,7 @@
 #'@param data An 's2dv_cube' object as provided function \code{CST_Load} in
 #'  package CSTools.
 #'@param threshold A single scalar or vector indicating the relative 
-#'  threshold(s).
+#'  threshold(s). It must contain values between 0 and 1.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to selectfrom the data by providing a list of two elements: the initial date
 #'  of the period and the initial month of the period. By default it is set to 
@@ -40,13 +40,14 @@
 #'threshold <- 0.9
 #'exp <- NULL
 #'exp$data <- array(rnorm(5 * 3 * 214 * 2),
-#'                    c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+#'                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 #'exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 #'class(exp) <- 's2dv_cube'
 #'exp_probs <- CST_Threshold(exp, threshold, start = list(21, 4), end = list(21, 6))
 #' 
@@ -74,6 +75,7 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL,
                      time_dim = time_dim, memb_dim = memb_dim,
                      sdate_dim = sdate_dim, na.rm = na.rm, ncores = ncores)
   data$data <- thres
+  data$dims <- dim(thres)
   if (!is.null(start) && !is.null(end)) {
     data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
                                             start = start, end = end, 
@@ -91,10 +93,10 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL,
 #'
 #'@param data A multidimensional array with named dimensions.
 #'@param threshold A single scalar or vector indicating the relative 
-#'  threshold(s).
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'  threshold(s). It must contain values between 0 and 1.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -105,7 +107,7 @@ CST_Threshold <- function(data, threshold, start = NULL, end = NULL,
 #'  the period and the final month of the period. By default it is set to NULL 
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the temporal 
-#'  dimension. By default, it is set to 'ftime'. More than one dimension name 
+#'  dimension. By default, it is set to 'time'. More than one dimension name 
 #'  matching the dimensions provided in the object \code{data$data} can be 
 #'  specified. This dimension is required to subset the data in a requested 
 #'  period.
@@ -147,7 +149,7 @@ Threshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
     names(dim(data)) <- c(memb_dim, sdate_dim)
   }
   if (is.null(threshold)) {
-      stop("Parameter 'threshold' cannot be NULL.")
+    stop("Parameter 'threshold' cannot be NULL.")
   }
    if (!is.numeric(threshold)) {
     stop("Parameter 'threshold' must be numeric.")
@@ -155,14 +157,24 @@ Threshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
   if (is.null(names(dim(data)))) {
     stop("Parameter 'data' must have named dimensions.")
   }
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
       }
-      data <- SelectPeriodOnData(data, dates, start, end, 
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }   
   if (!is.null(memb_dim)) {
diff --git a/R/TotalSpellTimeExceedingThreshold.R b/R/TotalSpellTimeExceedingThreshold.R
index 3ee22a27b21114a970476974067366adca2ffc7c..38a3ff4d897848113afb644de9adf764222295d3 100644
--- a/R/TotalSpellTimeExceedingThreshold.R
+++ b/R/TotalSpellTimeExceedingThreshold.R
@@ -48,7 +48,13 @@
 #'  computation.
 #'
 #'@return An 's2dv_cube' object containing the number of days that are part of a
-#'spell within a threshold in element \code{data}.
+#'spell within a threshold in element \code{data} with dimensions of the input 
+#'parameter 'data' except the dimension where the indicator has been computed.
+#'The 'Dates' array is updated to the dates corresponding to the beginning of 
+#'the aggregated time period. A new element called 'time_bounds' will be added  
+#'into the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+#'containing two elements, the start and end dates of the aggregated period with 
+#'the same dimensions of 'Dates' element.
 #'
 #'@examples
 #'exp <- NULL
@@ -60,10 +66,14 @@
 #'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 #'class(exp) <- 's2dv_cube'
-#'TTSET <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 23, spell = 3)
+#'TTSET <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 23, spell = 3,
+#'                                              start = list(21, 4), 
+#'                                              end = list(21, 6))
 #' 
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>',
                                                  start = NULL, end = NULL,
@@ -95,19 +105,41 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
       threshold[[2]] <- threshold[[2]]$data
     }
   }
+
+  Dates <- data$attrs$Dates
   
-  total <- TotalSpellTimeExceedingThreshold(data$data, data$attrs$Dates,
+  total <- TotalSpellTimeExceedingThreshold(data$data, Dates,
                                             threshold = threshold, spell = spell, 
                                             op = op, start = start, end = end, 
                                             time_dim = time_dim, 
                                             ncores = ncores)
   data$data <- total
-  if (!is.null(start) && !is.null(end)) {
-    data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
-                                            start = start, end = end,
-                                            time_dim = time_dim, 
-                                            ncores = ncores)
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
   }
+
   return(data)
 }
 #'Total Spell Time Exceeding Threshold
@@ -143,9 +175,9 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
 #'  are used it has to be a vector of a pair of two logical operators: 
 #'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
 #'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -171,9 +203,19 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
 #'values by values exceeding the threshold.
 
 #'@examples
-#'data <- array(rnorm(120), c(member = 1, sdate = 2, ftime = 20, lat = 4))
-#'threshold <- array(rnorm(4), c(lat = 4))
-#'total <- TotalSpellTimeExceedingThreshold(data, threshold, spell = 6)
+#'data <- array(1:100, c(member = 5, sdate = 3, ftime = 214, lon = 2))
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'
+#'threshold <- array(1:4, c(lat = 4))
+#'total <- TotalSpellTimeExceedingThreshold(data, threshold, dates = Dates, 
+#'                                          spell = 6, start = list(21, 4), 
+#'                                          end = list(21, 6))
 #' 
 #'@import multiApply
 #'@export
@@ -311,8 +353,11 @@ TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>',
     }
   }
   # dates
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
@@ -334,6 +379,14 @@ TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>',
           }
         }
       }
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
       data <- SelectPeriodOnData(data, dates, start, end, 
                                  time_dim = time_dim, ncores = ncores)
     }
diff --git a/R/TotalTimeExceedingThreshold.R b/R/TotalTimeExceedingThreshold.R
index ceda1eee6ee47ef0e08ad0a1649d93e986286df2..94d2c53096c3d15fc1dac378550a175dd95d47dd 100644
--- a/R/TotalTimeExceedingThreshold.R
+++ b/R/TotalTimeExceedingThreshold.R
@@ -54,11 +54,17 @@
 #'
 #'@return An 's2dv_cube' object containing in element \code{data} the total 
 #'number of the corresponding units of the data frequency that a variable is 
-#'exceeding a threshold during a period.
+#'exceeding a threshold during a period with dimensions of the input parameter 
+#''data' except the dimension where the indicator has been computed. The 
+#''Dates' array is updated to the dates corresponding to the beginning of the 
+#'aggregated time period. A new element called 'time_bounds' will be added into 
+#'the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+#'containing two elements, the start and end dates of the aggregated period with 
+#'the same dimensions of 'Dates' element.
 #' 
 #'@examples
 #'exp <- NULL
-#'exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*280),
+#'exp$data <- array(rnorm(5 * 3 * 214 * 2)*23,
 #'                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 #'exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
@@ -66,10 +72,13 @@
 #'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
 #'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
 #'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 #'class(exp) <- 's2dv_cube'
-#'DOT <- CST_TotalTimeExceedingThreshold(exp, threshold = 280)
+#'DOT <- CST_TotalTimeExceedingThreshold(exp, threshold = 23, start = list(21, 4), 
+#'                                       end = list(21, 6))
 #'
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
                                             start = NULL, end = NULL,
@@ -101,17 +110,38 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
       threshold[[2]] <- threshold[[2]]$data
     }
   }
-  total <- TotalTimeExceedingThreshold(data$data, dates = data$attrs$Dates,
+  
+  Dates <- data$attrs$Dates
+  total <- TotalTimeExceedingThreshold(data = data$data, dates = Dates,
                                        threshold = threshold, op = op,
                                        start = start, end = end, 
                                        time_dim = time_dim, na.rm = na.rm, 
                                        ncores = ncores)
   data$data <- total
-  if (!is.null(start) && !is.null(end)) {
-     data$attrs$Dates <- SelectPeriodOnDates(dates = data$attrs$Dates,
-                                             start = start, end = end,
-                                             time_dim = time_dim, 
-                                             ncores = ncores)
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
   }
   return(data)
 }
@@ -151,9 +181,9 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
 #'  are used it has to be a vector of a pair of two logical operators: 
 #'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
 #'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is 
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to define the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
@@ -173,12 +203,21 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
 #'
 #'@return A multidimensional array with named dimensions containing the total 
 #'number of the corresponding units of the data frequency that a variable is 
-#'exceeding a threshold during a period.
+#'exceeding a threshold during a period with dimensions of the input parameter 
+#''data' except the dimension where the indicator has been computed.
 #'
 #'@examples
-#'exp <- array(abs(rnorm(5 * 3 * 214 * 2)*280),
-#'             c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
-#'DOT <- TotalTimeExceedingThreshold(exp, threshold = 300, time_dim = 'ftime')
+#'data <- array(rnorm(5 * 3 * 214 * 2)*23,
+#'              c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'DOT <- TotalTimeExceedingThreshold(data, threshold = 23, dates = Dates,
+#'                                   start = list(21, 4), end = list(21, 6))
 #'
 #'@import multiApply
 #'@export
@@ -313,8 +352,11 @@ TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
     }
   }
   # dates
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
@@ -336,8 +378,14 @@ TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
           }
         }
       }
-      data <- SelectPeriodOnData(data, dates, start, end, 
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }
 
diff --git a/R/WindCapacityFactor.R b/R/WindCapacityFactor.R
index 8ed20844370f7036d2cee074df64eab79308ce18..76092dd7498b64f5c9da9db3cd8f0c352fa8e6e1 100644
--- a/R/WindCapacityFactor.R
+++ b/R/WindCapacityFactor.R
@@ -42,14 +42,23 @@
 #'@examples
 #'wind <- NULL
 #'wind$data <- array(rweibull(n = 100, shape = 2, scale = 6), 
-#'                   c(member = 10, lat = 2, lon = 5))
+#'                   c(member = 5, sdate = 3, ftime = 214, lon = 2, lat = 5))
 #'wind$coords <- list(lat =  c(40, 41), lon = 1:5)
 #'variable <- list(varName = 'sfcWind', 
 #'                 metadata = list(sfcWind = list(level = 'Surface')))
 #'wind$attrs <- list(Variable = variable, Datasets = 'synthetic', 
 #'                   when = Sys.time(), Dates = '1990-01-01 00:00:00')
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'wind$attrs$Dates <- Dates
 #'class(wind) <- 's2dv_cube'
-#'WCF <- CST_WindCapacityFactor(wind, IEC_class = "III")
+#'WCF <- CST_WindCapacityFactor(wind, IEC_class = "III", 
+#'                              start = list(21, 4), end = list(21, 6))
 #'
 #'@export
 CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III", "III"),
@@ -71,8 +80,11 @@ CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/II
   
   WindCapacity <- WindCapacityFactor(wind = wind$data, IEC_class = IEC_class, 
                                      dates = wind$attrs$Dates, start = start, 
-                                     end = end, ncores = ncores)
+                                     end = end, time_dim = time_dim, 
+                                     ncores = ncores)
   wind$data <- WindCapacity
+  wind$dims <- dim(WindCapacity)
+
   if ('Variable' %in% names(wind$attrs)) {
     if ('varName' %in% names(wind$attrs$Variable)) {
       wind$attrs$Variable$varName <- 'WindCapacityFactor'
@@ -111,9 +123,9 @@ CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/II
 #'  respectively. Classes \code{'I/II'} and \code{'II/III'} indicate
 #'  intermediate turbines that fit both classes. More details of the five
 #'  turbines and a plot of its power curves can be found in Lledó et al. (2019).
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period
 #'  to select from the data by providing a list of two elements: the initial
 #'  date of the period and the initial month of the period. By default it is set
@@ -124,7 +136,7 @@ CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/II
 #'  the period and the final month of the period. By default it is set to NULL
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the dimension to
-#'  compute the indicator. By default, it is set to 'ftime'. More than one
+#'  compute the indicator. By default, it is set to 'time'. More than one
 #'  dimension name matching the dimensions provided in the object
 #'  \code{data$data} can be specified.
 #'@param ncores An integer indicating the number of cores to use in parallel
@@ -134,8 +146,19 @@ CST_WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/II
 #'  Capacity Factor (unitless).
 #'
 #'@examples
-#'wind <- rweibull(n = 100, shape = 2, scale = 6)
-#'WCF <- WindCapacityFactor(wind, IEC_class = "III")
+#'wind <- array(rweibull(n = 32100, shape = 2, scale = 6), 
+#'              c(member = 5, sdate = 3, time = 214, lon = 2, lat = 5))
+#'
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, time = 214)
+#'
+#'WCF <- WindCapacityFactor(wind, IEC_class = "III", dates = Dates, 
+#'                          start = list(21, 4), end = list(21, 6))
 #'
 #'@importFrom stats approxfun
 #'@importFrom utils read.delim
@@ -153,14 +176,24 @@ WindCapacityFactor <- function(wind, IEC_class = c("I", "I/II", "II", "II/III",
 	)
 	pc_file <- system.file("power_curves", pc_files[IEC_class], package = "CSIndicators", mustWork = T)
 	pc <- read_pc(pc_file)
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
         stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
       }
-      wind <- SelectPeriodOnData(wind, dates, start, end,
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        wind <- SelectPeriodOnData(data = wind, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'wind' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }
 
diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R
index 357820999cb83f3c53bfdb08795d8265094026a9..3944e3d99c45e6ed289e7f6820fa20e684d6a841 100644
--- a/R/WindPowerDensity.R
+++ b/R/WindPowerDensity.R
@@ -6,8 +6,8 @@
 #'@description It is computed as 0.5*ro*wspd^3. As this function is non-linear, 
 #'it will give inaccurate results if used with period means.
 #'
-#'@param wind An s2dv_cube object with instantaneous wind speeds expressed in m/s 
-#'  obtained from CST_Load or s2dv_cube functions from CSTools pacakge.
+#'@param wind An 's2dv_cube' object with instantaneous wind speeds expressed in 
+#'  m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge.
 #'@param ro A scalar, or alternatively a multidimensional array with the same
 #'  dimensions as wind, with the air density expressed in kg/m^3. By default it
 #'  takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.
@@ -32,14 +32,23 @@
 #'@examples
 #'wind <- NULL
 #'wind$data <- array(rweibull(n = 100, shape = 2, scale = 6), 
-#'                   c(member = 10, lat = 2, lon = 5))
+#'                   c(member = 5, sdate = 3, ftime = 214, lon = 2, lat = 5))
 #'wind$coords <- list(lat =  c(40, 41), lon = 1:5)
 #'variable <- list(varName = 'sfcWind', 
 #'                 metadata = list(sfcWind = list(level = 'Surface')))
 #'wind$attrs <- list(Variable = variable, Datasets = 'synthetic', 
 #'                   when = Sys.time(), Dates = '1990-01-01 00:00:00')
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, ftime = 214)
+#'wind$attrs$Dates <- Dates
 #'class(wind) <- 's2dv_cube'
-#'WCF <- CST_WindPowerDensity(wind)
+#'WPD <- CST_WindPowerDensity(wind, start = list(21, 4), 
+#'                            end = list(21, 6))
 #'
 #'@export
 CST_WindPowerDensity <- function(wind, ro = 1.225, start = NULL, end = NULL, 
@@ -59,8 +68,10 @@ CST_WindPowerDensity <- function(wind, ro = 1.225, start = NULL, end = NULL,
   }
   WindPower <- WindPowerDensity(wind = wind$data, ro = ro, 
                                 dates = wind$attrs$Dates, start = start, 
-                                end = end, ncores = ncores)
+                                end = end, time_dim = time_dim, 
+                                ncores = ncores)
   wind$data <- WindPower
+  wind$dims <- dim(WindPower)
   if ('Variable' %in% names(wind$attrs)) {
     if ('varName' %in% names(wind$attrs$Variable)) {
       wind$attrs$Variable$varName <- 'WindPowerDensity'
@@ -88,9 +99,9 @@ CST_WindPowerDensity <- function(wind, ro = 1.225, start = NULL, end = NULL,
 #'@param ro A scalar, or alternatively a multidimensional array with the same
 #'  dimensions as wind, with the air density expressed in kg/m^3. By default it
 #'  takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.
-#'@param dates A vector of dates or a multidimensional array of dates with named
-#'  dimensions matching the dimensions on parameter 'data'. By default it is
-#'  NULL, to select a period this parameter must be provided.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
 #'@param start An optional parameter to defined the initial date of the period
 #'  to select from the data by providing a list of two elements: the initial
 #'  date of the period and the initial month of the period. By default it is set
@@ -101,7 +112,7 @@ CST_WindPowerDensity <- function(wind, ro = 1.225, start = NULL, end = NULL,
 #'  the period and the final month of the period. By default it is set to NULL
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the dimension to
-#'  compute the indicator. By default, it is set to 'ftime'. More than one
+#'  compute the indicator. By default, it is set to 'time'. More than one
 #'  dimension name matching the dimensions provided in the object
 #'  \code{data$data} can be specified.
 #'@param ncores An integer indicating the number of cores to use in parallel
@@ -111,20 +122,39 @@ CST_WindPowerDensity <- function(wind, ro = 1.225, start = NULL, end = NULL,
 #'Density expressed in W/m^2.
 #'
 #'@examples
-#'wind <- rweibull(n = 100, shape = 2, scale = 6)
-#'WPD <- WindPowerDensity(wind)
+#'wind <- array(rweibull(n = 32100, shape = 2, scale = 6), 
+#'              c(member = 5, sdate = 3, time = 214, lon = 2, lat = 5))
+#'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+#'           seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+#'               as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+#'dim(Dates) <- c(sdate = 3, time = 214)
+#'WPD <- WindPowerDensity(wind, dates = Dates, start = list(21, 4), 
+#'                        end = list(21, 6))
 #'
 #'@export
 WindPowerDensity <- function(wind, ro = 1.225, dates = NULL, start = NULL, 
                              end = NULL, time_dim = 'time', ncores = NULL) {
-  if (!is.null(dates)) {
-    if (!is.null(start) && !is.null(end)) {
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
       if (!any(c(is.list(start), is.list(end)))) {
        stop("Parameter 'start' and 'end' must be lists indicating the ",
              "day and the month of the period start and end.")
       }
-      wind <- SelectPeriodOnData(wind, dates, start, end,
-                                 time_dim = time_dim, ncores = ncores)
+      if (!is.null(dim(dates))) {
+        wind <- SelectPeriodOnData(data = wind, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'wind' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
     }
   }
 	return(0.5 * ro * wind^3)
diff --git a/man/AccumulationExceedingThreshold.Rd b/man/AccumulationExceedingThreshold.Rd
index 172592c89c091e771f5a13f2084f594dc25b51e2..0f9a64ab6d658b844d7067ab1dcab9c6c9a4d662 100644
--- a/man/AccumulationExceedingThreshold.Rd
+++ b/man/AccumulationExceedingThreshold.Rd
@@ -39,9 +39,9 @@ c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 between data and threshold (TRUE) or not (FALSE by default). It can only be
 TRUE if a unique threshold is used.}
 
-\item{dates}{A vector of dates or a multidimensional array with of dates with 
-named dimensions matching the dimensions on parameter 'data'. By default it
-is NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to define the initial date of the period
 to select from the data by providing a list of two elements: the initial 
@@ -86,12 +86,6 @@ function:
 # Assuming data is already (tasmax + tasmin)/2 - 10
 data <- array(rnorm(5 * 3 * 214 * 2, mean = 25, sd = 3),
                    c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
-Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
-              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
-          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
-              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
-          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
-              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
 GDD <- AccumulationExceedingThreshold(data, threshold = 0, start = list(1, 4),
                                      end = list(31, 10))
 }
diff --git a/man/CST_AccumulationExceedingThreshold.Rd b/man/CST_AccumulationExceedingThreshold.Rd
index bc0eb83f4d4d1b168b956c63facbbb6ea5de4fc4..ff02e4c7538e3fcd5986bbc973c9eab46c2e92cd 100644
--- a/man/CST_AccumulationExceedingThreshold.Rd
+++ b/man/CST_AccumulationExceedingThreshold.Rd
@@ -63,7 +63,12 @@ computation.}
 \value{
 An 's2dv_cube' object containing the aggregated values in the element
 \code{data} with dimensions of the input parameter 'data' except the dimension
-where the indicator has been computed.
+where the indicator has been computed. The 'Dates' array is updated to 
+the dates corresponding to the beginning of the aggregated time period. A new 
+element called 'time_bounds' will be added into the 'attrs' element in the 
+'s2dv_cube' object. It consists of a list containing two elements, the start 
+and end dates of the aggregated period with the same dimensions of 'Dates' 
+element.
 }
 \description{
 The accumulation (sum) of a variable in the days (or time steps) that the 
@@ -80,9 +85,19 @@ function:
 }
 \examples{
 exp <- NULL
-exp$data <- array(rnorm(216)*200, dim = c(dataset = 1, member = 2, sdate = 3, 
-                 ftime = 9, lat = 2, lon = 2))
+exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*100), 
+                 c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
 class(exp) <- 's2dv_cube'
-DOT <- CST_AccumulationExceedingThreshold(exp, threshold = 280)
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+exp$attrs$Dates <- Dates
+AT <- CST_AccumulationExceedingThreshold(data = exp, threshold = 100,
+                                        start = list(21, 4), 
+                                        end = list(21, 6))
 
 }
diff --git a/man/CST_PeriodAccumulation.Rd b/man/CST_PeriodAccumulation.Rd
index 39287052dc1e55d5bef9e24046e875bee38cdb73..77f4a389e28c8b6baeba40b2a063df57b2037a23 100644
--- a/man/CST_PeriodAccumulation.Rd
+++ b/man/CST_PeriodAccumulation.Rd
@@ -40,8 +40,14 @@ not (FALSE).}
 computation.}
 }
 \value{
-A 's2dv_cube' object containing the indicator in the element
-\code{data}.
+An 's2dv_cube' object containing the indicator in the element
+\code{data} with dimensions of the input parameter 'data' except the dimension
+where the accumulation has been computed (specified with 'time_dim'). The 
+'Dates' array is updated to the dates corresponding to the beginning of the 
+aggregated time period. A new element called 'time_bounds' will be added into 
+the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+containing two elements, the start and end dates of the aggregated period with 
+the same dimensions of 'Dates' element.
 }
 \description{
 Period Accumulation computes the sum (accumulation) of a given variable in a 
@@ -60,14 +66,16 @@ exp$data <- array(rnorm(216)*200, dim = c(dataset = 1, member = 2, sdate = 3,
                  ftime = 9, lat = 2, lon = 2))
 class(exp) <- 's2dv_cube'
 TP <- CST_PeriodAccumulation(exp)
-exp$data <- array(rnorm(5 * 3 * 214 * 2),
-                   c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
-exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
-                        as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
-                    seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
-                        as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
-                    seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
-                        as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+exp$data <- array(rnorm(5 * 3 * 214 * 2), 
+                 c(memb = 5, sdate = 3, ftime = 214, lon = 2)) 
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+exp$attrs$Dates <- Dates
 SprR <- CST_PeriodAccumulation(exp, start = list(21, 4), end = list(21, 6))
 dim(SprR$data)
 head(SprR$attrs$Dates)
diff --git a/man/CST_PeriodMean.Rd b/man/CST_PeriodMean.Rd
index b1004ad56896fe2aa943677debd1aeec1590888f..0aa4aa359645d70ba850f556586488e8980a431d 100644
--- a/man/CST_PeriodMean.Rd
+++ b/man/CST_PeriodMean.Rd
@@ -41,7 +41,13 @@ computation.}
 }
 \value{
 An 's2dv_cube' object containing the indicator in the element 
- \code{data}.
+\code{data} with dimensions of the input parameter 'data' except the 
+dimension where the mean has been computed (specified with 'time_dim'). The 
+'Dates' array is updated to the dates corresponding to the beginning of the 
+aggregated time period. A new element called 'time_bounds' will be added into 
+the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+containing two elements, the start and end dates of the aggregated period with 
+the same dimensions of 'Dates' element.
 }
 \description{
 Period Mean computes the average (mean) of a given variable in a period.
@@ -56,12 +62,19 @@ this function:
 }
 \examples{
 exp <- NULL
-exp$data <- array(rnorm(45), dim = c(member = 7, ftime = 8))
+exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, ftime = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, ftime = 3)
+exp$attrs$Dates <- Dates
 class(exp) <- 's2dv_cube'
-exp$attrs$Dates <- c(seq(as.Date("01-07-1993", "\%d-\%m-\%Y", tz = 'UTC'),
-                        as.Date("01-08-1993","\%d-\%m-\%Y", tz = 'UTC'), "day"),
-                    seq(as.Date("01-07-1994", "\%d-\%m-\%Y", tz = 'UTC'),
-                        as.Date("01-08-1994","\%d-\%m-\%Y", tz = 'UTC'), "day"))
-SA <- CST_PeriodMean(exp)
+
+SA <- CST_PeriodMean(exp, start = list(01, 12), end = list(01, 01))
 
 }
diff --git a/man/CST_QThreshold.Rd b/man/CST_QThreshold.Rd
index eda0fd1ced67dd2d87d4b6924f343e4e27304ea1..5a68bc3823bad195196c00fccf5c3024202078ee 100644
--- a/man/CST_QThreshold.Rd
+++ b/man/CST_QThreshold.Rd
@@ -83,14 +83,18 @@ exp <- NULL
 exp$data <- array(abs(rnorm(112)*26), dim = c(member = 7, sdate = 8, ftime = 2))
 class(exp) <- 's2dv_cube'
 exp_probs <- CST_QThreshold(exp, threshold)
-exp$data <- array(rnorm(5 * 3 * 214 * 2),
-                   c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+
+exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+                 c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
-exp_probs <- CST_QThreshold(exp, threshold)
+dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
+class(exp) <- 's2dv_cube'
+exp_probs <- CST_QThreshold(exp, threshold, start = list(21, 4), 
+                           end = list(21, 6))
 
 }
diff --git a/man/CST_Threshold.Rd b/man/CST_Threshold.Rd
index ffe06000d67a6c1371749a872a3d48d1d0dc4be8..e513ec058b650f776fffb569cdcd7ceb6c0547d6 100644
--- a/man/CST_Threshold.Rd
+++ b/man/CST_Threshold.Rd
@@ -21,7 +21,7 @@ CST_Threshold(
 package CSTools.}
 
 \item{threshold}{A single scalar or vector indicating the relative 
-threshold(s).}
+threshold(s). It must contain values between 0 and 1.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to selectfrom the data by providing a list of two elements: the initial date
@@ -67,13 +67,14 @@ given a dataset.
 threshold <- 0.9
 exp <- NULL
 exp$data <- array(rnorm(5 * 3 * 214 * 2),
-                   c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+                 c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 class(exp) <- 's2dv_cube'
 exp_probs <- CST_Threshold(exp, threshold, start = list(21, 4), end = list(21, 6))
 
diff --git a/man/CST_TotalSpellTimeExceedingThreshold.Rd b/man/CST_TotalSpellTimeExceedingThreshold.Rd
index e2f7d26c15bd92c2ef82bcf21d250ce8aac514df..0715414fcc4f22dee1f4ad506675c9b395170759 100644
--- a/man/CST_TotalSpellTimeExceedingThreshold.Rd
+++ b/man/CST_TotalSpellTimeExceedingThreshold.Rd
@@ -56,7 +56,13 @@ computation.}
 }
 \value{
 An 's2dv_cube' object containing the number of days that are part of a
-spell within a threshold in element \code{data}.
+spell within a threshold in element \code{data} with dimensions of the input 
+parameter 'data' except the dimension where the indicator has been computed.
+The 'Dates' array is updated to the dates corresponding to the beginning of 
+the aggregated time period. A new element called 'time_bounds' will be added  
+into the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+containing two elements, the start and end dates of the aggregated period with 
+the same dimensions of 'Dates' element.
 }
 \description{
 The number of days (when daily data is provided) that are part of a spell 
@@ -84,8 +90,11 @@ exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"),
                         as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 class(exp) <- 's2dv_cube'
-TTSET <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 23, spell = 3)
+TTSET <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 23, spell = 3,
+                                             start = list(21, 4), 
+                                             end = list(21, 6))
 
 }
 \seealso{
diff --git a/man/CST_TotalTimeExceedingThreshold.Rd b/man/CST_TotalTimeExceedingThreshold.Rd
index b09ae53bb10e5b46fdd792a7d523bd70300c4462..4a1f7366bfe3b7ba4fb266fe5259f556ae0d60df 100644
--- a/man/CST_TotalTimeExceedingThreshold.Rd
+++ b/man/CST_TotalTimeExceedingThreshold.Rd
@@ -58,7 +58,13 @@ computation.}
 \value{
 An 's2dv_cube' object containing in element \code{data} the total 
 number of the corresponding units of the data frequency that a variable is 
-exceeding a threshold during a period.
+exceeding a threshold during a period with dimensions of the input parameter 
+'data' except the dimension where the indicator has been computed. The 
+'Dates' array is updated to the dates corresponding to the beginning of the 
+aggregated time period. A new element called 'time_bounds' will be added into 
+the 'attrs' element in the 's2dv_cube' object. It consists of a list 
+containing two elements, the start and end dates of the aggregated period with 
+the same dimensions of 'Dates' element.
 }
 \description{
 The Total Time of a variable exceeding (or not) a Threshold. It returns the  
@@ -83,7 +89,7 @@ indices for heat stress can be obtained by using this function:
 }
 \examples{
 exp <- NULL
-exp$data <- array(abs(rnorm(5 * 3 * 214 * 2)*280),
+exp$data <- array(rnorm(5 * 3 * 214 * 2)*23,
                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
 exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
@@ -91,7 +97,9 @@ exp$attrs$Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"),
                         as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
                     seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
                         as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(exp$attrs$Dates) <- c(sdate = 3, ftime = 214)
 class(exp) <- 's2dv_cube'
-DOT <- CST_TotalTimeExceedingThreshold(exp, threshold = 280)
+DOT <- CST_TotalTimeExceedingThreshold(exp, threshold = 23, start = list(21, 4), 
+                                      end = list(21, 6))
 
 }
diff --git a/man/CST_WindCapacityFactor.Rd b/man/CST_WindCapacityFactor.Rd
index 638f5b858ff97d9dcc7302bfda05c49a5b921959..84c057d3277410cddfce394fc6696dd49eb8f1c8 100644
--- a/man/CST_WindCapacityFactor.Rd
+++ b/man/CST_WindCapacityFactor.Rd
@@ -60,14 +60,23 @@ below).
 \examples{
 wind <- NULL
 wind$data <- array(rweibull(n = 100, shape = 2, scale = 6), 
-                  c(member = 10, lat = 2, lon = 5))
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2, lat = 5))
 wind$coords <- list(lat =  c(40, 41), lon = 1:5)
 variable <- list(varName = 'sfcWind', 
                 metadata = list(sfcWind = list(level = 'Surface')))
 wind$attrs <- list(Variable = variable, Datasets = 'synthetic', 
                   when = Sys.time(), Dates = '1990-01-01 00:00:00')
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+                    seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+                    seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+wind$attrs$Dates <- Dates
 class(wind) <- 's2dv_cube'
-WCF <- CST_WindCapacityFactor(wind, IEC_class = "III")
+WCF <- CST_WindCapacityFactor(wind, IEC_class = "III", 
+                             start = list(21, 4), end = list(21, 6))
 
 }
 \references{
diff --git a/man/CST_WindPowerDensity.Rd b/man/CST_WindPowerDensity.Rd
index c33bd8d99a3cb0f0bb37e1872e962f1878ab7674..4b04aed2026256d78dcf1c5798c60b0160d5f2e1 100644
--- a/man/CST_WindPowerDensity.Rd
+++ b/man/CST_WindPowerDensity.Rd
@@ -14,8 +14,8 @@ CST_WindPowerDensity(
 )
 }
 \arguments{
-\item{wind}{An s2dv_cube object with instantaneous wind speeds expressed in m/s 
-obtained from CST_Load or s2dv_cube functions from CSTools pacakge.}
+\item{wind}{An 's2dv_cube' object with instantaneous wind speeds expressed in 
+m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge.}
 
 \item{ro}{A scalar, or alternatively a multidimensional array with the same
 dimensions as wind, with the air density expressed in kg/m^3. By default it
@@ -53,14 +53,23 @@ it will give inaccurate results if used with period means.
 \examples{
 wind <- NULL
 wind$data <- array(rweibull(n = 100, shape = 2, scale = 6), 
-                  c(member = 10, lat = 2, lon = 5))
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2, lat = 5))
 wind$coords <- list(lat =  c(40, 41), lon = 1:5)
 variable <- list(varName = 'sfcWind', 
                 metadata = list(sfcWind = list(level = 'Surface')))
 wind$attrs <- list(Variable = variable, Datasets = 'synthetic', 
                   when = Sys.time(), Dates = '1990-01-01 00:00:00')
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+                    seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+                    seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+                        as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+wind$attrs$Dates <- Dates
 class(wind) <- 's2dv_cube'
-WCF <- CST_WindPowerDensity(wind)
+WPD <- CST_WindPowerDensity(wind, start = list(21, 4), 
+                           end = list(21, 6))
 
 }
 \author{
diff --git a/man/PeriodAccumulation.Rd b/man/PeriodAccumulation.Rd
index 99033211af06dc55c429d728d69c2fbee910c756..614b65cbe112398eb62975b4c76495dd6709e859 100644
--- a/man/PeriodAccumulation.Rd
+++ b/man/PeriodAccumulation.Rd
@@ -17,9 +17,9 @@ PeriodAccumulation(
 \arguments{
 \item{data}{A multidimensional array with named dimensions.}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -64,14 +64,16 @@ exp <- array(rnorm(216)*200, dim = c(dataset = 1, member = 2, sdate = 3,
 TP <- PeriodAccumulation(exp, time_dim = 'ftime')
 data <- array(rnorm(5 * 3 * 214 * 2),
              c(memb = 5, sdate = 3, time = 214, lon = 2)) 
-# ftime tested
 Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
               as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
           seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
               as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
           seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
               as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
-SprR <- PeriodAccumulation(data, dates = Dates, start = list(21, 4), end = list(21, 6))
-HarR <- PeriodAccumulation(data, dates = Dates, start = list(21, 8), end = list(21, 10))
+dim(Dates) <- c(sdate = 3, time = 214)
+SprR <- PeriodAccumulation(data, dates = Dates, start = list(21, 4), 
+                          end = list(21, 6))
+HarR <- PeriodAccumulation(data, dates = Dates, start = list(21, 8), 
+                          end = list(21, 10))
 
 }
diff --git a/man/PeriodMean.Rd b/man/PeriodMean.Rd
index fffb332299a35a3a67e7fe274059628d592e7d5a..9637d58f54a4cc5b2362438d348cc44b9f4accf0 100644
--- a/man/PeriodMean.Rd
+++ b/man/PeriodMean.Rd
@@ -17,9 +17,9 @@ PeriodMean(
 \arguments{
 \item{data}{A multidimensional array with named dimensions.}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -33,7 +33,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the dimension to 
-compute the indicator. By default, it is set to 'ftime'. More than one 
+compute the indicator. By default, it is set to 'time'. More than one 
 dimension name matching the dimensions provided in the object 
 \code{data$data} can be specified.}
 
@@ -59,7 +59,16 @@ this function:
 }
 }
 \examples{
-exp <- array(rnorm(56), dim = c(member = 7, ftime = 8))
-SA <- PeriodMean(exp, time_dim = 'ftime')
+data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+SA <- PeriodMean(data, dates = Dates, start = list(01, 12), end = list(01, 01))
 
 }
diff --git a/man/QThreshold.Rd b/man/QThreshold.Rd
index 2af6e5f93df4ae8e374bdb6211f531d94546ff47..efc48cf5c844fc53977430a22758d2dab177e968 100644
--- a/man/QThreshold.Rd
+++ b/man/QThreshold.Rd
@@ -23,9 +23,9 @@ QThreshold(
 units as parameter 'data' and with the common dimensions of the element 
 'data' of the same length.}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -39,7 +39,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the temporal 
-dimension. By default, it is set to 'ftime'. More than one dimension name 
+dimension. By default, it is set to 'time'. More than one dimension name 
 matching the dimensions provided in the object \code{data$data} can be 
 specified. This dimension is required to subset the data in a requested 
 period.}
@@ -82,7 +82,17 @@ and memb_dim parameters:
 \examples{
 threshold = 25
 data <- array(rnorm(5 * 3 * 20 * 2, mean = 26), 
-             c(member = 5, sdate = 3, time = 20, lon = 2)) 
-thres_q <- QThreshold(data, threshold)
+             c(member = 5, sdate = 3, time = 214, lon = 2)) 
+
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, time = 214)
+
+thres_q <- QThreshold(data, threshold, dates = Dates, time_dim = 'time',
+                     start = list(21, 4), end = list(21, 6))
 
 }
diff --git a/man/Threshold.Rd b/man/Threshold.Rd
index db5981753e68504c918674c79c85f87fd2847450..dc9d2a24426af50029823aeeb1f9944df42c66d9 100644
--- a/man/Threshold.Rd
+++ b/man/Threshold.Rd
@@ -21,11 +21,11 @@ Threshold(
 \item{data}{A multidimensional array with named dimensions.}
 
 \item{threshold}{A single scalar or vector indicating the relative 
-threshold(s).}
+threshold(s). It must contain values between 0 and 1.}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -39,7 +39,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the temporal 
-dimension. By default, it is set to 'ftime'. More than one dimension name 
+dimension. By default, it is set to 'time'. More than one dimension name 
 matching the dimensions provided in the object \code{data$data} can be 
 specified. This dimension is required to subset the data in a requested 
 period.}
diff --git a/man/TotalSpellTimeExceedingThreshold.Rd b/man/TotalSpellTimeExceedingThreshold.Rd
index 276423bb14092674ee655377d78a89e3a9f8796a..10124de096bb539121ccb1db3444656f279be839 100644
--- a/man/TotalSpellTimeExceedingThreshold.Rd
+++ b/man/TotalSpellTimeExceedingThreshold.Rd
@@ -36,9 +36,9 @@ are used it has to be a vector of a pair of two logical operators:
 c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
 c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -86,9 +86,19 @@ different behaviour consider to modify the 'data' input by substituting NA
 values by values exceeding the threshold.
 }
 \examples{
-data <- array(rnorm(120), c(member = 1, sdate = 2, ftime = 20, lat = 4))
-threshold <- array(rnorm(4), c(lat = 4))
-total <- TotalSpellTimeExceedingThreshold(data, threshold, spell = 6)
+data <- array(1:100, c(member = 5, sdate = 3, ftime = 214, lon = 2))
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+
+threshold <- array(1:4, c(lat = 4))
+total <- TotalSpellTimeExceedingThreshold(data, threshold, dates = Dates, 
+                                         spell = 6, start = list(21, 4), 
+                                         end = list(21, 6))
 
 }
 \seealso{
diff --git a/man/TotalTimeExceedingThreshold.Rd b/man/TotalTimeExceedingThreshold.Rd
index 206847574f25df8eabde09914d03bbd2d116bc46..4dc00d019c6c31f60214ebb963ab6963bfceaf43 100644
--- a/man/TotalTimeExceedingThreshold.Rd
+++ b/man/TotalTimeExceedingThreshold.Rd
@@ -34,9 +34,9 @@ are used it has to be a vector of a pair of two logical operators:
 c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
 c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is 
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to define the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
@@ -62,7 +62,8 @@ computation.}
 \value{
 A multidimensional array with named dimensions containing the total 
 number of the corresponding units of the data frequency that a variable is 
-exceeding a threshold during a period.
+exceeding a threshold during a period with dimensions of the input parameter 
+'data' except the dimension where the indicator has been computed.
 }
 \description{
 The Total Time of a variable exceeding (or not) a Threshold. It returns the  
@@ -86,8 +87,16 @@ indices for heat stress can be obtained by using this function:
 }
 }
 \examples{
-exp <- array(abs(rnorm(5 * 3 * 214 * 2)*280),
-            c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
-DOT <- TotalTimeExceedingThreshold(exp, threshold = 300, time_dim = 'ftime')
+data <- array(rnorm(5 * 3 * 214 * 2)*23,
+             c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, ftime = 214)
+DOT <- TotalTimeExceedingThreshold(data, threshold = 23, dates = Dates,
+                                  start = list(21, 4), end = list(21, 6))
 
 }
diff --git a/man/WindCapacityFactor.Rd b/man/WindCapacityFactor.Rd
index 69549a817e3ebd951e3e47e553205699276f95ce..0b6b958109a22efced98cb67b78459c580fed689 100644
--- a/man/WindCapacityFactor.Rd
+++ b/man/WindCapacityFactor.Rd
@@ -25,9 +25,9 @@ respectively. Classes \code{'I/II'} and \code{'II/III'} indicate
 intermediate turbines that fit both classes. More details of the five
 turbines and a plot of its power curves can be found in Lledó et al. (2019).}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period
 to select from the data by providing a list of two elements: the initial
@@ -41,7 +41,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the dimension to
-compute the indicator. By default, it is set to 'ftime'. More than one
+compute the indicator. By default, it is set to 'time'. More than one
 dimension name matching the dimensions provided in the object
 \code{data$data} can be specified.}
 
@@ -65,8 +65,19 @@ different power curves that span different IEC classes can be selected (see
 below).
 }
 \examples{
-wind <- rweibull(n = 100, shape = 2, scale = 6)
-WCF <- WindCapacityFactor(wind, IEC_class = "III")
+wind <- array(rweibull(n = 32100, shape = 2, scale = 6), 
+             c(member = 5, sdate = 3, time = 214, lon = 2, lat = 5))
+
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, time = 214)
+
+WCF <- WindCapacityFactor(wind, IEC_class = "III", dates = Dates, 
+                         start = list(21, 4), end = list(21, 6))
 
 }
 \references{
diff --git a/man/WindPowerDensity.Rd b/man/WindPowerDensity.Rd
index 8e3c8e3d3d147ad0b3cb88c593ad6f390c51e189..9ca3234c3005f61683d6896e1825346179cac9ce 100644
--- a/man/WindPowerDensity.Rd
+++ b/man/WindPowerDensity.Rd
@@ -22,9 +22,9 @@ speeds expressed in m/s.}
 dimensions as wind, with the air density expressed in kg/m^3. By default it
 takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.}
 
-\item{dates}{A vector of dates or a multidimensional array of dates with named
-dimensions matching the dimensions on parameter 'data'. By default it is
-NULL, to select a period this parameter must be provided.}
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
 
 \item{start}{An optional parameter to defined the initial date of the period
 to select from the data by providing a list of two elements: the initial
@@ -38,7 +38,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the dimension to
-compute the indicator. By default, it is set to 'ftime'. More than one
+compute the indicator. By default, it is set to 'time'. More than one
 dimension name matching the dimensions provided in the object
 \code{data$data} can be specified.}
 
@@ -57,8 +57,17 @@ It is computed as 0.5*ro*wspd^3. As this function is non-linear,
 it will give inaccurate results if used with period means.
 }
 \examples{
-wind <- rweibull(n = 100, shape = 2, scale = 6)
-WPD <- WindPowerDensity(wind)
+wind <- array(rweibull(n = 32100, shape = 2, scale = 6), 
+             c(member = 5, sdate = 3, time = 214, lon = 2, lat = 5))
+Dates <- c(seq(as.Date("01-05-2000", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2000", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2001", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2001", format = "\%d-\%m-\%Y"), by = 'day'),
+          seq(as.Date("01-05-2002", format = "\%d-\%m-\%Y"), 
+              as.Date("30-11-2002", format = "\%d-\%m-\%Y"), by = 'day'))
+dim(Dates) <- c(sdate = 3, time = 214)
+WPD <- WindPowerDensity(wind, dates = Dates, start = list(21, 4), 
+                       end = list(21, 6))
 
 }
 \author{
diff --git a/tests/testthat/test-AccumulationExceedingThreshold.R b/tests/testthat/test-AccumulationExceedingThreshold.R
index 0331e8efae5e2f444bed43ccf3ca6891c93fd9c1..baa2c50b2936287d1d901159b05a25ae10279079 100644
--- a/tests/testthat/test-AccumulationExceedingThreshold.R
+++ b/tests/testthat/test-AccumulationExceedingThreshold.R
@@ -1,4 +1,7 @@
 ##############################################
+
+library(CSTools)
+
 # dat1
 dat1 <- 1:20
 
@@ -234,7 +237,6 @@ test_that("4. Output checks", {
 })
 
 ##############################################
-library(CSTools)
 
 test_that("5. Seasonal forecasts", {
   exp <- CSTools::lonlat_temp$exp
@@ -260,24 +262,25 @@ test_that("5. Seasonal forecasts", {
                  as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
              seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
                  as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+  dim(Dates) <- c(sdate = 3, ftime = 214)
   GDD <- AccumulationExceedingThreshold(exp - 17, threshold = 0, dates = Dates, time_dim = 'ftime', 
                                         start = list(1, 4), end = list(31, 10), na.rm = TRUE)
   
   expect_equal(
     round(GDD[,1,1,1]),
-    c(538, 367, 116, 519, 219, 282)
+    c(549, 387, 125, 554, 245, 282)
   )
   expect_equal(
     dim(GDD),
-    c(member = 6, sdate = 3, lat =4, lon = 4)
+    c(member = 6, sdate = 3, lat = 4, lon = 4)
   )
   expect_error(
     AccumulationExceedingThreshold(exp - 17, threshold = 0, dates = Dates, start = list(1, 4), end = list(31, 10), time_dim = 'time'),
     "Parameter 'time_dim' is not found in 'data' dimension."
   )
   expect_equal(
-    all(is.na(AccumulationExceedingThreshold(exp - 17, threshold = 0, dates = Dates, time_dim = 'ftime',start = list(1, 4), end = list(31, 10)))),
-    all(is.na(c(NA, NA)))
+    !any(is.na(AccumulationExceedingThreshold(exp - 17, threshold = 0, dates = Dates, time_dim = 'ftime',start = list(1, 4), end = list(31, 10)))),
+    !any(is.na(c(1, 1)))
   )
 
   # test the 'diff'
@@ -294,7 +297,6 @@ test_that("5. Seasonal forecasts", {
     AccumulationExceedingThreshold(input_1, threshold_1),
     204
   )
-
   expect_equal(
     AccumulationExceedingThreshold(input_2, threshold_2, op = '<'),
     -105
@@ -304,3 +306,47 @@ test_that("5. Seasonal forecasts", {
     -55
   )
 })
+
+##############################################
+test_that("6. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_AccumulationExceedingThreshold(data = exp, threshold = mean(exp$data),
+                                            time_dim = 'ftime', start = list(10, 03), 
+                                            end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
diff --git a/tests/testthat/test-PeriodAccumulation.R b/tests/testthat/test-PeriodAccumulation.R
index 81718f5d9337fc6f39158bcc8d94a8d5e0a02d0f..6898a937be04267b5bb4b0deb708d11fe080468b 100644
--- a/tests/testthat/test-PeriodAccumulation.R
+++ b/tests/testthat/test-PeriodAccumulation.R
@@ -1,5 +1,7 @@
+library(CSTools)
+
 ##############################################
-test_that("Sanity Checks", {
+test_that("1. Sanity Checks", {
   expect_error(
     PeriodAccumulation('x'), 
     "Parameter 'data' must be numeric."
@@ -25,18 +27,33 @@ test_that("Sanity Checks", {
     PeriodAccumulation(1:10), 
     55
   )
-  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
+  data <- array(1:24, c(sdate = 2, ftime = 3, lon = 4))
   expect_equal(
     PeriodAccumulation(data),
     array(c(9, 12, 27, 30, 45, 48, 63, 66), c(sdate = 2, lon = 4))
   )
-
+  # Test dates warning
+  expect_warning(
+    PeriodAccumulation(array(1:10, c(ftime = 10)), 
+                       dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                                   as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
+                       start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and 'end' ", 
+           "are not NULL. All data will be used.")
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    PeriodAccumulation(array(1:10, c(ftime = 10)), 
+                       start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
+
 ##############################################
-library(CSTools)
 
-test_that("seasonal", {
-  exp <- CSTools::lonlat_prec
+test_that("2. Seasonal", {
+  exp <- NULL
   exp$data <- array(1:(1 * 3 * 214 * 2),
                     c(memb = 1, sdate = 3, ftime = 214, lon = 2))
   exp$dims <- dim(exp$data)
@@ -47,6 +64,7 @@ test_that("seasonal", {
                        seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
                            as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
   dim(exp$attrs$Dates) <- c(ftime = 214, sdate = 3)
+  class(exp) <- 's2dv_cube'
   output <- exp
   output$data <- array(c(sum(exp$data[1,1,21:82,1]), sum(exp$data[1,2,21:82,1]),
                          sum(exp$data[1,3,21:82,1]), sum(exp$data[1,1,21:82,2]),
@@ -58,3 +76,60 @@ test_that("seasonal", {
     output$data
   )
 })
+
+##############################################
+
+test_that("3. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodAccumulation(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                                start = list(10, 03), end = list(20, 03))
+  res2 <- CST_PeriodAccumulation(data = CSTools::lonlat_prec, time_dim = 'ftime')
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  expect_equal(
+    dim(res2$data), 
+    dim(exp$data)[-which(names(dim(exp$data)) == 'ftime')]
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  expect_equal(
+    dim(res2$data)['sdate'], 
+    dim(res2$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  expect_equal(
+    res2$attrs$Dates, 
+    res2$attrs$time_bounds$start
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
\ No newline at end of file
diff --git a/tests/testthat/test-PeriodMean.R b/tests/testthat/test-PeriodMean.R
index 41866399e21aa251dfd3018607f6219ef7de85a0..cd9f5fe5ddca8b82c54fbbbc661bf68e0e8490a0 100644
--- a/tests/testthat/test-PeriodMean.R
+++ b/tests/testthat/test-PeriodMean.R
@@ -1,5 +1,7 @@
+library(CSTools)
+
 ##############################################
-test_that("Sanity Checks", {
+test_that("1. Sanity Checks", {
   expect_error(
     PeriodMean('x'), 
     "Parameter 'data' must be numeric."
@@ -16,14 +18,15 @@ test_that("Sanity Checks", {
   )
   expect_error(
     PeriodMean(1, dates = '2000-01-01', end = 3, start = 4),
-    "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end.")
+    "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end."
+  )
   suppressWarnings(
     expect_equal(
-      PeriodMean(array(1:10, c(time = 10))), 
+      PeriodMean(array(1:10, c(ftime = 10))), 
       5.5
     )
   )
-  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
+  data <- array(1:24, c(sdate = 2, ftime = 3, lon = 4))
   suppressWarnings(
     expect_equal(
       PeriodMean(data),
@@ -31,13 +34,28 @@ test_that("Sanity Checks", {
       c(sdate = 2, lon = 4))
     )
   )
+  # Test dates warning
+  expect_warning(
+    PeriodMean(array(1:10, c(ftime = 10)), 
+               dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                           as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and 'end' ", 
+           "are not NULL. All data will be used.")
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    PeriodMean(array(1:10, c(ftime = 10)), 
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
 ##############################################
-library(CSTools)
 
-test_that("seasonal", {
-  exp <- CSTools::lonlat_prec
+test_that("2. Seasonal", {
+  exp <- NULL
   exp$data <- array(1:(1 * 3 * 214 * 2),
                     c(memb = 1, sdate = 3, ftime = 214, lon = 2))
   exp$dims <- dim(exp$data)
@@ -48,6 +66,7 @@ test_that("seasonal", {
                        seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
                            as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
   dim(exp$attrs$Dates) <- c(ftime = 214, sdate = 3)
+  class(exp) <- 's2dv_cube'
   output <- exp
   output$data <- array(c(mean(exp$data[1,1,21:82,1]), mean(exp$data[1,2,21:82,1]),
                          mean(exp$data[1,3,21:82,1]), mean(exp$data[1,1,21:82,2]),
@@ -58,3 +77,47 @@ test_that("seasonal", {
     output$data
   )
 })
+
+##############################################
+test_that("3. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodMean(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
+
diff --git a/tests/testthat/test-QThreshold.R b/tests/testthat/test-QThreshold.R
index deb35df4d97d1263568b621cfac5d46ae4d63795..57883b15a8286fbbce6165e34f8ee07b092acaf6 100644
--- a/tests/testthat/test-QThreshold.R
+++ b/tests/testthat/test-QThreshold.R
@@ -1,5 +1,21 @@
+
+library(CSTools)
+
+# dat1
+threshold <- 26
+dat1 <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+dates0 <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+dates1 <- dates0
+dim(dates1) <- c(sdate = 3, ftime = 214)
+
 ##############################################
-test_that("Sanity checks", {
+test_that("1. Sanity checks", {
   expect_error(
     QThreshold(NULL),
     "Parameter 'data' cannot be NULL."
@@ -104,13 +120,30 @@ test_that("Sanity checks", {
     dim(res),
     c(sdate = 3, ftime = 52)
   )
-
+  # test start and end
+  expect_warning(
+    QThreshold(dat1, threshold = 26, dates = dates0, start = list(21, 4), 
+               end = list(21, 6)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and ",
+           "'end' are not NULL. All data will be used.")
+  )
+  expect_equal(
+    dim(QThreshold(dat1, threshold = 26, dates = dates1, start = list(21, 4), 
+                   end = list(21, 6))),
+    c(sdate = 3, member = 5, ftime = 52, lon = 2)
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    QThreshold(array(1:61, dim = c(ftime = 61, sdate = 3)), threshold = 25,
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
 ##############################################
-library(CSTools)
 
-test_that("Seasonal forecasts", {
+test_that("2. Seasonal forecasts", {
   obs <- CSTools::lonlat_temp$obs$data - 248
   obs_percentile <- QThreshold(obs, threshold = 35)
   expect_equal(
@@ -142,3 +175,5 @@ test_that("Seasonal forecasts", {
     c(0.6, 0.4, 0.6, 0.6, 0.4, 0.4)
   )
 })
+
+##############################################
diff --git a/tests/testthat/test-Threshold.R b/tests/testthat/test-Threshold.R
index 84e4e0d02ba41aa9bc488d3e83a5c4c45c0ca315..7cb83cff6a06e375a5f96fafb1ae10b4a73fb1c1 100644
--- a/tests/testthat/test-Threshold.R
+++ b/tests/testthat/test-Threshold.R
@@ -1,5 +1,21 @@
+
+library(CSTools)
+
+# dat1
+threshold <- 0.9
+dat1 <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+dates0 <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+dates1 <- dates0
+dim(dates1) <- c(sdate = 3, ftime = 214)
+
 ##############################################
-test_that("Sanity checks", {
+test_that("1. Sanity checks", {
   expect_error(
     Threshold(NULL),
     "Parameter 'data' cannot be NULL."
@@ -60,9 +76,30 @@ test_that("Sanity checks", {
     dim(Threshold(data, threshold, memb_dim = 'x', sdate_dim = 'ftime')),
     c(probs = 2)
   )
+  # test start and end
+  expect_warning(
+    Threshold(dat1, threshold = 0.9, dates = dates0, start = list(21, 4), 
+              end = list(21, 6)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and ",
+           "'end' are not NULL. All data will be used.")
+  )
+  expect_equal(
+    dim(Threshold(dat1, threshold = 0.8, dates = dates1, start = list(21, 4), 
+                  end = list(21, 6))),
+    c(ftime = 52, lon = 2)
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    Threshold(array(1:366, dim = c(ftime = 61, sdate = 3, member = 2)), 
+              threshold = 0.8, start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
-test_that("Seasonal forecasts", {
+##############################################
+
+test_that("2. Seasonal forecasts", {
   exp <- CSTools::lonlat_temp$exp$data
   thresholdP <- Threshold(exp, threshold = 0.9)
   expect_equal(
diff --git a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
index 36e46cf88e3a199019c5a7edfffd6b60a9ec8e94..b66c5ae745d81910e802b4092678f985a93bf5ac 100644
--- a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
@@ -142,6 +142,13 @@ test_that("1. Sanity checks", {
     paste0("Parameter 'start' and 'end' must be lists indicating the ",
            "day and the month of the period start and end.")
   )
+  # start and end when dates is not provided
+  expect_warning(
+    TotalSpellTimeExceedingThreshold(array(1:10, c(ftime = 10)), threshold = 5, spell = 2, 
+                                     start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
 ###########################################################################
@@ -240,7 +247,7 @@ test_that("4. Output checks", {
 
 })
 
-###########################################################################
+##################################################################
 
 test_that("5. Seasonal Forecasts", {
   exp <- CSTools::lonlat_temp$exp
@@ -263,3 +270,47 @@ test_that("5. Seasonal Forecasts", {
     WSDI1$data[3,3,3,],
     c(rep(0, 53)))
 })
+
+##################################################################
+
+test_that("6. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodMean(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
\ No newline at end of file
diff --git a/tests/testthat/test-TotalTimeExceedingThreshold.R b/tests/testthat/test-TotalTimeExceedingThreshold.R
index 1c3ca19d6e1fbf82e4c212d6eebea2f093c3c3e2..c025c3ed844cd99701e8053b36e4ec3fbb64a7e8 100644
--- a/tests/testthat/test-TotalTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalTimeExceedingThreshold.R
@@ -40,7 +40,7 @@ thres4_7 <- rnorm(5, 20)
 set.seed(2)
 thres4_8 <- rnorm(5, 25)
 
-###########################################################################
+#############################################################
 
 test_that("1. Sanity checks", {
   # data
@@ -134,10 +134,17 @@ test_that("1. Sanity checks", {
     paste0("Parameter 'start' and 'end' must be lists indicating the ",
            "day and the month of the period start and end.")
   )
+  # start and end when dates is not provided
+  expect_warning(
+    TotalTimeExceedingThreshold(array(1:10, c(ftime = 10)), threshold = 5,  
+                                start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
 
-###########################################################################
+#######################################################
 
 test_that("2. Output checks", {
   expect_equal(
@@ -231,7 +238,7 @@ test_that("4. Output checks", {
 
 ###########################################################################
 
-test_that("Seasonal forecasts", {
+test_that("5. Seasonal forecasts", {
   # compare with scalar fixed threshold
   exp <- CSTools::lonlat_temp$exp
   obs <- CSTools::lonlat_temp$obs
@@ -252,3 +259,47 @@ test_that("Seasonal forecasts", {
     c(3, 3, 3, 3, 3, 3)
   )
 })
+
+##################################################################
+
+test_that("6. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodMean(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
diff --git a/tests/testthat/test-WindCapacityFactor.R b/tests/testthat/test-WindCapacityFactor.R
index bcecdbce63a017d91e05c5326168d806b01ab334..3afa27e4de0e8b9b8074b084851937c17a88b40e 100644
--- a/tests/testthat/test-WindCapacityFactor.R
+++ b/tests/testthat/test-WindCapacityFactor.R
@@ -12,7 +12,19 @@ wind$attrs <- list(Variable = variable, Datasets = 'synthetic',
 class(wind) <- 's2dv_cube'
 WCF <- CST_WindCapacityFactor(wind, IEC_class = "III")
 
-###########################################################################
+# dat2
+dat2 <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+dates0 <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+dates2 <- dates0
+dim(dates2) <- c(sdate = 3, ftime = 214)
+
+###################################################
 test_that("1. Input checks", {
   # Check 's2dv_cube'
   expect_error(
@@ -33,9 +45,16 @@ test_that("1. Input checks", {
     paste0("Parameter 'start' and 'end' must be lists indicating the ",
            "day and the month of the period start and end.")
   )
+  # start and end when dates is not provided
+  expect_warning(
+    WindCapacityFactor(array(1:10, c(ftime = 10)), 
+                       start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
-###########################################################################
+####################################################
 test_that("2. Output checks", {
   expect_equal(
     CST_WindCapacityFactor(wind = wind)$attrs$Variable$varName, 
@@ -45,5 +64,17 @@ test_that("2. Output checks", {
     dim(CST_WindCapacityFactor(wind = wind)$data), 
     c(member = 10, lat = 2, lon = 5)
   )
+  # test start and end
+  expect_warning(
+    WindCapacityFactor(wind = dat2, IEC_class = "III", dates = dates0, 
+                       start = list(21, 4), end = list(21, 6)),
+    paste0("Parameter 'wind' must have named dimensions if 'start' and ",
+           "'end' are not NULL. All data will be used.")
+  )
+  expect_equal(
+    dim(WindCapacityFactor(wind = dat2, IEC_class = "III", dates = dates2, 
+                           start = list(21, 4), end = list(21, 6))),
+    c(member = 5, sdate = 3, ftime = 52, lon = 2)
+  )
 })
 
diff --git a/tests/testthat/test-WindPowerDensity.R b/tests/testthat/test-WindPowerDensity.R
index 184b06296759fa49dd237b001812fb63f8aecf9a..999235a3efe93b8bffb550ddf25ad55c84373581 100644
--- a/tests/testthat/test-WindPowerDensity.R
+++ b/tests/testthat/test-WindPowerDensity.R
@@ -11,6 +11,18 @@ wind$attrs <- list(Variable = variable, Datasets = 'synthetic',
                    when = Sys.time(), Dates = '1990-01-01 00:00:00')
 class(wind) <- 's2dv_cube'
 
+# dat2
+dat2 <- array(abs(rnorm(5 * 3 * 214 * 2)*50),
+                  c(member = 5, sdate = 3, ftime = 214, lon = 2)) 
+dates0 <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2000", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2001", format = "%d-%m-%Y"), by = 'day'),
+                     seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+                         as.Date("30-11-2002", format = "%d-%m-%Y"), by = 'day'))
+dates2 <- dates0
+dim(dates2) <- c(sdate = 3, ftime = 214)
+
 ###########################################################################
 test_that("1. Input checks", {
   # Check 's2dv_cube'
@@ -32,6 +44,13 @@ test_that("1. Input checks", {
     paste0("Parameter 'start' and 'end' must be lists indicating the ",
            "day and the month of the period start and end.")
   )
+  # start and end when dates is not provided
+  expect_warning(
+    WindPowerDensity(array(1:10, c(ftime = 10)), 
+                     start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
 })
 
 ###########################################################################
@@ -44,5 +63,17 @@ test_that("2. Output checks", {
     dim(CST_WindPowerDensity(wind = wind)$data), 
     c(member = 10, lat = 2, lon = 5)
   )
+  # test start and end
+  expect_warning(
+    WindPowerDensity(wind = dat2, dates = dates0, 
+                     start = list(21, 4), end = list(21, 6)),
+    paste0("Parameter 'wind' must have named dimensions if 'start' and ",
+           "'end' are not NULL. All data will be used.")
+  )
+  expect_equal(
+    dim(WindPowerDensity(wind = dat2, dates = dates2, 
+                         start = list(21, 4), end = list(21, 6))),
+    c(member = 5, sdate = 3, ftime = 52, lon = 2)
+  )
 })