diff --git a/R/AccumulationExceedingThreshold.R b/R/AccumulationExceedingThreshold.R
index f202dca791a07a96d523e8b5680f80159138b00e..5cb70f354e3df42ab95d56488a42a12b895248c0 100644
--- a/R/AccumulationExceedingThreshold.R
+++ b/R/AccumulationExceedingThreshold.R
@@ -12,47 +12,55 @@
 #'                  temperatures and 10°C between April 1st and October 31st}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in 
-#'  package CSTools.
-#'@param threshold An 's2dv_cube' object as output of a 'CST_' function in the
-#'  same units as parameter 'data' and with the common dimensions of the element 
-#'  'data' of the same length. A single scalar is also possible.
-#'@param op An operator '>' (by default), '<', '>=' or '<='. 
-#'@param diff A logical value indicating whether to accumulate the difference
-#'  between data and threshold (TRUE) or not (FALSE by default).
-#'@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
-#'  to NULL and the indicator is computed using all the data provided in
+#'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in
+#'  package CSTools. 
+#'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
+#'  object or a multidimensional array with named dimensions. It must be in the 
+#'  same units and with the common dimensions of the same length as parameter 
+#'  'data'. It can also be a vector with the same length of 'time_dim' from 
+#'  'data' or a scalar. If we want to use two thresholds: it can be a vector 
+#'  of two scalars, a list of two vectors with the same length of 
+#'  'time_dim' from 'data' or a list of two multidimensional arrays with the 
+#'  common dimensions of the same length as parameter 'data'. If two thresholds
+#'  are used, parameter 'op' must be also a vector of two elements.
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  are used it has to be a vector of a pair of two logical operators: 
+#'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+#'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
+#'@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 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 to NULL and the indicator is computed using all the data provided in 
 #'  \code{data}.
-#'@param end An optional parameter to defined the final date of the period to
-#'  select from the data by providing a list of two elements: the final day of
+#'@param end An optional parameter to define the final date of the period to 
+#'  select from the data by providing a list of two elements: the final day of 
 #'  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
-#'  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
-#'  not (FALSE). 
+#'@param time_dim A character string indicating the name of the dimension to 
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
+#'@param na.rm A logical value indicating whether to ignore NA values (TRUE) 
+#'  or not (FALSE).
 #'@param ncores An integer indicating the number of cores to use in parallel
-#'  computation.
+#'computation.
 #'
-#'@return A 's2dv_cube' object containing the indicator in the element \code{data}.
-#'
-#'@examples
+#'@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.
+#'@examples 
 #'exp <- NULL
 #'exp$data <- array(rnorm(216)*200, dim = c(dataset = 1, member = 2, sdate = 3, 
 #'                  ftime = 9, lat = 2, lon = 2))
 #'class(exp) <- 's2dv_cube'
 #'DOT <- CST_AccumulationExceedingThreshold(exp, threshold = 280)
-#' 
+#'
 #'@import multiApply
 #'@export
-CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', 
-                                               diff = FALSE,
-                                               start = NULL, end = NULL,
-                                               time_dim = 'ftime',
+CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FALSE,
+                                               start = NULL, end = NULL, time_dim = 'ftime',
                                                na.rm = FALSE, ncores = NULL) {
   if (!inherits(data, 's2dv_cube')) {
     stop("Parameter 'data' must be of the class 's2dv_cube', ",
@@ -62,29 +70,39 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>',
   if (!is.null(start) && !is.null(end)) {
     if (is.null(dim(data$Dates$start))) {
       if (length(data$Dates$start) != dim(data$data)[time_dim]) {
-        if (length(data$Dates$start) == 
+        if (length(data$Dates$start) ==
             prod(dim(data$data)[time_dim] * dim(data$data)['sdate'])) {
           dim(data$Dates$start) <- c(dim(data$data)[time_dim],
                                      dim(data$data)['sdate'])
-        } else {
-          warning("Dimensions in 'data' element 'Dates$start' are missed and ",
-                  "all data would be used.")
         }
+      } else {
+        warning("Dimensions in 'data' element 'Dates$start' are missed and",
+                "all data would be used.")
       }
     }
   }
- if (inherits(threshold, 's2dv_cube')) {
-    threshold <- threshold$data
+  if (length(op) == 1) {
+    if (inherits(threshold, 's2dv_cube')) {
+        threshold <- threshold$data
+    }
+  } else if (length(op) == 2) {
+    if (inherits(threshold[[1]], 's2dv_cube')) {
+      threshold[[1]] <- threshold[[1]]$data
+    }
+    if (inherits(threshold[[2]], 's2dv_cube')) {
+      threshold[[2]] <- threshold[[2]]$data
+    }
   }
- total <- AccumulationExceedingThreshold(data$data, data$Dates[[1]],
-                                         threshold = threshold, op = op, diff = diff,
-                                         start = start, end = end, time_dim = time_dim,
-                                         na.rm = na.rm, ncores = ncores)
+
+  total <- AccumulationExceedingThreshold(data$data, dates = data$Dates[[1]],
+                                          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$Dates <- SelectPeriodOnDates(dates = data$Dates$start,
-                                       start = start, end = end,
-                                       time_dim = time_dim, ncores = ncores)
+    data$Dates <- SelectPeriodOnDates(dates = data$Dates$start,
+                                      start = start, end = end,
+                                      time_dim = time_dim, ncores = ncores)
   }
   return(data)
 }
@@ -103,41 +121,50 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>',
 #'}
 #'
 #'@param data A multidimensional array with named dimensions.
-#'@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 op An operator '>' (by default), '<', '>=' or '<='.
-#'@param diff A logical value indicating whether to accumulate the difference
-#'  between data and threshold (TRUE) or not (FALSE by default).
-#'@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 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
-#'  to NULL and the indicator is computed using all the data provided in
+#'@param threshold If only one threshold is used: it can be a multidimensional 
+#'  array with named dimensions. It must be in the same units and with the 
+#'  common dimensions of the same length as parameter 'data'. It can also be a
+#'  vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+#'  want to use two thresholds: it can be a vector of two scalars, a list of 
+#'  two vectors with the same length of 'time_dim' from 'data' or a list of 
+#'  two multidimensional arrays with the common dimensions of the same length 
+#'  as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+#'  also a vector of two elements.
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  are used it has to be a vector of a pair of two logical operators: 
+#'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+#'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
+#'@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 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 to NULL and the indicator is computed using all the data provided in 
 #'  \code{data}.
-#'@param end An optional parameter to defined the final date of the period to
-#'  select from the data by providing a list of two elements: the final day of
+#'@param end An optional parameter to define the final date of the period to 
+#'  select from the data by providing a list of two elements: the final day of 
 #'  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
-#'  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 
-#'  not (FALSE). 
+#'@param time_dim A character string indicating the name of the dimension to 
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
+#'@param na.rm A logical value indicating whether to ignore NA values (TRUE) 
+#'  or not (FALSE).
 #'@param ncores An integer indicating the number of cores to use in parallel
-#'  computation.
+#'computation.
 #'
 #'@return A multidimensional array with named dimensions containing the 
-#'indicator in the element \code{data}.
+#'aggregated values with dimensions of the input parameter 'data' except the 
+#'dimension where the indicator has been computed.
 #'
-#'@import multiApply
 #'@examples
 #'# Assuming data is already (tasmax + tasmin)/2 - 10
 #'data <- array(rnorm(5 * 3 * 214 * 2, mean = 25, sd = 3),
-#'                    c(memb = 5, sdate = 3, time = 214, lon = 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"), 
@@ -146,12 +173,14 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>',
 #'               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
+#'@importFrom s2dv Reorder
 #'@export
-AccumulationExceedingThreshold <- function(data, threshold, op = '>',
-                                           diff = FALSE,
+AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FALSE,
                                            dates = NULL, start = NULL, end = NULL,
-                                           time_dim = 'time', na.rm = FALSE,
+                                           time_dim = 'ftime', na.rm = FALSE,
                                            ncores = NULL) {
+  # data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
@@ -162,82 +191,247 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>',
     dim(data) <- length(data)
     names(dim(data)) <- time_dim
   }
-  if (is.null(threshold)) {
+  if (is.null(names(dim(data)))) {
+    stop("Parameter 'data' must have named dimensions.")
+  }
+  # time_dim
+  if (!is.character(time_dim)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!all(time_dim %in% names(dim(data)))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+  if (length(time_dim) > 1) {
+    warning("Parameter 'time_dim' has length greater than 1 and ",
+            "only the first element will be used.")
+    time_dim <- time_dim[1]
+  }
+
+  # op 
+  if (!is.character(op)) {
+    stop("Parameter 'op' must be a character.")
+  }
+  if (length(op) == 1) {
+    if (!(op %in% c('>', '<', '>=', '<=', '='))) {
+      stop("Parameter 'op' must be a logical operator.")
+    }
+  } else if (length(op) == 2) {
+    op_list <- list(c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), 
+                    c('>', '<'), c('>', '<='), c('>=', '<'), c('>=', '<='))
+    if (!any(unlist(lapply(op_list, function(x) all(x == op))))) {
+      stop("Parameter 'op' is not an accepted pair of logical operators.")
+    }
+  } else {
+    stop("Parameter 'op' must be a logical operator with length 1 or 2.")
+  }
+  # threshold
+  if (is.null(unlist(threshold))) {
     stop("Parameter 'threshold' cannot be NULL.")
   }
-   if (!is.numeric(threshold)) {
+  if (!is.numeric(unlist(threshold))) {
     stop("Parameter 'threshold' must be numeric.")
   }
-  if (!is.array(threshold) && length(threshold) > 1) {
-    dim(threshold) <- length(threshold)
-    names(dim(threshold)) <- time_dim
-  } else if (length(threshold) == 1) {
-    dim(threshold) <- NULL
-  }
-  if (is.null(names(dim(data)))) {
-    stop("Parameter 'data' must have named dimensions.")
+  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."))
+    }
+    if (!is.numeric(threshold[[1]]) | !is.numeric(threshold[[2]])) {
+      stop("Parameter 'threshold' must be numeric.")
+    }
+    if (length(threshold[[1]]) != length(threshold[[2]])) {
+      stop("The pair of thresholds must have the same length.")
+    }
+
+    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.")
+      } else {
+        dim(threshold[[1]]) <- length(threshold[[1]])
+        dim(threshold[[2]]) <- length(threshold[[2]])
+        names(dim(threshold[[1]])) <- time_dim
+        names(dim(threshold[[2]])) <- time_dim
+      }
+    } else if (is.array(threshold[[1]]) && length(threshold[[1]]) > 1) {
+      if (is.null(names(dim(threshold[[1]])))) {
+        stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      if (!is.null(dim(threshold[[2]]))) {
+        if (!all(names(dim(threshold[[1]])) %in% names(dim(threshold[[2]])))) {
+          stop("The pair of thresholds must have the same dimension names.")
+        }
+      }
+      namedims <- names(dim(threshold[[1]]))
+      order <- match(namedims, names(dim(threshold[[2]])))
+      threshold[[2]] <- aperm(threshold[[2]], order)
+      if (!all(dim(threshold[[1]]) == dim(threshold[[2]]))) {
+        stop("The pair of thresholds must have the same dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold[[1]]) == 1) {
+      dim(threshold[[1]]) <- NULL
+      dim(threshold[[2]]) <- NULL
+    }
+  } 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.")
+      } else {
+        dim(threshold) <- length(threshold)
+        names(dim(threshold)) <- time_dim
+      }
+    } else if (is.array(threshold) && length(threshold) > 1) {
+      if (is.null(names(dim(threshold)))) {
+          stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold) == 1) {
+      dim(threshold) <- NULL
+    }
   }
-  if (is.null(names(dim(threshold))) && length(threshold) > 1) {
-    stop("Parameter 'threshold' must have named dimensions.")
+  # ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
   }
+  # dates
   if (!is.null(dates)) {
     if (!is.null(start) && !is.null(end)) {
       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.")
       }
-      if (all(time_dim %in% names(dim(threshold)))) {
-        if (dim(threshold)[time_dim] == dim(data)[time_dim]) {
+      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)
+          }
+        }
+      } 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)
+          }
         }
       }
       data <- SelectPeriodOnData(data, dates, start, end, 
                                  time_dim = time_dim, ncores = ncores)
     }
   }
