diff --git a/DESCRIPTION b/DESCRIPTION
index 2cb6f17e9a939c77e5f74706db509c3b78098a49..9123c589b1411638ea21e364a1592dc7be137f30 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: CSTools
 Title: Assessing Skill of Climate Forecasts on Seasonal-to-Decadal
     Timescales
-Version: 4.1.0
+Version: 4.1.1
 Authors@R: c(
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-8568-3071")),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "aut", comment = c(ORCID = "0000-0001-5221-0147")),
diff --git a/NEWS.md b/NEWS.md
index b1403786ad0eb1193f366ce3f5dbd61114ebcc74..72362e9bbbc6092b81aca5a1a95d0ce1428ffffd 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,10 @@
+### CSTools 4.1.1
+**Submission date to CRAN: 10-11-2022**
+- Fixes:
+    + CST_Analogs corrected input of ClimProjDiags::Subset()
+    + PlotCombinedMap corrected use of 'cex_bar_titles' parameter
+    + CST_Anomaly added 'memb_dim', 'dat_dim' and 'ftime_dim' and improved use for 'dim_anom' parameters
+
 ### CSTools 4.1.0
 **Submission date to CRAN: 25-10-2022**
 - New features:
diff --git a/R/CST_Analogs.R b/R/CST_Analogs.R
index 28e12c20771da4a46f50e5f1e5688ab5f5cbca38..6437868091a5b8c9e4f9731dc32953ffbcee159a 100644
--- a/R/CST_Analogs.R
+++ b/R/CST_Analogs.R
@@ -36,83 +36,82 @@
 #' from surface pressure using analogues.  Clim. Dyn., 41, 1419-1437. 
 #' \email{pascal.yiou@lsce.ipsl.fr}
 #'
-#'@param expL an 's2dv_cube' object containing the experimental field on the 
-#'large scale for which the analog is aimed. This field is used to in all the
-#'criterias. If parameter 'expVar' is not provided, the function will return 
-#'the expL analog. The element 'data' in the 's2dv_cube' object must have, at
-#'least, latitudinal and longitudinal dimensions. The object is expect to be 
-#'already subset for the desired large scale region.
-#'@param obsL an 's2dv_cube' object containing the observational field on the
-#'large scale. The element 'data' in the 's2dv_cube' object must have the same 
-#'latitudinal and longitudinal dimensions as parameter 'expL' and a temporal 
-#'dimension with the maximum number of available observations.
-#'@param expVar an 's2dv_cube' object containing the experimental field on the
-#'local scale, usually a different variable to the parameter 'expL'. If it is 
-#'not NULL (by default, NULL), the returned field by this function will be the 
-#'analog of parameter 'expVar'.
-#'@param obsVar an 's2dv_cube' containing the field of the same variable as the 
-#'passed in parameter 'expVar' for the same region.
-#'@param region a vector of length four indicating the minimum longitude, the 
-#'maximum longitude, the minimum latitude and the maximum latitude.
-#'@param criteria a character string indicating the criteria to be used for the
-#'selection of analogs:
-#'\itemize{
-#'\item{Large_dist} minimum Euclidean distance in the large scale pattern; 
-#'\item{Local_dist} minimum Euclidean distance in the large scale pattern 
-#'and minimum Euclidean distance in the local scale pattern; and
-#'\item{Local_cor} minimum Euclidean distance in the large scale pattern, 
-#'minimum Euclidean distance in the local scale pattern and highest 
-#'correlation in the local variable to downscale.}
-#'Criteria 'Large_dist' is recommended for CST_Analogs, for an advanced use of 
-#'the criterias 'Local_dist' and 'Local_cor' use 'Analogs' function.
-#'@param excludeTime an array of N named dimensions (coinciding with time 
-#'dimensions in expL)of character string(s) indicating the date(s) of the 
-#'observations in the format "dd/mm/yyyy" to be excluded during the search of 
-#'analogs. It can be NULL but if expL is not a forecast (time_expL contained in
-#'time_obsL), by default time_expL will be removed during the search of analogs.
-#'@param time_expL a character string indicating the date of the experiment 
-#'in the same format than time_obsL (i.e. "yyyy-mm-dd"). By default it is NULL 
-#'and dates are taken from element \code{$Dates$start} from expL.
-#'@param time_obsL a character string indicating the date of the observations 
-#'in the date format (i.e. "yyyy-mm-dd"). By default it is NULL and dates are 
-#'taken from element \code{$Dates$start} from obsL.
-#'@param region a vector of length four indicating the minimum longitude, 
-#'the maximum longitude, the minimum latitude and the maximum latitude.
-#'@param nAnalogs number of Analogs to be selected to apply the criterias 
-#''Local_dist' or 'Local_cor'. This is not the necessary the number of analogs 
-#'that the user can get, but the number of events with minimum distance in 
-#'which perform the search of the best Analog. The default value for the 
-#''Large_dist' criteria is 1, for 'Local_dist' and 'Local_cor' criterias must
-#' be greater than 1 in order to match with the first criteria, if nAnalogs is
-#' NULL for 'Local_dist' and 'Local_cor' the default value will be set at the 
-#' length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
-#' the best analog.
-#'@param AnalogsInfo TRUE to get a list with two elements: 1) the downscaled 
-#'field and 2) the AnalogsInfo which contains: a) the number of the best 
-#'analogs, b) the corresponding value of the metric used in the selected 
-#'criteria (distance values for Large_dist and Local_dist,correlation values 
-#'for Local_cor), c)dates of the analogs). The analogs are listed in decreasing
-#'order, the first one is the best analog (i.e if the selected criteria is 
-#'Local_cor the best analog will be the one with highest correlation, while for
-#'Large_dist criteria the best analog will be the day with minimum Euclidean 
-#'distance). Set to FALSE to get a single analog, the best analog, for instance
-#'for downscaling. 
+#'@param expL An 's2dv_cube' object containing the experimental field on the 
+#'  large scale for which the analog is aimed. This field is used to in all the
+#'  criterias. If parameter 'expVar' is not provided, the function will return 
+#'  the expL analog. The element 'data' in the 's2dv_cube' object must have, at
+#'  least, latitudinal and longitudinal dimensions. The object is expect to be 
+#'  already subset for the desired large scale region.
+#'@param obsL An 's2dv_cube' object containing the observational field on the
+#'  large scale. The element 'data' in the 's2dv_cube' object must have the same 
+#'  latitudinal and longitudinal dimensions as parameter 'expL' and a temporal 
+#'  dimension with the maximum number of available observations.
+#'@param expVar An 's2dv_cube' object containing the experimental field on the
+#'  local scale, usually a different variable to the parameter 'expL'. If it is 
+#'  not NULL (by default, NULL), the returned field by this function will be the 
+#'  analog of parameter 'expVar'.
+#'@param obsVar An 's2dv_cube' containing the field of the same variable as the 
+#'  passed in parameter 'expVar' for the same region.
+#'@param region A vector of length four indicating the minimum longitude, the 
+#'  maximum longitude, the minimum latitude and the maximum latitude.
+#'@param criteria A character string indicating the criteria to be used for the
+#'  selection of analogs:
+#'  \itemize{
+#'  \item{Large_dist} minimum Euclidean distance in the large scale pattern; 
+#'  \item{Local_dist} minimum Euclidean distance in the large scale pattern 
+#'  and minimum Euclidean distance in the local scale pattern; and
+#'  \item{Local_cor} minimum Euclidean distance in the large scale pattern, 
+#'  minimum Euclidean distance in the local scale pattern and highest 
+#'  correlation in the local variable to downscale.}
+#'  Criteria 'Large_dist' is recommended for CST_Analogs, for an advanced use of 
+#'  the criterias 'Local_dist' and 'Local_cor' use 'Analogs' function.
+#'@param excludeTime An array of N named dimensions (coinciding with time 
+#'  dimensions in expL)of character string(s) indicating the date(s) of the 
+#'  observations in the format "dd/mm/yyyy" to be excluded during the search of 
+#'  analogs. It can be NULL but if expL is not a forecast (time_expL contained in
+#'  time_obsL), by default time_expL will be removed during the search of analogs.
+#'@param time_expL A character string indicating the date of the experiment 
+#'  in the same format than time_obsL (i.e. "yyyy-mm-dd"). By default it is NULL 
+#'  and dates are taken from element \code{$Dates$start} from expL.
+#'@param time_obsL A character string indicating the date of the observations 
+#'  in the date format (i.e. "yyyy-mm-dd"). By default it is NULL and dates are 
+#'  taken from element \code{$Dates$start} from obsL.
+#'@param region A vector of length four indicating the minimum longitude, 
+#'  the maximum longitude, the minimum latitude and the maximum latitude.
+#'@param nAnalogs Number of Analogs to be selected to apply the criterias 
+#'  'Local_dist' or 'Local_cor'. This is not the necessary the number of analogs 
+#'  that the user can get, but the number of events with minimum distance in 
+#'  which perform the search of the best Analog. The default value for the 
+#'  'Large_dist' criteria is 1, for 'Local_dist' and 'Local_cor' criterias must
+#'  be greater than 1 in order to match with the first criteria, if nAnalogs is
+#'  NULL for 'Local_dist' and 'Local_cor' the default value will be set at the 
+#'  length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
+#'  the best analog.
+#'@param AnalogsInfo A logical value. TRUE to get a list with two elements: 
+#'  1) the downscaled field and 2) the AnalogsInfo which contains: 
+#'  a) the number of the best analogs, b) the corresponding value of the metric 
+#'  used in the selected criteria (distance values for Large_dist and Local_dist, 
+#'  correlation values for Local_cor), c)dates of the analogs). The analogs are 
+#'  listed in decreasing order, the first one is the best analog (i.e if the 
+#'  selected criteria is Local_cor the best analog will be the one with highest 
+#'  correlation, while for Large_dist criteria the best analog will be the day 
+#'  with minimum Euclidean distance). Set to FALSE to get a single analog, the 
+#'  best analog, for instance for downscaling. 
 #'@param ncores The number of cores to use in parallel computation
-#'@import multiApply
-#'@import abind
-#'@importFrom ClimProjDiags SelBox Subset
 #'
 #'@seealso code{\link{CST_Load}}, \code{\link[s2dv]{Load}} and 
 #'\code{\link[s2dv]{CDORemap}}
 #'
-#'@return An 'array' object containing the dowscaled values of the best 
-#'analogs. 
+#'@return An 's2dv_cube' object containing an array with the dowscaled values of 
+#'the best analogs in element 'data'. If 'AnalogsInfo' is TRUE, 'data' is a list 
+#'with an array of the downscaled fields and the analogs information in 
+#'elements 'analogs', 'metric' and 'dates'.
 #'@examples
 #'expL <- rnorm(1:200)
-#'dim(expL) <- c(member=10,lat = 4, lon = 5)
+#'dim(expL) <- c(member = 10,lat = 4, lon = 5)
 #'obsL <- c(rnorm(1:180),expL[1,,]*1.2)
 #'dim(obsL) <- c(time = 10,lat = 4, lon = 5)
-#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-")
+#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994:2003, sep = "-")
 #'time_expL <- time_obsL[1]
 #'lon <- seq(-1,5,1.5)
 #'lat <- seq(30,35,1.5)
@@ -122,12 +121,16 @@
 #'                  Dates = list(start = time_obsL, end = time_obsL))
 #'region <- c(min(lon), max(lon), min(lat), max(lat))
 #'downscaled_field <- CST_Analogs(expL = expL, obsL = obsL, region = region)
+#' 
+#'@import multiApply
+#'@import abind
+#'@importFrom ClimProjDiags SelBox Subset
 #'@export
 CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL,
                         criteria = 'Large_dist', excludeTime = NULL,
                         time_expL = NULL, time_obsL = NULL,
                         nAnalogs = NULL, AnalogsInfo = FALSE,
-                        ncores = 1) {
+                        ncores = NULL) {
   if (!inherits(expL, "s2dv_cube") || !inherits(obsL, "s2dv_cube")) {
     stop("Parameter 'expL' and 'obsL' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
@@ -237,75 +240,77 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL,
 #'and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column 
 #'from surface pressure using analogues.  Clim. Dyn., 41, 1419-1437. 
 #'\email{pascal.yiou@lsce.ipsl.fr}
-#'
-#'@param expL an array of N named dimensions containing the experimental field
-#'on the large scale for which the analog is aimed. This field is used to in 
-#'all the criterias. If parameter 'expVar' is not provided, the function will
-#'return the expL analog. The element 'data' in the 's2dv_cube' object must 
-#'have, at least, latitudinal and longitudinal dimensions. The object is expect
-#'to be already subset for the desired large scale region.
-#'@param obsL an array of N named dimensions containing the observational field 
-#'on the large scale. The element 'data' in the 's2dv_cube' object must have 
-#'the same latitudinal and longitudinal dimensions as parameter 'expL' and a 
-#' single temporal dimension with the maximum number of available observations.
-#'@param time_obsL a character string indicating the date of the observations 
-#'in the format "dd/mm/yyyy". Reference time to search for analogs. 
-#'@param time_expL an array of N named dimensions (coinciding with time 
-#'dimensions in expL) of character string(s) indicating the date(s) of the 
-#'experiment in the format "dd/mm/yyyy". Time(s) to find the analogs.
-#'@param lonL a vector containing the longitude of parameter 'expL'.
-#'@param latL a vector containing the latitude of parameter 'expL'.
-#'@param excludeTime an array of N named dimensions (coinciding with time 
-#'dimensions in expL) of character string(s) indicating the date(s) of the 
-#'observations in the format "dd/mm/yyyy" to be excluded during the search of 
-#'analogs. It can be NULL but if expL is not a forecast (time_expL contained in 
-#'time_obsL),by default time_expL will be removed during the search of analogs.
-#'@param expVar an array of N named dimensions containing the experimental 
-#'field on the local scale, usually a different variable to the parameter 
-#''expL'. If it is not NULL (by default, NULL), the returned field by this 
-#'function will be the analog of parameter 'expVar'.
-#'@param obsVar an array of N named dimensions containing the field of the 
-#'same variable as the passed in parameter 'expVar' for the same region.
-#'@param AnalogsInfo TRUE to get a list with two elements: 1) the downscaled 
-#'field and 2) the AnalogsInfo which contains: a) the number of the best 
-#'analogs, b) the corresponding value of the metric used in the selected 
-#'criteria (distance values for Large_dist and Local_dist,correlation values 
-#'for Local_cor), c)dates of the analogs). The analogs are listed in decreasing
-#'order, the first one is the best analog (i.e if the selected criteria is 
-#'Local_cor the best analog will be the one with highest correlation, while for
-#'Large_dist criteria the best analog will be the day with minimum Euclidean 
-#'distance). Set to FALSE to get a single analog, the best analog, for instance
-#'for downscaling. 
+#' 
+#'@param expL An array of N named dimensions containing the experimental field
+#'  on the large scale for which the analog is aimed. This field is used to in 
+#'  all the criterias. If parameter 'expVar' is not provided, the function will
+#'  return the expL analog. The element 'data' in the 's2dv_cube' object must 
+#'  have, at least, latitudinal and longitudinal dimensions. The object is 
+#'  expect to be already subset for the desired large scale region.
+#'@param obsL An array of N named dimensions containing the observational field 
+#'  on the large scale. The element 'data' in the 's2dv_cube' object must have 
+#'  the same latitudinal and longitudinal dimensions as parameter 'expL' and a 
+#'  single temporal dimension with the maximum number of available observations.
+#'@param time_obsL A character string indicating the date of the observations 
+#'  in the format "dd/mm/yyyy". Reference time to search for analogs. 
+#'@param time_expL An array of N named dimensions (coinciding with time 
+#'  dimensions in expL) of character string(s) indicating the date(s) of the 
+#'  experiment in the format "dd/mm/yyyy". Time(s) to find the analogs.
+#'@param lonL A vector containing the longitude of parameter 'expL'.
+#'@param latL A vector containing the latitude of parameter 'expL'.
+#'@param excludeTime An array of N named dimensions (coinciding with time 
+#'  dimensions in expL) of character string(s) indicating the date(s) of the 
+#'  observations in the format "dd/mm/yyyy" to be excluded during the search of 
+#'  analogs. It can be NULL but if expL is not a forecast (time_expL contained 
+#'  in time_obsL), by default time_expL will be removed during the search of 
+#'  analogs.
+#'@param expVar An array of N named dimensions containing the experimental 
+#'  field on the local scale, usually a different variable to the parameter 
+#'  'expL'. If it is not NULL (by default, NULL), the returned field by this 
+#'  function will be the analog of parameter 'expVar'.
+#'@param obsVar An array of N named dimensions containing the field of the 
+#'  same variable as the passed in parameter 'expVar' for the same region.
+#'@param AnalogsInfo A logical value. If it is TRUE it returns a list
+#'  with two elements: 1) the downscaled field and
+#'  2) the AnalogsInfo which contains: a) the number of the best 
+#'  analogs, b) the corresponding value of the metric used in the selected 
+#'  criteria (distance values for Large_dist and Local_dist,correlation values 
+#'  for Local_cor), c)dates of the analogs). The analogs are listed in 
+#'  decreasing order, the first one is the best analog (i.e if the selected 
+#'  criteria is Local_cor the best analog will be the one with highest 
+#'  correlation, while for Large_dist criteria the best analog will be the day 
+#'  with minimum Euclidean distance). Set to FALSE to get a single analog, the 
+#'  best analog, for instance for downscaling. 
 #'@param criteria a character string indicating the criteria to be used for the
-#'selection of analogs:
-#'\itemize{
-#'\item{Large_dist} minimum Euclidean distance in the large scale pattern; 
-#'\item{Local_dist} minimum Euclidean distance in the large scale pattern 
-#'and minimum Euclidean distance in the local scale pattern; and
-#'\item{Local_cor} minimum Euclidean distance in the large scale pattern, 
-#'minimum Euclidean distance in the local scale pattern and highest 
-#'correlation in the local variable to downscale.} 
+#' selection of analogs:
+#'  \itemize{\item{Large_dist} minimum Euclidean distance in the large scale pattern; 
+#'           \item{Local_dist} minimum Euclidean distance in the large scale pattern 
+#'            and minimum Euclidean distance in the local scale pattern; and
+#'           \item{Local_cor} minimum Euclidean distance in the large scale pattern, 
+#'            minimum Euclidean distance in the local scale pattern and highest 
+#'            correlation in the local variable to downscale.} 
 #'@param lonVar a vector containing the longitude of parameter 'expVar'.
 #'@param latVar a vector containing the latitude of parameter 'expVar'.
 #'@param region a vector of length four indicating the minimum longitude, 
-#'the maximum longitude, the minimum latitude and the maximum latitude.
+#'  the maximum longitude, the minimum latitude and the maximum latitude.
 #'@param nAnalogs number of Analogs to be selected to apply the criterias 
-#''Local_dist' or 'Local_cor'. This is not the necessary the number of analogs 
-#'that the user can get, but the number of events with minimum distance in 
-#'which perform the search of the best Analog. The default value for the 
-#''Large_dist' criteria is 1, for 'Local_dist' and 'Local_cor' criterias must
-#' be greater than 1 in order to match with the first criteria, if nAnalogs is
-#' NULL for 'Local_dist' and 'Local_cor' the default value will be set at the 
-#' length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
-#' the best analog.
+#'  'Local_dist' or 'Local_cor'. This is not the necessary the number of analogs 
+#'  that the user can get, but the number of events with minimum distance in 
+#'  which perform the search of the best Analog. The default value for the 
+#'  'Large_dist' criteria is 1, for 'Local_dist' and 'Local_cor' criterias must
+#'  be greater than 1 in order to match with the first criteria, if nAnalogs is
+#'  NULL for 'Local_dist' and 'Local_cor' the default value will be set at the 
+#'  length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
+#'  the best analog.
 #'@param ncores the number of cores to use in parallel computation.
 #'@import multiApply
 #'@import abind
 #'@importFrom ClimProjDiags SelBox Subset
 #'
-#'@return AnalogsFields, dowscaled values of the best analogs for the criteria 
-#'selected. If AnalogsInfo is set to TRUE the function also returns a 
-#'list with the dowsncaled field and the Analogs Information. 
+#'@return An array with the dowscaled values of the best analogs for the criteria 
+#'selected. If 'AnalogsInfo' is set to TRUE it returns a list with an array 
+#'of the dowsncaled field and the analogs information in elements 'analogs', 
+#''metric' and 'dates'.
 #'
 #'@examples
 #'# Example 1:Downscaling using criteria 'Large_dist' and a single variable:
@@ -327,7 +332,7 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL,
 #'obsSLP <- c(rnorm(1:1980), expSLP * 1.5)
 #'dim(obsSLP) <- c(lat = 4, lon = 5, time = 100)
 #'time_obsSLP <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-")
-#'downscale_field<- Analogs(expL = expSLP, obsL = obsSLP, time_obsSLP,
+#'downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, time_obsSLP,
 #'                          nAnalogs = 5, time_expL = "01-01-2003", 
 #'                          AnalogsInfo = TRUE, excludeTime = "01-01-2003")
 #'
@@ -341,7 +346,7 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL,
 #'
 #'# Example 5: Downscaling using criteria 'Local_dist' and 2 variables:
 #'# analogs of local scale using criteria 2
-#'region=c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33)
+#'region = c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33)
 #'Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP,
 #'                       obsVar = obs.pr, criteria = "Local_dist", 
 #'                       lonL = seq(-1, 5, 1.5),latL = seq(30, 35, 1.5), 
@@ -410,11 +415,10 @@ CST_Analogs <- function(expL, obsL, expVar = NULL, obsVar = NULL, region = NULL,
 #'@export
 Analogs <- function(expL, obsL, time_obsL, time_expL = NULL, 
                     lonL = NULL, latL = NULL, expVar = NULL,
-                    obsVar = NULL,
-                    criteria = "Large_dist",excludeTime = NULL,
-                    lonVar = NULL, latVar = NULL, region = NULL,
-                    nAnalogs = NULL, AnalogsInfo = FALSE,
-                    ncores = 1) {
+                    obsVar = NULL, criteria = "Large_dist", 
+                    excludeTime = NULL, lonVar = NULL, latVar = NULL, 
+                    region = NULL, nAnalogs = NULL, 
+                    AnalogsInfo = FALSE, ncores = NULL) {
   if (!all(c('lon', 'lat') %in% names(dim(expL)))) {
     stop("Parameter 'expL' must have the dimensions 'lat' and 'lon'.")
   }
@@ -610,53 +614,53 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
                                               names(dim(obsL)))
   if (!is.null(expVar)) {
     names(dim(expVar)) <- replace_repeat_dimnames(names(dim(expVar)),
-                                                names(dim(obsVar)))
+                                                  names(dim(obsVar)))
   }
   
   if (is.null(excludeTime)) {
     excludeTime <- vector(mode="character", length=length(time_expL))
   }
- if(length(time_expL)==length(excludeTime)){
- if (any(names(dim(expL)) %in% c('sdate_exp'))) {
-    dim(time_expL) <- c(dim(expL)['sdate_exp'], dim(expL)['time_exp'])
- } else if (any(names(dim(expL)) %in% c('sdate'))) {
-   if (any(names(dim(expL)) %in% c('time_exp'))) {
-     dim(time_expL) <- c(dim(expL)['sdate'], dim(expL)['time_exp'])
-     dim(excludeTime) <- c(dim(expL)['sdate'], dim(expL)['time_exp'])
-   } else if (any(names(dim(expL)) %in% c('time'))) {
-      dim(time_expL) <- c(dim(expL)['sdate'], dim(expL)['time'])
-      dim(excludeTime) <- c(dim(expL)['sdate'], dim(expL)['time'])
-   } else {
-      dim(time_expL) <- c(dim(expL)['sdate'])
-      dim(excludeTime) <- c(dim(expL)['sdate'])
-   } 
-  } else if (any(names(dim(expL)) %in% c('time'))) {
+  if (length(time_expL) == length(excludeTime)) {
+    if (any(names(dim(expL)) %in% c('sdate_exp'))) {
+      dim(time_expL) <- c(dim(expL)['sdate_exp'], dim(expL)['time_exp'])
+    } else if (any(names(dim(expL)) %in% c('sdate'))) {
+      if (any(names(dim(expL)) %in% c('time_exp'))) {
+        dim(time_expL) <- c(dim(expL)['sdate'], dim(expL)['time_exp'])
+        dim(excludeTime) <- c(dim(expL)['sdate'], dim(expL)['time_exp'])
+      } else if (any(names(dim(expL)) %in% c('time'))) {
+        dim(time_expL) <- c(dim(expL)['sdate'], dim(expL)['time'])
+        dim(excludeTime) <- c(dim(expL)['sdate'], dim(expL)['time'])
+      } else {
+        dim(time_expL) <- c(dim(expL)['sdate'])
+        dim(excludeTime) <- c(dim(expL)['sdate'])
+      } 
+    } else if (any(names(dim(expL)) %in% c('time'))) {
       dim(time_expL) <- c(dim(expL)['time'])
       dim(excludeTime) <- c(dim(expL)['time'])
-  } else if (any(names(dim(expL)) %in% c('time_exp'))) {
-    dim(time_expL) <- c(dim(expL)['time_exp'])
-    dim(excludeTime) <- c(dim(expL)['time_exp'])
+    } else if (any(names(dim(expL)) %in% c('time_exp'))) {
+      dim(time_expL) <- c(dim(expL)['time_exp'])
+      dim(excludeTime) <- c(dim(expL)['time_exp'])
+    }
   }
- }
   if (!AnalogsInfo) {
     if (is.null(obsVar)) {
       res <- Apply(list(expL, obsL),
-                     target_dims = list(c('lat', 'lon'), c('time','lat','lon')),
-                     fun = .analogs, time_obsL, expVar = expVar,
-                     time_expL=time_expL, excludeTime=excludeTime,
-                     obsVar = obsVar, criteria = criteria,
-                     lonL = lonL, latL = latL,
-                     lonVar = lonVar, latVar = latVar, region = region,
-                     nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo,
-                     output_dims = c('nAnalogs', 'lat', 'lon'), 
-                     ncores = ncores)$output1
+                   target_dims = list(c('lat', 'lon'), c('time','lat','lon')),
+                   fun = .analogs, time_obsL, expVar = expVar,
+                   time_expL = time_expL, excludeTime = excludeTime,
+                   obsVar = obsVar, criteria = criteria,
+                   lonL = lonL, latL = latL,
+                   lonVar = lonVar, latVar = latVar, region = region,
+                   nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo,
+                   output_dims = c('nAnalogs', 'lat', 'lon'), 
+                   ncores = ncores)$output1
 
     } else if (!is.null(obsVar) && is.null(expVar)) {
         res <- Apply(list(expL, obsL, obsVar), 
                      target_dims = list(c('lat', 'lon'), c('time','lat','lon'),
                                         c('time', 'lat', 'lon')),
-                     fun = .analogs,time_obsL, 
-                     time_expL=time_expL, excludeTime=excludeTime,
+                     fun = .analogs, time_obsL, 
+                     time_expL = time_expL, excludeTime = excludeTime,
                      expVar = expVar, criteria = criteria,
                      lonL = lonL, latL = latL,
                      lonVar = lonVar, latVar = latVar, region = region,
@@ -669,8 +673,8 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
                      target_dims = list(c('lat', 'lon'), c('time','lat','lon'),
                                         c('time','lat','lon'), c('lat','lon')),
                      fun = .analogs,
-                     criteria = criteria,time_obsL,
-                     time_expL=time_expL, excludeTime=excludeTime,
+                     criteria = criteria, time_obsL,
+                     time_expL = time_expL, excludeTime = excludeTime,
                      lonL = lonL, latL = latL,
                      lonVar = lonVar, latVar = latVar, region = region,
                      nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo,
@@ -682,7 +686,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
         res <- Apply(list(expL, obsL),
                      target_dims = list(c('lat', 'lon'), c('time','lat','lon')),
                      fun = .analogs, time_obsL, expVar = expVar,
-                     time_expL=time_expL, excludeTime=excludeTime,
+                     time_expL = time_expL, excludeTime = excludeTime,
                      obsVar = obsVar, criteria = criteria,
                      lonL = lonL, latL = latL,
                      lonVar = lonVar, latVar = latVar, region = region,
@@ -697,7 +701,7 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
                      target_dims = list(c('lat', 'lon'), c('time','lat','lon'),
                                         c('time', 'lat', 'lon')),
                      fun = .analogs,time_obsL,
-                     time_expL=time_expL, excludeTime=excludeTime,
+                     time_expL = time_expL, excludeTime = excludeTime,
                      expVar = expVar, criteria = criteria,
                      lonL = lonL, latL = latL,
                      lonVar = lonVar, latVar = latVar, region = region,
@@ -712,9 +716,9 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
       res <- Apply(list(expL, obsL, obsVar, expVar),
                    target_dims = list(c('lat', 'lon'), c('time', 'lat', 'lon'),
                                       c('time', 'lat', 'lon'), c('lat', 'lon')),
-                   fun = .analogs,time_obsL, 
+                   fun = .analogs, time_obsL, 
                    criteria = criteria,
-                   time_expL=time_expL, excludeTime=excludeTime,
+                   time_expL = time_expL, excludeTime = excludeTime,
                    lonL = lonL, latL = latL,
                    lonVar = lonVar, latVar = latVar, region = region,
                    nAnalogs = nAnalogs, AnalogsInfo = AnalogsInfo,
@@ -727,14 +731,14 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
   }
   return(res)
 }
-.analogs <- function(expL, obsL, time_expL,  excludeTime = NULL,
-                    obsVar = NULL, expVar = NULL,
-                    time_obsL, criteria = "Large_dist",
-                    lonL = NULL, latL = NULL,
-                    lonVar = NULL, latVar = NULL, region = NULL, 
-                    nAnalogs = NULL, AnalogsInfo = FALSE) {
+.analogs <- function (expL, obsL, time_expL,  excludeTime = NULL,
+                      obsVar = NULL, expVar = NULL,
+                      time_obsL, criteria = "Large_dist",
+                      lonL = NULL, latL = NULL,
+                      lonVar = NULL, latVar = NULL, region = NULL, 
+                      nAnalogs = NULL, AnalogsInfo = FALSE) {
 
-  if (all(excludeTime=="")) {
+  if (all(excludeTime == "")) {
     excludeTime = NULL
   }
 
@@ -744,15 +748,15 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
       if (is.null(excludeTime)) {
         excludeTime <- time_expL
         warning("Parameter 'excludeTime' is NULL, time_obsL contains
-               time_expL, so, by default, the date of 
-               time_expL will be excluded in the search of analogs")
+                 time_expL, so, by default, the date of 
+                 time_expL will be excluded in the search of analogs")
       } else {
         `%!in%` = Negate(`%in%`)
         if(any(time_expL %!in% excludeTime)) {
           excludeTime <- c(excludeTime, time_expL)
           warning("Parameter 'excludeTime' is not NULL, time_obsL contains
-               time_expL, so, by default, the date of 
-               time_expL will be excluded in the search of analogs")
+                   time_expL, so, by default, the date of 
+                   time_expL will be excluded in the search of analogs")
         }
       }
       time_ref <- time_obsL[-c(which(time_obsL %in% excludeTime))]
@@ -771,10 +775,10 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
       if (is.null(excludeTime)) {
         if (!is.null(obsVar)) {
           warning("Parameter 'excludeTime' is NULL, time_obsL does not contain 
-                  time_expL, obsVar not NULL")
+                   time_expL, obsVar not NULL")
         } else {
           warning("Parameter 'excludeTime' is NULL, time_obsL does not contain 
-                  time_expL")
+                   time_expL")
         }
       } else {
         time_ref <- time_obsL[-c(which(time_obsL %in% excludeTime))]
@@ -791,17 +795,17 @@ Analogs <- function(expL, obsL, time_obsL, time_expL = NULL,
         }
         if (!is.null(obsVar)) {
           warning("Parameter 'excludeTime' has a value and time_obsL does not 
-                    contain time_expL, obsVar not NULL")
+                   contain time_expL, obsVar not NULL")
         } else {
           warning("Parameter 'excludeTime' has a value and time_obsL does not 
-                    contain time_expL")
+                   contain time_expL")
         }
       }
     }
   } else {
     stop("parameter 'obsL' cannot be NULL")
   }
-  if(length(time_obsL)==0){
+  if (length(time_obsL)==0) {
     stop("Parameter 'time_obsL' can not be length 0")
   }
     Analog_result <- FindAnalog(expL = expL, obsL = obsL, time_obsL = time_obsL,
@@ -845,13 +849,13 @@ FindAnalog <- function(expL, obsL, time_obsL, expVar, obsVar, criteria,
       Analogs_fields <- Subset(obsVar, 
                                along = which(names(dim(obsVar)) == 'time'),
                                indices = best)
-      warning("Parameter 'obsVar' is NULL and the returned field", 
+      warning("Parameter 'obsVar' is NULL and the returned field ", 
               "will be computed from 'obsL' (same variable).")
       
     } else {
       obslocal <- SelBox(obsVar, lon = lonVar, lat = latVar, 
-                         region = region)
-      Analogs_fields <- Subset(obslocal$data, 
+                         region = region)$data
+      Analogs_fields <- Subset(obslocal, 
                                along = which(names(dim(obslocal)) == 'time'),
                                indices = best)
     }
@@ -922,7 +926,7 @@ BestAnalog <- function(position, nAnalogs = nAnalogs, AnalogsInfo = FALSE,
     pos1 <- pos1[1 : nAnalogs]
     pos2 <- pos2[1 : nAnalogs]
     best <- match(pos1, pos2)
-    if(length(best)==1){
+    if(length(best)==1) {
       warning("Just 1 best analog matching Large_dist and ", 
               "Local_dist criteria")
     } 
@@ -940,7 +944,7 @@ BestAnalog <- function(position, nAnalogs = nAnalogs, AnalogsInfo = FALSE,
     pos1 <- pos1[1 : nAnalogs]
     pos2 <- pos2[1 : nAnalogs]
     best <- match(pos1, pos2)
-    if(length(best)==1){
+    if (length(best)==1) {
       warning("Just 1 best analog matching Large_dist and ", 
               "Local_dist criteria")
     } 
@@ -985,7 +989,7 @@ Select <- function(expL, obsL,  expVar = NULL, obsVar = NULL,
     margins <- c(1 : (length(dim(metric1))))[-dim_time_obs]
     pos1 <- apply(metric1, margins, order)      
     names(dim(pos1))[1] <- 'time'
-    metric1.original=metric1
+    metric1.original = metric1
     metric1 <-  apply(metric1, margins, sort)
     names(dim(metric1))[1] <- 'time'
     names(dim(metric1.original))=names(dim(metric1))
@@ -1071,8 +1075,8 @@ Select <- function(expL, obsL,  expVar = NULL, obsVar = NULL,
   } 
   result
 }
-.time_ref<- function(time_obsL,time_expL,excludeTime){
-  sameTime=which(time_obsL %in% time_expL)
+.time_ref <- function(time_obsL,time_expL,excludeTime){
+  sameTime = which(time_obsL %in% time_expL)
   result<- c(time_obsL[1:(sameTime-excludeTime-1)],
              time_obsL[(sameTime+excludeTime+1):length(time_obsL)])
   result
@@ -1105,15 +1109,17 @@ replace_time_dimnames <- function(dataL, time_name = 'time',
     stop("Parameter 'names_obs' must be a vector of characters.")
   }
   time_dim_obs <- which(names_obs == time_name | 
-                          names_obs == stdate_name | names_obs == ftime_name)
-  if(length(time_dim_obs) >1){
+                        names_obs == stdate_name | names_obs == ftime_name)
+  if (length(time_dim_obs) > 1) {
     stop ("more than 1 time dimension, please give just 1")
   }
-  if(length(time_dim_obs) == 0){
+  if (length(time_dim_obs) == 0) {
     warning ("name of time dimension is not 'ftime' or 'time' or 'stdate' 
-             or time dimension is null")
+              or time dimension is null")
+  }
+  if (length(time_dim_obs) != 0) {
+    names_obs[time_dim_obs]= time_name
   }
-  if(length(time_dim_obs)!=0){ names_obs[time_dim_obs]= time_name}
-  names(dim(dataL))=names_obs
+  names(dim(dataL)) = names_obs
   return(dataL)
 }
diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R
index e4410185f14a53ecb40f0091a8503b7eb7e589e4..a84b6fc8538b03f4113b96aa8b3126189a0bdee9 100644
--- a/R/CST_Anomaly.R
+++ b/R/CST_Anomaly.R
@@ -2,23 +2,43 @@
 #'
 #'@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.
+#'@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 filter_span a numeric value indicating the degree of smoothing. This option is only available if parameter \code{cross} is set to FALSE.  
-#'@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.
+#'@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 dim_anom A character string indicating the name of the dimension 
+#'  along which the climatology will be computed. The default value is 'sdate'.
+#'@param cross A logical value indicating whether cross-validation should be 
+#'  applied or not. Default = FALSE.
+#'@param memb_dim A character string indicating the name of the member 
+#'  dimension. It must be one dimension in 'exp' and 'obs'. If there is no 
+#'  member dimension, set NULL. The default value is 'member'.
+#'@param memb A logical value indicating whether to subtract the climatology 
+#'  based on the individual members (TRUE) or the ensemble mean over all
+#'  members (FALSE) when calculating the anomalies. The default value is TRUE.
+#'@param dat_dim A character vector indicating the name of the dataset and 
+#'  member dimensions. If there is no dataset dimension, it can be NULL.
+#'  The default value is "c('dataset', 'member')".
+#'@param filter_span A numeric value indicating the degree of smoothing. This 
+#'  option is only available if parameter \code{cross} is set to FALSE.
+#'@param ftime_dim A character string indicating the name of the temporal 
+#'  dimension where the smoothing with 'filter_span' will be applied. It cannot
+#'  be NULL if 'filter_span' is provided. The default value is 'ftime'.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL. It will be used only when 
+#'  'filter_span' is not NULL.
 #'
-#' @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.
-#'
-#'@importFrom s2dv InsertDim Clim Ano_CrossValid
+#'@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.
 #'
 #'@examples
 #'# Example 1:
@@ -34,78 +54,69 @@
 #'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)
-#'
 #'anom5 <- CST_Anomaly(lonlat_temp$exp)
-#'
 #'anom6 <- CST_Anomaly(obs = lonlat_temp$obs)
 #'
 #'@seealso \code{\link[s2dv]{Ano_CrossValid}}, \code{\link[s2dv]{Clim}} and \code{\link{CST_Load}}
 #'
-#'
+#'@import multiApply
+#'@importFrom s2dv InsertDim Clim Ano_CrossValid Reorder
 #'@export
-CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
-                        filter_span = NULL, dim_anom = 3) {
- 
+CST_Anomaly <- function(exp = NULL, obs = NULL, dim_anom = 'sdate', cross = FALSE, 
+                        memb_dim = 'member', memb = TRUE, dat_dim = c('dataset', 'member'), 
+                        filter_span = NULL, ftime_dim = 'ftime', ncores = NULL) {
+  # s2dv_cube
   if (!inherits(exp, 's2dv_cube') & !is.null(exp) || 
       !inherits(obs, 's2dv_cube') & !is.null(obs)) {
     stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
   }
-
- if (!is.null(obs)) {
-    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.")
-    }
-  }
-  case_exp = case_obs = 0
-  if (is.null(exp) & is.null(obs)) {
+  # exp and obs
+  if (is.null(exp$data) & is.null(obs$data)) {
     stop("One of the parameter 'exp' or 'obs' cannot be NULL.")
   }
+  case_exp = case_obs = 0
   if (is.null(exp)) {
     exp <- obs 
     case_obs = 1
-    warning("Parameter 'exp' is not provided and will be recycled.")
+    warning("Parameter 'exp' is not provided and 'obs' will be used instead.")
   }  
   if (is.null(obs)) {
     obs <- exp 
     case_exp = 1
-    warning("Parameter 'obs' is not provided and will be recycled.")
+    warning("Parameter 'obs' is not provided and 'exp' will be used instead.")
   }
-
-
-  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.")
+  if(any(is.null(names(dim(exp$data))))| any(nchar(names(dim(exp$data))) == 0) |
+     any(is.null(names(dim(obs$data))))| any(nchar(names(dim(obs$data))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names in element 'data'.")
+  }
+  if(!all(names(dim(exp$data)) %in% names(dim(obs$data))) |
+     !all(names(dim(obs$data)) %in% names(dim(exp$data)))) {
+    stop("Parameter 'exp' and 'obs' must have same dimension names in element 'data'.")
   }
   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) {
+  # dim_anom
+  if (is.numeric(dim_anom) & length(dim_anom) == 1) {
+    warning("Parameter 'dim_anom' must be a character string and a numeric value will not be ",
+            "accepted in the next release. The corresponding dimension name is assigned.")
+    dim_anom <- dimnames_data[dim_anom]
+  }
+  if (!is.character(dim_anom)) {
+    stop("Parameter 'dim_anom' must be a character string.")
+  }
+  if (!dim_anom %in% names(dim_exp) | !dim_anom %in% names(dim_obs)) {
+    stop("Parameter 'dim_anom' is not found in 'exp' or in 'obs' dimension in element 'data'.")
+  }
+  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'.")
-  } 
+  # cross
   if (!is.logical(cross) | !is.logical(memb) ) {
     stop("Parameters 'cross' and 'memb' must be logical.")
   }
@@ -114,89 +125,97 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
     warning("Parameter 'cross' has length greater than 1 and only the first element",
             "will be used.")
   }
+  # memb
   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)
+  # memb_dim
+  if (!is.null(memb_dim)) {
+    if (!is.character(memb_dim) | length(memb_dim) > 1) {
+      stop("Parameter 'memb_dim' must be a character string.")
+    }
+    if (!memb_dim %in% names(dim_exp) | !memb_dim %in% names(dim_obs)) {
+      stop("Parameter 'memb_dim' is not found in 'exp' or in 'obs' dimension.")
+    }
+  }
+  # dat_dim
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim)) {
+      stop("Parameter 'dat_dim' must be a character vector.")
+    }
+    if (!all(dat_dim %in% names(dim_exp)) | !all(dat_dim %in% names(dim_obs))) {
+      stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension in element 'data'.",
+           " Set it as NULL if there is no dataset dimension.")
+    }
+  }
+  # filter_span
+  if (!is.null(filter_span)) {
+    if (!is.numeric(filter_span)) {
+      warning("Paramater 'filter_span' is not numeric and any filter",
+                " is being applied.")
+      filter_span <- NULL
+    }
+    # ncores
+    if (!is.null(ncores)) {
+      if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+        length(ncores) > 1) {
+        stop("Parameter 'ncores' must be a positive integer.")
+      }
+    }
+    # ftime_dim
+    if (!is.character(ftime_dim)) {
+      stop("Parameter 'ftime_dim' must be a character string.")
+    }
+    if (!ftime_dim %in% names(dim_exp) | !memb_dim %in% names(dim_obs)) {
+      stop("Parameter 'ftime_dim' is not found in 'exp' or in 'obs' dimension in element 'data'.")
+    }
   }
-  
   
   # Computating anomalies
   #----------------------
   
   # With cross-validation
   if (cross) {
-    ano <- Ano_CrossValid(exp = exp$data, obs = obs$data, memb = memb)
-    # reorder dimension back
-    ano$exp <- aperm(ano$exp, match(names(dim(exp$data)), names(dim(ano$exp))))
-    ano$obs <- aperm(ano$obs, match(names(dim(obs$data)), names(dim(ano$obs))))
+    ano <- Ano_CrossValid(exp = exp$data, obs = obs$data, time_dim = dim_anom, memb_dim = memb_dim, memb = memb, dat_dim = dat_dim)
    
     #  Without cross-validation 
   } else {
-    tmp <- Clim(exp = exp$data, obs = obs$data, memb = memb)
+    tmp <- Clim(exp = exp$data, obs = obs$data, time_dim = dim_anom, memb_dim = memb_dim, memb = memb, dat_dim = dat_dim)
     if (!is.null(filter_span)) {
-        if (is.numeric(filter_span)) {
-            pos_dims <- names(dim(tmp$clim_exp))
-            reorder <- match(pos_dims, c('ftime',
-                             pos_dims[-which(pos_dims == 'ftime')]))
-            tmp$clim_obs <- aperm(apply(tmp$clim_obs, c(1 : 
-               length(dim(tmp$clim_obs)))[-which(names(dim(tmp$clim_obs)) == 'ftime')],
-               .Loess, loess_span = filter_span), reorder)
-            tmp$clim_exp <- aperm(apply(tmp$clim_exp, c(1 :
-               length(dim(tmp$clim_exp)))[-which(names(dim(tmp$clim_exp)) == 'ftime')],
-               .Loess, loess_span = filter_span), reorder)
-        } else {
-            warning("Paramater 'filter_span' is not numeric and any filter",
-                    " is being applied.")
-        }
+      tmp$clim_exp <- Apply(tmp$clim_exp, 
+                            target_dims = c(ftime_dim),
+                            output_dims = c(ftime_dim),
+                            fun = .Loess,
+                            loess_span = filter_span, 
+                            ncores = ncores)$output1
+      tmp$clim_obs <- Apply(tmp$clim_obs, 
+                            target_dims = c(ftime_dim),
+                            output_dims = c(ftime_dim),
+                            fun = .Loess,
+                            loess_span = filter_span, 
+                            ncores = ncores)$output1
     }
     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(tmp$clim_exp, 1,  dim_exp[memb_dim]) 
+      clim_obs <- InsertDim(tmp$clim_obs, 1,  dim_obs[memb_dim]) 
     }
-   
-    clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) 
-    clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) 
-    ano <- NULL    
+    clim_exp <- InsertDim(clim_exp, 1, dim_exp[dim_anom]) 
+    clim_obs <- InsertDim(clim_obs, 1, dim_obs[dim_anom])
+    ano <- NULL
+
+    # Permuting back dimensions to original order
+    clim_exp <- Reorder(clim_exp, dimnames_data)
+    clim_obs <- Reorder(clim_obs, dimnames_data)
+
     ano$exp <- exp$data - clim_exp
     ano$obs <- obs$data - clim_obs 
   }
 
-  # Permuting back dimensions to original order
-  if  (dim_anom != 3) {
-    
-    if (case_obs == 0) {
-      ano$exp <- aperm(ano$exp, perm = dimperm)
-    } 
-    if (case_exp == 0) {
-      ano$obs <- aperm(ano$obs, perm = dimperm)
-    }
-        
-    #Updating back permuted dimensions
-    dim_exp <- dim(exp$data)
-    dim_obs <- dim(obs$data)
-  }
-  
-  # Adding dimensions names
-  attr(ano$exp, 'dimensions') <- dimnames_data
-  attr(ano$obs, 'dimensions') <- dimnames_data
   exp$data <- ano$exp
   obs$data <- ano$obs
   
@@ -212,6 +231,7 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
     return(list(exp = exp, obs = obs)) 
   }
 }
+
 .Loess <- function(clim, loess_span) {
   data <- data.frame(ensmean = clim, day = 1 : length(clim))
   loess_filt <- loess(ensmean ~ day, data, span = loess_span)
diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R
index b6a52cb3c870d5c8cfab6ec526e5e818c5d066b2..7d639a49016b44752231d17458efac1324e15d1d 100644
--- a/R/PlotCombinedMap.R
+++ b/R/PlotCombinedMap.R
@@ -1,37 +1,77 @@
 #'Plot Multiple Lon-Lat Variables In a Single Map According to a Decision Function
-#'@description Plot a number a two dimensional matrices with (longitude, latitude) dimensions on a single map with the cylindrical equidistant latitude and longitude projection.
+#' 
+#'@description Plot a number a two dimensional matrices with (longitude, 
+#'latitude) dimensions on a single map with the cylindrical equidistant 
+#'latitude and longitude projection.
+#' 
 #'@author Nicolau Manubens, \email{nicolau.manubens@bsc.es}
 #'@author Veronica Torralba, \email{veronica.torralba@bsc.es}
 #'
-#'@param maps List of matrices to plot, each with (longitude, latitude) dimensions, or 3-dimensional array with the dimensions (longitude, latitude, map). Dimension names are required.
-#'@param lon Vector of longitudes. Must match the length of the corresponding dimension in 'maps'.
-#'@param lat Vector of latitudes. Must match the length of the corresponding dimension in 'maps'.
-#'@param map_select_fun Function that selects, for each grid point, which value to take among all the provided maps. This function receives as input a vector of values for a same grid point for all the provided maps, and must return a single selected value (not its index!) or NA. For example, the \code{min} and \code{max} functions are accepted.
-#'@param display_range Range of values to be displayed for all the maps. This must be a numeric vector c(range min, range max). The values in the parameter 'maps' can go beyond the limits specified in this range. If the selected value for a given grid point (according to 'map_select_fun') falls outside the range, it will be coloured with 'col_unknown_map'.
-#'@param map_dim Optional name for the dimension of 'maps' along which the multiple maps are arranged. Only applies when 'maps' is provided as a 3-dimensional array. Takes the value 'map' by default.
-#'@param brks Colour levels to be sent to PlotEquiMap. This parameter is optional and adjusted automatically by the function.
-#'@param cols List of vectors of colours to be sent to PlotEquiMap for the colour bar of each map. This parameter is optional and adjusted automatically by the function (up to 5 maps). The colours provided for each colour bar will be automatically interpolated to match the number of breaks. Each item in this list can be named, and the name will be used as title for the corresponding colour bar (equivalent to the parameter 'bar_titles').
-#'@param col_unknown_map Colour to use to paint the grid cells for which a map is not possible to be chosen according to 'map_select_fun' or for those values that go beyond 'display_range'. Takes the value 'white' by default.
-#'@param mask Optional numeric array with dimensions (latitude, longitude), with values in the range [0, 1], indicating the opacity of the mask over each grid point. Cells with a 0 will result in no mask, whereas cells with a 1 will result in a totally opaque superimposed pixel coloured in 'col_mask'.
-#'@param col_mask Colour to be used for the superimposed mask (if specified in 'mask'). Takes the value 'grey' by default.
+#'@param maps List of matrices to plot, each with (longitude, latitude) 
+#'  dimensions, or 3-dimensional array with the dimensions (longitude, latitude, 
+#'  map). Dimension names are required.
+#'@param lon Vector of longitudes. Must match the length of the corresponding 
+#'  dimension in 'maps'.
+#'@param lat Vector of latitudes. Must match the length of the corresponding 
+#'  dimension in 'maps'.
+#'@param map_select_fun Function that selects, for each grid point, which value 
+#'  to take among all the provided maps. This function receives as input a 
+#'  vector of values for a same grid point for all the provided maps, and must 
+#'  return a single selected value (not its index!) or NA. For example, the 
+#'  \code{min} and \code{max} functions are accepted.
+#'@param display_range Range of values to be displayed for all the maps. This 
+#'  must be a numeric vector c(range min, range max). The values in the 
+#'  parameter 'maps' can go beyond the limits specified in this range. If the 
+#'  selected value for a given grid point (according to 'map_select_fun') falls 
+#'  outside the range, it will be coloured with 'col_unknown_map'.
+#'@param map_dim Optional name for the dimension of 'maps' along which the 
+#'  multiple maps are arranged. Only applies when 'maps' is provided as a 
+#'  3-dimensional array. Takes the value 'map' by default.
+#'@param brks Colour levels to be sent to PlotEquiMap. This parameter is 
+#'  optional and adjusted automatically by the function.
+#'@param cols List of vectors of colours to be sent to PlotEquiMap for the 
+#'  colour bar of each map. This parameter is optional and adjusted 
+#'  automatically by the function (up to 5 maps). The colours provided for each 
+#'  colour bar will be automatically interpolated to match the number of breaks. 
+#'  Each item in this list can be named, and the name will be used as title for 
+#'  the corresponding colour bar (equivalent to the parameter 'bar_titles').
+#'@param col_unknown_map Colour to use to paint the grid cells for which a map 
+#'  is not possible to be chosen according to 'map_select_fun' or for those 
+#'  values that go beyond 'display_range'. Takes the value 'white' by default.
+#'@param mask Optional numeric array with dimensions (latitude, longitude), with 
+#'  values in the range [0, 1], indicating the opacity of the mask over each 
+#'  grid point. Cells with a 0 will result in no mask, whereas cells with a 1 
+#'  will result in a totally opaque superimposed pixel coloured in 'col_mask'.
+#'@param col_mask Colour to be used for the superimposed mask (if specified in 
+#'  'mask'). Takes the value 'grey' by default.
 #'@param dots Array of same dimensions as 'var' or with dimensions 
 #'  c(n, dim(var)), where n is the number of dot/symbol layers to add to the 
 #'  plot. A value of TRUE at a grid cell will draw a dot/symbol on the 
 #'  corresponding square of the plot. By default all layers provided in 'dots' 
 #'  are plotted with dots, but a symbol can be specified for each of the 
 #'  layers via the parameter 'dot_symbol'.
-#'@param bar_titles Optional vector of character strings providing the titles to be shown on top of each of the colour bars.
-#'@param legend_scale Scale factor for the size of the colour bar labels. Takes 1 by default.
-#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 by default.
+#'@param bar_titles Optional vector of character strings providing the titles to 
+#'  be shown on top of each of the colour bars.
+#'@param legend_scale Scale factor for the size of the colour bar labels. Takes 
+#'  1 by default.
+#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 
+#'  by default.
 #'@param plot_margin Numeric vector of length 4 for the margin sizes in the 
 #'  following order: bottom, left, top, and right. If not specified, use the 
 #'  default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in 
 #'  s2dv::PlotEquiMap.
-#'@param fileout File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff
-#'@param width File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default.
-#'@param height File height, in the units specified in the parameter size_units (inches by default). Takes 5 by default.
-#'@param size_units Units of the size of the device (file or window) to plot in. Inches ('in') by default. See ?Devices and the creator function of the corresponding device.
-#'@param res Resolution of the device (file or window) to plot in. See ?Devices and the creator function of the corresponding device.
+#'@param fileout File where to save the plot. If not specified (default) a 
+#'  graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp 
+#'  and tiff
+#'@param width File width, in the units specified in the parameter size_units 
+#'  (inches by default). Takes 8 by default.
+#'@param height File height, in the units specified in the parameter size_units 
+#'  (inches by default). Takes 5 by default.
+#'@param size_units Units of the size of the device (file or window) to plot in. 
+#'  Inches ('in') by default. See ?Devices and the creator function of the 
+#'  corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. See ?Devices 
+#'  and the creator function of the corresponding device.
 #'@param drawleg Where to draw the common colour bar. Can take values TRUE, 
 #'  FALSE or:\cr
 #'  'up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N'\cr
@@ -39,13 +79,7 @@
 #'  'right', 'r', 'R', 'east', 'e', 'E'\cr
 #'  'left', 'l', 'L', 'west', 'w', 'W' 
 #'@param ... Additional parameters to be passed on to \code{PlotEquiMap}.
-
-#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap}
 #' 
-#'@importFrom s2dv PlotEquiMap ColorBar
-#'@importFrom maps map 
-#'@importFrom graphics box image layout mtext par plot.new
-#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off hcl jpeg pdf png postscript svg tiff
 #'@examples
 #'# Simple example
 #'x <- array(1:(20 * 10), dim = c(lat = 10, lon = 20)) / 200
@@ -60,7 +94,7 @@
 #'                map_select_fun = max,
 #'                display_range = c(0, 1),
 #'                bar_titles = paste('% of belonging to', c('a', 'b', 'c')), 
-#'                brks = 20, width = 10, height = 8)
+#'                brks = 20, width = 12, height = 10)
 #'}
 #'
 #'Lon <- c(0:40, 350:359)
@@ -72,9 +106,16 @@
 #'\dontrun{
 #'PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max,
 #'                display_range = range(data), mask = mask,
-#'                width = 12, height = 8) 
+#'                width = 14, height = 10) 
 #'}
 #'
+#'@seealso \code{PlotCombinedMap} and \code{PlotEquiMap}
+#' 
+#'@importFrom s2dv PlotEquiMap ColorBar
+#'@importFrom maps map 
+#'@importFrom graphics box image layout mtext par plot.new
+#'@importFrom grDevices adjustcolor bmp colorRampPalette dev.cur dev.new dev.off 
+#'  hcl jpeg pdf png postscript svg tiff
 #'@export
 PlotCombinedMap <- function(maps, lon, lat, 
                             map_select_fun, display_range, 
@@ -218,7 +259,7 @@ PlotCombinedMap <- function(maps, lon, lat,
                                    brks = brks, cols = cols, vertical = FALSE,
                                    subsampleg = NULL, bar_limits = display_range, var_limits = NULL,
                                    triangle_ends = NULL, plot = FALSE, draw_separators = TRUE,
-                                   bar_titles = bar_titles, title_scale = 1, label_scale = legend_scale * 1.5,
+                                   bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5,
                                    extra_margin = c(2, 0, 2, 0))
   
   # Check legend_scale
@@ -394,7 +435,7 @@ PlotCombinedMap <- function(maps, lon, lat,
                          brks = colorbar$brks, cols = colorbar$cols, vertical = FALSE,
                          subsampleg = NULL, bar_limits = display_range, var_limits = NULL,
                          triangle_ends = NULL, plot = TRUE, draw_separators = TRUE,
-                         bar_titles = bar_titles, title_scale = 1, label_scale = legend_scale * 1.5,
+                         bar_titles = bar_titles, title_scale = cex_bar_titles, label_scale = legend_scale * 1.5,
                          extra_margin = c(2, 0, 2, 0))
   }
   
diff --git a/man/Analogs.Rd b/man/Analogs.Rd
index fc26a5523fe0dd78c4755b03f106118d99b193b3..a7addc73e5bf4936c06e29be410841aae5519619 100644
--- a/man/Analogs.Rd
+++ b/man/Analogs.Rd
@@ -20,56 +20,56 @@ Analogs(
   region = NULL,
   nAnalogs = NULL,
   AnalogsInfo = FALSE,
-  ncores = 1
+  ncores = NULL
 )
 }
 \arguments{
-\item{expL}{an array of N named dimensions containing the experimental field
+\item{expL}{An array of N named dimensions containing the experimental field
 on the large scale for which the analog is aimed. This field is used to in 
 all the criterias. If parameter 'expVar' is not provided, the function will
 return the expL analog. The element 'data' in the 's2dv_cube' object must 
-have, at least, latitudinal and longitudinal dimensions. The object is expect
-to be already subset for the desired large scale region.}
+have, at least, latitudinal and longitudinal dimensions. The object is 
+expect to be already subset for the desired large scale region.}
 
-\item{obsL}{an array of N named dimensions containing the observational field 
+\item{obsL}{An array of N named dimensions containing the observational field 
 on the large scale. The element 'data' in the 's2dv_cube' object must have 
 the same latitudinal and longitudinal dimensions as parameter 'expL' and a 
 single temporal dimension with the maximum number of available observations.}
 
-\item{time_obsL}{a character string indicating the date of the observations 
+\item{time_obsL}{A character string indicating the date of the observations 
 in the format "dd/mm/yyyy". Reference time to search for analogs.}
 
-\item{time_expL}{an array of N named dimensions (coinciding with time 
+\item{time_expL}{An array of N named dimensions (coinciding with time 
 dimensions in expL) of character string(s) indicating the date(s) of the 
 experiment in the format "dd/mm/yyyy". Time(s) to find the analogs.}
 
-\item{lonL}{a vector containing the longitude of parameter 'expL'.}
+\item{lonL}{A vector containing the longitude of parameter 'expL'.}
 
-\item{latL}{a vector containing the latitude of parameter 'expL'.}
+\item{latL}{A vector containing the latitude of parameter 'expL'.}
 
-\item{expVar}{an array of N named dimensions containing the experimental 
+\item{expVar}{An array of N named dimensions containing the experimental 
 field on the local scale, usually a different variable to the parameter 
 'expL'. If it is not NULL (by default, NULL), the returned field by this 
 function will be the analog of parameter 'expVar'.}
 
-\item{obsVar}{an array of N named dimensions containing the field of the 
+\item{obsVar}{An array of N named dimensions containing the field of the 
 same variable as the passed in parameter 'expVar' for the same region.}
 
 \item{criteria}{a character string indicating the criteria to be used for the
 selection of analogs:
-\itemize{
-\item{Large_dist} minimum Euclidean distance in the large scale pattern; 
-\item{Local_dist} minimum Euclidean distance in the large scale pattern 
-and minimum Euclidean distance in the local scale pattern; and
-\item{Local_cor} minimum Euclidean distance in the large scale pattern, 
-minimum Euclidean distance in the local scale pattern and highest 
-correlation in the local variable to downscale.}}
-
-\item{excludeTime}{an array of N named dimensions (coinciding with time 
+ \itemize{\item{Large_dist} minimum Euclidean distance in the large scale pattern; 
+          \item{Local_dist} minimum Euclidean distance in the large scale pattern 
+           and minimum Euclidean distance in the local scale pattern; and
+          \item{Local_cor} minimum Euclidean distance in the large scale pattern, 
+           minimum Euclidean distance in the local scale pattern and highest 
+           correlation in the local variable to downscale.}}
+
+\item{excludeTime}{An array of N named dimensions (coinciding with time 
 dimensions in expL) of character string(s) indicating the date(s) of the 
 observations in the format "dd/mm/yyyy" to be excluded during the search of 
-analogs. It can be NULL but if expL is not a forecast (time_expL contained in 
-time_obsL),by default time_expL will be removed during the search of analogs.}
+analogs. It can be NULL but if expL is not a forecast (time_expL contained 
+in time_obsL), by default time_expL will be removed during the search of 
+analogs.}
 
 \item{lonVar}{a vector containing the longitude of parameter 'expVar'.}
 
@@ -88,23 +88,25 @@ NULL for 'Local_dist' and 'Local_cor' the default value will be set at the
 length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
 the best analog.}
 
-\item{AnalogsInfo}{TRUE to get a list with two elements: 1) the downscaled 
-field and 2) the AnalogsInfo which contains: a) the number of the best 
+\item{AnalogsInfo}{A logical value. If it is TRUE it returns a list
+with two elements: 1) the downscaled field and
+2) the AnalogsInfo which contains: a) the number of the best 
 analogs, b) the corresponding value of the metric used in the selected 
 criteria (distance values for Large_dist and Local_dist,correlation values 
-for Local_cor), c)dates of the analogs). The analogs are listed in decreasing
-order, the first one is the best analog (i.e if the selected criteria is 
-Local_cor the best analog will be the one with highest correlation, while for
-Large_dist criteria the best analog will be the day with minimum Euclidean 
-distance). Set to FALSE to get a single analog, the best analog, for instance
-for downscaling.}
+for Local_cor), c)dates of the analogs). The analogs are listed in 
+decreasing order, the first one is the best analog (i.e if the selected 
+criteria is Local_cor the best analog will be the one with highest 
+correlation, while for Large_dist criteria the best analog will be the day 
+with minimum Euclidean distance). Set to FALSE to get a single analog, the 
+best analog, for instance for downscaling.}
 
 \item{ncores}{the number of cores to use in parallel computation.}
 }
 \value{
-AnalogsFields, dowscaled values of the best analogs for the criteria 
-selected. If AnalogsInfo is set to TRUE the function also returns a 
-list with the dowsncaled field and the Analogs Information.
+An array with the dowscaled values of the best analogs for the criteria 
+selected. If 'AnalogsInfo' is set to TRUE it returns a list with an array 
+of the dowsncaled field and the analogs information in elements 'analogs', 
+'metric' and 'dates'.
 }
 \description{
 This function perform a downscaling using Analogs. To compute 
@@ -156,7 +158,7 @@ downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, obsVar = obs.pr,
 obsSLP <- c(rnorm(1:1980), expSLP * 1.5)
 dim(obsSLP) <- c(lat = 4, lon = 5, time = 100)
 time_obsSLP <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-")
-downscale_field<- Analogs(expL = expSLP, obsL = obsSLP, time_obsSLP,
+downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, time_obsSLP,
                          nAnalogs = 5, time_expL = "01-01-2003", 
                          AnalogsInfo = TRUE, excludeTime = "01-01-2003")
 
@@ -170,7 +172,7 @@ downscale_field <- Analogs(expL = expSLP, obsL = obsSLP, obsVar = obs.pr,
 
 # Example 5: Downscaling using criteria 'Local_dist' and 2 variables:
 # analogs of local scale using criteria 2
-region=c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33)
+region = c(lonmin = -1 ,lonmax = 2, latmin = 30, latmax = 33)
 Local_scale <- Analogs(expL = expSLP, obsL = obsSLP, time_obsL = time_obsSLP,
                       obsVar = obs.pr, criteria = "Local_dist", 
                       lonL = seq(-1, 5, 1.5),latL = seq(30, 35, 1.5), 
diff --git a/man/CST_Analogs.Rd b/man/CST_Analogs.Rd
index 7a67b0fb2facf261acb5e28f91321ecdcee816ba..b0242f584023df4d157937e7f92d4c972791fa9e 100644
--- a/man/CST_Analogs.Rd
+++ b/man/CST_Analogs.Rd
@@ -16,34 +16,34 @@ CST_Analogs(
   time_obsL = NULL,
   nAnalogs = NULL,
   AnalogsInfo = FALSE,
-  ncores = 1
+  ncores = NULL
 )
 }
 \arguments{
-\item{expL}{an 's2dv_cube' object containing the experimental field on the 
+\item{expL}{An 's2dv_cube' object containing the experimental field on the 
 large scale for which the analog is aimed. This field is used to in all the
 criterias. If parameter 'expVar' is not provided, the function will return 
 the expL analog. The element 'data' in the 's2dv_cube' object must have, at
 least, latitudinal and longitudinal dimensions. The object is expect to be 
 already subset for the desired large scale region.}
 
-\item{obsL}{an 's2dv_cube' object containing the observational field on the
+\item{obsL}{An 's2dv_cube' object containing the observational field on the
 large scale. The element 'data' in the 's2dv_cube' object must have the same 
 latitudinal and longitudinal dimensions as parameter 'expL' and a temporal 
 dimension with the maximum number of available observations.}
 
-\item{expVar}{an 's2dv_cube' object containing the experimental field on the
+\item{expVar}{An 's2dv_cube' object containing the experimental field on the
 local scale, usually a different variable to the parameter 'expL'. If it is 
 not NULL (by default, NULL), the returned field by this function will be the 
 analog of parameter 'expVar'.}
 
-\item{obsVar}{an 's2dv_cube' containing the field of the same variable as the 
+\item{obsVar}{An 's2dv_cube' containing the field of the same variable as the 
 passed in parameter 'expVar' for the same region.}
 
-\item{region}{a vector of length four indicating the minimum longitude, 
+\item{region}{A vector of length four indicating the minimum longitude, 
 the maximum longitude, the minimum latitude and the maximum latitude.}
 
-\item{criteria}{a character string indicating the criteria to be used for the
+\item{criteria}{A character string indicating the criteria to be used for the
 selection of analogs:
 \itemize{
 \item{Large_dist} minimum Euclidean distance in the large scale pattern; 
@@ -55,21 +55,21 @@ correlation in the local variable to downscale.}
 Criteria 'Large_dist' is recommended for CST_Analogs, for an advanced use of 
 the criterias 'Local_dist' and 'Local_cor' use 'Analogs' function.}
 
-\item{excludeTime}{an array of N named dimensions (coinciding with time 
+\item{excludeTime}{An array of N named dimensions (coinciding with time 
 dimensions in expL)of character string(s) indicating the date(s) of the 
 observations in the format "dd/mm/yyyy" to be excluded during the search of 
 analogs. It can be NULL but if expL is not a forecast (time_expL contained in
 time_obsL), by default time_expL will be removed during the search of analogs.}
 
-\item{time_expL}{a character string indicating the date of the experiment 
+\item{time_expL}{A character string indicating the date of the experiment 
 in the same format than time_obsL (i.e. "yyyy-mm-dd"). By default it is NULL 
 and dates are taken from element \code{$Dates$start} from expL.}
 
-\item{time_obsL}{a character string indicating the date of the observations 
+\item{time_obsL}{A character string indicating the date of the observations 
 in the date format (i.e. "yyyy-mm-dd"). By default it is NULL and dates are 
 taken from element \code{$Dates$start} from obsL.}
 
-\item{nAnalogs}{number of Analogs to be selected to apply the criterias 
+\item{nAnalogs}{Number of Analogs to be selected to apply the criterias 
 'Local_dist' or 'Local_cor'. This is not the necessary the number of analogs 
 that the user can get, but the number of events with minimum distance in 
 which perform the search of the best Analog. The default value for the 
@@ -79,22 +79,24 @@ NULL for 'Local_dist' and 'Local_cor' the default value will be set at the
 length of 'time_obsL'. If AnalogsInfo is FALSE the function returns just 
 the best analog.}
 
-\item{AnalogsInfo}{TRUE to get a list with two elements: 1) the downscaled 
-field and 2) the AnalogsInfo which contains: a) the number of the best 
-analogs, b) the corresponding value of the metric used in the selected 
-criteria (distance values for Large_dist and Local_dist,correlation values 
-for Local_cor), c)dates of the analogs). The analogs are listed in decreasing
-order, the first one is the best analog (i.e if the selected criteria is 
-Local_cor the best analog will be the one with highest correlation, while for
-Large_dist criteria the best analog will be the day with minimum Euclidean 
-distance). Set to FALSE to get a single analog, the best analog, for instance
-for downscaling.}
+\item{AnalogsInfo}{A logical value. TRUE to get a list with two elements: 
+1) the downscaled field and 2) the AnalogsInfo which contains: 
+a) the number of the best analogs, b) the corresponding value of the metric 
+used in the selected criteria (distance values for Large_dist and Local_dist, 
+correlation values for Local_cor), c)dates of the analogs). The analogs are 
+listed in decreasing order, the first one is the best analog (i.e if the 
+selected criteria is Local_cor the best analog will be the one with highest 
+correlation, while for Large_dist criteria the best analog will be the day 
+with minimum Euclidean distance). Set to FALSE to get a single analog, the 
+best analog, for instance for downscaling.}
 
 \item{ncores}{The number of cores to use in parallel computation}
 }
 \value{
-An 'array' object containing the dowscaled values of the best 
-analogs.
+An 's2dv_cube' object containing an array with the dowscaled values of 
+the best analogs in element 'data'. If 'AnalogsInfo' is TRUE, 'data' is a list 
+with an array of the downscaled fields and the analogs information in 
+elements 'analogs', 'metric' and 'dates'.
 }
 \description{
 This function perform a downscaling using Analogs. To compute 
@@ -124,10 +126,10 @@ function within 'CSTools' package.
 }
 \examples{
 expL <- rnorm(1:200)
-dim(expL) <- c(member=10,lat = 4, lon = 5)
+dim(expL) <- c(member = 10,lat = 4, lon = 5)
 obsL <- c(rnorm(1:180),expL[1,,]*1.2)
 dim(obsL) <- c(time = 10,lat = 4, lon = 5)
-time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-")
+time_obsL <- paste(rep("01", 10), rep("01", 10), 1994:2003, sep = "-")
 time_expL <- time_obsL[1]
 lon <- seq(-1,5,1.5)
 lat <- seq(30,35,1.5)
@@ -137,6 +139,7 @@ obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon,
                  Dates = list(start = time_obsL, end = time_obsL))
 region <- c(min(lon), max(lon), min(lat), max(lat))
 downscaled_field <- CST_Analogs(expL = expL, obsL = obsL, region = region)
+
 }
 \references{
 Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard,
diff --git a/man/CST_Anomaly.Rd b/man/CST_Anomaly.Rd
index 06c78c8958a07baeffc6676eafb8fadd5a85fd71..3af85b5fc6ca2d742cf29bd6e0a1159f54dc088e 100644
--- a/man/CST_Anomaly.Rd
+++ b/man/CST_Anomaly.Rd
@@ -7,35 +7,65 @@
 CST_Anomaly(
   exp = NULL,
   obs = NULL,
+  dim_anom = "sdate",
   cross = FALSE,
+  memb_dim = "member",
   memb = TRUE,
+  dat_dim = c("dataset", "member"),
   filter_span = NULL,
-  dim_anom = 3
+  ftime_dim = "ftime",
+  ncores = NULL
 )
 }
 \arguments{
-\item{exp}{an object of class \code{s2dv_cube} as returned by \code{CST_Load} function, containing the seasonal forecast experiment data in the element named \code{$data}.}
+\item{exp}{An object of class \code{s2dv_cube} as returned by \code{CST_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{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{dim_anom}{A character string indicating the name of the dimension 
+along which the climatology will be computed. The default value is 'sdate'.}
 
-\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{cross}{A logical value indicating whether cross-validation should be 
+applied or not. Default = FALSE.}
 
-\item{filter_span}{a numeric value indicating the degree of smoothing. This option is only available if parameter \code{cross} is set to FALSE.}
+\item{memb_dim}{A character string indicating the name of the member 
+dimension. It must be one dimension in 'exp' and 'obs'. If there is no 
+member dimension, set NULL. The default value is 'member'.}
 
-\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.}
+\item{memb}{A logical value indicating whether to subtract the climatology 
+based on the individual members (TRUE) or the ensemble mean over all
+members (FALSE) when calculating the anomalies. The default value is TRUE.}
+
+\item{dat_dim}{A character vector indicating the name of the dataset and 
+member dimensions. If there is no dataset dimension, it can be NULL.
+The default value is "c('dataset', 'member')".}
+
+\item{filter_span}{A numeric value indicating the degree of smoothing. This 
+option is only available if parameter \code{cross} is set to FALSE.}
+
+\item{ftime_dim}{A character string indicating the name of the temporal 
+dimension where the smoothing with 'filter_span' will be applied. It cannot
+be NULL if 'filter_span' is provided. The default value is 'ftime'.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL. It will be used only when 
+'filter_span' is not NULL.}
 }
 \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.
+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.
+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:
@@ -51,18 +81,10 @@ 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)
-
 anom5 <- CST_Anomaly(lonlat_temp$exp)
-
 anom6 <- CST_Anomaly(obs = lonlat_temp$obs)
 
 }
diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd
index cf08c8f623698f1c545d53f91c7727cde937cd46..d013f80e5aa5bf42c70bfdb8b4faa34990c24813 100644
--- a/man/PlotCombinedMap.Rd
+++ b/man/PlotCombinedMap.Rd
@@ -31,27 +31,53 @@ PlotCombinedMap(
 )
 }
 \arguments{
-\item{maps}{List of matrices to plot, each with (longitude, latitude) dimensions, or 3-dimensional array with the dimensions (longitude, latitude, map). Dimension names are required.}
-
-\item{lon}{Vector of longitudes. Must match the length of the corresponding dimension in 'maps'.}
-
-\item{lat}{Vector of latitudes. Must match the length of the corresponding dimension in 'maps'.}
-
-\item{map_select_fun}{Function that selects, for each grid point, which value to take among all the provided maps. This function receives as input a vector of values for a same grid point for all the provided maps, and must return a single selected value (not its index!) or NA. For example, the \code{min} and \code{max} functions are accepted.}
-
-\item{display_range}{Range of values to be displayed for all the maps. This must be a numeric vector c(range min, range max). The values in the parameter 'maps' can go beyond the limits specified in this range. If the selected value for a given grid point (according to 'map_select_fun') falls outside the range, it will be coloured with 'col_unknown_map'.}
-
-\item{map_dim}{Optional name for the dimension of 'maps' along which the multiple maps are arranged. Only applies when 'maps' is provided as a 3-dimensional array. Takes the value 'map' by default.}
-
-\item{brks}{Colour levels to be sent to PlotEquiMap. This parameter is optional and adjusted automatically by the function.}
-
-\item{cols}{List of vectors of colours to be sent to PlotEquiMap for the colour bar of each map. This parameter is optional and adjusted automatically by the function (up to 5 maps). The colours provided for each colour bar will be automatically interpolated to match the number of breaks. Each item in this list can be named, and the name will be used as title for the corresponding colour bar (equivalent to the parameter 'bar_titles').}
-
-\item{col_unknown_map}{Colour to use to paint the grid cells for which a map is not possible to be chosen according to 'map_select_fun' or for those values that go beyond 'display_range'. Takes the value 'white' by default.}
-
-\item{mask}{Optional numeric array with dimensions (latitude, longitude), with values in the range [0, 1], indicating the opacity of the mask over each grid point. Cells with a 0 will result in no mask, whereas cells with a 1 will result in a totally opaque superimposed pixel coloured in 'col_mask'.}
-
-\item{col_mask}{Colour to be used for the superimposed mask (if specified in 'mask'). Takes the value 'grey' by default.}
+\item{maps}{List of matrices to plot, each with (longitude, latitude) 
+dimensions, or 3-dimensional array with the dimensions (longitude, latitude, 
+map). Dimension names are required.}
+
+\item{lon}{Vector of longitudes. Must match the length of the corresponding 
+dimension in 'maps'.}
+
+\item{lat}{Vector of latitudes. Must match the length of the corresponding 
+dimension in 'maps'.}
+
+\item{map_select_fun}{Function that selects, for each grid point, which value 
+to take among all the provided maps. This function receives as input a 
+vector of values for a same grid point for all the provided maps, and must 
+return a single selected value (not its index!) or NA. For example, the 
+\code{min} and \code{max} functions are accepted.}
+
+\item{display_range}{Range of values to be displayed for all the maps. This 
+must be a numeric vector c(range min, range max). The values in the 
+parameter 'maps' can go beyond the limits specified in this range. If the 
+selected value for a given grid point (according to 'map_select_fun') falls 
+outside the range, it will be coloured with 'col_unknown_map'.}
+
+\item{map_dim}{Optional name for the dimension of 'maps' along which the 
+multiple maps are arranged. Only applies when 'maps' is provided as a 
+3-dimensional array. Takes the value 'map' by default.}
+
+\item{brks}{Colour levels to be sent to PlotEquiMap. This parameter is 
+optional and adjusted automatically by the function.}
+
+\item{cols}{List of vectors of colours to be sent to PlotEquiMap for the 
+colour bar of each map. This parameter is optional and adjusted 
+automatically by the function (up to 5 maps). The colours provided for each 
+colour bar will be automatically interpolated to match the number of breaks. 
+Each item in this list can be named, and the name will be used as title for 
+the corresponding colour bar (equivalent to the parameter 'bar_titles').}
+
+\item{col_unknown_map}{Colour to use to paint the grid cells for which a map 
+is not possible to be chosen according to 'map_select_fun' or for those 
+values that go beyond 'display_range'. Takes the value 'white' by default.}
+
+\item{mask}{Optional numeric array with dimensions (latitude, longitude), with 
+values in the range [0, 1], indicating the opacity of the mask over each 
+grid point. Cells with a 0 will result in no mask, whereas cells with a 1 
+will result in a totally opaque superimposed pixel coloured in 'col_mask'.}
+
+\item{col_mask}{Colour to be used for the superimposed mask (if specified in 
+'mask'). Takes the value 'grey' by default.}
 
 \item{dots}{Array of same dimensions as 'var' or with dimensions 
 c(n, dim(var)), where n is the number of dot/symbol layers to add to the 
@@ -60,26 +86,36 @@ corresponding square of the plot. By default all layers provided in 'dots'
 are plotted with dots, but a symbol can be specified for each of the 
 layers via the parameter 'dot_symbol'.}
 
-\item{bar_titles}{Optional vector of character strings providing the titles to be shown on top of each of the colour bars.}
+\item{bar_titles}{Optional vector of character strings providing the titles to 
+be shown on top of each of the colour bars.}
 
-\item{legend_scale}{Scale factor for the size of the colour bar labels. Takes 1 by default.}
+\item{legend_scale}{Scale factor for the size of the colour bar labels. Takes 
+1 by default.}
 
-\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 by default.}
+\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 
+by default.}
 
 \item{plot_margin}{Numeric vector of length 4 for the margin sizes in the 
 following order: bottom, left, top, and right. If not specified, use the 
 default of par("mar"), c(5.1, 4.1, 4.1, 2.1). Used as 'margin_scale' in 
 s2dv::PlotEquiMap.}
 
-\item{fileout}{File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff}
+\item{fileout}{File where to save the plot. If not specified (default) a 
+graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp 
+and tiff}
 
-\item{width}{File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default.}
+\item{width}{File width, in the units specified in the parameter size_units 
+(inches by default). Takes 8 by default.}
 
-\item{height}{File height, in the units specified in the parameter size_units (inches by default). Takes 5 by default.}
+\item{height}{File height, in the units specified in the parameter size_units 
+(inches by default). Takes 5 by default.}
 
-\item{size_units}{Units of the size of the device (file or window) to plot in. Inches ('in') by default. See ?Devices and the creator function of the corresponding device.}
+\item{size_units}{Units of the size of the device (file or window) to plot in. 
+Inches ('in') by default. See ?Devices and the creator function of the 
+corresponding device.}
 
-\item{res}{Resolution of the device (file or window) to plot in. See ?Devices and the creator function of the corresponding device.}
+\item{res}{Resolution of the device (file or window) to plot in. See ?Devices 
+and the creator function of the corresponding device.}
 
 \item{drawleg}{Where to draw the common colour bar. Can take values TRUE, 
 FALSE or:\cr
@@ -91,7 +127,9 @@ FALSE or:\cr
 \item{...}{Additional parameters to be passed on to \code{PlotEquiMap}.}
 }
 \description{
-Plot a number a two dimensional matrices with (longitude, latitude) dimensions on a single map with the cylindrical equidistant latitude and longitude projection.
+Plot a number a two dimensional matrices with (longitude, 
+latitude) dimensions on a single map with the cylindrical equidistant 
+latitude and longitude projection.
 }
 \examples{
 # Simple example
@@ -107,7 +145,7 @@ PlotCombinedMap(list(a, b, c), lons, lats,
                map_select_fun = max,
                display_range = c(0, 1),
                bar_titles = paste('\% of belonging to', c('a', 'b', 'c')), 
-               brks = 20, width = 10, height = 8)
+               brks = 20, width = 12, height = 10)
 }
 
 Lon <- c(0:40, 350:359)
@@ -119,7 +157,7 @@ dim(mask) <- c(lat = 26, lon = 51)
 \dontrun{
 PlotCombinedMap(data, lon = Lon, lat = Lat, map_select_fun = max,
                display_range = range(data), mask = mask,
-               width = 12, height = 8) 
+               width = 14, height = 10) 
 }
 
 }
