diff --git a/DESCRIPTION b/DESCRIPTION
index d76276602a8a3d7997e752b9937a613bfbd46167..7fd6c5d40393d199acca72430f517865536fbfa5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -47,7 +47,8 @@ Description: Exploits dynamical seasonal forecasts in order to provide
     Yiou et al. (2013) <doi:10.1007/s00382-012-1626-3>.
 Depends:
     R (>= 3.4.0),
-    maps
+    maps,
+    easyVerification
 Imports:
     s2dverification,
     s2dv,
diff --git a/NAMESPACE b/NAMESPACE
index 983a892656c807ef4ad315c44301d6aca22c2894..29dc0811bb0d77faefcd9658edd5343a3bb2b944 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -28,6 +28,7 @@ export(Calibration)
 export(EnsClustering)
 export(MergeDims)
 export(MultiEOF)
+export(MultiMetric)
 export(PlotCombinedMap)
 export(PlotForecastPDF)
 export(PlotMostLikelyQuantileMap)
@@ -49,7 +50,6 @@ import(multiApply)
 import(ncdf4)
 import(qmap)
 import(rainfarmr)
-import(s2dverification)
 import(stats)
 importFrom(ClimProjDiags,SelBox)
 importFrom(RColorBrewer,brewer.pal)
@@ -57,6 +57,8 @@ importFrom(abind,abind)
 importFrom(data.table,CJ)
 importFrom(data.table,data.table)
 importFrom(data.table,setkey)
+importFrom(easyVerification,climFairRpss)
+importFrom(easyVerification,veriApply)
 importFrom(grDevices,adjustcolor)
 importFrom(grDevices,bmp)
 importFrom(grDevices,colorRampPalette)
@@ -86,7 +88,22 @@ importFrom(maps,map)
 importFrom(plyr,.)
 importFrom(plyr,dlply)
 importFrom(reshape2,melt)
+importFrom(s2dv,ColorBar)
+importFrom(s2dv,Corr)
+importFrom(s2dv,InsertDim)
+importFrom(s2dv,MeanDims)
+importFrom(s2dv,PlotEquiMap)
+importFrom(s2dv,RMS)
+importFrom(s2dv,RMSSS)
 importFrom(s2dv,Reorder)
+importFrom(s2dverification,ACC)
+importFrom(s2dverification,Ano_CrossValid)
+importFrom(s2dverification,Clim)
+importFrom(s2dverification,EOF)
+importFrom(s2dverification,Eno)
+importFrom(s2dverification,InsertDim)
+importFrom(s2dverification,Load)
+importFrom(s2dverification,Mean1Dim)
 importFrom(s2dverification,Subset)
 importFrom(utils,glob2rx)
 importFrom(utils,head)
diff --git a/NEWS.md b/NEWS.md
index e57b835859743a5c64e74f3ef5367139517de885..8889b807ad9beb26c625dea9cf0b33f64d0c4fd1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -6,11 +6,13 @@
     + CST_RFWeights accepts s2dv_cube objects as input and new 'ncores' paramenter allows parallel computation
     + RFWeights is exposed to users
     + CST_RainFARM accepts multi-dimensional slopes and weights and handless missing values in sample dimensions.
+    + CST_MultiMetric includes 'rpss' metric and it is addapted to s2dv.
 
 - Fixes:
     + PlotForecastPDF correctly displays terciles labels 
     + CST_SaveExp correctly save time units
     + CST_SplitDims returns ordered output following ascending order provided in indices when it is numeric
+
 ### CSTools 3.1.0
 **Submission date to CRAN: 02-07-2020**
 
diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R
index 6e33c3411d1000abb0c15e4994aa9d7fbd473b23..ef616225fe3d8bf827325982937a3c13428f2283 100644
--- a/R/CST_Anomaly.R
+++ b/R/CST_Anomaly.R
@@ -18,7 +18,7 @@
 #'
 #' @return A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools.
 #'
-#'@import s2dverification
+#'@importFrom s2dverification Clim Ano_CrossValid InsertDim
 #'
 #'@examples
 #'# Example 1:
@@ -165,12 +165,12 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
       clim_exp <- tmp$clim_exp
       clim_obs <- tmp$clim_obs
     } else {
-      clim_exp <- InsertDim(tmp$clim_exp, 2, dim_exp[2]) 
-      clim_obs <- InsertDim(tmp$clim_obs, 2, dim_obs[2]) 
+      clim_exp <- s2dverification::InsertDim(tmp$clim_exp, 2, dim_exp[2]) 
+      clim_obs <- s2dverification::InsertDim(tmp$clim_obs, 2, dim_obs[2]) 
     }
    
-    clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) 
-    clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) 
+    clim_exp <- s2dverification::InsertDim(clim_exp, 3, dim_exp[3]) 
+    clim_obs <- s2dverification::InsertDim(clim_obs, 3, dim_obs[3]) 
     ano <- NULL    
     ano$ano_exp <- exp$data - clim_exp
     ano$ano_obs <- obs$data - clim_obs 
diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R
index 1da7fb5be8a8679d022a7d5179dcbdc2830acd26..9baf897bb680eff59d5b3bd4b3169e4a880ee65b 100644
--- a/R/CST_BiasCorrection.R
+++ b/R/CST_BiasCorrection.R
@@ -11,7 +11,6 @@
 #'
 #'@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)
 #'
-#'@import s2dverification
 #'@import multiApply
 #'@examples
 #'
diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R
index ca29039764d5cbe108f1c340cf9188963972141e..fd381790f936697fff109766e69b4aa847a77b73 100644
--- a/R/CST_Calibration.R
+++ b/R/CST_Calibration.R
@@ -13,7 +13,7 @@
 #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
 #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object.
 #'
-#'@import s2dverification
+#'@importFrom s2dverification InsertDim
 #'@import abind
 #'
 #'@seealso \code{\link{CST_Load}}
@@ -84,7 +84,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min",
 #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
 #'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array.
 #'
-#'@import s2dverification
+#'@importFrom s2dverification InsertDim
 #'@import abind
 #'
 #'@seealso \code{\link{CST_Load}}
@@ -251,7 +251,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 
 .calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations
   amt.mbr <- dim(fc)[1]
-  obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr)
+  obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr)
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs")
   obs.av <- mean(obs, na.rm = TRUE)
@@ -271,7 +271,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 
 .calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations
   amt.mbr <- dim(fc)[1]
-  obs.per.ens <- InsertDim(var = obs, posdim = 1, lendim = amt.mbr)
+  obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr)
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs")
   obs.av <- mean(obs, na.rm = TRUE)