-  if (diff == TRUE) {
+  # diff
+  if (length(op) == 2 & diff == TRUE) {
+    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.")
+    }
     data <- Apply(list(data, threshold),
                   target_dims = list(time_dim, NULL),
                   fun = function(x, y) {x - y}, ncores = ncores)$output1
     dim(data) <- dim(data)[-length(dim(data))]
     threshold <- 0
-  } 
-  if (is.null(dim(threshold))) {
-    total <- Apply(list(data), target_dims = time_dim,
-                   fun = .sumexceedthreshold,
-                   y = threshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
-  } else if (all(time_dim %in% names(dim(threshold)))) {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim, time_dim),
-                   fun = .sumexceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
-  } else if (any(time_dim %in% names(dim(threshold)))) {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim,
-                                      time_dim[time_dim %in% names(dim(threshold))]),
-                   fun = .sumexceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
+  }
+
+  ###
+
+  if (length(op) > 1) {
+    thres1 <- threshold[[1]]
+    thres2 <- threshold[[2]]
+    if (is.null(dim(thres1))) {
+      total <- Apply(list(data), target_dims = time_dim,
+                     fun = .sumexceedthreshold,
+                     y = thres1, y2 = thres2,
+                     op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(thres1)))) {
+        total <- Apply(list(data, thres1, thres2),
+                       target_dims = list(time_dim,
+                                          time_dim[time_dim %in% names(dim(thres1))],
+                                          time_dim[time_dim %in% names(dim(thres2))]),
+                       fun = .sumexceedthreshold, op = op, na.rm = na.rm,
+                       ncores = ncores)$output1
+    } else {
+        total <- Apply(list(data, thres1, thres2),
+                       target_dims = list(time_dim, thres1 = NULL, thres2 = NULL),
+                       fun = .sumexceedthreshold, op = op, na.rm = na.rm,
+                       ncores = ncores)$output1
+    }
   } else {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim, NULL),
-                   fun = .sumexceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
+    if (is.null(dim(threshold))) {
+        total <- Apply(list(data), target_dims = time_dim,
+                       fun = .sumexceedthreshold,
+                       y = threshold,
+                       op = op, na.rm = na.rm,
+                       ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(threshold)))) {
+      total <- Apply(list(data, threshold),
+                     target_dims = list(time_dim,
+                                        time_dim[time_dim %in% names(dim(threshold))]),
+                     fun = .sumexceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    } else {
+      total <- Apply(list(data, threshold),
+                     target_dims = list(time_dim, NULL),
+                     fun = .sumexceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    }
   }
   return(total)
 }
 
-.sumexceedthreshold <- function(x, y, op, na.rm) {
-  if (op == '>') {
-    res <- sum(x[x > y], na.rm = na.rm)
-  } else if (op == '<') {
-    res <- sum(x[x < y], na.rm = na.rm)
-  } else if (op == '<=') {
-    res <- sum(x[x <= y], na.rm = na.rm)
+
+.sumexceedthreshold <- function(x, y, y2 = NULL, op, na.rm) {
+  y <- as.vector(y)
+  y2 <- as.vector(y2)
+  if (is.null(y2)) {
+    if (op == '>') {
+      res <- sum(x[x > y], na.rm = na.rm)
+    } else if (op == '<') {
+      res <- sum(x[x < y], na.rm = na.rm)
+    } else if (op == '<=') {
+      res <- sum(x[x <= y], na.rm = na.rm)
+    } else if (op == '>=') {
+      res <- sum(x[x >= y], na.rm = na.rm)
+    } else {
+      res <- sum(x[x = y], na.rm = na.rm)
+    }
   } else {
-    res <- sum(x[x >= y], na.rm = na.rm)
+    if (all(op == c('<', '>'))) {
+      res <- sum(x[x < y & x > y2], na.rm = na.rm)
+    } else if (all(op == c('<', '>='))) {
+      res <- sum(x[x < y & x >= y2], na.rm = na.rm)
+    } else if (all(op == c('<=', '>'))) {
+      res <- sum(x[x <= y & x > y2], na.rm = na.rm)
+    } else if (all(op == c('<=', '>='))) {
+      res <- sum(x[x <= y & x >= y2], na.rm = na.rm)
+    } else if (all(op == c('>', '<'))) {
+      res <- sum(x[x > y & x < y2], na.rm = na.rm)
+    } else if (all(op == c('>', '<='))) {
+      res <- sum(x[x > y & x <= y2], na.rm = na.rm)
+    } else if (all(op ==  c('>=', '<'))) {
+      res <- sum(x[x >= y & x < y2], na.rm = na.rm)
+    } else if (all(op == c('>=', '<='))) {
+      res <- sum(x[x >= y & x <= y2], na.rm = na.rm)
+    }
   }
+  
   return(res)
-}
-
+}
\ No newline at end of file
diff --git a/R/TotalSpellTimeExceedingThreshold.R b/R/TotalSpellTimeExceedingThreshold.R
index ac2261ae3380e12245fbbb705ac2535fa95fcc8e..7531dee293c83db3b963ffcba639611cf7a54f62 100644
--- a/R/TotalSpellTimeExceedingThreshold.R
+++ b/R/TotalSpellTimeExceedingThreshold.R
@@ -18,14 +18,21 @@
 #'
 #'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in 
 #'  package CSTools.
-#'@param threshold An 's2dv_cube' object as output of a 'CST_' function in the 
-#'  same units as parameter 'data' and with the common dimensions of the element
-#'  'data' of the same length. A single scalar is also possible. If 
-#'  \code{timd_dim} is in the dimension (with the same length as \code{data}), 
-#'  the comparison will be done day by day. 
+#'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
+#'  object or a multidimensional array with named dimensions. It must be in the 
+#'  same units and with the common dimensions of the same length as parameter 
+#'  'data'. It can also be a vector with the same length of 'time_dim' from 
+#'  'data' or a scalar. If we want to use two thresholds: it can be a vector 
+#'  of two scalars, a list of two vectors with the same length of 
+#'  'time_dim' from 'data' or a list of two multidimensional arrays with the 
+#'  common dimensions of the same length as parameter 'data'. If two thresholds
+#'  are used, parameter 'op' must be also a vector of two elements.
 #'@param spell A scalar indicating the minimum length of the spell.
-#'@param op An operator '>' (by default), '<', '>=' or '<='. 
-#'@param start An optional parameter to defined the initial date of the period 
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  are used it has to be a vector of a pair of two logical operators: 
+#'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+#'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
+#'@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
 #'  to NULL and the indicator is computed using all the data provided in 
@@ -35,14 +42,13 @@
 #'  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 
-#'  dimension name matching the dimensions provided in the object 
-#'  \code{data$data} can be specified.
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
 #'@param ncores An integer indicating the number of cores to use in parallel 
 #'  computation.
 #'
-#'@return An 's2dv_cube' object containing the indicator in the element 
-#'\code{data}.
+#'@return An 's2dv_cube' object containing the number of days that are part of a
+#'spell within a threshold in element \code{data}.
 #'
 #'@examples
 #'exp <- NULL
@@ -82,18 +88,29 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
       }
     }
   }
-  if (inherits(threshold, 's2dv_cube')) {
-    threshold <- threshold$data
+  if (length(op) == 1) {
+    if (inherits(threshold, 's2dv_cube')) {
+        threshold <- threshold$data
+    }
+  } else if (length(op) == 2) {
+    if (inherits(threshold[[1]], 's2dv_cube')) {
+      threshold[[1]] <- threshold[[1]]$data
+    }
+    if (inherits(threshold[[2]], 's2dv_cube')) {
+      threshold[[2]] <- threshold[[2]]$data
+    }
   }
+
+
   total <- TotalSpellTimeExceedingThreshold(data$data, data$Dates[[1]],
                                             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$Dates <- SelectPeriodOnDates(dates = data$Dates$start,
-                                       start = start, end = end,
-                                       time_dim = time_dim, ncores = ncores)
+    data$Dates <- SelectPeriodOnDates(dates = data$Dates$start,
+                                      start = start, end = end,
+                                      time_dim = time_dim, ncores = ncores)
   }
   return(data)
 }
@@ -116,12 +133,20 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
 #'@seealso [Threshold()] and [AbsToProbs()].
 #'
 #'@param data A multidimensional array with named dimensions.
-#'@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. If \code{timd_dim} is in the dimension (with the 
-#'  same length as \code{data}), the comparison will be done day by day.
+#'@param threshold If only one threshold is used: it can be a multidimensional 
+#'  array with named dimensions. It must be in the same units and with the 
+#'  common dimensions of the same length as parameter 'data'. It can also be a
+#'  vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+#'  want to use two thresholds: it can be a vector of two scalars, a list of 
+#'  two vectors with the same length of 'time_dim' from 'data' or a list of 
+#'  two multidimensional arrays with the common dimensions of the same length 
+#'  as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+#'  also a vector of two elements.
 #'@param spell A scalar indicating the minimum length of the spell.
-#'@param op An operator '>' (by default), '<', '>=' or '<='. 
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  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.
@@ -130,34 +155,36 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
 #'  date of the period and the initial 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 end An optional parameter to defined the final date of the period to 
+#'@param end An optional parameter to define the final date of the period to 
 #'  select from the data by providing a list of two elements: the final day of 
 #'  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 
-#'  dimension name matching the dimensions provided in the object 
-#'  \code{data$data} can be specified.
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
 #'@param ncores An integer indicating the number of cores to use in parallel 
 #'  computation.
 #'
-#'@return A multidimensional array with named dimensions containing the indicator
-#'in the element \code{data}.
+#'@return A multidimensional array with named dimensions containing the number 
+#'of days that are part of a spell within a threshold with dimensions of the 
+#'input parameter 'data' except the dimension where the indicator has been 
+#'computed.
 #'
 #'@details This function considers NA values as the end of the spell. For a 
 #'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, time = 20, lat = 4))
+#'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)
 #' 
 #'@import multiApply
 #'@export
-TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>', dates = NULL,
-                                             start = NULL, end = NULL, time_dim = 'time',
-                                             ncores = NULL) {
+TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>',
+                                             dates = NULL, start = NULL, end = NULL, 
+                                             time_dim = 'ftime', ncores = NULL) {
+  # data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
@@ -168,62 +195,232 @@ TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>', d
     dim(data) <- length(data)
     names(dim(data)) <- time_dim
   }
-  if (is.null(threshold)) {
+  if (is.null(names(dim(data)))) {
+    stop("Parameter 'data' must have named dimensions.")
+  }
+  # time_dim
+  if (!is.character(time_dim)) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!all(time_dim %in% names(dim(data)))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+  if (length(time_dim) > 1) {
+    warning("Parameter 'time_dim' has length greater than 1 and ",
+            "only the first element will be used.")
+    time_dim <- time_dim[1]
+  }
+  # op
+  if (!is.character(op)) {
+    stop("Parameter 'op' must be a character.")
+  }
+  if (length(op) == 1) {
+    if (!(op %in% c('>', '<', '>=', '<=', '='))) {
+      stop("Parameter 'op' must be a logical operator.")
+    }
+  } else if (length(op) == 2) {
+    op_list <- list(c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), 
+                    c('>', '<'), c('>', '<='), c('>=', '<'), c('>=', '<='))
+    if (!any(unlist(lapply(op_list, function(x) all(x == op))))) {
+      stop("Parameter 'op' is not an accepted pair of logical operators.")
+    }
+  } else {
+    stop("Parameter 'op' must be a logical operator with length 1 or 2.")
+  }
+  # threshold
+  if (is.null(unlist(threshold))) {
     stop("Parameter 'threshold' cannot be NULL.")
   }
-  if (!is.numeric(threshold)) {
+  if (!is.numeric(unlist(threshold))) {
     stop("Parameter 'threshold' must be numeric.")
   }
-  if (is.null(names(dim(data)))) {
-    stop("Parameter 'data' must have named dimensions.")
+  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."))
+    }
+    if (!is.numeric(threshold[[1]]) | !is.numeric(threshold[[2]])) {
+      stop("Parameter 'threshold' must be numeric.")
+    }
+    if (length(threshold[[1]]) != length(threshold[[2]])) {
+      stop("The pair of thresholds must have the same length.")
+    }
+    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.")
+      } else {
+        dim(threshold[[1]]) <- length(threshold[[1]])
+        dim(threshold[[2]]) <- length(threshold[[2]])
+        names(dim(threshold[[1]])) <- time_dim
+        names(dim(threshold[[2]])) <- time_dim
+      }
+    } else if (is.array(threshold[[1]]) && length(threshold[[1]]) > 1) {
+      if (is.null(names(dim(threshold[[1]])))) {
+        stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      if (!is.null(dim(threshold[[2]]))) {
+        if (!all(names(dim(threshold[[1]])) %in% names(dim(threshold[[2]])))) {
+          stop("The pair of thresholds must have the same dimension names.")
+        }
+      }
+      namedims <- names(dim(threshold[[1]]))
+      order <- match(namedims, names(dim(threshold[[2]])))
+      threshold[[2]] <- aperm(threshold[[2]], order)
+      if (!all(dim(threshold[[1]]) == dim(threshold[[2]]))) {
+        stop("The pair of thresholds must have the same dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold[[1]]) == 1) {
+      dim(threshold[[1]]) <- NULL
+      dim(threshold[[2]]) <- NULL
+    }
+  } 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.")
+      } else {
+        dim(threshold) <- length(threshold)
+        names(dim(threshold)) <- time_dim
+      }
+    } else if (is.array(threshold) && length(threshold) > 1) {
+      if (is.null(names(dim(threshold)))) {
+          stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold) == 1) {
+      dim(threshold) <- NULL
+    }
   }
