diff --git a/NAMESPACE b/NAMESPACE
index 0a16d4d2dfb5888adf0007be1556694a7f87b899..d92da379f48f9ca97e2c65c7d476b85c73a73a60 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,7 +6,10 @@ export(CST_AbsToProbs)
 export(CST_AccumulationExceedingThreshold)
 export(CST_MergeRefToExp)
 export(CST_PeriodAccumulation)
+export(CST_PeriodMax)
 export(CST_PeriodMean)
+export(CST_PeriodMin)
+export(CST_PeriodVariance)
 export(CST_QThreshold)
 export(CST_SelectPeriodOnData)
 export(CST_Threshold)
@@ -16,7 +19,10 @@ export(CST_WindCapacityFactor)
 export(CST_WindPowerDensity)
 export(MergeRefToExp)
 export(PeriodAccumulation)
+export(PeriodMax)
 export(PeriodMean)
+export(PeriodMin)
+export(PeriodVariance)
 export(QThreshold)
 export(SelectPeriodOnData)
 export(SelectPeriodOnDates)
diff --git a/R/PeriodMax.R b/R/PeriodMax.R
new file mode 100644
index 0000000000000000000000000000000000000000..0806c51d1dd4b6a9b45131a5a750fd580aae7666
--- /dev/null
+++ b/R/PeriodMax.R
@@ -0,0 +1,214 @@
+#'Period Max on 's2dv_cube' objects
+#'
+#'Period Max computes the maximum (max) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO5}{(Providing temperature data) Max Temperature of Warmest 
+#'                   Month. The maximum monthly temperature occurrence over a 
+#'                   given year (time-series) or averaged span of years 
+#'                   (normal).}
+#'  \item\code{BIO13}{(Providing precipitation data) Precipitation of Wettest 
+#'                    Month. This index identifies the total precipitation 
+#'                    that prevails during the wettest month.}
+#'}
+#'
+#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
+#'  package CSTools.
+#'@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 
+#'  \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 
+#'  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.
+#'@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} with dimensions of the input parameter 'data' except the 
+#'dimension where the max has been computed (specified with 'time_dim'). A new 
+#'element called 'time_bounds' will be added into the 'attrs' element in the 
+#''s2dv_cube' object. It consists of a list containing two elements, the start 
+#'and end dates of the aggregated period with the same dimensions of 'Dates' 
+#'element.
+#'
+#'@examples
+#'exp <- NULL
+#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'exp$attrs$Dates <- Dates
+#'class(exp) <- 's2dv_cube'
+#'
+#'res <- CST_PeriodMax(exp, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@importFrom ClimProjDiags Subset
+#'@export
+CST_PeriodMax <- function(data, start = NULL, end = NULL,
+                           time_dim = 'time', na.rm = FALSE,
+                           ncores = NULL) {
+  # Check 's2dv_cube'
+  if (!inherits(data, 's2dv_cube')) {
+    stop("Parameter 'data' must be of the class 's2dv_cube'.")
+  }
+
+  # Dates subset
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dim(data$attrs$Dates))) {
+      warning("Dimensions in 'data' element 'attrs$Dates' are missed and ",
+              "all data would be used.")
+      start <- NULL
+      end <- NULL
+    }
+  }
+
+  Dates <- data$attrs$Dates
+  total <- PeriodMax(data = data$data, dates = Dates, start = start, end = end,
+                      time_dim = time_dim, na.rm = na.rm, ncores = ncores)
+  
+  data$data <- total
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
+  }
+  return(data)
+}
+
+#'Period max on multidimensional array objects
+#'
+#'Period max computes the average (max) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO5}{(Providing temperature data) Max Temperature of Warmest 
+#'                   Month. The maximum monthly temperature occurrence over a 
+#'                   given year (time-series) or averaged span of years 
+#'                   (normal).}
+#'  \item\code{BIO13}{(Providing precipitation data) Precipitation of Wettest 
+#'                    Month. This index identifies the total precipitation 
+#'                    that prevails during the wettest month.}
+#'}
+#'
+#'@param data A multidimensional array with named dimensions.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
+#'@param start An optional parameter to defined the initial date of the period 
+#'  to select from the data by providing a list of two elements: the initial 
+#'  date of the period and the initial month of the period. By default it is set
+#'  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 
+#'  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.
+#'@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}.
+#'
+#'@examples
+#'data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'res <- PeriodMax(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@export
+PeriodMax <- function(data, dates = NULL, start = NULL, end = NULL,
+                      time_dim = 'time', na.rm = FALSE, ncores = NULL) {
+  # Initial checks
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be numeric.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!is.array(data)) {
+    dim(data) <- length(data)
+    names(dim(data)) <- time_dim
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
+      if (!any(c(is.list(start), is.list(end)))) {
+      stop("Parameter 'start' and 'end' must be lists indicating the ",
+            "day and the month of the period start and end.")
+      }
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
+    }
+  }
+  total <- Apply(list(data), target_dims = time_dim, fun = max,
+                 na.rm = na.rm, ncores = ncores)$output1
+  return(total)
+}
+
+
diff --git a/R/PeriodMean.R b/R/PeriodMean.R
index db6a78f273965f19ac056e3b4f5d577e7b239963..7066fdd270bc0edf7a2cf47c9bd0489b3b72342e 100644
--- a/R/PeriodMean.R
+++ b/R/PeriodMean.R
@@ -22,7 +22,7 @@
 #'  the period and the final month of the period. By default it is set to NULL 
 #'  and the indicator is computed using all the data provided in \code{data}.
 #'@param time_dim A character string indicating the name of the dimension to 