@@ -296,7 +296,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE)
   fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE)
-  fc.ens.av.per.ens <- InsertDim(var = fc.ens.av, posdim = 1, lendim = amt.mbr)
+  fc.ens.av.per.ens <- s2dverification::InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr)
   fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE)
   fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE))
   fc.dev <- fc - fc.ens.av.per.ens
@@ -322,10 +322,10 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 .calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts
 
   amt.mbr <- dim(fc)[1]
-  repmat1.tmp <- InsertDim(var = fc, posdim = 1, lendim = amt.mbr)
+  repmat1.tmp <- s2dverification::InsertDim(fc, posdim = 1, lendim = amt.mbr)
   repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3))
   spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE)
-  spr.abs.per.ens <- InsertDim(var = spr.abs, posdim = 1, lendim = amt.mbr)
+  spr.abs.per.ens <- s2dverification::InsertDim(spr.abs, posdim = 1, lendim = amt.mbr)
 
   return(
     append(.calc.fc.quant(fc),
diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R
index 5b70d068449d58a4c1f56b00484d7ce06b604a1f..26ac54492112884c660ec3f4495aaed0205a8090 100644
--- a/R/CST_CategoricalEnsCombination.R
+++ b/R/CST_CategoricalEnsCombination.R
@@ -55,7 +55,7 @@
 #'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.
 #'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).
 #'
-#'@import s2dverification
+#'@importFrom s2dverification InsertDim
 #'@import abind
 #'@examples 
 #'
@@ -89,7 +89,7 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me
   exp$data <- .CategoricalEnsCombination.wrap(fc = exp$data, obs = obs$data, cat.method = cat.method, eval.method = eval.method, amt.cat = amt.cat, ...)
   names.dim.tmp[which(names.dim.tmp == "member")] <- "category"
   names(dim(exp$data)) <- names.dim.tmp
-  exp$data <- InsertDim(var = exp$data, lendim = 1, posdim = 2)
+  exp$data <- s2dverification::InsertDim(exp$data, lendim = 1, posdim = 2)
   names(dim(exp$data))[2] <- "member"
   exp$Datasets <- c(exp$Datasets, obs$Datasets)
   exp$source_files <- c(exp$source_files, obs$source_files)
@@ -277,8 +277,8 @@ comb.dims <- function(arr.in, dims.to.combine){
                           ui = constr.mtrx, ci = constr.vec,
                           freq.per.mdl.at.obs = freq.per.mdl.at.obs) 
       init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01)))          