+  # spell
+  if (!is.numeric(spell) | length(spell) != 1) {
+    stop("Parameter 'spell' must be a scalar.")
+  }
+  # ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  # dates
   if (!is.null(dates)) {
     if (!is.null(start) && !is.null(end)) {
       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.")
       }
-      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)
+      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)
+          }
+        }
+      } 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)
+          }
         }
       }
       data <- SelectPeriodOnData(data, dates, start, end, 
                                  time_dim = time_dim, ncores = ncores)
     }
   }
-  if (is.null(dim(threshold))) {
-    total <- Apply(list(data), target_dims = time_dim,
-                   fun = .totalspellthres,
-                   threshold = threshold, spell = spell, op = op,
-                   ncores = ncores)$output1
-  } else if (any(time_dim %in% names(dim(threshold)))) {
-    total <- Apply(list(data, threshold),
-                 target_dims = list(time_dim, 
-                                    time_dim[time_dim %in% names(dim(threshold))]),
-                 fun = .totalspellthres, spell = spell, op = op,
-                 ncores = ncores)$output1
 
+  if (length(op) > 1) {
+    thres1 <- threshold[[1]]
+    thres2 <- threshold[[2]]
+    if (is.null(dim(thres1))) {
+      total <- Apply(list(data), target_dims = time_dim,
+                    fun = .totalspellthres, y = thres1, y2 = thres2,
+                    spell = spell, op = op,
+                    ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(thres1)))) {
+      total <- Apply(list(data, thres1, thres2),
+                     target_dims = list(time_dim, 
+                                        time_dim[time_dim %in% names(dim(thres1))],
+                                        time_dim[time_dim %in% names(dim(thres2))]),
+                     fun = .totalspellthres, spell = spell, op = op,
+                     ncores = ncores)$output1
+
+    } else {
+      total <- Apply(list(data, thres1, thres2), 
+                     target_dims = list(time_dim, thres1 = NULL, thres2 = NULL), 
+                     fun = .totalspellthres, spell = spell, op = op,
+                     ncores = ncores)$output1
+    }
   } else {
-    total <- Apply(list(data, threshold), target_dims = list(time_dim, NULL), 
-                 fun = .totalspellthres, spell = spell, op = op,
-                 ncores = ncores)$output1
+    if (is.null(dim(threshold))) {
+      total <- Apply(list(data), target_dims = time_dim,
+                    fun = .totalspellthres,
+                    y = threshold, spell = spell, op = op,
+                    ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(threshold)))) {
+      total <- Apply(list(data, threshold),
+                     target_dims = list(time_dim, 
+                                        time_dim[time_dim %in% names(dim(threshold))]),
+                     fun = .totalspellthres, spell = spell, op = op,
+                     ncores = ncores)$output1
+
+    } else {
+      total <- Apply(list(data, threshold), 
+                     target_dims = list(time_dim, NULL), 
+                     fun = .totalspellthres, spell = spell, op = op,
+                     ncores = ncores)$output1
+    }
   }
   return(total) 
 }  
 
-.totalspellthres <- function(data, threshold, spell, op = '>') {
-  # data a time serie, threshold single value:
-  if (op == '>') {
-    exceed <- data > threshold
-  } else if (op == '<') {
-    exceed <- data < threshold
-  } else if (op == '<=') {
-    exceed <- data <= threshold
+.totalspellthres <- function(x, y, y2 = NULL, spell, op = '>') {
+  y <- as.vector(y)
+  y2 <- as.vector(y2)
+  if (is.null(y2)) {
+    if (op == '>') {
+      exceed <- x > y
+    } else if (op == '<') {
+      exceed <- x < y
+    } else if (op == '<=') {
+      exceed <- x <= y
+    } else {
+      exceed <- x >= y
+    }
   } else {
-    exceed <- data >= threshold
+    if (all(op == c('<', '>'))) {
+      exceed <- x < y & x > y2
+    } else if (all(op == c('<', '>='))) {
+      exceed <- x < y & x >= y2
+    } else if (all(op == c('<=', '>'))) {
+      exceed <- x <= y & x > y2
+    } else if (all(op == c('<=', '>='))) {
+      exceed <- x <= y & x >= y2
+    } else if (all(op == c('>', '<'))) {
+      exceed <- x > y & x < y2
+    } else if (all(op == c('>', '<='))) {
+      exceed <- x > y & x <= y2
+    } else if (all(op ==  c('>=', '<'))) {
+      exceed <- x >= y & x < y2
+    } else if (all(op == c('>=', '<='))) {
+      exceed <- x >= y & x <= y2
+    }
   }
+
   spells_exceed <- sequence(rle(as.character(exceed))$lengths)
   spells_exceed[exceed == FALSE] <- NA
   pos_spells <- which(spells_exceed == spell)
@@ -237,4 +434,4 @@ TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>', d
                         return(days)
                         })))
   return(total)
-}
+}
\ No newline at end of file
diff --git a/R/TotalTimeExceedingThreshold.R b/R/TotalTimeExceedingThreshold.R
index ec25244d879adf3cd6d561e3b4aff5a6bc31bfbd..01d7823faf0d5848f41cb6f1b8ff30e5e5b58146 100644
--- a/R/TotalTimeExceedingThreshold.R
+++ b/R/TotalTimeExceedingThreshold.R
@@ -1,12 +1,12 @@
 #'Total Time of a variable Exceeding (not exceeding) a Threshold
 #'
-#'The Total Time of a variable exceeding (or not) a Threshold returns the total 
-#'number of days (if the data provided is daily, or the corresponding units to 
-#'the data frequency provided) that a variable is exceeding a threshold during a 
-#'period. The threshold provided must be in the same units than the variable 
-#'units, i.e. to use a percentile as a scalar, 
-#'the function \code{AbsToProbs} or \code{QThreshold}  may be needed (see 
-#'examples). Providing maximum temperature daily data, the following agriculture 
+#'The Total Time of a variable exceeding (or not) a Threshold. It returns the  
+#'total number of days (if the data provided is daily, or the corresponding  
+#'units of the data frequency) that a variable is exceeding a threshold 
+#'during a period. The threshold provided must be in the same units as the 
+#'variable units, i.e. to use a percentile as a scalar, the function 
+#'\code{AbsToProbs} or \code{QThreshold}  may be needed (see examples). 
+#'Providing maximum temperature daily data, the following agriculture 
 #'indices for heat stress can be obtained by using this function:
 #'\itemize{
 #'  \item\code{SU35}{Total count of days when daily maximum temperatures exceed 
@@ -22,32 +22,39 @@
 #'
 #'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in 
 #'  package CSTools.
-#'@param threshold An 's2dv_cube' object as output of a 'CST_' function in the 
-#'  same units as parameter \code{data} and with the common dimensions of the 
-#'  element \code{data} of the same length (e.g. an array with the same lengths 
-#'  of longitude and latitude). A single scalar is also possible (for the case 
-#'  of comparing all grid points with the same scalar).
-#'@param op An operator '>' (by default), '<', '>=' or '<='. 
-#'@param start An optional parameter to defined the initial date of the period 
+#'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
+#'  object or a multidimensional array with named dimensions. It must be in the 
+#'  same units and with the common dimensions of the same length as parameter 
+#'  'data'. It can also be a vector with the same length of 'time_dim' from 
+#'  'data' or a scalar. If we want to use two thresholds: it can be a vector 
+#'  of two scalars, a list of two vectors with the same length of 
+#'  'time_dim' from 'data' or a list of two multidimensional arrays with the 
+#'  common dimensions of the same length as parameter 'data'. If two thresholds
+#'  are used, parameter 'op' must be also a vector of two elements.
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  are used it has to be a vector of a pair of two logical operators: 
+#'  c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+#'  c('>', '<='), c('>=', '<'),c('>=', '<=')).
+#'@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
 #'  to NULL and the indicator is computed using all the data provided in 
 #'  \code{data}.
-#'@param end An optional parameter to defined the final date of the period to 
+#'@param end An optional parameter to define the final date of the period to 
 #'  select from the data by providing a list of two elements: the final day of 
 #'  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 
-#'  dimension name matching the dimensions provided in the object 
-#'  \code{data$data} can be specified.
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
 #'  not (FALSE). 
 #'@param ncores An integer indicating the number of cores to use in parallel 
 #'  computation.
 #'
-#'@return An 's2dv_cube' object containing the indicator in the element 
-#'\code{data}.
+#'@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.
 #' 
 #'@examples
 #'exp <- NULL
@@ -87,8 +94,17 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
       }
     }
   }
-  if (inherits(threshold, 's2dv_cube')) {
-    threshold <- threshold$data
+  if (length(op) == 1) {
+    if (inherits(threshold, 's2dv_cube')) {
+        threshold <- threshold$data
+    }
+  } else if (length(op) == 2) {
+    if (inherits(threshold[[1]], 's2dv_cube')) {
+      threshold[[1]] <- threshold[[1]]$data
+    }
+    if (inherits(threshold[[2]], 's2dv_cube')) {
+      threshold[[2]] <- threshold[[2]]$data
+    }
   }
   total <- TotalTimeExceedingThreshold(data$data, data$Dates[[1]],
                                        threshold = threshold, op = op,
@@ -104,17 +120,18 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
 }
 #'Total Time of a variable Exceeding (not exceeding) a Threshold
 #'
-#'The Total Time of a variable exceeding (or not) a Threshold returns the total 
-#'number of days (if the data provided is daily, or the corresponding units to 
-#'the data frequency provided) that a variable is exceeding a threshold during a 
-#'period. The threshold provided must be in the same units than the variable 
-#'units, i.e. to use a percentile as a threshold, the function \code{Threshold} 
-#'or \code{QThreshold} may be needed (see examples). Providing maximum 
-#'temperature daily data, the following agriculture indices for heat stress can 
-#'be obtained by using this function:
+#'The Total Time of a variable exceeding (or not) a Threshold. It returns the  
+#'total number of days (if the data provided is daily, or the corresponding  
+#'units of the data frequency) that a variable is exceeding a threshold 
+#'during a period. The threshold provided must be in the same units as the 
+#'variable units, i.e. to use a percentile as a scalar, the function 
+#'\code{AbsToProbs} or \code{QThreshold}  may be needed (see examples). 
+#'Providing maximum temperature daily data, the following agriculture 
+#'indices for heat stress can be obtained by using this function:
 #'\itemize{
 #'  \item\code{SU35}{Total count of days when daily maximum temperatures exceed 
-#'                   35°C}
+#'                   35°C in the seven months from the start month given (e.g. 
+#'                   from April to October for start month of April).}
 #'  \item\code{SU36}{Total count of days when daily maximum temperatures exceed 
 #'                   36 between June 21st and September 21st}
 #'  \item\code{SU40}{Total count of days when daily maximum temperatures exceed 
@@ -124,35 +141,42 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
 #'}
 #'
 #'@param data A multidimensional array with named dimensions.
-#'@param threshold A multidimensional array with named dimensions in the same 
-#'  units as parameter \code{data} and with the common dimensions of the element 
-#'  \code{data} of the same length (e.g. an array with the same lengths of 
-#'  longitude and latitude). A single scalar is also possible (for the case of 
-#'  comparing all grid points with the same scalar).
-#'@param op A operator '>' (by default), '<', '>=' or '<='. 
+#'@param threshold If only one threshold is used: it can be a multidimensional 
+#'  array with named dimensions. It must be in the same units and with the 
+#'  common dimensions of the same length as parameter 'data'. It can also be a
+#'  vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+#'  want to use two thresholds: it can be a vector of two scalars, a list of 
+#'  two vectors with the same length of 'time_dim' from 'data' or a list of 
+#'  two multidimensional arrays with the common dimensions of the same length 
+#'  as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+#'  also a vector of two elements.
+#'@param op An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+#'  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 start An optional parameter to defined the initial date of the period 
+#'@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
 #'  to NULL and the indicator is computed using all the data provided in 
 #'  \code{data}.
-#'@param end An optional parameter to defined the final date of the period to 
+#'@param end An optional parameter to define the final date of the period to 
 #'  select from the data by providing a list of two elements: the final day of 
 #'  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 'time'. More than one 
-#'  dimension name matching the dimensions provided in the object 
-#'  \code{data$data} can be specified.
+#'  compute the indicator. By default, it is set to 'ftime'. It can only
+#'  indicate one time dimension.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
 #'  not (FALSE). 
 #'@param ncores An integer indicating the number of cores to use in parallel 
 #'  computation.
 #'
-#'@return A multidimensional array with named dimensions containing the 
-#'indicator in the element \code{data}.
+#'@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.
 #'
 #'@examples
 #'exp <- array(abs(rnorm(5 * 3 * 214 * 2)*280),
@@ -163,8 +187,9 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
 #'@export
 TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
                                         dates = NULL, start = NULL, end = NULL,
-                                        time_dim = 'time', na.rm = FALSE,
+                                        time_dim = 'ftime', na.rm = FALSE,
                                         ncores = NULL) {
+  # data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
@@ -175,81 +200,225 @@ TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
     dim(data) <- length(data)
     names(dim(data)) <- time_dim
   }