-#'  compute the indicator. By default, it is set to 'ftime'. More than one 
+#'  compute the indicator. By default, it is set to 'time'. More than one 
 #'  dimension name matching the dimensions provided in the object 
 #'  \code{data$data} can be specified.
 #'@param na.rm A logical value indicating whether to ignore NA values (TRUE) or 
@@ -41,7 +41,7 @@
 #'
 #'@examples
 #'exp <- NULL
-#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, ftime = 3))
+#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
 #'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
 #'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
 #'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
@@ -50,7 +50,7 @@
 #'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
 #'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
 #'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
-#'dim(Dates) <- c(sdate = 4, ftime = 3)
+#'dim(Dates) <- c(sdate = 4, time = 3)
 #'exp$attrs$Dates <- Dates
 #'class(exp) <- 's2dv_cube'
 #'
@@ -60,12 +60,13 @@
 #'@importFrom ClimProjDiags Subset
 #'@export
 CST_PeriodMean <- function(data, start = NULL, end = NULL,
-                           time_dim = 'ftime', na.rm = FALSE,
+                           time_dim = 'time', na.rm = FALSE,
                            ncores = NULL) {
-   # Check 's2dv_cube'
+  # Check 's2dv_cube'
   if (!inherits(data, 's2dv_cube')) {
     stop("Parameter 'data' must be of the class 's2dv_cube'.")
   }
+
   # Dates subset
   if (!is.null(start) && !is.null(end)) {
     if (is.null(dim(data$attrs$Dates))) {
@@ -163,16 +164,24 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
 #'@export
 PeriodMean <- function(data, dates = NULL, start = NULL, end = NULL,
                        time_dim = 'time', na.rm = FALSE, ncores = NULL) {
-
+  # Initial checks
+  ## data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
   if (!is.numeric(data)) {
     stop("Parameter 'data' must be numeric.")
   }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
   if (!is.array(data)) {
     dim(data) <- length(data)
-    names(data) <- time_dim
+    names(dim(data)) <- time_dim
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
   }
 
   if (!is.null(start) && !is.null(end)) {
diff --git a/R/PeriodMin.R b/R/PeriodMin.R
new file mode 100644
index 0000000000000000000000000000000000000000..842e2e8077673e973f11f1bf340efc03bf83f569
--- /dev/null
+++ b/R/PeriodMin.R
@@ -0,0 +1,214 @@
+#'Period Min on 's2dv_cube' objects
+#'
+#'Period Min computes the average (min) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO6}{(Providing temperature data) Min Temperature of Coldest  
+#'                   Month. The minimum monthly temperature occurrence over a 
+#'                   given year (time-series) or averaged span of years 
+#'                   (normal).}
+#'  \item\code{BIO14}{(Providing precipitation data) Precipitation of Driest  
+#'                    Month. This index identifies the total precipitation 
+#'                    that prevails during the driest month.}
+#'}
+#'
+#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
+#'  package CSTools.
+#'@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 
+#'  \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 
+#'  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.
+#'@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} with dimensions of the input parameter 'data' except the 
+#'dimension where the min has been computed (specified with 'time_dim'). A new 
+#'element called 'time_bounds' will be added into the 'attrs' element in the 
+#''s2dv_cube' object. It consists of a list containing two elements, the start 
+#'and end dates of the aggregated period with the same dimensions of 'Dates' 
+#'element.
+#'
+#'@examples
+#'exp <- NULL
+#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'exp$attrs$Dates <- Dates
+#'class(exp) <- 's2dv_cube'
+#'
+#'res <- CST_PeriodMin(exp, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@importFrom ClimProjDiags Subset
+#'@export
+CST_PeriodMin <- function(data, start = NULL, end = NULL,
+                           time_dim = 'time', na.rm = FALSE,
+                           ncores = NULL) {
+  # Check 's2dv_cube'
+  if (!inherits(data, 's2dv_cube')) {
+    stop("Parameter 'data' must be of the class 's2dv_cube'.")
+  }
+
+  # Dates subset
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dim(data$attrs$Dates))) {
+      warning("Dimensions in 'data' element 'attrs$Dates' are missed and ",
+              "all data would be used.")
+      start <- NULL
+      end <- NULL
+    }
+  }
+
+  Dates <- data$attrs$Dates
+  total <- PeriodMin(data = data$data, dates = Dates, start = start, end = end,
+                      time_dim = time_dim, na.rm = na.rm, ncores = ncores)
+  
+  data$data <- total
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
+  }
+  return(data)
+}
+
+#'Period Min on multidimensional array objects
+#'
+#'Period Min computes the average (min) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO6}{(Providing temperature data) Min Temperature of Coldest  
+#'                   Month. The minimum monthly temperature occurrence over a 
+#'                   given year (time-series) or averaged span of years 
+#'                   (normal).}
+#'  \item\code{BIO14}{(Providing precipitation data) Precipitation of Driest  
+#'                    Month. This index identifies the total precipitation 
+#'                    that prevails during the driest month.}
+#'}
+#'
+#'@param data A multidimensional array with named dimensions.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
+#'@param start An optional parameter to defined the initial date of the period 
+#'  to select from the data by providing a list of two elements: the initial 
+#'  date of the period and the initial month of the period. By default it is set
+#'  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 
+#'  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.
+#'@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}.
+#'
+#'@examples
+#'data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'res <- PeriodMin(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@export
+PeriodMin <- function(data, dates = NULL, start = NULL, end = NULL,
+                      time_dim = 'time', na.rm = FALSE, ncores = NULL) {
+  # Initial checks
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be numeric.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!is.array(data)) {
+    dim(data) <- length(data)
+    names(dim(data)) <- time_dim
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
+      if (!any(c(is.list(start), is.list(end)))) {
+      stop("Parameter 'start' and 'end' must be lists indicating the ",
+            "day and the month of the period start and end.")
+      }
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
+    }
+  }
+  total <- Apply(list(data), target_dims = time_dim, fun = min,
+                 na.rm = na.rm, ncores = ncores)$output1
+  return(total)
+}
+
+
diff --git a/R/PeriodVariance.R b/R/PeriodVariance.R
new file mode 100644
index 0000000000000000000000000000000000000000..b702981bbccbf1f447c5bb19b2bc2f02d1d0c078
--- /dev/null
+++ b/R/PeriodVariance.R
@@ -0,0 +1,228 @@
+#'Period Variance on 's2dv_cube' objects
+#'
+#'Period Variance computes the average (var) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO4}{(Providing temperature data) Temperature Seasonality 
+#'                   (Standard Deviation). The amount of temperature variation  
+#'                   over a given year (or averaged years) based on the standard 
+#'                   deviation (variation) of monthly temperature averages. }
+#'  \item\code{BIO15}{(Providing precipitation data) Precipitation Seasonality 
+#'                    (CV). This is a measure of the variation in 
+#'                    monthly precipitation totals over the course of the year. 
+#'                    This index is the ratio of the standard deviation of the 
+#'                    monthly total precipitation to the mean monthly total 
+#'                    precipitation (also known as the coefficient of variation) 
+#'                    and is expressed as a percentage}
+#'}
+#'
+#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
+#'  package CSTools.
+#'@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 
+#'  \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 
+#'  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.
+#'@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} with dimensions of the input parameter 'data' except the 
+#'dimension where the var has been computed (specified with 'time_dim'). A new 
+#'element called 'time_bounds' will be added into the 'attrs' element in the 
+#''s2dv_cube' object. It consists of a list containing two elements, the start 
+#'and end dates of the aggregated period with the same dimensions of 'Dates' 
+#'element.
+#'
+#'@examples
+#'exp <- NULL
+#'exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'exp$attrs$Dates <- Dates
+#'class(exp) <- 's2dv_cube'
+#'
+#'res <- CST_PeriodVariance(exp, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@importFrom ClimProjDiags Subset
+#'@export
+CST_PeriodVariance <- function(data, start = NULL, end = NULL,
+                           time_dim = 'time', na.rm = FALSE,
+                           ncores = NULL) {
+  # Check 's2dv_cube'
+  if (!inherits(data, 's2dv_cube')) {
+    stop("Parameter 'data' must be of the class 's2dv_cube'.")
+  }
+
+  # Dates subset
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dim(data$attrs$Dates))) {
+      warning("Dimensions in 'data' element 'attrs$Dates' are missed and ",
+              "all data would be used.")
+      start <- NULL
+      end <- NULL
+    }
+  }
+
+  Dates <- data$attrs$Dates
+  total <- PeriodVariance(data = data$data, dates = Dates, start = start, end = end,
+                          time_dim = time_dim, na.rm = na.rm, ncores = ncores)
+  
+  data$data <- total
+  data$dims <- dim(total)
+
+  if (!is.null(Dates)) {
+    if (!is.null(start) && !is.null(end)) {
+      Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
+                                   time_dim = time_dim, ncores = ncores)
+    }
+    if (is.null(dim(Dates))) {
+      warning("Element 'Dates' has NULL dimensions. They will not be ", 
+              "subset and 'time_bounds' will be missed.")
+      data$attrs$Dates <- Dates
+    } else {
+      # Create time_bounds
+      time_bounds <- NULL
+      time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                                 indices = 1, drop = 'selected')
+      time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
+                                               indices = dim(Dates)[time_dim], 
+                                               drop = 'selected')
+
+      # Add Dates in attrs
+      data$attrs$Dates <- time_bounds$start
+      data$attrs$time_bounds <- time_bounds
+    }
+  }
+  return(data)
+}
+
+#'Period Variance on multidimensional array objects
+#'
+#'Period Variance computes the average (var) of a given variable in a period.
+#'Two bioclimatic indicators can be obtained by using this function:
+#'\itemize{
+#'  \item\code{BIO4}{(Providing temperature data) Temperature Seasonality 
+#'                   (Standard Deviation). The amount of temperature variation  
+#'                   over a given year (or averaged years) based on the standard 
+#'                   deviation (variation) of monthly temperature averages. }
+#'  \item\code{BIO15}{(Providing precipitation data) Precipitation Seasonality 
+#'                    (CV). This is a measure of the variation in 
+#'                    monthly precipitation totals over the course of the year. 
+#'                    This index is the ratio of the standard deviation of the 
+#'                    monthly total precipitation to the mean monthly total 
+#'                    precipitation (also known as the coefficient of variation) 
+#'                    and is expressed as a percentage}
+#'}
+#'
+#'@param data A multidimensional array with named dimensions.
+#'@param dates A multidimensional array of dates with named dimensions matching 
+#'  the temporal dimensions on parameter 'data'. By default it is NULL, to  
+#'  select aperiod this parameter must be provided.
+#'@param start An optional parameter to defined the initial date of the period 
+#'  to select from the data by providing a list of two elements: the initial 
+#'  date of the period and the initial month of the period. By default it is set
+#'  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 
+#'  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.
+#'@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}.
+#'
+#'@examples
+#'data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+#'Dates <- c(seq(as.Date("2000-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2001-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2001-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2002-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2002-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2003-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"),
+#'           seq(as.Date("2003-11-01", "%Y-%m-%d", tz = "UTC"), 
+#'               as.Date("2004-01-01", "%Y-%m-%d", tz = "UTC"), by = "month"))
+#'dim(Dates) <- c(sdate = 4, time = 3)
+#'res <- PeriodVariance(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+#'
+#'@import multiApply
+#'@export
+PeriodVariance <- function(data, dates = NULL, start = NULL, end = NULL,
+                           time_dim = 'time', na.rm = FALSE, ncores = NULL) {
+  # Initial checks
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be numeric.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!is.array(data)) {
+    dim(data) <- length(data)
+    names(dim(data)) <- time_dim
+  }
+  if (!time_dim %in% names(dim(data))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+
+  if (!is.null(start) && !is.null(end)) {
+    if (is.null(dates)) {
+      warning("Parameter 'dates' is NULL and the average of the ",
+              "full data provided in 'data' is computed.")
+    } else {
+      if (!any(c(is.list(start), is.list(end)))) {
+      stop("Parameter 'start' and 'end' must be lists indicating the ",
+            "day and the month of the period start and end.")
+      }
+      if (!is.null(dim(dates))) {
+        data <- SelectPeriodOnData(data = data, dates = dates, start = start, 
+                                   end = end, time_dim = time_dim, 
+                                   ncores = ncores)
+      } else {
+        warning("Parameter 'dates' must have named dimensions if 'start' and ",
+                "'end' are not NULL. All data will be used.")
+      }
+    }
+  }
+  total <- Apply(list(data), target_dims = time_dim, 
+                 fun = .periodvariance,
+                 na.rm = na.rm, ncores = ncores)$output1
+  return(total)
+}
+
+.periodvariance <- function(data, na.rm) {
+  var <- sum((data - mean(data, na.rm = na.rm))^2) / (length(data)-1)
+  return(var)
+}
+
+
diff --git a/man/CST_PeriodMax.Rd b/man/CST_PeriodMax.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..7b01d145482423cf55a375ba9bc9a31fb7f6d4eb
--- /dev/null
+++ b/man/CST_PeriodMax.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodMax.R
+\name{CST_PeriodMax}
+\alias{CST_PeriodMax}
+\title{Period Max on 's2dv_cube' objects}
+\usage{
+CST_PeriodMax(
+  data,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
+package CSTools.}
+
+\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 
+\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 
+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.}
+
+\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{
+An 's2dv_cube' object containing the indicator in the element 
+\code{data} with dimensions of the input parameter 'data' except the 
+dimension where the max has been computed (specified with 'time_dim'). A new 
+element called 'time_bounds' will be added into the 'attrs' element in the 
+'s2dv_cube' object. It consists of a list containing two elements, the start 
+and end dates of the aggregated period with the same dimensions of 'Dates' 
+element.
+}
+\description{
+Period Max computes the maximum (max) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO5}{(Providing temperature data) Max Temperature of Warmest 
+                  Month. The maximum monthly temperature occurrence over a 
+                  given year (time-series) or averaged span of years 
+                  (normal).}
+ \item\code{BIO13}{(Providing precipitation data) Precipitation of Wettest 
+                   Month. This index identifies the total precipitation 
+                   that prevails during the wettest month.}
+}
+}
+\examples{
+exp <- NULL
+exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+exp$attrs$Dates <- Dates
+class(exp) <- 's2dv_cube'
+
+res <- CST_PeriodMax(exp, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/man/CST_PeriodMean.Rd b/man/CST_PeriodMean.Rd
index 0aa4aa359645d70ba850f556586488e8980a431d..ab420665dabdb4fa8f94a40f8c0ed2d65552bda0 100644
--- a/man/CST_PeriodMean.Rd
+++ b/man/CST_PeriodMean.Rd
@@ -8,7 +8,7 @@ CST_PeriodMean(
   data,
   start = NULL,
   end = NULL,
-  time_dim = "ftime",
+  time_dim = "time",
   na.rm = FALSE,
   ncores = NULL
 )
@@ -29,7 +29,7 @@ the period and the final month of the period. By default it is set to NULL
 and the indicator is computed using all the data provided in \code{data}.}
 
 \item{time_dim}{A character string indicating the name of the dimension to 
-compute the indicator. By default, it is set to 'ftime'. More than one 
+compute the indicator. By default, it is set to 'time'. More than one 
 dimension name matching the dimensions provided in the object 
 \code{data$data} can be specified.}
 
@@ -62,7 +62,7 @@ this function:
 }
 \examples{
 exp <- NULL
-exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, ftime = 3))
+exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
 Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
               as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
           seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
