diff --git a/NAMESPACE b/NAMESPACE
index 03fa0de4a2997108d8994eac872d50842c898dab..9b7ff9058f6bbd726800fbf151286383cc2c3f15 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(AnoAgree)
+export(BioIndicators)
 export(Climdex)
 export(CombineIndices)
 export(DTRIndicator)
@@ -8,6 +9,7 @@ export(DTRRef)
 export(DailyAno)
 export(Extremes)
 export(Lon2Index)
+export(QThreshold)
 export(SeasonSelect)
 export(SelBox)
 export(Subset)
diff --git a/R/BioIndicators.R b/R/BioIndicators.R
new file mode 100644
index 0000000000000000000000000000000000000000..0ee567504146083c860dc15c6f1c38b523600002
--- /dev/null
+++ b/R/BioIndicators.R
@@ -0,0 +1,145 @@
+#' Bioclimate indicators
+#'
+#'Different bioclimate indicators computation for agriculture:
+#'\itemize{
+#'  \item\code{GDD}{Growing Degree Days is defined as the summation of daily differences between daily average temperatures and 10°C between April 1st and October 31st}
+#'  \item\code{GST}{Growing Season Temperature is the average of the daily average temperatures for a given period. The daily average temperature here is computed by (daily temperature maximum +  daily temperature minimum) / 2.}
+#'  \item\code{SPRTX}{Sprtx is defined as the average of daily temperature maximums between April 1st and May 31st.}
+#'  \item\code{TP}{Total Precipitation sum of daily precipitation during a certain period. It allows to compute the specific indices \code{SprR} (Spring Rain takes the days from April 21st to June 21st) and \code{HarvestR} (Harvest Rain takes the days from August 21st to October 21st).}
+#'  \item\code{SU}{Total count of days when daily maximum temperatures exceed a threshold. It allows to compute specific indices \code{SU36} (when 36°C are exceeded between June 21st and September 21st) and \code{SU40} (when 40°C are exceeded between June 21st and September 21st).}
+#'}
+#' The user must provide the data for the proper variable and period depending on the desired index.
+#'
+#'@param data a n-dimensional array or a list of n-dimensional arrays with named dimensions containing daily data depending on each specific index:
+#'  \item\code{GDD}{data must be a list of two n-dimensional arrays for maximum and minimum daily temperatures in this order}
+#'  \item\code{GST}{data must be a list of two n-dimensional arrays for maximum and minimum daily temperatures in this order}
+#'  \item\code{SPRTX}{data must be a n-dimensional array for maximum daily temperatures}
+#'  \item\code{TP}{data must be a n-dimensional array fir daily precipitation.
+#'  \item\code{SU}{data must be a list of one n-dimensinoal array containing temperature data (usually daily maximum or minimum) and a scalar, vector or n-dimensional array containing the threshold}
+#'} 
+#'@param index a character string indicating the index to be computed among \code{GDD}, \code{GST}, \code{SPRTX}, \code{TP} and \code{SU}.
+#'@param time_dim a character string indicting the name of the temporal dimension.
+#'@param na.rm a logical value indicting if missing values should be removed.
+#'@param ncores an integer indicating the number of cores to use in parallel computation.
+#'@import multiApply
+#'@details 
+#' Requires extra definition?
+#'@examples
+#'data <- array(1:10, c(ftime = 5, x = 2))
+#'BioIndicators(data, index = 'TP') 
+#'@seealso \code{QThreshold}
+#'@export
+BioIndicators <- function(data, index, time_dim = 'ftime', na.rm = TRUE,
+                          ncores = NULL) {
+  if (!is.character(time_dim)) {
+    stop("Parameter 'time_dim' must be a character string indicating the",
+ " name of the dimension to compute the index.")
+  }
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores)) {
+      ncores <- 1
+      warning("Parameter 'ncores' must be numeric. Execution using 1 core.")
+    }
+    if (length(ncores) > 1) {
+      ncores <- ncores[1]
+      warning("Parameter 'ncores' has length greater than one and only",
+              " the first element will be used.")
+    }
+  }
+  if (!is.logical(na.rm)) {
+    na.rm <- TRUE
+    warning("Parameter 'na.rm' must be logical and it has been set to TRUE.")
+  }
+  if (length(na.rm) > 1) {
+    na.rm <- na.rm[1]
+    warning("Parameter 'na.rm' has length greater than one and only",
+            " the first element will be used.")
+  }
+  if (is.null(index)) {
+    stop("Parameter 'index' must be specified.")
+  }
+  if (length(index) > 1) {
+    index <- index[1]
+    warning("Parameter 'index' has length greater than one and only",
+            " the first element will be used.")
+  }  
+  if (!is.character(index)) {
+    stop("Parameter 'index' must be a character string indicating one of the ",
+         "available indices for computation.")
+  }
+  if (index == "GDD") {
+    if (!is.list(data) || length(data) != 2) {
+      stop("Parameter 'data' must be a list of two n-dimensional arrays.")
+    }
+    res <- Apply(data, target_dims = list(time_dim, time_dim), fun = .gdd,
+                 na.rm = na.rm, ncores = ncores)$output1
+  } else if (index == "GST") {
+    if (!is.list(data) || length(data) != 2) {
+      stop("Parameter 'data' must be a list of two n-dimensional arrays.")
+    }
+     res <- Apply(data, target_dims = list(time_dim, time_dim), 
+                  fun = .meantxtn, na.rm = na.rm, ncores = ncores)$output1
+  } else if (index == "SPRTX") {
+    if (is.list(data)) {
+      data <- data[[1]]
+    }
+    if (is.null(dim(data))) {
+      stop("Parameter 'data' must be a n-dimensional array.")
+    }
+    if (is.null(names(dim(data)))) {
+      stop("Parameter 'data' must be a n-dimensional array with named dimensions.")
+    }
+    res <-  Apply(list(data), target_dims = time_dim, fun = mean, 
+                  na.rm = na.rm, ncores = ncores)$output1
+  } else if (index == "TP") {
+    if (is.list(data)) {
+      data <- data[[1]]
+    }
+    if (is.null(dim(data))) {
+      stop("Parameter 'data' must be a n-dimensional array.")
+    }
+    if (is.null(names(dim(data)))) {
+      stop("Parameter 'data' must be a n-dimensional array with named dimensions.")
+    }
+    res <- Apply(list(data), target_dims = time_dim, fun = sum,
+                 na.rm = na.rm, ncores = ncores)$output1
+  } else if (index == "SU") {
+    if (!is.list(data) || length(data) != 2) {
+      stop("Parameter 'data' must be a list of two elements: one n-dimensional",
+           " array and one scalar, vector or n-dimensional array.")
+    }
+    threshold <- data[[2]]
+    data <- data[[1]]
+    if (is.null(dim(data))) {
+      dim(data) <- c(length(data))
+      names(dim(data)) <- time_dim
+    }
+    if (is.null(dim(threshold)) && length(threshold) > 1) {
+      dim(threshold) <- c(length(threshold))
+      names(dim(threshold)) <- time_dim 
+    }
+    threshold <- drop(threshold) 
+    if (!is.null(dim(threshold)) && !all(dim(threshold) == 1)) {
+      res <- Apply(list(data, threshold), target_dims = list(time_dim, time_dim),
+                   fun = .su, na.rm = na.rm, ncores = ncores)$output1
+    } else {
+      res <- Apply(list(data), target_dims = list(time_dim),
+                   fun = .su, threshold = threshold, na.rm = na.rm,
+                   ncores = ncores)$output1
+    }
+  } else {
+    stop("Parameter 'index' must be a character string among:",
+         " 'GDD', 'GST', 'SPRTX', 'TP' or 'SU'.")
+  }  
+  return(res)
+}    
+.gdd <- function(tx, tn, na.rm = na.rm) {
+  temp <- ((tx + tn) / 2) - 10
+  res <- sum(temp[temp > 0], na.rm = na.rm)
+}
+.meantxtn <- function(tx, tn, na.rm = na.rm) {
+  ta <- mean((tx + tn) / 2, na.rm = na.rm)
+}
+.su <- function(data, threshold, na.rm) {
+  result <- sum(data > threshold, na.rm = na.rm)
+}
diff --git a/R/CombineIndices.R b/R/CombineIndices.R
index 05f04c98b02c4ce59580296ab1c56da500743405..8728c08201a6453e72f66f639b4842aeca236470 100644
--- a/R/CombineIndices.R
+++ b/R/CombineIndices.R
@@ -85,4 +85,4 @@ CombineIndices <- function(indices, weights = NULL, operation = "mean") {
     }
   }
   comb_ind