-  if (is.null(threshold)) {
-        stop("Parameter 'threshold' cannot be NULL.")
+  if (is.null(names(dim(data)))) {
+    stop("Parameter 'data' must have named dimensions.")
   }
-   if (!is.numeric(threshold)) {
-    stop("Parameter 'threshold' must be numeric.")
+  # time_dim 
+  if (!is.character(time_dim)) {
+    stop("Parameter 'time_dim' must be a character string.")
   }
-  if (!is.array(threshold) && length(threshold) > 1) {
-    dim(threshold) <- length(threshold)
-    names(dim(threshold)) <- time_dim
-  } else if (length(threshold) == 1) {
-    dim(threshold) <- NULL
+  if (!all(time_dim %in% names(dim(data)))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
   }
-  if (is.null(names(dim(data)))) {
-    stop("Parameter 'data' must have named dimensions.")
+  if (length(time_dim) > 1) {
+    warning("Parameter 'time_dim' has length greater than 1 and ",
+            "only the first element will be used.")
+    time_dim <- time_dim[1]
+  }
+  # op
+  if (!is.character(op)) {
+    stop("Parameter 'op' must be a character.")
   }
-  if (is.null(names(dim(threshold))) && length(threshold) > 1) {
-    stop("Parameter 'threshold' must have named dimensions.")
+  if (length(op) == 1) {
+    if (!(op %in% c('>', '<', '>=', '<=', '='))) {
+      stop("Parameter 'op' must be a logical operator.")
+    }
+  } else if (length(op) == 2) {
+    op_list <- list(c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), 
+                    c('>', '<'), c('>', '<='), c('>=', '<'), c('>=', '<='))
+    if (!any(unlist(lapply(op_list, function(x) all(x == op))))) {
+      stop("Parameter 'op' is not an accepted pair of logical operators.")
+    }
+  } else {
+    stop("Parameter 'op' must be a logical operator with length 1 or 2.")
+  }
+  # threshold
+  if (is.null(unlist(threshold))) {
+    stop("Parameter 'threshold' cannot be NULL.")
+  }
+  if (!is.numeric(unlist(threshold))) {
+    stop("Parameter 'threshold' must be numeric.")
   }
-  common_dims <- which(names(dim(data)) %in% names(dim(threshold)))
-  if (length(threshold) > 1) {    
-    if (any(dim(data)[common_dims] != 
-            dim(threshold)[which(names(dim(threshold)) %in% names(dim(data)))])) {
-    stop("Parameter 'data' and 'threshold' must have the same length on common dimensions.")
+  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."))
     }
-  } 
+    if (!is.numeric(threshold[[1]]) | !is.numeric(threshold[[2]])) {
+      stop("Parameter 'threshold' must be numeric.")
+    }
+    if (length(threshold[[1]]) != length(threshold[[2]])) {
+      stop("The pair of thresholds must have the same length.")
+    }
+    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.")
+      } else {
+        dim(threshold[[1]]) <- length(threshold[[1]])
+        dim(threshold[[2]]) <- length(threshold[[2]])
+        names(dim(threshold[[1]])) <- time_dim
+        names(dim(threshold[[2]])) <- time_dim
+      }
+    } else if (is.array(threshold[[1]]) && length(threshold[[1]]) > 1) {
+      if (is.null(names(dim(threshold[[1]])))) {
+        stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      if (!is.null(dim(threshold[[2]]))) {
+        if (!all(names(dim(threshold[[1]])) %in% names(dim(threshold[[2]])))) {
+          stop("The pair of thresholds must have the same dimension names.")
+        }
+      }
+      namedims <- names(dim(threshold[[1]]))
+      order <- match(namedims, names(dim(threshold[[2]])))
+      threshold[[2]] <- aperm(threshold[[2]], order)
+      if (!all(dim(threshold[[1]]) == dim(threshold[[2]]))) {
+        stop("The pair of thresholds must have the same dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold[[1]]) == 1) {
+      dim(threshold[[1]]) <- NULL
+      dim(threshold[[2]]) <- NULL
+    }
+  } 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.")
+      } else {
+        dim(threshold) <- length(threshold)
+        names(dim(threshold)) <- time_dim
+      }
+    } else if (is.array(threshold) && length(threshold) > 1) {
+      if (is.null(names(dim(threshold)))) {
+          stop("If parameter 'threshold' is an array it must have named dimensions.")
+      }
+      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."))
+        }
+      }
+    } else if (length(threshold) == 1) {
+      dim(threshold) <- NULL
+    }
+  }
+  # ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  # dates
   if (!is.null(dates)) {
     if (!is.null(start) && !is.null(end)) {
       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.")
       }
-      if (time_dim %in% names(dim(threshold))) {
-        if (dim(threshold)[time_dim] == dim(data)[time_dim]) {
+      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)
+                                            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)
+          }
         }
       }
       data <- SelectPeriodOnData(data, dates, start, end, 
                                  time_dim = time_dim, ncores = ncores)
     }
   }