-      var.cat.fc[ , eval.dexes] <- apply(InsertDim(var = 
-        InsertDim(var = optim.tmp$par, lendim = amt.cat, posdim = 2),
+      var.cat.fc[ , eval.dexes] <- apply(s2dverification::InsertDim( 
+        s2dverification::InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2),
         lendim =  amt.sdate.ev, posdim = 3) *
         freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE)
     } else if (cat.method == "comb") {
@@ -363,7 +363,8 @@ comb.dims <- function(arr.in, dims.to.combine){
 
 .funct.optim.grad <- function(par, freq.per.mdl.at.obs){
   amt.model <- dim(freq.per.mdl.at.obs)[1]
-  return(-apply(freq.per.mdl.at.obs/InsertDim(var = drop(par %*% freq.per.mdl.at.obs),
+  return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% 
+freq.per.mdl.at.obs),
    lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE))
 }
 
@@ -373,7 +374,7 @@ comb.dims <- function(arr.in, dims.to.combine){
   amt.mdl <- mdl.feat$amt.mdl
   mdl.msk.tmp <- mdl.feat$mdl.msk
   amt.coeff <- amt.mdl + 1
-  msk.fc.obs <- (cat.fc == InsertDim(var = cat.obs, posdim = 1, lendim = amt.mbr))
+  msk.fc.obs <- (cat.fc == s2dverification::InsertDim(cat.obs, posdim = 1, lendim = amt.mbr))
   freq.per.mdl.at.obs <- array(NA, c(amt.coeff, amt.sdate))
   for (i.mdl in seq(1, amt.mdl)){
     freq.per.mdl.at.obs[i.mdl, ] <- apply(msk.fc.obs[mdl.msk.tmp[i.mdl, ], , drop = FALSE],
diff --git a/R/CST_Load.R b/R/CST_Load.R
index 76d5deb3335941151a227f174122027828eb7eee..65b695cdcf11a3bf0e5ecb06a7a36c24c4cec1d2 100644
--- a/R/CST_Load.R
+++ b/R/CST_Load.R
@@ -9,7 +9,7 @@
 #' @param ... Parameters that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`.
 #' @return A list with one or two S3 objects, named 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational data, respectively. These 's2dv_cube's can be ingested by other functions in CSTools. If the parameter `exp` in the call to `CST_Load` is set to `NULL`, then only the 'obs' component is returned, and viceversa.
 #' @author Nicolau Manubens, \email{nicolau.manubens@bsc.es}
-#' @import s2dverification
+#' @importFrom s2dverification Load
 #' @importFrom utils glob2rx
 #' @export
 #' @examples
diff --git a/R/CST_MergeDims.R b/R/CST_MergeDims.R
index a9923c581b95e6c48c37d2efa2e32d7588a947f0..720448c7fc02d4fe7fea20655d761473654275d6 100644
--- a/R/CST_MergeDims.R
+++ b/R/CST_MergeDims.R
@@ -10,7 +10,7 @@
 #'@param na.rm a logical indicating if the NA values should be removed or not.
 #'
 #'@import abind
-#'@import s2dverification
+#'@importFrom s2dverification Subset
 #'@examples
 #'
 #'data <- 1 : c(2 * 3 * 4 * 5 * 6 * 7)
@@ -49,7 +49,7 @@ CST_MergeDims <- function(data, merge_dims = c('ftime', 'monthly'), rename_dim =
 #'@param na.rm a logical indicating if the NA values should be removed or not.
 #'
 #'@import abind
-#'@import s2dverification
+#'@importFrom s2dverification Subset
 #'@examples
 #'
 #'data <- 1 : 20
diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R
index e85c8b68d7b205e492dac752d8186f7a5ec85215..d25ed5e49d888071e13755bcad69c93a6d5c2843 100644
--- a/R/CST_MultiMetric.R
+++ b/R/CST_MultiMetric.R
@@ -4,18 +4,24 @@
 #'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
 #'@description This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations.
 #'
-#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiment data in the element named \code{$data}.
+#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiments data in the element named \code{$data}.
 #'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}.
-#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.
+#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms', 'rmsss' and 'rpss'. If 'rpss' is chossen the terciles probabilities are evaluated.
 #'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed.
 #'
-#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0.
-#'@seealso \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}}
+#'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.
+#'@param memb_dim name of the member dimension. It can be NULL, the default value is 'member'.
+#'@param sdate_dim name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.
+#'@return an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the first position in the first 'nexp' dimension correspons to the Multi-Model Mean. 
+#'@seealso \code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}}
 #'@references 
 #'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z}
 #' 
-#'@import s2dverification
+#'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim
+#'@import abind
+#'@importFrom easyVerification climFairRpss veriApply
 #'@import stats
+#'@import multiApply
 #'@examples
 #'library(zeallot)
 #'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
@@ -31,21 +37,59 @@
 #'c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
 #'a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs)
 #'str(a)
+#'\donttest{
+#'exp <- lonlat_data$exp
+#'obs <- lonlat_data$obs
+#'a <- CST_MultiMetric(exp, obs, metric = 'rpss', multimodel = FALSE)
+#'a <- CST_MultiMetric(exp, obs, metric = 'correlation')
+#'a <- CST_MultiMetric(exp, obs, metric = 'rms')
+#'a <- CST_MultiMetric(exp, obs, metric = 'rmsss')
+#'}
 #'@export
-CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) {
+CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE,
+                            time_dim = 'ftime', memb_dim = 'member',
+                            sdate_dim = 'sdate') {
   if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) {
     stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
   }
-
-  if (dim(obs$data)['member'] != 1) {
-    stop("The length of the dimension 'member' in the component 'data' ",
-         "of the parameter 'obs' must be equal to 1.")
-  }
-
-  if (!is.null(names(dim(exp$data))) & !is.null(names(dim(obs$data)))) {
-    if (all(names(dim(exp$data)) %in% names(dim(obs$data)))) {
-      dimnames <- names(dim(exp$data))
+  result <- MultiMetric(exp$data, obs$data, metric = metric, multimodel = multimodel,
+                        time_dim = time_dim, memb_dim = memb_dim, sdate_dim = sdate_dim)
+  exp$data <- result
+  return(exp)
+}
+#'Multiple Metrics applied in Multiple Model Anomalies
+#'
+#'@author Mishra Niti, \email{niti.mishra@bsc.es}
+#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+#'@description This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations on arrays with named dimensions.
+#'
+#'@param exp a multidimensional array with named dimensions.
+#'@param obs a multidimensional array with named dimensions.
+#'@param metric a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.
+#'@param multimodel a logical value indicating whether a Multi-Model Mean should be computed.
+#'
+#'@param time_dim name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.
+#'@param memb_dim name of the member dimension. It can be NULL, the default value is 'member'.
+#'@param sdate_dim name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.
+#'@return a list of arrays containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean. 
+#'@seealso \code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}}
+#'@references 
+#'Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z}
+#' 
+#'@importFrom s2dv MeanDims Reorder Corr RMS RMSSS InsertDim
+#'@import abind
+#'@importFrom easyVerification climFairRpss veriApply
+#'@import stats
+#'@import multiApply
+#'@examples
+#'res <- MultiMetric(lonlat_data$exp$data, lonlat_data$obs$data)
+#'@export
+MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE,
+                        time_dim = 'ftime', memb_dim = 'member', sdate_dim = 'sdate') {
+  if (!is.null(names(dim(exp))) & !is.null(names(dim(obs)))) {
+    if (all(names(dim(exp)) %in% names(dim(obs)))) {
+      dimnames <- names(dim(exp))
     } else {
       stop("Dimension names of element 'data' from parameters 'exp'",
            " and 'obs' should have the same name dimmension.")
@@ -54,7 +98,6 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE)
     stop("Element 'data' from parameters 'exp' and 'obs'",
          " should have dimmension names.")
   }
-  
   if (!is.logical(multimodel)) {
     stop("Parameter 'multimodel' must be a logical value.")
   }
@@ -68,49 +111,110 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE)
     warning("Parameter 'multimodel' has length > 1 and only the first ",
             "element will be used.")
   }
-  
-  # seasonal average of anomalies per model
-  AvgExp <- MeanListDim(exp$data, narm = T, c(2, 4))
-  AvgObs <- MeanListDim(obs$data, narm = T, c(2, 4))
-  
-  # indv model correlation
-  if (metric == 'correlation') {
-    corr <- Corr(AvgExp, AvgObs, posloop = 1, poscor = 2)
-  } else if (metric == 'rms') {
-    corr <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2)
-  } else if (metric == 'rmsss') {
-    corr <- RMSSS(AvgExp, AvgObs, posloop = 1, posRMS = 2)
+  if (is.null(time_dim) | !is.character(time_dim)) {
+    time_dim <- 'time'
+  }
+  if (is.null(memb_dim) | !is.character(memb_dim)) {
+    memb_dim <- 'memb'
+  }
+  if( is.null(sdate_dim) | !is.character(sdate_dim)) {
+    sdate_dim <- 'sdate'
+  }
+  exp_dims <- dim(exp)
+  obs_dims <- dim(obs) 
+  if (!is.null(names(exp_dims)) & !is.null(names(obs_dims))) {
+    if (all(names(exp_dims) == names(obs_dims))) {
+      if (!(time_dim %in% names(exp_dims))) {
+        warning("Parameter 'time_dim' does not match with a dimension name in 'exp'",
+                " and 'obs'. A 'time_dim' of length 1 is added.")
+        dim(exp) <- c(exp_dims, time_dim = 1)
+        names(dim(exp))[length(dim(exp))] <- time_dim
+        dim(obs) <- c(obs_dims, time_dim = 1)
+        names(dim(obs))[length(dim(obs))] <- time_dim
+        exp_dims <- dim(exp)
+        obs_dims <- dim(obs) 
+      }
+      if (!(memb_dim %in% names(exp_dims))) {
+        warning("Parameter 'memb_dim' does not match with a dimension name in ",
+                "'exp' and 'obs'. A 'memb_dim' of length 1 is added.")
+        dim(exp) <- c(exp_dims, memb_dim = 1)
+        names(dim(exp))[length(dim(exp))] <- memb_dim
+        dim(obs) <- c(obs_dims, memb_dim = 1)
+        names(dim(obs))[length(dim(obs))] <- memb_dim
+        exp_dims <- dim(exp)
+        obs_dims <- dim(obs)
+      }
+      if (!(sdate_dim %in% names(exp_dims))) {
+        warning("Parameter 'sdate_dim' does not match with a dimension name in ",
+                "'exp' and 'obs'. A 'sdate_dim' of length 1 is added.")
+        dim(exp) <- c(exp_dims, sdate_dim = 1)
+        names(dim(exp))[length(dim(exp))] <- sdate_dim
+        dim(obs) <- c(obs_dims, sdate_dim = 1)
+        names(dim(obs))[length(dim(obs))] <- sdate_dim
+        exp_dims <- dim(exp)
+        obs_dims <- dim(obs)
+      }
+    } else {
+      stop("Dimension names of element 'data' from parameters 'exp'",
+           " and 'obs' should be the same and in the same order.")
+    }
   } else {
-    stop("Parameter 'metric' must be a character string indicating ",
-         "one of the options: 'correlation', 'rms' or 'rmse'.")
+    stop("Element 'data' from parameters 'exp' and 'obs'",
+         " should have dimmension names.")
   }
-  if (multimodel == TRUE) {
-    # seasonal avg of anomalies for multi-model
-    AvgExp_MMM <- MeanListDim(AvgExp, narm = TRUE, 1)
-    AvgObs_MMM <- MeanListDim(AvgObs, narm = TRUE, 1)
-    # multi model correlation
+  if (metric == 'rpss') {
+    if (multimodel == TRUE) {
+      warning("A probabilistic metric cannot be use to evaluate a multimodel mean.")
+    }
+    AvgExp <- MeanDims(exp, time_dim, na.rm = TRUE)
+    AvgObs <- MeanDims(obs, time_dim, na.rm = TRUE)
+    dif_dims <- which(dim(AvgExp) != dim(AvgObs))
+    dif_dims <- names(dif_dims[-which(names(dif_dims) == memb_dim)])
+    lapply(dif_dims, function(x) {
+      names(dim(AvgExp))[which(names(dim(AvgExp)) == x)] <<- paste0(dif_dims, '_exp')})
+
+    pos_memb <- which(names(dim(AvgExp)) == memb_dim)
+    dim(AvgObs) <- dim(AvgObs)[-pos_memb]
+    AvgExp <- Reorder(AvgExp, c(names(dim(AvgExp))[-pos_memb], memb_dim))
+    pos_memb <- which(names(dim(AvgExp)) == memb_dim)
+    pos_sdate <- which(names(dim(AvgExp)) == sdate_dim)
+    corr <- Apply(list(AvgExp, AvgObs),
+       target_dims = list(c(sdate_dim, 'lat', 'lon', memb_dim), c(sdate_dim, 'lat', 'lon')),
+                  fun = function(x, y) {
+                    veriApply('FairRpss', fcst = x, obs = y,
+                              ensdim = which(names(dim(x)) ==  'member'),
+                              tdim = which(names(dim(x)) == 'sdate'),
+                              prob = c(1/3, 2/3))})
+    
+  } else if (metric %in% c('correlation', 'rms', 'rmsss')) {
+    AvgExp <- MeanDims(exp, c(memb_dim, time_dim), na.rm = TRUE)
+    AvgObs <- MeanDims(obs, c(memb_dim, time_dim), na.rm = TRUE)
+    dataset_dim <- c('data', 'dataset', 'datsets', 'models')
+    if (any(dataset_dim %in% names(exp_dims))) {
+        dataset_dim <- dataset_dim[dataset_dim %in% names(dim(AvgExp))]
+    } else {
+        warning("Parameter 'exp' and 'obs' does not have a dimension 'dataset'.")
+    }
+    if (multimodel == TRUE) {
+        # seasonal avg of anomalies for multi-model
+        AvgExp_MMM <- MeanDims(AvgExp, c(dataset_dim), na.rm = TRUE)
+        pos_dataset <- which(names(dim(AvgExp)) == dataset_dim)
+        AvgExp_MMM <- s2dv::InsertDim(AvgExp_MMM, posdim = pos_dataset, lendim = 1,
+                                        name = dataset_dim)
+        AvgExp <- abind(AvgExp_MMM, AvgExp, along = pos_dataset)
+        names(dim(AvgExp)) <- names(dim(AvgExp_MMM))
+    }
     if (metric == 'correlation') {
-      corr_MMM <- Corr(var_exp = InsertDim(AvgExp_MMM, 1, 1),
-                       var_obs = InsertDim(AvgObs_MMM, 1, 1),
-                       posloop = 1, poscor = 2)
+      corr <- s2dv::Corr(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
     } else if (metric == 'rms') {
-      corr_MMM <- RMS(var_exp = InsertDim(AvgExp_MMM, 1, 1),
-                      var_obs = InsertDim(AvgObs_MMM, 1, 1),
-                      posloop = 1, posRMS = 2)
+      corr <- s2dv::RMS(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
     } else if (metric == 'rmsss') {
-      corr_MMM <- RMSSS(var_exp = InsertDim(AvgExp_MMM, 1, 1),
-                        var_obs = InsertDim(AvgObs_MMM, 1, 1), 
-                        posloop = 1, posRMS = 2)
-    } 
-    corr <- abind::abind(corr, corr_MMM, along = 1)
+      corr <- s2dv::RMSSS(AvgExp, AvgObs, dat_dim = dataset_dim, time_dim = sdate_dim)
+    }
+  } else {
+    stop("Parameter 'metric' must be a character string indicating ",
+         "one of the options: 'correlation', 'rms', 'rmsss' or 'rpss'.")
   }
   
-  names(dim(corr)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6])
-  
-  #exp$data <- ano$ano_exp
-  #obs$data <- ano$ano_obs
-  exp$data <- corr
-  exp$Datasets <- c(exp$Datasets, obs$Datasets)
-  exp$source_files <- c(exp$source_files, obs$source_files)
-  return(exp)
+  return(corr)
 }
diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R
index c5ab53b5291a690d50a139b317db39b77b3ccad7..70c88c09c8d88c6ed785e4e60eb2b6fd35d56080 100644
--- a/R/CST_MultivarRMSE.R
+++ b/R/CST_MultivarRMSE.R
@@ -10,8 +10,8 @@
 #'
 #'@return an object of class \code{s2dv_cube} containing the RMSE in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs.  An array with dimensions: c(number of exp, number of obs, 1 (the multivariate RMSE value), number of lat, number of lon)
 #'
-#'@seealso \code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}}
-#'@import s2dverification
+#'@seealso \code{\link[s2dv]{RMS}} and \code{\link{CST_Load}}
+#'@importFrom s2dv RMS MeanDims
 #'@examples
 #'# Creation of sample s2dverification objects. These are not complete
 #'# s2dverification objects though. The Load function returns complete objects.
@@ -108,17 +108,18 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) {
   sumweights <- 0
   for (j in 1 : nvar) {
     # seasonal average of anomalies
-    AvgExp <- MeanListDim(exp[[j]]$data, narm = TRUE, c(2, 4))
-    AvgObs <- MeanListDim(obs[[j]]$data, narm = TRUE, c(2, 4))
+    AvgExp <- MeanDims(exp[[j]]$data, c('member', 'ftime'), na.rm = TRUE)
+    AvgObs <- MeanDims(obs[[j]]$data, c('member', 'ftime'), na.rm = TRUE)
     # multivariate RMSE (weighted) 
-    rmse <- RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, conf = FALSE)
+    rmse <- s2dv::RMS(AvgExp, AvgObs, dat_dim = 'dataset', time_dim = 'sdate',
+                conf = FALSE)$rms
     stdev <- sd(AvgObs)
     mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j]))
     sumweights <- sumweights + as.numeric(weight[j])
   }
   mvrmse <- mvrmse / sumweights 
   
