diff --git a/NAMESPACE b/NAMESPACE
index 0a16d4d2dfb5888adf0007be1556694a7f87b899..7bed2f618f8cea4a56254e396088581a629f37c7 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,6 +7,8 @@ export(CST_AccumulationExceedingThreshold)
 export(CST_MergeRefToExp)
 export(CST_PeriodAccumulation)
 export(CST_PeriodMean)
+export(CST_PeriodPET)
+export(CST_PeriodSPEI)
 export(CST_QThreshold)
 export(CST_SelectPeriodOnData)
 export(CST_Threshold)
@@ -17,6 +19,8 @@ export(CST_WindPowerDensity)
 export(MergeRefToExp)
 export(PeriodAccumulation)
 export(PeriodMean)
+export(PeriodPET)
+export(PeriodSPEI)
 export(QThreshold)
 export(SelectPeriodOnData)
 export(SelectPeriodOnDates)
@@ -25,7 +29,15 @@ export(TotalSpellTimeExceedingThreshold)
 export(TotalTimeExceedingThreshold)
 export(WindCapacityFactor)
 export(WindPowerDensity)
+import(CSTools)
+import(ClimProjDiags)
+import(SPEI)
+import(TLMoments)
+import(lmom)
+import(lmomco)
+import(lubridate)
 import(multiApply)
+import(zoo)
 importFrom(ClimProjDiags,Subset)
 importFrom(stats,approxfun)
 importFrom(stats,ecdf)
