diff --git a/example_scripts/test_decadal.R b/example_scripts/test_decadal.R
index 900b26223d1263fdca0032378cb878509485b20b..12daa540012d05d876c1da2fe1f142c11b58d66a 100644
--- a/example_scripts/test_decadal.R
+++ b/example_scripts/test_decadal.R
@@ -4,7 +4,7 @@
 ##              do not modify.
 ###############################################################################
 
-source("modules/Loading/Loading_decadal.R")
+source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
@@ -15,21 +15,18 @@ recipe <- prepare_outputs(recipe_file)
 # archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$archive
 
 # Load datasets
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 
 # Calibrate datasets
-calibrated_data <- calibrate_datasets(recipe, data)
+calibrated_data <- Calibration(recipe, data)
 
 # Compute skill metrics
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 
 # Compute percentiles and probability bins
-probabilities <- compute_probabilities(recipe, calibrated_data)
-
-# Export all data to netCDF
-save_data(recipe, calibrated_data, skill_metrics, probabilities)
+probabilities <- Probabilities(recipe, calibrated_data)
 
 # Plot data
-plot_data(recipe, calibrated_data, skill_metrics, probabilities,
-	  significance = T)
+Visualization(recipe, calibrated_data, skill_metrics, probabilities,
+	      significance = T)
 
diff --git a/example_scripts/test_parallel_workflow.R b/example_scripts/test_parallel_workflow.R
index 8cf79c25dc2fe386bf7e129b16317d9c428ebf99..5f0265135cde51a0cdcc41bb2133e4d684da1ef3 100644
--- a/example_scripts/test_parallel_workflow.R
+++ b/example_scripts/test_parallel_workflow.R
@@ -9,15 +9,15 @@ args = commandArgs(trailingOnly = TRUE)
 recipe_file <- args[1]
 recipe <- read_atomic_recipe(recipe_file)
 # Load datasets
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 # Calibrate datasets
-data <- calibrate_datasets(recipe, data)
+data <- Calibration(recipe, data)
 # Compute anomalies
-data <- compute_anomalies(recipe, data)
+data <- Anomalies(recipe, data)
 # Compute skill metrics 
-skill_metrics <- compute_skill_metrics(recipe, data)
+skill_metrics <- Skill(recipe, data)
 # Compute percentiles and probability bins
-probabilities <- compute_probabilities(recipe, data)
+probabilities <- Probabilities(recipe, data)
 # Plot data
-plot_data(recipe, data, skill_metrics, probabilities,
-	  significance = T)
+Visualization(recipe, data, skill_metrics, probabilities,
+	      significance = T)
diff --git a/example_scripts/test_seasonal.R b/example_scripts/test_seasonal.R
index 9734aeafb6629664e18a2eeb907e0d42bec5a78a..2c7d673a68ebbb6d7116e5845da78a955b28dda4 100644
--- a/example_scripts/test_seasonal.R
+++ b/example_scripts/test_seasonal.R
@@ -18,19 +18,16 @@ recipe_file <- "recipes/atomic_recipes/recipe_test_multivar.yml"
 recipe <- prepare_outputs(recipe_file)
 
 # Load datasets
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 # Change units
 data <- Units(recipe, data)
 # Calibrate datasets
-data <- calibrate_datasets(recipe, data)
+data <- Calibration(recipe, data)
 # Compute anomalies
-data <- compute_anomalies(recipe, data)
+data <- Anomalies(recipe, data)
 # Compute skill metrics 
-skill_metrics <- compute_skill_metrics(recipe, data)
+skill_metrics <- Skill(recipe, data)
 # Compute percentiles and probability bins
-probabilities <- compute_probabilities(recipe, data)
-# Export all data to netCDF
-## TODO: Fix plotting
-# save_data(recipe, data, skill_metrics, probabilities)
+probabilities <- Probabilities(recipe, data)
 # Plot data