-  if (is.null(dim(threshold))) {
-    total <- Apply(list(data), target_dims = time_dim,
-                   fun = .exceedthreshold,
-                   y = threshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
-  } else if (all(time_dim %in% names(dim(threshold)))) {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim, time_dim),
-                   fun = .exceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
-  } else if (any(time_dim %in% names(dim(threshold)))) {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim, 
-                                      time_dim[time_dim %in% names(dim(threshold))]),
-                   fun = .exceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1   
+
+  if (length(op) > 1) {
+    thres1 <- threshold[[1]]
+    thres2 <- threshold[[2]]
+    if (is.null(dim(thres1))) {
+      total <- Apply(list(data), target_dims = time_dim,
+                    fun = .exceedthreshold, y = thres1, y2 = thres2,
+                    op = op, na.rm = na.rm,
+                    ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(thres1)))) {
+      total <- Apply(list(data, thres1, thres2),
+                     target_dims = list(time_dim, 
+                                        time_dim[time_dim %in% names(dim(thres1))],
+                                        time_dim[time_dim %in% names(dim(thres2))]),
+                     fun = .exceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+
+    } else {
+      total <- Apply(list(data, thres1, thres2), 
+                     target_dims = list(time_dim, thres1 = NULL, thres2 = NULL), 
+                     fun = .exceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    }
   } else {
-    total <- Apply(list(data, threshold),
-                   target_dims = list(time_dim, NULL),
-                   fun = .exceedthreshold, op = op, na.rm = na.rm,
-                   ncores = ncores)$output1
+    if (is.null(dim(threshold))) {
+      total <- Apply(list(data), target_dims = time_dim,
+                     fun = .exceedthreshold,
+                     y = threshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    } else if (any(time_dim %in% names(dim(threshold)))) {
+      total <- Apply(list(data, threshold),
+                     target_dims = list(time_dim, 
+                                        time_dim[time_dim %in% names(dim(threshold))]),
+                     fun = .exceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1   
+    } else {
+      total <- Apply(list(data, threshold),
+                     target_dims = list(time_dim, NULL),
+                     fun = .exceedthreshold, op = op, na.rm = na.rm,
+                     ncores = ncores)$output1
+    }
   }
   return(total)
 }
-.exceedthreshold <- function(x, y, op, na.rm) {
-  if (op == '>') {
-    res <- sum(x > y, na.rm = na.rm)
-  } else if (op == '<') {
-    res <- sum(x < y, na.rm = na.rm)
-  } else if (op == '<=') {
-    res <- sum(x <= y, na.rm = na.rm)
+
+.exceedthreshold <- function(x, y, y2 = NULL, op = '>', na.rm) {
+  y <- as.vector(y)
+  y2 <- as.vector(y2)
+  if (is.null(y2)) {
+    if (op == '>') {
+      res <- sum(x > y, na.rm = na.rm)
+    } else if (op == '<') {
+      res <- sum(x < y, na.rm = na.rm)
+    } else if (op == '<=') {
+      res <- sum(x <= y, na.rm = na.rm)
+    } else {
+      res <- sum(x >= y, na.rm = na.rm)
+    }
   } else {
-    res <- sum(x >= y, na.rm = na.rm)
+    if (all(op == c('<', '>'))) {
+      res <- sum(x < y & x > y2, na.rm = na.rm)
+    } else if (all(op == c('<', '>='))) {
+      res <- sum(x < y & x >= y2, na.rm = na.rm)
+    } else if (all(op == c('<=', '>'))) {
+      res <- sum(x <= y & x > y2, na.rm = na.rm)
+    } else if (all(op == c('<=', '>='))) {
+      res <- sum(x <= y & x >= y2, na.rm = na.rm)
+    } else if (all(op == c('>', '<'))) {
+      res <- sum(x > y & x < y2, na.rm = na.rm)
+    } else if (all(op == c('>', '<='))) {
+      res <- sum(x > y & x <= y2, na.rm = na.rm)
+    } else if (all(op ==  c('>=', '<'))) {
+      res <- sum(x >= y & x < y2, na.rm = na.rm)
+    } else if (all(op == c('>=', '<='))) {
+      res <- sum(x >= y & x <= y2, na.rm = na.rm)
+    }
   }
   return(res)
-}
-
+}
\ No newline at end of file
diff --git a/man/AccumulationExceedingThreshold.Rd b/man/AccumulationExceedingThreshold.Rd
index e646dc627b325dece44c37eb14f72842f54998a8..172592c89c091e771f5a13f2084f594dc25b51e2 100644
--- a/man/AccumulationExceedingThreshold.Rd
+++ b/man/AccumulationExceedingThreshold.Rd
@@ -12,7 +12,7 @@ AccumulationExceedingThreshold(
   dates = NULL,
   start = NULL,
   end = NULL,
-  time_dim = "time",
+  time_dim = "ftime",
   na.rm = FALSE,
   ncores = NULL
 )
@@ -20,44 +20,54 @@ AccumulationExceedingThreshold(
 \arguments{
 \item{data}{A multidimensional array with named dimensions.}
 
-\item{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.}
+\item{threshold}{If only one threshold is used: it can be a multidimensional 
+array with named dimensions. It must be in the same units and with the 
+common dimensions of the same length as parameter 'data'. It can also be a
+vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+want to use two thresholds: it can be a vector of two scalars, a list of 
+two vectors with the same length of 'time_dim' from 'data' or a list of 
+two multidimensional arrays with the common dimensions of the same length 
+as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+also a vector of two elements.}
 
-\item{op}{An operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+are used it has to be a vector of a pair of two logical operators: 
+c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{diff}{A logical value indicating whether to accumulate the difference
-between data and threshold (TRUE) or not (FALSE by default).}
+\item{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.}
 
-\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 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{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
-to NULL and the indicator is computed using all the data provided in
+\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 
+date of the period and the initial 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{end}{An optional parameter to defined the final date of the period to
-select from the data by providing a list of two elements: the final day of
+\item{end}{An optional parameter to define the final date of the period to 
+select from the data by providing a list of two elements: the final day of 
 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
-dimension name matching the dimensions provided in the object
-\code{data$data} can be specified.}
+\item{time_dim}{A character string indicating the name of the dimension to 
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
-\item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
-not (FALSE).}
+\item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) 
+or not (FALSE).}
 
 \item{ncores}{An integer indicating the number of cores to use in parallel
 computation.}
 }
 \value{
 A multidimensional array with named dimensions containing the 
-indicator in the element \code{data}.
+aggregated values with dimensions of the input parameter 'data' except the 
+dimension where the indicator has been computed.
 }
 \description{
 The accumulation (sum) of a variable in the days (or time steps) that the 
@@ -75,7 +85,7 @@ function:
 \examples{
 # Assuming data is already (tasmax + tasmin)/2 - 10
 data <- array(rnorm(5 * 3 * 214 * 2, mean = 25, sd = 3),
-                   c(memb = 5, sdate = 3, time = 214, lon = 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"), 
diff --git a/man/CST_AccumulationExceedingThreshold.Rd b/man/CST_AccumulationExceedingThreshold.Rd
index 9c0a521ae8bda6a57e102060bb91804bdaf496e6..bc0eb83f4d4d1b168b956c63facbbb6ea5de4fc4 100644
--- a/man/CST_AccumulationExceedingThreshold.Rd
+++ b/man/CST_AccumulationExceedingThreshold.Rd
@@ -17,42 +17,53 @@ CST_AccumulationExceedingThreshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in 
+\item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in
 package CSTools.}
 
-\item{threshold}{An 's2dv_cube' object as output of a 'CST_' function in the
-same units as parameter 'data' and with the common dimensions of the element 
-'data' of the same length. A single scalar is also possible.}
+\item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
+object or a multidimensional array with named dimensions. It must be in the 
+same units and with the common dimensions of the same length as parameter 
+'data'. It can also be a vector with the same length of 'time_dim' from 
+'data' or a scalar. If we want to use two thresholds: it can be a vector 
+of two scalars, a list of two vectors with the same length of 
+'time_dim' from 'data' or a list of two multidimensional arrays with the 
+common dimensions of the same length as parameter 'data'. If two thresholds
+are used, parameter 'op' must be also a vector of two elements.}
 
-\item{op}{An operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+are used it has to be a vector of a pair of two logical operators: 
+c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{diff}{A logical value indicating whether to accumulate the difference
-between data and threshold (TRUE) or not (FALSE by default).}
+\item{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.}
 
-\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
-date of the period and the initial month of the period. By default it is set
-to NULL and the indicator is computed using all the data provided in
+\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 
+date of the period and the initial 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{end}{An optional parameter to defined the final date of the period to
-select from the data by providing a list of two elements: the final day of
+\item{end}{An optional parameter to define the final date of the period to 
+select from the data by providing a list of two elements: the final day of 
 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
-dimension name matching the dimensions provided in the object
-\code{data$data} can be specified.}
+\item{time_dim}{A character string indicating the name of the dimension to 
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
-\item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or
-not (FALSE).}
+\item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) 
+or not (FALSE).}
 
 \item{ncores}{An integer indicating the number of cores to use in parallel
 computation.}
 }
 \value{
-A 's2dv_cube' object containing the indicator in the element \code{data}.
+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.
 }
 \description{
 The accumulation (sum) of a variable in the days (or time steps) that the 
diff --git a/man/CST_TotalSpellTimeExceedingThreshold.Rd b/man/CST_TotalSpellTimeExceedingThreshold.Rd
index 847fed25a2a4a66c22b918f3e7c23b2ea84833b6..75a2d1e2c55dceffc2e3b708e562b35d7e219435 100644
--- a/man/CST_TotalSpellTimeExceedingThreshold.Rd
+++ b/man/CST_TotalSpellTimeExceedingThreshold.Rd
@@ -19,17 +19,24 @@ CST_TotalSpellTimeExceedingThreshold(
 \item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in 
 package CSTools.}
 
-\item{threshold}{An 's2dv_cube' object as output of a 'CST_' function in the 
-same units as parameter 'data' and with the common dimensions of the element
-'data' of the same length. A single scalar is also possible. If 
-\code{timd_dim} is in the dimension (with the same length as \code{data}), 
-the comparison will be done day by day.}
+\item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
+object or a multidimensional array with named dimensions. It must be in the 
+same units and with the common dimensions of the same length as parameter 
+'data'. It can also be a vector with the same length of 'time_dim' from 
+'data' or a scalar. If we want to use two thresholds: it can be a vector 
+of two scalars, a list of two vectors with the same length of 
+'time_dim' from 'data' or a list of two multidimensional arrays with the 
+common dimensions of the same length as parameter 'data'. If two thresholds
+are used, parameter 'op' must be also a vector of two elements.}
 
 \item{spell}{A scalar indicating the minimum length of the spell.}
 
-\item{op}{An operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+are used it has to be a vector of a pair of two logical operators: 
+c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{start}{An optional parameter to defined the initial date of the period 
+\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 
 date of the period and the initial month of the period. By default it is set
 to NULL and the indicator is computed using all the data provided in 
@@ -41,16 +48,15 @@ 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 
-dimension name matching the dimensions provided in the object 
-\code{data$data} can be specified.}
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
 \item{ncores}{An integer indicating the number of cores to use in parallel 
 computation.}
 }
 \value{
-An 's2dv_cube' object containing the indicator in the element 
-\code{data}.
+An 's2dv_cube' object containing the number of days that are part of a
+spell within a threshold in element \code{data}.
 }
 \description{
 The number of days (when daily data is provided) that are part of a spell 
diff --git a/man/CST_TotalTimeExceedingThreshold.Rd b/man/CST_TotalTimeExceedingThreshold.Rd
index bbd05e04b7091c524e7468a754644d5c21ba3e1b..5dea9647c56d673e353f7357bfb1018a3d7656ef 100644
--- a/man/CST_TotalTimeExceedingThreshold.Rd
+++ b/man/CST_TotalTimeExceedingThreshold.Rd
@@ -19,29 +19,35 @@ CST_TotalTimeExceedingThreshold(
 \item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in 
 package CSTools.}
 
-\item{threshold}{An 's2dv_cube' object as output of a 'CST_' function in the 
-same units as parameter \code{data} and with the common dimensions of the 
-element \code{data} of the same length (e.g. an array with the same lengths 
-of longitude and latitude). A single scalar is also possible (for the case 
-of comparing all grid points with the same scalar).}
+\item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
+object or a multidimensional array with named dimensions. It must be in the 
+same units and with the common dimensions of the same length as parameter 
+'data'. It can also be a vector with the same length of 'time_dim' from 
+'data' or a scalar. If we want to use two thresholds: it can be a vector 
+of two scalars, a list of two vectors with the same length of 
+'time_dim' from 'data' or a list of two multidimensional arrays with the 
+common dimensions of the same length as parameter 'data'. If two thresholds
+are used, parameter 'op' must be also a vector of two elements.}
 
-\item{op}{An operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+are used it has to be a vector of a pair of two logical operators: 
+c('<', '>'), c('<', '>='), c('<=', '>'), c('<=', '>='), c('>', '<'), 
+c('>', '<='), c('>=', '<'),c('>=', '<=')).}
 
-\item{start}{An optional parameter to defined the initial date of the period 
+\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 
 date of the period and the initial 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{end}{An optional parameter to defined the final date of the period to 
+\item{end}{An optional parameter to define the final date of the period to 
 select from the data by providing a list of two elements: the final day of 
 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 
-dimension name matching the dimensions provided in the object 
-\code{data$data} can be specified.}
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
 \item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
 not (FALSE).}
@@ -50,17 +56,18 @@ not (FALSE).}
 computation.}
 }
 \value{
-An 's2dv_cube' object containing the indicator in the element 
-\code{data}.
+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.
 }
 \description{
-The Total Time of a variable exceeding (or not) a Threshold returns the total 
-number of days (if the data provided is daily, or the corresponding units to 
-the data frequency provided) that a variable is exceeding a threshold during a 
-period. The threshold provided must be in the same units than the variable 
-units, i.e. to use a percentile as a scalar, 
-the function \code{AbsToProbs} or \code{QThreshold}  may be needed (see 
-examples). Providing maximum temperature daily data, the following agriculture 
+The Total Time of a variable exceeding (or not) a Threshold. It returns the  
+total number of days (if the data provided is daily, or the corresponding  
+units of the data frequency) that a variable is exceeding a threshold 
+during a period. The threshold provided must be in the same units as the 
+variable units, i.e. to use a percentile as a scalar, the function 
+\code{AbsToProbs} or \code{QThreshold}  may be needed (see examples). 
+Providing maximum temperature daily data, the following agriculture 
 indices for heat stress can be obtained by using this function:
 \itemize{
  \item\code{SU35}{Total count of days when daily maximum temperatures exceed 
diff --git a/man/TotalSpellTimeExceedingThreshold.Rd b/man/TotalSpellTimeExceedingThreshold.Rd
index 37fd6cc0999f4bcf793bfa8da81248aa3758e87a..276423bb14092674ee655377d78a89e3a9f8796a 100644
--- a/man/TotalSpellTimeExceedingThreshold.Rd
+++ b/man/TotalSpellTimeExceedingThreshold.Rd
@@ -12,21 +12,29 @@ TotalSpellTimeExceedingThreshold(
   dates = NULL,
   start = NULL,
   end = NULL,
-  time_dim = "time",
+  time_dim = "ftime",
   ncores = NULL
 )
 }
 \arguments{
 \item{data}{A multidimensional array with named dimensions.}
 
-\item{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. If \code{timd_dim} is in the dimension (with the 
-same length as \code{data}), the comparison will be done day by day.}
+\item{threshold}{If only one threshold is used: it can be a multidimensional 
+array with named dimensions. It must be in the same units and with the 
+common dimensions of the same length as parameter 'data'. It can also be a
+vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+want to use two thresholds: it can be a vector of two scalars, a list of 
+two vectors with the same length of 'time_dim' from 'data' or a list of 
+two multidimensional arrays with the common dimensions of the same length 
+as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+also a vector of two elements.}
 
 \item{spell}{A scalar indicating the minimum length of the spell.}
 
-\item{op}{An operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+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 
@@ -38,22 +46,23 @@ date of the period and the initial 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{end}{An optional parameter to defined the final date of the period to 
+\item{end}{An optional parameter to define the final date of the period to 
 select from the data by providing a list of two elements: the final day of 
 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 
-dimension name matching the dimensions provided in the object 
-\code{data$data} can be specified.}
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
 \item{ncores}{An integer indicating the number of cores to use in parallel 
 computation.}
 }
 \value{
-A multidimensional array with named dimensions containing the indicator
-in the element \code{data}.
+A multidimensional array with named dimensions containing the number 
+of days that are part of a spell within a threshold with dimensions of the 
+input parameter 'data' except the dimension where the indicator has been 
+computed.
 }
 \description{
 The number of days (when daily data is provided) that are part of a spell 
@@ -77,7 +86,7 @@ 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, time = 20, lat = 4))
+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)
 
diff --git a/man/TotalTimeExceedingThreshold.Rd b/man/TotalTimeExceedingThreshold.Rd
index f874b50a3eac4e58a3d15ca287fe782c7e8c8176..206847574f25df8eabde09914d03bbd2d116bc46 100644
--- a/man/TotalTimeExceedingThreshold.Rd
+++ b/man/TotalTimeExceedingThreshold.Rd
@@ -11,7 +11,7 @@ TotalTimeExceedingThreshold(
   dates = NULL,
   start = NULL,
   end = NULL,
-  time_dim = "time",
+  time_dim = "ftime",
   na.rm = FALSE,
   ncores = NULL
 )
@@ -19,33 +19,39 @@ TotalTimeExceedingThreshold(
 \arguments{
 \item{data}{A multidimensional array with named dimensions.}
 
-\item{threshold}{A multidimensional array with named dimensions in the same 
-units as parameter \code{data} and with the common dimensions of the element 
-\code{data} of the same length (e.g. an array with the same lengths of 
-longitude and latitude). A single scalar is also possible (for the case of 
-comparing all grid points with the same scalar).}
+\item{threshold}{If only one threshold is used: it can be a multidimensional 
+array with named dimensions. It must be in the same units and with the 
+common dimensions of the same length as parameter 'data'. It can also be a
+vector with the same length of 'time_dim' from 'data' or a scalar. If we 
+want to use two thresholds: it can be a vector of two scalars, a list of 
+two vectors with the same length of 'time_dim' from 'data' or a list of 
+two multidimensional arrays with the common dimensions of the same length 
+as parameter 'data'. If two thresholds are used, parameter 'op' must be 
+also a vector of two elements.}
 
-\item{op}{A operator '>' (by default), '<', '>=' or '<='.}
+\item{op}{An operator '>' (by default), '<', '>=' or '<='. If  two thresholds
+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{start}{An optional parameter to defined the initial date of the period 
+\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 
 date of the period and the initial 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{end}{An optional parameter to defined the final date of the period to 
+\item{end}{An optional parameter to define the final date of the period to 
 select from the data by providing a list of two elements: the final day of 
 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 'time'. More than one 
-dimension name matching the dimensions provided in the object 
-\code{data$data} can be specified.}
+compute the indicator. By default, it is set to 'ftime'. It can only
+indicate one time dimension.}
 
 \item{na.rm}{A logical value indicating whether to ignore NA values (TRUE) or 
 not (FALSE).}
@@ -54,21 +60,23 @@ not (FALSE).}
 computation.}
 }
 \value{
-A multidimensional array with named dimensions containing the 
-indicator in the element \code{data}.
+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.
 }
 \description{
-The Total Time of a variable exceeding (or not) a Threshold returns the total 
-number of days (if the data provided is daily, or the corresponding units to 
-the data frequency provided) that a variable is exceeding a threshold during a 
-period. The threshold provided must be in the same units than the variable 
-units, i.e. to use a percentile as a threshold, the function \code{Threshold} 
-or \code{QThreshold} may be needed (see examples). Providing maximum 
-temperature daily data, the following agriculture indices for heat stress can 
-be obtained by using this function:
+The Total Time of a variable exceeding (or not) a Threshold. It returns the  
+total number of days (if the data provided is daily, or the corresponding  
+units of the data frequency) that a variable is exceeding a threshold 
+during a period. The threshold provided must be in the same units as the 
+variable units, i.e. to use a percentile as a scalar, the function 
+\code{AbsToProbs} or \code{QThreshold}  may be needed (see examples). 
+Providing maximum temperature daily data, the following agriculture 
+indices for heat stress can be obtained by using this function:
 \itemize{
  \item\code{SU35}{Total count of days when daily maximum temperatures exceed 
-                  35°C}
+                  35°C in the seven months from the start month given (e.g. 
+                  from April to October for start month of April).}
  \item\code{SU36}{Total count of days when daily maximum temperatures exceed 
                   36 between June 21st and September 21st}
  \item\code{SU40}{Total count of days when daily maximum temperatures exceed 
diff --git a/tests/testthat/test-AccumulationExceedingThreshold.R b/tests/testthat/test-AccumulationExceedingThreshold.R
index a0ad194cf45113454508980eb774d0276b440237..926fe1977e618beaaf5d1c240faec7396072284c 100644
--- a/tests/testthat/test-AccumulationExceedingThreshold.R
+++ b/tests/testthat/test-AccumulationExceedingThreshold.R
@@ -1,100 +1,306 @@
-context("Generic tests")
-test_that("Sanity checks", {
-  #source("csindicators/R/AccumulationExceedingThreshold.R")
-  expect_error(AccumulationExceedingThreshold(NULL), 
-    "Parameter 'data' cannot be NULL.")
-  expect_error(AccumulationExceedingThreshold('x'),
-   "Parameter 'data' must be numeric.")
-  data <- 1:20
-  expect_error(AccumulationExceedingThreshold(data, NULL),
-    "Parameter 'threshold' cannot be NULL.")
-  expect_error(AccumulationExceedingThreshold(data, 'x'),
-    "Parameter 'threshold' must be numeric.")
-  threshold <- 10
-  expect_equal(AccumulationExceedingThreshold(data, threshold), 155)
-  dim(data) <- c(2, 10)
-  expect_error(AccumulationExceedingThreshold(data, threshold),
-    "Parameter 'data' must have named dimensions.")
-  names(dim(data)) <- c('lat', 'time')
-  threshold <- array(1:2, 2)
-  expect_error(AccumulationExceedingThreshold(data, threshold),
-    "Parameter 'threshold' must have named dimensions.")
-  dim(threshold) <- c(time = 2)
-  data <- array(1:40, c(x = 2, ftime = 20))
-  expect_error(AccumulationExceedingThreshold(data, threshold),
-    "Could not find dimension 'time' in 1th object provided in 'data'.")    
-  threshold <- 10
-  expect_equal(AccumulationExceedingThreshold(data, threshold, time_dim = 'ftime'), 
-               array(c(375, 390), c(x = 2)))
-  dim(threshold) <- c(member = 1, ftime = 1)
-  expect_equal(AccumulationExceedingThreshold(data, threshold, time_dim = 'ftime'), 
-               array(c(375, 390), c(x = 2)))
-  expect_equal(AccumulationExceedingThreshold(data, threshold, time_dim = 'x'),
-               array(c(rep(0,5), seq(23, 79, 4)), c(ftime = 20)))
-  expect_error(AccumulationExceedingThreshold(data, threshold,
-                               time_dim = 'x', ncores = 'Z'),
-         "Parameter 'ncores' must be numeric")
-
-  expect_equal(AccumulationExceedingThreshold(data, threshold, time_dim = 2),
-         array(c(375, 390), c(x = 2)))
-  # dimensions:
-  data <- array(1:20, c(time = 5, sdate = 2, lat = 2))
-  # does this case made sense?
-  threshold <- array(1:5, c(time = 5))
-  expect_equal(dim(AccumulationExceedingThreshold(data, threshold)),
-               c(sdate = 2, lat = 2))
-  threshold <- array(1:2, c(lat = 2))
-  expect_equal(dim(AccumulationExceedingThreshold(data, threshold)),
-               c(sdate = 2, lat = 2))
-  data <- array(1:60, c(time = 5, fyear = 3, sdate = 2, lat = 2))
-  expect_equal(dim(AccumulationExceedingThreshold(data, threshold, 
-                      time_dim = c('time', 'fyear'))),
-               c(sdate = 2, lat = 2))
+context("CSIndicators::AccumulationExceedingThreshold tests")
+
+# dat1
+dat1 <- 1:20
+
+# dat2
+dat2_1 <- array(1:40, c(x = 2, ftime = 20))
+thres2_1 <- array(10, dim = c(member = 1, ftime = 1))
+dat2_3 <- array(1:20, c(ftime = 5, sdate = 2, lat = 2))
+thres2_3 <- array(1:5, c(ftime = 5))
+dat2_4 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+thres2_4 <- array(1:2, c(lat = 2))
+
+# dat3
+dat3_1 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+dat3_2 <- array(1:40, c(x = 2, ftime = 20))
+
+# dat4
+set.seed(1)
+dat4 <- array(rnorm(60, 23), c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+set.seed(1)
+thres4_1 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(2)
+thres4_2 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(1)
+thres4_3 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2))
+set.seed(2)
+thres4_4 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2))
+set.seed(1)
+thres4_5 <- array(rnorm(5, 20), c(ftime = 5))
+set.seed(2)
+thres4_6 <- array(rnorm(5, 25), c(ftime = 5))
+set.seed(1)
+thres4_7 <- rnorm(5, 20)
+set.seed(2)
+thres4_8 <- rnorm(5, 25)
+
+##############################################
+test_that("1. Input checks", {
+  # data
+  expect_error(
+    AccumulationExceedingThreshold(NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    AccumulationExceedingThreshold('x'),
+   "Parameter 'data' must be numeric."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(array(dat1, dim = c(2, 10)), 10),
+    "Parameter 'data' must have named dimensions."
+  )
+  # time_dim
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(array(dat1, dim = c('sdate' = 20)), 10),
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  # op
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = 1),
+    "Parameter 'op' must be a character."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = 'a'),
+    "Parameter 'op' must be a logical operator."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = c('=','=')),
+    "Parameter 'op' is not an accepted pair of logical operators."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = c('=','<','>')),
+    "Parameter 'op' must be a logical operator with length 1 or 2."
+  )
+  # threshold
+  expect_error(
+    AccumulationExceedingThreshold(dat1, NULL),
+    "Parameter 'threshold' cannot be NULL."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 'x'),
+    "Parameter 'threshold' must be numeric."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = c("<",">")),
+    "If 'op' is a  pair of logical operators parameter 'threshold' also has to be a pair of values."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, list(1:10,1:20), op = c("<",">")),
+    "The pair of thresholds must have the same length."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 5, lon = 2))), op = c("<",">")),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 2, lon = 5))), op = c("<",">")),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, list(array(1:10, c(time = 2)),array(1:10, c(time = 2))), op = c("<",">"), time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 1:10, op = "<"),
+    "If parameter 'threshold' is a vector it must have the same length as data time dimension."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, array(rnorm(10)), op = "<"),
+    "If parameter 'threshold' is an array it must have named dimensions."
+  )
+  expect_error(
+    AccumulationExceedingThreshold(dat1, array(20, dim = c(time = 2)), op = "<", time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  # ncores
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, time_dim = 'x', ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+  # dates
+  expect_error(
+    AccumulationExceedingThreshold(dat1, 10, op = "<", dates = 2, start = 'a', end = 'b'),
+    paste0("Parameter 'start' and 'end' must be lists indicating the ",
+           "day and the month of the period start and end.")
+  )
+  # diff
+  expect_error(
+    AccumulationExceedingThreshold(dat2_3, thres2_3, diff = T),
+    paste0("Parameter 'diff' can't be TRUE if the parameter 'threshold' is not a scalar.")
+  )
+})
+
+##############################################
+test_that("2. Output checks", {
+  expect_equal(
+    AccumulationExceedingThreshold(dat1, 10),
+    155
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat2_1, 10, time_dim = 'ftime'), 
+    array(c(375, 390), c(x = 2))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat2_1, thres2_1, time_dim = 'ftime'), 
+    array(c(375, 390), c(x = 2))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat2_1, thres2_1, time_dim = 'x'),
+    array(c(rep(0,5), seq(23, 79, 4)), c(ftime = 20))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat2_1, 10, time_dim = 'ftime'),
+    array(c(375, 390), c(x = 2))
+  )
+  # dimensions
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat2_3, thres2_3)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat2_3, thres2_4)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat2_4, thres2_4, time_dim = 'ftime')),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+
+})
+
+##############################################
+test_that("3. Output checks", {
+
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat3_1, c(55,58), c('<', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat3_1, c(55,58), c(">", "<")),
+    array(c(rep(0,11),113), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat3_1, c(55,58), c(">=", "<=")),
+    array(c(rep(0,10),55,171), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat3_2, c(46, 35), op = c("<", ">"), time_dim = 'ftime'), 
+    array(c(76, 114), c(x = 2))
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(dat3_2, c(7,11), op = c('>=', '<='), time_dim = 'ftime'), 
+    array(c(27, 18), c(x = 2))
+  )
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), op = c('>=', '<'))),
+    c(sdate = 2, lat = 2)
+  )
 
 })
 