-  names(dim(mvrmse)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6])
+ # names(dim(mvrmse)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6])
   exp_Datasets <- unlist(lapply(exp, function(x) {
                          x[[which(names(x) == 'Datasets')]]}))
   exp_source_files <- unlist(lapply(exp, function(x) {
diff --git a/R/CST_RegimesAssign.R b/R/CST_RegimesAssign.R
index 60f0d7045a622bb70526b7f7248b12130cd87637..72ae69f1e16ce301f194794d55ad5f9e55251201 100644
--- a/R/CST_RegimesAssign.R
+++ b/R/CST_RegimesAssign.R
@@ -27,7 +27,7 @@
 #'         that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'),
 #'         \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) ,
 #'         \code{$frequency} (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.),
-#'@import s2dverification
+#'@importFrom s2dverification ACC Mean1Dim InsertDim
 #'@import multiApply
 #'@examples
 #'\dontrun{
@@ -102,7 +102,7 @@ CST_RegimesAssign <- function(data,  ref_maps,
 #'         \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) ,
 #'         \code{$frequency} (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.),
 #'  
-#'@import s2dverification
+#'@importFrom s2dverification ACC Mean1Dim Eno InsertDim
 #'@import multiApply   
 #'@examples 
 #'\dontrun{
@@ -190,9 +190,9 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
     
     if (any(is.na(index))) {
       recon <-list(
-          composite = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
+          composite = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
                                 postime, dim(ref_maps)['composite.cluster']),
-          pvalue = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
+          pvalue = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
                              postime, dim(ref_maps)['composite.cluster']))
     } else {
       if (memb) {
@@ -255,7 +255,7 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
                     which(names(dim(target)) == 'lon')))
   
   # weights are defined
-  latWeights <- InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3])
+  latWeights <- s2dverification::InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3])
   
   
   rmsdiff <- function(x, y) {
@@ -278,11 +278,11 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
     corr <- rep(NA, nclust)
     for (i in 1:nclust) {
       corr[i] <-
-        ACC(InsertDim(InsertDim(
-          InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1
+        ACC(s2dverification::InsertDim(InsertDim(
+          s2dverification::InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1
         ), 3, 1),
-        InsertDim(InsertDim(
-          InsertDim(target * latWeights, 1, 1), 2, 1
+        s2dverification::InsertDim(InsertDim(
+          s2dverification::InsertDim(target * latWeights, 1, 1), 2, 1
         ), 3, 1))$ACC[2]
     }
     assign <- which(corr == max(corr))
diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R
index 1bdd4779b6ae20158815a660e7ff61cc2f8845e4..4ae560bff5757c03a91a3e251a9442dca13ed2fe 100644
--- a/R/CST_SaveExp.R
+++ b/R/CST_SaveExp.R
@@ -18,6 +18,7 @@
 #'
 #'@import ncdf4
 #'@importFrom s2dv Reorder
+#'@importFrom s2dverification InsertDim
 #'@import multiApply
 #'
 #'@examples
@@ -79,6 +80,7 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
 #' 
 #'@import ncdf4
 #'@importFrom s2dv Reorder
+#'@importFrom s2dverification InsertDim
 #'@import multiApply
 #'
 #'@examples
@@ -142,7 +144,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
   if (length(dataset_pos) == 0) {
     warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ",
             "All data is stored in the same 'dataset' folder.")
-    data$data <- InsertDim(var = data, posdim = 1, lendim = 1)
+    data$data <- s2dverification::InsertDim(data, posdim = 1, lendim = 1)
     names(dim(data))[1] <- "dataset"
     dimname <- c("dataset", dimname)
     dataset_pos = 1
diff --git a/R/CST_SplitDim.R b/R/CST_SplitDim.R
index 3326d435af6ecf535bc66ebe09d7d05a14cd499b..6332684f75ebc404f2c772f9cd0a3edfebcb8482 100644
--- a/R/CST_SplitDim.R
+++ b/R/CST_SplitDim.R
@@ -11,7 +11,7 @@
 #'@param new_dim_name a character string indicating the name of the new dimension.
 #'
 #'@import abind
-#'@import s2dverification
+#'@importFrom s2dverification Subset
 #'@examples
 #'
 #'data <- 1 : 20
@@ -73,7 +73,7 @@ CST_SplitDim <- function(data, split_dim = 'time', indices = NULL,
 #'@param freq a character string indicating the frequency: by 'day', 'month' and 'year' or 'monthly' (by default). 'month' identifies months between 1 and 12 independetly of the year they belong to, while 'monthly' differenciates months from different years. Parameter 'freq' can also be numeric indicating the length in which to subset the dimension.
 #'@param new_dim_name a character string indicating the name of the new dimension.
 #'@import abind
-#'@import s2dverification
+#'@importFrom s2dverification Subset
 #'@examples
 #'
 #'data <- 1 : 20
diff --git a/R/CST_WeatherRegimes.R b/R/CST_WeatherRegimes.R
index 72ab3987e38e2ab4b7f45bcc2b10016f3eecbb92..e7cc925c3d9bc9d11180c33dbd01ab722aecafdc 100644
--- a/R/CST_WeatherRegimes.R
+++ b/R/CST_WeatherRegimes.R
@@ -34,7 +34,7 @@
 #'          \code{cluster} (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.),
 #'         \code{persistence} (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)),
 #'         \code{frequency} (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).),
-#'@import s2dverification
+#'@importFrom s2dverification EOF
 #'@import multiApply
 #'@examples
 #'\dontrun{
@@ -109,7 +109,7 @@ CST_WeatherRegimes <- function(data, ncenters = NULL,
 #'         \code{cluster} (A matrix or vector with integers (from 1:k) indicating the cluster to which each time step is allocated.),
 #'         \code{persistence} (Percentage of days in a month/season before a cluster is replaced for a new one (only if method=’kmeans’ has been selected.)),
 #'         \code{frequency} (Percentage of days in a month/season belonging to each cluster (only if method=’kmeans’ has been selected).),
-#'@import s2dverification
+#'@importFrom s2dverification EOF
 #'@import multiApply
 #'@examples
 #'\dontrun{
diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R
index 7f8f5d64e0191ab40cd7e88b25d23fd95f9269e9..2aee71b8ef7a3861a3834e05f0110ea10bc1751d 100644
--- a/R/PlotCombinedMap.R
+++ b/R/PlotCombinedMap.R
@@ -25,7 +25,7 @@
 
 #'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap}
 #' 
-#'@import s2dverification
+#'@importFrom s2dv PlotEquiMap ColorBar
 #'@importFrom maps map 
 #'@importFrom graphics box image layout mtext par plot.new
 #'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff
diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R
index 9cd4dfebcad4c7b037daacd4d681ed534289f55e..8c4b1068b9e96893018f29c378c04ec30fdae71b 100644
--- a/R/PlotForecastPDF.R
+++ b/R/PlotForecastPDF.R
@@ -23,8 +23,7 @@
 #'@importFrom reshape2 melt
 #'@importFrom plyr .
 #'@importFrom plyr dlply
-#'@import s2dverification
-#'
+#'@importFrom s2dverification InsertDim
 #'@examples
 #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), 
 #'                    fcst3 = rnorm(10, -0.5, 0.9))
@@ -109,7 +108,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         if (length(tercile.limits) != 2) {
             stop("Provide two tercile limits")
         } 
-        tercile.limits <- InsertDim(tercile.limits, 1, npanels)
+        tercile.limits <- s2dverification::InsertDim(tercile.limits, 1, npanels)
     } else if (is.array(tercile.limits)) {
         if (length(dim(tercile.limits)) == 2) {
             if (dim(tercile.limits)[2] != 2) {
@@ -136,7 +135,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
             if (length(extreme.limits) != 2) {
                 stop("Provide two extreme limits")
             } 
-            extreme.limits <- InsertDim(extreme.limits, 1, npanels)
+            extreme.limits <- s2dverification::InsertDim(extreme.limits, 1, npanels)
         } else if (is.array(extreme.limits)) {
             if (length(dim(extreme.limits)) == 2) {
                 if (dim(extreme.limits)[2] != 2) {
diff --git a/R/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R
index 57be86435ddf526ebe051cb8b2cfd3fd4c394fa6..b4e974a88a70b39fac07e339e1fe7bdae7fe65c6 100644
--- a/R/PlotMostLikelyQuantileMap.R
+++ b/R/PlotMostLikelyQuantileMap.R
@@ -12,7 +12,6 @@
 #'@param ... additional parameters to be sent to \code{PlotCombinedMap} and \code{PlotEquiMap}.
 #'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap}
 #' 
-#'@import s2dverification
 #'@importFrom maps map 
 #'@importFrom graphics box image layout mtext par plot.new
 #'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff
diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R
index fcfb36bbb68b7b58019634d9417ce1ad0609523a..14d578a54b912c6a358cb1ac913ce63901f8c2dc 100644
--- a/R/PlotTriangles4Categories.R
+++ b/R/PlotTriangles4Categories.R
@@ -72,6 +72,7 @@
 #'@importFrom grDevices dev.new dev.off dev.cur 
 #'@importFrom graphics plot points polygon text title axis
 #'@importFrom RColorBrewer brewer.pal
+#'@importFrom s2dv ColorBar
 #'@export
 PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL,
                                      toptitle = NULL, sig_data = NULL,
diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd
index 8e3ce593b7024138e191c4fc818516358503be7d..e8e047dc6f6fba4e4dcaf6f30af4be5e3e0f4636 100644
--- a/man/CST_MultiMetric.Rd
+++ b/man/CST_MultiMetric.Rd
@@ -4,19 +4,33 @@
 \alias{CST_MultiMetric}
 \title{Multiple Metrics applied in Multiple Model Anomalies}
 \usage{
-CST_MultiMetric(exp, obs, metric = "correlation", multimodel = TRUE)
+CST_MultiMetric(
+  exp,
+  obs,
+  metric = "correlation",
+  multimodel = TRUE,
+  time_dim = "ftime",
+  memb_dim = "member",
+  sdate_dim = "sdate"
+)
 }
 \arguments{
-\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiment data in the element named \code{$data}.}
+\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of the seasonal forecast experiments data in the element named \code{$data}.}
 
 \item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the anomaly of observed data in the element named \code{$data}.}
 
-\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.}
+\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms', 'rmsss' and 'rpss'. If 'rpss' is chossen the terciles probabilities are evaluated.}
 
 \item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.}
+
+\item{time_dim}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.}
+
+\item{memb_dim}{name of the member dimension. It can be NULL, the default value is 'member'.}
+
+\item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.}
 }
 \value{
-an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is an array with two datset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest first dimension correspons to the Multi-Model Mean. The third dimension contains the statistics selected. For metric \code{correlation}, the third dimension is of length four and they corresponds to the lower limit of the 95\% confidence interval, the statistics itselfs, the upper limit of the 95\% confidence interval and the 95\% significance level. For metric \code{rms}, the third dimension is length three and they corresponds to the lower limit of the 95\% confidence interval, the RMSE and the upper limit of the 95\% confidence interval. For metric \code{rmsss}, the third dimension is length two and they corresponds to the statistics itselfs and the p-value of the one-sided Fisher test with Ho: RMSSS = 0.
+an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean.
 }
 \description{
 This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations.
@@ -36,12 +50,20 @@ attr(obs, 'class') <- 's2dv_cube'
 c(ano_exp, ano_obs) \%<-\% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
 a <- CST_MultiMetric(exp = ano_exp, obs = ano_obs)
 str(a)
+\donttest{
+exp <- lonlat_data$exp
+obs <- lonlat_data$obs
+a <- CST_MultiMetric(exp, obs, metric = 'rpss', multimodel = FALSE)
+a <- CST_MultiMetric(exp, obs, metric = 'correlation')
+a <- CST_MultiMetric(exp, obs, metric = 'rms')
+a <- CST_MultiMetric(exp, obs, metric = 'rmsss')
+}
 }
 \references{
 Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z}
 }
 \seealso{
-\code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link{CST_Load}}
+\code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}}
 }
 \author{
 Mishra Niti, \email{niti.mishra@bsc.es}
diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd
index 24af608c8360985de73763434a16f1a11fd35ffe..c318c105ee312313426535a35152684e4479a5fe 100644
--- a/man/CST_MultivarRMSE.Rd
+++ b/man/CST_MultivarRMSE.Rd
@@ -57,7 +57,7 @@ a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight)
 str(a)
 }
 \seealso{
-\code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}}
+\code{\link[s2dv]{RMS}} and \code{\link{CST_Load}}
 }
 \author{
 Deborah Verfaillie, \email{deborah.verfaillie@bsc.es}
diff --git a/man/MultiMetric.Rd b/man/MultiMetric.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f688f0a31050bb7bf3be5735454b3b1e41a263af
--- /dev/null
+++ b/man/MultiMetric.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_MultiMetric.R
+\name{MultiMetric}
+\alias{MultiMetric}
+\title{Multiple Metrics applied in Multiple Model Anomalies}
+\usage{
+MultiMetric(
+  exp,
+  obs,
+  metric = "correlation",
+  multimodel = TRUE,
+  time_dim = "ftime",
+  memb_dim = "member",
+  sdate_dim = "sdate"
+)
+}
+\arguments{
+\item{exp}{a multidimensional array with named dimensions.}
+
+\item{obs}{a multidimensional array with named dimensions.}
+
+\item{metric}{a character string giving the metric for computing the maximum skill. This must be one of the strings 'correlation', 'rms' or 'rmsss.}
+
+\item{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.}
+
+\item{time_dim}{name of the temporal dimension where a mean will be applied. It can be NULL, the default value is 'ftime'.}
+
+\item{memb_dim}{name of the member dimension. It can be NULL, the default value is 'member'.}
+
+\item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.}
+}
+\value{
+a list of arrays containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean.
+}
+\description{
+This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations on arrays with named dimensions.
+}
+\examples{
+res <- MultiMetric(lonlat_data$exp$data, lonlat_data$obs$data)
+}
+\references{
+Mishra, N., Prodhomme, C., & Guemas, V. (n.d.). Multi-Model Skill Assessment of Seasonal Temperature and Precipitation Forecasts over Europe, 29-31.\url{http://link.springer.com/10.1007/s00382-018-4404-z}
+}
+\seealso{
+\code{\link[s2dv]{Corr}}, \code{\link[s2dv]{RMS}}, \code{\link[s2dv]{RMSSS}} and \code{\link{CST_Load}}
+}
+\author{
+Mishra Niti, \email{niti.mishra@bsc.es}
+
+Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+}
diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R
index 1a41a83a53dfb9236d0c0a1eec513973dfd0fdbb..9d048ab3e86d0667c0e9c8a56961283ca6ea0e70 100644
--- a/tests/testthat/test-CST_MultiMetric.R
+++ b/tests/testthat/test-CST_MultiMetric.R
@@ -11,30 +11,38 @@ test_that("basic use case", {
   attr(exp, 'class') <- 's2dv_cube'
   attr(obs, 'class') <- 's2dv_cube'
   
-  result <- list(data = array(rep(c(rep(1, 9), rep(0, 3)), 48), 
-                dim = c(dataset = 3, dataset = 1, statistics = 4, 
-                       lat = 6, lon = 8)),
-                lat = lat, lon = lon) 
+  result <- list(data = list(corr = array(rep(1, 3* 48), 
+                                   dim = c(nexp = 3, nobs = 1,
+                                   lat = 6, lon = 8)),
+                             p.val = array(rep(0, 3 * 48), dim = c(nexp = 3, nobs = 1,
+                                                      lat = 6, lon = 8)),
+                             conf.lower = array(rep(1, 3* 48),
+                                   dim = c(nexp = 3, nobs = 1,     
+                                   lat = 6, lon = 8)),
+                             conf.upper = array(rep(1, 3* 48),
+                                   dim = c(nexp = 3, nobs = 1,     
+                                   lat = 6, lon = 8))),    
+                 lat = lat, lon = lon) 
   attr(result, 'class') <- 's2dv_cube'
   expect_equal(CST_MultiMetric(exp = exp, obs = obs), result)
   
   exp2 <- exp
   exp2$data[1, 1, 1, 2, 1, 1] = NA
-  CST_MultiMetric(exp = exp2, obs = obs) 
   res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms')
   expect_equal(length(res), 3)
-  expect_equal(dim(res$data), 
-        c(dataset = 3, dataset = 1, statistics = 3, lat = 6, lon = 8))
+  expect_equal(dim(res$data$rms), 
+        c(nexp = 3, nobs = 1, lat = 6, lon = 8))
   res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms', 
                          multimodel = FALSE)
-  expect_equal(dim(res$data), 
-        c(dataset = 2, dataset = 1, statistics = 3, lat = 6, lon = 8))
+  expect_equal(dim(res$data$rms), 
+        c(nexp = 2, nobs = 1, lat = 6, lon = 8))
+  expect_equal(length(res$data), 3)
   res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss')
-  expect_equal(dim(res$data), 
-        c(dataset = 3, dataset = 1, statistics = 2, lat = 6, lon = 8))
+  expect_equal(dim(res$data$rmsss), 
+        c(nexp = 3, nobs = 1, lat = 6, lon = 8))
   res <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss', multimodel = FALSE)
-  expect_equal(dim(res$data), 
-        c(dataset = 2, dataset = 1, statistics = 2, lat = 6, lon =8))
+  expect_equal(dim(res$data$rmsss), 
+        c(nexp = 2, nobs = 1, lat = 6, lon = 8))
   })
 
 
@@ -58,7 +66,7 @@ test_that("Sanity checks", {
   expect_error(
     CST_MultiMetric(exp = exp, obs = obs, metric = 1),
     paste0("Parameter 'metric' must be a character string indicating one ",
-           "of the options: 'correlation', 'rms' or 'rmse'"))
+           "of the options: 'correlation', 'rms', 'rmsss' or 'rpss'"))
   expect_error(
     CST_MultiMetric(exp = exp, obs = obs, metric = NA),
     "missing value where TRUE/FALSE needed")
@@ -69,4 +77,15 @@ test_that("Sanity checks", {
     CST_MultiMetric(exp = exp, obs = obs, metric = "correlation", 
                     multimodel = NULL),
     "Parameter 'multimodel' must be a logical value.")
+  expect_error(
+    MultiMetric(exp = lonlat_data$exp, obs = lonlat_data$obs, metric = "rpss",
+                    multimodel = TRUE),
+    "Element 'data' from parameters 'exp' and 'obs' should have dimmension names.")
+exp <- lonlat_data$exp$data[1,,,,,]
+obs <- lonlat_data$obs$data[1,,,,,]
+  expect_error(
+    MultiMetric(exp = exp, obs = obs, metric = "rpss",
+                    multimodel = TRUE),
+  paste0("Dimension names of element 'data' from parameters 'exp' and ",
+         "'obs' should have the same name dimmension."))
 })  