-}
\ No newline at end of file
+}
diff --git a/R/GDD.R b/R/GDD.R
deleted file mode 100644
index 47995d947fb351b8ea537eb26a425065a0cf8e77..0000000000000000000000000000000000000000
--- a/R/GDD.R
+++ /dev/null
@@ -1,31 +0,0 @@
-#'GDD, Growing Degree Days
-#'
-#'
-#'@description This function compute GDD which is defined as the summation of daily differences between daily average temperatures and 10°C between April 1st and October 31st
-#'@details This function requires users to provide the daily average temperature data with the correct period/days for the indicator of interst, and the days must be given in the dimension of 'time_dim'(e.g. 'ftime' in the example).
-#'
-#'@param tx A numeric n-dimensional array containing daily temperature maximum
-#'@param tn A numeric n-dimensional array containing daily temperature minimum
-#'
-#'@return a list of an array with n-1 dimensions
-#'@seealso https://rdrr.io/cran/multiApply/ for multiApply
-#'@import multiApply
-#'
-#'@examples
-#'##Example synthetic data:
-#'tx <- array(rnorm(365*2*3*2*3)*20, c(member = 2, sdate = 3, ftime = 12, lat = 2, lon = 3))
-#'tn <- array(rnorm(365*2*3*2*3)*2, c(member = 2, sdate = 3, ftime = 12, lat = 2, lon = 3))
-#'res <- GDD(tx, tn)
-#'tx1 <- tx[1,1,,1,1]; tn1 <- tn[1,1,,1,1]
-#'temp <- ((tx1 + tn1) / 2) - 10
-#'sum(temp[temp > 0])
-#'res[1,1,1,1]
-#'@export
-GDD <- function(tx, tn, time_dim = 'ftime') {
-  res <- Apply(list(tx, tn), target_dims = list(time_dim, time_dim), fun = .gdd)$output1
-}
-.gdd <- function(tx, tn) {
-  temp <- ((tx + tn) / 2) - 10
-  res <- sum(temp[temp > 0])
-}
-
diff --git a/R/MeanTx.R b/R/MeanTx.R
deleted file mode 100644
index 5cc379d3f552993c9ab199b9980f7b9549061e54..0000000000000000000000000000000000000000
--- a/R/MeanTx.R
+++ /dev/null
@@ -1,26 +0,0 @@
-#'MeanTx (for Sprtx)
-#'
-#'
-#'@description This function compute the average of the daily temperature maximum for a given period. For example, Sprtx is defined as the average of daily temperature maximums between April 1st and May 31st.
-#'@details This function requires users to provide the daily temperature data with the correct period/days for the indicator of interest, and the days must be given in the dimension of 'time_dim'(e.g. 'ftime' in the example).
-#'
-#'@param tx A numeric n-dimensional array containing daily temperature maximum
-#'
-#'@return a list of an array with n-1 dimensions
-#'@seealso https://rdrr.io/cran/multiApply/ for multiApply
-#'@import multiApply
-#'
-#'@examples
-#'##Example synthetic data:
-#'tx <- array(rnorm(365*2*3*2*3) * 5, c(member = 2, sdate = 3, ftime = 365, lat = 2, lon = 3))
-#'res <- MeanTx(tx)
-#'dim(res)
-#'mean(tx[1,1,,1,1])
-#'res[1,1,1,1]
-#'@export
-MeanTx <- function(tx, time_dim = 'ftime') {
-  res <-  Apply(list(tx), target_dims = time_dim, fun = .meantx)$output1
-}
-.meantx <- function(tx) {
-  ta <- mean(tx, na.rm = TRUE)
-}
diff --git a/R/MeanTxTn.R b/R/MeanTxTn.R
deleted file mode 100644
index 2715784491db4336e07ffc7bf17aeec8855b30fd..0000000000000000000000000000000000000000
--- a/R/MeanTxTn.R
+++ /dev/null
@@ -1,29 +0,0 @@
-#'MeanTxTn (GST, Growing Season Temperature)
-#'
-#'
-#'@description This function compute the average of the daily average temperatures for a given period. The daily average temperature here is computed by (daily temperature maximum +  daily temperature minimum) / 2.
-#'@details This function requires users to provide the daily precipitation data with the correct period/days for the indicator of interst, and the days must be given in the dimension of 'time_dim'(e.g. 'ftime' in the example).
-#'
-#'@param tx A numeric n-dimensional array containing daily temperature maximum
-#'@param tn A numeric n-dimensional array containing daily temperature minimum
-#'
-#'@return a list of an array with n-1 dimensions
-#'@seealso https://rdrr.io/cran/multiApply/ for multiApply
-#'@import multiApply
-#'
-#'@examples
-#'##Example synthetic data:
-#'tx <- array(rnorm(365*2*3*2*3) * 5, c(member = 2, sdate = 3, ftime = 365, lat = 2, lon = 3))
-#'tn <- array(rnorm(365*2*3*2*3) * 0.8, c(member = 2, sdate = 3, ftime = 365, lat = 2, lon = 3))
-#'res <- MeanTxTn(tx, tn)
-#'dim(res)
-#'mean((tx[1,1,,1,1]+tn[1,1,,1,1])/2)
-#'res[1,1,1,1]
-#'@export
-MeanTxTn <- function(tx, tn, time_dim = 'ftime') {
-  res <- Apply(list(tx, tn), target_dims = list(time_dim, time_dim), fun = .meantxtn)$output1
-}
-.meantxtn <- function(tx, tn) {
-  ta <- mean((tx + tn) / 2, na.rm = TRUE)
-}
-
diff --git a/R/QThreshold.R b/R/QThreshold.R
index 7d28b2a5a210be81f6a306054aad9efcd42228b7..f89f4a74bff882c0f1ffe841e09e56c2eb12d3d4 100644
--- a/R/QThreshold.R
+++ b/R/QThreshold.R
@@ -1,10 +1,11 @@
 #'Daily thresholds based on values for n-dimensional arrays
 #'