-test_that("Seasonal forecasts", {
+##############################################
+
+test_that("4. Output checks", {
+
+  expect_equal(
+    dim(AccumulationExceedingThreshold(dat4, list(thres4_2, thres4_1), c('<=', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    as.vector(AccumulationExceedingThreshold(dat4, list(thres4_1, thres4_2), c(">", "<="))[1:3]),
+    c(91.05107, 115.67568, 69.89353),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AccumulationExceedingThreshold(dat4, list(thres4_3, thres4_4), c(">", "<="), time_dim = 'ftime'))[1:5],
+    c(91.05107, 115.67568, 69.89353, 117.29783, 115.40615),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AccumulationExceedingThreshold(dat4, list(thres4_6, thres4_5), op = c("<", ">="), time_dim = 'ftime'))[1:5], 
+    c(91.05107, 115.67568, 69.89353, 117.29783, 94.39550),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AccumulationExceedingThreshold(dat4, list(thres4_7, thres4_8), op = c('>=', '<='), time_dim = 'ftime'))[4:10], 
+    c(117.29783, 94.39550, 113.25711, 90.85402, 91.89458, 115.14699, 116.19438),
+    tolerance = 0.0001
+  )
+
+})
 
+##############################################
+test_that("5. Seasonal forecasts", {
+  library(CSTools)
   exp <- CSTools::lonlat_temp$exp
-  exp$data <- exp$data[,1:4,1:2,,,]
-  res <- CST_AccumulationExceedingThreshold(exp, threshold = 280)
-  expect_equal(round(res$data[,2,2,2]),
-               c(0, 280, 281, 281))
+  exp$data <- exp$data[ , 1:4, 1:2, , , ]
+  res <- CST_AccumulationExceedingThreshold(exp, threshold = 280, time_dim = 'ftime')
+
+  expect_equal(
+    round(res$data[, 2, 2, 2]),
+    c(0, 280, 281, 281)
+  )
+
   # GDD
-  exp <- array(NA, dim = c(member = 6, sdate = 3, ftime = 214, lat =4, lon = 4))
+  exp <- array(NA, dim = c(member = 6, sdate = 3, ftime = 214, lat = 4, lon = 4))
   exp1 <- drop(CSTools::lonlat_prec$data) * 86400000
-  exp[,,1:31,,] <- exp1 + 10; exp[,,32:62,,] <- exp1 + 11
-  exp[,,63:93,,] <- exp1 + 12; exp[,,94:124,,] <- exp1 + 13
-  exp[,,125:155,,] <- exp1 + 14; exp[,,156:186,,] <- exp1 + 15
-  exp[,,187:214,,] <- exp1[,,1:28,,] + 16
+  exp[, , 1:31, , ] <- exp1 + 10; exp[, , 32:62, , ] <- exp1 + 11
+  exp[, , 63:93, , ] <- exp1 + 12; exp[, , 94:124, , ] <- exp1 + 13
+  exp[, , 125:155, , ] <- exp1 + 14; exp[, , 156:186, , ] <- exp1 + 15
+  exp[, , 187:214, , ] <- exp1[, , 1:28, , ] + 16
+
   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'))
+                 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(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))
-  expect_equal(dim(GDD),
-              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)),
-               "Could not find dimension 'time' in 1th object provided in 'data'.")
-  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))))
+  
+  expect_equal(
+    round(GDD[,1,1,1]),
+    c(538, 367, 116, 519, 219, 282)
+  )
+  expect_equal(
+    dim(GDD),
+    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)))
+  )
 
   # test the 'diff'
