diff --git a/modules/Anomalies/Anomalies.R b/modules/Anomalies/Anomalies.R
new file mode 100644
index 0000000000000000000000000000000000000000..2fe1da0e94467c2dfc0b48cfd29035b112a8b1fb
--- /dev/null
+++ b/modules/Anomalies/Anomalies.R
@@ -0,0 +1,96 @@
+source("modules/Anomalies/tmp/CST_Anomaly.R")
+
+# Compute the hcst, obs and fcst anomalies with or without cross-validation
+# and return them, along with the hcst and obs climatologies.
+
+compute_anomalies <- function(recipe, data) {
+
+  if (recipe$Analysis$Workflow$Anomalies$compute) {
+    if (recipe$Analysis$Workflow$Anomalies$cross_validation) {
+      cross <- TRUE
+      cross_msg <- "with"
+    } else {
+      cross <- FALSE
+      cross_msg <- "without"
+    }
+    original_dims <- dim(data$hcst$data)
+
+    # Compute anomalies
+    anom <- CST_Anomaly(data$hcst, data$obs,
+			cross = cross,
+			memb = TRUE,
+			memb_dim = 'ensemble',
+			dim_anom = 'syear',
+			dat_dim = c('dat', 'ensemble'),
+			ftime_dim = 'time')
+    # Reorder dims
+    anom$exp$data <- Reorder(anom$exp$data, names(original_dims))
+    anom$obs$data <- Reorder(anom$obs$data, names(original_dims))
+
+    # Save full fields
+    hcst_fullvalue <- data$hcst
+    obs_fullvalue <- data$obs
+    
+    # Hindcast climatology 
+
+    data$hcst <- anom$exp
+    data$obs <- anom$obs
+    remove(anom)
+    # Change variable metadata
+    # data$hcst$Variable$varName <- paste0(data$hcst$Variable$varName, "anomaly")
+    attr(data$hcst$Variable, "variable")$long_name <- 
+      paste(attr(data$hcst$Variable, "variable")$long_name, "anomaly")
+    # data$obs$Variable$varName <- paste0(data$obs$Variable$varName, "anomaly")
+    attr(data$obs$Variable, "variable")$long_name <-
+      paste(attr(data$obs$Variable, "variable")$long_name, "anomaly")
+
+    # Compute forecast anomaly field
+    if (!is.null(data$fcst)) {
+      # Compute hindcast climatology ensemble mean
+      clim <- s2dv::Clim(hcst_fullvalue$data, obs_fullvalue$data,
+          	         time_dim = "syear",
+        	         dat_dim = c("dat", "ensemble"),
+        	         memb = FALSE,
+        	         memb_dim = "ensemble",
+        	         ftime_dim = "time")
+      clim_hcst <- InsertDim(clim$clim_exp, posdim = 1, lendim = 1,
+                             name = "syear")
+      dims <- dim(clim_hcst)
+      clim_hcst <- rep(clim_hcst, dim(data$fcst$data)[['ensemble']])
+      dim(clim_hcst) <- c(dims, ensemble = dim(data$fcst$data)[['ensemble']])
+      # Get fcst anomalies
+      data$fcst$data <- data$fcst$data - clim_hcst
+      # Change metadata
+      # data$fcst$Variable$varName <- paste0(data$fcst$Variable$varName, "anomaly")
+      attr(data$fcst$Variable, "variable")$long_name <- 
+        paste(attr(data$fcst$Variable, "variable")$long_name, "anomaly")
+    }
+
+    info(recipe$Run$logger,
+  	 paste("The anomalies have been computed,", cross_msg,
+	       "cross-validation. The original full fields are returned as",
+	       "$hcst.full_val and $obs.full_val."))
+
+    info(recipe$Run$logger, "##### ANOMALIES COMPUTED SUCCESSFULLY #####")
+
+  } else {
+    warn(recipe$Run$logger, paste("The Anomalies module has been called, but",
+	 "recipe parameter Analysis:Variables:anomaly is set to FALSE.",
+	 "The full fields will be returned."))
+    hcst_fullvalue <- NULL
+    obs_fullvalue <- NULL
+    info(recipe$Run$logger, "##### ANOMALIES NOT COMPUTED #####")
+  }
+
+  ## TODO: Return fcst full value?
+
+  return(list(hcst = data$hcst, obs = data$obs, fcst = data$fcst,
+	      hcst.full_val = hcst_fullvalue, obs.full_val = obs_fullvalue))
+
+}
+
+
+
+
+
+
diff --git a/modules/Anomalies/tmp/CST_Anomaly.R b/modules/Anomalies/tmp/CST_Anomaly.R
new file mode 100644
index 0000000000000000000000000000000000000000..a84b6fc8538b03f4113b96aa8b3126189a0bdee9
--- /dev/null
+++ b/modules/Anomalies/tmp/CST_Anomaly.R
@@ -0,0 +1,241 @@
+#'Anomalies relative to a climatology along selected dimension with or without cross-validation
+#'
+#'@author Perez-Zanon Nuria, \email{nuria.perez@bsc.es}
+#'@author Pena Jesus, \email{jesus.pena@bsc.es}
+#'@description This function computes the anomalies relative to a climatology 
+#'computed along the selected dimension (usually starting dates or forecast 
+#'time) allowing the application or not of crossvalidated climatologies. The 
+#'computation is carried out independently for experimental and observational 
+#'data products.
+#'
+#'@param exp An object of class \code{s2dv_cube} as returned by \code{CST_Load} 
+#'  function, containing the seasonal forecast experiment data in the element 
+#'  named \code{$data}.
+#'@param obs An object of class \code{s2dv_cube} as returned by \code{CST_Load} 
+#'  function, containing the observed data in the element named \code{$data}.
+#'@param 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.
+#'
+#'@examples
+#'# Example 1:
+#'mod <- 1 : (2 * 3 * 4 * 5 * 6 * 7)
+#'dim(mod) <- c(dataset = 2, member = 3, sdate = 4, ftime = 5, lat = 6, lon = 7)
+#'obs <- 1 : (1 * 1 * 4 * 5 * 6 * 7)
+#'dim(obs) <- c(dataset = 1, member = 1, sdate = 4, ftime = 5, lat = 6, lon = 7)
+#'lon <- seq(0, 30, 5)
+#'lat <- seq(0, 25, 5)
+#'exp <- list(data = mod, lat = lat, lon = lon)
+#'obs <- list(data = obs, lat = lat, lon = lon)
+#'attr(exp, 'class') <- 's2dv_cube'
+#'attr(obs, 'class') <- 's2dv_cube'
+#'
+#'anom1 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = TRUE)
+#'anom2 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = TRUE)
+#'anom3 <- CST_Anomaly(exp = exp, obs = obs, cross = TRUE, memb = FALSE)
+#'anom4 <- CST_Anomaly(exp = exp, obs = obs, cross = FALSE, memb = FALSE)
+#'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, 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.")
+  }
+  # 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 'obs' will be used instead.")
+  }  
+  if (is.null(obs)) {
+    obs <- exp 
+    case_exp = 1
+    warning("Parameter 'obs' is not provided and 'exp' will be used instead.")
+  }
+  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)
+  # 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.")
+  }
+  # cross
+  if (!is.logical(cross) | !is.logical(memb) ) {
+    stop("Parameters 'cross' and 'memb' must be logical.")
+  }
+  if (length(cross) > 1 | length(memb) > 1 ) {
+    cross <- cross[1]
+    warning("Parameter 'cross' has length greater than 1 and only the first element",
+            "will be used.")
+  }
+  # memb
+  if (length(memb) > 1) {
+     memb <- memb[1]
+     warning("Parameter 'memb' has length greater than 1 and only the first element",
+             "will be used.")
+  }
+  # 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, 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, time_dim = dim_anom, memb_dim = memb_dim, memb = memb, dat_dim = dat_dim)
+    if (!is.null(filter_span)) {
+      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, 1,  dim_exp[memb_dim]) 
+      clim_obs <- InsertDim(tmp$clim_obs, 1,  dim_obs[memb_dim]) 
+    }
+    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 
+  }
+
+  exp$data <- ano$exp
+  obs$data <- ano$obs
+  
+  #  Outputs
+  # ~~~~~~~~~
+    if (case_obs == 1) {
+      return(obs)
+    } 
+    else if (case_exp == 1) {
+      return(exp)
+    }
+    else {
+    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)
+  output <- predict(loess_filt)
+  return(output)
+}
+
diff --git a/modules/Calibration/Calibration.R b/modules/Calibration/Calibration.R
index 59e5451a437d2ed414c6c7fd9060f300ed9bd029..899b12913bfade3e1b3955a8236b22fe387e33f1 100644
--- a/modules/Calibration/Calibration.R
+++ b/modules/Calibration/Calibration.R
@@ -4,19 +4,27 @@ calibrate_datasets <- function(recipe, data) {
   # recipe. If the forecast is not null, it calibrates it as well.
   #
   # data: list of s2dv_cube objects containing the hcst, obs and fcst.
+  #       Optionally, it may also have hcst.full_val and obs.full_val.
   # recipe: object obtained when passing the .yml recipe file to read_yaml()
   
   method <- tolower(recipe$Analysis$Workflow$Calibration$method)
 
   if (method == "raw") {
-    warn(recipe$Run$logger, "The Calibration module has been called, 
-	 but the calibration method in the recipe is 'raw'. 
-	 The hcst and fcst will not be calibrated.")
+    warn(recipe$Run$logger,
+	 paste("The Calibration module has been called, but the calibration",
+	       "method in the recipe is 'raw'. The hcst and fcst will not be",
+	       "calibrated."))
     fcst_calibrated <- data$fcst
     hcst_calibrated <- data$hcst
+    if (!is.null(data$hcst.full_val)) {
+      hcst_full_calibrated <- data$hcst.full_val
+    } else {
+      hcst_full_calibrated <- NULL
+    }
     CALIB_MSG <- "##### NO CALIBRATION PERFORMED #####"
 
   } else {
+    ## TODO: Calibrate full fields when present
     # Calibration function params
     mm <- recipe$Analysis$Datasets$Multimodel
     if (is.null(recipe$Analysis$ncores)) {
@@ -32,6 +40,7 @@ calibrate_datasets <- function(recipe, data) {
     
     CALIB_MSG <- "##### CALIBRATION COMPLETE #####"
     # Replicate observation array for the multi-model case
+    ## TODO: Implement for obs.full_val
     if (mm) {
       obs.mm <- obs$data
       for(dat in 1:(dim(data$hcst$data)['dat'][[1]]-1)) {
@@ -47,13 +56,12 @@ calibrate_datasets <- function(recipe, data) {
     
       CST_CALIB_METHODS <- c("bias", "evmos", "mse_min", "crps_min", "rpc-based")
       ## TODO: implement other calibration methods
-      ## TODO: Restructure the code?
       if (!(method %in% CST_CALIB_METHODS)) {
-	error(recipe$Run$logger, "Calibration method in the recipe is not
-	      available for monthly data.")
+	error(recipe$Run$logger, 
+	      paste("Calibration method in the recipe is not available for",
+	            "monthly data."))
         stop()
       } else {
-        ## Alba's version of CST_Calibration (pending merge) is being used
         # Calibrate the hindcast
         hcst_calibrated <- CST_Calibration(data$hcst, data$obs,
                                            cal.method = method,
@@ -66,8 +74,25 @@ calibrate_datasets <- function(recipe, data) {
                                            memb_dim = "ensemble",
                                            sdate_dim = "syear",
                                            ncores = ncores)
+	# In the case where anomalies have been computed, calibrate full values
+	if (!is.null(data$hcst.full_val)) {
+	  hcst_full_calibrated <- CST_Calibration(data$hcst.full_val,
+						  data$obs.full_val,
+						  cal.method = method,
+						  eval.method = "leave-one-out",
+						  multi.model = mm,
+						  na.fill = TRUE,
+						  na.rm = na.rm,
+						  apply_to = NULL,
+						  memb_dim = "ensemble",
+						  sdate_dim = "syear",
+						  ncores = ncores)
+	} else {
+	  hcst_full_calibrated <- NULL
+	}
+
+	# Calibrate the forecast
         if (!is.null(data$fcst)) {
-          # Calibrate the forecast
           fcst_calibrated <- CST_Calibration(data$hcst, data$obs, data$fcst,
                                              cal.method = method,
                                              eval.method = "leave-one-out",
@@ -86,9 +111,9 @@ calibrate_datasets <- function(recipe, data) {
     } else if (recipe$Analysis$Variables$freq == "daily_mean") {
       # Daily data calibration using Quantile Mapping
       if (!(method %in% c("qmap"))) {
-	error(recipe$Run$logger, "Calibration method in the recipe is not
-	      available for daily data. Only quantile mapping 'qmap is
-	      implemented.")
+	error(recipe$Run$logger,
+	      paste("Calibration method in the recipe is not available for",
+		    "daily data. Only quantile mapping 'qmap is implemented."))
         stop()
       }
       # Calibrate the hindcast
@@ -104,7 +129,21 @@ calibrate_datasets <- function(recipe, data) {
 					     wet.day = F)
       # Restore dimension order
       hcst_calibrated$data <- Reorder(hcst_calibrated$data, dim_order)
-  
+      # In the case where anomalies have been computed, calibrate full values
+      if (!is.null(data$hcst.full_val)) {
+	  hcst_full_calibrated <- CST_QuantileMapping(data$hcst.full_val,
+						      data$obs.full_val,
+						      exp_cor = NULL,
+						      sdate_dim = "syear",
+						      memb_dim = "ensemble",
+						      method = "QUANT",
+						      ncores = ncores,
+						      na.rm = na.rm,
+						      wet.day = F)
+      } else {
+	hcst_full_calibrated <- NULL
+      }
+
       if (!is.null(data$fcst)) {
         # Calibrate the forecast
         fcst_calibrated <- CST_QuantileMapping(data$hcst, data$obs,
@@ -124,5 +163,14 @@ calibrate_datasets <- function(recipe, data) {
     }
   }
   info(recipe$Run$logger, CALIB_MSG)
-  return(list(hcst = hcst_calibrated, fcst = fcst_calibrated))
+  ## TODO: Sort out returns
+  return_list <- list(hcst = hcst_calibrated,
+		      obs = data$obs,
+		      fcst = fcst_calibrated)
+  if (!is.null(hcst_full_calibrated)) {
+    return_list <- append(return_list,
+			  list(hcst.full_val = hcst_full_calibrated,
+			       obs.full_val = data$obs.full_val))
+  }
+  return(return_list)
 }
diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R
index 8d54d63def1acb94639e0bf7732a040a7f98b206..66a53451091a85836183e3074fca50d81ba78ed5 100644
--- a/modules/Loading/Loading.R
+++ b/modules/Loading/Loading.R
@@ -3,9 +3,9 @@ source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R")
 # Load required libraries/funs
 source("modules/Loading/dates2load.R")
 source("modules/Loading/check_latlon.R")
+## TODO: Move to prepare_outputs.R
 source("tools/libs.R")
 
-
 load_datasets <- function(recipe) {
 
   # -------------------------------------------
@@ -343,6 +343,7 @@ load_datasets <- function(recipe) {
     }
   }
 
+  # Compute anomalies if requested
   # Print a summary of the loaded data for the user, for each object
   if (recipe$Run$logger$threshold <= 2) {
     data_summary(hcst, recipe)
diff --git a/modules/Loading/testing_recipes/recipe_circular-sort-test.yml b/modules/Loading/testing_recipes/recipe_circular-sort-test.yml
deleted file mode 100644
index 700fd3b2a7fa6a41158c4fe51ea89f2ff32f639e..0000000000000000000000000000000000000000
--- a/modules/Loading/testing_recipes/recipe_circular-sort-test.yml
+++ /dev/null
@@ -1,43 +0,0 @@
-Description:
-  Author: V. Agudetse
-  Info: For testing the behavior of the loading module when loading data
-        that crosses the date line or the Greenwich meridian.
-Analysis:
-  Horizon: Seasonal
-  Variables:
-    name: tas
-    freq: monthly_mean
-  Datasets:
-    System:
-      name: system7c3s
-    Multimodel: False
-    Reference:
-      name: era5
-  Time:
-    sdate: '1101'
-    fcst_year: 
-    hcst_start: '1993'
-    hcst_end: '2003'
-    leadtimemin: 2
-    leadtimemax: 2
-  Region:
-    latmin: -10
-    latmax: 10
-    lonmin: 320
-    lonmax: 350
-  Regrid:
-    method: bilinear
-    type: to_system
-  Workflow:
-    Calibration:
-      method: mse_min
-    Skill:
-      metric: BSS90
-    Indicators:
-      index: no
-  Output_format: S2S4E
-Run:
-  Loglevel: INFO
-  Terminal: yes
-  output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/
-  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
diff --git a/modules/Loading/testing_recipes/recipe_decadal.yml b/modules/Loading/testing_recipes/recipe_decadal.yml
index de4f959138c0c448107032efb9211242e5b1370a..986578f7cc7a74e44604c83fb6080ada63637406 100644
--- a/modules/Loading/testing_recipes/recipe_decadal.yml
+++ b/modules/Loading/testing_recipes/recipe_decadal.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system #to_reference
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: bias
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_decadal_daily.yml b/modules/Loading/testing_recipes/recipe_decadal_daily.yml
index f362329e387a7e0db259cac2d220d8f62d3b9b1a..9d404bfa45da3d70620b5f22920067d16f12afe5 100644
--- a/modules/Loading/testing_recipes/recipe_decadal_daily.yml
+++ b/modules/Loading/testing_recipes/recipe_decadal_daily.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system #to_reference
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: qmap
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_decadal_monthly_2.yml b/modules/Loading/testing_recipes/recipe_decadal_monthly_2.yml
index f9130e16d3dc97f6509c6e7668754f678cddbdd7..38b25d42295e198d714956e11cdf8ffdf006ed12 100644
--- a/modules/Loading/testing_recipes/recipe_decadal_monthly_2.yml
+++ b/modules/Loading/testing_recipes/recipe_decadal_monthly_2.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system #to_reference
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: bias
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_era5land.yml b/modules/Loading/testing_recipes/recipe_era5land.yml
deleted file mode 100644
index c99906b68c37e65c19c98ab062947d6f2d078e37..0000000000000000000000000000000000000000
--- a/modules/Loading/testing_recipes/recipe_era5land.yml
+++ /dev/null
@@ -1,46 +0,0 @@
-Description:
-  Author: V. Agudetse
-
-Analysis:
-  Horizon: Seasonal
-  Variables:
-    name: tas
-    freq: monthly_mean
-  Datasets:
-    System:
-      name: system7c3s
-    Multimodel: False
-    Reference:
-      name: era5land
-  Time:
-    sdate: '1101'
-    fcst_year: '2020'
-    hcst_start: '1993'
-    hcst_end: '2016'
-    ftime_min: 1
-    ftime_max: 3
-  Region:
-    latmin: -10
-    latmax: 10
-    lonmin: 0
-    lonmax: 20
-  Regrid:
-    method: bilinear
-    type: to_system
-  Workflow:
-    Calibration:
-      method: mse_min
-    Skill:
-      metric: RPS RPSS CRPS CRPSS FRPSS BSS10 BSS90 EnsCorr Corr
-    Probabilities:
-      percentiles: [[1/4, 2/4, 3/4]]
-    Indicators:
-      index: no
-  ncores: 1
-  remove_NAs: yes
-  Output_format: S2S4E
-Run:
-  Loglevel: INFO
-  Terminal: yes
-  output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/
-  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
diff --git a/modules/Loading/testing_recipes/recipe_system2c3s-prlr-nofcst.yml b/modules/Loading/testing_recipes/recipe_system2c3s-prlr-nofcst.yml
deleted file mode 100644
index 3bfad08ea575ad8545371cd1a8384221ec369265..0000000000000000000000000000000000000000
--- a/modules/Loading/testing_recipes/recipe_system2c3s-prlr-nofcst.yml
+++ /dev/null
@@ -1,42 +0,0 @@
-Description:
-  Author: V. Agudetse
-
-Analysis:
-  Horizon: Seasonal
-  Variables:
-    name: prlr
-    freq: monthly_mean
-  Datasets:
-    System:
-      name: system2c3s
-    Multimodel: False
-    Reference:
-      name: era5
-  Time:
-    sdate: '0301'
-    fcst_year: #
-    hcst_start: '1993'
-    hcst_end: '2016'
-    ftime_min: 1
-    ftime_max: 1
-  Region:
-    latmin: -10
-    latmax: 10
-    lonmin: 0
-    lonmax: 20
-  Regrid:
-    method: bilinear
-    type: to_system
-  Workflow:
-    Calibration:
-      method: evmos
-    Skill:
-      metric: FRPS
-    Indicators:
-      index: no
-  Output_format: S2S4E
-Run:
-  Loglevel: INFO
-  Terminal: yes
-  output_dir: /esarchive/scratch/vagudets/repos/auto-s2s/out-logs/
-  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-rsds.yml b/modules/Loading/testing_recipes/recipe_system5c3s-rsds.yml
index ca52e8cc161ece1eb0a23e1be391df1b2318bf5f..94fc716cfcb32db9b361266fd8fc1724cad04aee 100644
--- a/modules/Loading/testing_recipes/recipe_system5c3s-rsds.yml
+++ b/modules/Loading/testing_recipes/recipe_system5c3s-rsds.yml
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: mse_min
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml
index 27cdccdc89c56ea8ef0e3322063573c2017be744..3a2bc72e7f022fbbfd792e6fbac5dae6838780db 100644
--- a/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml
+++ b/modules/Loading/testing_recipes/recipe_system5c3s-tas.yml
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: raw
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_system7c3s-prlr.yml b/modules/Loading/testing_recipes/recipe_system7c3s-prlr.yml
index 197c109c0e5dbb7f2419131cc2f4a328995dc34f..23b630b5b0194f280dfc3ce857ab52ea818d7ecf 100644
--- a/modules/Loading/testing_recipes/recipe_system7c3s-prlr.yml
+++ b/modules/Loading/testing_recipes/recipe_system7c3s-prlr.yml
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: mse_min
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_system7c3s-tas.yml b/modules/Loading/testing_recipes/recipe_system7c3s-tas.yml
index f83f91bb44f645ddbc84f6541f57fb82a09caa7b..df82c349eab2eb46f785ce5cdc0b0b46d4ed2601 100644
--- a/modules/Loading/testing_recipes/recipe_system7c3s-tas.yml
+++ b/modules/Loading/testing_recipes/recipe_system7c3s-tas.yml
@@ -16,9 +16,9 @@ Analysis:
     sdate: '1101'
     fcst_year: '2020'
     hcst_start: '1993'
-    hcst_end: '2016'
+    hcst_end: '2010'
     ftime_min: 1
-    ftime_max: 6
+    ftime_max: 2
   Region:
     latmin: -10
     latmax: 10
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: yes # yes/no, default yes
+      cross_validation: yes # yes/no, default yes
     Calibration:
       method: mse_min
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-reference.yml b/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-reference.yml
index 71e386fae2420e0bb4cdc01c241b40169d787f68..364d3dd60c6386eadb578fb383fb72c35e0c6ae7 100644
--- a/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-reference.yml
+++ b/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-reference.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear # Mandatory, str: Interpolation method. See docu.
     type: to_reference # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
   Workflow:
+    Anomalies:
+      compute: no # Whether to compute the anomalies and use them for skill metrics
+      cross_validation: # whether they should be computed in cross-validation
     Calibration:
       method: qmap # Mandatory, str: Calibration method. See docu.
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-system.yml b/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-system.yml
index 233c14eb4a74c7dad982ddbdf3e335c16cad601e..244a5654979ed02cf7cecd1ae462e264d6aeec32 100644
--- a/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-system.yml
+++ b/modules/Loading/testing_recipes/recipe_tas-daily-regrid-to-system.yml
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross_validation:
     Calibration:
       method: qmap
     Skill:
diff --git a/modules/Loading/testing_recipes/recipe_system7c3s-tas-specs.yml b/modules/Loading/testing_recipes/recipe_test_anomalies.yml
similarity index 72%
rename from modules/Loading/testing_recipes/recipe_system7c3s-tas-specs.yml
rename to modules/Loading/testing_recipes/recipe_test_anomalies.yml
index b1829d22e76b8e1e94a3218af3b91d46715600e1..cdf5e3ca885bff21d13e9cdeeed18618b48582af 100644
--- a/modules/Loading/testing_recipes/recipe_system7c3s-tas-specs.yml
+++ b/modules/Loading/testing_recipes/recipe_test_anomalies.yml
@@ -8,17 +8,17 @@ Analysis:
     freq: monthly_mean
   Datasets:
     System:
-      name: system7c3s
+      name: system5c3s
     Multimodel: False
     Reference:
       name: era5
   Time:
     sdate: '1101'
     fcst_year: '2020'
-    hcst_start: '1993'
-    hcst_end: '2016'
+    hcst_start: '1999'
+    hcst_end: '2010'
     ftime_min: 1
-    ftime_max: 3
+    ftime_max: 2
   Region:
     latmin: -10
     latmax: 10
@@ -29,13 +29,18 @@ Analysis:
     type: to_system
   Workflow:
     Calibration:
-      method: mse_min
+      method: raw
+    Anomalies:
+      compute: yes
+      cross_validation: yes
     Skill:
-      metric: CRPS CRPSS FCRPS FCRPSS FRPS_Specs
+      metric: RPS RPSS CRPS CRPSS BSS10 BSS90 EnsCorr mean_bias mean_bias_SS
     Probabilities:
       percentiles: [[1/3, 2/3], [1/10, 9/10]]
     Indicators:
       index: no
+  ncores: 7
+  remove_NAs: yes
   Output_format: S2S4E
 Run:
   Loglevel: INFO
diff --git a/modules/Saving/Saving.R b/modules/Saving/Saving.R
index 961baceebdfeec3de44635a034b1ba080d476893..1d72b30e18794f625375587184e79cc70fb0b55b 100644
--- a/modules/Saving/Saving.R
+++ b/modules/Saving/Saving.R
@@ -3,27 +3,28 @@
 source("modules/Saving/paths2save.R")
 
 save_data <- function(recipe, data,
-		      calibrated_data = NULL,
 		      skill_metrics = NULL,
 		      probabilities = NULL,
 		      archive = NULL) {
-
   # Wrapper for the saving functions.
   # recipe: The auto-s2s recipe
   # archive: The auto-s2s archive
   # data: output of load_datasets()
-  # calibrated_data: output of calibrate_datasets()
+  # data: output of calibrate_datasets()
   # skill_metrics: output of compute_skill_metrics()
   # probabilities: output of compute_probabilities()
   # mean_bias: output of compute_mean_bias()
 
   if (is.null(recipe)) {
-    stop("The 'recipe' parameter is mandatory.")
+    error(recipe$Run$logger, "The 'recipe' parameter is mandatory.")
+    stop()
   }
 
   if (is.null(data)) {
-    stop("The 'data' parameter is mandatory. It should be the output of",
-	 "load_datasets().")
+    error(recupe$Run$logger,
+	  paste("The 'data' parameter is mandatory. It should be a list",
+		"of at least two s2dv_cubes containing the hcst and obs."))
+    stop()
   }
   if (is.null(archive)) {
     if (tolower(recipe$Analysis$Horizon) == "seasonal") {
@@ -41,15 +42,13 @@ save_data <- function(recipe, data,
   dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
 
   # Export hindcast, forecast and observations onto outfile
-  if (!is.null(calibrated_data)) {
-    save_forecast(calibrated_data$hcst, recipe, dict, outdir,
-		  archive = archive, type = 'hcst')
-    if (!is.null(calibrated_data$fcst)) {
-      save_forecast(calibrated_data$fcst, recipe, dict, outdir,
-		    archive = archive, type = 'fcst')
-    }
-  save_observations(data$obs, recipe, dict, outdir, archive = archive)
+  save_forecast(data$hcst, recipe, dict, outdir, archive = archive,
+		type = 'hcst')
+  if (!is.null(data$fcst)) {
+    save_forecast(data$fcst, recipe, dict, outdir,
+		  archive = archive, type = 'fcst')
   }
+  save_observations(data$obs, recipe, dict, outdir, archive = archive)
 
   # Separate ensemble correlation from the rest of the metrics, as it has one
   # extra dimension "ensemble" and must be saved to a different file
@@ -106,7 +105,6 @@ get_global_attributes <- function(recipe, archive) {
 
 get_times <- function(store.freq, fcst.horizon, leadtimes, sdate, calendar) {
   # Generates time dimensions and the corresponding metadata.
-  ## TODO: Add calendar
   ## TODO: Subseasonal
 
   switch(fcst.horizon,
@@ -176,6 +174,7 @@ save_forecast <- function(data_cube,
   fcst.horizon <- tolower(recipe$Analysis$Horizon)
   store.freq <- recipe$Analysis$Variables$freq
   calendar <- archive$System[[global_attributes$system]]$calendar
+
 #  if (fcst.horizon == "seasonal") {
 #    calendar <- attr(data_cube$Variable, "variable")$dim$time$calendar
 #  } else {
@@ -186,12 +185,13 @@ save_forecast <- function(data_cube,
   dates <- as.PCICt(ClimProjDiags::Subset(data_cube$Dates$start, 'syear', 1),
 		    cal = calendar)
   if (fcst.horizon == 'decadal') {
-    # Method 1: Use the first date as init_date. But it may be better to use the real initialized date (ask users)
+    ## Method 1: Use the first date as init_date. But it may be better to use 
+    ## the real initialized date (ask users)
 #    init_date <- as.Date(data_cube$Dates$start[1], format = '%Y%m%d')
-    # Method 2: use initial month
+    ## Method 2: use initial month
     init_month <- archive$System[[recipe$Analysis$Datasets$System$name]]$initial_month
     if (type == 'hcst') {
-#PROBLEM for fcst!!!!!!!!!!!!
+    ## PROBLEM for fcst!!!!!!!!!!!!
       init_date <- as.PCICt(paste0(recipe$Analysis$Time$hcst_start, '-', 
                                    sprintf('%02d', init_month), '-01'),
 			    cal = calendar)
@@ -290,7 +290,8 @@ save_forecast <- function(data_cube,
       ArrayToNc(vars, outfile)
     }
   }
-  info(recipe$Run$logger, "##### FCST SAVED TO NETCDF FILE #####")
+  info(recipe$Run$logger, paste("#####", toupper(type),
+				"SAVED TO NETCDF FILE #####"))
 }
 
 
@@ -321,9 +322,10 @@ save_observations <- function(data_cube,
   dates <- as.PCICt(ClimProjDiags::Subset(data_cube$Dates$start, 'syear', 1),
 		    cal = calendar)
   if (fcst.horizon == 'decadal') {
-    # Method 1: Use the first date as init_date. But it may be better to use the real initialized date (ask users)
+    ## Method 1: Use the first date as init_date. But it may be better to use 
+    ## the real initialized date (ask users)
 #    init_date <- as.Date(data_cube$Dates$start[1], format = '%Y%m%d')
-    # Method 2: use initial month
+    ##  Method 2: use initial month
     init_month <- archive$System[[recipe$Analysis$Datasets$System$name]]$initial_month
     init_date <- as.PCICt(paste0(recipe$Analysis$Time$hcst_start, '-', 
                                  sprintf('%02d', init_month), '-01'),
@@ -457,6 +459,13 @@ save_metrics <- function(skill,
 
   # Add global and variable attributes
   global_attributes <- get_global_attributes(recipe, archive)
+  if (recipe$Analysis$Workflow$Anomalies$compute) {
+    global_attributes <- c(list(from_anomalies = "Yes"),
+			   global_attributes)
+  } else {
+    global_attributes <- c(list(from_anomalies = "No"),
+			   global_attributes)
+  }
   attr(skill[[1]], 'global_attrs') <- global_attributes
 
   for (i in 1:length(skill)) {
@@ -468,7 +477,7 @@ save_metrics <- function(skill,
       sdname <- paste0(metric, " region-aggregated metric")
       dims <- c('Country', 'time')
     } else {
-      sdname <- paste0(metric, " grid point metric")
+      sdname <- paste0(metric) #, " grid point metric")
       dims <- c(lalo, 'time')
     }
     metadata <- list(metric = list(name = metric,
@@ -564,6 +573,13 @@ save_corr <- function(skill,
 
   # Add global and variable attributes
   global_attributes <- get_global_attributes(recipe, archive)
+  if (recipe$Analysis$Workflow$Anomalies$compute) {
+    global_attributes <- c(global_attributes, 
+			   list(from_anomalies = "Yes"))
+  } else {
+    global_attributes <- c(global_attributes,
+			   list(from_anomalies = "No"))
+  }
   attr(skill[[1]], 'global_attrs') <- global_attributes
 
   for (i in 1:length(skill)) {
@@ -575,7 +591,7 @@ save_corr <- function(skill,
       sdname <- paste0(metric, " region-aggregated metric")
       dims <- c('Country', 'ensemble', 'time')
     } else {
-      sdname <- paste0(metric, " grid point metric") # formerly names(metric)
+      sdname <- paste0(metric) #, " grid point metric") # formerly names(metric)
       dims <- c(lalo, 'ensemble', 'time')
     }
     metadata <- list(metric = list(name = metric,
@@ -670,6 +686,13 @@ save_percentiles <- function(percentiles,
 
   # Add global and variable attributes
   global_attributes <- get_global_attributes(recipe, archive)
+  if (recipe$Analysis$Workflow$Anomalies$compute) {
+    global_attributes <- c(list(from_anomalies = "Yes"),
+			   global_attributes)
+  } else {
+    global_attributes <- c(list(from_anomalies = "No"),
+			   global_attributes)
+  }
   attr(percentiles[[1]], 'global_attrs') <- global_attributes
 
   for (i in 1:length(percentiles)) {
@@ -769,6 +792,14 @@ save_probabilities <- function(probs,
   variable <- data_cube$Variable$varName
   var.longname <- attr(data_cube$Variable, 'variable')$long_name
   global_attributes <- get_global_attributes(recipe, archive)
+  # Add anomaly computation to global attributes
+  if (recipe$Analysis$Workflow$Anomalies$compute) {
+    global_attributes <- c(list(from_anomalies = "Yes"),
+			   global_attributes)
+  } else {
+    global_attributes <- c(list(from_anomalies = "No"),
+			   global_attributes)
+  }
   fcst.horizon <- tolower(recipe$Analysis$Horizon)
   store.freq <- recipe$Analysis$Variables$freq
   calendar <- archive$System[[global_attributes$system]]$calendar
diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R
index 76b08942c805cf2a0e5d889cb3f84761076d7d72..7e43469787abf17e1d0718cc6d573dd75771c196 100644
--- a/modules/Skill/Skill.R
+++ b/modules/Skill/Skill.R
@@ -52,12 +52,32 @@ source("modules/Skill/tmp/RMSSS.R")
 #                    " running Skill module ", "\n",
 #                    " it can call ", metric_fun ))
 
-compute_skill_metrics <- function(recipe, exp, obs) {
-  # exp: s2dv_cube containing the hindcast
+# compute_skill_metrics <- function(recipe, data$hcst, obs,
+# 				  clim_data$hcst = NULL,
+# 				  clim_obs = NULL) {
+compute_skill_metrics <- function(recipe, data) {
+
+  # data$hcst: s2dv_cube containing the hindcast
   # obs: s2dv_cube containing the observations
   # recipe: auto-s2s recipe as provided by read_yaml
   
   ## TODO: Adapt time_dims to subseasonal case
+  ## TODO: Add dat_dim
+  ## TODO: Refine error messages
+  ## TODO: Add check to see if anomalies are provided (info inside s2dv_cube)
+
+#   if (recipe$Analysis$Workflow$Anomalies$compute) {
+#     if (is.null(clim_data$hcst) || is.null(clim_obs)) {
+#       warn(recipe$Run$logger, "Anomalies have been requested in the recipe,
+# 	   but the climatologies have not been provided in the 
+# 	   compute_skill_metrics call. Be aware that some metrics like the 
+# 	   Mean Bias may not be correct.")
+#     }
+#   } else {
+#     warn(recipe$Run$logger, "Anomaly computation was not requested in the
+# 	 recipe. Be aware that some metrics, such as the CRPSS may not be 
+# 	 correct.")
+#   }
   time_dim <- 'syear'
   memb_dim <- 'ensemble'
   metrics <- tolower(recipe$Analysis$Workflow$Skill$metric)
@@ -88,22 +108,31 @@ compute_skill_metrics <- function(recipe, exp, obs) {
     } 
     # Ranked Probability Score and Fair version
     if (metric %in% c('rps', 'frps')) {
-      skill <- RPS(exp$data, obs$data, time_dim = time_dim, memb_dim = memb_dim,
-          	   Fair = Fair, ncores = ncores)
+      skill <- RPS(data$hcst$data, data$obs$data,
+		   time_dim = time_dim,
+		   memb_dim = memb_dim,
+          	   Fair = Fair,
+		   ncores = ncores)
       skill <- .drop_dims(skill)
       skill_metrics[[ metric ]] <- skill
     # Ranked Probability Skill Score and Fair version
     } else if (metric %in% c('rpss', 'frpss')) {
-      skill <- RPSS(exp$data, obs$data, time_dim = time_dim, memb_dim = memb_dim,
-          	    Fair = Fair, ncores = ncores)
+      skill <- RPSS(data$hcst$data, data$obs$data,
+		    time_dim = time_dim,
+		    memb_dim = memb_dim,
+          	    Fair = Fair,
+		    ncores = ncores)
       skill <- lapply(skill, function(x) {
 		      .drop_dims(x)})
       skill_metrics[[ metric ]] <- skill$rpss
       skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign
     # Brier Skill Score - 10th percentile
     } else if (metric == 'bss10') {
-      skill <- RPSS(exp$data, obs$data, time_dim = time_dim, 
-		    memb_dim = memb_dim, prob_thresholds = 0.1, Fair = Fair,
+      skill <- RPSS(data$hcst$data, data$obs$data,
+		    time_dim = time_dim, 
+		    memb_dim = memb_dim,
+		    prob_thresholds = 0.1,
+		    Fair = Fair,
 		    ncores = ncores)
       skill <- lapply(skill, function(x) {
 		      .drop_dims(x)})
@@ -111,8 +140,11 @@ compute_skill_metrics <- function(recipe, exp, obs) {
       skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign
     # Brier Skill Score - 90th percentile
     } else if (metric == 'bss90') {
-      skill <- RPSS(exp$data, obs$data, time_dim = time_dim,
-		    memb_dim = memb_dim, prob_thresholds = 0.9, Fair = Fair,
+      skill <- RPSS(data$hcst$data, data$obs$data,
+		    time_dim = time_dim,
+		    memb_dim = memb_dim,
+		    prob_thresholds = 0.9,
+		    Fair = Fair,
 		    ncores = ncores)
       skill <- lapply(skill, function(x) {
 		      .drop_dims(x)})
@@ -120,28 +152,56 @@ compute_skill_metrics <- function(recipe, exp, obs) {
       skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign
     # CRPS and FCRPS
     } else if (metric %in% c('crps', 'fcrps')) {
-      skill <- CRPS(exp$data, obs$data, time_dim = time_dim,
-		    memb_dim = memb_dim, Fair = Fair, ncores = ncores)
+      skill <- CRPS(data$hcst$data, data$obs$data,
+		    time_dim = time_dim,
+		    memb_dim = memb_dim,
+		    Fair = Fair,
+		    ncores = ncores)
       skill <- .drop_dims(skill)
       skill_metrics[[ metric ]] <- skill
     # CRPSS and FCRPSS
     } else if (metric %in% c('crpss', 'fcrpss')) {
-      skill <- CRPSS(exp$data, obs$data, time_dim = time_dim,
-		     memb_dim = memb_dim, Fair = Fair, ncores = ncores)
+      skill <- CRPSS(data$hcst$data, data$obs$data,
+		     time_dim = time_dim,
+		     memb_dim = memb_dim,
+		     Fair = Fair,
+		     ncores = ncores)
       skill <- lapply(skill, function(x) {
 		      .drop_dims(x)})
       skill_metrics[[ metric ]] <- skill$crpss
       skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign
-   # Mean bias (climatology) 
-   } else if (metric == 'mean_bias') {
-      skill <- Bias(exp$data, obs$data, time_dim = time_dim,
-		    memb_dim = memb_dim, ncores = ncores)
+    # Mean bias (climatology) 
+    } else if (metric == 'mean_bias') {
+      ## TODO: Eliminate option to compute from anomalies
+      # Compute from full field
+      if ((!is.null(data$hcst.full_val)) && (!is.null(data$obs.full_val)) && 
+	  (recipe$Analysis$Workflow$Anomalies$compute)) {
+	skill <- Bias(data$hcst.full_val$data, data$obs.full_val$data,
+		      time_dim = time_dim,
+		      memb_dim = memb_dim,
+		      ncores = ncores)
+      } else {
+        skill <- Bias(data$hcst$data, data$obs$data,
+                      time_dim = time_dim,
+           	      memb_dim = memb_dim,
+		      ncores = ncores)
+      }
       skill <- .drop_dims(skill)
       skill_metrics[[ metric ]] <- skill
     # Mean bias skill score
     } else if (metric == 'mean_bias_ss') {
-      skill <- AbsBiasSS(exp$data, obs$data, time_dim = time_dim,
-		         memb_dim = memb_dim, ncores = ncores)
+      if ((!is.null(data$hcst.full_val)) && (!is.null(data$obs.full_val)) && 
+	  (recipe$Analysis$Workflow$Anomalies$compute)) {
+	skill <- AbsBiasSS(data$hcst.full_val$data, data$obs.full_val$data,
+		           time_dim = time_dim,
+		           memb_dim = memb_dim,
+		           ncores = ncores)
+      } else {
+        skill <- AbsBiasSS(data$hcst$data, data$obs$data,
+	  		   time_dim = time_dim,
+		           memb_dim = memb_dim,
+			   ncores = ncores)
+      }
       skill <- lapply(skill, function(x) {
 		      .drop_dims(x)})
       skill_metrics[[ metric ]] <- skill$biasSS
@@ -150,7 +210,7 @@ compute_skill_metrics <- function(recipe, exp, obs) {
     } else if (metric %in% c('enscorr', 'corr')) {
       ## TODO: Return significance
       ## TODO: Implement option for Kendall and Spearman methods?
-      skill <- s2dv::Corr(exp$data, obs$data, 
+      skill <- s2dv::Corr(data$hcst$data, data$obs$data, 
           		  dat_dim = 'dat', 
           		  time_dim = time_dim, 
           		  method = 'pearson', 
@@ -181,14 +241,14 @@ compute_skill_metrics <- function(recipe, exp, obs) {
       skill_metrics[[ paste0(metric, "_significance") ]] <- skill$sign
     } else if (metric == 'enssprerr') {
       # Remove ensemble dim from obs to avoid veriApply warning
-      obs_noensdim <- ClimProjDiags::Subset(obs$data, "ensemble", 1,
+      obs_noensdim <- ClimProjDiags::Subset(data$obs$data, "ensemble", 1,
 					    drop = "selected")
       capture.output(
       skill <- easyVerification::veriApply(verifun = 'EnsSprErr',
-					   fcst = exp$data,
+					   fcst = data$hcst$data,
 					   obs = obs_noensdim,
-					   tdim = which(names(dim(exp$data))==time_dim),
-					   ensdim = which(names(dim(exp$data))==memb_dim),
+					   tdim = which(names(dim(data$hcst$data))==time_dim),
+					   ensdim = which(names(dim(data$hcst$data))==memb_dim),
 					   na.rm = na.rm,
 					   ncpus = ncores)
       )
@@ -202,12 +262,11 @@ compute_skill_metrics <- function(recipe, exp, obs) {
       metric_name <- (strsplit(metric, "_"))[[1]][1] # Get metric name
       if (!(metric_name %in% c('frpss', 'frps', 'bss10', 'bss90', 'enscorr',
 			       'rpss'))) {
-	## TODO: Test this scenario
 	warn(recipe$Run$logger,
-	     "Some of the requested metrics are not available.")
+	     "Some of the requested SpecsVerification metrics are not available.")
       }
       capture.output(
-      skill <- Compute_verif_metrics(exp$data, obs$data,
+      skill <- Compute_verif_metrics(data$hcst$data, data$obs$data,
           			     skill_metrics = metric_name,
           			     verif.dims=c("syear", "sday", "sweek"),
           			     na.rm = na.rm,
@@ -226,6 +285,7 @@ compute_skill_metrics <- function(recipe, exp, obs) {
 }
 
 compute_probabilities <- function(recipe, data) {
+  ## TODO: Do hcst and fcst at the same time 
 
   if (is.null(recipe$Analysis$ncores)) {
     ncores <- 1
@@ -293,7 +353,7 @@ compute_probabilities <- function(recipe, data) {
   if (!("time" %in% names(dim(metric_array)))) {
     dim(metric_array) <- c("time" = 1, dim(metric_array))
   }
-  # If array has memb_exp (Corr case), change name to 'ensemble'
+  # If array has memb dim (Corr case), change name to 'ensemble'
   if ("exp_memb" %in% names(dim(metric_array))) {
     names(dim(metric_array))[which(names(dim(metric_array)) == 
 				   "exp_memb")] <- "ensemble"
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 2d50a1a0a4f8971087f29e64e7bda832ff63043c..4a4039a940ff536cdbd6b580ebad117b6be24d5a 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -8,12 +8,10 @@ source("modules/Visualization/tmp/PlotCombinedMap.R")
 
 plot_data <- function(recipe,
 		      data,
-		      calibrated_data = NULL,
 		      skill_metrics = NULL,
 		      probabilities = NULL,
 		      archive = NULL,
 		      significance = F) {
-
   # Try to produce and save several basic plots.
   # recipe: the auto-s2s recipe as read by read_yaml()
   # archive: the auto-s2s archive as read by read_yaml()
@@ -26,10 +24,11 @@ plot_data <- function(recipe,
   outdir <- paste0(get_dir(recipe), "/plots/")
   dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
 
-  if ((is.null(skill_metrics)) && (is.null(calibrated_data)) && 
-      (is.null(data$fcst))) {
-    stop("The Visualization module has been called, but there is no data ",
-	 "that can be plotted.")
+  if ((is.null(skill_metrics)) && (is.null(data$fcst))) { 
+    error(recipe$Run$logger, "The Visualization module has been called,
+	  but there is no fcst in 'data', and 'skill_metrics' is NULL
+	  so there is no data that can be plotted.")
+    stop()
   }
 
   if (is.null(archive)) {
@@ -38,7 +37,7 @@ plot_data <- function(recipe,
 				  "conf/archive.yml"))$archive
     } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
       archive <- read_yaml(paste0(recipe$Run$code_dir,
-				   "conf/archive_decadal.yml"))$archive
+				  "conf/archive_decadal.yml"))$archive
     }
   }
 
@@ -49,21 +48,12 @@ plot_data <- function(recipe,
   }
 
   # Plot forecast ensemble mean
-  if (!is.null(calibrated_data$fcst)) {
-    plot_ensemble_mean(recipe, archive, calibrated_data$fcst, outdir)
-  } else if (!is.null(data$fcst)) {
-    warn(recipe$Run$logger, "Only the uncalibrated forecast was provided. 
-	 Using this data to plot the forecast ensemble mean.")
+  if (!is.null(data$fcst)) {
     plot_ensemble_mean(recipe, archive, data$fcst, outdir)
   }
 
   # Plot Most Likely Terciles
-  if ((!is.null(probabilities)) && (!is.null(calibrated_data$fcst))) {
-    plot_most_likely_terciles(recipe, archive, calibrated_data$fcst,
-			      probabilities$percentiles, outdir)
-  } else if ((!is.null(probabilities)) && (!is.null(data$fcst))) {
-    warn(recipe$Run$logger, "Only the uncalibrated forecast was provided.
-	 Using this data to plot the most likely terciles.")
+  if ((!is.null(probabilities)) && (!is.null(data$fcst))) {
     plot_most_likely_terciles(recipe, archive, data$fcst,
 			      probabilities$percentiles, outdir)
   }
@@ -74,7 +64,9 @@ plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics,
 
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq == "daily_mean") {
-    stop("Visualization functions not yet implemented for daily data.")
+    error(recipe$Run$logger, "Visualization functions not yet implemented
+	  for daily data.")
+    stop()
   }
   # Abort if skill_metrics is not list
   if (!is.list(skill_metrics) || is.null(names(skill_metrics))) {
@@ -196,6 +188,7 @@ plot_skill_metrics <- function(recipe, archive, data_cube, skill_metrics,
 
 plot_ensemble_mean <- function(recipe, archive, fcst, outdir) {
 
+  ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq == "daily_mean") {
     stop("Visualization functions not yet implemented for daily data.")
@@ -273,6 +266,7 @@ plot_most_likely_terciles <- function(recipe, archive,
 				      percentiles, 
 				      outdir) {
 
+  ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq == "daily_mean") {
     stop("Visualization functions not yet implemented for daily data.")
diff --git a/modules/test_decadal.R b/modules/test_decadal.R
index 80304f978052889c19eccdf85ba9f1e99488dfe8..270e38005688c34a7f41128e73e08fdc42b8f0da 100644
--- a/modules/test_decadal.R
+++ b/modules/test_decadal.R
@@ -16,15 +16,15 @@ data <- load_datasets(recipe)
 calibrated_data <- calibrate_datasets(recipe, data)
 
 # Compute skill metrics
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 
 # Compute percentiles and probability bins
 probabilities <- compute_probabilities(recipe, calibrated_data$hcst)
 
 # Export all data to netCDF
-save_data(recipe, data, calibrated_data, skill_metrics, probabilities)
+save_data(recipe, calibrated_data, skill_metrics, probabilities)
 
 # Plot data
-plot_data(recipe, data, calibrated_data, skill_metrics, probabilities,
+plot_data(recipe, calibrated_data, skill_metrics, probabilities,
 	  significance = T)
 
diff --git a/modules/test_seasonal.R b/modules/test_seasonal.R
index fba75bfeefe3753ab331737c7d1e5361a730d554..9347ea0b7f78e7328188519e9cfdda86be91b072 100644
--- a/modules/test_seasonal.R
+++ b/modules/test_seasonal.R
@@ -1,5 +1,6 @@
 source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
+source("modules/Anomalies/Anomalies.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
 source("modules/Visualization/Visualization.R")
@@ -12,12 +13,14 @@ recipe <- prepare_outputs(recipe_file)
 data <- load_datasets(recipe)
 # Calibrate datasets
 calibrated_data <- calibrate_datasets(recipe, data)
-# Compute skill metrics
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+# Compute anomalies
+calibrated_data <- compute_anomalies(recipe, calibrated_data)
+## TODO: Turn arguments into (recipe, data)? 
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 # Compute percentiles and probability bins
 probabilities <- compute_probabilities(recipe, calibrated_data$hcst)
 # Export all data to netCDF
-save_data(recipe, data, calibrated_data, skill_metrics, probabilities)
+save_data(recipe, calibrated_data, skill_metrics, probabilities)
 # Plot data
-plot_data(recipe, data, calibrated_data, skill_metrics, probabilities,
+plot_data(recipe, calibrated_data, skill_metrics, probabilities,
 	  significance = T)
diff --git a/tests/recipes/recipe-decadal_daily_1.yml b/tests/recipes/recipe-decadal_daily_1.yml
index 90048009bf356ba79f87306efbd75c4914d78289..ab0fb9a625d7eb3d182678f138435e90d057672b 100644
--- a/tests/recipes/recipe-decadal_daily_1.yml
+++ b/tests/recipes/recipe-decadal_daily_1.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system #to_reference
   Workflow:
+    Anomalies:
+      compute: no
+      cross-validation: 
     Calibration:
       method: qmap
     Skill:
diff --git a/tests/recipes/recipe-decadal_monthly_1.yml b/tests/recipes/recipe-decadal_monthly_1.yml
index cedb07f002cfb46ed50e980ead642ba7e1884414..35b55b1a1985142d0bb6833ccd925da5142cfd3e 100644
--- a/tests/recipes/recipe-decadal_monthly_1.yml
+++ b/tests/recipes/recipe-decadal_monthly_1.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system #to_reference
   Workflow:
+    Anomalies:
+      compute: no
+      cross-validation: 
     Calibration:
       method: bias
     Skill:
diff --git a/tests/recipes/recipe-decadal_monthly_2.yml b/tests/recipes/recipe-decadal_monthly_2.yml
index 6c33d30270e41536dfa6a82afc31c412a23ca807..49824f42d326240039cf2eb40dfe5d5b66511fc4 100644
--- a/tests/recipes/recipe-decadal_monthly_2.yml
+++ b/tests/recipes/recipe-decadal_monthly_2.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross-validation:
     Calibration:
       method: raw
     Skill:
diff --git a/tests/recipes/recipe-decadal_monthly_3.yml b/tests/recipes/recipe-decadal_monthly_3.yml
index 87197af62ba40df07bed0da3243bed8d440678cf..1e2daa704725120f5b7450160d824e55a012e3e5 100644
--- a/tests/recipes/recipe-decadal_monthly_3.yml
+++ b/tests/recipes/recipe-decadal_monthly_3.yml
@@ -29,6 +29,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies: 
+      compute: no
+      cross-validation: 
     Calibration:
       method: 'evmos'
     Skill:
diff --git a/tests/recipes/recipe-seasonal_daily_1.yml b/tests/recipes/recipe-seasonal_daily_1.yml
index fc1bc58c49b2747999100a098aaf000e9282b926..637b53710504e14fba299a1feefa74b51554401a 100644
--- a/tests/recipes/recipe-seasonal_daily_1.yml
+++ b/tests/recipes/recipe-seasonal_daily_1.yml
@@ -28,6 +28,9 @@ Analysis:
     method: conservative
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross-validation:
     Calibration:
       method: qmap
     Skill:
diff --git a/tests/recipes/recipe-seasonal_monthly_1.yml b/tests/recipes/recipe-seasonal_monthly_1.yml
index 0e1b6f4b023c39e01841ffa98353e16a0bca690d..e75ccad5e6e12941719f11845d35bf46d6912ef2 100644
--- a/tests/recipes/recipe-seasonal_monthly_1.yml
+++ b/tests/recipes/recipe-seasonal_monthly_1.yml
@@ -28,6 +28,9 @@ Analysis:
     method: bilinear
     type: to_system
   Workflow:
+    Anomalies:
+      compute: no
+      cross-validation: 
     Calibration:
       method: mse_min
     Skill:
diff --git a/tests/testthat/test-decadal_monthly_1.R b/tests/testthat/test-decadal_monthly_1.R
index 5cf1922e73e1ee4bc1a8c46f9998cc64ecbc145b..184f6f43b22af2ce1a2ba0edf190eff959b71f5a 100644
--- a/tests/testthat/test-decadal_monthly_1.R
+++ b/tests/testthat/test-decadal_monthly_1.R
@@ -24,7 +24,7 @@ suppressWarnings({invisible(capture.output(
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
 probs <- compute_probabilities(recipe, calibrated_data$hcst)
@@ -32,15 +32,14 @@ probs <- compute_probabilities(recipe, calibrated_data$hcst)
 
 # Saving
 suppressWarnings({invisible(capture.output(
-save_data(recipe = recipe, data = data, calibrated_data = calibrated_data, 
+save_data(recipe = recipe, data =  calibrated_data, 
           skill_metrics = skill_metrics, probabilities = probs, archive = archive)
 ))})
 
 # Plotting
 suppressWarnings({invisible(capture.output(
-plot_data(recipe = recipe, archive = archive, data = data,
-	  calibrated_data = calibrated_data, skill_metrics = skill_metrics,
-	  probabilities = probs, significance = T)
+plot_data(recipe = recipe, archive = archive, data = calibrated_data,
+	  skill_metrics = skill_metrics, probabilities = probs, significance = T)
 ))})
 
 
@@ -137,7 +136,7 @@ TRUE
 )
 expect_equal(
 names(calibrated_data),
-c("hcst", "fcst")
+c("hcst", "obs", "fcst")
 )
 expect_equal(
 class(calibrated_data$hcst),
diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R
index 4dd72ebf13c632479bbd34772a54509805b7c057..6549ce0e763b33035ac12711c91d1627e3a25906 100644
--- a/tests/testthat/test-decadal_monthly_2.R
+++ b/tests/testthat/test-decadal_monthly_2.R
@@ -22,7 +22,7 @@ suppressWarnings({invisible(capture.output(
 
 # Compute skill metrics
 suppressMessages({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
 probs <- compute_probabilities(recipe, calibrated_data$hcst)
@@ -132,12 +132,13 @@ test_that("2. Calibration", {
 
 expect_equal(
 names(calibrated_data),
-c("hcst", "fcst")
-)
-expect_equal(
-calibrated_data,
-data[1:2]
+c("hcst", "obs", "fcst")
 )
+## TODO: Ask An-Chi about this test
+# expect_equal(
+# calibrated_data,
+# data[1:2]
+# )
 
 })
 
diff --git a/tests/testthat/test-decadal_monthly_3.R b/tests/testthat/test-decadal_monthly_3.R
index 7535e8dc6fcdedeb2eb7a69fa1f13e917e6178bc..70e98160132d5d515ed75a2c3fae8962bed385ac 100644
--- a/tests/testthat/test-decadal_monthly_3.R
+++ b/tests/testthat/test-decadal_monthly_3.R
@@ -23,7 +23,7 @@ suppressWarnings({invisible(capture.output(
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
 probs <- compute_probabilities(recipe, calibrated_data$hcst)
@@ -110,7 +110,7 @@ test_that("2. Calibration", {
 
 expect_equal(
 names(calibrated_data),
-c("hcst", "fcst")
+c("hcst", "obs", "fcst")
 )
 expect_equal(
 as.vector(aperm(drop(calibrated_data$hcst$data), c(5, 1:4))[3, , 2, 2, 2]),
diff --git a/tests/testthat/test-seasonal_daily.R b/tests/testthat/test-seasonal_daily.R
index 5b771d77fd4bf22eb74626ec7b2170602234dd0b..ddcca22fd93750647b02ecfc5290591edfc5167d 100644
--- a/tests/testthat/test-seasonal_daily.R
+++ b/tests/testthat/test-seasonal_daily.R
@@ -19,7 +19,7 @@ calibrated_data <- calibrate_datasets(recipe, data)
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 ))})
 
 test_that("1. Loading", {
@@ -106,7 +106,7 @@ TRUE
 )
 expect_equal(
 names(calibrated_data),
-c("hcst", "fcst")
+c("hcst", "obs", "fcst")
 )
 expect_equal(
 class(calibrated_data$hcst),
@@ -163,3 +163,5 @@ c(0.7509920, 0.6514916, 0.5118371),
 tolerance=0.0001
 )
 })
+
+unlink(recipe$Run$output_dir)
diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R
index 86feedfbb14eb550cdf8cb26ee00283964d3df4d..d46d20dcccbf4f17a7279b9437fea62c471a840e 100644
--- a/tests/testthat/test-seasonal_monthly.R
+++ b/tests/testthat/test-seasonal_monthly.R
@@ -22,7 +22,7 @@ calibrated_data <- calibrate_datasets(recipe, data)
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data$hcst, data$obs)
+skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
 ))})
 
 suppressWarnings({invisible(capture.output(
@@ -31,13 +31,13 @@ probs <- compute_probabilities(recipe, calibrated_data$hcst)
 
 # Saving
 suppressWarnings({invisible(capture.output(
-save_data(recipe = recipe, data = data, calibrated_data = calibrated_data,
+save_data(recipe = recipe, data = calibrated_data,
 	  skill_metrics = skill_metrics, probabilities = probs)
 ))})
 
 # Plotting
 suppressWarnings({invisible(capture.output(
-plot_data(recipe = recipe,  data = data, calibrated_data = calibrated_data,
+plot_data(recipe = recipe,  data = calibrated_data,
 	  skill_metrics = skill_metrics, probabilities = probs,
 	  significance = T)
 ))})
@@ -133,7 +133,7 @@ TRUE
 )
 expect_equal(
 names(calibrated_data),
-c("hcst", "fcst")
+c("hcst", "obs", "fcst")
 )
 expect_equal(
 class(calibrated_data$hcst),