-#'Function to compute the corresponding quantile for a given value as required on function SU for the "Percentile_StressT"
+#'@description Function to compute the corresponding quantile for a given value as required on function SU for the "Percentile_StressT"
 #'@param data A multidimensional array with named dimensions. 
 #'@param threshold How many types of thresholds? When Experiments with several members are provided, threshold is not used. When Observations are provided it can be a single scalar. Should a user be able to provide an array?
-#'@param ncores
-#'@importFrom s2dv Reorder
+#'@param ncores integer indicating the number of cores to use in parallel computation
+#'@return description missed
+#'
 #'@import multiApply
 #'@examples
 #'tasmax <- array(rnorm(15*4, mean = 26, sd = 3),
@@ -18,9 +19,27 @@
 QThreshold <- function(data, threshold = NULL, memb_dim = 'member',
                        sdate_dim = 'sdate', ncores = NULL) {
   dims <- dim(data)
+  if (is.null(dims)) {
+    stop("Parameter 'data' must have dimension names.")
+  }
+  if (!(sdate_dim %in% names(dims))) {
+    stop("Parameter 'data' must have at least one dimension matching the name ",
+         "provided in parameter 'sdate_dim'.")
+  }
+  if (is.null(threshold)) {
+    if (length(dims) == 1) {
+      stop("Parameter 'data' has one dimension and parameter 'threshold'",
+           " cannot be not NULL.")
+    } else {
+      if (!(memb_dim %in% names(dims))) {
+        stop("Neither parameter 'threshold' provided or parameter 'memb_dim' ",
+           "matching names in parameter 'data'.")
+      }
+    }
+  }
   if (memb_dim %in% names(dims) && dims[which(names(dims) == memb_dim)] > 1) {
     res <- Apply(list(data), target_dims = c(memb_dim, sdate_dim),
-                 fun = .qthreshold, ncores = ncores)
+                 fun = .qthreshold, threshold = threshold, ncores = ncores)
   } else {
      res <- Apply(list(data), target_dims = sdate_dim,
                  fun = .qthreshold, threshold = threshold, ncores = ncores) 
@@ -46,14 +65,21 @@ QThreshold <- function(data, threshold = NULL, memb_dim = 'member',
   # order dims: member, sdate (done by Apply)
   dims <- dim(data)
   # For observations (no 'member' involving)
-  if (length(dims) == 1 && any(names(dims) != memb_dim)) {
+  if (length(dims) == 1 && !is.null(threshold)) {
     qres <- unlist(lapply(1:dims, function(x) { 
                              ecdf(data[-x])(threshold)}))
     dim(qres) <- c(dims)
-  }
+  } else if (length(dims) == 2 && !is.null(threshold)) {
+    qres <- unlist(
+    lapply(1:(dim(data)[1]), function(x) { # dim 1: member
+              lapply(1:(dim(data)[2]), function(y) { # dim 2: sdate
+                       ecdf(as.vector(data[,-y]))(threshold)
+                       })
+          }))
+    dim(qres) <- c(dim(data)[2], dim(data)[1])
+    qres <- s2dv::Reorder(qres, c(2, 1))
+  } else if (length(dims) == 2 && is.null(threshold)) {
   # For predictions (all 'member' used)
-
-  if (length(dims) == 2 && memb_dim %in% names(dims)) {
     qres <- unlist(
     lapply(1:(dim(data)[1]), function(x) { # dim 1: member
               lapply(1:(dim(data)[2]), function(y) { # dim 2: sdate
diff --git a/R/SU.R b/R/SU.R
deleted file mode 100644
index f26b8f037a104f23b54b237eb6238d9b6997eff4..0000000000000000000000000000000000000000
--- a/R/SU.R
+++ /dev/null
@@ -1,65 +0,0 @@
-#' Index SU on multidimensional arrays
-#' SU35, SU36 and SU40 (Number of Heat Stress Days above a certain threshold) 
-#' SU is the total count of days when the daily maximum temperatures exceed a certain threshold.
-#'@param data a multidimensional array with named dimensions
-#'@param threshold a single scalar, vector or multidimensional array with named dimensions
-#'@param time_dim a character string indicating the name of the temporal dimension to compute SU index
-#'@param ncores a integer indicating the number of cores to use in parallel computation
-#'@details list of indices (threshold and period considered)
-#' SU35: 35C, seven months
-#' SU36: 36C, June 21 - September 21
-#' SU40: 40C, June 21 - September 21
-#' Spr32: 32C, April 21 - June 21
-###
-#'@import multiApply
-#'@examples
-#'tasmaxexp <- array(rnorm(120),
-#'                   c(member = 3, sdate = 4, ftime = 5, lat = 2))
-#'threshold_exp <- QThreshold(tasmaxexp)
-#'tasmaxobs <- array(rnorm(40),
-#'                   c(member = 1, sdate = 4, ftime = 5, lat = 2))
-#'threshold_obs <- QThreshold(tasmaxobs, 26)
-#'res <- SU(tasmaxexp, threshold_exp)
-#'res <- SU(tasmaxobs, threshold_obs)
-#'res <- SU(threshold_exp, threshold_obs)
-#'dim(res$su_exp)
-#'dim(res$su_obs)
-#'@export
-SU <- function(data, threshold, time_dim = 'ftime', ncores = NULL) {
-  #Checks missed
-  if (is.null(dim(data))) {
-    dim(data) <- c(length(data))
-    names(dim(data)) <- time_dim
-  }
-  if (is.null(dim(threshold)) && length(threshold) > 1) {
-    dim(threshold) <- c(length(threshold))
-    names(dim(threshold)) <- time_dim 
-  }
-  if (!is.character(time_dim)) {
-    stop("Parameter 'time_dim' must be a character string indicating the",
- " name of the dimension to compute the index.")
-  }
-  if (!is.null(ncores)) {
-    if (!is.numeric(ncores)) {
-      ncores <- 1
-      warning("Parameter 'ncores' must be numeric. Execution using 1 core.")
-    }
-    if (length(ncores) > 1) {
-      ncores <- ncores[1]
-      warning("Parameter 'ncores' has length greater than one and only",
-              " the first element will be used.")
-    }
-  }
-  threshold <- drop(threshold) 
-  if (!is.null(dim(threshold)) && !all(dim(threshold) == 1)) {
-    data <- Apply(list(data, threshold), target_dims = list(time_dim, time_dim),
-                 fun = .su, ncores = ncores)$output1
-  } else {
-    data <- Apply(list(data), target_dims = list(time_dim),
-                 fun = .su, threshold = threshold, ncores = ncores)$output1
-  }
-  return(data)
-}
-.su <- function(data, threshold) {
-  result <- sum(data > threshold)
-}
diff --git a/R/TP.R b/R/TP.R
deleted file mode 100644
index 3b89409cf312cb3b769eaa4a0db1917890ef19c3..0000000000000000000000000000000000000000
--- a/R/TP.R
+++ /dev/null
@@ -1,27 +0,0 @@
-#'Total Precipitation (e.g. Spring Rain, SprR / Harvest Rain, HarvestR)
-#'
-#'@description This function computes sum of daily precipitation during a certain period. For example, Spring Rain takes the days from April 21st to June 21st and Harvest Rain takes the days from August 21st to October 21st
-#'
-#'@details This function requires users to provide the daily precipitation data with the correct period/days for the indicator of interst, and the days must be given in the dimension of 'time_dim'(e.g. 'ftime' in the example).
-#'
-#'@param data A numeric n-dimensional array containing daily precipitation
-#'
-#'@return a list of an array with n-1 dimensions
-#'@seealso https://rdrr.io/cran/multiApply/ for multiApply
-#'@import multiApply
-#'
-#'@examples
-#'##Example synthetic data:
-#'prlrexp <- array(rnorm(365*2*3*2), c(member = 2, sdate = 3, ftime = 365, lat = 2))
-#'res <- TP(prlrexp)
-#'dim(prlrexp); dim(res)
-#'sum(prlrexp[1,1,,1]); res[1,1,1]
-#'@export
-TP <- function(data, time_dim = 'ftime') {
-  res <- Apply(list(data), target_dims = time_dim, fun = .tp)$output1
-}
-
-.tp <- function(prlr) {
-  res <- sum(prlr)
-}
-
diff --git a/man/BioIndicators.Rd b/man/BioIndicators.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d47af5167361c5901026893f27d5cb1acb254fef
--- /dev/null
+++ b/man/BioIndicators.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/BioIndicators.R
+\name{BioIndicators}
+\alias{BioIndicators}
+\title{Bioclimate indicators}
+\usage{
+BioIndicators(data, index, time_dim = "ftime", na.rm = TRUE,
+  ncores = NULL)
+}
+\arguments{
+\item{data}{a n-dimensional array or a list of n-dimensional arrays with named dimensions containing daily data depending on each specific index:
+ \item\code{GDD}{data must be a list of two n-dimensional arrays for maximum and minimum daily temperatures in this order}
+ \item\code{GST}{data must be a list of two n-dimensional arrays for maximum and minimum daily temperatures in this order}
+ \item\code{SPRTX}{data must be a n-dimensional array for maximum daily temperatures}
+ \item\code{TP}{data must be a n-dimensional array fir daily precipitation.
+ \item\code{SU}{data must be a list of one n-dimensinoal array containing temperature data (usually daily maximum or minimum) and a scalar, vector or n-dimensional array containing the threshold}
+}}
+
+\item{index}{a character string indicating the index to be computed among \code{GDD}, \code{GST}, \code{SPRTX}, \code{TP} and \code{SU}.}
+
+\item{time_dim}{a character string indicting the name of the temporal dimension.}
+
+\item{na.rm}{a logical value indicting if missing values should be removed.}
+
+\item{ncores}{an integer indicating the number of cores to use in parallel computation.}
+}
+\description{
+Different bioclimate indicators computation for agriculture:
+\itemize{
+ \item\code{GDD}{Growing Degree Days is defined as the summation of daily differences between daily average temperatures and 10°C between April 1st and October 31st}
+ \item\code{GST}{Growing Season Temperature is the average of the daily average temperatures for a given period. The daily average temperature here is computed by (daily temperature maximum +  daily temperature minimum) / 2.}
+ \item\code{SPRTX}{Sprtx is defined as the average of daily temperature maximums between April 1st and May 31st.}
+ \item\code{TP}{Total Precipitation sum of daily precipitation during a certain period. It allows to compute the specific indices \code{SprR} (Spring Rain takes the days from April 21st to June 21st) and \code{HarvestR} (Harvest Rain takes the days from August 21st to October 21st).}
+ \item\code{SU}{Total count of days when daily maximum temperatures exceed a threshold. It allows to compute specific indices \code{SU36} (when 36°C are exceeded between June 21st and September 21st) and \code{SU40} (when 40°C are exceeded between June 21st and September 21st).}
+}
+The user must provide the data for the proper variable and period depending on the desired index.
+}
+\details{
+Requires extra definition?
+}
+\examples{
+data <- array(1:10, c(ftime = 5, x = 2))
+BioIndicators(data, index = 'TP') 
+}
+\seealso{
+\code{QThreshold}
+}
+
diff --git a/man/QThreshold.Rd b/man/QThreshold.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c6134a0d357a077cb593a0a2d43fa89263464ca3
--- /dev/null
+++ b/man/QThreshold.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/QThreshold.R
+\name{QThreshold}
+\alias{QThreshold}
+\title{Daily thresholds based on values for n-dimensional arrays}
+\usage{
+QThreshold(data, threshold = NULL, memb_dim = "member",
+  sdate_dim = "sdate", ncores = NULL)
+}
+\arguments{
+\item{data}{A multidimensional array with named dimensions.}
+
+\item{threshold}{How many types of thresholds? When Experiments with several members are provided, threshold is not used. When Observations are provided it can be a single scalar. Should a user be able to provide an array?}
+
+\item{ncores}{integer indicating the number of cores to use in parallel computation}
+}
+\value{
+description missed
+}
+\description{
+Function to compute the corresponding quantile for a given value as required on function SU for the "Percentile_StressT"
+}
+\examples{
+tasmax <- array(rnorm(15*4, mean = 26, sd = 3),
+     c(dataset = 1, member = 1, sdate = 5, ftime = 3, lon = 2, lat = 2))
+res <- QThreshold(tasmax, 26)
+tasmax <- array(rnorm(15*4*2, mean = 26, sd = 3),
+     c(dataset = 1, member = 2, sdate = 5, ftime = 3, lon = 2, lat = 2))
+res <- QThreshold(tasmax)
+}
+
diff --git a/tests/testthat/test-QThreshold.R b/tests/testthat/test-QThreshold.R
new file mode 100644
index 0000000000000000000000000000000000000000..bfa474968c979fe18dc15325f9fb45c93e892628
--- /dev/null
+++ b/tests/testthat/test-QThreshold.R
@@ -0,0 +1,27 @@
+context("Generic tests")
+test_that("Sanity checks", {
+  data <- 1:10
+  expect_error(QThreshold(data),
+    "Parameter 'data' must have dimension names.")
+  dim(data) <- c(x = 10)
+  expect_error(QThreshold(data), 
+      paste("Parameter 'data' must have at least one dimension",
+            "matching the name provided in parameter 'sdate_dim'."))
+  expect_error(QThreshold(data, sdate = 'x'), 
+      paste("Parameter 'data' has one dimension and parameter",
+       "'threshold' cannot be not NULL."))
+  expect_equal(round(QThreshold(data, threshold = 2, sdate_dim = 'x'),2),
+               array(c(rep(0.11, 2), rep(0.22, 8)), c(x = 10)))
+  dim(data) <- c(x = 5, ftime = 2)
+  expect_error(QThreshold(data), 
+               paste("Parameter 'data' must have at least one dimension",
+            "matching the name provided in parameter 'sdate_dim'."))
+  expect_error(QThreshold(data, sdate_dim = 'ftime'),
+              paste("Neither parameter 'threshold' provided or parameter",
+                    "'memb_dim' matching names in parameter 'data'."))
+  expect_equal(QThreshold(data, sdate_dim = 'ftime', memb_dim = 'x'),
+               array(c(rep(0,5), rep(1,5)), c(x = 5, ftime = 2)))
+  expect_equal(QThreshold(data, sdate_dim = 'ftime', threshold = 2),
+               array(c(0, 1, 0, 1, rep(0,6)), c(ftime = 2, x = 5)))
+          QThreshold(data, sdate_dim = 'ftime', threshold = 2, memb_dim = 'x') 
+}) 
diff --git a/tests/testthat/test-SU.R b/tests/testthat/test-SU.R
index eebfc122e7e43d43ee069cafc355fd61a48df9de..76dc6cb833ec802db16173b9fe6ea3eade399d8a 100644
--- a/tests/testthat/test-SU.R
+++ b/tests/testthat/test-SU.R
@@ -2,20 +2,23 @@ context("Generic tests")
 test_that("Sanity checks", {
   data <- 1:20
   threshold <- 10
-  expect_equal(SU(data, threshold), 10)
+  expect_equal(BioIndicators(list(data, threshold), index = 'SU'), 10)
   data <- array(1:40, c(x = 2, ftime = 20))
-  expect_equal(SU(data, threshold), array(c(15, 15), c(x = 2)))
+  expect_equal(BioIndicators(list(data, threshold), index = 'SU'), 
+               array(c(15, 15), c(x = 2)))
   dim(threshold) <- c(member = 1, ftime = 1)
-  expect_equal(SU(data, threshold), array(c(15, 15), c(x = 2)))
-  expect_equal(SU(data, threshold, time_dim = 'x'),
+  expect_equal(BioIndicators(list(data, threshold), index = 'SU'), 
+               array(c(15, 15), c(x = 2)))
+  expect_equal(BioIndicators(list(data, threshold), index = 'SU', time_dim = 'x'),
                array(c(rep(0,5), rep(2,15)), c(ftime = 20)))
-  expect_warning(SU(data, threshold, time_dim = 'x', ncores = 'Z'),
+  expect_warning(BioIndicators(list(data, threshold), index = 'SU',
+                               time_dim = 'x', ncores = 'Z'),
          "Parameter 'ncores' must be numeric. Execution using 1 core.")
-  expect_warning(SU(data, threshold, time_dim = 'x', ncores = c(1,2)),
+  expect_warning(BioIndicators(list(data, threshold), index = 'SU',
+                              time_dim = 'x', ncores = c(1,2)),
          paste("Parameter 'ncores' has length greater than one and only",
                "the first element will be used."))
-         
-  expect_error(SU(data, threshold, time_dim = 3),
+  expect_error(BioIndicators(list(data, threshold), index = 'SU', time_dim = 3),
          paste("Parameter 'time_dim' must be a character string",
                 "indicating the name of the dimension to compute the index."))