diff --git a/vignettes/MultiModelSkill_vignette.Rmd b/vignettes/MultiModelSkill_vignette.Rmd
index 734652fe67fdd3416825c40a7cf9d81673dd740f..e44be8529629d822878c0cf1987049433fe53533 100644
--- a/vignettes/MultiModelSkill_vignette.Rmd
+++ b/vignettes/MultiModelSkill_vignette.Rmd
@@ -145,9 +145,12 @@ While other relevant data is being stored in the corresponding element of the ob
 
 
 ```r
-> dim(AnomDJF$data)
-   dataset    dataset statistics        lat        lon 
-         4          1          4         35         43 
+> str(AnomDJF$data)
+List of 4
+ $ corr      : num [1:4, 1, 1:35, 1:64] 0.586 0.614 0.143 0.501 0.419 ...
+ $ p.val     : num [1:4, 1, 1:35, 1:64] 0.0026 0.00153 0.26805 0.01036 0.02931 ...
+ $ conf.lower: num [1:4, 1, 1:35, 1:64] 0.2073 0.2485 -0.3076 0.0883 -0.0154 ...
+ $ conf.upper: num [1:4, 1, 1:35, 1:64] 0.812 0.827 0.541 0.767 0.72 ...
 > names(AnomDJF)
 [1] "data"            "lon"             "lat"             "Variable"        "Datasets"        "Dates"          
 [7] "when"            "source_files"    "load_parameters"
@@ -155,26 +158,19 @@ While other relevant data is being stored in the corresponding element of the ob
 [1] "glosea5"                "ecmwf/system4_m1"       "meteofrance/system5_m1" "erainterim"    
 ```
  