@@ -71,7 +71,7 @@ Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"),
               as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
           seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
               as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
-dim(Dates) <- c(sdate = 4, ftime = 3)
+dim(Dates) <- c(sdate = 4, time = 3)
 exp$attrs$Dates <- Dates
 class(exp) <- 's2dv_cube'
 
diff --git a/man/CST_PeriodMin.Rd b/man/CST_PeriodMin.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..939569986386893c495d8789bdbaca904393f1d4
--- /dev/null
+++ b/man/CST_PeriodMin.Rd
@@ -0,0 +1,81 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodMin.R
+\name{CST_PeriodMin}
+\alias{CST_PeriodMin}
+\title{Period Min on 's2dv_cube' objects}
+\usage{
+CST_PeriodMin(
+  data,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
+package CSTools.}
+
+\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 
+\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 
+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.}
+
+\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{
+An 's2dv_cube' object containing the indicator in the element 
+\code{data} with dimensions of the input parameter 'data' except the 
+dimension where the min has been computed (specified with 'time_dim'). A new 
+element called 'time_bounds' will be added into the 'attrs' element in the 
+'s2dv_cube' object. It consists of a list containing two elements, the start 
+and end dates of the aggregated period with the same dimensions of 'Dates' 
+element.
+}
+\description{
+Period Min computes the average (min) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO6}{(Providing temperature data) Min Temperature of Coldest  
+                  Month. The minimum monthly temperature occurrence over a 
+                  given year (time-series) or averaged span of years 
+                  (normal).}
+ \item\code{BIO14}{(Providing precipitation data) Precipitation of Driest  
+                   Month. This index identifies the total precipitation 
+                   that prevails during the driest month.}
+}
+}
+\examples{
+exp <- NULL
+exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+exp$attrs$Dates <- Dates
+class(exp) <- 's2dv_cube'
+
+res <- CST_PeriodMin(exp, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/man/CST_PeriodVariance.Rd b/man/CST_PeriodVariance.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e28bf9551727659ffb53e257ec2c29fe7a8829b2
--- /dev/null
+++ b/man/CST_PeriodVariance.Rd
@@ -0,0 +1,85 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodVariance.R
+\name{CST_PeriodVariance}
+\alias{CST_PeriodVariance}
+\title{Period Variance on 's2dv_cube' objects}
+\usage{
+CST_PeriodVariance(
+  data,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
+package CSTools.}
+
+\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 
+\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 
+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.}
+
+\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{
+An 's2dv_cube' object containing the indicator in the element 
+\code{data} with dimensions of the input parameter 'data' except the 
+dimension where the var has been computed (specified with 'time_dim'). A new 
+element called 'time_bounds' will be added into the 'attrs' element in the 
+'s2dv_cube' object. It consists of a list containing two elements, the start 
+and end dates of the aggregated period with the same dimensions of 'Dates' 
+element.
+}
+\description{
+Period Variance computes the average (var) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO4}{(Providing temperature data) Temperature Seasonality 
+                  (Standard Deviation). The amount of temperature variation  
+                  over a given year (or averaged years) based on the standard 
+                  deviation (variation) of monthly temperature averages. }
+ \item\code{BIO15}{(Providing precipitation data) Precipitation Seasonality 
+                   (CV). This is a measure of the variation in 
+                   monthly precipitation totals over the course of the year. 
+                   This index is the ratio of the standard deviation of the 
+                   monthly total precipitation to the mean monthly total 
+                   precipitation (also known as the coefficient of variation) 
+                   and is expressed as a percentage}
+}
+}
+\examples{
+exp <- NULL
+exp$data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+exp$attrs$Dates <- Dates
+class(exp) <- 's2dv_cube'
+
+res <- CST_PeriodVariance(exp, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/man/PeriodMax.Rd b/man/PeriodMax.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..cb776d24cb9d308153cb1f230cab81ce3831828a
--- /dev/null
+++ b/man/PeriodMax.Rd
@@ -0,0 +1,76 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodMax.R
+\name{PeriodMax}
+\alias{PeriodMax}
+\title{Period max on multidimensional array objects}
+\usage{
+PeriodMax(
+  data,
+  dates = NULL,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A multidimensional array with named dimensions.}
+
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
+
+\item{start}{An optional parameter to defined the initial date of the period 
+to select from the data by providing a list of two elements: the initial 
+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 
+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.}
+
+\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}.
+}
+\description{
+Period max computes the average (max) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO5}{(Providing temperature data) Max Temperature of Warmest 
+                  Month. The maximum monthly temperature occurrence over a 
+                  given year (time-series) or averaged span of years 
+                  (normal).}
+ \item\code{BIO13}{(Providing precipitation data) Precipitation of Wettest 
+                   Month. This index identifies the total precipitation 
+                   that prevails during the wettest month.}
+}
+}
+\examples{
+data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+res <- PeriodMax(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/man/PeriodMin.Rd b/man/PeriodMin.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..154acf0872b04ca1e2de993ceeaa94cf7b36975a
--- /dev/null
+++ b/man/PeriodMin.Rd
@@ -0,0 +1,76 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodMin.R
+\name{PeriodMin}
+\alias{PeriodMin}
+\title{Period Min on multidimensional array objects}
+\usage{
+PeriodMin(
+  data,
+  dates = NULL,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A multidimensional array with named dimensions.}
+
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
+
+\item{start}{An optional parameter to defined the initial date of the period 
+to select from the data by providing a list of two elements: the initial 
+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 
+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.}
+
+\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}.
+}
+\description{
+Period Min computes the average (min) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO6}{(Providing temperature data) Min Temperature of Coldest  
+                  Month. The minimum monthly temperature occurrence over a 
+                  given year (time-series) or averaged span of years 
+                  (normal).}
+ \item\code{BIO14}{(Providing precipitation data) Precipitation of Driest  
+                   Month. This index identifies the total precipitation 
+                   that prevails during the driest month.}
+}
+}
+\examples{
+data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+res <- PeriodMin(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/man/PeriodVariance.Rd b/man/PeriodVariance.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e1e8d7cd6137accd0d8fdfaf704470246b8a272d
--- /dev/null
+++ b/man/PeriodVariance.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodVariance.R
+\name{PeriodVariance}
+\alias{PeriodVariance}
+\title{Period Variance on multidimensional array objects}
+\usage{
+PeriodVariance(
+  data,
+  dates = NULL,
+  start = NULL,
+  end = NULL,
+  time_dim = "time",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A multidimensional array with named dimensions.}
+
+\item{dates}{A multidimensional array of dates with named dimensions matching 
+the temporal dimensions on parameter 'data'. By default it is NULL, to  
+select aperiod this parameter must be provided.}
+
+\item{start}{An optional parameter to defined the initial date of the period 
+to select from the data by providing a list of two elements: the initial 
+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 
+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.}
+
+\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}.
+}
+\description{
+Period Variance computes the average (var) of a given variable in a period.
+Two bioclimatic indicators can be obtained by using this function:
+\itemize{
+ \item\code{BIO4}{(Providing temperature data) Temperature Seasonality 
+                  (Standard Deviation). The amount of temperature variation  
+                  over a given year (or averaged years) based on the standard 
+                  deviation (variation) of monthly temperature averages. }
+ \item\code{BIO15}{(Providing precipitation data) Precipitation Seasonality 
+                   (CV). This is a measure of the variation in 
+                   monthly precipitation totals over the course of the year. 
+                   This index is the ratio of the standard deviation of the 
+                   monthly total precipitation to the mean monthly total 
+                   precipitation (also known as the coefficient of variation) 
+                   and is expressed as a percentage}
+}
+}
+\examples{
+data <- array(rnorm(45), dim = c(member = 7, sdate = 4, time = 3))
+Dates <- c(seq(as.Date("2000-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2001-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2001-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2002-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2002-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2003-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"),
+          seq(as.Date("2003-11-01", "\%Y-\%m-\%d", tz = "UTC"), 
+              as.Date("2004-01-01", "\%Y-\%m-\%d", tz = "UTC"), by = "month"))
+dim(Dates) <- c(sdate = 4, time = 3)
+res <- PeriodVariance(data, dates = Dates, start = list(01, 12), end = list(01, 01))
+
+}
diff --git a/tests/testthat/test-PeriodMax.R b/tests/testthat/test-PeriodMax.R
new file mode 100644
index 0000000000000000000000000000000000000000..97907e6cf86a1c3614fcaeb8597faea8d93b94c8
--- /dev/null
+++ b/tests/testthat/test-PeriodMax.R
@@ -0,0 +1,133 @@
+library(CSTools)
+
+##############################################
+test_that("1. Sanity Checks", {
+  # data
+  expect_error(
+    PeriodMax(data = NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    PeriodMax('x'), 
+    "Parameter 'data' must be numeric."
+  )
+  # time_dim
+  expect_error(
+    PeriodMax(array(1:10), time_dim = 1), 
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodMax(array(1:10, dim = c(ftime = 10)), time_dim = 'time'), 
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  suppressWarnings(
+    expect_equal(
+      PeriodMax(array(1, c(x = 1)), time_dim = 'x'), 
+      1
+    )
+  )
+  expect_error(
+    PeriodMax(1, dates = '2000-01-01', end = 3, start = 4),
+    "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end."
+  )
+  suppressWarnings(
+    expect_equal(
+      PeriodMax(array(1:10, c(time = 10))), 
+      10
+    )
+  )
+  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
+  suppressWarnings(
+    expect_equal(
+      PeriodMax(data),
+      array(c(5, 6, 11, 12, 17, 18, 23, 24), 
+      c(sdate = 2, lon = 4))
+    )
+  )
+  # Test dates warning
+  expect_warning(
+    PeriodMax(array(1:10, c(time = 10)), 
+               dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                           as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and 'end' ", 
+           "are not NULL. All data will be used.")
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    PeriodMax(array(1:10, c(time = 10)), 
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
+})
+
+##############################################
+
+test_that("2. Seasonal", {
+  exp <- NULL
+  exp$data <- array(1:(1 * 3 * 214 * 2),
+                    c(memb = 1, sdate = 3, time = 214, lon = 2))
+  exp$dims <- dim(exp$data)
+  exp$attrs$Dates <- c(seq(as.Date("01-04-2000", format = "%d-%m-%Y"),
+                           as.Date("31-10-2000", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2001", format = "%d-%m-%Y"),
+                           as.Date("31-10-2001", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
+                           as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
+  dim(exp$attrs$Dates) <- c(time = 214, sdate = 3)
+  class(exp) <- 's2dv_cube'
+  output <- exp
+  output$data <- array(c(max(exp$data[1,1,21:82,1]), max(exp$data[1,2,21:82,1]),
+                         max(exp$data[1,3,21:82,1]), max(exp$data[1,1,21:82,2]),
+                         max(exp$data[1,2,21:82,2]), max(exp$data[1,3,21:82,2])),
+                       c(memb = 1, sdate = 3, lon = 2))
+  expect_equal(
+    CST_PeriodMax(exp, start = list(21, 4), end = list(21, 6))$data,
+    output$data
+  )
+})
+
+##############################################
+test_that("3. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodMax(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
+
diff --git a/tests/testthat/test-PeriodMean.R b/tests/testthat/test-PeriodMean.R
index cd9f5fe5ddca8b82c54fbbbc661bf68e0e8490a0..f51f7a17e579070fb2a250ac7118cc84e5893290 100644
--- a/tests/testthat/test-PeriodMean.R
+++ b/tests/testthat/test-PeriodMean.R
@@ -2,6 +2,7 @@ library(CSTools)
 
 ##############################################
 test_that("1. Sanity Checks", {
+  # data
   expect_error(
     PeriodMean('x'), 
     "Parameter 'data' must be numeric."
@@ -16,17 +17,26 @@ test_that("1. Sanity Checks", {
     PeriodMean(data = NULL), 
     "Parameter 'data' cannot be NULL."
   )
+  # time_dim
+  expect_error(
+    PeriodMean(array(1:10), time_dim = 1), 
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodMean(array(1:10, dim = c(ftime = 10)), time_dim = 'time'), 
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
   expect_error(
     PeriodMean(1, dates = '2000-01-01', end = 3, start = 4),
     "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end."
   )
   suppressWarnings(
     expect_equal(
-      PeriodMean(array(1:10, c(ftime = 10))), 
+      PeriodMean(array(1:10, c(time = 10))), 
       5.5
     )
   )
-  data <- array(1:24, c(sdate = 2, ftime = 3, lon = 4))
+  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
   suppressWarnings(
     expect_equal(
       PeriodMean(data),
@@ -36,7 +46,7 @@ test_that("1. Sanity Checks", {
   )
   # Test dates warning
   expect_warning(
-    PeriodMean(array(1:10, c(ftime = 10)), 
+    PeriodMean(array(1:10, c(time = 10)), 
                dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
                            as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
                start = list(05, 02), end = list(05, 09)),
@@ -45,7 +55,7 @@ test_that("1. Sanity Checks", {
   )
   # start and end when dates is not provided
   expect_warning(
-    PeriodMean(array(1:10, c(ftime = 10)), 
+    PeriodMean(array(1:10, c(time = 10)), 
                start = list(05, 02), end = list(05, 09)),
     paste0("Parameter 'dates' is NULL and the average of the full data ", 
            "provided in 'data' is computed.")
@@ -57,7 +67,7 @@ test_that("1. Sanity Checks", {
 test_that("2. Seasonal", {
   exp <- NULL
   exp$data <- array(1:(1 * 3 * 214 * 2),
-                    c(memb = 1, sdate = 3, ftime = 214, lon = 2))
+                    c(memb = 1, sdate = 3, time = 214, lon = 2))
   exp$dims <- dim(exp$data)
   exp$attrs$Dates <- c(seq(as.Date("01-04-2000", format = "%d-%m-%Y"),
                            as.Date("31-10-2000", format = "%d-%m-%Y"), by = 'day'),
@@ -65,7 +75,7 @@ test_that("2. Seasonal", {
                            as.Date("31-10-2001", format = "%d-%m-%Y"), by = 'day'),
                        seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
                            as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
-  dim(exp$attrs$Dates) <- c(ftime = 214, sdate = 3)
+  dim(exp$attrs$Dates) <- c(time = 214, sdate = 3)
   class(exp) <- 's2dv_cube'
   output <- exp
   output$data <- array(c(mean(exp$data[1,1,21:82,1]), mean(exp$data[1,2,21:82,1]),
diff --git a/tests/testthat/test-PeriodMin.R b/tests/testthat/test-PeriodMin.R
new file mode 100644
index 0000000000000000000000000000000000000000..5ed6c1f6e8014964a5963553b08fae9c37b9f2ca
--- /dev/null
+++ b/tests/testthat/test-PeriodMin.R
@@ -0,0 +1,129 @@
+library(CSTools)
+
+##############################################
+test_that("1. Sanity Checks", {
+  # data
+  expect_error(
+    PeriodMin('x'), 
+    "Parameter 'data' must be numeric."
+  )
+  suppressWarnings(
+    expect_equal(
+      PeriodMin(array(1, c(x = 1)), time_dim = 'x'), 
+      1
+    )
+  )
+  expect_error(
+    PeriodMin(data = NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  # time_dim
+  expect_error(
+    PeriodMin(array(1:10), time_dim = 1), 
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodMin(array(1:10, dim = c(ftime = 10)), time_dim = 'time'), 
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  expect_error(
+    PeriodMin(1, dates = '2000-01-01', end = 3, start = 4),
+    "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end."
+  )
+  expect_equal(
+    PeriodMin(array(1:10, c(time = 10))), 
+    1
+  )
+  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
+  expect_equal(
+    PeriodMin(data),
+    array(c(1, 2, 7, 8, 13, 14, 19, 20), 
+    c(sdate = 2, lon = 4))
+  )
+  # Test dates warning
+  expect_warning(
+    PeriodMin(array(1:10, c(time = 10)), 
+               dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                           as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and 'end' ", 
+           "are not NULL. All data will be used.")
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    PeriodMin(array(1:10, c(time = 10)), 
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
+})
+
+##############################################
+
+test_that("2. Seasonal", {
+  exp <- NULL
+  exp$data <- array(1:(1 * 3 * 214 * 2),
+                    c(memb = 1, sdate = 3, time = 214, lon = 2))
+  exp$dims <- dim(exp$data)
+  exp$attrs$Dates <- c(seq(as.Date("01-04-2000", format = "%d-%m-%Y"),
+                           as.Date("31-10-2000", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2001", format = "%d-%m-%Y"),
+                           as.Date("31-10-2001", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
+                           as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
+  dim(exp$attrs$Dates) <- c(time = 214, sdate = 3)
+  class(exp) <- 's2dv_cube'
+  output <- exp
+  output$data <- array(c(min(exp$data[1,1,21:82,1]), min(exp$data[1,2,21:82,1]),
+                         min(exp$data[1,3,21:82,1]), min(exp$data[1,1,21:82,2]),
+                         min(exp$data[1,2,21:82,2]), min(exp$data[1,3,21:82,2])),
+                       c(memb = 1, sdate = 3, lon = 2))
+  expect_equal(
+    CST_PeriodMin(exp, start = list(21, 4), end = list(21, 6))$data,
+    output$data
+  )
+})
+
+##############################################
+test_that("3. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodMin(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
+
diff --git a/tests/testthat/test-PeriodVariance.R b/tests/testthat/test-PeriodVariance.R
new file mode 100644
index 0000000000000000000000000000000000000000..e52564490b827b6e71d03c028605905f13499b92
--- /dev/null
+++ b/tests/testthat/test-PeriodVariance.R
@@ -0,0 +1,128 @@
+library(CSTools)
+
+##############################################
+test_that("1. Sanity Checks", {
+  # data
+  expect_error(
+    PeriodVariance('x'), 
+    "Parameter 'data' must be numeric."
+  )
+  expect_equal(
+    PeriodVariance(array(1:2, c(x = 2)), time_dim = 'x'), 
+    0.5
+  )
+  expect_error(
+    PeriodVariance(data = NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  # time_dim
+  expect_error(
+    PeriodVariance(array(1:10), time_dim = 1), 
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodVariance(array(1:10, dim = c(ftime = 10)), time_dim = 'time'), 
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  expect_error(
+    PeriodVariance(1, dates = '2000-01-01', end = 3, start = 4),
+    "Parameter 'start' and 'end' must be lists indicating the day and the month of the period start and end."
+  )
+  expect_equal(
+    PeriodVariance(array(1:10, c(time = 10))), 
+    9.166667,
+    tolerance = 0.001
+  )
+  data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
+  expect_equal(
+    PeriodVariance(data),
+    array(rep(4, 8), 
+    c(sdate = 2, lon = 4))
+  )
+  # Test dates warning
+  expect_warning(
+    PeriodVariance(array(1:10, c(time = 10)), 
+               dates = seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                           as.Date("10-05-2000", format = "%d-%m-%Y"), by = 'day'),
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' must have named dimensions if 'start' and 'end' ", 
+           "are not NULL. All data will be used.")
+  )
+  # start and end when dates is not provided
+  expect_warning(
+    PeriodVariance(array(1:10, c(time = 10)), 
+               start = list(05, 02), end = list(05, 09)),
+    paste0("Parameter 'dates' is NULL and the average of the full data ", 
+           "provided in 'data' is computed.")
+  )
+})
+
+##############################################
+
+test_that("2. Seasonal", {
+  exp <- NULL
+  exp$data <- array(1:(1 * 3 * 214 * 2),
+                    c(memb = 1, sdate = 3, ftime = 214, lon = 2))
+  exp$dims <- dim(exp$data)
+  exp$attrs$Dates <- c(seq(as.Date("01-04-2000", format = "%d-%m-%Y"),
+                           as.Date("31-10-2000", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2001", format = "%d-%m-%Y"),
+                           as.Date("31-10-2001", format = "%d-%m-%Y"), by = 'day'),
+                       seq(as.Date("01-04-2002", format = "%d-%m-%Y"),
+                           as.Date("31-10-2002", format = "%d-%m-%Y"), by = 'day'))
+  dim(exp$attrs$Dates) <- c(ftime = 214, sdate = 3)
+  class(exp) <- 's2dv_cube'
+  output <- exp
+  output$data <- array(c(var(exp$data[1,1,21:82,1]), var(exp$data[1,2,21:82,1]),
+                         var(exp$data[1,3,21:82,1]), var(exp$data[1,1,21:82,2]),
+                         var(exp$data[1,2,21:82,2]), var(exp$data[1,3,21:82,2])),
+                       c(memb = 1, sdate = 3, lon = 2))
+  expect_equal(
+    CST_PeriodVariance(exp, start = list(21, 4), end = list(21, 6))$data,
+    output$data
+  )
+})
+
+##############################################
+test_that("3. Subset Dates and check time_bounds", {
+  exp <- CSTools::lonlat_prec
+  res <- CST_PeriodVariance(data = CSTools::lonlat_prec, time_dim = 'ftime',
+                        start = list(10, 03), end = list(20, 03))
+  # Check dims
+  expect_equal(
+    dim(res$data), 
+    res$dims
+  )
+  # Check Dates
+  expect_equal(
+    dim(res$data)['sdate'], 
+    dim(res$attrs$Dates)
+  )
+  # Check time_bounds
+  expect_equal(
+    res$attrs$Dates, 
+    res$attrs$time_bounds$start
+  )
+  expect_equal(
+    dim(res$attrs$time_bounds$start), 
+    dim(res$attrs$time_bounds$end)
+  )
+  # Check 'sdates'
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$start) == 3),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$start) == 10),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::month(res$attrs$time_bounds$end) == 03),
+    TRUE
+  )
+  expect_equal(
+    all(lubridate::day(res$attrs$time_bounds$end) == 20),
+    TRUE
+  )
+})
+