-  input <- c(1:20)
-  threshold <- 3
-  expect_equal(AccumulationExceedingThreshold(input, threshold, diff = TRUE),
-               153)
-  expect_equal(AccumulationExceedingThreshold(input, threshold),
-               204)
-  input1 <- -input[1:15]
-  threshold <- -5
-  expect_equal(AccumulationExceedingThreshold(input1, threshold, op = '<'),
-               -105)
-  expect_equal(AccumulationExceedingThreshold(input1, threshold, op = '<', diff = TRUE),
-               -55)
+  input_1 <- c(1:20)
+  threshold_1 <- 3
+  input_2 <- -input_1[1:15]
+  threshold_2 <- -5
+
+  expect_equal(
+    AccumulationExceedingThreshold(input_1, threshold_1, diff = TRUE),
+    153
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(input_1, threshold_1),
+    204
+  )
+
+  expect_equal(
+    AccumulationExceedingThreshold(input_2, threshold_2, op = '<'),
+    -105
+  )
+  expect_equal(
+    AccumulationExceedingThreshold(input_2, threshold_2, op = '<', diff = TRUE),
+    -55
+  )
 })
diff --git a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
index 4c55098dac48293136dcb6a8fbae1ce1230a1256..d2155298df51c082450d75c7d449293020ff3c07 100644
--- a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
@@ -1,55 +1,265 @@
-context("Generic tests")
-test_that("Sanity checks", {
-  #source("csindicators/R/TotalSpellTimeExceedingThreshold.R")
-  expect_error(TotalSpellTimeExceedingThreshold(NULL), 
-    "Parameter 'data' cannot be NULL.")
-  expect_error(TotalSpellTimeExceedingThreshold('x'),
-   "Parameter 'data' must be numeric.")
-  data <- 1:20
-  expect_error(TotalSpellTimeExceedingThreshold(data, NULL),
-    "Parameter 'threshold' cannot be NULL.")
-  expect_error(TotalSpellTimeExceedingThreshold(data, 'x'),
-    "Parameter 'threshold' must be numeric.")
-  threshold <- 10
-  expect_error(TotalSpellTimeExceedingThreshold(data, threshold),
-    paste("argument", '"spell"', "is missing, with no default"))
-  dim(data) <- c(2, 10)
-  expect_error(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2),
-    "Parameter 'data' must have named dimensions.")
-  names(dim(data)) <- c('time', 'lat')
-  threshold <- array(1:2, 2)
-  dim(threshold) <- c(time = 2)
-  expect_equal(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2),  
-    array(c(0,rep(2,9)), c(lat = 10)))
-  data <- array(1:40, c(x = 2, ftime = 20))
-  expect_error(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2),
-    "Could not find dimension 'time' in 1th object provided in 'data'.")    
-  threshold <- 10
-  expect_equal(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2, time_dim = 'ftime'), 
-               array(c(15, 15), c(x = 2)))
-  threshold <- rep(10, 20)
-  dim(threshold) <- c(member = 1, ftime = 20)
-  expect_equal(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2, time_dim = 'ftime'), 
-               array(c(15, 15), c(x = 2, member = 1)))
-  expect_error(TotalSpellTimeExceedingThreshold(data, threshold, spell = 2, time_dim = 'y'),
-     paste("Could not find dimension 'y' in 1th object provided in 'data'"))
+context("CSIndicators::TotalSpellTimeExceedingThreshold tests")
+
+# dat1
+dat <- array(1:20, dim = c(2, 10))
+thres <- 10
+dat1 <- array(1:20, dim = c(time = 2, lat = 10))
+thres1 <- array(1:2, dim = c(time = 2))
+dat1_2 <- array(1:40, c(x = 2, ftime = 20))
+threshold1_2 <- array(rep(10, 20), dim = c(member = 1, ftime = 20))
+
+# dat2
+dat2_1 <- array(1:40, c(x = 2, ftime = 20))
+thres2_1 <- array(10, dim = c(member = 1, ftime = 1))
+dat2_3 <- array(1:20, c(ftime = 5, sdate = 2, lat = 2))
+thres2_3 <- array(1:5, c(ftime = 5))
+dat2_4 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+thres2_4 <- array(1:2, c(lat = 2))
+
+# dat3
+dat3_1 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+dat3_2 <- array(1:40, c(x = 2, ftime = 20))
+
+# dat4
+set.seed(1)
+dat4 <- array(rnorm(60, 23), c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+set.seed(1)
+thres4_1 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(2)
+thres4_2 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(1)
+thres4_3 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2))
+set.seed(2)
+thres4_4 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2))
+set.seed(1)
+thres4_5 <- array(rnorm(5, 20), c(ftime = 5))
+set.seed(2)
+thres4_6 <- array(rnorm(5, 25), c(ftime = 5))
+set.seed(1)
+thres4_7 <- rnorm(5, 20)
+set.seed(2)
+thres4_8 <- rnorm(5, 25)
+
+###########################################################################
+test_that("1. Sanity checks", {
+  # data
+  expect_error(
+    TotalSpellTimeExceedingThreshold(NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold('x'),
+    "Parameter 'data' must be numeric."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(array(dat1, dim = c(2, 10)), 10),
+    "Parameter 'data' must have named dimensions."
+  )
+  # time_dim
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(array(dat1, dim = c('sdate' = 20)), 10),
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  # op
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = 1, time_dim = 'time'),
+    "Parameter 'op' must be a character."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = 'a', time_dim = 'time'),
+    "Parameter 'op' must be a logical operator."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = c('=','='), time_dim = 'time'),
+    "Parameter 'op' is not an accepted pair of logical operators."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = c('=','<','>'), time_dim = 'time'),
+    "Parameter 'op' must be a logical operator with length 1 or 2."
+  )
+  # threshold
+  expect_error(
+    TotalSpellTimeExceedingThreshold(1:20, NULL),
+    "Parameter 'threshold' cannot be NULL."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(1:20, 'x'),
+    "Parameter 'threshold' must be numeric."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = c("<",">"), spell = 6,time_dim = 'time'),
+    "If 'op' is a  pair of logical operators parameter 'threshold' also has to be a pair of values."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, list(1:10,1:20), op = c("<",">"), spell = 6, time_dim = 'time'),
+    "The pair of thresholds must have the same length."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 5, lon = 2))), op = c("<",">"), spell = 6, time_dim = 'time'),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 2, lon = 5))), op = c("<",">"), spell = 6, time_dim = 'time'),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, list(array(1:3, c(time = 3)),array(1:3, c(time = 3))), op = c("<",">"), spell = 6, time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 1:10, spell = 6, op = "<", time_dim = 'time'),
+    "If parameter 'threshold' is a vector it must have the same length as data time dimension."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, array(rnorm(10)), spell = 6, op = "<", time_dim = 'time'),
+    "If parameter 'threshold' is an array it must have named dimensions."
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, array(3, dim = c(time = 3)), op = "<", spell = 6, time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  # spell
+  expect_error(
+    TotalSpellTimeExceedingThreshold(1:20, thres),
+    paste("argument", '"spell"', "is missing, with no default")
+  )
+  expect_error(
+    TotalSpellTimeExceedingThreshold(1:20, thres, spell = 1:2),
+    "Parameter 'spell' must be a scalar."
+  )
+  # ncores
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, time_dim = 'time', spell = 6, ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+  # dates
+  expect_error(
+    TotalSpellTimeExceedingThreshold(dat1, 10, op = "<", dates = 2, spell = 6, start = 'a', end = 'b', time_dim = 'time'),
+    paste0("Parameter 'start' and 'end' must be lists indicating the ",
+           "day and the month of the period start and end.")
+  )
+})
+
+###########################################################################
+
+test_that("2. Output checks", {
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat1, thres1, spell = 2, time_dim = 'time'),  
+    array(c(0,rep(2,9)), c(lat = 10))
+  )    
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat1_2, 10, spell = 2), 
+    array(c(15, 15), c(x = 2))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat1_2,threshold1_2, spell = 2), 
+    array(c(15, 15), c(x = 2, member = 1))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat2_1, thres2_1, spell = 10), 
+    array(c(15, 15), c(x = 2))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat2_1, thres2_1, spell = 2, time_dim = 'x'),
+    array(c(rep(0,5), rep(2,15)), c(ftime = 20))
+  )
+  # dimensions
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat2_3, thres2_3, spell = 3)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat2_3, thres2_4, spell = 3)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat2_4, thres2_4, spell = 3, time_dim = 'ftime')),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+})
+
+##############################################
+
+test_that("3. Output checks", {
+
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat3_1, c(55,58), spell = 3, c('<', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat3_1, c(30,60), spell = 3, c(">", "<")),
+    array(c(rep(0,6),rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat3_1, c(55,58), spell = 3, c(">=", "<=")),
+    array(c(rep(0,11),3), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat3_2, c(46, 35), spell = 3, op = c("<", ">"), time_dim = 'ftime'), 
+    array(c(0, 3), c(x = 2))
+  )
+  expect_equal(
+    TotalSpellTimeExceedingThreshold(dat3_2, c(7,11), spell = 3, op = c('>=', '<='), time_dim = 'ftime'), 
+    array(c(3, 0), c(x = 2))
+  )
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), spell = 3, op = c('>=', '<'))),
+    c(sdate = 2, lat = 2)
+  )
+
+})
+
+##############################################
+
+test_that("4. Output checks", {
+
+  expect_equal(
+    dim(TotalSpellTimeExceedingThreshold(dat4, list(thres4_2, thres4_1), spell = 3, c('<=', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    as.vector(TotalSpellTimeExceedingThreshold(dat4, list(thres4_1, thres4_2), spell = 3, c(">", "<="))[1:3]),
+    c(3, 5, 0)
+  )
+  expect_equal(
+    as.vector(TotalSpellTimeExceedingThreshold(dat4, list(thres4_3, thres4_4), spell = 4, c(">", "<="), time_dim = 'ftime'))[1:5],
+    c(0, 5, 0, 5, 5)
+  )
+  expect_equal(
+    as.vector(TotalSpellTimeExceedingThreshold(dat4, list(thres4_6, thres4_5), spell = 3, op = c("<", ">="), time_dim = 'ftime'))[1:5], 
+    c(3, 5, 0, 5, 3)
+  )
+  expect_equal(
+    as.vector(TotalSpellTimeExceedingThreshold(dat4, list(thres4_7, thres4_8), spell = 3, op = c('>=', '<='), time_dim = 'ftime'))[4:10], 
+    c(5, 3, 5, 4, 3, 5, 5)
+  )
 
 })
 
-test_that("Seasonal Forecasts", {
+###########################################################################
 
+test_that("5. Seasonal Forecasts", {
   exp <- CSTools::lonlat_temp$exp
   exp$data <- exp$data[1,1:3,1:3,,,]
   res <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 260, spell = 2)
-  expect_equal(res$data[,,1,1],
-               array(c(3,3,3,3,3,3,3,3,3), c(member = 3, sdate = 3)))
+  expect_equal(
+    res$data[,,1,1],
+    array(c(3,3,3,3,3,3,3,3,3), c(member = 3, sdate = 3))
+  )
   # compare with percentile
   thresholdP <- Threshold(exp$data, threshold = 0.9)
   WSDI <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = thresholdP, spell = 2)
-  expect_equal(WSDI$data[3,3,3,],
-               c(rep(0, 6), rep(2, 2), 0, 2, 0, 2, rep(0, 29), 2, rep(0, 11)))
+  expect_equal(
+    WSDI$data[3,3,3,],
+    c(rep(0, 6), rep(2, 2), 0, 2, 0, 2, rep(0, 29), 2, rep(0, 11))
+  )
   thresholdP1 <- thresholdP[1,,]
   WSDI1 <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = thresholdP1, spell = 2)
-  expect_equal(WSDI1$data[3,3,3,],
-               c(rep(0, 53)))
+  expect_equal(
+    WSDI1$data[3,3,3,],
+    c(rep(0, 53)))
 })
diff --git a/tests/testthat/test-TotalTimeExceedingThreshold.R b/tests/testthat/test-TotalTimeExceedingThreshold.R
index a1104b461f7e5a4cb443de833517a743abbcbb2e..4a5a365d65439e56f321e74a694e661fe2377e92 100644
--- a/tests/testthat/test-TotalTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalTimeExceedingThreshold.R
@@ -1,74 +1,252 @@
-context("Generic tests")
-test_that("Sanity checks", {
-  #source("csindicators/R/TotalTimeExceedingThreshold.R")
-  expect_error(TotalTimeExceedingThreshold(NULL), 
-    "Parameter 'data' cannot be NULL.")
-  expect_error(TotalTimeExceedingThreshold('x'),
-   "Parameter 'data' must be numeric.")
-  data <- 1:20
-  expect_error(TotalTimeExceedingThreshold(data, NULL),
-    "Parameter 'threshold' cannot be NULL.")
-  expect_error(TotalTimeExceedingThreshold(data, 'x'),
-    "Parameter 'threshold' must be numeric.")
-  threshold <- 10
-  expect_equal(TotalTimeExceedingThreshold(data, threshold), 10)
-  dim(data) <- c(2, 10)
-  expect_error(TotalTimeExceedingThreshold(data, threshold),
-    "Parameter 'data' must have named dimensions.")
-  names(dim(data)) <- c('lat', 'time')
-  threshold <- array(1:2, 2)
-  expect_error(TotalTimeExceedingThreshold(data, threshold),
-    "Parameter 'threshold' must have named dimensions.")
-  dim(threshold) <- c(time = 2)
-  expect_error(TotalTimeExceedingThreshold(data, threshold),  
-    "Parameter 'data' and 'threshold' must have the same length on common dimensions.")
-  data <- array(1:40, c(x = 2, ftime = 20))
-  expect_error(TotalTimeExceedingThreshold(data, threshold),
-    "Could not find dimension 'time' in 1th object provided in 'data'.")    
-  threshold <- 10
-  expect_equal(TotalTimeExceedingThreshold(data, threshold, time_dim = 'ftime'), 
-               array(c(15, 15), c(x = 2)))
-  dim(threshold) <- c(member = 1, ftime = 1)
-  expect_equal(TotalTimeExceedingThreshold(data, threshold, time_dim = 'ftime'), 
-               array(c(15, 15), c(x = 2)))
-  expect_equal(TotalTimeExceedingThreshold(data, threshold, time_dim = 'x'),
-               array(c(rep(0,5), rep(2,15)), c(ftime = 20)))
-  expect_error(TotalTimeExceedingThreshold(data, threshold,
-                               time_dim = 'x', ncores = 'Z'),
-         "Parameter 'ncores' must be numeric")
+context("CSIndicators::TotalTimeExceedingThreshold tests")
 
-  expect_equal(TotalTimeExceedingThreshold(data, threshold, time_dim = 2),
-         array(c(15,15), c(x = 2)))
-  # dimensions:
-  data <- array(1:20, c(time = 5, sdate = 2, lat = 2))
-  # does this case made sense?
-  threshold <- array(1:5, c(time = 5))
-  expect_equal(dim(TotalTimeExceedingThreshold(data, threshold)),
-               c(sdate = 2, lat = 2))
-  threshold <- array(1:2, c(lat = 2))
-  expect_equal(dim(TotalTimeExceedingThreshold(data, threshold)),
-               c(sdate = 2, lat = 2))
-  data <- array(1:60, c(time = 5, fyear = 3, sdate = 2, lat = 2))
-  expect_equal(dim(TotalTimeExceedingThreshold(data, threshold, 
-                      time_dim = c('time', 'fyear'))),
-               c(sdate = 2, lat = 2))
+# dat1
+dat <- array(1:20, dim = c(2, 10))
+thres <- 10
+dat1 <- array(1:20, dim = c(time = 2, lat = 10))
+thres1 <- array(1:2, dim = c(time = 2))
+dat1_2 <- array(1:40, c(x = 2, ftime = 20))
+threshold1_2 <- array(rep(10, 20), dim = c(member = 1, ftime = 20))
 
+# dat2
+dat2_1 <- array(1:40, c(x = 2, ftime = 20))
+thres2_1 <- array(10, dim = c(member = 1, ftime = 1))
+dat2_3 <- array(1:20, c(ftime = 5, sdate = 2, lat = 2))
+thres2_3 <- array(1:5, c(ftime = 5))
+dat2_4 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+thres2_4 <- array(1:2, c(lat = 2))
+
+# dat3
+dat3_1 <- array(1:60, c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+dat3_2 <- array(1:40, c(x = 2, ftime = 20))
+
+# dat4
+set.seed(1)
+dat4 <- array(rnorm(60, 23), c(ftime = 5, fyear = 3, sdate = 2, lat = 2))
+set.seed(1)
+thres4_1 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(2)
+thres4_2 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2, lat = 2))
+set.seed(1)
+thres4_3 <- array(rnorm(20, 20), c(ftime = 5, sdate = 2))
+set.seed(2)
+thres4_4 <- array(rnorm(20, 25), c(ftime = 5, sdate = 2))
+set.seed(1)
+thres4_5 <- array(rnorm(5, 20), c(ftime = 5))
+set.seed(2)
+thres4_6 <- array(rnorm(5, 25), c(ftime = 5))
+set.seed(1)
+thres4_7 <- rnorm(5, 20)
+set.seed(2)
+thres4_8 <- rnorm(5, 25)
+
+###########################################################################
+
+test_that("1. Sanity checks", {
+  # data
+  expect_error(
+    TotalTimeExceedingThreshold(NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold('x'),
+    "Parameter 'data' must be numeric."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(array(dat1, dim = c(2, 10)), 10),
+    "Parameter 'data' must have named dimensions."
+  )
+  # time_dim
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(array(dat1, dim = c('sdate' = 20)), 10),
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  # op
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = 1, time_dim = 'time'),
+    "Parameter 'op' must be a character."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = 'a', time_dim = 'time'),
+    "Parameter 'op' must be a logical operator."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = c('=','='), time_dim = 'time'),
+    "Parameter 'op' is not an accepted pair of logical operators."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = c('=','<','>'), time_dim = 'time'),
+    "Parameter 'op' must be a logical operator with length 1 or 2."
+  )
+  # threshold
+  expect_error(
+    TotalTimeExceedingThreshold(1:20, NULL),
+    "Parameter 'threshold' cannot be NULL."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(1:20, 'x'),
+    "Parameter 'threshold' must be numeric."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = c("<",">"), time_dim = 'time'),
+    "If 'op' is a pair of logical operators parameter 'threshold' also has to be a pair of values."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, list(1:10,1:20), op = c("<",">"), time_dim = 'time'),
+    "The pair of thresholds must have the same length."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 5, lon = 2))), op = c("<",">"), time_dim = 'time'),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, list(array(1:10, c(lat = 2, lon = 5)),array(1:10, c(time = 2, lon = 5))), op = c("<",">"), time_dim = 'time'),
+    "The pair of thresholds must have the same dimension names."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, list(array(1:3, c(time = 3)),array(1:3, c(time = 3))), op = c("<",">"), time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 1:10, op = "<", time_dim = 'time'),
+    "If parameter 'threshold' is a vector it must have the same length as data time dimension."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, array(rnorm(10)), op = "<", time_dim = 'time'),
+    "If parameter 'threshold' is an array it must have named dimensions."
+  )
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, array(3, dim = c(time = 3)), op = "<", time_dim = 'time'),
+    "Parameter 'data' and 'threshold' must have same length of all common dimensions."
+  )
+  # ncores
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, time_dim = 'time', ncores = 1.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+  # dates
+  expect_error(
+    TotalTimeExceedingThreshold(dat1, 10, op = "<", dates = 2, start = 'a', end = 'b', time_dim = 'time'),
+    paste0("Parameter 'start' and 'end' must be lists indicating the ",
+           "day and the month of the period start and end.")
+  )
+})
+
+
+###########################################################################
+
+test_that("2. Output checks", {
+  expect_equal(
+    TotalTimeExceedingThreshold(dat1, thres1, time_dim = 'time'),  
+    array(c(0,rep(2,9)), c(lat = 10))
+  )    
+  expect_equal(
+    TotalTimeExceedingThreshold(dat1_2, 10), 
+    array(c(15, 15), c(x = 2))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat1_2,threshold1_2), 
+    array(c(15, 15), c(x = 2, member = 1))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat2_1, thres2_1), 
+    array(c(15, 15), c(x = 2))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat2_1, thres2_1, time_dim = 'x'),
+    array(c(rep(0,5), rep(2,15)), c(ftime = 20))
+  )
+  # dimensions
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat2_3, thres2_3)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat2_3, thres2_4)),
+    c(sdate = 2, lat = 2)
+  )
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat2_4, thres2_4, time_dim = 'ftime')),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
 })
 