-In the element $data of the `AnomDJF` object, the third dimension contains the lower limit of the 95% confidence interval, the correlation, the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test. 
-
-
-```r
-corre <- AnomDJF$data[ , , 2, , ]
-names(dim(corre)) <- c("maps", "lat", "lon")
-```
-
+In the element $data of the `AnomDJF` object is a list of object for the metric and its statistics: correlation, p-value, the lower limit of the 95% confidence interval and the upper limit of the 95% confidence interval and the 95% significance level given by a one-sided T-test. 
 
 To obtain a spatial plot with a scale from -1 to 1 value of correlation for the model with the highest correlation for each grid point, the following lines should be run:
 
 ```r
-PlotCombinedMap(corre, lon = Lon, lat = Lat, map_select_fun = max,
-                display_range = c(0, 1), map_dim = 'maps',
+PlotCombinedMap(AnomDJF$data$corr[,1,,], lon = Lon, lat = Lat, map_select_fun = max,
+                display_range = c(0, 1), map_dim = 'nexp',
                 legend_scale = 0.5, brks = 11,
-                cols = list(c('white', 'darkblue'),
+                cols = list(c('white', 'black'),
+                            c('white', 'darkblue'),
                             c('white', 'darkred'),
-                            c('white', 'darkorange'),
-                            c('white', 'black')),
-                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"),
+                            c('white', 'darkorange')),
+                bar_titles = c("MMM", names(AnomDJF$Datasets)),
                 fileout = "./vignettes/Figures/MultiModelSkill_cor_tas_1992-2012.png",
                 width = 14, height = 8)                
 ```