-plot_data(recipe, data, skill_metrics, probabilities, significance = T)
+Visualization(recipe, data, skill_metrics, probabilities, significance = T)
diff --git a/modules/Anomalies/Anomalies.R b/modules/Anomalies/Anomalies.R
index 470d4706b68f4cfc2d0349e824a8db36fef86d66..386317fbea79fbb637aa82dc6a1f514861b76ec8 100644
--- a/modules/Anomalies/Anomalies.R
+++ b/modules/Anomalies/Anomalies.R
@@ -1,7 +1,10 @@
+## TODO: Remove in the next release
+source("modules/Anomalies/compute_anomalies.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) {
+Anomalies <- function(recipe, data) {
 
   if (is.null(recipe$Analysis$Workflow$Anomalies$compute)) {
     error(recipe$Run$logger, 
diff --git a/modules/Anomalies/compute_anomalies.R b/modules/Anomalies/compute_anomalies.R
new file mode 100644
index 0000000000000000000000000000000000000000..2ef36a34086a9bc270f09740a337438370a43950
--- /dev/null
+++ b/modules/Anomalies/compute_anomalies.R
@@ -0,0 +1,7 @@
+compute_anomalies <- function(recipe, data) {
+  warning(paste0("The function compute_anomalies() has been renamed to: ",
+                 "'Anomalies()'. The name 'compute_anomalies()' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Anomalies(recipe, data))
+}
diff --git a/modules/Calibration/Calibration.R b/modules/Calibration/Calibration.R
index 5811fa0578fdf4b7ec52e668cb7c8c9d01f23f6c..171b22cf39eeecf8edde3843277a7556db721e94 100644
--- a/modules/Calibration/Calibration.R
+++ b/modules/Calibration/Calibration.R
@@ -1,5 +1,7 @@
+## TODO: Remove in the next release
+source("modules/Calibration/calibrate_datasets.R")
 
-calibrate_datasets <- function(recipe, data) {
+Calibration <- function(recipe, data) {
   # Function that calibrates the hindcast using the method stated in the
   # recipe. If the forecast is not null, it calibrates it as well.
   #
diff --git a/modules/Calibration/calibrate_datasets.R b/modules/Calibration/calibrate_datasets.R
new file mode 100644
index 0000000000000000000000000000000000000000..8264992f949788b40e6bd6d5bfd5471b1536f1e5
--- /dev/null
+++ b/modules/Calibration/calibrate_datasets.R
@@ -0,0 +1,7 @@
+calibrate_datasets <- function(recipe, data) {
+  warning(paste0("The function calibrate_datasets() has been renamed to: ",
+                 "'Calibration()'. The name 'calibrate_datasets' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Calibration(recipe, data))
+}
diff --git a/modules/Loading/Loading.R b/modules/Loading/Loading.R
index 8cc46dc7f5e3fd25b82d5846eca2254b997962b5..cf97efc20d55a9e82e0825284163aad01136fe0c 100644
--- a/modules/Loading/Loading.R
+++ b/modules/Loading/Loading.R
@@ -1,452 +1,19 @@
-## TODO: remove paths to personal scratchs 
-source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R")
-# Load required libraries/funs
-source("modules/Loading/R/dates2load.R")
-source("modules/Loading/R/get_timeidx.R")
-source("modules/Loading/R/check_latlon.R")
-## TODO: Move to prepare_outputs.R
 source("tools/libs.R")
+## TODO: Remove with the next release
+source("modules/Loading/load_datasets.R")
 
-load_datasets <- function(recipe) {
-  
-  # -------------------------------------------
-  # Set params -----------------------------------------
-  
-  hcst.inityear <- recipe$Analysis$Time$hcst_start
-  hcst.endyear <- recipe$Analysis$Time$hcst_end
-  lats.min <- recipe$Analysis$Region$latmin
-  lats.max <- recipe$Analysis$Region$latmax
-  lons.min <- recipe$Analysis$Region$lonmin
-  lons.max <- recipe$Analysis$Region$lonmax
-  ref.name <- recipe$Analysis$Datasets$Reference$name
-  exp.name <- recipe$Analysis$Datasets$System$name
-  
-  variable <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]]
-  store.freq <- recipe$Analysis$Variables$freq
-  
-  # get sdates array
-  ## LOGGER: Change dates2load to extract logger from recipe?
-  sdates <- dates2load(recipe, recipe$Run$logger)
-  
-  idxs <- NULL
-  idxs$hcst <- get_timeidx(sdates$hcst,
-                           recipe$Analysis$Time$ftime_min,
-                           recipe$Analysis$Time$ftime_max,
-                           time_freq=store.freq)
-  
-  if (!(is.null(sdates$fcst))) {
-    idxs$fcst <- get_timeidx(sdates$fcst,
-                             recipe$Analysis$Time$ftime_min,
-                             recipe$Analysis$Time$ftime_max,
-                             time_freq=store.freq)
-  }
-  
-  ## TODO: Examine this verifications part, verify if it's necessary
-  # stream <- verifications$stream
-  # sdates <- verifications$fcst.sdate
-  
-  ## TODO: define fcst.name
-  ##fcst.name <- recipe$Analysis$Datasets$System[[sys]]$name
-  
-  # get esarchive datasets dict:
-  ## TODO: Adapt to 'filesystem' option in recipe
-  archive <- read_yaml("conf/archive.yml")$esarchive
-  exp_descrip <- archive$System[[exp.name]]
-  
-  freq.hcst <- unlist(exp_descrip[[store.freq]][variable[1]])
-  reference_descrip <- archive$Reference[[ref.name]]
-  freq.obs <- unlist(reference_descrip[[store.freq]][variable[1]])
-  obs.dir <- reference_descrip$src
-  fcst.dir <- exp_descrip$src
-  hcst.dir <- exp_descrip$src
-  fcst.nmember <- exp_descrip$nmember$fcst
-  hcst.nmember <- exp_descrip$nmember$hcst
-  
-  ## TODO: it is necessary?
-  ##if ("accum" %in% names(reference_descrip)) {
-  ##  accum <- unlist(reference_descrip$accum[store.freq][[1]])
-  ##} else {
-  ##  accum <- FALSE
-  ##}
-
-  var_dir_obs <- reference_descrip[[store.freq]][variable]
-  var_dir_exp <- exp_descrip[[store.freq]][variable]
-  
-  # -----------
-  obs.path <- paste0(archive$src,
-                     obs.dir, store.freq, "/$var$", "$var_dir$",
-                     "/$var$_$file_date$.nc")
-  
-  hcst.path <- paste0(archive$src,
-                      hcst.dir, store.freq, "/$var$", "$var_dir$",
-                      "$var$_$file_date$.nc")
-  
-  fcst.path <- paste0(archive$src,
-                      hcst.dir, store.freq, "/$var$", "$var_dir$",
-                      "/$var$_$file_date$.nc")
-  
-  # Define regrid parameters:
-  #-------------------------------------------------------------------
-  regrid_params <- get_regrid_params(recipe, archive)
-  
-  # Longitude circular sort and latitude check
-  #-------------------------------------------------------------------
-  circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max)
-  
-  if (recipe$Analysis$Variables$freq == "monthly_mean"){
-    split_multiselected_dims = TRUE
-  } else {
-    split_multiselected_dims = FALSE
-  }
-  
-  # Load hindcast 
-  #-------------------------------------------------------------------
-  hcst <- Start(dat = hcst.path,
-                var = variable,
-                var_dir = var_dir_exp,
-                file_date = sdates$hcst,
-                time = idxs$hcst,
-                var_dir_depends = 'var',
-                latitude = values(list(lats.min, lats.max)),
-                latitude_reorder = Sort(),
-                longitude = values(list(lons.min, lons.max)),
-                longitude_reorder = circularsort,
-                transform = regrid_params$fcst.transform,
-                transform_params = list(grid = regrid_params$fcst.gridtype,
-                                        method = regrid_params$fcst.gridmethod),
-                transform_vars = c('latitude', 'longitude'),
-                synonims = list(latitude = c('lat', 'latitude'),
-                                longitude = c('lon', 'longitude'),
-                                ensemble = c('member', 'ensemble')),
-                ensemble = indices(1:hcst.nmember),
-                metadata_dims = 'var', # change to just 'var'?
-                return_vars = list(latitude = 'dat',
-                                   longitude = 'dat',
-                                   time = 'file_date'),
-                split_multiselected_dims = split_multiselected_dims,
-                retrieve = TRUE)
-  
-  # Remove var_dir dimension
-  if ("var_dir" %in% names(dim(hcst))) {
-    hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected")
-  }
-  
-  if (store.freq %in% c("daily_mean", "daily")) {
-    # Adjusts dims for daily case, could be removed if startR allows
-    # multidim split
-    names(dim(hcst))[which(names(dim(hcst)) == 'file_date')] <- "syear"
-    default_dims <- c(dat = 1, var = 1, sday = 1, 
-                      sweek = 1, syear = 1, time = 1, 
-                      latitude = 1, longitude = 1, ensemble = 1)
-    default_dims[names(dim(hcst))] <- dim(hcst)
-    dim(hcst) <- default_dims
-    # Change time attribute dimensions
-    default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-    names(dim(attr(hcst, "Variables")$common$time))[which(names(
-      dim(attr(hcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-    default_time_dims[names(dim(attr(hcst, "Variables")$common$time))] <- 
-      dim(attr(hcst, "Variables")$common$time)
-    dim(attr(hcst, "Variables")$common$time) <- default_time_dims
-  }
-  
-  # Convert hcst to s2dv_cube object
-  ## TODO: Give correct dimensions to $Dates
-  ## (sday, sweek, syear instead of file_date)
-  hcst <- as.s2dv_cube(hcst)
-  # Adjust dates for models where the time stamp goes into the next month
-  if (recipe$Analysis$Variables$freq == "monthly_mean") {
-    hcst$attrs$Dates[] <- hcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag)
-  }
-  
-  # Load forecast
-  #-------------------------------------------------------------------
-  if (!is.null(recipe$Analysis$Time$fcst_year)) {
-    # the call uses file_date instead of fcst_syear so that it can work
-    # with the daily case and the current version of startR not allowing
-    # multiple dims split
-    
-    fcst <- Start(dat = fcst.path,
-                  var = variable,
-                  var_dir = var_dir_exp,
-                  var_dir_depends = 'var',
-                  file_date = sdates$fcst,
-                  time = idxs$fcst,
-                  latitude = values(list(lats.min, lats.max)),
-                  latitude_reorder = Sort(),
-                  longitude = values(list(lons.min, lons.max)),
-                  longitude_reorder = circularsort,
-                  transform = regrid_params$fcst.transform,
-                  transform_params = list(grid = regrid_params$fcst.gridtype,
-                                          method = regrid_params$fcst.gridmethod),
-                  transform_vars = c('latitude', 'longitude'),
-                  synonims = list(latitude = c('lat', 'latitude'),
-                                  longitude = c('lon', 'longitude'),
-                                  ensemble = c('member', 'ensemble')),
-                  ensemble = indices(1:fcst.nmember),
-                  metadata_dims = 'var',
-                  return_vars = list(latitude = 'dat',
-                                     longitude = 'dat',
-                                     time = 'file_date'),
-                  split_multiselected_dims = split_multiselected_dims,
-                  retrieve = TRUE)
-    
-    if ("var_dir" %in% names(dim(fcst))) {
-      fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected")
-    }
-    
-    if (store.freq %in% c("daily_mean", "daily")) {
-      # Adjusts dims for daily case, could be removed if startR allows
-      # multidim split
-      names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear"
-      default_dims <- c(dat = 1, var = 1, sday = 1,
-                        sweek = 1, syear = 1, time = 1,
-                        latitude = 1, longitude = 1, ensemble = 1)
-      default_dims[names(dim(fcst))] <- dim(fcst)
-      dim(fcst) <- default_dims
-      # Change time attribute dimensions
-      default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-      names(dim(attr(fcst, "Variables")$common$time))[which(names(
-        dim(attr(fcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-      default_time_dims[names(dim(attr(fcst, "Variables")$common$time))] <- 
-        dim(attr(fcst, "Variables")$common$time)
-      dim(attr(fcst, "Variables")$common$time) <- default_time_dims
-    }
-    
-    # Convert fcst to s2dv_cube
-    fcst <- as.s2dv_cube(fcst)
-    # Adjust dates for models where the time stamp goes into the next month
-    if (recipe$Analysis$Variables$freq == "monthly_mean") {
-      fcst$attrs$Dates[] <- 
-        fcst$attrs$Dates - seconds(exp_descrip$time_stamp_lag)
-    }
-    
+Loading <- function(recipe) {
+  # Source correct function depending on time horizon
+  ## TODO: Add condition for GRIB files
+  time_horizon <- tolower(recipe$Analysis$Horizon)
+  if (time_horizon == "seasonal") {
+    source("modules/Loading/R/load_seasonal.R")
+    data <- load_seasonal(recipe)
+  } else if (time_horizon == "decadal") {
+    source("modules/Loading/R/load_decadal.R")
+    data <- load_decadal(recipe)
   } else {
-    fcst <- NULL
-  }
-  
-  # Load reference
-  #-------------------------------------------------------------------
-  
-  # Obtain dates and date dimensions from the loaded hcst data to make sure
-  # the corresponding observations are loaded correctly.
-  dates <- hcst$attrs$Dates
-  dim(dates) <- hcst$dims[c("sday", "sweek", "syear", "time")]
-  
-  # Separate Start() call for monthly vs daily data
-  if (store.freq == "monthly_mean") {
-    
-    dates_file <- format(as.Date(dates, '%Y%m%d'), "%Y%m")
-    dim(dates_file) <- dim(dates)
-    
-    obs <- Start(dat = obs.path,
-                 var = variable,
-                 var_dir = var_dir_obs,
-                 var_dir_depends = 'var',
-                 file_date = dates_file,
-                 latitude = values(list(lats.min, lats.max)),
-                 latitude_reorder = Sort(),
-                 longitude = values(list(lons.min, lons.max)),
-                 longitude_reorder = circularsort,
-                 transform = regrid_params$obs.transform,
-                 transform_params = list(grid = regrid_params$obs.gridtype,
-                                         method = regrid_params$obs.gridmethod),
-                 transform_vars = c('latitude', 'longitude'),
-                 synonims = list(latitude = c('lat','latitude'),
-                                 longitude = c('lon','longitude')),
-                 metadata_dims = 'var',
-                 return_vars = list(latitude = 'dat',
-                                    longitude = 'dat',
-                                    time = 'file_date'),
-                 split_multiselected_dims = TRUE,
-                 retrieve = TRUE)
-    
-  } else if (store.freq %in% c("daily_mean", "daily")) {
-    
-    # Get year and month for file_date
-    dates_file <- sapply(dates, format, '%Y%m')
-    dim(dates_file) <- dim(dates)
-    # Set hour to 12:00 to ensure correct date retrieval for daily data
-    lubridate::hour(dates) <- 12 
-    lubridate::minute(dates) <- 00
-    # Restore correct dimensions
-    dim(dates) <- dim(dates_file)
-    
-    obs <- Start(dat = obs.path,
-                 var = variable,
-                 var_dir = var_dir_obs,
-                 var_dir_depends = 'var',
-                 file_date = sort(unique(dates_file)),
-                 time = dates,
-                 time_var = 'time',
-                 time_across = 'file_date',
-                 merge_across_dims = TRUE,
-                 merge_across_dims_narm = TRUE,
-                 latitude = values(list(lats.min, lats.max)),
-                 latitude_reorder = Sort(),
-                 longitude = values(list(lons.min, lons.max)),
-                 longitude_reorder = circularsort,
-                 transform = regrid_params$obs.transform,
-                 transform_params = list(grid = regrid_params$obs.gridtype,
-                                         method = regrid_params$obs.gridmethod),
-                 transform_vars = c('latitude', 'longitude'),
-                 synonims = list(latitude = c('lat','latitude'),
-                                 longitude = c('lon','longitude')),
-                 metadata_dims = 'var',
-                 return_vars = list(latitude = 'dat',
-                                    longitude = 'dat',
-                                    time = 'file_date'),
-                 split_multiselected_dims = TRUE,
-                 retrieve = TRUE)
-  }
-  
-  
-  # Remove var_dir dimension
-  if ("var_dir" %in% names(dim(obs))) {
-    obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected")
-  }
-  # Adds ensemble dim to obs (for consistency with hcst/fcst)
-  default_dims <- c(dat = 1, var = 1, sday = 1, 
-                    sweek = 1, syear = 1, time = 1, 
-                    latitude = 1, longitude = 1, ensemble = 1)
-  default_dims[names(dim(obs))] <- dim(obs)
-  dim(obs) <- default_dims
-  
-  # Convert obs to s2dv_cube
-  obs <- as.s2dv_cube(obs)
-  
-  # Check for consistency between hcst and obs grid
-  if (!(recipe$Analysis$Regrid$type == 'none')) {
-    if (!isTRUE(all.equal(as.vector(hcst$lat), as.vector(obs$lat)))) {
-      lat_error_msg <- paste("Latitude mismatch between hcst and obs.",
-                             "Please check the original grids and the",
-                             "regrid parameters in your recipe.")
-      error(recipe$Run$logger, lat_error_msg)
-      hcst_lat_msg <- paste0("First hcst lat: ", hcst$lat[1],
-                             "; Last hcst lat: ", hcst$lat[length(hcst$lat)])
-      info(recipe$Run$logger, hcst_lat_msg)
-      obs_lat_msg <- paste0("First obs lat: ", obs$lat[1],
-                            "; Last obs lat: ", obs$lat[length(obs$lat)])
-      info(recipe$Run$logger, obs_lat_msg)
-      stop("hcst and obs don't share the same latitudes.")
-    }
-    if (!isTRUE(all.equal(as.vector(hcst$lon), as.vector(obs$lon)))) {
-      lon_error_msg <- paste("Longitude mismatch between hcst and obs.",
-                             "Please check the original grids and the",
-                             "regrid parameters in your recipe.")
-      error(recipe$Run$logger, lon_error_msg)
-      hcst_lon_msg <- paste0("First hcst lon: ", hcst$lon[1],
-                             "; Last hcst lon: ", hcst$lon[length(hcst$lon)])
-      info(recipe$Run$logger, hcst_lon_msg)
-      obs_lon_msg <- paste0("First obs lon: ", obs$lon[1],
-                            "; Last obs lon: ", obs$lon[length(obs$lon)])
-      info(recipe$Run$logger, obs_lon_msg)
-      stop("hcst and obs don't share the same longitudes.")
-      
-    }
-  }
-  
-  # Remove negative values in accumulative variables
-  dictionary <- read_yaml("conf/variable-dictionary.yml")
-  for (var_idx in 1:length(variable)) {
-    var_name <- variable[var_idx]
-    if (dictionary$vars[[var_name]]$accum) {
-      info(recipe$Run$logger,
-           paste0("Accumulated variable ", var_name,
-                  ": setting negative values to zero."))
-      # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data,
-      #                                                along = "var",
-      #                                                indices = var_idx, F), 0)
-      obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0
-      hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0
-      if (!is.null(fcst)) {
-        fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0
-      }
-    }
-    
-    # Convert prlr from m/s to mm/day
-    ## TODO: Make a unit conversion function
-#    if (variable[[var_idx]] == "prlr") {
-#      # Verify that the units are m/s and the same in obs and hcst
-#      if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-#           (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
-#          ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-#           (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
-#        info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
-#        obs$data[, var_idx, , , , , , , ] <-
-#          obs$data[, var_idx, , , , , , , ]*86400*1000
-#        obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-#        hcst$data[, var_idx, , , , , , , ] <-
-#          hcst$data[, var_idx, , , , , , , ]*86400*1000
-#        hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-#        if (!is.null(fcst)) {
-#          fcst$data[, var_idx, , , , , , , ] <- 
-#            fcst$data[, var_idx, , , , , , , ]*86400*1000
-#          fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-#        }
-#      }
-#    }
-  }
-  # 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)
-    data_summary(obs, recipe)
-    if (!is.null(fcst)) {
-      data_summary(fcst, recipe)
-    }
+    stop("Incorrect time horizon.")
   }
-  
-  info(recipe$Run$logger,
-       "##### DATA LOADING COMPLETED SUCCESSFULLY #####")
-  
-  ############################################################################
-  #
-  #     CHECKS ON MISSING FILES
-  #
-  ############################################################################
-  
-  #obs.NA_dates.ind <- Apply(obs, 
-  #                          fun=(function(x){ all(is.na(x))}), 
-  #                          target_dims=c('time', 'latitude', 'longitude'))[[1]]
-  #obs.NA_dates <- dates_file[obs.NA_dates.ind]
-  #obs.NA_dates <- obs.NA_dates[order(obs.NA_dates)]
-  #obs.NA_files <- paste0(obs.dir, store.freq,"/",variable,"_",
-  #                       freq.obs,"obs.grid","/",variable,"_",obs.NA_dates,".nc")
-  #
-  #if (any(is.na(hcst))){
-  #  fatal(recipe$Run$logger,
-  #    paste(" ERROR: MISSING HCST VALUES FOUND DURING LOADING # ",
-  #          " ################################################# ",
-  #          " ######           MISSING FILES               #### ",
-  #          " ################################################# ",
-  #          "hcst files:",
-  #          hcst.NA_files,
-  #          " ################################################# ",
-  #          " ################################################# ",
-  #          sep="\n"))
-  #  quit(status = 1)
-  #}
-  #
-  #if (any(is.na(obs)) && !identical(obs.NA_dates,character(0))){
-  #  fatal(recipe$logger,
-  #    paste(" ERROR: MISSING OBS VALUES FOUND DURING  LOADING # ",
-  #          " ################################################# ",
-  #          " ######           MISSING FILES               #### ",
-  #          " ################################################# ",
-  #          "obs files:",
-  #          obs.NA_files,
-  #          " ################################################# ",
-  #          " ################################################# ",
-  #          sep="\n"))
-  #  quit(status=1)
-  #}
-  #
-  #info(recipe$logger,
-  #     "######### DATA LOADING COMPLETED SUCCESFULLY ##############")
-  
-  ############################################################################
-  ############################################################################
-  
-  return(list(hcst = hcst, fcst = fcst, obs = obs))
-  
+  return(data)
 }
diff --git a/modules/Loading/helper_loading_decadal.R b/modules/Loading/R/helper_loading_decadal.R
similarity index 100%
rename from modules/Loading/helper_loading_decadal.R
rename to modules/Loading/R/helper_loading_decadal.R
diff --git a/modules/Loading/R/load_decadal.R b/modules/Loading/R/load_decadal.R
index 74c97e299f7e7e232fbf2fb3e261bec9673a394e..ba4e1386571c83563bf7453367c8f0b5476bbce8 100644
--- a/modules/Loading/R/load_decadal.R
+++ b/modules/Loading/R/load_decadal.R
@@ -7,19 +7,17 @@
 ## TODO: remove paths to personal scratchs 
 source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R")
 # Load required libraries/funs
-source("modules/Loading/helper_loading_decadal.R")
+source("modules/Loading/R/helper_loading_decadal.R")
 source("modules/Loading/R/dates2load.R")
 source("modules/Loading/R/check_latlon.R")
 source("modules/Loading/R/get_timeidx.R")
-source("tools/libs.R")
 
 #====================================================================
 
 # recipe_file <- "recipes/atomic_recipes/recipe_decadal.yml"
 # recipe_file <- "recipes/atomic_recipes/recipe_decadal_daily.yml"
 
-load_datasets <- function(recipe) {
-
+load_decadal <- function(recipe) {
   ## 
   archive <- read_yaml(paste0("conf/archive_decadal.yml"))$esarchive
 
@@ -34,7 +32,7 @@ load_datasets <- function(recipe) {
   #-------------------------
   exp.name <- recipe$Analysis$Datasets$System$name #'HadGEM3'
   ref.name <- recipe$Analysis$Datasets$Reference$name #'era5'
-  member <- strsplit(recipe$Analysis$Datasets$System$member, ',')[[1]] #c("r1i1p1f2", "r2i1p1f2")
+  member <- strsplit(recipe$Analysis$Datasets$System$member, ', | |,')[[1]] #c("r1i1p1f2", "r2i1p1f2")
 #  variable <- recipe$Analysis$Variables$name #'tas'
   variable <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]]
   store.freq <- recipe$Analysis$Variables$freq #monthly_mean
@@ -491,47 +489,47 @@ load_datasets <- function(recipe) {
 #-------------------------------------------
 # Step 5. Tune data
 #-------------------------------------------
-  # Remove negative values in accumulative variables
-  dictionary <- read_yaml("conf/variable-dictionary.yml")
-  for (var_idx in 1:length(variable)) {
-    var_name <- variable[var_idx]
-    if (dictionary$vars[[var_name]]$accum) {
-      info(recipe$Run$logger,
-           paste0("Accumulated variable ", var_name,
-                  ": setting negative values to zero."))
-      # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data,
-      #                                                along = "var",
-      #                                                indices = var_idx, F), 0)
-      obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0
-      hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0
-      if (!is.null(fcst)) {
-        fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0
-      }
-    }
-    
-    # Convert prlr from m/s to mm/day
-    ## TODO: Make a unit conversion function
-    if (variable[[var_idx]] == "prlr") {
-      # Verify that the units are m/s and the same in obs and hcst
-      if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-           (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
-          ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
-           (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
-        info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
-        obs$data[, var_idx, , , , , , , ] <-
-          obs$data[, var_idx, , , , , , , ]*86400*1000
-        obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        hcst$data[, var_idx, , , , , , , ] <-
-          hcst$data[, var_idx, , , , , , , ]*86400*1000
-        hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        if (!is.null(fcst)) {
-          fcst$data[, var_idx, , , , , , , ] <- 
-            fcst$data[, var_idx, , , , , , , ]*86400*1000
-          fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
-        }
-      }
-    }
-  }
+  # # Remove negative values in accumulative variables
+  # dictionary <- read_yaml("conf/variable-dictionary.yml")
+  # for (var_idx in 1:length(variable)) {
+  #   var_name <- variable[var_idx]
+  #   if (dictionary$vars[[var_name]]$accum) {
+  #     info(recipe$Run$logger,
+  #          paste0("Accumulated variable ", var_name,
+  #                 ": setting negative values to zero."))
+  #     # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data,
+  #     #                                                along = "var",
+  #     #                                                indices = var_idx, F), 0)
+  #     obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0
+  #     hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0
+  #     if (!is.null(fcst)) {
+  #       fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0
+  #     }
+  #   }
+  #   
+  #   # Convert prlr from m/s to mm/day
+  #   ## TODO: Make a unit conversion function
+  #   if (variable[[var_idx]] == "prlr") {
+  #     # Verify that the units are m/s and the same in obs and hcst
+  #     if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+  #          (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
+  #         ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+  #          (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
+  #       info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
+  #       obs$data[, var_idx, , , , , , , ] <-
+  #         obs$data[, var_idx, , , , , , , ]*86400*1000
+  #       obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+  #       hcst$data[, var_idx, , , , , , , ] <-
+  #         hcst$data[, var_idx, , , , , , , ]*86400*1000
+  #       hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+  #       if (!is.null(fcst)) {
+  #         fcst$data[, var_idx, , , , , , , ] <- 
+  #           fcst$data[, var_idx, , , , , , , ]*86400*1000
+  #         fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+  #       }
+  #     }
+  #   }
+  # }
 
 #-------------------------------------------
 # Step 6. Print summary
diff --git a/modules/Loading/R/load_seasonal.R b/modules/Loading/R/load_seasonal.R
index 8cc46dc7f5e3fd25b82d5846eca2254b997962b5..12a19c6ce7a2b5a77506928c2fc39a219a9aa3c1 100644
--- a/modules/Loading/R/load_seasonal.R
+++ b/modules/Loading/R/load_seasonal.R
@@ -4,11 +4,9 @@ source("/esarchive/scratch/vagudets/repos/csoperational/R/get_regrid_params.R")
 source("modules/Loading/R/dates2load.R")
 source("modules/Loading/R/get_timeidx.R")
 source("modules/Loading/R/check_latlon.R")
-## TODO: Move to prepare_outputs.R
-source("tools/libs.R")
 
-load_datasets <- function(recipe) {
-  
+load_seasonal <- function(recipe) {
+
   # -------------------------------------------
   # Set params -----------------------------------------
   
diff --git a/modules/Loading/load_datasets.R b/modules/Loading/load_datasets.R
new file mode 100644
index 0000000000000000000000000000000000000000..86010780c1667b65406f14b57f7ff03ac1821809
--- /dev/null
+++ b/modules/Loading/load_datasets.R
@@ -0,0 +1,7 @@
+load_datasets <- function(recipe) {
+  warning(paste0("The function load_datasets() has been renamed to: ",
+                 "'Loading()'. The name 'load_datasets' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Loading(recipe))
+}
diff --git a/modules/Saving/Saving.R b/modules/Saving/Saving.R
index 2d067ba394fdc4a3bc0ea47f6033fb0326c8c6a2..7d0520c239faaffeddcd4fe479b6c2a3fbe771cd 100644
--- a/modules/Saving/Saving.R
+++ b/modules/Saving/Saving.R
@@ -14,11 +14,12 @@ source("modules/Saving/R/get_time.R")
 source("modules/Saving/R/get_latlon.R")
 source("modules/Saving/R/get_global_attributes.R")
 source("modules/Saving/R/drop_dims.R")
-save_data <- function(recipe, data,
-		      skill_metrics = NULL,
-		      probabilities = NULL,
-                      agg = 'global',
-		      archive = NULL) {
+
+Saving <- function(recipe, data,
+		   skill_metrics = NULL,
+		   probabilities = NULL,
+                   agg = 'global',
+		   archive = NULL) {
   # Wrapper for the saving functions.
   # recipe: The auto-s2s recipe
   # archive: The auto-s2s archive
diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R
index 6c5c10e654f4b9cc0df093a4af3b754a9433ad88..95b6856a5336d367b7c367de36521744a930889f 100644
--- a/modules/Skill/Skill.R
+++ b/modules/Skill/Skill.R
@@ -20,47 +20,11 @@ source("modules/Skill/R/tmp/MSE.R")
 source("modules/Skill/R/tmp/MSSS.R")
 source("modules/Skill/R/RPS_clim.R")
 source("modules/Skill/R/CRPS_clim.R")
+## TODO: Remove in the next release
+source("modules/Skill/compute_skill_metrics.R")
+source("modules/Skill/compute_probabilities.R")
 
-## TODO: Implement this in the future
-## Which parameter are required?
-# if (!("obs" %in% ls()) || is.null(obs)) {
-#   error(logger,
-#         "There is no object 'obs' in the global environment or it is NULL")
-# }
-# if (stream == "fcst" && (!("fcst" %in% ls()) || is.null(fcst))) {
-#   error(logger,
-#         "There is no object 'fcst' in the global environment or it is NULL")
-# }
-# if (!("hcst" %in% ls()) || is.null(hcst)) {
-#   error(logger,
-#         "There is no object 'hcst' in the global environment or it is NULL")
-# }
-# if (!("metric" %in% ls()) || is.null(metric)) {
-#   warn(logger,
-#        "Verification metric not found and it is set as 'EnsCorr'.")
-#   metric <- 'EnsCorr'
-# }
-# if (metric %in% c('FRPSS', 'RPSS')) {
-#   metric_fun <- "veriApply"
-#   metric_method <- "FairRpss"
-# } else if (metric %in% c("FCRPSS", "CRPSS")) {
-#   metric_fun <- "veriApply"
-# } else if (metric %in% c("EnsCorr", "EnsCor")) {
-#   metric_fun <- "veriApply"
-#   metric_method <- "EnsCorr"
-# #...
-# } else {
-#   error(logger, "Unknown verification metric defined in the recipe.")
-#   metric_fun <- 'NotFound'
-# }
-# info(logger, paste("#-------------------------- ", "\n",
-#                    " running Skill module ", "\n",
-#                    " it can call ", metric_fun ))
-
-# compute_skill_metrics <- function(recipe, data$hcst, obs,
-# 				  clim_data$hcst = NULL,
-# 				  clim_obs = NULL) {
-compute_skill_metrics <- function(recipe, data, agg = 'global') {
+Skill <- function(recipe, data, agg = 'global') {
 
   # data$hcst: s2dv_cube containing the hindcast
   # obs: s2dv_cube containing the observations
@@ -372,9 +336,8 @@ compute_skill_metrics <- function(recipe, data, agg = 'global') {
   return(skill_metrics)
 }
 
-compute_probabilities <- function(recipe, data) {
-  ## TODO: Do hcst and fcst at the same time 
-
+Probabilities <- function(recipe, data) {
+  ## TODO: Do hcst and fcst at the same time
   if (is.null(recipe$Analysis$ncores)) {
     ncores <- 1
   } else {
diff --git a/modules/Skill/compute_probabilities.R b/modules/Skill/compute_probabilities.R
new file mode 100644
index 0000000000000000000000000000000000000000..44a91b96e8573d4dac77cb1f02b59fa9a3c3b3cc
--- /dev/null
+++ b/modules/Skill/compute_probabilities.R
@@ -0,0 +1,7 @@
+compute_probabilities <- function(recipe, data) {
+  warning(paste0("The function compute_probabilities() has been renamed to: ",
+                 "'Probabilities()'. The name 'compute_probabilities' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Probabilities(recipe, data))
+}
diff --git a/modules/Skill/compute_skill_metrics.R b/modules/Skill/compute_skill_metrics.R
new file mode 100644
index 0000000000000000000000000000000000000000..98f0f2ebd2f0975fa483bc808e3f1a9df8ecfec4
--- /dev/null
+++ b/modules/Skill/compute_skill_metrics.R
@@ -0,0 +1,7 @@
+compute_skill_metrics <- function(recipe, data) {
+  warning(paste0("The function compute_skill_metrics() has been renamed to: ",
+                 "'Skill()'. The name 'compute_skill_metrics' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Skill(recipe, data))
+}
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 1ddd6ec747cf78092fb5f93db02fa475b6a0f25c..7df72086c6338f30e4e7431dde8c9c70c2c066c4 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -10,13 +10,15 @@ source("modules/Visualization/R/get_proj_code.R")
 source("modules/Visualization/R/tmp/PlotRobinson.R")
 source("modules/Visualization/R/plot_most_likely_terciles_map.R")
 source("modules/Visualization/R/plot_ensemble_mean.R")
+## TODO: Remove in the next release
+source("modules/Visualization/plot_data.R")
 
-plot_data <- function(recipe,
-                      data,
-                      skill_metrics = NULL,
-                      probabilities = NULL,
-                      significance = F,
-                      output_conf = NULL) {
+Visualization <- function(recipe,
+                          data,
+                          skill_metrics = NULL,
+                          probabilities = NULL,
+                          significance = F,
+                          output_conf = NULL) {
   # Try to produce and save several basic plots.
   # recipe: the auto-s2s recipe as read by read_yaml()
   # data: list containing the hcst, obs and (optional) fcst s2dv_cube objects
diff --git a/modules/Visualization/plot_data.R b/modules/Visualization/plot_data.R
new file mode 100644
index 0000000000000000000000000000000000000000..910b40fa16c6c912649101c9358fbe1a1973a932
--- /dev/null
+++ b/modules/Visualization/plot_data.R
@@ -0,0 +1,11 @@
+plot_data <- function(recipe,
+                      data,
+                      skill_metrics = NULL,
+                      probabilities = NULL,
+                      significance = F) {
+  warning(paste0("The function plot_data() has been renamed to: ",
+                 "'Visualization()'. The name 'plot_data()' will be ",
+                 "deprecated in the next release. Please change your scripts ",
+                 "accordingly."))
+  return(Visualization(recipe, data, skill_metrics, probabilities, significance))
+}
diff --git a/recipes/recipe_decadal_split.yml b/recipes/recipe_decadal_split.yml
index a94ffd0307e5fe9db4be48d5a53f04c9d81f61da..70764027b20671618ac8f9df64eb5e90679adeeb 100644
--- a/recipes/recipe_decadal_split.yml
+++ b/recipes/recipe_decadal_split.yml
@@ -8,7 +8,7 @@ Analysis:
     - {name: pr, freq: monthly_mean}
   Datasets:
     System: 
-      name: EC-Earth3-i4
+      name: EC-Earth3-i4,
       member: r1i4p1f1,r2i4p1f1
     Multimodel: no
     Reference:
@@ -41,6 +41,7 @@ Analysis:
       save: 'all'
     Visualization:
       plots: skill_metrics
+      multi_panel: yes
     Indicators:
       index: FALSE
   ncores: 8 # Optional, int: number of cores, defaults to 1
@@ -50,9 +51,9 @@ Run:
   Loglevel: INFO
   Terminal: yes
   filesystem: esarchive
-  output_dir: /esarchive/scratch/cdelgado/auto-s2s_logs/
+  output_dir: /esarchive/scratch/vagudets/auto-s2s_logs/
   code_dir: /esarchive/scratch/cdelgado/gitlab/auto-s2s/
-  autosubmit: yes
+  autosubmit: no
   # fill only if using autosubmit
   auto_conf:
     script: /esarchive/scratch/cdelgado/gitlab/cdelgado_copernicus/ESS_evaluation_tool/main_decadal.R
diff --git a/recipes/tests/recipe_seasonal_example.yml b/recipes/tests/recipe_seasonal_example.yml
index 787484ead1627b230a3ba820efc45835aa2646ed..ccd2c7f29ec443713bf2e6f8c8ca8f2d4b9eac7b 100644
--- a/recipes/tests/recipe_seasonal_example.yml
+++ b/recipes/tests/recipe_seasonal_example.yml
@@ -13,7 +13,7 @@ Description:
 Analysis:
   Horizon: Seasonal
   Variables:
-    - {name: tas, freq: monthly_mean}
+    - {name: tas prlr, freq: monthly_mean}
     - {name: prlr, freq: monthly_mean}
   Datasets:
     System:
diff --git a/tests/recipes/recipe-seasonal_daily_1.yml b/tests/recipes/recipe-seasonal_daily_1.yml
index f70f0c0369d0ff16dc45b00a29a2fbbd63eaaee2..42603c2cdb0bd7d421f9d1c48e6a7070e7ed4bf4 100644
--- a/tests/recipes/recipe-seasonal_daily_1.yml
+++ b/tests/recipes/recipe-seasonal_daily_1.yml
@@ -4,14 +4,13 @@ Description:
 Analysis:
   Horizon: Seasonal
   Variables:
-    name: tas
-    freq: daily_mean
+    - {name: tas, freq: daily_mean}
   Datasets:
     System:
-      name: ECMWF-SEAS5
+      - name: ECMWF-SEAS5
     Multimodel: False
     Reference:
-      name: ERA5
+      - name: ERA5
   Time:
     sdate: '1201'
     fcst_year:
@@ -20,10 +19,7 @@ Analysis:
     ftime_min: 1
     ftime_max: 1
   Region:
-    latmin: 17
-    latmax: 20
-    lonmin: 12
-    lonmax: 15
+    - {latmin: 17, latmax: 20, lonmin: 12, lonmax: 15}
   Regrid:
     method: conservative
     type: to_system
diff --git a/tests/testthat/test-decadal_daily_1.R b/tests/testthat/test-decadal_daily_1.R
index c26b89785a8c31afcac87b98441a3b0d758e9421..c847fd10074ee8c457a32b83a2b96f5d44dc9749 100644
--- a/tests/testthat/test-decadal_daily_1.R
+++ b/tests/testthat/test-decadal_daily_1.R
@@ -2,7 +2,7 @@ context("Decadal daily data - 1")
 
 ###########################################
 
-source("modules/Loading/Loading_decadal.R")
+source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
@@ -13,18 +13,18 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 ## Calibrate datasets
 #suppressWarnings({invisible(capture.output(
-#  calibrated_data <- calibrate_datasets(data, recipe)
+#  calibrated_data <- Calibration(data, recipe)
 #))})
 #
 ## Compute skill metrics
 #suppressWarnings({invisible(capture.output(
-#skill_metrics <- compute_skill_metrics(calibrated_data$hcst, data$obs,
-#                                       recipe, na.rm = T, ncores = 4)
+#skill_metrics <- Skill(calibrated_data$hcst, data$obs,
+#                       recipe, na.rm = T, ncores = 4)
 #))})
 
 #======================================
diff --git a/tests/testthat/test-decadal_monthly_1.R b/tests/testthat/test-decadal_monthly_1.R
index a5c95e0cd5330f50691e01b7e3d997ddc951003c..3346e529ae76f9e58c08371f95cf0ded95fb6533 100644
--- a/tests/testthat/test-decadal_monthly_1.R
+++ b/tests/testthat/test-decadal_monthly_1.R
@@ -2,7 +2,7 @@ context("Decadal monthly data - 1")
 
 ###########################################
 
-source("modules/Loading/Loading_decadal.R")
+source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
@@ -14,27 +14,27 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Calibrate datasets
 suppressWarnings({invisible(capture.output(
-  calibrated_data <- calibrate_datasets(recipe, data)
+  calibrated_data <- Calibration(recipe, data)
 ))})
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
-probs <- compute_probabilities(recipe, calibrated_data)
+probs <- Probabilities(recipe, calibrated_data)
 ))})
 
 # Plotting
 suppressWarnings({invisible(capture.output(
-plot_data(recipe = recipe, data = calibrated_data,
-	  skill_metrics = skill_metrics, probabilities = probs,
-          significance = T)
+Visualization(recipe = recipe, data = calibrated_data,
+	      skill_metrics = skill_metrics, probabilities = probs,
+              significance = T)
 ))})
 
 #======================================
@@ -297,20 +297,20 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data_b <- load_datasets(recipe)
+data_b <- Loading(recipe)
 ))})
 
 # Calibrate datasets
 suppressWarnings({invisible(capture.output(
-  calibrated_data_b <- calibrate_datasets(recipe, data_b)
+  calibrated_data_b <- Calibration(recipe, data_b)
 ))})
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics_b <- compute_skill_metrics(recipe, calibrated_data_b)
+skill_metrics_b <- Skill(recipe, calibrated_data_b)
 ))})
 suppressWarnings({invisible(capture.output(
-probs_b <- compute_probabilities(recipe, calibrated_data_b)
+probs_b <- Probabilities(recipe, calibrated_data_b)
 ))})
 
 
diff --git a/tests/testthat/test-decadal_monthly_2.R b/tests/testthat/test-decadal_monthly_2.R
index 2b6051097abc9807266fe7920211f31e26d57f2a..9adc16b65deb3b84f9a825d93c175504f6304f2b 100644
--- a/tests/testthat/test-decadal_monthly_2.R
+++ b/tests/testthat/test-decadal_monthly_2.R
@@ -2,7 +2,7 @@ context("Decadal monthly data - 2")
 
 ###########################################
 
-source("modules/Loading/Loading_decadal.R")
+source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
@@ -13,27 +13,27 @@ recipe <- prepare_outputs(recipe_file)
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Calibrate datasets
 suppressWarnings({invisible(capture.output(
-  calibrated_data <- calibrate_datasets(recipe, data)
+  calibrated_data <- Calibration(recipe, data)
 ))})
 
 # Compute skill metrics
 suppressMessages({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
-probs <- compute_probabilities(recipe, calibrated_data)
+probs <- Probabilities(recipe, calibrated_data)
 ))})
 
 # Plotting
 suppressWarnings({invisible(capture.output(
-plot_data(recipe = recipe, data = calibrated_data,
-	  skill_metrics = skill_metrics,
-          probabilities = probs, significance = T)
+Visualization(recipe = recipe, data = calibrated_data,
+	      skill_metrics = skill_metrics,
+              probabilities = probs, significance = T)
 ))})
 
 
diff --git a/tests/testthat/test-decadal_monthly_3.R b/tests/testthat/test-decadal_monthly_3.R
index 7232ebfdcc8ac6ce38da9e9dbc8328b79fe38703..d1b42fd5c763b5f35863afca1695aee3fbbd13ee 100644
--- a/tests/testthat/test-decadal_monthly_3.R
+++ b/tests/testthat/test-decadal_monthly_3.R
@@ -2,7 +2,7 @@ context("Decadal monthly data - 3")
 
 ###########################################
 
-source("modules/Loading/Loading_decadal.R")
+source("modules/Loading/Loading.R")
 source("modules/Calibration/Calibration.R")
 source("modules/Skill/Skill.R")
 source("modules/Saving/Saving.R")
@@ -13,20 +13,20 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive_decadal.yml"))$ar
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Calibrate datasets
 suppressWarnings({invisible(capture.output(
-  calibrated_data <- calibrate_datasets(recipe, data)
+  calibrated_data <- Calibration(recipe, data)
 ))})
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 ))})
 suppressWarnings({invisible(capture.output(
-probs <- compute_probabilities(recipe, calibrated_data)
+probs <- Probabilities(recipe, calibrated_data)
 ))})
 
 #======================================
diff --git a/tests/testthat/test-seasonal_NAO.R b/tests/testthat/test-seasonal_NAO.R
index 7a03ecbff94e6915540eb1518b0a6af4868196cf..7a9c9de88fc9ea6d2f2fde5cef617241071f2532 100644
--- a/tests/testthat/test-seasonal_NAO.R
+++ b/tests/testthat/test-seasonal_NAO.R
@@ -12,12 +12,12 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Compute anomalies
 suppressWarnings({invisible(capture.output(
-ano_data <- compute_anomalies(recipe, data)
+ano_data <- Anomalies(recipe, data)
 ))})
 
 # Compute index NAO
@@ -27,7 +27,7 @@ nao_data <- Indices(recipe = recipe, data = ano_data)
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, nao_data, agg = 'region')
+skill_metrics <- Skill(recipe, nao_data, agg = 'region')
 ))})
 
 
diff --git a/tests/testthat/test-seasonal_daily.R b/tests/testthat/test-seasonal_daily.R
index 8a51cc86af6f2c91e92b70da1a2847a069e33af0..6cfa4384ee93c907789046c1de9955bd554e012d 100644
--- a/tests/testthat/test-seasonal_daily.R
+++ b/tests/testthat/test-seasonal_daily.R
@@ -9,17 +9,17 @@ recipe_file <- "tests/recipes/recipe-seasonal_daily_1.yml"
 recipe <- prepare_outputs(recipe_file, disable_checks = F)
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Calibrate data
 suppressWarnings({invisible(capture.output(
-calibrated_data <- calibrate_datasets(recipe, data)
+calibrated_data <- Calibration(recipe, data)
 ))})
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 ))})
 
 test_that("1. Loading", {
diff --git a/tests/testthat/test-seasonal_monthly.R b/tests/testthat/test-seasonal_monthly.R
index 0e311424b5ef4db3d49f80e9a980a273e42fa8e0..d112acf9a9187bb073de5c1c84778f4d5f74ae8a 100644
--- a/tests/testthat/test-seasonal_monthly.R
+++ b/tests/testthat/test-seasonal_monthly.R
@@ -12,34 +12,34 @@ archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive
 
 # Load datasets
 suppressWarnings({invisible(capture.output(
-data <- load_datasets(recipe)
+data <- Loading(recipe)
 ))})
 
 # Calibrate data
 suppressWarnings({invisible(capture.output(
-calibrated_data <- calibrate_datasets(recipe, data)
+calibrated_data <- Calibration(recipe, data)
 ))})
 
 # Compute skill metrics
 suppressWarnings({invisible(capture.output(
-skill_metrics <- compute_skill_metrics(recipe, calibrated_data)
+skill_metrics <- Skill(recipe, calibrated_data)
 ))})
 
 suppressWarnings({invisible(capture.output(
-probs <- compute_probabilities(recipe, calibrated_data)
+probs <- Probabilities(recipe, calibrated_data)
 ))})
 
 # Saving
 suppressWarnings({invisible(capture.output(
-save_data(recipe = recipe, data = calibrated_data,
-	  skill_metrics = skill_metrics, probabilities = probs)
+Saving(recipe = recipe, data = calibrated_data,
+       skill_metrics = skill_metrics, probabilities = probs)
 ))})
 
 # Plotting
 suppressWarnings({invisible(capture.output(
-plot_data(recipe = recipe,  data = calibrated_data,
-	  skill_metrics = skill_metrics, probabilities = probs,
-	  significance = T)
+Visualization(recipe = recipe,  data = calibrated_data,
+	      skill_metrics = skill_metrics, probabilities = probs,
+	      significance = T)
 ))})
 outdir <- get_dir(recipe = recipe, variable = data$hcst$attrs$Variable$varName)
 
diff --git a/tools/divide_recipe.R b/tools/divide_recipe.R
index e9309be89148bc1cc7d178323f8c1ad39b52f73a..f4a56fb1c56ba4305fbd777fd2a9e27f19aa70ff 100644
--- a/tools/divide_recipe.R
+++ b/tools/divide_recipe.R
@@ -35,6 +35,13 @@ divide_recipe <- function(recipe) {
   
   # duplicate recipe by Datasets:
   # check Systems
+
+  if (any(c("name", "member") %in% names(recipe$Analysis$Datasets$System))) {
+    system <- recipe$Analysis$Datasets$System
+    recipe$Analysis$Datasets$System <- NULL
+    recipe$Analysis$Datasets$System[[1]] <- system
+  }
+  
   if (recipe$Analysis$Datasets$Multimodel %in% c(TRUE, 'both')) {
     for (reci in 1:length(all_recipes)) {
       all_recipes[[reci]]$Analysis$Datasets <- 
@@ -43,34 +50,17 @@ divide_recipe <- function(recipe) {
              Reference = NULL)
     }
   } else {
-    if (tolower(recipe$Analysis$Horizon) == 'seasonal') {
-      for (sys in 1:length(recipe$Analysis$Datasets$System)) {
-        for (reci in 1:length(all_recipes)) {
-          all_recipes[[reci]]$Analysis$Datasets <- 
-            list(System = recipe$Analysis$Datasets$System[[sys]],
-                 Multimodel = recipe$Analysis$Datasets$Multimodel,
-                 Reference = NULL)
-        }
-        if (sys == 1) {
-          recipes <- all_recipes
-        } else {
-          recipes <- append(recipes, all_recipes)
-        }
+    for (sys in 1:length(recipe$Analysis$Datasets$System)) {
+      for (reci in 1:length(all_recipes)) {
+        all_recipes[[reci]]$Analysis$Datasets <- 
+          list(System = recipe$Analysis$Datasets$System[[sys]],
+               Multimodel = recipe$Analysis$Datasets$Multimodel,
+               Reference = NULL)
       }
-    } else if (tolower(recipe$Analysis$Horizon) == 'decadal') {
-      for (sys in 1:length(recipe$Analysis$Datasets$System$name)) {
-        for (reci in 1:length(all_recipes)) {
-          all_recipes[[reci]]$Analysis$Datasets <- 
-            list(System = list(name = recipe$Analysis$Datasets$System$name[[sys]],
-                               member = recipe$Analysis$Datasets$System$member[[sys]]),
-                 Multimodel = recipe$Analysis$Datasets$Multimodel,
-                 Reference = NULL)
-        }
-        if (sys == 1) {
-          recipes <- all_recipes
-        } else {
-          recipes <- append(recipes, all_recipes)
-        }
+      if (sys == 1) {
+        recipes <- all_recipes
+      } else {
+        recipes <- append(recipes, all_recipes)
       }
     } # Rest of horizons
     all_recipes <- recipes
@@ -80,7 +70,7 @@ divide_recipe <- function(recipe) {
   for (ref in 1:length(recipe$Analysis$Datasets$Reference)) {
     for (reci in 1:length(all_recipes)) {
       all_recipes[[reci]]$Analysis$Datasets$Reference <- 
-	recipe$Analysis$Datasets$Reference[[ref]]
+	recipe$Analysis$Datasets$Reference[ref]
     }
     if (ref == 1) {
         recipes <- all_recipes
@@ -91,6 +81,12 @@ divide_recipe <- function(recipe) {
   all_recipes <- recipes
   # Duplicate recipe by Region
   recipes <- list()
+  if (any(c("latmin", "latmax", "lonmin", "lonmax") %in% 
+          names(recipe$Analysis$Region))) {
+    region <- recipe$Analysis$Region
+    recipe$Analysis$Region <- NULL
+    recipe$Analysis$Region[[1]] <- region
+  }
   for (reg in 1:length(recipe$Analysis$Region)) {
     for (reci in 1:length(all_recipes)) {
       all_recipes[[reci]]$Analysis$Region <- recipe$Analysis$Region[[reg]]
@@ -101,6 +97,7 @@ divide_recipe <- function(recipe) {
       recipes <- append(recipes, all_recipes)
     }
   }
+
   all_recipes <- recipes
   rm(list = 'recipes')
   # Duplicate recipe by start date
diff --git a/tools/libs.R b/tools/libs.R
index a67f9549fe2236d8c014ebf8f4af0bca111b6500..ec34ad2e91387998375289cbef2ffc91b1873547 100644
--- a/tools/libs.R
+++ b/tools/libs.R
@@ -29,6 +29,7 @@ source("tools/data_summary.R")
 source("tools/read_atomic_recipe.R")
 source("tools/write_autosubmit_conf.R")
 source("tools/get_archive.R")
+source("tools/restructure_recipe.R")
 # source("tools/add_dims.R") # Not sure if necessary yet
 
 # Settings
diff --git a/tools/prepare_outputs.R b/tools/prepare_outputs.R
index f9e71132705a91aca152c0a05390460d63fbfe59..8f9968f45dd7e0d4b14770539f154a552c9571f6 100644
--- a/tools/prepare_outputs.R
+++ b/tools/prepare_outputs.R
@@ -97,5 +97,7 @@ prepare_outputs <- function(recipe_file,
   } else {
     check_recipe(recipe)
   }
+  # Restructure the recipe to make the atomic recipe more readable
+  recipe <- restructure_recipe(recipe)
   return(recipe)
 }
diff --git a/tools/read_atomic_recipe.R b/tools/read_atomic_recipe.R
index de2ad5b554e81dec54146afeefd2724cf8b0c069..fb26cb11a4e71b4b7059a48acccd0ec4c85d235d 100644
--- a/tools/read_atomic_recipe.R
+++ b/tools/read_atomic_recipe.R
@@ -53,5 +53,7 @@ read_atomic_recipe <- function(recipe_file) {
   }
   recipe$Run$logger <- logger
   recipe$Run$logfile <- logfile
+  # Restructure recipe to flatten redundant lists
+  recipe <- restructure_recipe(recipe)
   return(recipe)
 }
diff --git a/tools/restructure_recipe.R b/tools/restructure_recipe.R
new file mode 100644
index 0000000000000000000000000000000000000000..0a9f28f42f1e8b27cb2957242320b920ba925b61
--- /dev/null
+++ b/tools/restructure_recipe.R
@@ -0,0 +1,24 @@
+restructure_recipe <- function(recipe) {
+  # Flattens the structure of the recipe to improve the readability of the code
+  # System
+  if ((length(recipe$Analysis$Datasets$System) == 1) &&
+      (is.list(recipe$Analysis$Datasets$System[[1]]))) {
+    recipe$Analysis$Datasets$System <- recipe$Analysis$Datasets$System[[1]]
+  }
+  # Reference
+  if ((length(recipe$Analysis$Datasets$Reference) == 1) &&
+      (is.list(recipe$Analysis$Datasets$Reference[[1]]))) {
+    recipe$Analysis$Datasets$Reference <- recipe$Analysis$Datasets$Reference[[1]]
+  }
+  # Variable
+  if ((length(recipe$Analysis$Variables) == 1) &&
+      (is.list(recipe$Analysis$Variables[[1]]))) {
+    recipe$Analysis$Variables <- recipe$Analysis$Variables[[1]]
+  }
+  # Region
+  if ((length(recipe$Analysis$Region) == 1) &&
+      (is.list(recipe$Analysis$Region[[1]]))) {
+    recipe$Analysis$Region <- recipe$Analysis$Region[[1]]
+  }
+  return(recipe)
+}