+##############################################
+
+test_that("3. Output checks", {
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat3_1, c(55,58), c('<', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat3_1, c(30,60), c(">", "<")),
+    array(c(rep(0, 6), rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat3_1, c(55, 58), c(">=", "<=")),
+    array(c(rep(0, 10), 1, 3), dim = c(fyear = 3, sdate = 2, lat = 2))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat3_2, c(46, 35), op = c("<", ">"), time_dim = 'ftime'), 
+    array(c(2, 3), c(x = 2))
+  )
+  expect_equal(
+    TotalTimeExceedingThreshold(dat3_2, c(7, 11), op = c('>=', '<='), time_dim = 'ftime'), 
+    array(c(3, 2), c(x = 2))
+  )
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), op = c('>=', '<'))),
+    c(sdate = 2, lat = 2)
+  )
+})
+
+##############################################
+
+test_that("4. Output checks", {
+  expect_equal(
+    dim(TotalTimeExceedingThreshold(dat4, list(thres4_2, thres4_1), c('<=', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2)
+  )
+  expect_equal(
+    as.vector(TotalTimeExceedingThreshold(dat4, list(thres4_1, thres4_2), c(">", "<="))[1:3]),
+    c(4, 5, 3)
+  )
+  expect_equal(
+    as.vector(TotalTimeExceedingThreshold(dat4, list(thres4_3, thres4_4), c(">", "<="), time_dim = 'ftime'))[1:5],
+    c(4, 5, 3, 5, 5)
+  )
+  expect_equal(
+    as.vector(TotalTimeExceedingThreshold(dat4, list(thres4_6, thres4_5), op = c("<", ">="), time_dim = 'ftime'))[1:5], 
+    c(4, 5, 3, 5, 4)
+  )
+  expect_equal(
+    as.vector(TotalTimeExceedingThreshold(dat4, list(thres4_7, thres4_8), op = c('>=', '<='), time_dim = 'ftime'))[4:10], 
+    c(5, 4, 5, 4, 4, 5, 5)
+  )
+})
+
+###########################################################################
+
 test_that("Seasonal forecasts", {
   # compare with scalar fixed threshold
   exp <- CSTools::lonlat_temp$exp
-  exp$data <- exp$data[1,1:3,,,,] - 247
-  SU35_NoP <- CST_TotalTimeExceedingThreshold(exp, threshold = 35)$data  
-  expect_equal(SU35_NoP[1,,15,3], c(0,1,1,1,0,0))
-  # convert to percentile
   obs <- CSTools::lonlat_temp$obs
-  exp_percentile <- AbsToProbs(exp$data)
-  obs_percentile <- drop(QThreshold(obs$data, threshold = 35))
+  exp$data <- exp$data[1, 1:3, , , , ] - 247
+  SU35_NoP <- CST_TotalTimeExceedingThreshold(exp, threshold = 35)$data  
+  expect_equal(
+    SU35_NoP[1, , 15, 3], c(0, 1, 1, 1, 0, 0))
+    # convert to percentile
+    exp_percentile <- AbsToProbs(exp$data)
+    obs_percentile <- drop(QThreshold(obs$data, threshold = 35)
+  )
   data <- exp
   data$data <- exp_percentile
   SU35_P <- CST_TotalTimeExceedingThreshold(data, threshold = obs_percentile)$data
-  expect_equal(SU35_P[,2,5,5],  c(3,3,3,3,3,3))
+  expect_equal(
+    SU35_P[ ,2, 5, 5], c(3, 3, 3, 3, 3, 3)
+  )
 })
-
-