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/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/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