diff --git a/NAMESPACE b/NAMESPACE
index a2a12c100f95cd4ccbe4b1ff6688c7384e7f940c..872506f75a52612871a1f7a21a1cd06876fec13f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(CST_Anomaly)
 export(CST_BiasCorrection)
 export(CST_Calibration)
 export(CST_Load)
diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R
new file mode 100644
index 0000000000000000000000000000000000000000..e76c46ab67164a77196460bf0e7dbb890d785819
--- /dev/null
+++ b/R/CST_Anomaly.R
@@ -0,0 +1,188 @@
+#'Anomalies relative to a climatology along selected dimension with or without cross-validation
+#'
+#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+#'@author Pena Jesus, \email{jesus.pena@bsc.es}
+#'@description This function computes the anomalies relative to a climatology computed along the 
+#'selected dimension (usually starting dates or forecast time) allowing the application or not of
+#'crossvalidated climatologies. The computation is carried out independently for experimental and 
+#'observational data products.
+#'
+#'@param exp an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.
+#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.'
+#'@param cross A logical value indicating whether cross-validation should be applied or not. Default = FALSE.
+#'@param memb A logical value indicating whether Clim() computes one climatology for each experimental data 
+#'product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE.
+#'
+#'@param dim_anom An integer indicating the dimension along which the climatology will be computed. It 
+#'usually corresponds to 3 (sdates) or 4 (ftime). Default = 3.
+#'
+#' @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
+#'
+#'@examples
+#'# Example 1:
+#'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
+#'dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
+#'obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
+#'dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
+#'lon <- seq(0, 30, 5)
+#'lat <- seq(0, 25, 5)
+#'exp <- list(data = mod, lat = lat, lon = lon)
+#'obs <- list(data = obs, lat = lat, lon = lon)
+#'attr(exp, 'class') <- 's2dv_cube'
+#'attr(obs, 'class') <- 's2dv_cube'
+#'
+#'anom1 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE)
+#'str(anom1)
+#'anom2 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
+#'str(anom2)
+#'
+#'anom3 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = FALSE)
+#'str(anom3)
+#'
+#'anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE)
+#'str(anom4)
+#'
+#'@seealso \code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}}
+#'
+#'
+#'@export
+CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE, dim_anom = 3) {
+ 
+  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 (!is.null(obs) & 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.")
+  }
+
+  case_exp = case_obs = 0
+  if (is.null(exp) & is.null(obs)) {
+    stop("One of the parameter 'exp' or 'obs' cannot be NULL.")
+  }
+  if (is.null(exp)) {
+    exp <- obs 
+    case_obs = 1
+    warning("Parameter 'exp' is not provided and will be recycled.")
+  }  
+  if (is.null(obs)) {
+    obs <- exp 
+    case_exp = 1
+    warning("Parameter 'obs' is not provided and will be recycled.")
+  }
+
+  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))
+    } else {
+      stop("Dimension names of element 'data' from parameters 'exp'",
+           " and 'obs' should have the same name dimmension.")
+    }
+  } else {
+    stop("Element 'data' from parameters 'exp' and 'obs'",
+         " should have dimmension names.")
+  }
+  dim_exp <- dim(exp$data)
+  dim_obs <- dim(obs$data)
+  dimnames_data <- names(dim_exp)
+ 
+  if (dim_exp[dim_anom] == 1 | dim_obs[dim_anom] == 1) {
+    stop("The length of dimension 'dim_anom' in label 'data' of the parameter ",
+         "'exp' and 'obs' must be greater than 1.")
+  }
+  if (!any(names(dim_exp)[dim_anom] == c('sdate', 'time', 'ftime'))) {
+    warning("Parameter 'dim_anom' correspond to a position name different ",
+            "than 'sdate', 'time' or 'ftime'.")
+  } 
+  if (!is.logical(cross) | !is.logical(memb) ) {
+    stop("Parameters 'cross' and 'memb' must be logical.")
+  }
+  if (length(cross) > 1 | length(memb) > 1 ) {
+    cross <- cross[1]
+    warning("Parameter 'cross' has length greater than 1 and only the first element",
+            "will be used.")
+  }
+  if (length(memb) > 1) {
+     memb <- memb[1]
+     warning("Parameter 'memb' has length greater than 1 and only the first element",
+             "will be used.")
+  }
+  
+  # Selecting time dimension through dimensions permutation
+  if (dim_anom != 3) {
+    dimperm <- 1 : length(dim_exp)
+    dimperm[3] <- dim_anom
+    dimperm[dim_anom] <- 3
+    
+    var_exp <- aperm(exp$data, perm = dimperm)
+    var_obs <- aperm(obs$data, perm = dimperm)
+    
+    #Updating permuted dimensions
+    dim_exp <- dim(exp$data)
+    dim_obs <- dim(obs$data)
+  }
+  
+  
+  # Computating anomalies
+  #----------------------
+  
+  # With cross-validation
+  if (cross) {
+    ano <- Ano_CrossValid(var_exp = exp$data, var_obs = obs$data, memb = memb)
+   
+    #  Without cross-validation 
+  } else {
+    tmp <- Clim(var_exp = exp$data, var_obs = obs$data, memb = memb)
+    
+    if (memb) { 
+      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 <- InsertDim(clim_exp, 3, dim_exp[3]) 
+    clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) 
+    ano <- NULL    
+    ano$ano_exp <- exp$data - clim_exp
+    ano$ano_obs <- obs$data - clim_obs 
+  }
+  
+  # Permuting back dimensions to original order
+  if  (dim_anom != 3) {
+    
+    if (case_obs == 0) {
+      ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm)
+    } 
+    if (case_exp == 0) {
+      ano$ano_obs <- aperm(ano$ano_obs, perm = dimperm)
+    }
+        
+    #Updating back permuted dimensions
+    dim_exp <- dim(exp$data)
+    dim_obs <- dim(obs$data)
+  }
+  
+  # Adding dimensions names
+  attr(ano$ano_exp, 'dimensions') <- dimnames_data
+  attr(ano$ano_obs, 'dimensions') <- dimnames_data
+  exp$data <- ano$ano_exp
+  obs$data <- ano$ano_obs
+  
+  #  Outputs
+  # ~~~~~~~~~
+    if (case_obs == 1) {
+      return(obs)
+    } 
+    else if (case_exp == 1) {
+      return(exp)
+    }
+    else {
+    return(list(exp = exp, obs = obs)) 
+  }
+}
diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R
index 596b6152c0de7b406582e4a2b1149536c776b65d..e7bfc729f7b8a87092254fe2bbd32ac8e981e525 100644
--- a/R/CST_MultiMetric.R
+++ b/R/CST_MultiMetric.R
@@ -4,8 +4,8 @@
 #'@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_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.
