diff --git a/.Rbuildignore b/.Rbuildignore
index 43690b7004542a97fd35b30d4f2dd4bb6e43bf62..62a63bb483658fb6c751ec2ebd39bb69a8d15f2c 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -3,5 +3,5 @@
 .*\.tar.gz$
 .*\.pdf$
 ./.nc$
-.*\.RData$
+.*^(?!data)\.RData$
 .*\.gitlab-ci.yml$
diff --git a/.gitignore b/.gitignore
index 7937bb7b02565bf89478ef96048dda30674562fd..2f6c062a8dacc06b63d3a1d570bf3bc74adb5cd3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,5 +15,5 @@ master_pull.txt
 *.ps
 Rplots.pdf
 .nfs*
-.RData
-data/
+*.RData
+!data/*.RData
diff --git a/DESCRIPTION b/DESCRIPTION
index efdbbd988fa1c9725f566b890e39d03c2851f0cb..6b64acde7b126e43fcf2358e2fcabe09cdb5f895 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,6 @@
 Package: CSTools
-Title: Set of Tools to Assess Skills Climate Forecasts on Seasonal-to-Decadal
-    Timescales
-Version: 0.0.1
+Title: Set of Tools to Assess Skill of Climate Forecasts on Seasonal-to-Decadal Timescales
+Version: 1.0.0
 Authors@R: c(
     person("BSC-CNS", role = c("cph")),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")),
@@ -13,13 +12,9 @@ Authors@R: c(
     person("Veronica", "Torralba", , "veronica.torralba@bsc.es", role = "aut"),
     person("Deborah", "Verfaillie", , "deborah.verfaillie@bsc.es", role = "aut"),
     person("Lauriane", "Batte", , "lauriane.batte@meteo.fr", role = "ctb"),
+    person("Jesus", "Peña", , "jesus.pena@bsc.es", role = "ctb"),
     person("Bert", "van Schaeybroeck", , "bertvs@meteo.be", role = "ctb"))
-Description: This package provides tools to exploit dynamical seasonal forecasts
-    in order to provide information relevant to stakeholders at the seasonal 
-    timescale. This toolbox contains process-based methods for forecast 
-    calibration, bias correction, statistical and stochastic downscaling, 
-    optimal forecast combination and multivariate verification, as well as 
-    basic and advanced tools to obtain tailored products.
+Description: Tools to exploit dynamical seasonal forecasts in order to provide information relevant to stakeholders at the seasonal timescale. The package contains process-based methods for forecast calibration, bias correction, statistical and stochastic downscaling, optimal forecast combination and multivariate verification, as well as basic and advanced tools to obtain tailored products.
 Depends:
     R (>= 3.2.0),
     maps
@@ -27,6 +22,9 @@ Imports:
     s2dverification,
     abind,
     stats,
+    utils,
+    grDevices,
+    graphics,
     JuliaCall,
     ncdf4,
     ggplot2,
@@ -35,7 +33,10 @@ Imports:
     reshape2,
     multiApply
 Suggests:
-    testthat
+    zeallot,
+    testthat,
+    R.rsp
+VignetteBuilder: R.rsp
 License: Apache License 2.0 | file LICENSE
 Encoding: UTF-8
 LazyData: true
diff --git a/NAMESPACE b/NAMESPACE
index dad92f85818b81a230118c79df6a4233c96702c6..872506f75a52612871a1f7a21a1cd06876fec13f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,9 @@
 # Generated by roxygen2: do not edit by hand
 
+export(CST_Anomaly)
 export(CST_BiasCorrection)
 export(CST_Calibration)
+export(CST_Load)
 export(CST_MultiMetric)
 export(CST_MultivarRMSE)
 export(CST_RFSlope)
@@ -20,7 +22,29 @@ import(stats)
 importFrom(data.table,CJ)
 importFrom(data.table,data.table)
 importFrom(data.table,setkey)
+importFrom(grDevices,adjustcolor)
+importFrom(grDevices,bmp)
+importFrom(grDevices,colorRampPalette)
+importFrom(grDevices,dev.cur)
+importFrom(grDevices,dev.new)
+importFrom(grDevices,dev.off)
+importFrom(grDevices,hcl)
+importFrom(grDevices,jpeg)
+importFrom(grDevices,pdf)
+importFrom(grDevices,png)
+importFrom(grDevices,postscript)
+importFrom(grDevices,svg)
+importFrom(grDevices,tiff)
+importFrom(graphics,box)
+importFrom(graphics,image)
+importFrom(graphics,layout)
+importFrom(graphics,mtext)
+importFrom(graphics,par)
+importFrom(graphics,plot.new)
 importFrom(maps,map)
 importFrom(plyr,.)
 importFrom(plyr,dlply)
 importFrom(reshape2,melt)
+importFrom(s2dverification,Load)
+importFrom(utils,glob2rx)
+importFrom(utils,installed.packages)
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_BiasCorrection.R b/R/CST_BiasCorrection.R
index 5fc35696a547f54ad538dfe697f3619c1ccb143f..c186936b705b7aa3a4db1350fe155f19530718ea 100644
--- a/R/CST_BiasCorrection.R
+++ b/R/CST_BiasCorrection.R
@@ -3,9 +3,10 @@
 #'@author Verónica Torralba, \email{veronica.torralba@bsc.es}
 #'@description This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset.
 #'
-#'@param data CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).
+#'@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}.
 #'
-#'@return a CSTools object (s2dverification object) with the bias corrected forecasts a element called \code{data$biascalibration}.
+#'@return an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data.
 #'
 #'@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)
 #'
@@ -22,45 +23,32 @@
 #'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
 #'lon <- seq(0, 30, 5)
 #'lat <- seq(0, 25, 5)
-#'data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-#'a <- CST_BiasCorrection(data1)
+#'exp <- list(data = mod1, lat = lat, lon = lon)
+#'obs <- list(data = obs1, lat = lat, lon = lon)
+#'attr(exp, 'class') <- 's2dv_cube'
+#'attr(obs, 'class') <- 's2dv_cube'
+#'a <- CST_BiasCorrection(exp = exp, obs = obs)
 #'str(a)
 #'@export
-CST_BiasCorrection <- function(data) {
-  is_s2dv_object <- function(x) {
-    if (all(c('mod', 'obs', 'lon', 'lat') %in% names(x))) { #&& length(x) == 11) {
-      TRUE
-    } else {
-      FALSE
-    }
-  }
-  wrong_input <- FALSE
-  if (!is.list(data)) {
-    wrong_input <- TRUE
-  }
-  
-  if (!is_s2dv_object(data)) {
-    wrong_input <- TRUE
-  }
-
-if (wrong_input) {
-    stop("Parameter 'data' must be a list as output of Load function
-    from s2dverification package")
+CST_BiasCorrection <- function(exp, obs) {
+  if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) {
+    stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
+         "as output by CSTools::CST_Load.")
   }
 
-  if (dim(data$obs)['member']!=1) {
-    stop("The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.")
+  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.")
   }
 
-  BiasCorrected <- BiasCorrection(data$mod, data$obs)
+  BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data)
   dimnames <- names(dim(BiasCorrected))
   BiasCorrected <- aperm(BiasCorrected, c(3, 1, 2, 4, 5, 6))
   names(dim(BiasCorrected)) <- dimnames[c(3, 1, 2, 4, 5, 6)]
-  data$biascorrection <- BiasCorrected
-  data$obs <- NULL
-  data$mod <- NULL
-  data <- data[c("biascorrection", names(data)[-which(names(data) == "biascorrection")])]
-  return(data)
+  exp$data <- BiasCorrected
+  exp$Datasets <- c(exp$Datasets, obs$Datasets)
+  exp$source_files <- c(exp$source_files, obs$source_files)
+  return(exp)
 }
 
 BiasCorrection <- function (exp, obs) {
diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R
index 238dfbd9e15e9fdbf381a6b25741a7a4f2608b2f..a9b9f78fb12d62eaf8f7cf72bf5d92bdd620400f 100644
--- a/R/CST_Calibration.R
+++ b/R/CST_Calibration.R
@@ -5,62 +5,45 @@
 #'
 #'@references Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting—II calibration and combination. Tellus A. 2005;57:234–252. doi:10.1111/j.1600-0870.2005.00104.x
 #'
-#'@param data a CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).
+#'@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}.
 #'
-#'@return a CSTools object (s2dverification object) with the calibrated forecasts in a element called \code{data$calibration}.
+#'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data.
 #'
 #'@import s2dverification
 #'@import multiApply
 #'
+#'@seealso \code{\link{CST_Load}}
+#'
 #'@examples
+#'\donttest{
 #'# Example
-#'# Creation of sample s2dverification objects. These are not complete
-#'# s2dverification objects though. The Load function returns complete objects.
-#'mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
-#'dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
-#'obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
-#'dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
-#'lon <- seq(0, 30, 5)
-#'lat <- seq(0, 25, 5)
-#'data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-#'a <- CST_Calibration(data1)
-#'str(a)
-#'
+#'# Load data using CST_Load or use the sample data provided:
+#'library(zeallot)
+#'c(exp, obs) %<-% lonlat_data
+#'exp_calibrated <- CST_Calibration(exp = exp, obs = obs)
+#'str(exp_calibrated)
+#'}
 #'@export
-CST_Calibration <- function(data) {
-  is_s2dv_object <- function(x) {
-    if (all(c('mod', 'obs', 'lon', 'lat') %in% names(x))) { #&& length(x) == 11) {
-      TRUE
-    } else {
-      FALSE
-    }
-  }
-  wrong_input <- FALSE
-  if (!is.list(data)) {
-    wrong_input <- TRUE
-  }
-  
-  if (!is_s2dv_object(data)) {
-    wrong_input <- TRUE
-  }
-  if (wrong_input) {
-    stop("Parameter 'data' must be a list as output of Load function
-from s2dverification package.")
+CST_Calibration <- function(exp, obs) {
+  if (!inherits(exp, 's2dv_cube') || !inherits(exp, 's2dv_cube')) {
+    stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
+         "as output by CSTools::CST_Load.")
   }
 
-  if (dim(data$obs)[which(names(dim(data$obs)) == 'member')] != 1) {
-    stop("The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.")
+  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.")
   }
-  
-  Calibrated <- Calibration(data$mod, data$obs)
+
+  Calibrated <- Calibration(exp = exp$data, obs = obs$data)
   dimnames <- names(dim(Calibrated))
   Calibrated <- aperm(Calibrated, c(3, 1, 2, 4, 5, 6))
   names(dim(Calibrated)) <- dimnames[c(3, 1, 2, 4, 5, 6)] 
-  data$calibration <- Calibrated
-  data$obs <- NULL
-  data$mod <- NULL
-  data <- data[c("calibration", names(data)[-which(names(data) == "calibration")])]
-  return(data)
+  exp$data <- Calibrated
+  exp$Datasets <- c(exp$Datasets, obs$Datasets)
+  exp$source_files <- c(exp$source_files, obs$source_files)
+  return(exp)
 }
 
 Calibration <- function(exp, obs) {
diff --git a/R/CST_Load.R b/R/CST_Load.R
new file mode 100644
index 0000000000000000000000000000000000000000..5e1dc3aaa5f20321e0ce148fd2f41c527438ca2a
--- /dev/null
+++ b/R/CST_Load.R
@@ -0,0 +1,111 @@
+#' CSTools Data Retreival Function
+#'
+#' This function aggregates, subsets and retrieves sub-seasonal, seasonal, decadal or climate projection data from NetCDF files in a local file system or on remote OPeNDAP servers, and arranges it for easy application of the CSTools functions.
+#'
+#' It receives any number of parameters (`...`) that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`.
+#'
+#' It is recommended to use this function in combination with the `zeallot::"%<-%"` operator, to directly assing the two returned 's2dv_cube's to two separate variables, which can then be sent independently to other functions in CSTools as needed. E.g.: `c(exp, obs) <- CST_Load(...)`.
+#'
+#' @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}
+#' @importFrom s2dverification Load
+#' @importFrom utils glob2rx
+#' @export
+#' @examples
+#' \dontrun{
+#' library(zeallot)
+#' startDates <- c('20001101', '20011101', '20021101', 
+#'                 '20031101', '20041101', '20051101')
+#' c(exp, obs) %<-% 
+#'   CST_Load(
+#'     var = 'tas', 
+#'     exp = 'system5c3s', 
+#'     obs = 'era5', 
+#'     nmember = 15,
+#'     sdates = startDates,
+#'     leadtimemax = 3,
+#'     latmin = 27, latmax = 48,
+#'     lonmin = -12, lonmax = 40,
+#'     output = 'lonlat',
+#'     nprocs = 1
+#'   )
+#' }
+#' \dontshow{
+#' exp <- CSTools::lonlat_data$exp
+#' obs <- CSTools::lonlat_data$obs
+#' }
+CST_Load <- function(...) {
+  exp <- s2dverification::Load(...)
+
+  if (is.null(exp) || (is.null(exp$mod) && is.null(exp$obs))) {
+    stop("The s2dverification::Load call did not return any data.")
+  }
+
+  obs <- exp
+  obs$mod <- NULL
+  exp$obs <- NULL
+  names(exp)[[1]] <- 'data'
+  names(obs)[[1]] <- 'data'
+
+  remove_matches <- function(v, patterns) {
+    if (length(v) > 0) {
+      matches <- c()
+      for (pattern in patterns) {
+        matches <- c(matches, which(grepl(pattern, v)))
+      }
+      if (length(matches) > 0) {
+        v <- v[-matches]
+      }
+    }
+    v
+  }
+
+  harmonize_patterns <- function(v) {
+    matches <- grepl('.*\\.nc$', v)
+    if (sum(!matches) > 0) {
+      match_indices <- which(!matches)
+      v[match_indices] <- sapply(v[match_indices], function(x) paste0(x, '*'))
+    }
+    v <- glob2rx(v)
+    v <- gsub('\\$.*\\$', '*', v)
+    v
+  }
+
+  if (!is.null(obs$data)) {
+    obs$Datasets$exp <- NULL
+    obs$Datasets <- obs$Datasets$obs
+    obs_path_patterns <- sapply(obs$Datasets, function(x) attr(x, 'source'))
+    obs_path_patterns <- harmonize_patterns(obs_path_patterns)
+  }
+
+  if (!is.null(exp$data)) {
+    exp$Datasets$obs <- NULL
+    exp$Datasets <- exp$Datasets$exp
+    exp_path_patterns <- sapply(exp$Datasets, function(x) attr(x, 'source'))
+    exp_path_patterns <- harmonize_patterns(exp_path_patterns)
+  }
+
+  if (!is.null(obs$data) && !is.null(exp$data)) {
+    obs$source_files <- remove_matches(obs$source_files,
+                                       exp_path_patterns)
+    obs$not_found_files <- remove_matches(obs$not_found_files,
+                                          exp_path_patterns)
+
+    exp$source_files <- remove_matches(exp$source_files,
+                                       obs_path_patterns)
+    exp$not_found_files <- remove_matches(exp$not_found_files,
+                                          obs_path_patterns)
+  }
+
+  result <- list()
+  if (!is.null(exp$data)) {
+    class(exp) <- 's2dv_cube'
+    result$exp <- exp
+  }
+  if (!is.null(obs$data)) {
+    class(obs) <- 's2dv_cube'
+    result$obs <- obs
+  }
+  result
+}
diff --git a/R/CST_MultiMetric.R b/R/CST_MultiMetric.R
index 1939622121a90f194ebd4255cd74e8fb11460ac1..e7bfc729f7b8a87092254fe2bbd32ac8e981e525 100644
--- a/R/CST_MultiMetric.R
+++ b/R/CST_MultiMetric.R
@@ -4,84 +4,74 @@
 #'@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 data a CSTools object (list) giving as output of \code{Load} function from S2dverification package.
+#'@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.
-#'@param names a character vector indicating the names of the models or experiments to be compared.
 #'
-#'@return A CSTools object as the input parameter \code{data}, without the elements \code{mod} and \code{obs} and including:
-#'\itemize{
-#'\item\code{$ano_exp} {an array with the same dimensions as the \code{mod} label in the input \code{data}.}
-#'\item\code{$ano_exp} {an array with the same dimensions as the \code{obs} label in the input \code{data}.}
-#'\item\code{$metric} {An array with two datset dimensions from the \code{mod} and \code{obs} label in the input \code{data}. 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.}}
-#'@details The output will include attributes in element \code{ano_exp} returning a vector with the names of the models or experiments.
-#'@seealso \code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{RMSSS}} and \code{\link[s2dverification]{Load}}
+#'@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}}
 #'@references 
-#' \url{http://link.springer.com/10.1007/s00382-018-4404-z}
+#'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
 #'@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)
 #'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)
-#'dat <- list(mod = mod, obs = obs, lat = lat, lon = lon)
-#'a <- CST_MultiMetric(dat)
+#'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'
+#'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(data, metric = "correlation", multimodel = TRUE, names = NULL) {
-  if (!is.list(data)) {
-    stop("Parameter 'data' must be a list as output of Load function
-from s2dverification package.")
+CST_MultiMetric <- function(exp, obs, metric = "correlation", multimodel = TRUE) {
+  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 (!(all(c('mod', 'obs', 'lon', 'lat') %in% names(data)) && length(data) >= 4)) {
-      stop("Parameter 'data' must contain at least 4 elements 'mod', 'obs', 
-           'lon', 'lat'.")
+
+  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(data$mod))) & !is.null(names(dim(data$obs)))) {
-    if (all(names(dim(data$mod)) %in% names(dim(data$obs)))) {
-      dimnames <- names(dim(data$mod))
+
+  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 elements 'mod' and 'obs' from parameter 'data' 
-           should have the same name dimmension.")
+      stop("Dimension names of element 'data' from parameters 'exp'",
+           " and 'obs' should have the same name dimmension.")
     }
   } else {
-    stop("Elements 'mod' and 'obs' from parameter 'data' should have 
-         dimmension names.")
+    stop("Element 'data' from parameters 'exp' and 'obs'",
+         " should have dimmension names.")
   }
   
-  if (is.null(names)) {
-    if (!is.null(names(data$Datasets$exp))) {
-      names <-  names(data$Datasets$exp)
-    } else {
-      names <- paste0('mod', 1 : dim(data$mod)[1])
-    }
-  }
-  if (!is.null(names)) {
-    if (length(names) != dim(data$mod)[1]) {
-      stop("Parameter 'names' must be the same length as the number of 
-           models in the element 'mod' of parameter 'data'.")
-    }
-  }
   if (!is.logical(multimodel)) {
     stop("Parameter 'multimodel' must be a logical value.")
   }
   if (length(multimodel) > 1) {
     multimodel <- multimodel[1]
-    warning("Parameter 'multimodel' has length > 1 and only the first element will be used.")
+    warning("Parameter 'multimodel' has length > 1 and only the first ",
+            "element will be used.")
   }
   if (length(metric) > 1) {
     metric <- metric[1]
-    warning("Parameter 'multimodel' has length > 1 and only the first element will be used.")
+    warning("Parameter 'multimodel' has length > 1 and only the first ",
+            "element will be used.")
   }
-  ano <- Ano_CrossValid(var_exp = data$mod, var_obs = data$obs)
   
   # 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') {
@@ -91,8 +81,8 @@ from s2dverification package.")
   } else if (metric == 'rmsss') {
     corr <- RMSSS(AvgExp, AvgObs, posloop = 1, posRMS = 2)
   } else {
-    stop("Parameter 'metric' must be a character string indicating one of the 
-         options: 'correlation', 'rms' or 'rmse'.")
+    stop("Parameter 'metric' must be a character string indicating ",
+         "one of the options: 'correlation', 'rms' or 'rmse'.")
   }
   if (multimodel == TRUE) {
     # seasonal avg of anomalies for multi-model
@@ -112,20 +102,15 @@ from s2dverification package.")
                         var_obs = InsertDim(AvgObs_MMM, 1, 1), 
                         posloop = 1, posRMS = 2)
     } 
-       corr <- abind::abind(corr, corr_MMM, along = 1)
-    names <- c(names, 'MMM')
+    corr <- abind::abind(corr, corr_MMM, along = 1)
   }
   
   names(dim(corr)) <- c(dimnames[1], dimnames[1], 'statistics', dimnames[5 : 6])
   
-  data$ano_exp <- ano$ano_exp
-  data$ano_obs <- ano$ano_obs
-  data$metric <- corr
-  data$obs <- NULL
-  data$mod <- NULL
-#  print(names)
- # attributes(data$exp)$names <- NULL
-  data <- data[order(names(data))]
-  attr(data[[1]], "experiments") <- names
-  return(data)
+  #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)
 }
diff --git a/R/CST_MultivarRMSE.R b/R/CST_MultivarRMSE.R
index 05138b1071bb29e06c4a14241f5ad21b142630e4..88967a134aca76502e3eb3bf1a32d18b9bf48217 100644
--- a/R/CST_MultivarRMSE.R
+++ b/R/CST_MultivarRMSE.R
@@ -3,17 +3,20 @@
 #'@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 data a list of s2dverification objects (lists) as output by the \code{Load} function from the s2dverification package, one for each variable.
+#'@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.
 #'
-#'@return \code{$mvrmse} {An array with dimensions: c(number of exp, number of obs, 1 (the multivariate RMSE value), number of lat, number of lon)}
+#'@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
 #'@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)
@@ -23,69 +26,113 @@
 #'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)
-#'data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-#'data2 <- list(mod = mod2, obs = obs2, lat = lat, lon = lon)
-#'data <- list(data1, data2)
+#'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'
+#'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'
+#'
+#'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(data, weight)
+#'a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight)
 #'str(a)
 #'@export
-CST_MultivarRMSE <- function(data, weight = NULL) {
-  is_s2dv_object <- function(x) {
-    if (all(c('mod', 'obs', 'lon', 'lat') %in% names(x))) { #&& length(x) == 11) {
-      TRUE
-    } else {
-      FALSE
-    }
+CST_MultivarRMSE <- function(exp, obs, weight = NULL) {
+  if (!is.list(exp) | !is.list(obs)) {
+    stop("Parameters 'exp' and 'obs' must be lists of 's2dv_cube' objects")
   }
 
-  wrong_input <- FALSE
-  if (!is.list(data)) {
-    wrong_input <- TRUE
+  if (!(all(sapply(exp, inherits, 's2dv_cube')))) {
+     stop("Elements of the list in parameter 'exp' must be of the class ",
+          "'s2dv_cube', as output by CSTools::CST_Load.")
   }
 
-  if (!(all(sapply(data, is_s2dv_object)))) {
-    wrong_input <- TRUE
+  if (!(all(sapply(obs, inherits, 's2dv_cube')))) {
+     stop("Elements of the list in parameter 'obs' must be of the class ",
+          "'s2dv_cube', as output by CSTools::CST_Load.")
   }
 
-  if (wrong_input) {
-    stop("Parameter 'data' must be a list of s2dverification objects as returned by ",
-         "the s2dverification::Load function.")
+  if (length(exp) != length(obs)) {
+    stop("Parameters 'exp' and 'obs' must be of the same length.")
   }
-
-  nvar <- length(data)
+  
+  nvar <- length(exp)
   
   if (nvar < 2) {
-    stop("Parameter 'data' must contain at least two s2dverification objects for two ",
-         "different variables.")
+    stop("Parameters 'exp' and 'obs'  must contain at least two", 
+         " s2dverification objects for two different variables.")
+  }
+  
+  for (j in 1 : nvar) {
+    if (!is.null(names(dim(exp[[j]]$data))) & !is.null(names(dim(obs[[j]]$data)))) {
+      if (all(names(dim(exp[[j]]$data)) %in% names(dim(obs[[j]]$data)))) {
+        dimnames <- names(dim(exp[[j]]$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.")
+    }
   }
   
   if (is.null(weight)) {
     weight <- c(rep(1, nvar))
   } else if (length(weight) != nvar) {
-    stop("Parameter 'weight' must have a length equal to the number of variables.")
+    stop("Parameter 'weight' must have a length equal to the number ",
+         "of variables.")
   }
+  obs_var <- unlist(lapply(obs, function(x) {
+                    x[[which(names(x) == 'Variable')]]}))
+
+  exp_var <- unlist(lapply(exp, function(x) {
+                    x[[which(names(x) == 'Variable')]]}))
   
+  if (all(exp_var != obs_var)) {
+    stop("Variables in parameters 'exp' and 'obs' must be in the same order.")
+  }
   mvrmse <- 0
   sumweights <- 0
-  for (j in 1:nvar) {
-    ano <- Ano_CrossValid(data[[j]]$mod, data[[j]]$obs)
+  for (j in 1 : nvar) {
     # 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)
+    rmse <- s2dverification::RMS(AvgExp, AvgObs, posloop = 1, posRMS = 2, 
+                                 conf = FALSE)
     stdev <- sd(AvgObs)
     mvrmse <- mvrmse + (rmse / stdev * as.numeric(weight[j]))
     sumweights <- sumweights + as.numeric(weight[j])
   }
   mvrmse <- mvrmse / sumweights 
-  data <- data[[1]]
-  data$mvrmse <- mvrmse
-  data$obs <- NULL
-  data$mod <- NULL
-  data <- data[c("mvrmse", names(data)[-which(names(data) == "mvrmse")])]
-  return(data)
+  
+  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) {
+                         x[[which(names(x) == 'source_files')]]}))
+  obs_Datasets <- unlist(lapply(obs, function(x) {
+                         x[[which(names(x) == 'Datasets')]]}))
+  obs_source_files <- unlist(lapply(obs, function(x) {
+                         x[[which(names(x) == 'source_files')]]}))
+
+  exp <- exp[[1]]
+  exp$data <- mvrmse
+  exp$Datasets <- c(exp_Datasets, obs_Datasets)
+  exp$source_files <- c(exp_source_files, obs_source_files)
+  exp$Variable <- c(exp_var)
+  return(exp)
 }
diff --git a/R/CST_RFSlope.R b/R/CST_RFSlope.R
index c1e10db42eb1f70ba3ab98ff0c04bf8fbc006d60..6d7ae8af6a6b001a1a1cd68f0880b6fbc772aadc 100644
--- a/R/CST_RFSlope.R
+++ b/R/CST_RFSlope.R
@@ -4,12 +4,12 @@
 #' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
 #' @description This function computes spatial spectral slopes from a CSTools object
-#' to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object in input.
+#' to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object (of the class 's2dv_cube') as input.
 #' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine.
 #'
-#' @param data CSTools object containing the spatial precipitation fields to downscale.
-#' The data object is expected to have an element named \code{exp} with at least two
+#' @param data An object of the class 's2dv_cube', containing the spatial precipitation fields to downscale.
+#' The data object is expected to have an element named \code{$data} with at least two
 #' spatial dimensions named "lon" and "lat" and one or more dimensions over which
 #' to average these slopes, which can be specified by parameter \code{time_dim}.
 #' @param kmin First wavenumber for spectral slope (default \code{kmin=1}).
@@ -21,17 +21,19 @@
 #' (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes).
 #' The returned array has the same dimensions as the \code{exp} element of the input object,
 #' minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}.
+#' importFrom utils installed.packages
 #' @export
 #' @examples
-#'
 #' #Example using CST_RFSlope for a CSTools object
+#'
+#' \donttest{
 #' exp <- 1 : (2 * 3 * 4 * 8 * 8)
 #' dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
 #' lon <- seq(10, 13.5, 0.5)
 #' dim(lon) <- c(lon = length(lon))
 #' lat <- seq(40, 43.5, 0.5)
 #' dim(lat) <- c(lat = length(lat))
-#' data <- list(exp = exp, lon = lon, lat = lat)
+#' data <- list(data = exp, lon = lon, lat = lat)
 #' slopes <- CST_RFSlope(data)
 #' dim(slopes)
 #' # dataset  member   sdate
@@ -40,11 +42,12 @@
 #' #         [,1]     [,2]     [,3]
 #' #[1,] 1.893503 1.893503 1.893503
 #' #[2,] 1.893503 1.893503 1.893503
+#' }
 
 CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 
-  slopes <- RFSlope(data$exp, kmin, time_dim,
-                     lon_dim = "lon", lat_dim = "lat")
+  slopes <- RFSlope(data$data, kmin, time_dim,
+                    lon_dim = "lon", lat_dim = "lat")
 
   return(slopes)
 }
@@ -81,6 +84,7 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 #' @examples
 #' # Example for the 'reduced' RFSlope function 
 #'
+#' \donttest{
 #' # Create a test array with dimension 8x8 and 20 timesteps, 
 #' # 3 starting dates and 20 ensemble members.
 #' pr <- 1:(4*3*8*8*20)
@@ -98,10 +102,24 @@ CST_RFSlope <- function(data, kmin = 1, time_dim = NULL) {
 #' #[2,] 1.893503 1.893503 1.893503
 #' #[3,] 1.893503 1.893503 1.893503
 #' #[4,] 1.893503 1.893503 1.893503
+#' }
 
 RFSlope <- function(data, kmin = 1, time_dim = NULL,
                     lon_dim = "lon", lat_dim = "lat") {
 
+  # Ensure input  grid is square and with even dimensions
+  if ( (dim(data)[lon_dim] != dim(data)[lat_dim]) |
+       (dim(data)[lon_dim] %% 2 == 1)) {
+    print("Warning: input data are expected to be on a square grid")
+    print("         with an even number of pixels per side")
+    nmin <- min(dim(data)[lon_dim], dim(data)[lat_dim])
+    nmin <- floor(nmin / 2) * 2
+    data <- .subset(data, lat_dim, 1:nmin)
+    data <- .subset(data, lon_dim, 1:nmin)
+    print(paste("The input data have been cut to a square of",
+                nmin, "pixels on each side"))
+  }
+
   # Check/detect time_dim
   if (is.null(time_dim)) {
     time_dim_names <- c("ftime", "sdate", "time")
@@ -134,14 +152,13 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
   # reorder and group time_dim together at the end
   cdim0 <- dim(data)
   imask <- names(cdim0) %in% time_dim
-  data <- aperm(data,c(which(!imask), which(imask)))
+  data <- aperm(data, c(which(!imask), which(imask)))
   cdim <- dim(data)
   ind <- 1:length(which(!imask))
   # compact (multiply) time_dim dimensions
   dim(data) <- c(cdim[ind], rainfarm_samples = prod(cdim[-ind]))
 
   # Repeatedly apply .RFSlope
-  print(dim(data))
   result <- Apply(data, c(lon_dim, lat_dim, "rainfarm_samples"),
                   .RFSlope, kmin)$output1
 
@@ -162,3 +179,20 @@ RFSlope <- function(data, kmin = 1, time_dim = NULL,
                    need_return = "R")
   return(sx)
 }
+
+# Function to generalize through do.call() n-dimensional array subsetting 
+# and array indexing. Derived from Stack Overflow issue
+# https://stackoverflow.com/questions/14500707/select-along-one-of-n-dimensions-in-array
+.subset <- function(field, dim_name, range, drop = FALSE) {
+
+  idim  <- which( names(dim(field)) %in% dim_name )
+  # Create list representing arguments supplied to [
+  # bquote() creates an object corresponding to a missing argument
+  indices <- rep(list(bquote()), length(dim(field)))
+  indices[[idim]] <- range
+
+  # do.call on the indices
+  field <- do.call("[",c(list(field), indices, list(drop = drop)))
+
+  return(field)
+}
diff --git a/R/CST_RainFARM.R b/R/CST_RainFARM.R
index 144cd67fdb5914a095317912487f076ed007fee3..1544a1e2ef314adae6ee941abaa24913df16c892 100644
--- a/R/CST_RainFARM.R
+++ b/R/CST_RainFARM.R
@@ -4,7 +4,8 @@
 #' @author Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 #'
 #' @description This function implements the RainFARM stochastic precipitation
-#' downscaling method and accepts a CSTools object in input.
+#' downscaling method and accepts a CSTools object (an object of the class 
+#' 's2dv_cube' as provided by `CST_Load`) as input.
 #' Adapted for climate downscaling and including orographic correction
 #' as described in Terzago et al. 2018.
 #' References: Terzago, S. et al. (2018). NHESS 18(11), 2825–2840.
@@ -13,11 +14,14 @@
 #' Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, 
 #' which takes a long time only the first time it is called on a machine.
 #'
-#' @param data CSTools object containing the spatial precipitation fields to downscale.
-#' The data object is expected to have an element named \code{exp} with at least two
+#' @param data An object of the class 's2dv_cube' as returned by `CST_Load`, 
+#' containing the spatial precipitation fields to downscale.
+#' The data object is expected to have an element named \code{$data} with at least two
 #' spatial dimensions named "lon" and "lat" and one or more dimensions over which
 #' to compute average spectral slopes (unless specified with parameter \code{slope}),
 #' which can be specified by parameter \code{time_dim}.
+#' The number of longitudes and latitudes in the input data is expected to be even and the same. If not
+#' the function will perform a subsetting to ensure this condition.
 #' @param weights Matrix with climatological weights which can be obtained using
 #' the \code{RFWeights} function. If \code{weights=1.} (default) no weights are used.
 #' The matrix should have dimensions (lon, lat) in this order.
@@ -46,13 +50,17 @@
 #'
 #' 3) if \code{nens>1} and a "member" dimension does not exist: the "realization" dimension is renamed to "member".
 #'
-#' @return CST_RainFARM() returns a downscaled CSTools object.
-#' If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array
-#' after the "member" dimension (unless \code{drop_realization_dim=TRUE} is specified).
-#' The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained.
+#' @return CST_RainFARM() returns a downscaled CSTools object (i.e., of the 
+#' class 's2dv_cube').
+#' If \code{nens>1} an additional dimension named "realization" is added to the 
+#' \code{$data} array after the "member" dimension (unless 
+#' \code{drop_realization_dim=TRUE} is specified).
+#' The ordering of the remaining dimensions in the \code{$data} element of the input object is maintained.
 #' @import multiApply
+#' @importFrom utils installed.packages
 #' @export
 #' @examples
+#' \donttest{
 #' #Example using CST_RainFARM for a CSTools object
 #' nf <- 8   # Choose a downscaling by factor 8
 #' exp <- 1 : (2 * 3 * 4 * 8 * 8)
@@ -61,30 +69,30 @@
 #' dim(lon) <- c(lon = length(lon))
 #' lat <- seq(40, 43.5, 0.5)
 #' dim(lat) <- c(lat = length(lat))
-#' data <- list(exp = exp, lon = lon, lat = lat)
+#' data <- list(data = exp, lon = lon, lat = lat)
 #' # Create a test array of weights
 #' ww <- array(1., dim = c(8 * nf, 8 * nf))
 #' res <- CST_RainFARM(data, nf, ww, nens=3)
 #' str(res)
 #' #List of 3
-#' # $ exp: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ...
+#' # $ data: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ...
 #' # $ lon : num [1:64] 9.78 9.84 9.91 9.97 10.03 ...
 #' # $ lat : num [1:64] 39.8 39.8 39.9 40 40 ...
-#' dim(res$exp)
+#' dim(res$data)
 #' # dataset  member   sdate   ftime     lat     lon    realization
 #' #       1       2       3       4      64      64     3
-#'
+#' }
 CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
                          nens = 1, fglob = FALSE, fsmooth = TRUE,
                          time_dim = NULL, verbose = FALSE,
                          drop_realization_dim = FALSE) {
 
-  res <- RainFARM(data$exp, data$lon, data$lat,
+  res <- RainFARM(data$data, data$lon, data$lat,
                   nf, weights, nens, slope, kmin, fglob, fsmooth,
                   time_dim, lon_dim = "lon", lat_dim = "lat",
                   drop_realization_dim, verbose)
 
-  data$exp <- res$data
+  data$data <- res$data
   data$lon <- res$lon
   data$lat <- res$lat
 
@@ -109,6 +117,8 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters)
 #' and one or more dimensions over which to average these slopes,
 #' which can be specified by parameter \code{time_dim}.
+#' The number of longitudes and latitudes in the input data is expected to be even and the same. If not
+#' the function will perform a subsetting to ensure this condition.
 #' @param lon Vector or array of longitudes.
 #' @param lat Vector or array of latitudes.
 #' @param weights Matrix with climatological weights which can be obtained using
@@ -150,6 +160,7 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' @import multiApply
 #' @export
 #' @examples
+#' \donttest{
 #' # Example for the 'reduced' RainFARM function 
 #' nf <- 8   # Choose a downscaling by factor 8
 #' nens <- 3 # Number of ensemble members
@@ -163,10 +174,8 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' ww <- array(1., dim = c(8 * nf, 8 * nf))
 #' # ... or create proper weights using an external fine-scale climatology file
 #' #     Specify a weightsfn filename if you wish to save the weights
-#' \donttest{
 #' ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, 
 #'                 weightsfn = "", fsmooth = TRUE)
-#' }
 #' # downscale using weights (ww=1. means do not use weights)
 #' res <- RainFARM(pr, lon_mat, lat_mat, nf, 
 #'                 fsmooth = TRUE, fglob = FALSE, 
@@ -179,12 +188,41 @@ CST_RainFARM <- function(data, nf, weights = 1., slope = 0, kmin = 1,
 #' dim(res$data)
 #' #  lon         lat       ftime realization 
 #' #   64          64          20           2 
+#' }
 #'
 RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
                      slope = 0, kmin = 1, fglob = FALSE, fsmooth = TRUE,
                      time_dim = NULL, lon_dim = "lon", lat_dim = "lat",
                      drop_realization_dim = FALSE, verbose = FALSE) {
 
+  # Ensure input  grid is square and with even dimensions
+  if ( (dim(data)[lon_dim] != dim(data)[lat_dim]) |
+       (dim(data)[lon_dim] %% 2 == 1)) {
+    print("Warning: input data are expected to be on a square grid")
+    print("         with an even number of pixels per side")
+    nmin <- min(dim(data)[lon_dim], dim(data)[lat_dim])
+    nmin <- floor(nmin / 2) * 2
+    data <- .subset(data, lat_dim, 1:nmin)
+    data <- .subset(data, lon_dim, 1:nmin)
+    if (length(dim(lon)) == 2) {
+      lon <- lon[1:nmin, 1:nmin]
+      lat <- lat[1:nmin, 1:nmin]
+    } else {
+      lon <- lon[1:nmin]
+      lat <- lat[1:nmin]
+    }
+    print("The input data have been cut to the range")
+    print(paste0("lon: [", lon[1], ", ", lon[length(lon)], "] ",
+                 " lat: [", lat[1], ", ", lat[length(lat)], "]"))
+  }
+  if (!(length(dim(weights)) == 0) &
+      !(all(dim(weights) == c(dim(data)[lon_dim] * nf,
+                              dim(data)[lat_dim] * nf)))) {
+    stop(paste("The dimensions of the weights matrix (", dim(weights)[1],
+               "x", dim(weights)[2] ,
+               ") are not consistent with the size of the data (",
+               dim(data)[lon_dim], ") and the refinement factor (", nf, ")"))
+  }
   # Check/detect time_dim
   if (is.null(time_dim)) {
     time_dim_names <- c("ftime", "sdate", "time")
@@ -201,7 +239,7 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
       stop("Could not automatically detect a target time dimension ",
            "in the provided data in 'data'.")
     }
-    print(paste("Selected time dim: ", time_dim))
+    print(paste("Selected time dim:", time_dim))
   }
 
   # Check availability of JuliaCall
@@ -236,21 +274,26 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
   # Expand back rainfarm_samples to compacted dims
   dim(result) <- c(dim(result)[1:2], cdim[-ind], dim(result)[-(1:3)])
 
-  # Reorder as it was in original data + realization dim after member if it exists
+  # Reorder as it was in original data
+  # + realization dim after member if it exists
   ienspos <- which(names(cdim0) == "member")
-  if( length(ienspos) == 0) ienspos <- length(names(cdim0))
-  iorder <- sapply(c(names(cdim0)[1:ienspos], "realization", names(cdim0)[-(1:ienspos)]), grep, names(dim(result)))
+  if ( length(ienspos) == 0) ienspos <- length(names(cdim0))
+  iorder <- sapply(c(names(cdim0)[1:ienspos], "realization",
+                     names(cdim0)[-(1:ienspos)]),
+                   grep, names(dim(result)))
   result <- aperm(result, iorder)
 
   if (drop_realization_dim) {
-    cdim = dim(result)
-    if (nens == 1) { 
+    cdim <- dim(result)
+    if (nens == 1) {
       # just drop it if only one member
       dim(result) <- cdim[-which(names(cdim) == "realization")[1]]
-    } else if("member" %in% names(cdim)) {
-      # compact member and realization dimension if member dim exists, else rename realization to member
+    } else if ("member" %in% names(cdim)) {
+      # compact member and realization dimension if member dim exists,
+      # else rename realization to member
       ind <- which(names(cdim) %in% c("member", "realization"))
-      dim(result) <- c(cdim[1:(ind[1]-1)], cdim[ind[1]]*cdim[ind[2]], cdim[(ind[2]+1):length(cdim)])
+      dim(result) <- c(cdim[1:(ind[1] - 1)], cdim[ind[1]] * cdim[ind[2]],
+                       cdim[(ind[2] + 1):length(cdim)])
     } else {
       ind <- which(names(cdim) %in% "realization")
       names(dim(result))[ind] <- "member"
@@ -296,3 +339,20 @@ RainFARM <- function(data, lon, lat, nf, weights = 1., nens = 1,
   }
   return(r)
 }
+
+# Function to generalize through do.call() n-dimensional array subsetting 
+# and array indexing. Derived from Stack Overflow issue
+# https://stackoverflow.com/questions/14500707/select-along-one-of-n-dimensions-in-array
+.subset <- function(field, dim_name, range, drop = FALSE) {
+
+  idim  <- which( names(dim(field)) %in% dim_name )
+  # Create list representing arguments supplied to [
+  # bquote() creates an object corresponding to a missing argument
+  indices <- rep(list(bquote()), length(dim(field)))
+  indices[[idim]] <- range
+
+  # do.call on the indices
+  field <- do.call("[",c(list(field), indices, list(drop = drop)))
+
+  return(field)
+}
diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R
index 83b4e239734d85662b73d7d92fe998e0a766e7c8..73afb63909de63565ec44845aa4049d43dc4e05b 100644
--- a/R/PlotForecastPDF.R
+++ b/R/PlotForecastPDF.R
@@ -25,17 +25,23 @@
 #'@importFrom plyr dlply
 #'
 #'@examples
+#'\donttest{
 #'fcsts <- data.frame(fcst1 = rnorm(50), fcst2 = rnorm(50, 0.5, 1.2), 
 #'                    fcst3 = rnorm(50, -0.5, 0.9))
 #'PlotForecastPDF(fcsts,c(-1,1))
 #'fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5))
 #'PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), 
 #'                fcst.names = paste0('random fcst ', 1 : 5), obs = 0.7)
+#'}
 #'@export
 PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = NULL, 
     plotfile = NULL, title = "Set a title", var.name = "Varname (units)", fcst.names = NULL, 
     add.ensmemb = c("above", "below", "no"), color.set = c("ggplot", "s2s4e", "hydro")) {
     value <- init <- extremes <- x <- ymin <- ymax <- tercile <- y <- xend <- yend <- yjitter <- MLT <- NULL
+    ggColorHue <- function(n) {
+      hues <- seq(15, 375, length = n + 1)
+      hcl(h = hues, l = 65, c = 100)[1:n]
+    }
     #------------------------
     # Define color sets
     #------------------------
@@ -53,10 +59,6 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         colorObs <- "purple"
         colorLab <- c("darkorange3", "blue")
     } else if (color.set == "ggplot") {
-        ggColorHue <- function(n) {
-            hues <- seq(15, 375, length = n + 1)
-            hcl(h = hues, l = 65, c = 100)[1:n]
-        }
         colorFill <- rev(ggColorHue(3))
         colorHatch <- c("deepskyblue3", "indianred1")
         colorMember <- c("#ffff7f")
diff --git a/R/PlotMostLikelyQuantileMap.R b/R/PlotMostLikelyQuantileMap.R
index 193bc712b00e1046a350541dbb3cc331b3307076..c8abb372122bd0ea915a39217c5aa84d8682d3ff 100644
--- a/R/PlotMostLikelyQuantileMap.R
+++ b/R/PlotMostLikelyQuantileMap.R
@@ -14,6 +14,8 @@
 #' 
 #'@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
 #'@examples
 #'# Simple example
 #'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
diff --git a/R/RFWeights.R b/R/RFWeights.R
index 66b592e187893a6f0d2f63805b8138f6c7dd345c..0459e285029489d210a4800ae491f27b082342a2 100644
--- a/R/RFWeights.R
+++ b/R/RFWeights.R
@@ -21,8 +21,10 @@
 #' can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org)
 #' websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org).
 #' @param nf Refinement factor for downscaling (the output resolution is increased by this factor).
-#' @param lon Vector or array of longitudes (can be omitted if \code{reffile} is used).
-#' @param lat Vector or array of latitudes (can be omitted if \code{reffile} is used).
+#' @param lon Vector of longitudes (can be omitted if \code{reffile} is used).
+#' @param lat Vector of latitudes (can be omitted if \code{reffile} is used).
+#' The number of longitudes and latitudes is expected to be even and the same. If not
+#' the function will perform a subsetting to ensure this condition.
 #' @param reffile Filename of reference file (e.g. the file to downscale).
 #' Not needed if \code{lon} and \code{lat} are passed.
 #' @param varname Name of the variable to be read from \code{reffile}.
@@ -37,7 +39,7 @@
 #'
 #' \donttest{
 #' # Specify lon and lat of the input
-#' lon_mat <- seq(10,13.5,0.5) # could also be a 2d matrix
+#' lon_mat <- seq(10,13.5,0.5)
 #' lat_mat <- seq(40,43.5,0.5)
 #'
 #' nf <- 8
@@ -57,15 +59,27 @@ RFWeights <- function(orofile, nf, lon=NULL, lat=NULL, reffile="",
      stop("Neither a reference filename nor lon and lat specified!\n")
   }
   if (reffile == "") {
-     xvar <- ncdim_def("lon", "degrees_east", lon, longname = "longitude")
-     yvar <- ncdim_def("lat", "degrees_north", lat, longname = "latitude")
-     tvar <- ncdim_def("time", "days since 01-01-2000", 0., longname = "time",
-                       unlim = T)
-     my_ncdf <- ncvar_def("grid", "m", list(xvar, yvar, tvar), -999,
-                          longname = "grid", prec = "single")
-     reffile <- tempfile()
-     ncfile <- nc_create(reffile, my_ncdf, force_v4 = FALSE)
-     nc_close(ncfile)
+    # Ensure input  grid is square and with even dimensions
+    if ((length(lat) != length(lon)) | (length(lon) %% 2 == 1)) {
+      print("Warning: input data are expected to be on a square grid")
+      print("         with an even number of pixels per side")
+      nmin <- min(length(lon), length(lat))
+      nmin <- floor(nmin / 2) * 2
+      lon <- lon[1:nmin]
+      lat <- lat[1:nmin]
+      print("The input data have been cut to the range")
+      print(paste0("lon: [", lon[1], ", ", lon[length(lon)], "] ",
+                   " lat: [", lat[1], ", ", lat[length(lat)], "]"))
+    }
+    xvar <- ncdim_def("lon", "degrees_east", lon, longname = "longitude")
+    yvar <- ncdim_def("lat", "degrees_north", lat, longname = "latitude")
+    tvar <- ncdim_def("time", "days since 01-01-2000", 0., longname = "time",
+                      unlim = T)
+    my_ncdf <- ncvar_def("grid", "m", list(xvar, yvar, tvar), -999,
+                         longname = "grid", prec = "single")
+    reffile <- tempfile()
+    ncfile <- nc_create(reffile, my_ncdf, force_v4 = FALSE)
+    nc_close(ncfile)
   }
 
   ww <- julia_call("rfweights", orofile, reffile, as.integer(nf),
diff --git a/R/sample_data.R b/R/sample_data.R
new file mode 100644
index 0000000000000000000000000000000000000000..3ceb7eedee75cd3beb1a7ae692b52866ea38c683
--- /dev/null
+++ b/R/sample_data.R
@@ -0,0 +1,77 @@
+#' Sample Of Experimental And Observational Climate Data In Function Of Longitudes And Latitudes
+#'
+#' This sample data set contains gridded seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). The data was generated on (or interpolated onto, for the reconstruction) a rectangular regular grid of size 360 by 181.
+#' 
+#' It is recommended to use the data set as follows:
+#' ```
+#' require(zeallot)
+#' c(exp, obs) %<-% CSTools::lonlat_data
+#' ```
+#'
+#' The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system.
+#' ```
+#' library(CSTools)
+#' require(zeallot)
+#'
+#' startDates <- c('20001101', '20011101', '20021101',
+#'                 '20031101', '20041101', '20051101')
+#'
+#' c(exp, obs) %<-% 
+#'   CST_Load(
+#'     var = 'tas', 
+#'     exp = 'system5c3s', 
+#'     obs = 'era5', 
+#'     nmember = 15,
+#'     sdates = startDates,
+#'     leadtimemax = 3,
+#'     latmin = 27, latmax = 48,
+#'     lonmin = -12, lonmax = 40, 
+#'     output = 'lonlat',
+#'     nprocs = 1
+#'   )
+#' ```
+#'
+#' @name lonlat_data
+#' @docType data
+#' @author Nicolau Manubens \email{nicolau.manubens@bsc.es}
+#' @keywords data
+NULL
+
+#' Sample Of Experimental And Observational Climate Data Averaged Over A Region
+#'
+#' This sample data set contains area-averaged seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E).
+#' 
+#' It is recommended to use the data set as follows:
+#' ```
+#' require(zeallot)
+#' c(exp, obs) %<-% CSTools::areave_data
+#' ```
+#'
+#' The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system.
+#' ```
+#' library(CSTools)
+#' require(zeallot)
+#'
+#' startDates <- c('20001101', '20011101', '20021101',
+#'                 '20031101', '20041101', '20051101')
+#'
+#' c(exp, obs) %<-% 
+#'   CST_Load(
+#'     var = 'tas', 
+#'     exp = 'system5c3s', 
+#'     obs = 'era5', 
+#'     nmember = 15,
+#'     sdates = startDates,
+#'     leadtimemax = 3,
+#'     latmin = 27, latmax = 48,
+#'     lonmin = -12, lonmax = 40, 
+#'     output = 'areave',
+#'     nprocs = 1
+#'   )
+#' ```
+#'
+#' @name areave_data
+#' @docType data
+#' @author Nicolau Manubens \email{nicolau.manubens@bsc.es}
+#' @keywords data
+NULL
diff --git a/data/areave_data.RData b/data/areave_data.RData
new file mode 100644
index 0000000000000000000000000000000000000000..426d4e165fada19e79a8a6912cc78124b24fd9f1
Binary files /dev/null and b/data/areave_data.RData differ
diff --git a/data/lonlat_data.RData b/data/lonlat_data.RData
new file mode 100644
index 0000000000000000000000000000000000000000..137786f0fd5c6b2825df17d451a85927ff54eb6f
Binary files /dev/null and b/data/lonlat_data.RData differ
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_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd
index 99602261c785d2b577d8988589d05addc5960611..a1b415fb525a9f2e3c72171e631c2633e0aefcd6 100644
--- a/man/CST_BiasCorrection.Rd
+++ b/man/CST_BiasCorrection.Rd
@@ -4,13 +4,15 @@
 \alias{CST_BiasCorrection}
 \title{Bias Correction based on the mean and standard deviation adjustment}
 \usage{
-CST_BiasCorrection(data)
+CST_BiasCorrection(exp, obs)
 }
 \arguments{
-\item{data}{CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).}
+\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}.}
 }
 \value{
-a CSTools object (s2dverification object) with the bias corrected forecasts a element called \code{data$biascalibration}.
+an object of class \code{s2dv_cube} containing the bias corrected forecasts in the element called \code{$data} with the same dimensions of the experimental data.
 }
 \description{
 This function applies the simple bias adjustment technique described in Torralba et al. (2017). The adjusted forecasts have an equivalent standard deviation and mean to that of the reference dataset.
@@ -26,8 +28,11 @@ obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
 dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
 lon <- seq(0, 30, 5)
 lat <- seq(0, 25, 5)
-data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-a <- CST_BiasCorrection(data1)
+exp <- list(data = mod1, lat = lat, lon = lon)
+obs <- list(data = obs1, lat = lat, lon = lon)
+attr(exp, 'class') <- 's2dv_cube'
+attr(obs, 'class') <- 's2dv_cube'
+a <- CST_BiasCorrection(exp = exp, obs = obs)
 str(a)
 }
 \references{
diff --git a/man/CST_Calibration.Rd b/man/CST_Calibration.Rd
index b324e465b2f975077aa8e7bcccfdba82b03c121b..e3e402a5a0c4eeea83c2c6e3799da290dd4539a0 100644
--- a/man/CST_Calibration.Rd
+++ b/man/CST_Calibration.Rd
@@ -4,35 +4,35 @@
 \alias{CST_Calibration}
 \title{Forecast Calibration based on the ensemble inflation}
 \usage{
-CST_Calibration(data)
+CST_Calibration(exp, obs)
 }
 \arguments{
-\item{data}{a CSTools object (an s2dverification object as output by the \code{Load} function from the s2dverification package).}
+\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}.}
 }
 \value{
-a CSTools object (s2dverification object) with the calibrated forecasts in a element called \code{data$calibration}.
+an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions of the experimental data.
 }
 \description{
 This function applies a variance inflation technique described in Doblas-Reyes et al. (2005) in leave-one-out cross-validation. This bias adjustment method produces calibrated forecasts with equivalent mean and variance to that of the reference dataset, but at the same time preserve reliability.
 }
 \examples{
+\donttest{
 # Example
-# Creation of sample s2dverification objects. These are not complete
-# s2dverification objects though. The Load function returns complete objects.
-mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
-dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
-obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
-dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
-lon <- seq(0, 30, 5)
-lat <- seq(0, 25, 5)
-data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-a <- CST_Calibration(data1)
-str(a)
-
+# Load data using CST_Load or use the sample data provided:
+library(zeallot)
+c(exp, obs) \%<-\% lonlat_data
+exp_calibrated <- CST_Calibration(exp = exp, obs = obs)
+str(exp_calibrated)
+}
 }
 \references{
 Doblas-Reyes F.J, Hagedorn R, Palmer T.N. The rationale behind the success of multi-model ensembles in seasonal forecasting—II calibration and combination. Tellus A. 2005;57:234–252. doi:10.1111/j.1600-0870.2005.00104.x
 }
+\seealso{
+\code{\link{CST_Load}}
+}
 \author{
 Verónica Torralba, \email{veronica.torralba@bsc.es}
 }
diff --git a/man/CST_Load.Rd b/man/CST_Load.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..bf03ba42d98810643a98bfca27f4563f6361162e
--- /dev/null
+++ b/man/CST_Load.Rd
@@ -0,0 +1,49 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_Load.R
+\name{CST_Load}
+\alias{CST_Load}
+\title{CSTools Data Retreival Function}
+\usage{
+CST_Load(...)
+}
+\arguments{
+\item{...}{Parameters that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`.}
+}
+\value{
+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.
+}
+\description{
+This function aggregates, subsets and retrieves sub-seasonal, seasonal, decadal or climate projection data from NetCDF files in a local file system or on remote OPeNDAP servers, and arranges it for easy application of the CSTools functions.
+}
+\details{
+It receives any number of parameters (`...`) that are automatically forwarded to the `s2dverification::Load` function. See details in `?s2dverification::Load`.
+
+It is recommended to use this function in combination with the `zeallot::"%<-%"` operator, to directly assing the two returned 's2dv_cube's to two separate variables, which can then be sent independently to other functions in CSTools as needed. E.g.: `c(exp, obs) <- CST_Load(...)`.
+}
+\examples{
+\dontrun{
+library(zeallot)
+startDates <- c('20001101', '20011101', '20021101', 
+                '20031101', '20041101', '20051101')
+c(exp, obs) \%<-\% 
+  CST_Load(
+    var = 'tas', 
+    exp = 'system5c3s', 
+    obs = 'era5', 
+    nmember = 15,
+    sdates = startDates,
+    leadtimemax = 3,
+    latmin = 27, latmax = 48,
+    lonmin = -12, lonmax = 40,
+    output = 'lonlat',
+    nprocs = 1
+  )
+}
+\dontshow{
+exp <- CSTools::lonlat_data$exp
+obs <- CSTools::lonlat_data$obs
+}
+}
+\author{
+Nicolau Manubens, \email{nicolau.manubens@bsc.es}
+}
diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd
index 70407776600012043ee1974ec28bcd4f67bdf570..99fdb400f1a82dd7267cd284a2b6597dfb314797 100644
--- a/man/CST_MultiMetric.Rd
+++ b/man/CST_MultiMetric.Rd
@@ -4,47 +4,44 @@
 \alias{CST_MultiMetric}
 \title{Multiple Metrics applied in Multiple Model Anomalies}
 \usage{
-CST_MultiMetric(data, metric = "correlation", multimodel = TRUE,
-  names = NULL)
+CST_MultiMetric(exp, obs, metric = "correlation", multimodel = TRUE)
 }
 \arguments{
-\item{data}{a CSTools object (list) giving as output of \code{Load} function from S2dverification package.}
+\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_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{multimodel}{a logical value indicating whether a Multi-Model Mean should be computed.}
-
-\item{names}{a character vector indicating the names of the models or experiments to be compared.}
 }
 \value{
-A CSTools object as the input parameter \code{data}, without the elements \code{mod} and \code{obs} and including:
-\itemize{
-\item\code{$ano_exp} {an array with the same dimensions as the \code{mod} label in the input \code{data}.}
-\item\code{$ano_exp} {an array with the same dimensions as the \code{obs} label in the input \code{data}.}
-\item\code{$metric} {An array with two datset dimensions from the \code{mod} and \code{obs} label in the input \code{data}. 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 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.
 }
 \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.
 }
-\details{
-The output will include attributes in element \code{ano_exp} returning a vector with the names of the models or experiments.
-}
 \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)
 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)
-dat <- list(mod = mod, obs = obs, lat = lat, lon = lon)
-a <- CST_MultiMetric(dat)
+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'
+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{
-\url{http://link.springer.com/10.1007/s00382-018-4404-z}
+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[s2dverification]{Load}}
+\code{\link[s2dverification]{Corr}}, \code{\link[s2dverification]{RMS}}, \code{\link[s2dverification]{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 65cdfa952c1c0bd431413f122eb43af08e7e1895..24af608c8360985de73763434a16f1a11fd35ffe 100644
--- a/man/CST_MultivarRMSE.Rd
+++ b/man/CST_MultivarRMSE.Rd
@@ -4,15 +4,17 @@
 \alias{CST_MultivarRMSE}
 \title{Multivariate Root Mean Square Error (RMSE)}
 \usage{
-CST_MultivarRMSE(data, weight = NULL)
+CST_MultivarRMSE(exp, obs, weight = NULL)
 }
 \arguments{
-\item{data}{a list of s2dverification objects (lists) as output by the \code{Load} function from the s2dverification package, one for each variable.}
+\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_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.}
 }
 \value{
-\code{$mvrmse} {An array with dimensions: c(number of exp, number of obs, 1 (the multivariate RMSE value), number of lat, number of lon)}
+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)
 }
 \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).
@@ -20,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)
@@ -30,15 +33,32 @@ 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)
-data1 <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-data2 <- list(mod = mod2, obs = obs2, lat = lat, lon = lon)
-data <- list(data1, data2)
+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'
+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'
+
+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(data, weight)
+a <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight = weight)
 str(a)
 }
+\seealso{
+\code{\link[s2dverification]{RMS}} and \code{\link{CST_Load}}
+}
 \author{
 Deborah Verfaillie, \email{deborah.verfaillie@bsc.es}
 }
diff --git a/man/CST_RFSlope.Rd b/man/CST_RFSlope.Rd
index a931374068f7ce9b762f76641faae54605e4420d..489e7a97ed29f57e308020e71cc0e330dba09761 100644
--- a/man/CST_RFSlope.Rd
+++ b/man/CST_RFSlope.Rd
@@ -7,8 +7,8 @@
 CST_RFSlope(data, kmin = 1, time_dim = NULL)
 }
 \arguments{
-\item{data}{CSTools object containing the spatial precipitation fields to downscale.
-The data object is expected to have an element named \code{exp} with at least two
+\item{data}{An object of the class 's2dv_cube', containing the spatial precipitation fields to downscale.
+The data object is expected to have an element named \code{$data} with at least two
 spatial dimensions named "lon" and "lat" and one or more dimensions over which
 to average these slopes, which can be specified by parameter \code{time_dim}.}
 
@@ -24,23 +24,25 @@ CST_RFSlope() returns spectral slopes using the RainFARM convention
 (the logarithmic slope of k*|A(k)|^2 where A(k) are the spectral amplitudes).
 The returned array has the same dimensions as the \code{exp} element of the input object,
 minus the dimensions specified by \code{lon_dim}, \code{lat_dim} and \code{time_dim}.
+importFrom utils installed.packages
 }
 \description{
 This function computes spatial spectral slopes from a CSTools object
-to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object in input.
+to be used for RainFARM stochastic precipitation downscaling method and accepts a CSTools object (of the class 's2dv_cube') as input.
 Please notice: The function uses the \code{julia_setup()} function from the JuliaCall library, 
 which takes a long time only the first time it is called on a machine.
 }
 \examples{
-
 #Example using CST_RFSlope for a CSTools object
+
+\donttest{
 exp <- 1 : (2 * 3 * 4 * 8 * 8)
 dim(exp) <- c(dataset = 1, member = 2, sdate = 3, ftime = 4, lat = 8, lon = 8)
 lon <- seq(10, 13.5, 0.5)
 dim(lon) <- c(lon = length(lon))
 lat <- seq(40, 43.5, 0.5)
 dim(lat) <- c(lat = length(lat))
-data <- list(exp = exp, lon = lon, lat = lat)
+data <- list(data = exp, lon = lon, lat = lat)
 slopes <- CST_RFSlope(data)
 dim(slopes)
 # dataset  member   sdate
@@ -50,6 +52,7 @@ slopes
 #[1,] 1.893503 1.893503 1.893503
 #[2,] 1.893503 1.893503 1.893503
 }
+}
 \author{
 Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 }
diff --git a/man/CST_RainFARM.Rd b/man/CST_RainFARM.Rd
index 2120194e1eee8ac44ef3525fbf01e9511b0cbd16..5343af7aaf8599ecd997cdb300576f2ee1fcf100 100644
--- a/man/CST_RainFARM.Rd
+++ b/man/CST_RainFARM.Rd
@@ -9,11 +9,14 @@ CST_RainFARM(data, nf, weights = 1, slope = 0, kmin = 1, nens = 1,
   drop_realization_dim = FALSE)
 }
 \arguments{
-\item{data}{CSTools object containing the spatial precipitation fields to downscale.
-The data object is expected to have an element named \code{exp} with at least two
+\item{data}{An object of the class 's2dv_cube' as returned by `CST_Load`, 
+containing the spatial precipitation fields to downscale.
+The data object is expected to have an element named \code{$data} with at least two
 spatial dimensions named "lon" and "lat" and one or more dimensions over which
 to compute average spectral slopes (unless specified with parameter \code{slope}),
-which can be specified by parameter \code{time_dim}.}
+which can be specified by parameter \code{time_dim}.
+The number of longitudes and latitudes in the input data is expected to be even and the same. If not
+the function will perform a subsetting to ensure this condition.}
 
 \item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).}
 
@@ -53,14 +56,17 @@ with the following behaviour if set to TRUE:
 3) if \code{nens>1} and a "member" dimension does not exist: the "realization" dimension is renamed to "member".}
 }
 \value{
-CST_RainFARM() returns a downscaled CSTools object.
-If \code{nens>1} an additional dimension named "realization" is added to the \code{exp} array
-after the "member" dimension (unless \code{drop_realization_dim=TRUE} is specified).
-The ordering of the remaining dimensions in the \code{exp} element of the input object is maintained.
+CST_RainFARM() returns a downscaled CSTools object (i.e., of the 
+class 's2dv_cube').
+If \code{nens>1} an additional dimension named "realization" is added to the 
+\code{$data} array after the "member" dimension (unless 
+\code{drop_realization_dim=TRUE} is specified).
+The ordering of the remaining dimensions in the \code{$data} element of the input object is maintained.
 }
 \description{
 This function implements the RainFARM stochastic precipitation
-downscaling method and accepts a CSTools object in input.
+downscaling method and accepts a CSTools object (an object of the class 
+'s2dv_cube' as provided by `CST_Load`) as input.
 Adapted for climate downscaling and including orographic correction
 as described in Terzago et al. 2018.
 References: Terzago, S. et al. (2018). NHESS 18(11), 2825–2840.
@@ -70,6 +76,7 @@ Please notice: The function uses the \code{julia_setup()} function from the Juli
 which takes a long time only the first time it is called on a machine.
 }
 \examples{
+\donttest{
 #Example using CST_RainFARM for a CSTools object
 nf <- 8   # Choose a downscaling by factor 8
 exp <- 1 : (2 * 3 * 4 * 8 * 8)
@@ -78,19 +85,19 @@ lon <- seq(10, 13.5, 0.5)
 dim(lon) <- c(lon = length(lon))
 lat <- seq(40, 43.5, 0.5)
 dim(lat) <- c(lat = length(lat))
-data <- list(exp = exp, lon = lon, lat = lat)
+data <- list(data = exp, lon = lon, lat = lat)
 # Create a test array of weights
 ww <- array(1., dim = c(8 * nf, 8 * nf))
 res <- CST_RainFARM(data, nf, ww, nens=3)
 str(res)
 #List of 3
-# $ exp: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ...
+# $ data: num [1:4, 1:64, 1:64, 1, 1:2, 1:3] 201 115 119 307 146 ...
 # $ lon : num [1:64] 9.78 9.84 9.91 9.97 10.03 ...
 # $ lat : num [1:64] 39.8 39.8 39.9 40 40 ...
-dim(res$exp)
+dim(res$data)
 # dataset  member   sdate   ftime     lat     lon    realization
 #       1       2       3       4      64      64     3
-
+}
 }
 \author{
 Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
diff --git a/man/PlotForecastPDF.Rd b/man/PlotForecastPDF.Rd
index 0af93c072220be143acb582921d3aff2a23a5330..8b8bae0f2feaa99cbe3e855bc7c81a9d526df6cc 100644
--- a/man/PlotForecastPDF.Rd
+++ b/man/PlotForecastPDF.Rd
@@ -38,6 +38,7 @@ a ggplot object containing the plot.
 This function plots the probability distribution function of several ensemble forecasts for the same event, either initialized at different moments or by different models. Probabilities for tercile categories are computed, plotted in colors and annotated. An asterisk marks the tercile with higher probabilities. Probabilities for extreme categories (above P90 and below P10) can also be included as hatched areas. Individual ensemble members can be plotted as jittered points. The observed value is optionally shown as a diamond.
 }
 \examples{
+\donttest{
 fcsts <- data.frame(fcst1 = rnorm(50), fcst2 = rnorm(50, 0.5, 1.2), 
                    fcst3 = rnorm(50, -0.5, 0.9))
 PlotForecastPDF(fcsts,c(-1,1))
@@ -45,6 +46,7 @@ fcsts2 <- array(rnorm(100), dim = c(members = 20, fcst = 5))
 PlotForecastPDF(fcsts2, c(-0.66, 0.66), extreme.limits = c(-1.2, 1.2), 
                fcst.names = paste0('random fcst ', 1 : 5), obs = 0.7)
 }
+}
 \author{
 Llorenç Lledó \email{llledo@bsc.es}
 }
diff --git a/man/RFSlope.Rd b/man/RFSlope.Rd
index 1f5689aac73e0a5b2489bf25eb0bf2cc98412e9d..628838b7b84aaa7de866c54a5c5671f2ff4af5fe 100644
--- a/man/RFSlope.Rd
+++ b/man/RFSlope.Rd
@@ -42,6 +42,7 @@ which takes a long time only the first time it is called on a machine.
 \examples{
 # Example for the 'reduced' RFSlope function 
 
+\donttest{
 # Create a test array with dimension 8x8 and 20 timesteps, 
 # 3 starting dates and 20 ensemble members.
 pr <- 1:(4*3*8*8*20)
@@ -60,6 +61,7 @@ slopes
 #[3,] 1.893503 1.893503 1.893503
 #[4,] 1.893503 1.893503 1.893503
 }
+}
 \author{
 Jost von Hardenberg - ISAC-CNR, \email{j.vonhardenberg@isac.cnr.it}
 }
diff --git a/man/RFWeights.Rd b/man/RFWeights.Rd
index 1f90924f2f5aaa0dd9197c11d4c5d01f799cf652..073cda877ce3926ef14d5ec03a0dfb396ec97186 100644
--- a/man/RFWeights.Rd
+++ b/man/RFWeights.Rd
@@ -20,9 +20,11 @@ websites. The latter data will need to be converted to NetCDF format before bein
 
 \item{nf}{Refinement factor for downscaling (the output resolution is increased by this factor).}
 
-\item{lon}{Vector or array of longitudes (can be omitted if \code{reffile} is used).}
+\item{lon}{Vector of longitudes (can be omitted if \code{reffile} is used).}
 
-\item{lat}{Vector or array of latitudes (can be omitted if \code{reffile} is used).}
+\item{lat}{Vector of latitudes (can be omitted if \code{reffile} is used).
+The number of longitudes and latitudes is expected to be even and the same. If not
+the function will perform a subsetting to ensure this condition.}
 
 \item{reffile}{Filename of reference file (e.g. the file to downscale).
 Not needed if \code{lon} and \code{lat} are passed.}
@@ -53,7 +55,7 @@ which takes a long time only the first time it is called on a machine
 
 \donttest{
 # Specify lon and lat of the input
-lon_mat <- seq(10,13.5,0.5) # could also be a 2d matrix
+lon_mat <- seq(10,13.5,0.5)
 lat_mat <- seq(40,43.5,0.5)
 
 nf <- 8
diff --git a/man/RainFARM.Rd b/man/RainFARM.Rd
index b02b4102fe3868513c41aaccf41ba0f70086d64b..34a4d418fcfb50c9bec63fc427c46177feb8b212 100644
--- a/man/RainFARM.Rd
+++ b/man/RainFARM.Rd
@@ -14,7 +14,9 @@ RainFARM(data, lon, lat, nf, weights = 1, nens = 1, slope = 0,
 The input array is expected to have at least two dimensions named "lon" and "lat" by default
 (these default names can be changed with the \code{lon_dim} and \code{lat_dim} parameters)
 and one or more dimensions over which to average these slopes,
-which can be specified by parameter \code{time_dim}.}
+which can be specified by parameter \code{time_dim}.
+The number of longitudes and latitudes in the input data is expected to be even and the same. If not
+the function will perform a subsetting to ensure this condition.}
 
 \item{lon}{Vector or array of longitudes.}
 
@@ -80,6 +82,7 @@ Please notice: The function uses the \code{julia_setup()} function from the Juli
 which takes a long time only the first time it is called on a machine.
 }
 \examples{
+\donttest{
 # Example for the 'reduced' RainFARM function 
 nf <- 8   # Choose a downscaling by factor 8
 nens <- 3 # Number of ensemble members
@@ -93,10 +96,8 @@ lat_mat <- seq(40, 43.5, 0.5)
 ww <- array(1., dim = c(8 * nf, 8 * nf))
 # ... or create proper weights using an external fine-scale climatology file
 #     Specify a weightsfn filename if you wish to save the weights
-\donttest{
 ww <- RFWeights("./worldclim.nc", nf, lon = lon_mat, lat = lat_mat, 
                 weightsfn = "", fsmooth = TRUE)
-}
 # downscale using weights (ww=1. means do not use weights)
 res <- RainFARM(pr, lon_mat, lat_mat, nf, 
                 fsmooth = TRUE, fglob = FALSE, 
@@ -109,6 +110,7 @@ str(res)
 dim(res$data)
 #  lon         lat       ftime realization 
 #   64          64          20           2 
+}
 
 }
 \author{
diff --git a/man/areave_data.Rd b/man/areave_data.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d42663944e89e94abe084df648badc329e2778de
--- /dev/null
+++ b/man/areave_data.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_data.R
+\docType{data}
+\name{areave_data}
+\alias{areave_data}
+\title{Sample Of Experimental And Observational Climate Data Averaged Over A Region}
+\description{
+This sample data set contains area-averaged seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E).
+}
+\details{
+It is recommended to use the data set as follows:
+```
+require(zeallot)
+c(exp, obs) %<-% CSTools::areave_data
+```
+
+The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system.
+```
+library(CSTools)
+require(zeallot)
+
+startDates <- c('20001101', '20011101', '20021101',
+                '20031101', '20041101', '20051101')
+
+c(exp, obs) %<-% 
+  CST_Load(
+    var = 'tas', 
+    exp = 'system5c3s', 
+    obs = 'era5', 
+    nmember = 15,
+    sdates = startDates,
+    leadtimemax = 3,
+    latmin = 27, latmax = 48,
+    lonmin = -12, lonmax = 40, 
+    output = 'areave',
+    nprocs = 1
+  )
+```
+}
+\author{
+Nicolau Manubens \email{nicolau.manubens@bsc.es}
+}
+\keyword{data}
diff --git a/man/lonlat_data.Rd b/man/lonlat_data.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..a835ad89ef8890fc136265f7ee721481baaa81a6
--- /dev/null
+++ b/man/lonlat_data.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sample_data.R
+\docType{data}
+\name{lonlat_data}
+\alias{lonlat_data}
+\title{Sample Of Experimental And Observational Climate Data In Function Of Longitudes And Latitudes}
+\description{
+This sample data set contains gridded seasonal forecast and corresponding observational data from the Copernicus Climate Change ECMWF-System 5 forecast system, and from the Copernicus Climate Change ERA-5 reconstruction. Specifically, for the 'tas' (2-meter temperature) variable, for the 15 first forecast ensemble members, monthly averaged, for the 3 first forecast time steps (lead months 1 to 4) of the November start dates of 2000 to 2005, for the Mediterranean region (27N-48N, 12W-40E). The data was generated on (or interpolated onto, for the reconstruction) a rectangular regular grid of size 360 by 181.
+}
+\details{
+It is recommended to use the data set as follows:
+```
+require(zeallot)
+c(exp, obs) %<-% CSTools::lonlat_data
+```
+
+The `CST_Load` call used to generate the data set in the infrastructure of the Earth Sciences Department of the Barcelona Supercomputing Center is shown next. Note that `CST_Load` internally calls `s2dverification::Load`, which would require a configuration file (not provided here) expressing the distribution of the 'system5c3s' and 'era5' NetCDF files in the file system.
+```
+library(CSTools)
+require(zeallot)
+
+startDates <- c('20001101', '20011101', '20021101',
+                '20031101', '20041101', '20051101')
+
+c(exp, obs) %<-% 
+  CST_Load(
+    var = 'tas', 
+    exp = 'system5c3s', 
+    obs = 'era5', 
+    nmember = 15,
+    sdates = startDates,
+    leadtimemax = 3,
+    latmin = 27, latmax = 48,
+    lonmin = -12, lonmax = 40, 
+    output = 'lonlat',
+    nprocs = 1
+  )
+```
+}
+\author{
+Nicolau Manubens \email{nicolau.manubens@bsc.es}
+}
+\keyword{data}
diff --git a/tests/testthat/test-CST_BiasCorrection.R b/tests/testthat/test-CST_BiasCorrection.R
index acc7425d51c29fa2958145400f223e6520ff4c18..39d5daa510a2c85e46daf49b7758bd60b5d055de 100644
--- a/tests/testthat/test-CST_BiasCorrection.R
+++ b/tests/testthat/test-CST_BiasCorrection.R
@@ -1,44 +1,44 @@
 context("Generic tests")
 test_that("Sanity checks", {
   expect_error(
-    CST_BiasCorrection(data = 1),
-    "Parameter 'data' must be a list as output of Load function
-    from s2dverification package")
+    CST_BiasCorrection(exp = 1),
+    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
+         "as output by CSTools::CST_Load."))
   
   mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
   obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
-  dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
-  dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
+  dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, 
+                 lat = 6, lon = 7)
+  dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, 
+                 lat = 6, lon = 7)
   lon <- seq(0, 30, 5)
   lat <- seq(0, 25, 5)
-  data <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
+  exp <- list(data = mod1, lat = lat, lon = lon)
+  obs <- list(data = obs1, lat = lat, lon = lon)
+  attr(exp, 'class') <- 's2dv_cube'
+  attr(obs, 'class') <- 's2dv_cube'
   
-  expect_equal(length(CST_BiasCorrection(data = data)), 3)
-  expect_equal(dim(CST_BiasCorrection(data = data)$biascorrection), 
-               c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7))
-  expect_equal(CST_BiasCorrection(data = data)$lat, lat)
-  expect_equal(CST_BiasCorrection(data = data)$lon, lon)
+  expect_equal(length(CST_BiasCorrection(exp = exp, obs = obs)), 3)
+  expect_equal(dim(CST_BiasCorrection(exp = exp, obs = obs)$data), 
+               c(dataset = 1, member = 3, sdate = 4, ftime = 5, 
+                 lat = 6, lon = 7))
+  expect_equal(CST_BiasCorrection(exp = exp, obs = obs)$lat, lat)
+  expect_equal(CST_BiasCorrection(exp = exp, obs = obs)$lon, lon)
   
-  data <- list(mod = mod1, obs = mod1, lat = lat, lon = lon)
-  expect_error(CST_BiasCorrection(data = data),
-               "The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.")
+  expect_error(CST_BiasCorrection(exp = exp, obs = exp),
+        paste0("The length of the dimension 'member' in the component 'data' ",
+         "of the parameter 'obs' must be equal to 1."))
   
-  mod2 <- mod1
-  mod2[1, 2, 1, 1, 1, 1] <- NA
-  data <- list(mod = mod2, obs = obs1, lat = lat, lon = lon)
-  expect_warning(
-    CST_BiasCorrection(data = data),
+  exp2 <- exp
+  exp2$data[1, 2, 1, 1, 1, 1] <- NA
+  expect_warning(CST_BiasCorrection(exp = exp2, obs = obs),
     "Parameter 'exp' contains NA values.")
   
-  obs2 <- obs1
-  obs2[1, 1, 2, 1, 1, 1] <- NA
-  data <- list(mod = mod1, obs = obs2, lat = lat, lon = lon)
-  expect_warning(
-    CST_BiasCorrection(data = data),
+  obs2 <- obs
+  obs2$data[1, 1, 2, 1, 1, 1] <- NA
+  expect_warning(CST_BiasCorrection(exp = exp, obs = obs2),
     "Parameter 'obs' contains NA values.")
   
-  data <- list( mod = mod2, obs = obs2, lat = lat, lon = lon)
-  expect_warning(
-    CST_BiasCorrection(data = data),
+  expect_warning(CST_BiasCorrection(exp = exp2, obs = obs2),
     "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.")
-})
\ No newline at end of file
+})
diff --git a/tests/testthat/test-CST_Calibration.R b/tests/testthat/test-CST_Calibration.R
index 382bc6345da068a2025461ab704ab27708ba736f..23cd4c21311c06e83b27c79c48a3856c0d8b5913 100644
--- a/tests/testthat/test-CST_Calibration.R
+++ b/tests/testthat/test-CST_Calibration.R
@@ -1,44 +1,35 @@
 context("Generic tests")
 test_that("Sanity checks", {
   expect_error(
-    CST_Calibration(data = 1),
-    "Parameter 'data' must be a list as output of Load function
-from s2dverification package")
-  
-  mod1 <- 1 : (1 * 3 * 4 * 5 * 6 * 7)
-  obs1 <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
-  dim(mod1) <- c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
-  dim(obs1) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
-  lon <- seq(0, 30, 5)
-  lat <- seq(0, 25, 5)
-  data <- list(mod = mod1, obs = obs1, lat = lat, lon = lon)
-  
-  expect_equal(length(CST_Calibration(data = data)), 3)
-  expect_equal(dim(CST_Calibration(data = data)$calibration), 
-            c(dataset = 1, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7))
-  expect_equal(CST_Calibration(data = data)$lat, lat)
-  expect_equal(CST_Calibration(data = data)$lon, lon)
-  
-  data <- list(mod = mod1, obs = mod1, lat = lat, lon = lon)
-  expect_error(CST_Calibration(data = data),
-  "The length of the dimension 'member' in label 'obs' of parameter 'data' must be equal to 1.")
+    CST_Calibration(exp = 1),
+    c("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
+         "as output by CSTools::CST_Load."))
+  expect_error(
+    CST_Calibration(obs = 1),
+    c("argument \"exp\" is missing, with no default"))
+ 
+  library(zeallot)  
+  c(exp, obs) %<-% lonlat_data
+  expect_equal(length(CST_Calibration(exp = exp, obs = obs)), 9)
+  expect_equal(dim(CST_Calibration(exp = exp, obs = obs)$data), dim(exp$data))
+  expect_equal(CST_Calibration(exp = exp, obs = obs)$lat, exp$lat)
+  expect_equal(CST_Calibration(exp = exp, obs = obs)$lat, obs$lat)
+  expect_equal(CST_Calibration(exp = exp, obs = obs)$lon, exp$lon)
+  expect_equal(CST_Calibration(exp = exp, obs = obs)$lon, obs$lon)
+  expect_error(CST_Calibration(exp = exp, obs = exp),
+  c("The length of the dimension 'member' in the component 'data' ",
+         "of the parameter 'obs' must be equal to 1."))
   
-  mod2 <- mod1
-  mod2[1, 2, 1, 1, 1, 1] <- NA
-  data <- list(mod = mod2, obs = obs1, lat = lat, lon = lon)
-  expect_warning(
-    CST_Calibration(data = data),
+  exp2 <- exp
+  exp2$data[1, 2, 1, 1, 1, 1] <- NA
+  expect_warning(CST_Calibration(exp = exp2, obs = obs),
     "Parameter 'exp' contains NA values.")
   
-  obs2 <- obs1
-  obs2[1, 1, 2, 1, 1, 1] <- NA
-  data <- list(mod = mod1, obs = obs2, lat = lat, lon = lon)
-  expect_warning(
-    CST_Calibration(data = data),
+  obs2 <- obs
+  obs2$data[1, 1, 2, 1, 1, 1] <- NA
+  expect_warning(CST_Calibration(exp = exp, obs = obs2),
     "Parameter 'obs' contains NA values.")
   
-  data <- list( mod = mod2, obs = obs2, lat = lat, lon = lon)
-  expect_warning(
-    CST_Calibration(data = data),
+  expect_warning(CST_Calibration(exp = exp2, obs = obs2),
     "Parameter 'obs' contains NA values", "Parameter 'exp' contains NA values.")
-})  
\ No newline at end of file
+})  
diff --git a/tests/testthat/test-CST_MultiMetric.R b/tests/testthat/test-CST_MultiMetric.R
index bdcc4d0f41ddf371eed3662246d946a68cdae386..b08b68ba9f941c2cb195b16d3511b3b995bd990f 100644
--- a/tests/testthat/test-CST_MultiMetric.R
+++ b/tests/testthat/test-CST_MultiMetric.R
@@ -6,73 +6,68 @@ test_that("basic use case", {
   dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 8)
   lon <- seq(0, 30, 5)
   lat <- seq(0, 30, 5)
-  dat <- list(mod = mod, obs = obs, lat = lat, lon = lon)
+  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'
   
-  result <- list(ano_exp = array(rep(
-                c(rep(-12, 6), rep(-4, 6), rep(4, 6), rep(12, 6)), 360),
-                dim = c(dataset = 2, member = 3, sdate = 4, ftime = 5, 
-                lat = 6, lon = 8)),
-                ano_obs = array(rep(c(-2, -2/3,  2/3,  2), 360),
-                dim = c(dataset = 1, member = 1, sdate = 4, ftime = 5,
-                lat = 6, lon = 8)),
-                lat = lat,
-                lon = lon,
-                metric = array(rep(c(rep(1, 9), rep(0.9, 3)), 48), 
+  result <- list(data = array(rep(c(rep(1, 9), 
+rep(0.89999999999999991118215802998747676610950, 3)), 48), 
                 dim = c(dataset = 3, dataset = 1, statistics = 4, 
-                       lat = 6, lon = 8))) 
-  attr(result[[1]], "experiments") <- c("mod1", "mod2", "MMM")
-  expect_equal(CST_MultiMetric(dat), result)
+                       lat = 6, lon = 8)),
+                lat = lat, lon = lon) 
+  attr(result, 'class') <- 's2dv_cube'
+  expect_equal(CST_MultiMetric(exp = exp, obs = obs), result)
   
-  res <- CST_MultiMetric(dat, metric = 'rms')
-  expect_equal(length(res), 5)
-  expect_equal(dim(res$ano_exp), 
-        c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 8))
-  expect_equal(dim(res$ano_obs), 
-        c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 8))
-  expect_equal(dim(res$metric), 
+  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))
-  res <- CST_MultiMetric(dat, metric = 'rms', multimodel = FALSE)
-  expect_equal(dim(res$metric), 
+  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))
-  res <- CST_MultiMetric(dat, metric = 'rmsss')
-  expect_equal(dim(res$metric), 
+  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))
-  res <- CST_MultiMetric(dat, metric = 'rmsss', multimodel = FALSE)
-  expect_equal(dim(res$metric), 
+  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))
-  res <- CST_MultiMetric(dat, names = c('A', 'B'))
-  expect_equal(attributes(res$ano_exp)$experiments, c('A', 'B', 'MMM'))
-})
+  })
 
 
 test_that("Sanity checks", {
   expect_error(
-    CST_MultiMetric(data = 1),
-      "Parameter 'data' must be a list as output of Load function
-from s2dverification package")
+    CST_MultiMetric(exp = 1),
+    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
+         "as output by CSTools::CST_Load."))
   mod <- 1 : (2 * 3 * 4 * 5 * 6 * 8)
   dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 8)
   obs <- 1 : (1 * 1 * 4 * 5 * 6 * 8)
   dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 8)
   lon <- seq(0, 30, 5)
   lat <- seq(0, 30, 5)
-  dat <- list(mod = mod, obs = obs, lat = lat, lon = lon)
+  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'
+
   
   expect_error(
-    CST_MultiMetric(data = dat, metric = 1),
-    "Parameter 'metric' must be a character string indicating one of the 
-         options: 'correlation', 'rms' or 'rmse'")
+    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'"))
   expect_error(
-    CST_MultiMetric(data = dat, metric = NA),
+    CST_MultiMetric(exp = exp, obs = obs, metric = NA),
     "missing value where TRUE/FALSE needed")
   expect_error(
-    CST_MultiMetric(data = dat, metric = NULL),
+    CST_MultiMetric(exp = exp, obs = obs, metric = NULL),
     "argument is of length zero")
   expect_error(
-    CST_MultiMetric(data = dat, metric = "correlation", multimodel = NULL),
+    CST_MultiMetric(exp = exp, obs = obs, metric = "correlation", 
+                    multimodel = NULL),
     "Parameter 'multimodel' must be a logical value.")
-  expect_error(
-    CST_MultiMetric(data = dat, metric = "correlation", multimodel = FALSE, names = 1),
-    "Parameter 'names' must be the same length as the number of 
-           models in the element 'mod' of parameter 'data'.")
 })  
diff --git a/vignettes/Figures/RainFARM_fig1.png b/vignettes/Figures/RainFARM_fig1.png
new file mode 100644
index 0000000000000000000000000000000000000000..1818c6e3e2ac0f264ad07b8b73e987492a5d3a7f
Binary files /dev/null and b/vignettes/Figures/RainFARM_fig1.png differ
diff --git a/vignettes/Figures/RainFARM_fig2.png b/vignettes/Figures/RainFARM_fig2.png
new file mode 100644
index 0000000000000000000000000000000000000000..a821f4ad0384de5de83b6475d6da76210568ea1c
Binary files /dev/null and b/vignettes/Figures/RainFARM_fig2.png differ
diff --git a/vignettes/MultiModelSkill_vignette.md b/vignettes/MultiModelSkill_vignette.md
index bd658ea0b5b6efbdd1c9a2e3d7aff280a889bada..6fc676109908add0004f31bdcfe2ab141b33cd4a 100644
--- a/vignettes/MultiModelSkill_vignette.md
+++ b/vignettes/MultiModelSkill_vignette.md
@@ -1,3 +1,13 @@
+---
+author: "Nuria Perez"
+date: "`r Sys.Date()`"
+output: html_document
+vignette: >
+  %\VignetteEngine{R.rsp::md}
+  %\VignetteIndexEntry{Multi-model Skill Assessment}
+  %\usepackage[utf8]{inputenc}
+---
+
 Multi-model Skill Assessment
 -----------------------------------------
 
@@ -11,6 +21,7 @@ A version of s2dverification under development should be loaded by running:
 devtools::install_git('https://earth.bsc.es/gitlab/es/s2dverification', 
                        branch = 'mostlikely')
 library(s2dverification)
+library(zeallot)
 ```
 
 
@@ -56,78 +67,99 @@ grid <- "256x128"
 obs <- "erainterim"
 ```
 
-Using the `Load` function from **s2dverification package**, the data available in our data store can be loaded. The following lines, shows how this function can be used. However, the data is loaded from a previous saved `.RData` file:
+Using the `CST_Load` function, the data available in our data store can be loaded. The following lines, shows how this function can be used. However, the data is loaded from a previous saved `.RData` file:
 Ask nuria.perez at bsc.es to achieve 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')
-#StartData <- Load(clim_var, exp = list(glosea5, list(name = 'ecmwf/system4_m1'),
-#                                  list(name = 'meteofrance/system5_m1')),
-#                  obs = obs, sdates = dateseq, leadtimemin = 2, leadtimemax = 4,
-#                  lonmin = -15, lonmax = 45, latmin = 25, latmax = 50, 
-#                  storefreq = "monthly", sampleperiod = 1,  nmember = 9, 
-#                  output = "lonlat", method = "bilinear", 
-#                  grid = paste("r", grid, sep = ""))
-#save(StartData, StartData, file = "tas_toydata.RData")         
-#load(file = "./tas_toydata.RData")
+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'),
+                                  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")   
+
+# Or use the following line to load the file provided in .RData format:
+load(file = "./tas_toydata.RData")
 ```
 
+There should be two new elements loaded in the R working environment: `exp` and `obs`, containing the experimental and the observed data for temperature. It's possible to check that they are of class `sd2v_cube` by running:
+
 
-In this example the first two elements of `StartData` corresponds to arrays containing the simulated and the observed data:
+```
+class(exp)
+class(obs)
+```
 
+The corresponding data is saved in the element `data` of each object, while other relevant information is saved in different elements, such as `lat` and `lon`:
 
 ```r
-> dim(StartData$mod)
+> dim(exp$data)
 dataset  member   sdate   ftime     lat     lon 
       3       9      21       3      35      64 
-> dim(StartData$obs)
+> dim(obs$data)
 dataset  member   sdate   ftime     lat     lon 
       1       1      21       3      35      64 
-```      
+Lat <- exp$lat
+Lon <- exp$lon
+```
 
 
-Latitudes and longitudes of the common grid can be saved:
+### 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). 
+
+
+First step is to compute the anomalies over the loaded data applying cross validation technique on individual members by running:
 
-```r
-Lat <- StartData$lat
-Lon <- StartData$lon
+```
+c(ano_exp, ano_obs) %<-% CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
 ```
 
+The dimensions are preserved:
 
-### 2.- Computing and plotting Anomaly Correlation Coefficient
+```
+> 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 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 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(data = StartData, metric = 'correlation', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'correlation', 
                           multimodel = TRUE)
 ```
 
 
-The first element of the list returned by the function `CST_MultiMetric` contains the anomaly data for the experiments and the observations, the element called 'metric' is the correlations (including the correlation for the MMM in the latest position).
-
-Here, it's a summary of the output:
+The output of the function `CST_MultiMetric` is a object of class `s2dv_cube`, it contains the result of the metric, in this case correlation, in the `data` element (including the correlation for the MMM in the latest position).
+While other relevant data is being stored in the corresponding element of the object:
 
 
 ```r
-> str(AnomDJF)
-List of 12
- $ ano_exp        : num [1:3, 1:9, 1:21, 1:3, 1:35, 1:64] -1.647 -0.478 -0.096 1.575 1.086 ...
-  ..- attr(*, "experiments")= chr [1:4] "exp1" "ecmwf/system4_m1" "meteofrance/system5_m1" "MMM"
- $ ano_obs        : num [1, 1, 1:21, 1:3, 1:35, 1:64] 0.0235 1.546 1.3885 -0.344 -5.972 ...
- $ Datasets       :List of 2
-(...)
- $ metric         : num [1:4, 1, 1:3, 1:35, 1:64] 1.02 1.35 1.17 1.05 1.3 ...
-(...)
+> dim(AnomDJF$data)
+   dataset    dataset statistics        lat        lon 
+         4          1          4         35         43 
+> names(AnomDJF)
+[1] "data"            "lon"             "lat"             "Variable"        "Datasets"        "Dates"          
+[7] "when"            "source_files"    "load_parameters"
+> names(AnomDJF$Datasets)
+[1] "glosea5"                "ecmwf/system4_m1"       "meteofrance/system5_m1" "erainterim"    
 ```
  
 In the second element of the list `AnomDJF`, `metric`, 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$metric[ , , 2, , ]
+corre <- AnomDJF$data[ , , 2, , ]
 names(dim(corre)) <- c("maps", "lat", "lon")
 ```
 
@@ -142,7 +174,7 @@ PlotCombinedMap(corre, lon = Lon, lat = Lat, map_select_fun = max,
                             c('white', 'darkred'),
                             c('white', 'darkorange'),
                             c('white', 'black')),
-                bar_titles = c('glosea5', attributes(AnomDJF$ano_exp)$experiments[-1]),
+                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"),
                 fileout = "./vignettes/Figures/MultiModelSkill_cor_tas_1992-2012.png",
                 width = 14, height = 8)                
 ```
@@ -159,9 +191,9 @@ 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(data = StartData, metric = 'rms', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rms', 
                           multimodel = TRUE)
-RMS <- AnomDJF$metric[ , , 2, , ]
+RMS <- AnomDJF$data[ , , 2, , ]
 ```
 
 The following lines are necessary to obtain the plot which visualizes the best model given this metric for each grid point. 
@@ -175,7 +207,7 @@ PlotCombinedMap(RMS, lon = Lon, lat = Lat, map_select_fun = min,
                             c('darkred', 'white'),
                             c('darkorange', 'white'),
                             c('black', 'white')),
-                bar_titles = c('glosea5', attributes(AnomDJF$ano_exp)$experiments[-1]), 
+                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), 
                 fileout = "./vignettes/Figures/MultiModelSkill_rms_tas_1992-2012.png",
                 width = 14, height = 8)              
 ```
@@ -192,9 +224,9 @@ Notice that the perfect RMSSS is 1 and the parameter  `map_select_fun` from `Plo
 
 
 ```r
-AnomDJF <- CST_MultiMetric(data = StartData, metric = 'rmsss', 
+AnomDJF <- CST_MultiMetric(exp = ano_exp, obs = ano_obs, metric = 'rmsss', 
                           multimodel = TRUE)
-RMSSS <- AnomDJF$metric[ , , 1, , ]
+RMSSS <- AnomDJF$data[ , , 1, , ]
 names(dim(RMSSS)) <- c("maps", "lat", "lon")
 
 PlotCombinedMap(RMSSS, lon = Lon, lat = Lat, 
@@ -206,7 +238,7 @@ PlotCombinedMap(RMSSS, lon = Lon, lat = Lat,
                             c('white', 'darkred'),
                             c('white', 'darkorange'),
                             c('white', 'black')),
-                bar_titles = c('glosea5', attributes(AnomDJF$ano_exp)$experiments[-1]), 
+                bar_titles = c(names(AnomDJF$Datasets)[-4], "MMM"), 
                 fileout = "./vignettes/Figures/MultiModelSkill_rmsss_tas_1992-2012.png",
                 width = 14, height = 8)                
 ```
diff --git a/vignettes/MultivarRMSE_vignette.md b/vignettes/MultivarRMSE_vignette.md
index ec52bc320a2551176735e4e67f3c50c11d0fcf7a..78af60807aa75815237fed804619e25f172de392 100644
--- a/vignettes/MultivarRMSE_vignette.md
+++ b/vignettes/MultivarRMSE_vignette.md
@@ -1,8 +1,13 @@
 ---
-output:
-  pdf_document: default
-  html_document: default
+author: "Deborah Verfaillie"
+date: "`r Sys.Date()`"
+output: html_document
+vignette: >
+  %\VignetteEngine{R.rsp::md}
+  %\VignetteIndexEntry{Multivariate RMSE}
+  %\usepackage[utf8]{inputenc}
 ---
+
 Multivariate Root Mean Square Error (RMSE)
 ------------------------------------------
 
@@ -12,14 +17,16 @@ To run this vignette, the next R packages should be installed and loaded:
 
 ```r
 library(s2dverification)
-library(startR)
 library(RColorBrewer)
+library(zeallot)
 ```
 
-The function hosted in the *CSTools* gitlab project should be loaded by running:
+
+Library *CSTools*, should be installed from CRAN and loaded:
+
 
 ```r
-devtools::install_git('https://earth.bsc.es/gitlab/external/cstools', branch = 'master')
+install.packages("CSTools")
 library(CSTools)
 ```
 
@@ -58,53 +65,85 @@ grid <- "256x128"
 obs <- "erainterim"
 ```
 
-Using the `Load` function from **s2dverification package**, the data available in our data store can be loaded. The following lines show how this function can be used. Here, the data is loaded from a previous saved `.RData` file:
-Ask deborah.verfaillie@bsc.es for the data to run the recipe.
+Using the `CST_Load` function from **CSTool package**, the data available in our data store can be loaded. The following lines show how this function can be used. Here, the data is loaded from a previous saved `.RData` file:
+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')
-# Temp <- Load(temp, exp = list(glosea5), 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 = ""))
-# Precip <- Load(precip, exp = list(glosea5), 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(Temp, Precip, file = "tas_prlr_toydata.RData")         
-load(file = "tas_prlr_toydata.RData")
+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 = ""))
+c(exp_P, obs_P) %<-% 
+   CST_Load(var = precip,  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 = ""))
+#save(exp_T, obs_T, exp_P, obs_P, file = "./tas_prlr_toydata.RData")         
+
+# Or use the following line to load the file provided in .RData format:
+load(file = "./tas_prlr_toydata.RData")
 ```
 
-In this example the first two elements of `Temp` and `Precip` correspond to arrays containing the simulated and the observed data:
+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:
 
 
-```r
-> dim(Temp$mod)
+```
+class(exp_T)
+class(obs_T)
+class(exp_P)
+class(obs_P)
+```
+
+Loading the data using `CST_Load` allows to obtain two lists, one for the experimental data and another for the observe data, with the same elements and compatible dimensions of the data element:
+
+
+```
+> dim(exp_T$data)
 dataset  member   sdate   ftime     lat     lon 
-      1       9      21       3      35      64 
-> dim(Temp$obs)
+      1       9      21       3      35      64
+> dim(obs_T$data)
 dataset  member   sdate   ftime     lat     lon 
-      1       1      21       3      35      64 
-```      
-
+      1       1      21       3      35      64
+```
 
 Latitudes and longitudes of the common grid can be saved:
 
 
 ```r
-Lat <- Temp$lat
-Lon <- Temp$lon
+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.
 
-A list containing the `Temp` and `Precip` 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
-StartData <- list(Temp, Precip)
+ano_exp <- list(ano_exp_T, ano_exp_P)
+ano_obs <- list(ano_obs_T, ano_obs_P)
 weight <- c(2, 1)
 ```
 
@@ -115,17 +154,20 @@ It is obtained by running the `CST_MultivarRMSE` function:
 
 
 ```r
-mvrmse <- CST_MultivarRMSE(StartData, weight)
+mvrmse <- CST_MultivarRMSE(exp = ano_exp, obs = ano_obs, weight)
 ```
 
 
-The function `CST_MultivarRMSE` returns the multivariate RMSE value for 2 or more variables. The output is a CSTool object containing a list of elements being the first one the mvrmse:
+The function `CST_MultivarRMSE` returns the multivariate RMSE value for 2 or more variables. The output is a CSTool object containing  the RMSE values in the `data` element and other relevant information:
 
 
 ```r
-> str(mvrmse)
-List of 10
-$ mvrmse: num [1, 1, 1, 1:35, 1:64] 0.764 0.8 0.67 0.662 0.615 ...
+> class(mvrmse)
+> str(mvrmse$data)
+ num [1, 1, 1, 1:35, 1:64] 0.764 0.8 0.67 0.662 0.615 ...
+> str(mvrmse$Variable)
+ Named chr [1:2] "tas" "prlr"
+ - attr(*, "names")= chr [1:2] "varName" "varName"
 ```
 
 
@@ -133,7 +175,7 @@ The following lines plot the multivariate RMSE
 
 
 ```r
-PlotEquiMap(mvrmse$mvrmse, lon = Lon, lat = Lat, filled.continents = FALSE, 
+PlotEquiMap(mvrmse$data, lon = Lon, lat = Lat, filled.continents = FALSE, 
             toptitle = "Multivariate RMSE tas, prlr 1992 - 2012", colNA = "white", 
             bar_limits = c(0,2.5), cols = brewer.pal(n=5,name='Reds'),
             fileout = "./MultivarRMSE_gloseas5_tas_prlr_1992-2012.png")
diff --git a/vignettes/RainFARM_vignette.md b/vignettes/RainFARM_vignette.md
new file mode 100644
index 0000000000000000000000000000000000000000..f8c286436caa9452ff826818339f25a360b4c337
--- /dev/null
+++ b/vignettes/RainFARM_vignette.md
@@ -0,0 +1,216 @@
+---
+title: "Rainfall Filtered Autoregressive Model (RainFARM) precipitation downscaling"
+author: "Jost von Hardenberg (ISAC-CNR)"
+date: "26/03/2019"
+output:
+  pdf_document: default
+  html_document: default
+vignette: >
+  %\VignetteEngine{R.rsp::md}
+  %\VignetteIndexEntry{RainFARM}
+  %\usepackage[utf8]{inputenc}
+---
+
+<!--
+```{r setup, include=FALSE}
+knitr::opts_chunk$set(echo = TRUE)
+```
+-->
+
+## Introduction
+
+It becomes often necessary to bridge the scale gap between climate scenarios of precipitation obtained from global and regional climate models and the small scales needed for impact studies, in order to measure uncertainty and reproduce extremes at these scales. An effective approach is provided by stochastic downscaling techniques (Ferraris et al. 2003a,b), in particular so-called full-field weather generators, which aim at generating synthetic spatiotemporal precipitation fields whose statistical properties are consistent with the small-scale statistics of observed precipitation, based only on knowledge of the large-scale precipitation field. These methods allow for example to estimate uncertainties at small scales in rainfall scenarios, by generating large ensembles of synthetic small-scale precipitation fields. 
+
+In particular the _Rainfall Filtered Autoregressive Model_ (**RainFARM**; Rebora et al. 2006a,b) 
+is a stochastic downscaling procedure based on the nonlinear transformation of a linearly correlated stochastic field, generated by small-scale extrapolation of the Fourier spectrum of a large-scale precipitation field. Developed originally for downscaling at weather timescales, the method has been adapted for downscaling at climate timescales by D'Onofrio et al. 2014 and recently improved for regions with complex orography in Terzago et al. 2018. 
+
+An efficient implementation of RainFARM is provided for CSTools by the `CST_RainFARM()` and `RainFARM()` functions. 
+
+## Downscaling seasonal precipitation forecasts with RainFARM
+
+### Preliminary setup
+
+In order to run the examples in this vignette, the *CSTools* package and some other support R packages need to be loaded by running:
+
+```r
+#install.packages('CSTools')
+library(CSTools)
+```
+
+We use the `CST_Load` function to load a seasonal precipitation forecast from the ecmwfS5 forecast. This function is forwarding all of its parameters to the `Load` data loading function in the *s2dverification* package, and returns data objects compatible with all functions in CSTools.
+
+There are 25 ensemble members available in the data set. We select a 4x4 pixel cutout at resolution 1 degree in a smaller area lon=[6,9], lat=[44,47] covering North-Western Italy and three starting dates on November 1st for years 2010, 2011 and 2012, selecting forecast times for the following month of March:
+
+```r  
+infile <- list(path = '../medscope/nwalps/data/$VAR_NAME$_$START_DATE$_nwalps.nc')
+exp <- CST_Load('prlr', exp = list(infile), obs = NULL,
+                sdates = c('20101101', '20111101', '20121101'), 
+                leadtimemin = 121, leadtimemax = 151, 
+                latmin = 44, latmax = 47, 
+                lonmin = 5, lonmax = 9, 
+                nmember = 25,
+                storefreq = "daily", sampleperiod = 1, 
+                output = "lonlat"
+               )$exp
+```
+
+This gives us a CSTools object `exp`, containing an element `exp$data` with dimensions:
+```r  
+dim(exp$data)
+dataset  member   sdate   ftime     lat     lon 
+      1      25       3      31       4       4 
+```
+
+Please notice that RainFARM (in this version) only accepts square domains, possibly with an even number of pixels on each side, so we will need to select an appropriate cutout.
+
+### Standard downscaling without climatological weights
+
+Our goal is to downscale with RainFARM these data from the resolution of 1 degree (about 100 km at these latitudes) to 0.05 degrees (about 5 km) using the `CST_RainFARM()` function. This means that we need to increase resolution by a factor `nf = 20`. 
+RainFARM can compute automatically its only free parameter, i.e. the spatial spectral slope, from the large-scale field (here only with size 4x4 pixel, but we suggest to use at least 8x8 pixel in a real application). 
+In this example we would like to compute this slope as an average over the _member_ and _ftime_ dimensions, while we will use different slopes for the remaining _dataset_ and _sdate_ dimensions (a different choice may be more appropriate in a real application). To obtain this we specify the parameter `time_dim = c("member", "ftime")`. The slope is computed starting from the wavenumber corresponding to the box, `kmin=1`. We create 5 stochastic realizations for each dataset, member, starting date and forecast time with `nens=5`. The command to donwscale and the resulting fields are:
+
+```r  
+exp_down <- CST_RainFARM(exp, nf=20, kmin = 1, nens = 5,
+                         time_dim = c("member", "ftime"))
+
+dim(exp_down$data)
+    dataset      member realization       sdate       ftime         lat        lon
+          1          25           5           3          31          80         80
+str(exp_down$lon)
+ num [1:80] 5.53 5.58 5.62 5.67 5.72 ...
+str(exp_down$lat)
+ num [1:80] 47.5 47.4 47.4 47.3 47.3 ...
+```
+The function returns an array `exp_down$data` with the additional "realization" dimension for the stochastic ensemble with 5 members. The longitudes and latitudes have been correspondingly interpolated to the finer resolution. 
+
+Alternatively we could have used the "reduced" function `RainFARM` which accepts directly a data array (with arbitrary dimensions, provided a longitude, a latitude and a "time" dimension exist) and two arrays to describe longitudes and latitudes:
+
+```r  
+downscaled <- RainFARM(exp$data, exp$lon, exp$lat, nf = 20, kmin = 1, nens = 5,
+                       time_dim = c("member", "ftime"))
+```
+The resulting array has the same dimensions as the `$data` element in the result of `CST_RainFARM`.
+
+Each instant and each realization will of course be different, but let's plot and compare the original and the downscaled data of the first ensemble member and for the first realization for a specific forecast date (Fig. 1):
+
+<!--  
+png("Figures/RainFARM_fig1.png", width = 640, height = 365)
+par(mfrow=c(1,2))
+-->
+```r  
+a <- exp$data[1, 1, 1, 17, , ] * 86400 * 1000
+a[a > 60] <- 60
+image(exp$lon, rev(exp$lat), t(apply(a, 2, rev)), xlab = "lon", ylab = "lat",
+      col = rev(terrain.colors(20)), zlim = c(0,60))
+map("world", add = TRUE)
+title(main = "pr 17/03/2010 original")
+a <- exp_down$data[1, 1, 1, 1, 17, , ] * 86400 * 1000
+a[a > 60] <- 60
+image(exp_down$lon, rev(exp_down$lat), t(apply(a, 2, rev)), xlab = "lon", ylab = "lat", 
+      col = rev(terrain.colors(20)), zlim = c(0, 60))
+map("world", add = TRUE)
+title(main = "pr 17/03/2010 downscaled")
+```
+<!-- dev.off() -->
+
+![Comparison between an original precipitation field at 1-degree resolution (left) and the same field downscaled with RainFARM to 0.05 degree resolution (right)](Figures/RainFARM_fig1.png){width=420px}
+
+RainFARM has downscaled the original field with a realistic fine-scale correlation structure. Precipitation is conserved in an average sense (in this case smoothing both the original and the downscaled fields with a circular kernel with a diameter equal to the original field grid spacing would lead to the same results). The downscaled field presents more extreme precipitation peaks.
+
+### Downscaling using climatological weights
+
+The area of interest in our example presents a complex orography, but the basic RainFARM algorithm used does not consider topographic elevation in deciding how to distribute fine-scale precipitation. A long term climatology of the downscaled fields would have a resolution comparable to that of the original coarse fields and would not resemble the fine-scale structure of an observed climatology. If an external fine-scale climatology of precipitation is available, we can use the method discussed in Terzago et al. (2018) to change the distribution of precipitation by RainFARM for each timestep, so that the long-term average is close to this reference climatology in terms of precipitation distribution (while the total precipitation amount of the original fields to downscale is preserved).
+
+Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (http://www.worldclim.org) or CHELSA (http://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://www.gdal.org).
+We will assume that a copy of the WORLDCLIM precipitation climatology at 30 arcseconds (about 1km resolution) is available in the local file `medscope.nc`. From this file we can derive suitable weights to be used with RainFARM using the `RFWeights` functions as follows:
+```r  
+ww <- RFWeights("./worldclim.nc", nf = 20, lon = data$lon, lat = data$lat)
+```
+The result is a two-dimensional weights matrix with the same `lon`and `lat` dimensions as requested.
+The weights (varying around an average value of 1) encode how to distribute differently precipitation in each stochastic realization of RainFARM.
+We call again `CST_RainFARM()`, this time passing climatological weights:
+```r  
+exp_down_weights <- CST_RainFARM(exp, nf = 20, kmin = 1, nens = 5, 
+                                 weights = ww, time_dim = c("member", "ftime"))
+```
+
+We plot again the same realization as before in Fig. 2 (left panel). The WORLDCLIM climatology data used in this case have missing values over the ocean, so that the resulting downscaled maps have valid points only over land. 
+
+From a single realization and time it is not possible to see that a more realistic precipitation has been achieved, but this becomes clear if we compare the climatological average over all 5 stochastic ensemble realizations, over all forecast times and over all forecast times. The resulting plot in panels Fig2b,c show that the RainFARM climatology downscaled without weights presents on average a very coarse structure, comparable with that of the original fields, while when using the weights a much more realistic distribution is achieved.
+
+<!--  
+```r
+exp_down1 <- exp_down$data[1, , , 1, , , ]
+exp_down_weights1 <- exp_down_weights$data[1, , , 1, , , ]
+dim(exp_down1) <- c(member = 25 * 5 * 31, lat = 80, lon = 80) 
+dim(exp_down_weights1) <- c(member = 25 * 5 * 31, lat = 80, lon = 80) 
+ad <- apply(exp_down1, c(2, 3), mean)
+adw <- apply(exp_down_weights1, c(2, 3), mean);
+
+png("Figures/RainFARM_fig2.png", width = 640, height = 243)
+par(mfrow = c(1,3))
+a <- exp_down_weights$data[1, 1, 1, 1, 17, , ] * 86400 * 1000
+a[a > 60] <- 60
+image(exp_down$lon, rev(exp_down$lat), t(apply(a, 2, rev)), xlab = "lon",
+      ylab = "lat", col = rev(terrain.colors(20)), zlim = c(0, 60))
+map("world", add = TRUE)
+title(main = "pr 17/03/2010 with weights")
+a <- ad * 86400 * 1000
+a[a > 5] <- 5
+image(exp_down$lon, rev(exp_down$lat), t(apply(a, 2, rev)), xlab = "lon",
+      ylab="lat", col = rev(terrain.colors(20)), zlim = c(0, 5))
+map("world", add = TRUE)
+title(main = "climatology no weights")
+a <- adw * 86400 * 1000
+a[a > 5] <- 5
+image(exp_down$lon, rev(exp_down$lat), t(apply(a, 2, rev)), xlab = "lon",
+      ylab = "lat", col = rev(terrain.colors(20)), zlim = c(0, 5))
+map("world", add = TRUE)
+title(main = "climatology with weights")
+dev.off()
+```
+-->
+
+![The same field as in Fig. 1 downscaled using climatological weights from the WORLDCLIM dataset. (left panel). The center and right panel show the climatology of the downscaled fields over all ensemble members, realizations and forecast times without and with climatological weights respectively. The center and right panels use a different colorscale than the left panel.](Figures/RainFARM_fig2.png){width=640px}
+
+### Determining the spectral slopes
+
+The only free parameter of the RainFARM method is represented by the spatial spectral slope used to extrapolate the Fourier spectrum to the unresolved scales. When fine-scale precipitation data for the study area are available this slope could be determined offline from a study of these data. Else this slope can be determined directly from the large-scale fields to downscale. The RainFARM functions compute automatically this parameter by default over the input dimensions specified by the `time_dim`parameter. In many practical applications a careful study of an appropriate average to be used may still be necessary. This can be achieved using the function `CST_RFSlope()`. If for example we wish to determine the average slopes over all member and forecast times, but computing separate slopes for each starting date and dataset, we can do so with:
+```r  
+slopes <- CST_RFSlope(exp, time_dim = c("member", "ftime"))
+dim(slopes)
+    dataset   sdate 
+          1       3   
+slopes
+         [,1]     [,2]     [,3]
+[1,] 1.532351 1.664028 1.459252
+```
+which return an array of spectral slopes, one for each "dataset" and starting date "sdate".
+
+### Compacting dimensions
+
+RainFARM creates an additional dimension in the `$data` array of the output data structure, named "realization". This is an additional ensemble dimension used to accomodate multiple stochastic realizations (if `nens > 1` is chosen), in addition to the "member" dimension of the input array (which instead is usually used to represent for example multiple forecast realizations with the same model from the same starting date). Having a separate ensemble dimension for the stochastic ensemble allows to evaluate independently uncertainty due to the small scales resolved with downscaling, compared to other sources of uncertainty such as the dependence to initial conditions gauged for example by the original "member" dimension. Since other CSTools functions support only one ensemble dimension named "member", it may be useful in some cases to compact together these two ensemble dimensions. In the output array the "realization" dimension is placed immediately after the "member" dimension, so the user could do this after calling the function reshaping the dimensions of the array. It is also possible to use the option `drop_realization_dim` of the `RainFARM` and `CST_RainFARM` functions which has the following behaviour when set to TRUE:
+
+1) if `nens == 1`: the "realization" dimension is dropped;
+
+2) if `nens > 1` and a "member" dimension exists: the
+          "realization" and "member" dimensions are compacted and the
+          resulting dimension is named "member";
+
+3) if `nens > 1` and a "member" dimension does not exist: the
+          "realization" dimension is renamed to "member".
+
+## Technical issues
+
+For maximum efficiency, the `CST_RainFARM()`, `RainFARM()`, `RFWeights()` and `RFSlope()` functions are based on an external library  written in the Julia language (https://github.com/jhardenberg/RainFARM.jl). The `JuliaCall` R package (https://CRAN.R-project.org/package=JuliaCall) is used to call the external library. 
+In particular the function uses the `julia_setup()` function from the JuliaCall library, which takes a long time (up to a few tens of seconds) only the very first time it is called. Further calls to these functions will present no delay.
+
+## Bibliography
+
+* Terzago, S., Palazzi, E., von Hardenberg, J. Stochastic downscaling of precipitation in complex orography: A simple method to reproduce a realistic fine-scale climatology. Natural Hazards and Earth System Sciences, 18(11), 2825–2840, doi:10.5194/nhess-18-2825-2018 (2018) 
+* D'Onofrio, D.; Palazzi, E., von Hardenberg, J., Provenzale A., Calmanti S. Stochastic Rainfall Downscaling of Climate Models. J of Hydrometeorology 15 (2), 830-843, doi:10.1175/JHM-D-13-096.1 (2014)
+* Rebora, N., L. Ferraris, J. von Hardenberg, A. Provenzale. Rainfall downscaling and flood forecasting: A case study in the Mediterranean area. Nat. Hazards Earth Syst. Sci., 6, 611–619, doi:10.5194/nhess-6-611-2006 (2006a)
+* Rebora, N., Ferraris, L., von Hardenberg, J., Provenzale, A. RainFARM: Rainfall downscaling by a filtered autoregressive model. J. Hydrometeor., 7, 724–738, doi:10.1175/JHM517.1 (2006b)
+* Ferraris, L., Gabellani, S., Parodi, U., Rebora, N., von Hardenberg, J., Provenzale, A. Revisiting multifractality in rainfall fields. J. Hydrometeor., 4, 544–551, doi:10.1175/ 1525-7541(2003)004,0544:RMIRF.2.0.CO;2 (2003a)
+* Ferraris, L., Gabellani, S., Rebora, N., and Provenzale, A.. A comparison of stochastic models for spatial rainfall downscaling, Water Resour. Res., 39, 1368, https://doi.org/10.1029/2003WR002504, (2003b)
+