diff --git a/tests/testthat/test-CST_Anomaly.R b/tests/testthat/test-CST_Anomaly.R
new file mode 100644
index 0000000000000000000000000000000000000000..b4137015bd494d3cf1b911e5f214e87b38d5db97
--- /dev/null
+++ b/tests/testthat/test-CST_Anomaly.R
@@ -0,0 +1,178 @@
+context("CSTools::CST_Anomaly tests")
+
+##############################################
+# dat
+set.seed(1)
+mod <- rnorm(2 * 3 * 4 * 5 * 6 * 7)
+dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
+set.seed(2)
+obs <- rnorm(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'
+
+# dat1
+exp1 <- exp
+exp1$data <- NULL
+
+# dat2
+exp2 <- exp
+exp2$data <- array(rnorm(2 * 3 * 4 * 5 * 6 * 7), dim = c(var = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7))
+
+obs2 <- obs
+obs2$data <- array(rnorm(2 * 3 * 2 * 5 * 6 * 1), dim = c(member = 1, sdate = 2, ftime = 5, lat = 6, lon = 7))
+
+# dat3
+exp3 <- exp
+obs3 <- obs
+exp3$data <- array(rnorm(2 * 3 * 1 * 5 * 6 * 7), dim = c(dataset = 2, member = 3, sdate = 1, ftime = 5, lat = 6, lon = 7))
+obs3$data <- array(rnorm(2 * 3 * 1 * 5 * 6 * 1), dim = c(dataset = 2, member = 1, sdate = 1, ftime = 5, lat = 6, lon = 7))
+
+# dat4
+set.seed(1)
+mod4 <- rnorm(2 * 3 * 4 * 5 * 6 * 7)
+dim(mod4) <- c(datasets = 2, members = 3, sdates = 4, ftimes = 5, lat = 6, lon = 7)
+set.seed(2)
+obs4 <- rnorm(1 * 1 * 4 * 5 * 6 * 7)
+dim(obs4) <- c(datasets = 1, members = 1, sdates = 4, ftimes = 5, lat = 6, lon = 7)
+lon <- seq(0, 30, 5)
+lat <- seq(0, 25, 5)
+exp4 <- list(data = mod4, lat = lat, lon = lon)
+obs4 <- list(data = obs4, lat = lat, lon = lon)
+attr(exp4, 'class') <- 's2dv_cube'
+attr(obs4, 'class') <- 's2dv_cube'
+##############################################
+
+test_that("1. Input checks", {
+  # s2dv_cube
+  expect_error(
+    CST_Anomaly(exp = 1, obs = 1),
+    "Parameter 'exp' and 'obs' must be of the class 's2dv_cube', as output by CSTools::CST_Load."
+  )
+  # exp and obs
+  expect_error(
+    CST_Anomaly(exp = exp1, obs = exp1),
+    "One of the parameter 'exp' or 'obs' cannot be NULL."
+  )
+  expect_error(
+    CST_Anomaly(exp = exp2, obs = obs),
+    "Parameter 'exp' and 'obs' must have same dimension names in element 'data'."
+  )
+  # dim_anom
+  expect_warning(
+    CST_Anomaly(exp = exp, obs = obs, dim_anom = 3),
+    paste0("Parameter 'dim_anom' must be a character string and a numeric value will not be ",
+            "accepted in the next release. The corresponding dimension name is assigned.")
+  )
+  expect_error(
+    CST_Anomaly(exp = exp3, obs = obs3),
+    "The length of dimension 'dim_anom' in label 'data' of the parameter 'exp' and 'obs' must be greater than 1."
+  )
+  expect_error(
+    CST_Anomaly(exp4, obs4),
+    "Parameter 'dim_anom' is not found in 'exp' or in 'obs' dimension in element 'data'."
+  )
+  # cross and memb
+  expect_error(
+    CST_Anomaly(exp = exp, obs = obs, cross = 1),
+    "Parameters 'cross' and 'memb' must be logical."
+  )
+  expect_error(
+    CST_Anomaly(exp = exp, obs = obs, memb = 1),
+    "Parameters 'cross' and 'memb' must be logical."
+  )
+  # memb_dim
+  expect_error(
+    CST_Anomaly(exp = exp, obs = obs, memb_dim = 1),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    names(CST_Anomaly(exp4, obs4, dim_anom = 'sdates')),
+    "Parameter 'memb_dim' is not found in 'exp' or in 'obs' dimension."
+  )
+  # filter_span
+  expect_warning(
+    CST_Anomaly(exp = exp, obs = obs, filter_span = 'a'),
+    "Paramater 'filter_span' is not numeric and any filter is being applied."
+  )
+  # dat_dim
+  expect_error(
+    names(CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members')),
+    "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension in element 'data'. Set it as NULL if there is no dataset dimension."
+  )
+  # ftime_dim
+  expect_error(
+    CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'), memb = FALSE, filter_span = 2)$exp$data[1,1,1,,1,1],
+    "Parameter 'ftime_dim' is not found in 'exp' or in 'obs' dimension in element 'data'.",
+  )
+})
+
+##############################################
+test_that("2. Output checks: dat", {
+
+  expect_equal(
+    names(CST_Anomaly(exp, obs)),
+    c("exp", "obs")
+  )
+  expect_equal(
+    dim(CST_Anomaly(exp, obs)$exp$data),
+    c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
+  )
+  expect_equal(
+    CST_Anomaly(exp, obs)$exp$data[1,1,1,,1,1],
+    c(-0.64169277, 0.04953107, -0.59168037, 0.42660289, -0.77391490),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    CST_Anomaly(exp, obs)$obs$data[1,1,1,,1,1],
+    c(-0.83326562, -0.21035806, 1.17320132, 0.09760576, 0.28872829),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    CST_Anomaly(exp, obs, memb = FALSE)$exp$data[1,1,1,,1,1],
+    c(-0.9385105, 0.5140613, -0.3985370, 0.2916146, -1.0413568),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    CST_Anomaly(exp, obs, memb = FALSE, filter_span = 2)$exp$data[1,1,1,,1,1],
+    c(-0.8645582, 0.3478374, -0.3914569, 0.4555659, -1.1119619),
+    tolerance = 0.0001
+  )
+})
+##############################################
+test_that("3. Output checks: dat4", {
+  
+  expect_equal(
+    names(CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'))),
+    c('exp', 'obs')
+  )
+  expect_equal(
+    dim(CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'))$exp$data),
+    c(datasets = 2, members = 3, sdates = 4, ftimes = 5, lat = 6, lon = 7)
+  )
+  expect_equal(
+    CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'))$exp$data[1,1,1,,1,1],
+    c(-0.64169277, 0.04953107, -0.59168037, 0.42660289, -0.77391490),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'))$obs$data[1,1,1,,1,1],
+    c(-0.83326562, -0.21035806, 1.17320132, 0.09760576, 0.28872829),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'), memb = FALSE)$exp$data[1,1,1,,1,1],
+    c(-0.9385105, 0.5140613, -0.3985370, 0.2916146, -1.0413568),
+    tolerance = 0.0001
+  )
+  
+  expect_equal(
+    CST_Anomaly(exp4, obs4, dim_anom = 'sdates', memb_dim = 'members', dat_dim = c('datasets', 'members'), ftime_dim = 'ftimes', memb = FALSE, filter_span = 2)$exp$data[1,1,1,,1,1],
+    c(-0.8645582, 0.3478374, -0.3914569, 0.4555659, -1.1119619),
+    tolerance = 0.0001
+  )
+})
\ No newline at end of file
diff --git a/vignettes/WeatherRegimes_vignette.Rmd b/vignettes/WeatherRegimes_vignette.Rmd
index 47d095ce6d424a07a9ce8ee627b1e245d9366eef..788b25761eaac40a729cd90132b2d2abf5335a23 100644
--- a/vignettes/WeatherRegimes_vignette.Rmd
+++ b/vignettes/WeatherRegimes_vignette.Rmd
@@ -30,7 +30,7 @@ library(zeallot)
 The data employed in this example are described below.  
 - Sea level pressure (psl): this has been selected as the circulation variable, however other variables such as geopotential at 500 hPa can be also used.
 - Region: Euro-Atlantic domain [85.5ºW-45ºE; 27-81ºN].
-- Datasets: seasonal predictions from ECMWF System 4 ([**Molteni et al. 2011**] (https://www.ecmwf.int/sites/default/files/elibrary/2011/11209-new-ecmwf-seasonal-forecast-system-system-4.pdf)) and ERA-Interim reanalysis ([**Dee et al. 2011**] (https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.828)) as a reference dataset.
+- Datasets: seasonal predictions from ECMWF System 4 ([**Molteni et al. 2011**](https://www.ecmwf.int/sites/default/files/elibrary/2011/11209-new-ecmwf-seasonal-forecast-system-system-4.pdf)) and ERA-Interim reanalysis ([**Dee et al. 2011**](https://doi.org/10.1002/qj.828)) as a reference dataset.
 - Period: 1991-2010. Only 20 years have been selected for illustrative purposes, but the full hindcast period could be used for the analysis.  
 
 
@@ -75,7 +75,7 @@ The LOESS filter has been applied to the climatology to remove the short-term va
 
 ### 4- Weather regimes in observations 
 
-`CST_WeatherRegimes()` function is used to define the clusters based on the sea level pressure anomalies from ERA-Interim. This function is based on the [*kmeans function*] (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/kmeans.html)
+`CST_WeatherRegimes()` function is used to define the clusters based on the sea level pressure anomalies from ERA-Interim. This function is based on the [*kmeans function*](https://stat.ethz.ch/R-manual/R-devel/library/stats/html/kmeans.html)
 from the stats R package. In this example we have made different assumptions: four clusters (`ncenters=4`) will be produced and the Empirical orthogonal functions are not used to filter the data (`EOFS=FALSE`) just to take into account the extreme values. More details about the methodology can be found in Cortesi et al. 2018 (submitted).