-#'@param obs an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed 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 experiment 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 multimodel a logical value indicating whether a Multi-Model Mean should be computed.
 #'
@@ -17,6 +17,7 @@
 #'@import s2dverification
 #'@import stats
 #'@examples
+#'library(zeallot)
 #'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
 #'dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
 #'obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
@@ -27,7 +28,8 @@
 #'obs <- list(data = obs, lat = lat, lon = lon)
 #'attr(exp, 'class') <- 's2dv_cube'
 #'attr(obs, 'class') <- 's2dv_cube'
-#'a <- CST_MultiMetric(exp = exp, obs = obs)
+#'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)
 #'@export
 CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) {
@@ -66,11 +68,10 @@ CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE)
     warning("Parameter 'multimodel' has length > 1 and only the first ",
             "element will be used.")
   }
-  ano <- Ano_CrossValid(var_exp = exp$data, var_obs = obs$data)
   
   # seasonal average of anomalies per model
-  AvgExp <- MeanListDim(ano$ano_exp, narm = T, c(2, 4))
-  AvgObs <- MeanListDim(ano$ano_obs, narm = T, c(2, 4))
+  AvgExp <- MeanListDim(exp$data, narm = T, c(2, 4))
+  AvgObs <- MeanListDim(obs$data, narm = T, c(2, 4))
   
   # indv model correlation
   if (metric == 'correlation') {
diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R
index d1fc4f2939454bc2b714e2c733bb598a84ef2977..88967a134aca76502e3eb3bf1a32d18b9bf48217 100644
--- a/R/CST_MultivarRMSE.R
+++ b/R/CST_MultivarRMSE.R
@@ -3,8 +3,8 @@
 #'@author Deborah Verfaillie, \email{deborah.verfaillie@bsc.es}
 #'@description This function calculates the RMSE from multiple variables, as the mean of each variable's RMSE scaled by its observed standard deviation. Variables can be weighted based on their relative importance (defined by the user).
 #'
-#'@param exp a list of objects, one for each variable, of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.
-#'@param obs a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.
+#'@param exp a list of objects, one for each variable, 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 obs a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the observed anomaly data in the element named \code{$data}.
 #'
 #'@param weight (optional) a vector of weight values to assign to each variable. If no weights are defined, a value of 1 is assigned to every variable.
 #'
@@ -15,7 +15,8 @@
 #'@examples
 #'# Creation of sample s2dverification objects. These are not complete
 #'# s2dverification objects though. The Load function returns complete objects.
-#'
+#'# using package zeallot is optional:
+#' library(zeallot)
 #'# Example with 2 variables
 #'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
 #'mod2 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
@@ -25,25 +26,27 @@
 #'obs2 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
 #'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
 #'dim(obs2) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
-#'lon <- seq(0, 25, 5)
-#'lat <- seq(0, 30, 5)
+#'lon <- seq(0, 30, 5)
+#'lat <- seq(0, 25, 5)
 #'exp1 <- list(data = mod1, lat = lat, lon = lon, Datasets = "EXP1", 
 #'             source_files = "file1", Variable = list('pre'))
 #'attr(exp1, 'class') <- 's2dv_cube'
 #'exp2 <- list(data = mod2, lat = lat, lon = lon, Datasets = "EXP2", 
 #'             source_files = "file2", Variable = list('tas'))
 #'attr(exp2, 'class') <- 's2dv_cube'
-#'exp <- list(exp1, exp2)
 #'obs1 <- list(data = obs1, lat = lat, lon = lon, Datasets = "OBS1", 
 #'             source_files = "file1", Variable = list('pre'))
 #'attr(obs1, 'class') <- 's2dv_cube'
 #'obs2 <- list(data = obs2, lat = lat, lon = lon, Datasets = "OBS2", 
 #'             source_files = "file2", Variable = list('tas'))
 #'attr(obs2, 'class') <- 's2dv_cube'
-
-#'obs <- list(obs1, obs2)
+#'
+#'c(ano_exp1, ano_obs1) %<-% CST_Anomaly(exp1, obs1, cross = TRUE, memb = TRUE)
+#'c(ano_exp2, ano_obs2) %<-% CST_Anomaly(exp2, obs2, cross = TRUE, memb = TRUE)
+#'ano_exp <- list(exp1, exp2)
+#'ano_obs <- list(ano_obs1, ano_obs2)
 #'weight <- c(1, 2)
-#'a <- CST_MultivarRMSE(exp = exp, obs = obs, weight = weight)
+#'a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight)
 #'str(a)
 #'@export
 CST_MultivarRMSE <- function(exp, obs, weight = NULL) {
@@ -104,10 +107,9 @@ CST_MultivarRMSE <- function(exp, obs, weight = NULL) {
   mvrmse <- 0
   sumweights <- 0
   for (j in 1 : nvar) {
-    ano <- Ano_CrossValid(exp[[j]]$data, obs[[j]]$data)
     # seasonal average of anomalies
-    AvgExp <- MeanListDim(ano$ano_exp, narm = TRUE, c(2, 4))
-    AvgObs <- MeanListDim(ano$ano_obs, narm = TRUE, c(2, 4))
+    AvgExp <- MeanListDim(exp[[j]]$data, narm = TRUE, c(2, 4))
+    AvgObs <- MeanListDim(obs[[j]]$data, narm = TRUE, c(2, 4))
     # multivariate RMSE (weighted) 
     rmse <- s2dverification::RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, 
                                  conf = FALSE)
diff --git a/man/CST_Anomaly.Rd b/man/CST_Anomaly.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..b214a31a299c6736e7704bec16378ba86ef7b0d4
--- /dev/null
+++ b/man/CST_Anomaly.Rd
@@ -0,0 +1,64 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_Anomaly.R
+\name{CST_Anomaly}
+\alias{CST_Anomaly}
+\title{Anomalies relative to a climatology along selected dimension with or without cross-validation}
+\usage{
+CST_Anomaly(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
+  dim_anom = 3)
+}
+\arguments{
+\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.}
+
+\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.'}
+
+\item{cross}{A logical value indicating whether cross-validation should be applied or not. Default = FALSE.}
+
+\item{memb}{A logical value indicating whether Clim() computes one climatology for each experimental data 
+product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE.}
+
+\item{dim_anom}{An integer indicating the dimension along which the climatology will be computed. It 
+usually corresponds to 3 (sdates) or 4 (ftime). Default = 3.}
+}
+\value{
+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.
+}
+\description{
+This function computes the anomalies relative to a climatology computed along the 
+selected dimension (usually starting dates or forecast time) allowing the application or not of
+crossvalidated climatologies. The computation is carried out independently for experimental and 
+observational data products.
+}
+\examples{
+# Example 1:
+mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
+dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
+obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
+dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
+lon <- seq(0, 30, 5)
+lat <- seq(0, 25, 5)
+exp <- list(data = mod, lat = lat, lon = lon)
+obs <- list(data = obs, lat = lat, lon = lon)
+attr(exp, 'class') <- 's2dv_cube'
+attr(obs, 'class') <- 's2dv_cube'
+
+anom1 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE)
+str(anom1)
+anom2 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
+str(anom2)
+
+anom3 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = FALSE)
+str(anom3)
+
+anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE)
+str(anom4)
+
+}
+\seealso{
+\code{\link[s2dverification]{Ano_CrossValid}}, \code{\link[s2dverification]{Clim}} and \code{\link{CST_Load}}
+}
+\author{
+Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+
+Pena Jesus, \email{jesus.pena@bsc.es}
+}
diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd
index 0f5f2acff9c6fd9a51bb121be18e99ba8c011f73..99fdb400f1a82dd7267cd284a2b6597dfb314797 100644
--- a/man/CST_MultiMetric.Rd
+++ b/man/CST_MultiMetric.Rd
@@ -7,9 +7,9 @@
 CST_MultiMetric(exp, obs, metric = "correlation", multimodel = TRUE)
 }
 \arguments{
-\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing 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 experiment data in the element named \code{$data}.}
 