@@ -200,14 +196,14 @@ The following lines are necessary to obtain the plot which visualizes the best m
 
 ```r
 names(dim(RMS)) <- c("maps", "lat", "lon")
-PlotCombinedMap(RMS, lon = Lon, lat = Lat, map_select_fun = min,
-                display_range = c(0, ceiling(max(abs(RMS)))), map_dim = 'maps',
+PlotCombinedMap(AnomDJF$data$rms[,1,,], lon = Lon, lat = Lat, map_select_fun = min,
+                display_range = c(0, ceiling(max(abs(AnomDJF$data$rms)))), map_dim = 'nexp',
                 legend_scale = 0.5, brks = 11,
-                cols = list(c('darkblue', 'white'),
+                cols = list(c('black', 'white'),
+                            c('darkblue', 'white'),
                             c('darkred', 'white'),
-                            c('darkorange', 'white'),
-                            c('black', 'white')),
-                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), 
+                            c('darkorange', 'white')),
+                bar_titles = c("MMM", names(AnomDJF$Datasets)), 
                 fileout = "./vignettes/Figures/MultiModelSkill_rms_tas_1992-2012.png",
                 width = 14, height = 8)              
 ```
@@ -226,19 +222,17 @@ Notice that the perfect RMSSS is 1 and the parameter  `map_select_fun` from `Plo
 ```r
 AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', 
                           multimodel = TRUE)
-RMSSS <- AnomDJF$data[ , , 1, , ]
-names(dim(RMSSS)) <- c("maps", "lat", "lon")
 
-PlotCombinedMap(RMSSS, lon = Lon, lat = Lat, 
+PlotCombinedMap(AnomDJF$data$rmsss[,1,,], lon = Lon, lat = Lat, 
                 map_select_fun = function(x) {x[which.min(abs(x - 1))]},
                 display_range = c(0, 
-                                  ceiling(max(abs(RMSSS)))), map_dim = 'maps',
+                                  ceiling(max(abs(AnomDJF$data$rmsss)))), map_dim = 'nexp',
                 legend_scale = 0.5, brks = 11,
-                cols = list(c('white', 'darkblue'),
+                cols = list(c('white', 'black'),
+                            c('white', 'darkblue'),
                             c('white', 'darkred'),
-                            c('white', 'darkorange'),
-                            c('white', 'black')),
-                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), 
+                            c('white', 'darkorange')),
+                bar_titles = c("MMM", names(AnomDJF$Datasets)), 
                 fileout = "./vignettes/Figures/MultiModelSkill_rmsss_tas_1992-2012.png",
                 width = 14, height = 8)                
 ```