diff --git a/NAMESPACE b/NAMESPACE
index e7d7c003b822cbb73b70d26b94fb43d55faab57e..bd5d0f1649b1d8529ffb6d16118bc8c87ec186ab 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -11,6 +11,7 @@ export(CST_Calibration)
 export(CST_CategoricalEnsCombination)
 export(CST_EnsClustering)
 export(CST_Load)
+export(CST_MergeDims)
 export(CST_MultiEOF)
 export(CST_MultiMetric)
 export(CST_MultivarRMSE)
@@ -21,6 +22,7 @@ export(CST_RainFARM)
 export(CST_SaveExp)
 export(CST_SplitDim)
 export(EnsClustering)
+export(MergeDims)
 export(MultiEOF)
 export(PlotCombinedMap)
 export(PlotForecastPDF)
diff --git a/R/CST_BEI_Weighting.R b/R/CST_BEI_Weighting.R
index 51c8ba32fdaede17917f610ecd685aa6991ebcda..de7470110c4f9d090bc1ad84b6c786448a5c63c4 100644
--- a/R/CST_BEI_Weighting.R
+++ b/R/CST_BEI_Weighting.R
@@ -19,6 +19,10 @@
 #' (time, member), when 'time' is the temporal dimension as default. 
 #' When 'aweights' parameter has any other dimensions (as e.g. 'lat') and 
 #' 'var_exp' parameter has also the same dimension, they must be equals.  
+#' @param terciles A numeric array with at least one dimension 'tercil' equal to 
+#' 2, the first element is the lower tercil for a hindcast period, and the second  
+#' element is the upper tercile. By default is NULL, the terciles are computed  
+#' from var_exp data.
 #' @param type A character string indicating the type of output. 
 #' If 'type' =  'probs', the function returns, in the element data from 
 #' 'var_exp' parameter, an array with at least two 
@@ -36,7 +40,7 @@
 #' with weighted members.
 #' @param time_dim_name A character string indicating the name of the 
 #' temporal dimension, by default 'time'.
-#'  
+#' 
 #' @return CST_BEI_Weighting() returns a CSTools object (i.e., of the
 #' class 's2dv_cube').
 #' This object has at least an element named \code{$data}
@@ -56,8 +60,10 @@
 #' # time     lat     lon   dataset
 #' #    2       3       2         2
 #' @export
-CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean', 
-                         time_dim_name = 'time') {
+
+CST_BEI_Weighting <- function(var_exp, aweights, terciles = NULL, 
+                              type = 'ensembleMean', time_dim_name = 'time') {
+  
   if (!is.character(time_dim_name)) {
     stop("Parameter 'time_dim_name' must be a character string indicating",
          " the name of the temporal dimension.")
@@ -80,6 +86,27 @@ CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean',
     stop("Parameter 'var_exp' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
   }
+  if (!is.null(terciles)){
+    if(!is.array(terciles)){
+      stop("Parameter 'terciles' must be an array.")
+    } 
+    if (is.null(names(dim(terciles)))) {
+      stop("Parameters 'terciles' should have dimmension names.")
+    }
+    if(!('tercil' %in% names(dim(terciles)))) {
+      stop("Parameter 'terciles' must have dimension 'tercil'.")
+    }
+    if (dim(terciles)['tercil'] != 2) {
+      stop("Length of dimension 'tercil' ",
+           "of parameter 'terciles' must be equal to 2.")
+    }
+    if(time_dim_name %in% names(dim(terciles))) {
+      stop("Parameter 'terciles' must not have temporal dimension.")
+    }
+    if('member' %in% names(dim(terciles))) {
+      stop("Parameter 'terciles' must not have dimension 'member'.")
+    }
+  }
   if (!is.array(aweights)) {
     stop("Parameter 'aweights' must be an array.")
   }
@@ -115,7 +142,10 @@ CST_BEI_Weighting <- function(var_exp, aweights, type = 'ensembleMean',
     em <- BEI_EMWeighting(var_exp$data, aweights, time_dim_name)
     var_exp$data <- em
   } else if (type == 'probs'){
-    probs <- BEI_ProbsWeighting(var_exp$data, aweights, time_dim_name)
+    if (is.null(terciles)){
+      terciles <- BEI_TercilesWeighting(var_exp$data, aweights, time_dim_name)
+    }
+    probs <- BEI_ProbsWeighting(var_exp$data, aweights, terciles, time_dim_name)
     var_exp$data <- probs
   } else {
       stop("Parameter 'type' must be a character string ('probs' or ",
@@ -265,6 +295,9 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
 #' variable, as 'time' the spatial dimension by default.
 #' @param aweights Normalized weights array with at least dimensions 
 #' (time, member), when 'time' is the temporal dimension as default.
+#' @param terciles A numeric array with at least one dimension 'tercil' equal to 
+#' 2, the first element is the lower tercil for a hindcast period, and the second  
+#' element is the upper tercile. 
 #' @param time_dim_name A character string indicating the name of the 
 #' temporal dimension, by default 'time'.
 #' 
@@ -284,7 +317,9 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
 #' dim(var_exp) <- c(time = 2, member = 4)
 #' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
 #' dim(aweights) <- c(time = 2, member = 4)
-#' res <- BEI_ProbsWeighting(var_exp, aweights)
+#' terciles <- c(2.5,5)
+#' dim(terciles) <- c(tercil = 2)
+#' res <- BEI_ProbsWeighting(var_exp, aweights, terciles)
 #' dim(res)
 #' # time tercil 
 #' #    2      3 
@@ -293,12 +328,15 @@ BEI_EMWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
 #' dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3)
 #' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
 #' dim(aweights) <- c(time = 2, member = 4)
-#' res <- BEI_ProbsWeighting(var_exp, aweights)
+#' terciles <- rep(c(48,50), 2*3)
+#' dim(terciles) <- c(tercil = 2, lat = 2, lon = 3)
+#' res <- BEI_ProbsWeighting(var_exp, aweights, terciles)
 #' dim(res)
 #' # time  tercil  lat     lon
 #' #    2       3    2       3
 #' @noRd
-BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
+BEI_ProbsWeighting <- function(var_exp, aweights, terciles, 
+                               time_dim_name = 'time') {
   
   if (!is.character(time_dim_name)) {
     stop("Parameter 'time_dim_name' must be a character string indicating",
@@ -309,6 +347,28 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
             "only the first element will be used.")
     time_dim_name <- time_dim_name[1]
   }
+  if (is.null(terciles)){
+    stop("Parameter 'terciles' is null")
+  }
+  if(!is.array(terciles)){
+    stop("Parameter 'terciles' must be an array.")
+  } 
+  if (is.null(names(dim(terciles)))) {
+    stop("Parameters 'terciles' should have dimmension names.")
+  }
+  if(!('tercil' %in% names(dim(terciles)))) {
+    stop("Parameter 'terciles' must have dimension 'tercil'.")
+  }
+  if (dim(terciles)['tercil'] != 2) {
+    stop("Length of dimension 'tercil' ",
+         "of parameter 'terciles' must be equal to 2.")
+  }
+  if(time_dim_name %in% names(dim(terciles))) {
+    stop("Parameter 'terciles' must not have temporal dimension.")
+  }
+  if('member' %in% names(dim(terciles))) {
+    stop("Parameter 'terciles' must not have dimension 'member'.")
+  }
   if (!is.array(var_exp)) {
     stop("Parameter 'var_exp' must be an array.")
   }
@@ -326,10 +386,10 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
     stop("Parameter 'aweights' must have temporal dimension.")
   }
   if(!('member' %in% names(dim(var_exp)))) {
-    stop("Parameter 'var_exp' must have temporal dimension.")
+    stop("Parameter 'var_exp' must have dimension 'member'.")
   }
   if(!('member' %in% names(dim(aweights)))) {
-    stop("Parameter 'aweights' must have temporal dimension.")
+    stop("Parameter 'aweights' must have dimension 'member'.")
   }
   if (dim(var_exp)[time_dim_name] != dim(aweights)[time_dim_name]) {
     stop("Length of temporal dimensions ",
@@ -340,8 +400,11 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
          "of parameter 'var_exp' and 'aweights' must be equals.")
   }
   
-  res <- Apply(list(var_exp, aweights),
-               target_dims = list(c(time_dim_name,'member'), c(time_dim_name,'member')),
+  
+  res <- Apply(list(var_exp, aweights, terciles),
+               target_dims = list(c(time_dim_name,'member'), 
+                                  c(time_dim_name,'member'),
+                                  c('tercil')),
                fun = .BEI_ProbsWeighting, time_dim_name)$output1
   return(res)
 }
@@ -352,8 +415,12 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
 #' by default 'time', and dimension 'member'.
 #' @param aweights Normalized weights array with a temporal dimension, 
 #' by default 'time', and dimension 'member'
+#' @param terciles A numeric array with one dimension 'tercil' equal to 2, 
+#' the first element is the lower tercil for a hindcast period, and the second  
+#' element is the upper tercile.
 #' @param time_dim_name A character string indicating the name of the 
 #' temporal dimension, by default 'time'.
+#' 
 #' @return .BEI_ProbsWeighting returns an array of with a temporal dimension,
 #' as default 'time', and 'tercil' dimension, containing the probabilities 
 #' for each tercile computing with weighted members.
@@ -363,48 +430,207 @@ BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
 #' # Example
 #' var_exp <- 1 : 8
 #' dim(var_exp) <- c(stime = 2, member = 4)
-#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.3, 0.3)
+#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
 #' dim(aweights) <- c(stime = 2, member = 4)
-#' res <- .BEI_ProbsWeighting(var_exp, aweights, time_dim_name = 'stime')
+#' terciles <- quantile(1:8, probs = c(1/3, 2/3))
+#' dim(terciles) <- c(tercil = 2)
+#' res <- .BEI_ProbsWeighting(var_exp, aweights, terciles, time_dim_name = 'stime')
 #' dim(res)
 #' # stime tercil
 #' #     2      3
 #' @noRd
-.BEI_ProbsWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
-  # computing terciles
-  terciles_exp <- WeightTerciles(var_exp, aweights, time_dim_name)
-  lowerTercile <- terciles_exp$lowerTercile
-  upperTercile <- terciles_exp$upperTercile
-  
-  # Probabilities
-  aTerciles <- Apply(list(var_exp), target_dims = list('member'),
-                     fun = Data2Tercil, lowerTercile, upperTercile)$output1
+.BEI_ProbsWeighting <- function(var_exp, aweights, terciles, 
+                                time_dim_name = 'time') {
+  if(any(is.na(var_exp)) || any(is.na(aweights))){
+    probTercile <- array(NA, dim = c(dim(var_exp)[time_dim_name], tercil = 3))
+  } else {
+    if(any(is.na(terciles))) stop("Terciles are NAs")
+    terciles_exp <- list(lowerTercile = terciles[1], 
+                         upperTercile = terciles[2])
+    
+    lowerTercile <- terciles_exp$lowerTercile
+    upperTercile <- terciles_exp$upperTercile
+    
+    # Probabilities
+    aTerciles <- Apply(list(var_exp), target_dims = list('member'),
+                       fun = Data2Tercil, lowerTercile, upperTercile)$output1
+    
+    pos <- match(names(dim(aTerciles)), c(time_dim_name,'member'))
+    aTerciles <- aperm(aTerciles,pos)
+    names(dim(aTerciles)) <- c(time_dim_name,'member')
+    
+    probTercile <- array(NA, dim = c(dim(var_exp)[time_dim_name], tercil = 3))
+    for (idTercil in 1:3){
+      probTercile[,idTercil] <- Apply(list(aTerciles, aweights),
+                                      target_dims = list('member','member'),
+                                      fun = WeightTercil2Prob, idTercil)$output1
+    }
+  }
+  return(probTercile)
+}
+
+#' Computing the weighted terciles for SFSs.
+#' @author Eroteida Sanchez-Garcia - AEMET, \email{esanchezg@aemet.es}
+#' 
+#' @description This function implements the computation to obtain the terciles
+#' for a weighted variable for SFSs using a normalized weights array,
+#' 
+#' @references Regionally improved seasonal forecast of precipitation through Best
+#' estimation of winter NAO, Sanchez-Garcia, E. et al.,
+#' Adv. Sci. Res., 16, 165174, 2019, https://doi.org/10.5194/asr-16-165-2019
+#' 
+#' @param var_exp Variable (e.g. precipitation, temperature, NAO index)
+#' array from a SFS with at least dimensions (time, member) for a spatially 
+#' aggregated variable or dimensions (time, member, lat, lon) for a spatial 
+#' variable, as 'time' the spatial dimension by default.
+#' @param aweights Normalized weights array with at least dimensions 
+#' (time, member), when 'time' is the temporal dimension as default.
+#' @param time_dim_name A character string indicating the name of the 
+#' temporal dimension, by default 'time'.
+#' 
+#' @return BEI_TercilesWeighting() returns an array with at least one 
+#' dimension depending if the variable is a spatially aggregated variable 
+#' (as e.g. NAO index)(tercil) or it is spatial variable (as e.g. 
+#' precipitation or temperature)(tercil, lat, lon), containing the 
+#' terciles computing with weighted members.
+#' The first tercil is the lower tercile, the second is the upper tercile.
+#'  
+#' @import multiApply
+#' 
+#' @examples
+#' # Example 1 
+#' var_exp <- 1 : (2 * 4)
+#' dim(var_exp) <- c(time = 2, member = 4)
+#' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
+#' dim(aweights) <- c(time = 2, member = 4)
+#' res <- BEI_TercilesWeighting(var_exp, aweights)
+#' dim(res)
+#' # tercil 
+#' #      2
+#' # Example 2 
+#' var_exp <- rnorm(48, 50, 9)
+#' dim(var_exp) <- c(time = 2, member = 4, lat = 2, lon = 3)
+#' aweights<- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
+#' dim(aweights) <- c(time = 2, member = 4)
+#' res <- BEI_TercilesWeighting(var_exp, aweights)
+#' dim(res)
+#' # tercil  lat     lon
+#' #      2    2       3
+#' @noRd
+BEI_TercilesWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
   
-  pos <- match(names(dim(aTerciles)), c(time_dim_name,'member'))
-  aTerciles <- aperm(aTerciles,pos)
-  names(dim(aTerciles)) <- c(time_dim_name,'member')
+  if (!is.character(time_dim_name)) {
+    stop("Parameter 'time_dim_name' must be a character string indicating",
+         " the name of the temporal dimension.")
+  }
+  if (length(time_dim_name) > 1) {
+    warning("Parameter 'time_dim_name' has length greater than 1 and ",
+            "only the first element will be used.")
+    time_dim_name <- time_dim_name[1]
+  }
+  if (!is.array(var_exp)) {
+    stop("Parameter 'var_exp' must be an array.")
+  }
+  if (!is.array(aweights)) {
+    stop("Parameter 'aweights' must be an array.")
+  }
+  if (is.null(names(dim(var_exp))) || is.null(names(dim(aweights)))) {
+    stop("Parameters 'var_exp' and 'aweights'",
+         " should have dimmension names.")
+  }
+  if(!(time_dim_name %in% names(dim(var_exp)))) {
+    stop("Parameter 'var_exp' must have temporal dimension.")
+  }
+  if(!(time_dim_name %in% names(dim(aweights)))) {
+    stop("Parameter 'aweights' must have temporal dimension.")
+  }
+  if(!('member' %in% names(dim(var_exp)))) {
+    stop("Parameter 'var_exp' must have temporal dimension.")
+  }
+  if(!('member' %in% names(dim(aweights)))) {
+    stop("Parameter 'aweights' must have temporal dimension.")
+  }
+  if (dim(var_exp)[time_dim_name] != dim(aweights)[time_dim_name]) {
+    stop("Length of temporal dimensions ",
+         "of parameter 'var_exp' and 'aweights' must be equals.")
+  }
+  if (dim(var_exp)['member'] != dim(aweights)['member']) {
+    stop("Length of dimension 'member' ",
+         "of parameter 'var_exp' and 'aweights' must be equals.")
+  }
   
-  probTercile <- array(NA, dim = c(dim(var_exp)[time_dim_name], tercil = 3))
-  for (idTercil in 1:3){
-    probTercile[,idTercil] <- Apply(list(aTerciles, aweights),
-                                    target_dims = list('member','member'),
-                                    fun = WeightTercil2Prob, idTercil)$output1
+  res <- Apply(list(var_exp, aweights),
+               target_dims = list(c(time_dim_name,'member'), c(time_dim_name,'member')),
+               fun = .BEI_TercilesWeighting, time_dim_name)$output1
+  return(res)
+}
+
+#' Atomic BEI_TercilesWeighting
+#' @param var_exp Variable (e.g. precipitation, temperature, NAO index)
+#' array from a SFS with a temporal dimension, 
+#' by default 'time', and dimension 'member'.
+#' @param aweights Normalized weights array with a temporal dimension, 
+#' by default 'time', and dimension 'member'
+#' @param time_dim_name A character string indicating the name of the 
+#' temporal dimension, by default 'time'.
+#' @return .BEI_TercilesWeighting returns a numeric array with dimension tercil
+#' equal to 2, the first is the lower tercil and the second the upper tercile, 
+#' computing with weighted members considering all members and all period.
+#' If any member value for any period is NA , the terciles are not computed, and
+#' the function return NA value as tercile upper and lower.
+#' @examples
+#' # Example
+#' var_exp <- 1 : 8
+#' dim(var_exp) <- c(stime = 2, member = 4)
+#' aweights <- c(0.2, 0.1, 0.3, 0.4, 0.1, 0.2, 0.4, 0.3)
+#' dim(aweights) <- c(stime = 2, member = 4)
+#' res <- .BEI_TercilesWeighting(var_exp, aweights, time_dim_name = 'stime')
+#' dim(res)
+#' # tercil
+#' #      2
+#' @noRd
+.BEI_TercilesWeighting <- function(var_exp, aweights, time_dim_name = 'time') {
+  if(any(is.na(var_exp)) || any(is.na(aweights))){
+    terciles_exp <- array(c(NA, NA), dim = c(tercil = 2))
+  } else {
+    l_terciles_exp <- WeightTerciles(var_exp, aweights, time_dim_name)
+    terciles_exp <- array(c(l_terciles_exp$lowerTercile, 
+                            l_terciles_exp$upperTercile), dim = c(tercil = 2))
   }
-  return(probTercile)
+  return(terciles_exp)
 }
 
+# Auxiliar function to compute in which tercile is a data value
+Data2Tercil_old <- function(x,lt,ut) {
+  if(is.na(lt) || is.na(ut)){
+    y <- rep(NA, length(x))
+  } else {
+    y <- rep(2,length(x))
+    y[x <= lt] <- 1
+    y[x >= ut] <- 3
+    if (lt == ut) {
+      warning("The upper and lower terciles are equals")
+    }
+  }
+  dim(y) <- c(member = length(x))
+  return (y)
+}
 # Auxiliar function to compute in which tercile is a data value
 Data2Tercil <- function(x,lt,ut) {
-  y <- rep(2,length(x))
-  y[x <= lt] <- 1
-  y[x >= ut] <- 3
-  if (lt == ut) {
-    warning("The upper and lower terciles are equals")
+  if(is.na(lt) || is.na(ut)){
+    y <- rep(NA, length(x))
+  } else {
+    y <- rep(2,length(x))
+    y[x <= lt] <- 1
+    y[x >= ut] <- 3
+    if (lt == ut) {
+      y <- rep(NA, length(x))
+    }
   }
   dim(y) <- c(member = length(x))
+  y[which(is.na(x))] <- NA
   return (y)
 }
-
 # Auxiliar function to convers weighted terciles to probabilities
 WeightTercil2Prob <- function(aTerciles, aWeights, idTercil) {
   return(sum(aWeights[which(aTerciles == idTercil)]))
@@ -419,10 +645,12 @@ WeightTerciles <- function(data, aweights, time_dim_name = 'time') {
   names(dim(aweights)) <- namesdimdata
   vectorData <- as.vector(data)
   vectorWeights <- as.vector(aweights/dim(aweights)[time_dim_name]) # normalized
-  lSortData <- sort(vectorData,index.return=TRUE)
-  indSort <- lSortData$ix # index asociated to weight
+  #lSortData <- sort(vectorData,index.return=TRUE)
+  indSort <- order(vectorData) # index asociated to weight
+  # indSort <- lSortData$ix # index asociated to weight
   # corresponding for this data
-  dataSort <- lSortData$x
+  dataSort <- vectorData[indSort]
+  # dataSort <- lSortData$x
   # Adding normalized weights. When 1/3 is reached, the data value
   # is lower tercile and when 2/3 is reached, it is the upper tercile.
   sumWeights <- 0
diff --git a/R/CST_MergeDims.R b/R/CST_MergeDims.R
index 560f6c97b56e651e7cee345d0fdfe7218827c426..a9923c581b95e6c48c37d2efa2e32d7588a947f0 100644
--- a/R/CST_MergeDims.R
+++ b/R/CST_MergeDims.R
@@ -6,7 +6,7 @@
 #'
 #'@param data a 's2dv_cube' object
 #'@param merge_dims a character vector indicating the names of the dimensions to merge
-#'@param remane_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.
+#'@param rename_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.
 #'@param na.rm a logical indicating if the NA values should be removed or not.
 #'
 #'@import abind
@@ -45,7 +45,7 @@ CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), rename_dim =
 #'
 #'@param data an n-dimensional array with named dimensions
 #'@param merge_dims a character vector indicating the names of the dimensions to merge
-#'@param remane_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.
+#'@param rename_dim a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.
 #'@param na.rm a logical indicating if the NA values should be removed or not.
 #'
 #'@import abind
diff --git a/man/CST_BEI_Weighting.Rd b/man/CST_BEI_Weighting.Rd
index 6b9a448ae742037887d2e058e679a08edccf71d0..0e60a356d6f3e541775aa55bd398710a855a7900 100644
--- a/man/CST_BEI_Weighting.Rd
+++ b/man/CST_BEI_Weighting.Rd
@@ -4,7 +4,7 @@
 \alias{CST_BEI_Weighting}
 \title{Weighting SFSs of a CSTools object.}
 \usage{
-CST_BEI_Weighting(var_exp, aweights, type = "ensembleMean",
+CST_BEI_Weighting(var_exp, aweights, terciles = NULL, type = "ensembleMean",
   time_dim_name = "time")
 }
 \arguments{
@@ -18,6 +18,11 @@ at least a temporal dimension and a dimension named 'member'.}
 When 'aweights' parameter has any other dimensions (as e.g. 'lat') and 
 'var_exp' parameter has also the same dimension, they must be equals.}
 
+\item{terciles}{A numeric array with at least one dimension 'tercil' equal to 
+2, the first element is the lower tercil for a hindcast period, and the second  
+element is the upper tercile. By default is NULL, the terciles are computed  
+from var_exp data.}
+
 \item{type}{A character string indicating the type of output. 
 If 'type' =  'probs', the function returns, in the element data from 
 'var_exp' parameter, an array with at least two 
diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd
index 485199eae252039266185188df10ee21cf1af52b..e8a82af0b9e67719a600e1e6b04e9fec2c0d96b6 100644
--- a/man/CST_BiasCorrection.Rd
+++ b/man/CST_BiasCorrection.Rd
@@ -41,4 +41,4 @@ Verónica Torralba, \email{veronica.torralba@bsc.es}
 \references{
 Torralba, V., F.J. Doblas-Reyes, D. MacLeod, I. Christel and M. Davis (2017). Seasonal climate prediction: a new source of information for the management of wind energy resources. Journal of Applied Meteorology and Climatology, 56, 1231-1247, doi:10.1175/JAMC-D-16-0204.1. (CLIM4ENERGY, EUPORIAS, NEWA, RESILIENCE, SPECS)
 }
-\encoding{UTF-8}
+
diff --git a/man/CST_Calibration.Rd b/man/CST_Calibration.Rd
index 36171dbde22451fd681d7dbb899f5f54379a410e..210c080fe5bb6489468c92f90c490579c8c10a3d 100644
--- a/man/CST_Calibration.Rd
+++ b/man/CST_Calibration.Rd
@@ -57,4 +57,4 @@ attr(obs, 'class') <- 's2dv_cube'
 a <- CST_Calibration(exp = exp, obs = obs, cal.method = "mse_min", eval.method = "in-sample")
 str(a)
 }
-\encoding{UTF-8}
+
diff --git a/man/CST_MergeDims.Rd b/man/CST_MergeDims.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..449e011ee7ea5088996f0255206974e01ba91319
--- /dev/null
+++ b/man/CST_MergeDims.Rd
@@ -0,0 +1,41 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_MergeDims.R
+\name{CST_MergeDims}
+\alias{CST_MergeDims}
+\title{Function to  Merge Dimensions}
+\usage{
+CST_MergeDims(data, merge_dims = c("ftime", "monthly"), rename_dim = NULL,
+  na.rm = FALSE)
+}
+\arguments{
+\item{data}{a 's2dv_cube' object}
+
+\item{merge_dims}{a character vector indicating the names of the dimensions to merge}
+
+\item{rename_dim}{a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.}
+
+\item{na.rm}{a logical indicating if the NA values should be removed or not.}
+}
+\description{
+This function merges two dimensions of the array \code{data} in a 's2dv_cube' object into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them.
+}
+\examples{
+
+data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7)
+dim(data) <- c(time = 7, lat = 2, lon = 3, monthly = 4, member = 6,
+              dataset = 5, var = 1)
+data[2,,,,,,] <- NA
+data[c(3,27)] <- NA
+data <-list(data = data)
+class(data) <- 's2dv_cube'
+new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'))
+dim(new_data$data)
+new_data <- CST_MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid')
+dim(new_data$data)
+new_data <- CST_MergeDims(data, merge_dims = c('time', 'monthly'), na.rm = TRUE)
+dim(new_data$data)
+}
+\author{
+Nuria Perez-Zanon, \email{nuria.perez@bsc.es}
+}
+
diff --git a/man/MergeDims.Rd b/man/MergeDims.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..585049e883c5be79772d01462444f820aea79aba
--- /dev/null
+++ b/man/MergeDims.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_MergeDims.R
+\name{MergeDims}
+\alias{MergeDims}
+\title{Function to Split Dimension}
+\usage{
+MergeDims(data, merge_dims = c("time", "monthly"), rename_dim = NULL,
+  na.rm = FALSE)
+}
+\arguments{
+\item{data}{an n-dimensional array with named dimensions}
+
+\item{merge_dims}{a character vector indicating the names of the dimensions to merge}
+
+\item{rename_dim}{a character string indicating the name of the output dimension. If left at NULL, the first dimension name provided in parameter \code{merge_dims} will be used.}
+
+\item{na.rm}{a logical indicating if the NA values should be removed or not.}
+}
+\description{
+This function merges two dimensions of an array into one. The user can select the dimensions to merge and provide the final name of the dimension. The user can select to remove NA values or keep them.
+}
+\examples{
+
+data <- 1 : 20
+dim(data) <- c(time = 10, lat = 2)
+new_data <- MergeDims(data, merge_dims = c('time', 'lat'))
+}
+\author{
+Nuria Perez-Zanon, \email{nuria.perez@bsc.es}
+}
+
diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd
index bed0bd3137789bc45216ab058c8d07d9b019f167..3d1d65c3335738ee1878723ec8c28f1a0e001fc9 100644
--- a/man/PlotForecastPDF.Rd
+++ b/man/PlotForecastPDF.Rd
@@ -49,4 +49,4 @@ PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2),
 \author{
 Llorenç Lledó \email{llledo@bsc.es}
 }
-\encoding{UTF-8}
+
diff --git a/vignettes/Figures/BestEstimateIndex_fig2.png b/vignettes/Figures/BestEstimateIndex_fig2.png
index e10333a577d1a2cba98d4efa702af3f6e8c14bab..24dbceb26d63f9c71765efd849130f0302780e04 100644
Binary files a/vignettes/Figures/BestEstimateIndex_fig2.png and b/vignettes/Figures/BestEstimateIndex_fig2.png differ