-\item{obs}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed 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.}
 
@@ -22,6 +22,7 @@ an object of class \code{s2dv_cube} containing the statistics of the selected me
 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.
 }
 \examples{
+library(zeallot)
 mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
 dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
 obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
@@ -32,7 +33,8 @@ exp <- list(data = mod, lat = lat, lon = lon)
 obs <- list(data = obs, lat = lat, lon = lon)
 attr(exp, 'class') <- 's2dv_cube'
 attr(obs, 'class') <- 's2dv_cube'
-a <- CST_MultiMetric(exp = exp, obs = obs)
+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)
 }
 \references{
diff --git a/man/CST_MultivarRMSE.Rd b/man/CST_MultivarRMSE.Rd
index 5db63421519e9e0e2d89a915d0f5ac942921c4d4..24af608c8360985de73763434a16f1a11fd35ffe 100644
--- a/man/CST_MultivarRMSE.Rd
+++ b/man/CST_MultivarRMSE.Rd
@@ -7,9 +7,9 @@
 CST_MultivarRMSE(exp, obs, weight = NULL)
 }
 \arguments{
-\item{exp}{a list of objects, one for each variable, of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.}
+\item{exp}{a list of objects, one for each variable, 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{obs}{a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the observed data in the element named \code{$data}.}
+\item{obs}{a list of objects, one for each variable (in the same order than the input in 'exp') of class \code{s2dv_cube} as returned by \code{CST_Anomaly} function, containing the observed anomaly data in the element named \code{$data}.}
 
 \item{weight}{(optional) a vector of weight values to assign to each variable. If no weights are defined, a value of 1 is assigned to every variable.}
 }
@@ -22,7 +22,8 @@ This function calculates the RMSE from multiple variables, as the mean of each v
 \examples{
 # Creation of sample s2dverification objects. These are not complete
 # s2dverification objects though. The Load function returns complete objects.
-
+# using package zeallot is optional:
+library(zeallot)
 # Example with 2 variables
 mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
 mod2 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
@@ -32,24 +33,27 @@ obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
 obs2 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
 dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
 dim(obs2) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
-lon <- seq(0, 25, 5)
-lat <- seq(0, 30, 5)
+lon <- seq(0, 30, 5)
+lat <- seq(0, 25, 5)
 exp1 <- list(data = mod1, lat = lat, lon = lon, Datasets = "EXP1", 
             source_files = "file1", Variable = list('pre'))
 attr(exp1, 'class') <- 's2dv_cube'
 exp2 <- list(data = mod2, lat = lat, lon = lon, Datasets = "EXP2", 
             source_files = "file2", Variable = list('tas'))
 attr(exp2, 'class') <- 's2dv_cube'
-exp <- list(exp1, exp2)
 obs1 <- list(data = obs1, lat = lat, lon = lon, Datasets = "OBS1", 
             source_files = "file1", Variable = list('pre'))
 attr(obs1, 'class') <- 's2dv_cube'
 obs2 <- list(data = obs2, lat = lat, lon = lon, Datasets = "OBS2", 
             source_files = "file2", Variable = list('tas'))
 attr(obs2, 'class') <- 's2dv_cube'
-obs <- list(obs1, obs2)
+
+c(ano_exp1, ano_obs1) \%<-\% CST_Anomaly(exp1, obs1, cross = TRUE, memb = TRUE)
+c(ano_exp2, ano_obs2) \%<-\% CST_Anomaly(exp2, obs2, cross = TRUE, memb = TRUE)
+ano_exp <- list(exp1, exp2)
+ano_obs <- list(ano_obs1, ano_obs2)
 weight <- c(1, 2)
-a <- CST_MultivarRMSE(exp = exp, obs = obs, weight = weight)
+a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight)
 str(a)
 }
 \seealso{
diff --git a/vignettes/MultiModelSkill_vignette.md b/vignettes/MultiModelSkill_vignette.md
index 6bebe138672b60eeb212ef24daac0c99c8eaefac..6fc676109908add0004f31bdcfe2ab141b33cd4a 100644
--- a/vignettes/MultiModelSkill_vignette.md
+++ b/vignettes/MultiModelSkill_vignette.md
@@ -71,16 +71,19 @@ Using the `CST_Load` function, the data available in our data store can be loade
 Ask nuria.perez at bsc.es to achieve the data to run the recipe.
 
 ```r  
-#glosea5 <- '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc'
+glosea5 <- '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc'
 c(exp, obs) %<-% 
-   CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5), list(name = 'ecmwf/system4_m1'),
+   CST_Load(var = clim_var, exp = list(list(name = 'glosea5', path = glosea5), 
+                                  list(name = 'ecmwf/system4_m1'),
                                   list(name = 'meteofrance/system5_m1')),
                   obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
                   lonmin = -20, lonmax = 70, latmin = 25, latmax = 75, 
                   storefreq = "monthly", sampleperiod = 1,  nmember = 9, 
                   output = "lonlat", method = "bilinear", 
                   grid = paste("r", grid, sep = ""))
-#save(exp, obs, file = "../tas_toydata.RData")         
+#save(exp, obs, file = "../tas_toydata.RData")   
+
+# Or use the following line to load the file provided in .RData format:
 load(file = "./tas_toydata.RData")
 ```
 
@@ -108,11 +111,31 @@ Lon <- exp$lon
 
 ### 2.- Computing and plotting Anomaly Correlation Coefficient
 
-The Anomaly Correlation Coefficient (ACC) is the most widely used skill metric for Seasonal Climate Forecast quality (Mishra et al., 2018). The ACC is obtained by running the `CST_MultiMetric` function defining the parameter 'metric' as correlation. Furthermore, the function includes the option of computing the Multi-Model Mean ensemble (MMM).
+The Anomaly Correlation Coefficient (ACC) is the most widely used skill metric for Seasonal Climate Forecast quality (Mishra et al., 2018). 
+
+
+First step is to compute the anomalies over the loaded data applying cross validation technique on individual members by running:
+
+```
+c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
+```
+
+The dimensions are preserved:
+
+```
+> str(ano_exp$data)
+ num [1:3, 1:9, 1:21, 1:3, 1:35, 1:64] -1.647 -0.478 -0.096 1.575 1.086 ...
+ - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ...
+> str(ano_obs$data)
+ num [1, 1, 1:21, 1:3, 1:35, 1:64] 0.0235 1.546 1.3885 -0.344 -5.972 ...
+ - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ...
+```
+
+The ACC is obtained by running the `CST_MultiMetric` function defining the parameter 'metric' as correlation. The function also includes the option of computing the Multi-Model Mean ensemble (MMM).
 
 
 ```r
-AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'correlation', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'correlation', 
                           multimodel = TRUE)
 ```
 
@@ -168,7 +191,7 @@ The next figure is the map of the maximum positive Anomaly Correlation Coefficie
 The same function can be used to compute the RMS error by defining the parameter `metric` as 'rms'.
 
 ```r
-AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rms', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rms', 
                           multimodel = TRUE)
 RMS <- AnomDJF$data[ , , 2, , ]
 ```
@@ -201,7 +224,7 @@ Notice that the perfect RMSSS is 1 and the parameter  `map_select_fun` from `Plo
 
 
 ```r
-AnomDJF <- CST_MultiMetric(exp = exp, obs = obs, metric = 'rmsss', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', 
                           multimodel = TRUE)
 RMSSS <- AnomDJF$data[ , , 1, , ]
 names(dim(RMSSS)) <- c("maps", "lat", "lon")
diff --git a/vignettes/MultivarRMSE_vignette.md b/vignettes/MultivarRMSE_vignette.md
index fc18b99fb9185282e6c88043d7e93f8b707f774c..78af60807aa75815237fed804619e25f172de392 100644
--- a/vignettes/MultivarRMSE_vignette.md
+++ b/vignettes/MultivarRMSE_vignette.md
@@ -17,7 +17,6 @@ To run this vignette, the next R packages should be installed and loaded:
 
 ```r
 library(s2dverification)
-library(startR)
 library(RColorBrewer)
 library(zeallot)
 ```
@@ -70,14 +69,13 @@ Using the `CST_Load` function from **CSTool package**, the data available in our
 Ask nuria.perez at bsc.es for the data to run the recipe.
 
 ```r  
-#glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc')
+glosea5 <- list(path = '/esnas/exp/glosea5/specs-seasonal_i1p1/$STORE_FREQ$_mean/$VAR_NAME$-allmemb/$VAR_NAME$_$START_DATE$.nc')
  c(exp_T, obs_T) %<-% 
    CST_Load(var = temp, exp = list(glosea5),
      obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
      latmin = 25, latmax = 75, lonmin = -20, lonmax = 70, output = 'lonlat',
      nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9,
      method = "bilinear", grid = paste("r", grid, sep = ""))
-#ecmwf <- list(path = '/esnas/exp/ecmwf/system4_m1/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc')
 c(exp_P, obs_P) %<-% 
    CST_Load(var = precip,  exp = list(glosea5),
      obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
@@ -85,7 +83,9 @@ c(exp_P, obs_P) %<-%
      nprocs = 1, storefreq = "monthly", sampleperiod = 1, nmember = 9,
      method = "bilinear", grid = paste("r", grid, sep = ""))
 #save(exp_T, obs_T, exp_P, obs_P, file = "./tas_prlr_toydata.RData")         
-load(file = "../tas_prlr_toydata.RData")
+
+# Or use the following line to load the file provided in .RData format:
+load(file = "./tas_prlr_toydata.RData")
 ```
 
 There should be four new elements loaded in the R working environment: `exp_T`, `obs_T`, `exp_P` and `obs_P`. The first two elements correspond to the experimental and observed data for temperature and the other are the equivalent for the precipitation data. It's possible to check that they are of class `sd2v_cube` by running:
@@ -118,14 +118,32 @@ Lat <- exp_T$lat
 Lon <- exp_T$lon
 ```
 
+The next step is to compute the anomalies of the experimental and observational data using `CST_Anomaly` function, which could be applied over data from each variable, and in this case it's compute applying cross validation technique over individual members:
+
+```
+c(ano_exp_T, ano_obs_T) %<-% CST_Anomaly(exp = exp_T, obs = obs_T, cross = TRUE, memb = TRUE)
+c(ano_exp_P, ano_obs_P) %<-% CST_Anomaly(exp = exp_P, obs = obs_P, cross = TRUE, memb = TRUE)
+```
+
+The original dimensions are preserved and the anomalies are stored in the `data` element of the correspondent object:
+
+```
+> str(ano_exp_T$data)
+ num [1, 1:9, 1:21, 1:3, 1:35, 1:64] -1.647 1.575 2.77 0.048 -1.886 ...
+ - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ...
+> str(ano_obs_T$data)
+ num [1, 1, 1:21, 1:3, 1:35, 1:64] 0.0235 1.546 1.3885 -0.344 -5.972 ...
+ - attr(*, "dimensions")= chr [1:6] "dataset" "member" "sdate" "ftime" ...
+```
+
+Two lists containing the experiment ,`ano_exp`, and the observation, `ano_obs`, lists should be put together to serve as input of the function to compute multivariate RMSEs.
 
-Two lists containing the `exp` and `obs` lists should be put together to serve as input of the function to compute multivariate RMSEs.
 Furthermore, some weights can be applied to the difference variables based on their relative importance (if no weights are given, a value of 1 is automatically assigned to each variable). For this example, we'll give a weight of 2 to the temperature dataset and a weight of 1 to the precipitation dataset:
 
 
 ```r
-exp <- list(exp_T, exp_P)
-obs <- list(obs_T, obs_P)
+ano_exp <- list(ano_exp_T, ano_exp_P)
+ano_obs <- list(ano_obs_T, ano_obs_P)
 weight <- c(2, 1)
 ```
 
@@ -136,7 +154,7 @@ It is obtained by running the `CST_MultivarRMSE` function:
 
 
 ```r
-mvrmse <- CST_MultivarRMSE(exp = exp, obs = obs, weight)
+mvrmse <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight)
 ```