diff --git a/R/PeriodPET.R b/R/PeriodPET.R
new file mode 100644
index 0000000000000000000000000000000000000000..f365d39b444bca4aca8108e7b4f33523c896ba4c
--- /dev/null
+++ b/R/PeriodPET.R
@@ -0,0 +1,409 @@
+#'Compute the Potential Evapotranspiration
+#' 
+#'Compute the Potential evapotranspiration (PET) that is the amount of 
+#'evaporation and transpiration that would occur if a sufficient water source 
+#'were available. This function calculate PET according to the Thornthwaite, 
+#'Hargreaves or Hargreaves-modified equations. 
+#' 
+#'This function is build to work be compatible with other tools in 
+#'that work with 's2dv_cube' object class. The input data must be this object 
+#'class. If you don't work with 's2dv_cube', see PeriodPET. For more information 
+#'on the SPEI calculation, see functions CST_PeriodStandardization and 
+#'CST_PeriodAccumulation.
+#'
+#'@param data A named list with the needed \code{s2dv_cube} objects containing 
+#'  the seasonal forecast experiment in the data element for each variable. 
+#'  Specific variables are needed for each method used in computing the 
+#'  Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+#'  variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for 
+#'  'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' 
+#'  'tmean' is required. The units for temperature variables
+#'  ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+#'  precipitation ('pr') need to be in mm/month. 
+#'@param dates An array of temporal dimensions containing the Dates of 
+#'  'data'. It must be of class 'Date' or 'POSIXct'.
+#'@param lat A numeric vector containing the latitude values of 'data'.
+#'@param pet_method A character string indicating the method used to compute 
+#'  the potential evapotranspiration. The accepted methods are:
+#'  'hargreaves' and 'hargreaves_modified', that require the data to have 
+#'  variables tmin and tmax; and 'thornthwaite', that requires variable 
+#'  'tmean'. 
+#'@param time_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'syear'. 
+#'@param leadtime_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'time'. 
+#'@param lat_dim A character string indicating the name of the latitudinal 
+#'  dimension. By default it is set by 'latitude'.
+#'@param na.rm A logical value indicating whether NA values should be removed 
+#'  from data. It is FALSE by default. 
+#'@param ncores An integer value indicating the number of cores to use in 
+#'  parallel computation.
+#' 
+#'@examples 
+#'dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1)
+#' 
+#'exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+#'exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+#'exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+#'end_year <- 2012
+#'dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), 
+#'                          paste0(2010:end_year, "-09-15"), 
+#'                          paste0(2010:end_year, "-10-16")), "UTC")
+#'dim(dates_exp) <- c(syear = 3, time = 3)
+#'
+#'lat <- c(40)
+#'
+#'exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr)
+#'
+#'res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp)
+#'
+#'@import SPEI
+#'@import lubridate
+#'@import multiApply
+#'@export 
+CST_PeriodPET <- function(data, pet_method = 'hargreaves', 
+                          time_dim = 'syear', leadtime_dim = 'time',
+                          lat_dim = 'latitude', na.rm = FALSE, 
+                          ncores = NULL) {
+  # Check 's2dv_cube'
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!all(sapply(data, function(x) inherits(x, 's2dv_cube')))) {
+    stop("Parameter 'data' must be a list of 's2dv_cube' class.")
+  }
+  # latitude
+  if (!any(names(data[[1]]$coords) %in% .KnownLatNames())) {
+    stop("Spatial coordinate names of parameter 'data' do not match any ",
+         "of the names accepted by the package.")
+  }
+  # Dates
+  dates_exp <- data[[1]]$attrs$Dates
+  if (!'Dates' %in% names(data[[1]]$attrs)) {
+    stop("Element 'Dates' is not found in 'attrs' list of 'data'. ", 
+         "See 's2dv_cube' object description in README file for more ", 
+         "information.")
+  }
+  lat_dim <- names(data[[1]]$coords)[[which(names(data[[1]]$coords) %in% .KnownLatNames())]]
+
+  res <- PeriodPET(data = lapply(data, function(x) x$data), 
+                   dates = data[[1]]$attrs$Dates, 
+                   lat = data[[1]]$coords[[lat_dim]], 
+                   pet_method = pet_method, time_dim = time_dim, 
+                   leadtime_dim = leadtime_dim, lat_dim = lat_dim, 
+                   na.rm = na.rm, ncores = ncores)
+  # Add metadata
+  source_files <- lapply(data, function(x) {x$attrs$source_files})
+  coords <- data[[1]]$coords
+  Dates <- data[[1]]$attrs$Dates
+  metadata <- data[[1]]$attrs$Variable$metadata
+  metadata_names <- intersect(names(dim(res)), names(metadata))
+  suppressWarnings(
+    res <- CSTools::s2dv_cube(data = res, coords = coords, 
+                              varName = paste0('PET'), 
+                              metadata = metadata[metadata_names],
+                              Dates = Dates, 
+                              source_files = source_files, 
+                              when = Sys.time())
+  )
+  return(res)
+}
+
+#'Compute the Potential Evapotranspiration
+#' 
+#'Compute the Potential evapotranspiration (PET) that is the amount of 
+#'evaporation and transpiration that would occur if a sufficient water source 
+#'were available. This function calculate PET according to the Thornthwaite, 
+#'Hargreaves or Hargreaves-modified equations. 
+#' 
+#'For more information on the SPEI calculation, see functions 
+#'PeriodStandardization and PeriodAccumulation.
+#'
+#'@param data A named list with the needed \code{s2dv_cube} objects containing 
+#'  the seasonal forecast experiment in the data element for each variable. 
+#'  Specific variables are needed for each method used in computing the 
+#'  Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+#'  variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for 
+#'  'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' 
+#'  'tmean' is required. The units for temperature variables
+#'  ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+#'  precipitation ('pr') need to be in mm/month. 
+#'@param dates An array of temporal dimensions containing the Dates of 
+#'  'data'. It must be of class 'Date' or 'POSIXct'.
+#'@param lat A numeric vector containing the latitude values of 'data'.
+#'@param pet_method A character string indicating the method used to compute 
+#'  the potential evapotranspiration. The accepted methods are:
+#'  'hargreaves' and 'hargreaves_modified', that require the data to have 
+#'  variables tmin and tmax; and 'thornthwaite', that requires variable 
+#'  'tmean'. 
+#'@param time_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'syear'. 
+#'@param leadtime_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'time'. 
+#'@param lat_dim A character string indicating the name of the latitudinal 
+#'  dimension. By default it is set by 'latitude'.
+#'@param na.rm A logical value indicating whether NA values should be removed 
+#'  from data. It is FALSE by default. 
+#'@param ncores An integer value indicating the number of cores to use in 
+#'  parallel computation.
+#' 
+#'@examples 
+#'dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1)
+#' 
+#'exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+#'exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+#'exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+#'end_year <- 2012
+#'dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), 
+#'                          paste0(2010:end_year, "-09-15"), 
+#'                          paste0(2010:end_year, "-10-16")), "UTC")
+#'dim(dates_exp) <- c(syear = 3, time = 3)
+#'
+#'lat <- c(40)
+#'
+#'exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr)
+#'
+#'res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp)
+#'
+#'@import SPEI
+#'@import lubridate
+#'@import multiApply
+#'@export 
+PeriodPET <- function(data, dates, lat, pet_method = 'hargreaves', 
+                      time_dim = 'syear', leadtime_dim = 'time',
+                      lat_dim = 'latitude', na.rm = FALSE, 
+                      ncores = NULL) {
+
+  # Initial checks
+  ## data
+  if (!inherits(data, 'list')) {
+    stop("Parameter 'data' needs to be a named list with the needed variables.")
+  }
+  if (is.null(names(data))) {
+    stop("Parameter 'data' needs to be a named list with the variable names.")
+  }
+  if (any(sapply(data, function(x) is.null(names(dim(x)))))) {
+    stop("Parameter 'data' needs to be a list of arrays with dimension names.")
+  }
+  dims <- lapply(data, function(x) dim(x))
+  first_dims <- dims[[1]]
+  all_equal <- all(sapply(dims[-1], function(x) identical(first_dims, x)))
+  if (!all_equal) {
+    stop("Parameter 'data' variables need to have the same dimensions.")
+  }
+  # lat
+  if (!is.numeric(lat)) {
+    stop("Parameter 'lat' must be numeric.")
+  }
+  if (!lat_dim %in% names(dims[[1]])) {
+    stop("Parameter 'data' must have 'lat_dim' dimension.")
+  }
+  if (any(sapply(dims, FUN = function(x) x[lat_dim] != length(lat)))) {
+    stop("Parameter 'lat' needs to have the same length of latitudinal", 
+         "dimension of all the variables arrays in 'data'.")
+  }
+
+  # data (2)
+  if (all(c('tmin', 'tmax', 'pr') %in% names(data))) {
+    # hargreaves modified: 'tmin', 'tmax', 'pr' and 'lat'
+    if (!(pet_method %in% c('hargreaves_modified', 'hargreaves'))) {
+      warning("Parameter 'pet_method' needs to be 'hargreaves' or ", 
+              "'hargreaves_modified'. It is set to 'hargreaves_modified'.")
+      pet_method <- 'hargreaves_modified'
+    }
+  } else if (all(c('tmin', 'tmax') %in% names(data))) { 
+    if (!(pet_method %in% c('hargreaves'))) {
+      warning("Parameter 'pet_method' will be set as 'hargreaves'.")
+      pet_method <- 'hargreaves'
+    }
+  } else if (c('tmean') %in% names(data)) {
+    # thornthwaite: 'tmean' (mean), 'lat'
+    if (!(pet_method == 'thornthwaite')) {
+      warning("Parameter 'pet_method' it is set to be 'thornthwaite'.")
+      pet_method <- 'thornthwaite'
+    }
+  } else {
+    stop("Parameter 'data' needs to be a named list with accepted ", 
+          "variable names. See documentation.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!all(sapply(data, function(x) time_dim %in% names(dim(x))))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+  ## leadtime_dim
+  if (!is.character(leadtime_dim) | length(leadtime_dim) != 1) {
+    stop("Parameter 'leadtime_dim' must be a character string.")
+  }
+  if (!all(sapply(data, function(x) leadtime_dim %in% names(dim(x))))) {
+    stop("Parameter 'leadtime_dim' is not found in 'data' dimension.")
+  }
+  ## lat_dim
+  if (!is.character(lat_dim) | length(lat_dim) != 1) {
+    stop("Parameter 'lat_dim' must be a character string.")
+  }
+  if (!all(sapply(data, function(x) lat_dim %in% names(dim(x))))) {
+    stop("Parameter 'lat_dim' is not found in 'data' dimension.")
+  }
+  # dates
+  if (is.null(dates)) {
+    stop("Parameter 'dates' is missing, dates must be provided.")
+  }
+  if (!(is.Date(dates)) & !(is.POSIXct(dates))) {
+    stop("Parameter 'dates' is not of the correct class, ", 
+         "only 'Date' and 'POSIXct' classes are accepted.")
+  }
+  if (!time_dim %in% names(dim(dates)) | !leadtime_dim %in% names(dim(dates))) {
+    stop("Parameter 'dates' must have 'time_dim' and 'leadtime_dim' ",
+         "dimension.")
+  }
+  if (!all(dim(data[[1]])[c(time_dim, leadtime_dim)] == 
+           dim(dates)[c(time_dim, leadtime_dim)])) {
+    stop("Parameter 'dates' needs to have the same length as 'time_dim' ", 
+         "and 'leadtime_dim' as 'data'.")
+  }
+  ## na.rm
+  if (!is.logical(na.rm) | length(na.rm) > 1) {
+    stop("Parameter 'na.rm' must be one logical value.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+
+  # complete dates
+  dates_monthly <- .datesmask(dates)
+
+  lat_mask <- array(lat, dim = c(1, length(lat)))
+  names(dim(lat_mask)) <- c('dat', lat_dim)
+
+  # extract mask of NA locations to return to NA the final result
+  mask_na <- array(1, dim = dim(data[[1]]))
+  if (pet_method == 'hargreaves') {
+    varnames <- c('tmax', 'tmin')
+    mask_na[which(is.na(data$tmax))] <- 0
+    mask_na[which(is.na(data$tmin))] <- 0
+  } else if (pet_method == 'hargreaves_modified') {
+    varnames <- c('tmax', 'tmin', 'pr')
+    mask_na[which(is.na(data$tmax))] <- 0
+    mask_na[which(is.na(data$tmin))] <- 0
+    mask_na[which(is.na(data$pr))] <- 0
+  } else if (pet_method == 'thornthwaite') {
+    varnames <- c('tmean')
+    mask_na[which(is.na(data$tmean))] <- 0
+  }
+
+  # replace NA with 0
+  for (dd in 1:length(data)) {
+    data[[dd]][which(is.na(data[[dd]]))] <- 0
+  }
+
+  # prepare data
+  target_dims_data <- lapply(data[varnames], function(x) rep(c(leadtime_dim, time_dim), 1))
+  pet <- Apply(data = c(list(lat_mask = lat_mask), data[varnames]), 
+               target_dims = c(list(lat_mask = 'dat'), target_dims_data), 
+               fun = .pet, 
+               dates_monthly = dates_monthly, pet_method = pet_method, 
+               leadtime_dim = leadtime_dim, time_dim = time_dim, 
+               output_dims = c(leadtime_dim, time_dim),
+               ncores = ncores)$output1
+  # reorder dims in pet_estimated
+  pos <- match(names(dim(data[[1]])), names(dim(pet)))
+  pet <- aperm(pet, pos)
+
+  # restore original NAs from mask_na
+  pet[which(mask_na == 0)] <- NA
+
+  return(pet)
+}
+
+.pet <- function(lat_mask, data2, data3 = NULL, data4 = NULL, 
+                 dates_monthly, pet_method = 'hargreaves', 
+                 leadtime_dim = 'time', time_dim = 'syear') {
+ 
+  dims <- dim(data2)
+
+  # create a vector from data but adding 0 to achive complete time series 
+  # of the considered period 
+  # (starting in January of the first year) so that the solar radiation 
+  # estimation is computed in each case for the correct month
+
+  if (!is.null(data2)) {
+    data_tmp <- as.vector(data2)
+    data2 <- array(0, dim = length(dates_monthly))
+    count <- 1
+    for (dd in 1:length(dates_monthly)) {
+      if (dates_monthly[dd] == 1) {
+        data2[dd] <- data_tmp[count]
+        count <- count + 1
+      }
+    }
+    rm(data_tmp)
+  }
+  if (!is.null(data3)) {
+    data_tmp <- as.vector(data3)
+    data3 <- array(0, dim = length(dates_monthly))
+    count <- 1
+    for (dd in 1:length(dates_monthly)) {
+      if (dates_monthly[dd] == 1) {
+        data3[dd] <- data_tmp[count]
+        count <- count + 1
+      }
+    }
+    rm(data_tmp)
+  }
+  if (!is.null(data4)) {
+    data_tmp <- as.vector(data4)
+    data4 <- array(0, dim = length(dates_monthly))
+    count <- 1
+    for (dd in 1:length(dates_monthly)) {
+      if (dates_monthly[dd] == 1) {
+        data4[dd] <- data_tmp[count]
+        count <- count + 1
+      }
+    }
+    rm(data_tmp)
+  }
+  if (pet_method == 'hargreaves') {
+    pet <- hargreaves(Tmin = as.vector(data3), Tmax = as.vector(data2), 
+                      lat = lat_mask, na.rm = FALSE, verbose = FALSE)
+    # line to return the vector to the size of the actual original data
+    pet <- array(pet[which(dates_monthly == 1)], dim = dims)
+  }
+
+  if (pet_method == 'hargreaves_modified') {
+    pet <- hargreaves(Tmin = as.vector(data3), Tmax = as.vector(data2), 
+                      lat = lat_mask, Pre = as.vector(data4), na.rm = FALSE,
+                      verbose = FALSE)
+    pet <- array(pet[which(dates_monthly == 1)], dim = dims)
+  }
+  
+  if (pet_method == 'thornthwaite') {
+    pet <- thornthwaite(as.vector(data2), lat = lat_mask, na.rm = TRUE,
+                        verbose = FALSE)
+    # line to return the vector to the size of the actual original data
+    pet <- array(pet[which(dates_monthly == 1)], dim = dims)
+  }
+  return(pet)
+}
+
+
+.datesmask <- function(dates) {
+  ini <- as.Date(paste(lubridate::year(min(dates)), 01, 01, sep = '-'))
+  end <- as.Date(paste(lubridate::year(max(dates)), 12, 31, sep = '-'))
+  daily <- as.Date(ini:end)
+  monthly <- daily[which(lubridate::day(daily) == 1)]
+  dates_mask <- array(0, dim = length(monthly))
+  for (dd in 1:length(dates)) {
+    ii <- which(monthly == as.Date(paste(lubridate::year(dates[dd]), 
+                                         lubridate::month(dates[dd]), 
+                                         01, sep = '-')))
+    dates_mask[ii] <- 1
+  }
+  return(dates_mask)
+}
\ No newline at end of file
diff --git a/R/PeriodSPEI.R b/R/PeriodSPEI.R
new file mode 100644
index 0000000000000000000000000000000000000000..c2a04ecc0300c6f20ecd74ee0aa81c406421967d
--- /dev/null
+++ b/R/PeriodSPEI.R
@@ -0,0 +1,1130 @@
+#'Compute the Standardised Precipitation-Evapotranspiration Index
+#' 
+#'Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) 
+#'that is a multiscalar drought index based on climatic data. It can be used for 
+#'determining the onset, duration and magnitude of drought conditions with 
+#'respect to normal conditions in a variety of natural and managed systems such 
+#'as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI 
+#'is to compare the highest possible evapotranspiration with the current water 
+#'availability. The SPEI uses the monthly (or weekly) difference between 
+#'precipitation and potential evapotranspiration. This represents a simple 
+#'climatic water balance which is calculated at different time scales to obtain 
+#'the SPEI. This function is build to work be compatible with other tools in 
+#'that work with 's2dv_cube' object class. The input data must be this object 
+#'class. If you don't work with 's2dv_cube', see PeriodSPEI.
+#' 
+#'Next, some specifications for the calculation of this indicator will be 
+#'discussed. On the one hand, the model to be used to calculate potential 
+#'evapotranspiration is specified with the pet_method parameter (hargreaves, 
+#'hargraves modified or thornwhite). On the other hand, to choose the time scale 
+#'in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the 
+#'accum parameter, where you must indicate the number of time steps you want to 
+#'accumulate throughout leadtime_dim. Since the accumulation is done for the 
+#'elapsed time steps, there will be no complete accumulations until reaching the 
+#'time instant equal to the value of the parameter. For this reason, in the 
+#'result, we will find that for the dimension where the accumulation has been 
+#'carried out, the values of the array will be NA since they do not include 
+#'complete accumulations. Also, there is a parameter to specify if the 
+#'standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a 
+#'probability distribution to transform the original values to standardized 
+#'units that are comparable in space and time and at different SPEI time scales. 
+#'The na.rm parameter is a logical parameter used to decide whether to remove 
+#'the NA values from the data before doing the calculation. It must be taken 
+#'into account that if na.rm == FALSE and there is some NA value in the specific 
+#'coordinates which the SPEI is computed, standardization cannot be carried out 
+#'for those coordinates and therefore, the result will be filled with NA for the 
+#'specific coordinates. However, when na.rm == TRUE,  if the amount of data for 
+#'those specific coordinates is smaller than 4, it will not be possible to carry 
+#'out because we will not have enough data and the result will be also filled 
+#'with NAs for that coordinates.
+#' 
+#'@param exp A named list with the needed \code{s2dv_cube} objects containing 
+#'  the seasonal forecast experiment in the data element for each variable. 
+#'  Specific variables are needed for each method used in computing the 
+#'  Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+#'  variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and 
+#'  'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. 
+#'  Variable 'pr' is always needed. The units for temperature variables
+#'  ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+#'  precipitation ('pr') need to be in mm/month. 
+#'@param exp_cor A named list with the needed \code{s2dv_cube} objects for each 
+#'  variable in which the quantile PeriodSPEI should be applied. If it is not 
+#'  specified, the PeriodSPEI is calculated from object 'exp'.
+#'@param time_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'syear'. 
+#'@param leadtime_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'time'. 
+#'@param memb_dim A character string indicating the name of the dimension in 
+#'  which the ensemble members are stored. When set it to NULL, threshold is 
+#'  computed for individual members.
+#'@param lat_dim A character string indicating the name of the latitudinal 
+#'  dimension. By default it is set by 'latitude'.
+#'@param accum An integer value indicating the number of months for the 
+#'  accumulation for each variable. When it is greater than 1, the result will 
+#'  be filled with NA until the accum time_dim dimension number due to the 
+#'  accumulation to previous months.
+#'@param ref_period A list with two numeric values with the starting and end 
+#'  points of the reference period used for computing the index. The default 
+#'  value is NULL indicating that the first and end values in data will be 
+#'  used as starting and end points.
+#'@param params A multidimensional array with named dimensions for computing the 
+#'  SPEI. This option overrides computation of fitting parameters. It needs 
+#'  to have the following dimensions: same leadtime dimension of 'exp' 
+#'  (specified in 'leadtime_dim'); time dimension of length 1 (specified in 
+#'  'time_dim'); and a dimension named 'coef' with the length of the 
+#'  coefficients needed for the used distribution (for 'Gamma' coef dimension is 
+#'  of lenght 2, for 'log-Logistic' or 'PearsonIII' is of length 3). It can't 
+#'  have member dimension (specified in 'memb_dim').
+#'@param pet_exp A multidimensional array containing the Potential  
+#'  EvapoTranspiration data of 'exp'. It must have the same dimensions of the 
+#'  variable 'pr' of 'exp'. If it is NULL it is calculated using the provided 
+#'  variables with specified 'pet_method'. It is NULL by default.
+#'@param pet_expcor A multidimensional array containing the Potential  
+#'  EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of 
+#'  the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the 
+#'  provided variables with specified 'pet_method'. It is NULL by default.
+#'@param standardization A logical value indicating wether the standardization 
+#'  is computed.
+#'@param cross_validation A logical value indicating if cross validation is 
+#'  done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and 
+#'  is not provided. It is FALSE by default.
+#'@param pet_method A character string indicating the method used to compute 
+#'  the potential evapotranspiration. The accepted methods are:
+#'  'hargreaves' and 'hargreaves_modified', that require the data to have 
+#'  variables tmin and tmax; and 'thornthwaite', that requires variable 
+#'  'tmean'.  
+#'@param method A character string indicating the standardization method used. 
+#'  If can be: 'parametric' or 'non-parametric'. It is set to 'parametric' by 
+#'  default.
+#'@param distribution A character string indicating the name of the distribution 
+#'  function to be used for computing the SPEI. The accepted names are: 
+#'  'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The 
+#'  'Gamma' method only works when only precipitation is provided and other 
+#'  variables are 0 because it is positive defined (SPI indicator).
+#'@param param_error A numeric value with the error accepted.
+#'@param handle_infinity A logical value wether to return Infinite values (TRUE)
+#'  or not (FALSE).
+#'@param na.rm A logical value indicating whether NA values should be removed 
+#'  from data. It is FALSE by default. If it is FALSE and there are NA values, 
+#'  (if standardization is TRUE) all values of other dimensions except time_dim 
+#'  and leadtime_dim will be set to NA directly. On the other hand, if it is 
+#'  TRUE, if the data from other dimensions except time_dim and leadtime_dim is 
+#'  not reaching 4 values, it is not enough values to estimate the parameters 
+#'  and the result will include NA.
+#'@param ncores An integer value indicating the number of cores to use in 
+#'  parallel computation.
+#' 
+#'@return An 's2dv_cube' object containing the SPEI multidimensional array in 
+#'element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' 
+#'will be provided. The parameters of the standardization will only be returned 
+#'if 'return_params' is TRUE. The SPEI will only be computed if 
+#''standardization' is TRUE. If 'standardization' is FALSE, only the climatic 
+#'water balance (precipitation minus evapotranspiration) will be returned. The 
+#'resultant arrays will have the same dimensions as the initial input data. The 
+#'other elements in the 's2dv_cube' will be updated with the combined 
+#'information of the input data arrays.
+#' 
+#'@examples 
+#'dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, 
+#'          latitude = 2, longitude = 1, ensemble = 25)
+#'dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, 
+#'             latitude = 2, longitude = 1, ensemble = 15)
+#'
+#'dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), 
+#'                          paste0(2010:2015, "-09-15"), 
+#'                          paste0(2010:2015, "-10-16")), 'UTC')
+#'dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3)
+#'dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), 
+#'                             paste0(2020, "-10-16")), 'UTC')
+#'dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3)
+#'lat <- c(40,40.1)
+#'
+#'exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims)
+#'exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims)
+#'exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims)
+#'
+#'expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor)
+#'expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor)
+#'expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor)
+#'
+#'exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+#'exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 
+#'                'pr' = expcor_pr)
+#'
+#'exp <- lapply(exp, CSTools::s2dv_cube, coords = list(latitude = lat), 
+#'              Dates = dates_exp)
+#'exp_cor <- lapply(exp_cor, CSTools::s2dv_cube, coords = list(latitude = lat), 
+#'                  Dates = dates_expcor)
+#'
+#'res <- CST_PeriodSPEI(exp = exp, exp_cor = exp_cor)
+#' 
+#'@import multiApply
+#'@import ClimProjDiags
+#'@import SPEI
+#'@import zoo
+#'@import TLMoments
+#'@import lmomco
+#'@import lmom
+#'@import lubridate
+#'@import CSTools
+#'@export
+CST_PeriodSPEI <- function(exp, exp_cor = NULL, 
+                           time_dim = 'syear', leadtime_dim = 'time', 
+                           memb_dim =  'ensemble', lat_dim = 'latitude',
+                           accum = 1, ref_period = NULL, params = NULL,
+                           pet_exp = NULL, pet_expcor = NULL, 
+                           standardization = TRUE, cross_validation = FALSE, 
+                           pet_method = 'hargreaves', method = 'parametric', 
+                           distribution = 'log-Logistic',
+                           param_error = -9999, handle_infinity = FALSE, 
+                           return_params = FALSE, na.rm = FALSE, 
+                           ncores = NULL) {
+                          
+  # Check 's2dv_cube'
+  if (is.null(exp)) {
+    stop("Parameter 'exp' cannot be NULL.")
+  }
+  if (!all(sapply(exp, function(x) inherits(x, 's2dv_cube')))) {
+    stop("Parameter 'exp' must be a list of 's2dv_cube' class.")
+  }
+  if (!is.null(exp_cor)) {
+    if (!all(sapply(exp_cor, function(x) inherits(x, 's2dv_cube')))) {
+      stop("Parameter 'exp_cor' must be a list of 's2dv_cube' class.")
+    }
+  }
+  # latitude
+  if (!any(names(exp[[1]]$coords) %in% .KnownLatNames())) {
+    stop("Spatial coordinate names of parameter 'exp' do not match any ",
+         "of the names accepted by the package.")
+  }
+  # Dates
+  dates_exp <- exp[[1]]$attrs$Dates
+  if (!'Dates' %in% names(exp[[1]]$attrs)) {
+    stop("Element 'Dates' is not found in 'attrs' list of 'exp'. ", 
+         "See 's2dv_cube' object description in README file for more ", 
+         "information.")
+  }
+
+  if (!is.null(exp_cor)) {
+    if (!'Dates' %in% names(exp_cor[[1]]$attrs)) {
+      stop("Element 'Dates' is not found in 'attrs' list of 'exp_cor'. ", 
+           "See 's2dv_cube' object description in README file for more ", 
+           "information.")
+    }
+  }
+
+  lat_dim <- names(exp[[1]]$coords)[[which(names(exp[[1]]$coords) %in% .KnownLatNames())]]
+
+  res <- PeriodSPEI(exp = lapply(exp, function(x) x$data), 
+                    dates_exp = exp[[1]]$attrs$Dates, 
+                    lat = exp[[1]]$coords[[lat_dim]], 
+                    exp_cor = if (is.null(exp_cor)) {NULL} else {
+                                  lapply(exp_cor, function(x) x$data)}, 
+                    dates_expcor = exp_cor[[1]]$attrs$Dates, 
+                    time_dim = time_dim, leadtime_dim = leadtime_dim, 
+                    memb_dim = memb_dim, lat_dim = lat_dim, 
+                    accum = accum, ref_period = ref_period, params = params, 
+                    pet_exp = pet_exp, pet_expcor = pet_expcor, 
+                    standardization = standardization, 
+                    cross_validation = cross_validation, 
+                    pet_method = pet_method, method = method, 
+                    distribution = distribution, 
+                    param_error = param_error, handle_infinity = handle_infinity, 
+                    return_params = return_params, na.rm = na.rm, 
+                    ncores = ncores)
+
+  if (!is.null(exp_cor)) {
+    source_files_expcor <- lapply(exp_cor, function(x) {x$attrs$source_files})
+    source_files <- lapply(exp, function(x) {x$attrs$source_files})
+    source_files <- c(exp = source_files, exp_cor = source_files_expcor)
+    coords <- exp_cor[[1]]$coords
+    Dates <- exp_cor[[1]]$attrs$Dates
+    metadata <- exp_cor[[1]]$attrs$Variable$metadata
+    metadata_names <- names(metadata)
+  } else {
+    source_files <- lapply(exp, function(x) {x$attrs$source_files})
+    coords <- exp[[1]]$coords
+    Dates <- exp[[1]]$attrs$Dates
+    metadata <- exp[[1]]$attrs$Variable$metadata
+    metadata_names <- names(metadata) 
+  }
+
+  if (standardization) {
+    varname <- 'SPEI'
+  } else {
+    varname <- 'Precipitation minus accumulated PET'
+  }
+
+  if (return_params & standardization) {
+    metadata_names <- intersect(names(dim(res[[1]])), metadata_names)
+    suppressWarnings(
+      res[[1]] <- CSTools::s2dv_cube(data = res[[1]], coords = coords, 
+                                     varName = varname, 
+                                     metadata = metadata[metadata_names],
+                                     Dates = Dates, 
+                                     source_files = source_files,
+                                     when = Sys.time())
+    )
+    return(list(spei = res[[1]], params = res[[2]]))
+  } else {
+    metadata_names <- intersect(names(dim(res)), metadata_names)
+    suppressWarnings(
+      res <- CSTools::s2dv_cube(data = res, coords = coords, 
+                                varName = varname, 
+                                metadata = metadata[metadata_names],
+                                Dates = Dates, 
+                                source_files = source_files, 
+                                when = Sys.time())
+    )
+    return(res)
+  }
+}
+
+
+#'Compute the Standardised Precipitation-Evapotranspiration Index
+#' 
+#'Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) 
+#'that is a multiscalar drought index based on climatic data. It can be used for 
+#'determining the onset, duration and magnitude of drought conditions with 
+#'respect to normal conditions in a variety of natural and managed systems such 
+#'as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI 
+#'is to compare the highest possible evapotranspiration with the current water 
+#'availability. The SPEI uses the monthly (or weekly) difference between 
+#'precipitation and potential evapotranspiration. This represents a simple 
+#'climatic water balance which is calculated at different time scales to obtain 
+#'the SPEI.
+#' 
+#'Next, some specifications for the calculation of this indicator will be 
+#'discussed. On the one hand, the model to be used to calculate potential 
+#'evapotranspiration is specified with the pet_method parameter (hargreaves, 
+#'hargraves modified or thornwhite). On the other hand, to choose the time scale 
+#'in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the 
+#'accum parameter, where you must indicate the number of time steps you want to 
+#'accumulate throughout leadtime_dim. Since the accumulation is done for the 
+#'elapsed time steps, there will be no complete accumulations until reaching the 
+#'time instant equal to the value of the parameter. For this reason, in the 
+#'result, we will find that for the dimension where the accumulation has been 
+#'carried out, the values of the array will be NA since they do not include 
+#'complete accumulations. Also, there is a parameter to specify if the 
+#'standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a 
+#'probability distribution to transform the original values to standardized 
+#'units that are comparable in space and time and at different SPEI time scales. 
+#'The na.rm parameter is a logical parameter used to decide whether to remove 
+#'the NA values from the data before doing the calculation. It must be taken 
+#'into account that if na.rm == FALSE and there is some NA value in the specific 
+#'coordinates which the SPEI is computed, standardization cannot be carried out 
+#'for those coordinates and therefore, the result will be filled with NA for the 
+#'specific coordinates. However, when na.rm == TRUE,  if the amount of data for 
+#'those specific coordinates is smaller than 4, it will not be possible to carry 
+#'out because we will not have enough data and the result will be also filled 
+#'with NAs for that coordinates.
+#' 
+#'@param exp A named list with multidimensional array objects containing 
+#'  the seasonal forecast experiment in the data element for each variable. 
+#'  Specific variables are needed for each method used in computing the 
+#'  Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+#'  variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and 
+#'  'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. 
+#'  Variable 'pr' is always needed. The units for temperature variables
+#'  ('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+#'  precipitation ('pr') need to be in mm/month. 
+#'@param dates_exp An array of temporal dimensions containing the Dates of 
+#'  'exp'. It must be of class 'Date' or 'POSIXct'.
+#'@param lat A numeric vector containing the latitude values of 'exp'.
+#'@param exp_cor A named list with multidimensional array objects for each 
+#'  variable in which the quantile PeriodSPEI should be applied. If it is not 
+#'  specified, the PeriodSPEI is calculated from object 'exp'.
+#'@param dates_expcor An array of temporal dimensions containing the Dates of 
+#'  'exp_cor'. It must be of class 'Date' or 'POSIXct'.
+#'@param time_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'syear'. 
+#'@param leadtime_dim A character string indicating the name of the temporal 
+#'  dimension. By default, it is set to 'time'. 
+#'@param memb_dim A character string indicating the name of the dimension in 
+#'  which the ensemble members are stored. When set it to NULL, threshold is 
+#'  computed for individual members.
+#'@param lat_dim A character string indicating the name of the latitudinal 
+#'  dimension. By default it is set by 'latitude'.
+#'@param accum accum An integer value indicating the number of months for the 
+#'  accumulation for each variable. When it is greater than 1, the result will 
+#'  be filled with NA until the accum time_dim dimension number due to the 
+#'  accumulation to previous months.
+#'@param ref_period A list with two numeric values with the starting and end 
+#'  points of the reference period used for computing the index. The default 
+#'  value is NULL indicating that the first and end values in data will be 
+#'  used as starting and end points.
+#'@param params A multidimensional array with named dimensions for computing the 
+#'  SPEI. This option overrides computation of fitting parameters. It needs 
+#'  to be of same time dimension (specified in 'time_dim') of 'exp' and a 
+#'  dimension named 'coef' with the length of the coefficients needed for the 
+#'  used distribution (for 'Gamma' coef dimension is of lenght 2, for 
+#'  'log-Logistic' or 'PearsonIII' is of length 3). It also needs to have a 
+#'  leadtime dimension (specified in 'leadtime_dim') of length 1.
+#'@param pet_exp A multidimensional array containing the Potential  
+#'  EvapoTranspiration data of 'exp'. It must have the same dimensions of the 
+#'  variable 'pr' of 'exp'. If it is NULL it is calculated using the provided 
+#'  variables with specified 'pet_method'. It is NULL by default.
+#'@param pet_expcor A multidimensional array containing the Potential  
+#'  EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of 
+#'  the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the 
+#'  provided variables with specified 'pet_method'. It is NULL by default.
+#'@param standardization A logical value indicating wether the standardization 
+#'  is computed.
+#'@param cross_validation A logical value indicating if cross validation is 
+#'  done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and 
+#'  is not provided. It is FALSE by default.
+#'@param pet_method A character string indicating the method used to compute 
+#'  the potential evapotranspiration. The accepted methods are:
+#'  'hargreaves' and 'hargreaves_modified', that require the data to have 
+#'  variables tmin and tmax; and 'thornthwaite', that requires variable 
+#'  'tmean'.  
+#'@param method A character string indicating the standardization method used. 
+#'  If can be: 'parametric' or 'non-parametric'.
+#'@param distribution A character string indicating the name of the distribution 
+#'  function to be used for computing the SPEI. The accepted names are: 
+#'  'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The 
+#'  'Gamma' method only works when only precipitation is provided and other 
+#'  variables are 0 because it is positive defined (SPI indicator).
+#'@param param_error A numeric value with the error accepted.
+#'@param handle_infinity A logical value wether to return Infinite values (TRUE)
+#'  or not (FALSE).
+#'@param return_params A logical value indicating wether to return parameters 
+#'  array (TRUE) or not (FALSE). It is FALSE by default.
+#'@param na.rm A logical value indicating whether NA values should be removed 
+#'  from data. It is FALSE by default. If it is FALSE and there are NA values, 
+#'  (if standardization is TRUE) all values of other dimensions except time_dim 
+#'  and leadtime_dim will be set to NA directly. On the other hand, if it is 
+#'  TRUE, if the data from other dimensions except time_dim and leadtime_dim is 
+#'  not reaching 4 values, it is not enough values to estimate the parameters 
+#'  and the result will include NA.
+#'@param ncores An integer value indicating the number of cores to use in 
+#'  parallel computation.
+#' 
+#'@return An 's2dv_cube' object containing the SPEI multidimensional array in 
+#'element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' 
+#'will be provided. The parameters of the standardization will only be returned 
+#'if 'return_params' is TRUE. The SPEI will only be computed if 
+#''standardization' is TRUE. If 'standardization' is FALSE, only the climatic 
+#'water balance (precipitation minus evapotranspiration) will be returned. The 
+#'resultant arrays will have the same dimensions as the initial input data.
+#' 
+#'@examples 
+#'dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, 
+#'          latitude = 2, longitude = 1, ensemble = 25)
+#'dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, 
+#'             latitude = 2, longitude = 1, ensemble = 15)
+#'exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims)
+#'exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims)
+#'exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims)
+#' 
+#'expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor)
+#'expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor)
+#'expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor)
+#' 
+#'dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), 
+#'                       paste0(2010:2015, "-09-15"), 
+#'                       paste0(2010:2015, "-10-16")))
+#'dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3)
+#'dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), 
+#'                             paste0(2020, "-09-15"), 
+#'                             paste0(2020, "-10-16")), 'UTC')
+#'dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3)
+#'lat <- c(40,40.1)
+#' 
+#'exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+#'exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 
+#'                'pr' = expcor_pr)
+#'res <- PeriodSPEI(exp = exp, exp_cor = exp_cor, lat = lat, 
+#'                  dates_exp = dates_exp, dates_expcor = dates_expcor)
+#' 
+#'@import multiApply
+#'@import ClimProjDiags
+#'@import SPEI
+#'@import zoo
+#'@import TLMoments
+#'@import lmomco
+#'@import lmom
+#'@import lubridate
+#'@export
+PeriodSPEI <- function(exp, dates_exp, lat, 
+                       exp_cor = NULL, dates_expcor = NULL,
+                       time_dim = 'syear', leadtime_dim = 'time', 
+                       memb_dim = 'ensemble', lat_dim = 'latitude',
+                       accum = 1, ref_period = NULL, params = NULL, 
+                       pet_exp = NULL, pet_expcor = NULL, 
+                       standardization = TRUE, cross_validation = FALSE,
+                       pet_method = 'hargreaves', method = 'parametric', 
+                       distribution = 'log-Logistic',
+                       param_error = -9999, handle_infinity = FALSE, 
+                       return_params = FALSE, na.rm = FALSE, ncores = NULL) {
+
+  # Initial checks
+  ## exp
+  if (!inherits(exp, 'list')) {
+    stop("Parameter 'exp' needs to be a named list with the needed variables.")
+  }
+  if (is.null(names(exp))) {
+    stop("Parameter 'exp' needs to be a named list with the variable names.")
+  }
+  if (any(sapply(exp, function(x) is.null(names(dim(x)))))) {
+    stop("Parameter 'exp' needs to be a list of arrays with dimension names.")
+  }
+  dims <- lapply(exp, function(x) dim(x))
+  first_dims <- dims[[1]]
+  all_equal <- all(sapply(dims[-1], function(x) identical(first_dims, x)))
+  if (!all_equal) {
+    stop("Parameter 'exp' variables need to have the same dimensions.")
+  }
+
+  ## exp_cor
+  if (!is.null(exp_cor)) {
+    if (!inherits(exp_cor, 'list')) {
+      stop("Parameter 'exp_cor' needs to be a named list with the needed ", 
+           "variables if it is not NULL.")
+    }
+    if (is.null(names(exp_cor))) {
+      stop("Parameter 'exp_cor' needs to be a named list with the variable names.")
+    }
+    if (any(sapply(exp_cor, function(x) is.null(names(dim(x)))))) {
+      stop("Parameter 'exp_cor' needs to be a list of arrays with dimension names.")
+    }
+    dimscor <- lapply(exp_cor, function(x) dim(x))
+    first_dimscor <- dimscor[[1]]
+    all_equal <- all(sapply(dimscor[-1], function(x) identical(first_dimscor, x)))
+    if (!all_equal) {
+      stop("Parameter 'exp_cor' variables need to have the same dimensions.")
+    }
+  }
+  # Variable checks
+  ## exp (2)
+  pet <- vector("list", 2)
+  if (!('pr' %in% names(exp))) {
+    stop("Variable 'pr' is not included in 'exp'.")
+  }
+  ## exp_cor (2)
+  if (!is.null(exp_cor)) {
+    if (!('pr' %in% names(exp_cor))) {
+      stop("Variable 'pr' is not included in 'exp_cor'.")
+    }
+    if (length(pet_method) == 1) {
+      pet_method <- rep(pet_method, 2)
+    }
+  }
+  
+  ## pet_exp
+  if (!is.null(pet_exp)) {
+    if (length(dim(exp[['pr']])) != length(dim(pet_exp))) {
+      stop("Parameter 'pet_exp' must have the same length of all the ", 
+            "dimensions as variable 'pr' in 'exp'.") 
+    }
+    if (!all(dim(exp[['pr']]) %in% dim(pet_exp))) {
+      stop("Parameter 'pet_exp' must have the same length of all the ", 
+           "dimensions as variable 'pr' in 'exp'.")
+    }
+    if (any(names(dim(exp[['pr']])) != names(dim(pet_exp)))) {
+      pos <- match(names(dim(exp[['pr']])), names(dim(pet_exp)))
+      pet_exp <- aperm(pet_exp, pos)
+    }
+    pet[[1]] <- pet_exp
+  } else if (is.null(dates_exp)) {
+    stop("Parameter 'dates_exp' must be provided.")
+  }
+  ## pet_expcor
+  if (!is.null(exp_cor)) {
+    if (!is.null(pet_expcor)) {
+      if (length(dim(exp_cor[['pr']])) != length(dim(pet_expcor))) {
+        stop("Parameter 'pet_expcor' must have the same length of all the ", 
+              "dimensions as variable 'pr' in 'exp_cor'.") 
+      }
+      if (!all(dim(exp_cor[['pr']]) %in% dim(pet_expcor))) {
+        stop("Parameter 'pet_expcor' must have the same length of all the ", 
+            "dimensions as variable 'pr' in 'exp_cor'.")
+      }
+      if (any(names(dim(exp_cor[['pr']])) != names(dim(pet_expcor)))) {
+        pos <- match(names(dim(exp_cor[['pr']])), names(dim(pet_expcor)))
+        pet_expcor <- aperm(pet_expcor, pos)
+      } 
+      pet[[2]] <- pet_expcor
+    } else if (is.null(dates_expcor)) {
+      stop("Parameter 'dates_expcor' must be provided.")
+    }
+  }
+
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!all(sapply(exp, function(x) time_dim %in% names(dim(x))))) {
+    stop("Parameter 'time_dim' is not found in 'exp' dimension.")
+  }
+  if (!is.null(exp_cor)) {
+    if (!all(sapply(exp_cor, function(x) time_dim %in% names(dim(x))))) {
+      stop("Parameter 'time_dim' is not found in 'exp_cor' dimension.")
+    }
+  }
+  ## leadtime_dim
+  if (!is.character(leadtime_dim) | length(leadtime_dim) != 1) {
+    stop("Parameter 'leadtime_dim' must be a character string.")
+  }
+  if (!all(sapply(exp, function(x) leadtime_dim %in% names(dim(x))))) {
+    stop("Parameter 'leadtime_dim' is not found in 'exp' dimension.")
+  }
+  if (!is.null(exp_cor)) {
+    if (!all(sapply(exp_cor, function(x) leadtime_dim %in% names(dim(x))))) {
+      stop("Parameter 'leadtime_dim' is not found in 'exp_cor' dimension.")
+    }
+  }
+  ## memb_dim
+  if (!is.character(memb_dim) | length(memb_dim) != 1) {
+    stop("Parameter 'memb_dim' must be a character string.")
+  }
+  if (!all(sapply(exp, function(x) memb_dim %in% names(dim(x))))) {
+    stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+  }
+  if (!is.null(exp_cor)) {
+    if (!all(sapply(exp_cor, function(x) memb_dim %in% names(dim(x))))) {
+      stop("Parameter 'memb_dim' is not found in 'exp_cor' dimension.")
+    }
+  }
+
+  ## lat_dim
+  if (!is.character(lat_dim) | length(lat_dim) != 1) {
+    stop("Parameter 'lat_dim' must be a character string.")
+  }
+  if (!all(sapply(exp, function(x) lat_dim %in% names(dim(x))))) {
+    stop("Parameter 'lat_dim' is not found in 'exp' dimension.")
+  }
+  if (!is.null(exp_cor)) {
+    if (!all(sapply(exp_cor, function(x) lat_dim %in% names(dim(x))))) {
+      stop("Parameter 'lat_dim' is not found in 'exp_cor' dimension.")
+    }
+  }
+  
+  ## dates
+  if (is.null(dates_exp)) {
+    stop("Parameter 'dates_exp' is missing, dates must be provided.")
+  }
+  if (!(is.Date(dates_exp)) & !(is.POSIXct(dates_exp))) {
+    stop("Parameter 'dates_exp' is not of the correct class, ", 
+         "only 'Date' and 'POSIXct' classes are accepted.")
+  }
+  if (!time_dim %in% names(dim(dates_exp)) | !leadtime_dim %in% names(dim(dates_exp))) {
+    stop("Parameter 'dates_exp' must have 'time_dim' and 'leadtime_dim' ",
+         "dimension.")
+  }
+  if (!all(dim(exp[[1]])[c(time_dim, leadtime_dim)] == 
+           dim(dates_exp)[c(time_dim, leadtime_dim)])) {
+    stop("Parameter 'dates_exp' needs to have the same length as 'time_dim' ", 
+         "and 'leadtime_dim' as 'exp'.")
+  }
+
+  if (!is.null(exp_cor)) {
+    if (is.null(dates_expcor)) {
+      stop("Parameter 'dates_expcor' is missing, dates for 'exp_cor' must be ", 
+           "provided if exp_cor is not NULL.")
+    }
+    if (!(is.Date(dates_expcor)) & !(is.POSIXct(dates_expcor))) {
+      stop("Element 'Dates' in 'exp_cor' is not of the correct class, ", 
+           "only 'Date' and 'POSIXct' classes are accepted.")
+    }
+    if (!time_dim %in% names(dim(dates_expcor)) | !leadtime_dim %in% names(dim(dates_expcor))) {
+      stop("Parameter 'dates_expcor' must have 'time_dim' and 'leadtime_dim' ",
+           "dimension.")
+    }
+    if (!all(dim(exp_cor[[1]])[c(time_dim, leadtime_dim)] == 
+             dim(dates_expcor)[c(time_dim, leadtime_dim)])) {
+      stop("Parameter 'dates_expcor' needs to have the same length as ", 
+           "'time_dim' and 'leadtime_dim' as 'exp_cor'.")
+    }
+  }
+
+  ## accum
+  if (accum > dim(exp[[1]])[leadtime_dim]) {
+    stop(paste0("Cannot compute accumulation of ", accum, " months because ", 
+                "loaded data has only ", dim(exp[[1]])[leadtime_dim], " months."))
+  }
+
+  ## ref_period
+  if (!is.null(ref_period)) {
+    if (length(ref_period) != 2) {
+      warning("Parameter 'ref_period' must be of length two indicating the ",
+              "first and end years of the reference period. It will not ", 
+              "be used.")
+      ref_period <- NULL
+    } else if (!all(sapply(ref_period, is.numeric))) {
+      warning("Parameter 'ref_period' must be a numeric vector indicating the ", 
+              "'start' and 'end' years of the reference period. It will not ",
+              "be used.")
+      ref_period <- NULL
+    } else if (ref_period[[1]] > ref_period[[2]]) {
+      warning("In parameter 'ref_period' 'start' cannot be after 'end'. It ",
+              "will not be used.")
+      ref_period <- NULL
+    } else if (!all(unlist(ref_period) %in% year(dates_exp))) {
+      warning("Parameter 'ref_period' contain years outside the dates. ", 
+              "It will not be used.")
+      ref_period <- NULL
+    } else {
+      years <- year(ClimProjDiags::Subset(dates_exp, along = leadtime_dim, 
+                                          indices = 1))
+      ref_period[[1]] <- which(ref_period[[1]] == years)
+      ref_period[[2]] <- which(ref_period[[2]] == years)
+    }
+  }
+
+  ## standardization
+  if (!is.logical(standardization)) {
+    stop("Parameter 'standardization' must be a logical value.")
+  }
+
+  ## param_error
+  if (!is.numeric(param_error)) {
+    stop("Parameter 'param_error' must be a numeric value.")
+  }
+
+  ## handle_infinity
+  if (!is.logical(handle_infinity)) {
+    stop("Parameter 'handle_infinity' must be a logical value.")
+  }
+
+  ## cross_validation
+  if (!is.logical(cross_validation)) {
+    stop("Parameter 'cross_validation' must be a logical value.")
+  }
+  if (!is.null(exp_cor)) {
+    if (cross_validation & standardization) {
+      warning("Detected 'cross_validation' = TRUE. It will be set as FALSE ", 
+              "since 'exp_cor' is provided.")
+      cross_validation <- FALSE
+    }
+  }
+
+  ## method
+  if (!(method %in% c('parametric', 'non-parametric'))) {
+    stop("Parameter 'method' must be a character string containing one of ", 
+         "the following methods: 'parametric' or 'non-parametric'.")
+  }
+
+  ## distribution
+  if (!(distribution %in% c('log-Logistic', 'Gamma', 'PearsonIII'))) {
+    stop("Parameter 'distribution' must be a character string containing one ", 
+         "of the following distributions: 'log-Logistic', 'Gamma' or ", 
+         "'PearsonIII'.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+
+  ## na.rm
+  if (!is.logical(na.rm)) {
+    stop("Parameter 'na.rm' must be logical.")
+  }
+
+  ## params
+  if (!is.null(params)) {
+    if (!is.numeric(params)) {
+      stop("Parameter 'params' must be numeric.")
+    }
+    if (!all(c(time_dim, leadtime_dim, 'coef') %in% names(dim(params)))) {
+      stop("Parameter 'params' must be a multidimensional array with named ", 
+           "dimensions: '", time_dim, "', '", leadtime_dim, "' and 'coef'.")
+    }
+    if (distribution == "Gamma") {
+      if (dim(params)['coef'] != 2) {
+        stop("For '", distribution, "' distribution, params array should have ",
+             "'coef' dimension of length 2.")
+      }
+    } else {
+      if (dim(params)['coef'] != 3) {
+        stop("For '", distribution, "' distribution, params array should have ",
+             "'coef' dimension of length 3.")
+      }
+    }
+  }
+
+  ## return_params
+  if (!is.logical(return_params)) {
+    stop("Parameter 'return_params' must be logical.")
+  }
+
+  # Complete dates
+  dates <- .return2list(dates_exp, dates_expcor)
+
+  # Compute PeriodSPEI
+  k = 0
+  spei_res <- NULL
+  computed_pet <- FALSE
+
+  for (data in .return2list(exp, exp_cor)) {
+    k = k + 1
+    # Evapotranspiration estimation (unless pet is already provided)
+    if (is.null(pet[[k]]) | computed_pet) {
+      pet[[k]] <- PeriodPET(data = data, dates = dates[[k]], 
+                            lat = lat, pet_method = pet_method[k], 
+                            time_dim = time_dim, leadtime_dim = leadtime_dim, 
+                            lat_dim = lat_dim, na.rm = na.rm, ncores = ncores)
+      computed_pet <- TRUE
+    }
+    if (any(is.na(pet[[k]]))) {
+      mask_na <- which(is.na(pet[[k]]))
+      warning("There are NAs in PET.")
+    }
+
+    # Accumulation
+    diff_p_pet <- data$pr - pet[[k]]
+    dates_monthly <- .datesmask(dates[[k]])
+    accumulated <- .Accumulation(data = diff_p_pet, 
+                                 dates_monthly = dates_monthly, accum = accum, 
+                                 time_dim = time_dim, leadtime_dim = leadtime_dim, 
+                                 ncores = ncores)
+
+    # Standardization
+    if (standardization) {
+      spei <- .Standardization(data = accumulated, params = params, 
+                               accum = accum, time_dim = time_dim,
+                               leadtime_dim = leadtime_dim,
+                               memb_dim = memb_dim, ref_period = ref_period,
+                               cross_validation = cross_validation,
+                               handle_infinity = handle_infinity,
+                               param_error = param_error,
+                               method = method, distribution = distribution, 
+                               na.rm = TRUE, ncores = ncores)
+
+      ref_period <- NULL
+      params <- spei$params
+      spei_res <- spei[[1]]
+    } else {
+      spei_res <- accumulated
+    }
+    pos <- match(names(dim(data[[1]])), names(dim(spei_res)))
+    spei_res <- aperm(spei_res, pos)
+  }
+
+  if (standardization) {
+    if (return_params) {
+      return(list(spei = spei_res, params = params))
+    } else {
+      return(spei_res)
+    }
+  } else {
+    return(spei_res)
+  }
+}
+
+.Standardization <- function(data, params = NULL, accum = 1, time_dim = 'syear', 
+                             leadtime_dim = 'time', memb_dim = 'ensemble', 
+                             ref_period = NULL, cross_validation = FALSE,  
+                             handle_infinity = FALSE, param_error = -9999, 
+                             method = 'parametric', distribution = 'log-Logistic', 
+                             na.rm = FALSE, ncores = NULL) {
+
+  nleadtime <- dim(data)[leadtime_dim]
+  ntime <- dim(data)[time_dim]
+
+  if (!cross_validation) {
+    ntime <- 1
+  }
+
+  coef = switch(distribution,
+                "Gamma" = array(NA, dim = 2, dimnames = list(c('alpha', 'beta'))),
+                "log-Logistic" = array(NA, dim = 3, dimnames = list(c('xi', 'alpha', 'kappa'))),
+                "PearsonIII" = array(NA, dim = 3, dimnames = list(c('mu', 'sigma', 'gamma'))))
+
+  if (is.null(params)) {
+    params <- array(NA, dim = c(ntime, nleadtime, length(coef)))
+    names(dim(params)) <- c(time_dim, leadtime_dim, 'coef')
+  } else if (length(dim(params)) < 2) {
+    params <- array(params, dim = c(length(params), ntime, nleadtime))
+    # dim(params): [time_dim, leadtime_dim, coef] 
+    # with the values repeated each time_dim and leadtime_dim
+    params <- aperm(params, c(2,3,1)) 
+    names(dim(params)) <- c(time_dim, leadtime_dim, 'coef')
+  }
+
+  spei <- Apply(data = list(data = data, params = params), 
+                target_dims = list(data = c(leadtime_dim, time_dim, memb_dim), 
+                                   params = c(time_dim, leadtime_dim, 'coef')),
+                fun = .standardization, 
+                coef = coef, 
+                leadtime_dim = leadtime_dim, 
+                time_dim = time_dim, memb_dim = memb_dim, 
+                ref_period = ref_period, handle_infinity = handle_infinity,
+                cross_validation = cross_validation, param_error = param_error, 
+                method = method, distribution = distribution, 
+                na.rm = na.rm, 
+                output_dims = list(spei = c(leadtime_dim, time_dim, memb_dim), 
+                                   params = c(time_dim, leadtime_dim, 'coef')), 
+                ncores = ncores)
+  return(spei)
+}
+
+.standardization <- function(data, params, coef, leadtime_dim = 'time', 
+                             time_dim = 'syear', memb_dim = 'ensemble', 
+                             ref_period = NULL, handle_infinity = FALSE, 
+                             cross_validation = FALSE, param_error = -9999, 
+                             method = 'parametric', distribution = 'log-Logistic',
+                             na.rm = FALSE) {
+  # data: [leadtime_dim, time_dim, memb_dim]
+  # params: [time_dim, leadtime_dim, 'coef']
+
+  # maximum number of parameters needed to define any of the considered distributions
+  ncoef <- length(coef) 
+  nleadtime <- as.numeric(dim(data)[leadtime_dim])
+  ntime <- as.numeric(dim(data)[time_dim])
+  nmemb <- as.numeric(dim(data)[memb_dim])
+
+  if (cross_validation) {
+    params_result <- array(data = NA, dim = c(ntime, nleadtime, ncoef))
+  } else {
+    params_result <- array(data = NA, dim = c(1, nleadtime, ncoef))
+  }
+  names(dim(params_result)) <- c(time_dim, leadtime_dim, 'coef')
+
+  if (all(is.na(data))) {
+    spei_mod <- array(NA, dim(data))
+    # if the data [time, sdate, memb] has no variability it will raise an error 
+    # further down, so we assign a value to the result and skip the step
+  } else if (anyNA(data) && !na.rm) {
+    spei_mod <- array(NA, dim(data))
+  } else if (var(data, na.rm = T) == 0) {                              
+    spei_mod <- array(param_error, dim(data)) # Add this?
+  } else {
+    if (is.null(ref_period)) {
+      ref.start <- NULL
+      ref.end <- NULL
+    } else {
+      ref.start <- ref_period[[1]]
+      ref.end <- ref_period[[2]]
+    }
+
+    spei_mod <- array(data = NA, dim = c(nleadtime, ntime, nmemb))
+    names(dim(spei_mod)) <- c(leadtime_dim, time_dim, memb_dim)
+
+    for (ff in 1:nleadtime) {
+      data_subset <- ClimProjDiags::Subset(data, along = leadtime_dim, 
+                                           indices = ff, drop = 'selected')
+      params_tmp <- if (all(is.na(params))) {NULL} else {params[, ff, ]}
+
+      spei_data <- .std(data = data_subset, coef = coef, ntime = ntime, 
+                        nmemb = nmemb, method = method, 
+                        distribution = distribution, na.rm = na.rm, 
+                        ref.start = ref.start, ref.end = ref.end, 
+                        params = params_tmp, handle_infinity = handle_infinity, 
+                        cross_validation = cross_validation)
+
+      spei_mod[ff, , ] <- spei_data[[1]]
+      params_ff <- spei_data[[2]]
+      # lengthen dimension coef of params_ff in case it doesn't match the 
+      # corresponding dimension of parms_months
+      if (!is.null(params_ff)) {
+        if (length(params_ff) < ncoef) {
+          params_ff <- append(params_ff, array(NA, dim = ncoef - length(params_ff)))
+        }
+        params_result[, ff, ] <- params_ff
+      }
+    }
+  }
+  return(list(spei = spei_mod, params = params_result))
+}
+
+.std <- function(data, coef, ntime, nmemb, method = 'parametric', 
+                 distribution = 'log-Logistic', na.rm = FALSE, 
+                 ref.start = NULL, ref.end = NULL, params = NULL, 
+                 handle_infinity = FALSE, cross_validation = FALSE) {
+  
+  # data: [time_dim, memb_dim]
+  # params: NULL or [(ntime), coef]
+
+  fit = 'ub-pwm' # hard-coded
+
+  if (method == 'non-parametric') { 
+    bp <- matrix(0, length(data), 1)
+    for (i in 1:length(data)) {
+      bp[i,1] = sum(data[] <= data[i], na.rm = na.rm);  # Writes the rank of the data
+    }
+    std_index <- qnorm((bp - 0.44)/(length(data) + 0.12))
+    dim(std_index) <- c(ntime, nmemb)
+    # it won't return params to be used in exp_cor; also it is not using 
+    # handle_infinity nor cross_validation
+    params_result <- NULL
+  } else {
+    std_index <- array(NA, c(ntime, nmemb))
+    # Select window if necessary
+    if (!is.null(ref.start) && !is.null(ref.end)) {
+      data.fit <- window(data, ref.start, ref.end)	
+    } else {
+      data.fit <- data
+    }
+
+    if (cross_validation) {  
+      loop_years <- ntime
+    } else {
+      loop_years <- 1
+    }
+
+    params_result <- array(NA, dim = c(loop_years, length(coef)))
+    colnames(params_result) <- names(coef)
+    
+    for (nsd in 1:loop_years) { # loop over years (in case of cross_validation)
+      # Cumulative series (acu)
+      if (cross_validation) {
+        acu <- as.vector(data.fit[-nsd, ])
+      } else {
+        acu <- as.vector(data.fit)
+      }
+
+      acu.sorted <- sort.default(acu, method = "quick")
+      # remove NAs (no need if(na.rm) because if there are NA and na.rm = F 
+      # we don't get to this point)
+      acu.sorted <- acu.sorted[!is.na(acu.sorted)] 
+      # else all acu was NA and we don't need to continue with this case   
+      
+      if (length(acu.sorted) != 0) { 
+        acu_sd <- sd(acu.sorted)
+        if (!is.na(acu_sd)) {
+          if (acu_sd != 0) {
+            if (distribution != "log-Logistic") {
+              pze <- sum(acu == 0) / length(acu)
+              acu.sorted = acu.sorted[acu.sorted > 0]
+            }
+            if (!is.null(params)) {
+              f_params <- as.vector(params)
+              params_result[nsd, ] <- f_params
+            } else {
+              # else coef will be NA
+              if (length(acu.sorted) >= 4) { 
+                # Calculate probability weighted moments based on fit with lmomco or TLMoments
+                pwm = switch(fit,
+                             "pp-pwm" = pwm.pp(acu.sorted, -0.35, 0, nmom = 3),
+                             pwm.ub(acu.sorted, nmom = 3)
+                             # TLMoments::PWM(acu.sorted, order = 0:2)
+                             )
+
+                # Check L-moments validity
+                lmom <- pwm2lmom(pwm)          
+                if (!(!are.lmom.valid(lmom) || anyNA(lmom[[1]]) || any(is.nan(lmom[[1]])))) {
+                  # lmom fortran functions need specific inputs L1, L2, T3
+                  # this is handled by lmomco internally with lmorph
+                  fortran_vec = c(lmom$lambdas[1:2], lmom$ratios[3])
+                  # Calculate parameters based on distribution with lmom then lmomco
+                  f_params = switch(distribution,
+                                    "log-Logistic" = tryCatch(lmom::pelglo(fortran_vec), 
+                                                              error = function(e){parglo(lmom)$para}),
+                                    "Gamma" = tryCatch(lmom::pelgam(fortran_vec), 
+                                                       error = function(e){pargam(lmom)$para}),
+                                    "PearsonIII" = tryCatch(lmom::pelpe3(fortran_vec), 
+                                                            error = function(e){parpe3(lmom)$para}))
+                  # Adjust if user chose log-Logistic and max-lik
+                  if (distribution == 'log-Logistic' && fit == 'max-lik') {
+                    f_params = parglo.maxlik(acu.sorted, f_params)$para
+                  }
+                  params_result[nsd, ] <- f_params
+                } # end for the case the L-moments are not valid (std_index will be NA)
+              } # end for case there are not enough values to estimate the parameters (std_index will be NA)
+            } # end estimation of f_param
+            # Calculate cdf based on distribution with lmom
+            if (all(is.na(params_result[nsd,]))) {
+              cdf_res <- NA
+            } else {
+              f_params <- params_result[nsd,]
+              f_params <- f_params[which(!is.na(f_params))]
+              cdf_res = switch(distribution,
+                               "log-Logistic" = lmom::cdfglo(data, f_params),
+                               "Gamma" = lmom::cdfgam(data, f_params),
+                               "PearsonIII" = lmom::cdfpe3(data, f_params))
+            }
+
+            std_index_cv <- array(qnorm(cdf_res), dim = c(ntime, nmemb))
+    
+            # Adjust if user chose Gamma or PearsonIII - Not tested: For future development
+            # if (distribution != 'log-Logistic') { 
+            #   std_index[ff,s] = qnorm(pze + (1-pze)*pnorm(std_index[ff,s])) # ff doesn't exist at this point
+            # }
+            if (cross_validation) {
+              std_index[nsd, ] <- std_index_cv[nsd, ]  
+            } else {
+              std_index <- std_index_cv
+            }
+          }
+        } # end if for the case there is no variability        
+      } # end if for the case all NA in acu
+    } # next year (in case of cross_validation or all done if cross_validation == F)
+
+    if (handle_infinity) { # could also use "param_error" ?; we are giving it the min/max value of the grid point
+      std_index[is.infinite(std_index) & std_index < 0] <- min(std_index[!is.infinite(std_index)])
+      std_index[is.infinite(std_index) & std_index > 0] <- max(std_index[!is.infinite(std_index)]) 
+    }
+    # f_params will be params only if cross_validation is FALSE
+    # (otherwise there will be one f_params per year;
+    # but the output params will be read only in the case that
+    # it is called with cross_validation FALSE)
+  }
+  return(list(std_index, params_result)) 
+}
+
+.Accumulation <- function(data, dates_monthly, accum = 2, 
+                          time_dim = 'syear', leadtime_dim = 'time', 
+                          ncores = NULL) {
+  
+  accumulated <- Apply(data = list(data), 
+                       target_dims = list(data = c(leadtime_dim, time_dim)), 
+                       dates_monthly = dates_monthly,
+                       accum = accum,
+                       output_dims = c(leadtime_dim, time_dim), 
+                       leadtime_dim = leadtime_dim, time_dim = time_dim,
+                       fun = .accumulation, 
+                       ncores = ncores)$output1
+
+  pos <- match(names(dim(accumulated)), names(dim(data)))
+  accumulated <- aperm(accumulated, pos)
+
+  return(accumulated)
+
+}
+
+.accumulation <- function(data, dates_monthly, accum = 1, 
+                          time_dim = 'syear', leadtime_dim = 'time') {
+
+  # data:[time, syear]
+  dims <- dim(data)
+  
+  data_vector <- array(NA, dim = length(dates_monthly))
+  count <- 1
+  for (dd in 1:length(dates_monthly)) {
+    if (dates_monthly[dd] == 1) {
+      data_vector[dd] <- as.vector(data)[count]
+      count <- count + 1
+    }
+  }
+
+  data_sum_x <- rollapply(data_vector, accum, sum) 
+  # adds as many NAs as needed at the begining to account for the months that 
+  # cannot be added (depends on accu) and so that the position in the vector 
+  # corresponds to the accumulated of the previous months (instead of the 
+  # accumulated of the next months)
+  data_sum_x <- c(rep(NA, accum-1), data_sum_x) 
+  # discard the months that don't appear in the original data
+  data_sum_x <- data_sum_x[which(dates_monthly == 1)] 
+
+  accum_result <- array(data_sum_x, dim = c(dims))
+  return(accum_result)
+}
+
+.datesmask <- function(dates) {
+  ini <- as.Date(paste(lubridate::year(min(dates)), 01, 01, sep = '-'))
+  end <- as.Date(paste(lubridate::year(max(dates)), 12, 31, sep = '-'))
+  daily <- as.Date(ini:end)
+  monthly <- daily[which(lubridate::day(daily) == 1)]
+  dates_mask <- array(0, dim = length(monthly))
+  for (dd in 1:length(dates)) {
+    ii <- which(monthly == as.Date(paste(lubridate::year(dates[dd]), 
+                                         lubridate::month(dates[dd]), 
+                                         01, sep = '-')))
+    dates_mask[ii] <- 1
+  }
+  return(dates_mask)
+}
\ No newline at end of file
diff --git a/R/zzz.R b/R/zzz.R
index cf9163970f76c2da74b3d811c1fa0b71445beeab..47d871d4d8d4f2a3d1dfe01b9a79c5aea8674e70 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -88,4 +88,22 @@ wind2CF <- function(wind, pc) {
   power <- wind2power(wind, pc)
 	CF <- power / pc$attr$RatedPower
 	return(CF)
+}
+
+.KnownLonNames <- function() {
+  known_lon_names <- c('lon', 'lons', 'longitude', 'x', 'i', 'nav_lon')
+}
+
+.KnownLatNames <- function() {
+  known_lat_names <- c('lat', 'lats', 'latitude', 'y', 'j', 'nav_lat')
+}
+
+.return2list <- function(data1, data2 = NULL) {
+  if (is.null(data1) & is.null(data2)) {
+    return(NULL)
+  } else if (is.null(data2)) {
+    return(list(data1))
+  } else {
+    return(list(data1, data2))
+  }
 }
\ No newline at end of file
diff --git a/man/CST_PeriodPET.Rd b/man/CST_PeriodPET.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..10383afcef72d310a54aae7dd30493237d029005
--- /dev/null
+++ b/man/CST_PeriodPET.Rd
@@ -0,0 +1,85 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodPET.R
+\name{CST_PeriodPET}
+\alias{CST_PeriodPET}
+\title{Compute the Potential Evapotranspiration}
+\usage{
+CST_PeriodPET(
+  data,
+  pet_method = "hargreaves",
+  time_dim = "syear",
+  leadtime_dim = "time",
+  lat_dim = "latitude",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A named list with the needed \code{s2dv_cube} objects containing 
+the seasonal forecast experiment in the data element for each variable. 
+Specific variables are needed for each method used in computing the 
+Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for 
+'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' 
+'tmean' is required. The units for temperature variables
+('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+precipitation ('pr') need to be in mm/month.}
+
+\item{pet_method}{A character string indicating the method used to compute 
+the potential evapotranspiration. The accepted methods are:
+'hargreaves' and 'hargreaves_modified', that require the data to have 
+variables tmin and tmax; and 'thornthwaite', that requires variable 
+'tmean'.}
+
+\item{time_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'syear'.}
+
+\item{leadtime_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'time'.}
+
+\item{lat_dim}{A character string indicating the name of the latitudinal 
+dimension. By default it is set by 'latitude'.}
+
+\item{na.rm}{A logical value indicating whether NA values should be removed 
+from data. It is FALSE by default.}
+
+\item{ncores}{An integer value indicating the number of cores to use in 
+parallel computation.}
+
+\item{dates}{An array of temporal dimensions containing the Dates of 
+'data'. It must be of class 'Date' or 'POSIXct'.}
+
+\item{lat}{A numeric vector containing the latitude values of 'data'.}
+}
+\description{
+Compute the Potential evapotranspiration (PET) that is the amount of 
+evaporation and transpiration that would occur if a sufficient water source 
+were available. This function calculate PET according to the Thornthwaite, 
+Hargreaves or Hargreaves-modified equations.
+}
+\details{
+This function is build to work be compatible with other tools in 
+that work with 's2dv_cube' object class. The input data must be this object 
+class. If you don't work with 's2dv_cube', see PeriodPET. For more information 
+on the SPEI calculation, see functions CST_PeriodStandardization and 
+CST_PeriodAccumulation.
+}
+\examples{
+dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1)
+
+exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+end_year <- 2012
+dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), 
+                         paste0(2010:end_year, "-09-15"), 
+                         paste0(2010:end_year, "-10-16")), "UTC")
+dim(dates_exp) <- c(syear = 3, time = 3)
+
+lat <- c(40)
+
+exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr)
+
+res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp)
+
+}
diff --git a/man/CST_PeriodSPEI.Rd b/man/CST_PeriodSPEI.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1500a648b3af77766ec6ad85ec3fbad343f24f00
--- /dev/null
+++ b/man/CST_PeriodSPEI.Rd
@@ -0,0 +1,212 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodSPEI.R
+\name{CST_PeriodSPEI}
+\alias{CST_PeriodSPEI}
+\title{Compute the Standardised Precipitation-Evapotranspiration Index}
+\usage{
+CST_PeriodSPEI(
+  exp,
+  exp_cor = NULL,
+  time_dim = "syear",
+  leadtime_dim = "time",
+  memb_dim = "ensemble",
+  lat_dim = "latitude",
+  accum = 1,
+  ref_period = NULL,
+  params = NULL,
+  pet_exp = NULL,
+  pet_expcor = NULL,
+  standardization = TRUE,
+  cross_validation = FALSE,
+  pet_method = "hargreaves",
+  method = "parametric",
+  distribution = "log-Logistic",
+  param_error = -9999,
+  handle_infinity = FALSE,
+  return_params = FALSE,
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named list with the needed \code{s2dv_cube} objects containing 
+the seasonal forecast experiment in the data element for each variable. 
+Specific variables are needed for each method used in computing the 
+Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and 
+'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. 
+Variable 'pr' is always needed. The units for temperature variables
+('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+precipitation ('pr') need to be in mm/month.}
+
+\item{exp_cor}{A named list with the needed \code{s2dv_cube} objects for each 
+variable in which the quantile PeriodSPEI should be applied. If it is not 
+specified, the PeriodSPEI is calculated from object 'exp'.}
+
+\item{time_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'syear'.}
+
+\item{leadtime_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'time'.}
+
+\item{memb_dim}{A character string indicating the name of the dimension in 
+which the ensemble members are stored. When set it to NULL, threshold is 
+computed for individual members.}
+
+\item{lat_dim}{A character string indicating the name of the latitudinal 
+dimension. By default it is set by 'latitude'.}
+
+\item{accum}{An integer value indicating the number of months for the 
+accumulation for each variable. When it is greater than 1, the result will 
+be filled with NA until the accum time_dim dimension number due to the 
+accumulation to previous months.}
+
+\item{ref_period}{A list with two numeric values with the starting and end 
+points of the reference period used for computing the index. The default 
+value is NULL indicating that the first and end values in data will be 
+used as starting and end points.}
+
+\item{params}{A multidimensional array with named dimensions for computing the 
+SPEI. This option overrides computation of fitting parameters. It needs 
+to have the following dimensions: same leadtime dimension of 'exp' 
+(specified in 'leadtime_dim'); time dimension of length 1 (specified in 
+'time_dim'); and a dimension named 'coef' with the length of the 
+coefficients needed for the used distribution (for 'Gamma' coef dimension is 
+of lenght 2, for 'log-Logistic' or 'PearsonIII' is of length 3). It can't 
+have member dimension (specified in 'memb_dim').}
+
+\item{pet_exp}{A multidimensional array containing the Potential
+EvapoTranspiration data of 'exp'. It must have the same dimensions of the 
+variable 'pr' of 'exp'. If it is NULL it is calculated using the provided 
+variables with specified 'pet_method'. It is NULL by default.}
+
+\item{pet_expcor}{A multidimensional array containing the Potential
+EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of 
+the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the 
+provided variables with specified 'pet_method'. It is NULL by default.}
+
+\item{standardization}{A logical value indicating wether the standardization 
+is computed.}
+
+\item{cross_validation}{A logical value indicating if cross validation is 
+done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and 
+is not provided. It is FALSE by default.}
+
+\item{pet_method}{A character string indicating the method used to compute 
+the potential evapotranspiration. The accepted methods are:
+'hargreaves' and 'hargreaves_modified', that require the data to have 
+variables tmin and tmax; and 'thornthwaite', that requires variable 
+'tmean'.}
+
+\item{method}{A character string indicating the standardization method used. 
+If can be: 'parametric' or 'non-parametric'. It is set to 'parametric' by 
+default.}
+
+\item{distribution}{A character string indicating the name of the distribution 
+function to be used for computing the SPEI. The accepted names are: 
+'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The 
+'Gamma' method only works when only precipitation is provided and other 
+variables are 0 because it is positive defined (SPI indicator).}
+
+\item{param_error}{A numeric value with the error accepted.}
+
+\item{handle_infinity}{A logical value wether to return Infinite values (TRUE)
+or not (FALSE).}
+
+\item{na.rm}{A logical value indicating whether NA values should be removed 
+from data. It is FALSE by default. If it is FALSE and there are NA values, 
+(if standardization is TRUE) all values of other dimensions except time_dim 
+and leadtime_dim will be set to NA directly. On the other hand, if it is 
+TRUE, if the data from other dimensions except time_dim and leadtime_dim is 
+not reaching 4 values, it is not enough values to estimate the parameters 
+and the result will include NA.}
+
+\item{ncores}{An integer value indicating the number of cores to use in 
+parallel computation.}
+}
+\value{
+An 's2dv_cube' object containing the SPEI multidimensional array in 
+element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' 
+will be provided. The parameters of the standardization will only be returned 
+if 'return_params' is TRUE. The SPEI will only be computed if 
+'standardization' is TRUE. If 'standardization' is FALSE, only the climatic 
+water balance (precipitation minus evapotranspiration) will be returned. The 
+resultant arrays will have the same dimensions as the initial input data. The 
+other elements in the 's2dv_cube' will be updated with the combined 
+information of the input data arrays.
+}
+\description{
+Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) 
+that is a multiscalar drought index based on climatic data. It can be used for 
+determining the onset, duration and magnitude of drought conditions with 
+respect to normal conditions in a variety of natural and managed systems such 
+as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI 
+is to compare the highest possible evapotranspiration with the current water 
+availability. The SPEI uses the monthly (or weekly) difference between 
+precipitation and potential evapotranspiration. This represents a simple 
+climatic water balance which is calculated at different time scales to obtain 
+the SPEI. This function is build to work be compatible with other tools in 
+that work with 's2dv_cube' object class. The input data must be this object 
+class. If you don't work with 's2dv_cube', see PeriodSPEI.
+}
+\details{
+Next, some specifications for the calculation of this indicator will be 
+discussed. On the one hand, the model to be used to calculate potential 
+evapotranspiration is specified with the pet_method parameter (hargreaves, 
+hargraves modified or thornwhite). On the other hand, to choose the time scale 
+in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the 
+accum parameter, where you must indicate the number of time steps you want to 
+accumulate throughout leadtime_dim. Since the accumulation is done for the 
+elapsed time steps, there will be no complete accumulations until reaching the 
+time instant equal to the value of the parameter. For this reason, in the 
+result, we will find that for the dimension where the accumulation has been 
+carried out, the values of the array will be NA since they do not include 
+complete accumulations. Also, there is a parameter to specify if the 
+standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a 
+probability distribution to transform the original values to standardized 
+units that are comparable in space and time and at different SPEI time scales. 
+The na.rm parameter is a logical parameter used to decide whether to remove 
+the NA values from the data before doing the calculation. It must be taken 
+into account that if na.rm == FALSE and there is some NA value in the specific 
+coordinates which the SPEI is computed, standardization cannot be carried out 
+for those coordinates and therefore, the result will be filled with NA for the 
+specific coordinates. However, when na.rm == TRUE,  if the amount of data for 
+those specific coordinates is smaller than 4, it will not be possible to carry 
+out because we will not have enough data and the result will be also filled 
+with NAs for that coordinates.
+}
+\examples{
+dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, 
+         latitude = 2, longitude = 1, ensemble = 25)
+dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, 
+            latitude = 2, longitude = 1, ensemble = 15)
+
+dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), 
+                         paste0(2010:2015, "-09-15"), 
+                         paste0(2010:2015, "-10-16")), 'UTC')
+dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3)
+dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), 
+                            paste0(2020, "-10-16")), 'UTC')
+dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3)
+lat <- c(40,40.1)
+
+exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims)
+exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims)
+exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims)
+
+expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor)
+expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor)
+expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor)
+
+exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 
+               'pr' = expcor_pr)
+
+exp <- lapply(exp, CSTools::s2dv_cube, coords = list(latitude = lat), 
+             Dates = dates_exp)
+exp_cor <- lapply(exp_cor, CSTools::s2dv_cube, coords = list(latitude = lat), 
+                 Dates = dates_expcor)
+
+res <- CST_PeriodSPEI(exp = exp, exp_cor = exp_cor)
+
+}
diff --git a/man/PeriodPET.Rd b/man/PeriodPET.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..7ca1073429f35717fe6a3be22bbf7cffe95e46b0
--- /dev/null
+++ b/man/PeriodPET.Rd
@@ -0,0 +1,84 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodPET.R
+\name{PeriodPET}
+\alias{PeriodPET}
+\title{Compute the Potential Evapotranspiration}
+\usage{
+PeriodPET(
+  data,
+  dates,
+  lat,
+  pet_method = "hargreaves",
+  time_dim = "syear",
+  leadtime_dim = "time",
+  lat_dim = "latitude",
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{data}{A named list with the needed \code{s2dv_cube} objects containing 
+the seasonal forecast experiment in the data element for each variable. 
+Specific variables are needed for each method used in computing the 
+Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+variables for each method are: for 'hargreaves', 'tmin' and 'tmax'; for 
+'hargreaves_modified' are 'tmin', 'tmax' and 'pr'; for method 'thornthwaite' 
+'tmean' is required. The units for temperature variables
+('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+precipitation ('pr') need to be in mm/month.}
+
+\item{dates}{An array of temporal dimensions containing the Dates of 
+'data'. It must be of class 'Date' or 'POSIXct'.}
+
+\item{lat}{A numeric vector containing the latitude values of 'data'.}
+
+\item{pet_method}{A character string indicating the method used to compute 
+the potential evapotranspiration. The accepted methods are:
+'hargreaves' and 'hargreaves_modified', that require the data to have 
+variables tmin and tmax; and 'thornthwaite', that requires variable 
+'tmean'.}
+
+\item{time_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'syear'.}
+
+\item{leadtime_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'time'.}
+
+\item{lat_dim}{A character string indicating the name of the latitudinal 
+dimension. By default it is set by 'latitude'.}
+
+\item{na.rm}{A logical value indicating whether NA values should be removed 
+from data. It is FALSE by default.}
+
+\item{ncores}{An integer value indicating the number of cores to use in 
+parallel computation.}
+}
+\description{
+Compute the Potential evapotranspiration (PET) that is the amount of 
+evaporation and transpiration that would occur if a sufficient water source 
+were available. This function calculate PET according to the Thornthwaite, 
+Hargreaves or Hargreaves-modified equations.
+}
+\details{
+For more information on the SPEI calculation, see functions 
+PeriodStandardization and PeriodAccumulation.
+}
+\examples{
+dims <- c(time = 3, syear = 3, ensemble = 1, latitude = 1)
+
+exp_tasmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+exp_tasmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+exp_prlr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+end_year <- 2012
+dates_exp <- as.POSIXct(c(paste0(2010:end_year, "-08-16"), 
+                         paste0(2010:end_year, "-09-15"), 
+                         paste0(2010:end_year, "-10-16")), "UTC")
+dim(dates_exp) <- c(syear = 3, time = 3)
+
+lat <- c(40)
+
+exp1 <- list('tmax' = exp_tasmax, 'tmin' = exp_tasmin, 'pr' = exp_prlr)
+
+res <- PeriodPET(data = exp1, lat = lat, dates = dates_exp)
+
+}
diff --git a/man/PeriodSPEI.Rd b/man/PeriodSPEI.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..049e084cd301704f4d157e77b6c96ee64feb4f3e
--- /dev/null
+++ b/man/PeriodSPEI.Rd
@@ -0,0 +1,215 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PeriodSPEI.R
+\name{PeriodSPEI}
+\alias{PeriodSPEI}
+\title{Compute the Standardised Precipitation-Evapotranspiration Index}
+\usage{
+PeriodSPEI(
+  exp,
+  dates_exp,
+  lat,
+  exp_cor = NULL,
+  dates_expcor = NULL,
+  time_dim = "syear",
+  leadtime_dim = "time",
+  memb_dim = "ensemble",
+  lat_dim = "latitude",
+  accum = 1,
+  ref_period = NULL,
+  params = NULL,
+  pet_exp = NULL,
+  pet_expcor = NULL,
+  standardization = TRUE,
+  cross_validation = FALSE,
+  pet_method = "hargreaves",
+  method = "parametric",
+  distribution = "log-Logistic",
+  param_error = -9999,
+  handle_infinity = FALSE,
+  return_params = FALSE,
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named list with multidimensional array objects containing 
+the seasonal forecast experiment in the data element for each variable. 
+Specific variables are needed for each method used in computing the 
+Potential Evapotranspiration. See parameter 'pet_method'. The accepted 
+variables are: 'tmin' and 'tmax', required for methods 'hargreaves' and 
+'hargreaves_modified' and 'tmean', required for method 'thornthwaite'. 
+Variable 'pr' is always needed. The units for temperature variables
+('tmin', 'tmax' and 'tmean') need to be in Celcius degrees; the units for 
+precipitation ('pr') need to be in mm/month.}
+
+\item{dates_exp}{An array of temporal dimensions containing the Dates of 
+'exp'. It must be of class 'Date' or 'POSIXct'.}
+
+\item{lat}{A numeric vector containing the latitude values of 'exp'.}
+
+\item{exp_cor}{A named list with multidimensional array objects for each 
+variable in which the quantile PeriodSPEI should be applied. If it is not 
+specified, the PeriodSPEI is calculated from object 'exp'.}
+
+\item{dates_expcor}{An array of temporal dimensions containing the Dates of 
+'exp_cor'. It must be of class 'Date' or 'POSIXct'.}
+
+\item{time_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'syear'.}
+
+\item{leadtime_dim}{A character string indicating the name of the temporal 
+dimension. By default, it is set to 'time'.}
+
+\item{memb_dim}{A character string indicating the name of the dimension in 
+which the ensemble members are stored. When set it to NULL, threshold is 
+computed for individual members.}
+
+\item{lat_dim}{A character string indicating the name of the latitudinal 
+dimension. By default it is set by 'latitude'.}
+
+\item{accum}{accum An integer value indicating the number of months for the 
+accumulation for each variable. When it is greater than 1, the result will 
+be filled with NA until the accum time_dim dimension number due to the 
+accumulation to previous months.}
+
+\item{ref_period}{A list with two numeric values with the starting and end 
+points of the reference period used for computing the index. The default 
+value is NULL indicating that the first and end values in data will be 
+used as starting and end points.}
+
+\item{params}{A multidimensional array with named dimensions for computing the 
+SPEI. This option overrides computation of fitting parameters. It needs 
+to be of same time dimension (specified in 'time_dim') of 'exp' and a 
+dimension named 'coef' with the length of the coefficients needed for the 
+used distribution (for 'Gamma' coef dimension is of lenght 2, for 
+'log-Logistic' or 'PearsonIII' is of length 3). It also needs to have a 
+leadtime dimension (specified in 'leadtime_dim') of length 1.}
+
+\item{pet_exp}{A multidimensional array containing the Potential
+EvapoTranspiration data of 'exp'. It must have the same dimensions of the 
+variable 'pr' of 'exp'. If it is NULL it is calculated using the provided 
+variables with specified 'pet_method'. It is NULL by default.}
+
+\item{pet_expcor}{A multidimensional array containing the Potential
+EvapoTranspiration data of 'exp_cor'. It must have the same dimensions of 
+the variable 'pr' of 'exp_cor'. If it is NULL it is calculated using the 
+provided variables with specified 'pet_method'. It is NULL by default.}
+
+\item{standardization}{A logical value indicating wether the standardization 
+is computed.}
+
+\item{cross_validation}{A logical value indicating if cross validation is 
+done (TRUE), or not (FALSE). It only will be used when 'exp_cor' and 
+is not provided. It is FALSE by default.}
+
+\item{pet_method}{A character string indicating the method used to compute 
+the potential evapotranspiration. The accepted methods are:
+'hargreaves' and 'hargreaves_modified', that require the data to have 
+variables tmin and tmax; and 'thornthwaite', that requires variable 
+'tmean'.}
+
+\item{method}{A character string indicating the standardization method used. 
+If can be: 'parametric' or 'non-parametric'.}
+
+\item{distribution}{A character string indicating the name of the distribution 
+function to be used for computing the SPEI. The accepted names are: 
+'log-Logistic' and 'Gamma'. It is set to 'log-Logistic' by default. The 
+'Gamma' method only works when only precipitation is provided and other 
+variables are 0 because it is positive defined (SPI indicator).}
+
+\item{param_error}{A numeric value with the error accepted.}
+
+\item{handle_infinity}{A logical value wether to return Infinite values (TRUE)
+or not (FALSE).}
+
+\item{return_params}{A logical value indicating wether to return parameters 
+array (TRUE) or not (FALSE). It is FALSE by default.}
+
+\item{na.rm}{A logical value indicating whether NA values should be removed 
+from data. It is FALSE by default. If it is FALSE and there are NA values, 
+(if standardization is TRUE) all values of other dimensions except time_dim 
+and leadtime_dim will be set to NA directly. On the other hand, if it is 
+TRUE, if the data from other dimensions except time_dim and leadtime_dim is 
+not reaching 4 values, it is not enough values to estimate the parameters 
+and the result will include NA.}
+
+\item{ncores}{An integer value indicating the number of cores to use in 
+parallel computation.}
+}
+\value{
+An 's2dv_cube' object containing the SPEI multidimensional array in 
+element \code{data}. If 'exp_cor' is provided, only results from 'exp_cor' 
+will be provided. The parameters of the standardization will only be returned 
+if 'return_params' is TRUE. The SPEI will only be computed if 
+'standardization' is TRUE. If 'standardization' is FALSE, only the climatic 
+water balance (precipitation minus evapotranspiration) will be returned. The 
+resultant arrays will have the same dimensions as the initial input data.
+}
+\description{
+Calculation of the Standardised Precipitation-Evapotranspiration Index (SPEI) 
+that is a multiscalar drought index based on climatic data. It can be used for 
+determining the onset, duration and magnitude of drought conditions with 
+respect to normal conditions in a variety of natural and managed systems such 
+as crops, ecosystems, rivers, water resources, etc. The idea behind the SPEI 
+is to compare the highest possible evapotranspiration with the current water 
+availability. The SPEI uses the monthly (or weekly) difference between 
+precipitation and potential evapotranspiration. This represents a simple 
+climatic water balance which is calculated at different time scales to obtain 
+the SPEI.
+}
+\details{
+Next, some specifications for the calculation of this indicator will be 
+discussed. On the one hand, the model to be used to calculate potential 
+evapotranspiration is specified with the pet_method parameter (hargreaves, 
+hargraves modified or thornwhite). On the other hand, to choose the time scale 
+in which you want to accumulate the SPEI (SPEI3, SPEI6...) is done using the 
+accum parameter, where you must indicate the number of time steps you want to 
+accumulate throughout leadtime_dim. Since the accumulation is done for the 
+elapsed time steps, there will be no complete accumulations until reaching the 
+time instant equal to the value of the parameter. For this reason, in the 
+result, we will find that for the dimension where the accumulation has been 
+carried out, the values of the array will be NA since they do not include 
+complete accumulations. Also, there is a parameter to specify if the 
+standardization that can be TRUE or FALSE. If it is TRUE, the data is fit to a 
+probability distribution to transform the original values to standardized 
+units that are comparable in space and time and at different SPEI time scales. 
+The na.rm parameter is a logical parameter used to decide whether to remove 
+the NA values from the data before doing the calculation. It must be taken 
+into account that if na.rm == FALSE and there is some NA value in the specific 
+coordinates which the SPEI is computed, standardization cannot be carried out 
+for those coordinates and therefore, the result will be filled with NA for the 
+specific coordinates. However, when na.rm == TRUE,  if the amount of data for 
+those specific coordinates is smaller than 4, it will not be possible to carry 
+out because we will not have enough data and the result will be also filled 
+with NAs for that coordinates.
+}
+\examples{
+dims <- c(var = 1, sday = 1, sweek = 1, syear = 6, time = 3, 
+         latitude = 2, longitude = 1, ensemble = 25)
+dimscor <- c(var = 1, sday = 1, sweek = 1, syear = 1, time = 3, 
+            latitude = 2, longitude = 1, ensemble = 15)
+exp_tmax <- array(rnorm(900, 27.73, 5.26), dim = dims)
+exp_tmin <- array(rnorm(900, 14.83, 3.86), dim = dims)
+exp_pr <- array(rnorm(900, 21.19, 25.64), dim = dims)
+
+expcor_tmax <- array(rnorm(540, 29.03, 5.67), dim = dimscor)
+expcor_tmin <- array(rnorm(540, 15.70, 4.40), dim = dimscor)
+expcor_pr <- array(rnorm(540, 15.62, 21.38), dim = dimscor)
+
+dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), 
+                      paste0(2010:2015, "-09-15"), 
+                      paste0(2010:2015, "-10-16")))
+dim(dates_exp) <- c(sday = 1, sweek = 1, syear = 6, time = 3)
+dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), 
+                            paste0(2020, "-09-15"), 
+                            paste0(2020, "-10-16")), 'UTC')
+dim(dates_expcor) <- c(sday = 1, sweek = 1, syear = 1, time = 3)
+lat <- c(40,40.1)
+
+exp <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+exp_cor <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 
+               'pr' = expcor_pr)
+res <- PeriodSPEI(exp = exp, exp_cor = exp_cor, lat = lat, 
+                 dates_exp = dates_exp, dates_expcor = dates_expcor)
+
+}
diff --git a/tests/testthat/test-PeriodPET.R b/tests/testthat/test-PeriodPET.R
new file mode 100644
index 0000000000000000000000000000000000000000..647d1e08727046185f58b1516d83e2ea920eda6c
--- /dev/null
+++ b/tests/testthat/test-PeriodPET.R
@@ -0,0 +1,197 @@
+##############################################
+
+# cube1
+cube1 <- NULL
+cube1$data <- array(rnorm(10), dim = c(lat = 2, time = 5))
+class(cube1) <- 's2dv_cube'
+
+# cube2
+cube2 <- NULL
+cube2$data <- array(rnorm(10), dim = c(lat = 2, time = 5))
+class(cube2) <- 's2dv_cube'
+cube2$coords <- list(lat = 1:2)
+
+# dat1
+dims <- c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10)
+
+set.seed(1)
+exp_tmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+set.seed(2)
+exp_tmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+set.seed(3)
+exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+
+dates_exp <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), 
+                       paste0(2010:2015, "-10-16")))
+dim(dates_exp) <- c(syear = 6, time = 3)
+
+lat <- c(40,40.1)
+
+exp1 <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+
+# dat2
+dims2 <- c(styear = 6, ftime = 3, lat = 2, lon = 1, member = 10)
+
+set.seed(1)
+exp_tmean <- array(rnorm(100, 17.34, 9.18), dim = dims2)
+set.seed(2)
+exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims2)
+
+dates_exp2 <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), 
+                       paste0(2010:2015, "-10-16")))
+dim(dates_exp2) <- c(sday = 1, sweek = 1, styear = 6, ftime = 3)
+
+lat <- c(40,40.1)
+
+exp2 <- list('tmean' = exp_tmean, 'pr' = exp_pr)
+
+# cube4
+cube4_exp <- lapply(exp1, function(x) {
+  suppressWarnings(
+    CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), 
+                       varName = 'test', Dates = as.POSIXct(dates_exp))
+  )
+})
+
+##############################################
+
+test_that("1. Initial checks CST_PeriodPET", {
+  # Check 's2dv_cube'
+  expect_error(
+    CST_PeriodPET(data = NULL), 
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    CST_PeriodPET(data = array(10)), 
+    "Parameter 'data' must be a list of 's2dv_cube' class."
+  )
+  # latitude
+  expect_error(
+    CST_PeriodPET(data = list(cube1)), 
+    paste0("Spatial coordinate names of parameter 'data' do not match any ",
+           "of the names accepted by the package.")
+  )
+  # Dates
+  expect_error(
+    CST_PeriodPET(data = list(cube2)),
+    paste0("Element 'Dates' is not found in 'attrs' list of 'data'. ", 
+           "See 's2dv_cube' object description in README file for more ", 
+           "information.")
+  )
+})
+
+##############################################
+
+test_that("1. Initial checks PeriodPET", {
+  # data
+  expect_error(
+    PeriodPET(data = NULL), 
+    "Parameter 'data' needs to be a named list with the needed variables."
+  )
+  expect_error(
+    PeriodPET(data = list(1)), 
+    "Parameter 'data' needs to be a named list with the variable names."
+  )
+  expect_error(
+    PeriodPET(data = list(tmax = array(10))), 
+    "Parameter 'data' needs to be a list of arrays with dimension names."
+  )
+  expect_error(
+    PeriodPET(data = list(tmax = array(10, c(time = 10)), 
+                          tmin = array(10, c(time = 11)))),
+    "Parameter 'data' variables need to have the same dimensions."
+  )
+  expect_error(
+    PeriodPET(data = list(tmax = array(10, c(time = 10)), 
+                          tmin = array(10, c(ftime = 10)))),
+    "Parameter 'data' variables need to have the same dimensions."
+  )
+  # lat
+  expect_error(
+    PeriodPET(data = exp1, lat = 'lat'),
+    "Parameter 'lat' must be numeric."
+  )
+  expect_error(
+    PeriodPET(data = list(tmax = array(10, c(time = 10)), 
+                          tmin = array(10, c(time = 10))), lat = 1:2),
+    "Parameter 'data' must have 'lat_dim' dimension."
+  )
+  # data (2)
+  expect_warning(
+    PeriodPET(data = exp1, pet_method = '1', dates = dates_exp, lat = lat), 
+    paste0("Parameter 'pet_method' needs to be 'hargreaves' or ", 
+           "'hargreaves_modified'. It is set to 'hargreaves_modified'.")
+  )
+  # time_dim
+  expect_error(
+    PeriodPET(data = exp1, time_dim = 1, dates = dates_exp, lat = lat),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodPET(data = exp2, lat = lat, dates = dates_exp2, 
+              lat_dim = 'lat', pet_method = 'thornthwaite'),
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  # leadtime_dim
+  expect_error(
+    PeriodPET(data = exp1, leadtime_dim = 1, dates = dates_exp, lat = lat),
+    "Parameter 'leadtime_dim' must be a character string."
+  )
+  expect_error(
+    PeriodPET(data = exp2, lat = lat, dates = dates_exp2, 
+               lat_dim = 'lat', time_dim = 'ftime', pet_method = 'thornthwaite'),
+    "Parameter 'leadtime_dim' is not found in 'data' dimension."
+  )
+  # lat_dim
+  expect_error(
+    PeriodPET(data = exp1, lat_dim = 1, dates = dates_exp, lat = lat)
+  )
+  expect_error(
+    PeriodPET(data = exp2, lat = lat, dates = dates_exp2),
+    "Parameter 'data' must have 'lat_dim' dimension."
+  )
+  # na.rm
+  expect_error(
+    PeriodPET(data = exp1, na.rm = 1.5, dates = dates_exp, lat = lat),
+    "Parameter 'na.rm' must be one logical value."
+  )
+  # ncores
+  expect_error(
+    PeriodPET(data = exp1, ncores = 1.5, dates = dates_exp, lat = lat),
+    "Parameter 'ncores' must be a positive integer."
+  ) 
+})
+
+##############################################
+
+test_that("2. Output checks", {
+  res01 <- CST_PeriodPET(data = cube4_exp)
+  res1 <- PeriodPET(data = exp1, lat = lat, dates = dates_exp)
+  res2 <- PeriodPET(data = exp2, lat = lat, dates = dates_exp2, 
+                    pet_method = c('thornthwaite'), 
+                    lat_dim = 'lat', time_dim = 'styear', 
+                    leadtime_dim = 'ftime')
+  # structure
+  expect_equal(
+    names(res01),
+    c('data', 'dims', 'coords', 'attrs')
+  )
+  # dims
+  expect_equal(
+    dim(res1),
+    c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10)
+  )
+  # values
+  expect_equal(
+    res1[1:4],
+    c(137.77342, 154.55548, 65.72859, 222.20438),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    res2[1:4],
+    c(77.76124, 118.94212, 66.57568, 185.67074),
+    tolerance = 0.0001
+  )
+})
+
+##############################################
\ No newline at end of file
diff --git a/tests/testthat/test-PeriodSPEI.R b/tests/testthat/test-PeriodSPEI.R
new file mode 100644
index 0000000000000000000000000000000000000000..2f1d7b0745bf4e507fb094965f6abf50ff4b7e10
--- /dev/null
+++ b/tests/testthat/test-PeriodSPEI.R
@@ -0,0 +1,481 @@
+##############################################
+# cube1
+cube1 <- NULL
+cube1$data <- array(rnorm(10), dim = c(lat = 2, time = 5))
+class(cube1) <- 's2dv_cube'
+
+# cube2
+cube2 <- NULL
+cube2$data <- array(rnorm(10), dim = c(lat = 2, time = 5))
+class(cube2) <- 's2dv_cube'
+cube2$coords <- list(lat = 1:2)
+
+# cube3
+cube3 <- NULL
+cube3$data <- array(rnorm(10), dim = c(lat = 2, time = 5))
+class(cube3) <- 's2dv_cube'
+cube3$coords <- list(lat = 1:2)
+cube3$attrs$Dates <- as.Date(c(paste0(2010:2014, "-08-16")))
+
+# dat1
+dims <- c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10)
+dimscor <- c(syear = 1, time = 3, latitude = 2, longitude = 1, ensemble = 15)
+
+set.seed(1)
+exp_tmax <- array(rnorm(360, 27.73, 5.26), dim = dims)
+set.seed(2)
+exp_tmin <- array(rnorm(360, 14.83, 3.86), dim = dims)
+set.seed(3)
+exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims)
+
+set.seed(1)
+expcor_tmax <- array(rnorm(60, 29.03, 5.67), dim = dimscor)
+set.seed(2)
+expcor_tmin <- array(rnorm(60, 15.70, 4.40), dim = dimscor)
+set.seed(3)
+expcor_pr <- array(rnorm(60, 15.62, 21.38), dim = dimscor)
+
+dates_exp <- as.POSIXct(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), 
+                          paste0(2010:2015, "-10-16")), "UTC")
+dim(dates_exp) <- c(syear = 6, time = 3)
+
+dates_expcor <- as.POSIXct(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), 
+                             paste0(2020, "-10-16")), "UTC")
+dim(dates_expcor) <- c(syear = 1, time = 3)
+
+lat <- c(40,40.1)
+
+exp1 <- list('tmax' = exp_tmax, 'tmin' = exp_tmin, 'pr' = exp_pr)
+exp_cor1 <- list('tmax' = expcor_tmax, 'tmin' = expcor_tmin, 'pr' = expcor_pr)
+params1 <- array(abs(rnorm(100)), dim = c(syear = 1, time = 3, latitude = 2, 
+                                          longitude = 1, coef = 3))
+
+# dat2
+dims2 <- c(styear = 6, ftime = 3, lat = 2, lon = 1, member = 10)
+dimscor2 <- c(styear = 1, ftime = 3, lat = 2, lon = 1, member = 15)
+
+set.seed(1)
+exp_tmean <- array(rnorm(100, 17.34, 9.18), dim = dims2)
+set.seed(2)
+exp_pr <- array(rnorm(360, 21.19, 25.64), dim = dims2)
+
+set.seed(1)
+expcor_tmean <- array(rnorm(100, 17.23, 9.19), dim = dimscor2)
+set.seed(2)
+expcor_pr <- array(rnorm(60, 15.62, 21.38), dim = dimscor2)
+
+dates_exp2 <- as.Date(c(paste0(2010:2015, "-08-16"), paste0(2010:2015, "-09-15"), 
+                        paste0(2010:2015, "-10-16")))
+dim(dates_exp2) <- c(sday = 1, sweek = 1, styear = 6, ftime = 3)
+
+dates_expcor2 <- as.Date(c(paste0(2020, "-08-16"), paste0(2020, "-09-15"), 
+                          paste0(2020, "-10-16")))
+dim(dates_expcor2) <- c(sday = 1, sweek = 1, styear = 1, ftime = 3)
+
+lat <- c(40,40.1)
+
+exp2 <- list('tmean' = exp_tmean, 'pr' = exp_pr)
+exp_cor2 <- list('tmean' = expcor_tmean, 'pr' = expcor_pr)
+
+# cube4
+cube4_exp <- lapply(exp1, function(x) {
+  suppressWarnings(
+    CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), 
+                       varName = 'test', Dates = dates_exp)
+  )
+})
+cube4_expcor <- lapply(exp_cor1, function(x) {
+  suppressWarnings(
+    CSTools::s2dv_cube(data = x, coords = list(latitude = c(40, 40.1)), 
+                       varName = 'test', Dates = dates_expcor)
+  )
+})
+
+##############################################
+
+test_that("1. Initial checks CST_PeriodSPEI", {
+  # Check 's2dv_cube'
+  expect_error(
+    CST_PeriodSPEI(exp = NULL), 
+    "Parameter 'exp' cannot be NULL."
+  )
+  expect_error(
+    CST_PeriodSPEI(exp = array(10)), 
+    "Parameter 'exp' must be a list of 's2dv_cube' class."
+  )
+  # latitude
+  expect_error(
+    CST_PeriodSPEI(exp = list(cube1)), 
+    paste0("Spatial coordinate names of parameter 'exp' do not match any ",
+           "of the names accepted by the package.")
+  )
+  # Dates
+  expect_error(
+    CST_PeriodSPEI(exp = list(cube2)),
+    paste0("Element 'Dates' is not found in 'attrs' list of 'exp'. ", 
+           "See 's2dv_cube' object description in README file for more ", 
+           "information.")
+  )
+  expect_error(
+    CST_PeriodSPEI(exp = list(cube3), exp_cor = list(cube2)),
+    paste0("Element 'Dates' is not found in 'attrs' list of 'exp_cor'. ", 
+           "See 's2dv_cube' object description in README file for more ", 
+           "information.")
+  )
+})
+
+##############################################
+
+test_that("1. Initial checks PeriodSPEI", {
+  # exp
+  expect_error(
+    PeriodSPEI(exp = NULL), 
+    "Parameter 'exp' needs to be a named list with the needed variables."
+  )
+  expect_error(
+    PeriodSPEI(exp = list(1)), 
+    "Parameter 'exp' needs to be a named list with the variable names."
+  )
+  expect_error(
+    PeriodSPEI(exp = list(tmax = array(10))), 
+    "Parameter 'exp' needs to be a list of arrays with dimension names."
+  )
+  expect_error(
+    PeriodSPEI(exp = list(tmax = array(10, c(time = 10)), 
+                          tmin = array(10, c(time = 11)))),
+    "Parameter 'exp' variables need to have the same dimensions."
+  )
+  expect_error(
+    PeriodSPEI(exp = list(tmax = array(10, c(time = 10)), 
+                          tmin = array(10, c(ftime = 10)))),
+    "Parameter 'exp' variables need to have the same dimensions."
+  )
+  # exp_cor
+  expect_error(
+    PeriodSPEI(exp = exp1, exp_cor = 1), 
+    paste0("Parameter 'exp_cor' needs to be a named list with the needed ", 
+           "variables if it is not NULL.")
+  )
+  expect_error(
+    PeriodSPEI(exp = exp1, exp_cor = list(1)), 
+    "Parameter 'exp_cor' needs to be a named list with the variable names."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp1, exp_cor = list('tmax' = array(10))), 
+    "Parameter 'exp_cor' needs to be a list of arrays with dimension names."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp1, exp_cor = list(tmax = array(10, c(time = 10)), 
+                                          tmin = array(10, c(time = 11)))),
+    "Parameter 'exp_cor' variables need to have the same dimensions."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp1, lat = 'lat', dates_exp = dates_exp),
+    "Parameter 'lat' must be numeric."
+  )
+  expect_error(
+    PeriodSPEI(exp = list(pr = array(10, c(time = 10, syear = 1, ensemble = 1))),  
+                          lat = 1:2, dates_exp = dates_exp),
+    "Parameter 'lat_dim' is not found in 'exp' dimension."
+  )
+  # exp (2)
+  expect_warning(
+    PeriodSPEI(exp = exp1, pet_method = '1', dates_exp = dates_exp, lat = lat), 
+    paste0("Parameter 'pet_method' needs to be 'hargreaves' or ", 
+           "'hargreaves_modified'. It is set to 'hargreaves_modified'.")
+  )
+  # time_dim
+  expect_error(
+    PeriodSPEI(exp = exp1, time_dim = 1, dates_exp = dates_exp, lat = lat),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, 
+               dates_exp = dates_exp2, dates_expcor = dates_expcor2,
+               lat_dim = 'lat', pet_method = 'thornthwaite'),
+    "Parameter 'time_dim' is not found in 'exp' dimension."
+  )
+  # leadtime_dim
+  expect_error(
+    PeriodSPEI(exp = exp1, leadtime_dim = 1, dates_exp = dates_exp, lat = lat),
+    "Parameter 'leadtime_dim' must be a character string."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, 
+               dates_exp = dates_exp2, dates_expcor = dates_expcor2,
+               lat_dim = 'lat', time_dim = 'ftime', pet_method = 'thornthwaite'),
+    "Parameter 'leadtime_dim' is not found in 'exp' dimension."
+  )
+  # memb_dim
+  expect_error(
+    PeriodSPEI(exp = exp1, memb_dim = 1, dates_exp = dates_exp, lat = lat),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, 
+               dates_exp = dates_exp2, dates_expcor = dates_expcor2,
+               lat_dim = 'lat', time_dim = 'ftime', leadtime_dim = 'styear',
+               pet_method = 'thornthwaite'),
+    "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  # lat_dim
+  expect_error(
+    PeriodSPEI(exp = exp1, lat_dim = 1, dates_exp = dates_exp, lat = lat)
+  )
+  expect_error(
+    PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, 
+               dates_exp = dates_exp2, dates_expcor = dates_expcor2),
+    "Parameter 'time_dim' is not found in 'exp' dimension."
+  )
+  # accum
+  expect_error(
+    PeriodSPEI(exp = exp1, accum = 10, dates_exp = dates_exp, lat = lat),
+    "Cannot compute accumulation of 10 months because loaded data has only 3 months."
+  ) 
+  # ref_period
+  expect_warning(
+    PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, 
+               dates_expcor = dates_expcor, lat = lat, ref_period = 1),
+    paste0("Parameter 'ref_period' must be of length two indicating the ",
+           "first and end years of the reference period. It will not ", 
+           "be used.")
+  )
+  expect_warning(
+    PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, 
+               dates_expcor = dates_expcor, lat = lat, ref_period = list('a', 1)),
+    paste0("Parameter 'ref_period' must be a numeric vector indicating the ", 
+           "'start' and 'end' years of the reference period. It will not ",
+           "be used.")
+  )
+  expect_warning(
+    PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, 
+               dates_expcor = dates_expcor, lat = lat, ref_period = list(2012, 2011)),
+    paste0("In parameter 'ref_period' 'start' cannot be after 'end'. It ",
+           "will not be used.")
+  )
+  expect_warning(
+    PeriodSPEI(exp = exp1, exp_cor = exp_cor1, dates_exp = dates_exp, 
+               dates_expcor = dates_expcor, lat = lat, ref_period = list(2008, 2021)),
+    paste0("Parameter 'ref_period' contain years outside the dates. ", 
+           "It will not be used.")
+  )
+  # standardization
+  expect_error(
+    PeriodSPEI(exp = exp1, standardization = 10, dates_exp = dates_exp, lat = lat),
+    "Parameter 'standardization' must be a logical value."
+  ) 
+  # param_error
+  expect_error(
+    PeriodSPEI(exp = exp1, param_error = TRUE, dates_exp = dates_exp, lat = lat),
+    "Parameter 'param_error' must be a numeric value."
+  ) 
+  # handle_infinity
+  expect_error(
+    PeriodSPEI(exp = exp1, handle_infinity = 1, dates_exp = dates_exp, lat = lat),
+    "Parameter 'handle_infinity' must be a logical value."
+  ) 
+  # cross_validation
+  expect_error(
+    PeriodSPEI(exp = exp1, cross_validation = 1, dates_exp = dates_exp, lat = lat),
+    "Parameter 'cross_validation' must be a logical value."
+  ) 
+  # method
+  expect_error(
+    PeriodSPEI(exp = exp1, method = 1, dates_exp = dates_exp, lat = lat),
+    paste0("Parameter 'method' must be a character string containing one ",
+           "of the following methods: 'parametric' or 'non-parametric'.")
+  ) 
+  # distribution
+  expect_error(
+    PeriodSPEI(exp = exp1, distribution = 1, dates_exp = dates_exp, lat = lat),
+    paste0("Parameter 'distribution' must be a character string containing one ", 
+           "of the following distributions: 'log-Logistic', 'Gamma' or ", 
+           "'PearsonIII'.")
+  ) 
+  # ncores
+  expect_error(
+    PeriodSPEI(exp = exp1, ncores = 1.5, dates_exp = dates_exp, lat = lat),
+    "Parameter 'ncores' must be a positive integer."
+  ) 
+})
+
+##############################################
+
+test_that("2. Output checks: CST_PeriodSPEI", {
+  res1 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = NULL)
+  res2 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = cube4_expcor)
+  res3 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = cube4_expcor, standardization = F)
+  res4 <- CST_PeriodSPEI(exp = cube4_exp, exp_cor = NULL, return_params = T)
+  expect_equal(
+    names(res1),
+    c("data", "dims", "coords", "attrs")
+  )
+  expect_equal(
+    res2$attrs$Variable$varName,
+    "SPEI"
+  )
+  expect_equal(
+    res3$attrs$Variable$varName,
+    "Precipitation minus accumulated PET"
+  )
+  expect_equal(
+    names(res3),
+    c("data", "dims", "coords", "attrs")
+  )
+  expect_equal(
+    names(res4),
+    c("spei", "params")
+  )
+})
+
+##############################################
+
+test_that("2. Output checks", {
+  res1 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+                     dates_exp = dates_exp, dates_expcor = dates_expcor,
+                     return_params = TRUE)
+  res2 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+                     dates_exp = dates_exp, dates_expcor = dates_expcor, 
+                     standardization = FALSE) # No info about accumulation
+  res3 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, 
+                     dates_exp = dates_exp, return_params = TRUE)
+  # output dims
+  expect_equal(
+    names(res1), 
+    c('spei', 'params')
+  )
+  expect_equal(
+    dim(res2), 
+    c(syear = 1, time = 3, latitude = 2, longitude = 1, ensemble = 15)
+  )
+  expect_equal(
+    names(res3), 
+    c('spei', 'params')
+  ) 
+  expect_equal(
+    dim(res1[[1]]),
+    dimscor
+  )
+  expect_equal(
+    dim(res1[[2]])[which(!names(dim(res1[[2]])) %in% c('coef', 'syear'))],
+    dims[which(!names(dims) %in% c('syear', 'ensemble'))]
+  )
+  expect_equal(
+    dim(res2),
+    dimscor
+  )
+  expect_equal(
+    dim(res3[[2]]),
+    c(syear = 1, time = 3, coef = 3, latitude = 2, longitude = 1)
+  )
+  # exp
+  # exp_cor
+  # pet
+  # accum
+  res11 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, accum = 2,
+                      dates_exp = dates_exp, na.rm = TRUE)
+  # expect_equal(
+  #   res11[1,3,1,1,][1:4], 
+  #   c(-0.4292409, -0.1375149, -0.5564081, -0.4273380),
+  #   tolerance = 0.0001
+  # )
+  # ref_period
+  res_ref <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, accum = 2,
+                        dates_exp = dates_exp, dates_expcor = dates_expcor,  
+                        na.rm = TRUE, ref_period = list(2011, 2013))
+  expect_equal(
+    !identical(res1[[1]], res_ref),
+    TRUE
+  )
+  # params
+  res5 <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, 
+                     params = params1, return_params = TRUE)
+  expect_error(
+    PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, 
+              params = array(abs(rnorm(100)), dim = dimscor)),
+    paste0("Parameter 'params' must be a multidimensional array with named ",
+           "dimensions: 'syear', 'time' and 'coef'.")
+  )
+  expect_equal(
+    dim(res5$params),
+    c(syear = 1, time = 3, coef = 3, latitude = 2, longitude = 1)
+  )
+  # standarization
+  res4 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, 
+                     dates_exp = dates_exp, standardization = FALSE)
+  expect_equal(
+    names(res4), 
+    NULL
+  )
+  expect_equal(
+    dim(res4),
+    c(syear = 6, time = 3, latitude = 2, longitude = 1, ensemble = 10)
+  )  
+  # cross_validation 
+  expect_warning(
+    PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+               dates_exp = dates_exp, dates_expcor = dates_expcor, 
+               cross_validation = TRUE),
+    paste0("Detected 'cross_validation' = TRUE. It will be set as FALSE ",
+           "since 'exp_cor' is provided.")
+  )
+  res_crossval_T <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, 
+                               cross_validation = TRUE, return_params = TRUE)
+  res_crossval_F <- PeriodSPEI(exp = exp1, lat = lat, dates_exp = dates_exp, 
+                               cross_validation = FALSE, return_params = TRUE)
+  # cross_validation = TRUE
+  expect_equal(
+    dim(res_crossval_T$spei),
+    dims
+  )
+  expect_equal(
+    dim(res_crossval_T$params),
+    c(syear = 6, time = 3, coef = 3, latitude = 2, longitude = 1)
+  )
+  # cross_validation = FALSE
+  expect_equal(
+    dim(res_crossval_F$params)[-which(names(dim(res_crossval_F$params)) == 'coef')],
+    dimscor[-which(names(dimscor) == 'ensemble')]
+  )
+
+  # pet_method - ok
+  res5 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+                     dates_exp = dates_exp, dates_expcor = dates_expcor,
+                     pet_method = c('hargreaves', 'hargreaves_modified'))
+  res6 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+                     dates_exp = dates_exp, dates_expcor = dates_expcor,
+                     pet_method = c('hargreaves_modified', 'hargreaves'))
+  expect_equal(
+    identical(res5, res6),
+    FALSE
+  )
+  expect_equal(
+    dim(res5),
+    dim(res6)
+  )
+
+  # time_dim, leadtime_dim, memb_dim, lat_dim
+  res7 <- PeriodSPEI(exp = exp2, exp_cor = exp_cor2, lat = lat, 
+                     dates_exp = dates_exp2, dates_expcor = dates_expcor2,
+                     pet_method = c('thornthwaite', 'thornthwaite'), 
+                     lat_dim = 'lat', time_dim = 'styear', 
+                     leadtime_dim = 'ftime', memb_dim = 'member')
+  # method - ok
+  res8 <- PeriodSPEI(exp = exp1, exp_cor = exp_cor1, lat = lat, 
+                     dates_exp = dates_exp, dates_expcor = dates_expcor,
+                     method = 'non-parametric')
+  # distribution - Only works for 'log-Logistic'
+  res9 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, 
+                     dates_exp = dates_exp, distribution = 'PearsonIII') # NA
+  res10 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, 
+                      dates_exp = dates_exp, distribution = 'Gamma') # NA
+  # param_error -
+  # handle_infinity - OK
+  res9 <- PeriodSPEI(exp = exp1, exp_cor = NULL, lat = lat, 
+                     dates_exp = dates_exp, handle_infinity = FALSE)
+  # na.rm -
+  
+  # ncores
+})
+
+##############################################
\ No newline at end of file