diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..d95b20cb549a163332b75a724f5fe33078f08724
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+out-logs/
diff --git a/OperationalCS.R b/OperationalCS.R
index 03bfe7d1f0bd39258eca5f96d18e0f5d48853657..1d79b4ff48cd51135ff949063e1d0ef352078721 100644
--- a/OperationalCS.R
+++ b/OperationalCS.R
@@ -6,31 +6,24 @@ args = commandArgs(trailingOnly = TRUE)
 # This code should process and check the recipe and build the workflow
 print(args)
 library(yaml)
-conf <- read_yaml(args[1]) # conf <- read_yaml("recipes/conf_user.yml")
-recipe <- read_yaml(args[2]) # recipe <- read_yaml("recipes/seasonal_oper.yml")
+
+# To test:
+# args <- NULL; args[1] <- "recipes/conf_user.yml"; args[2] <- "recipes/seasonal_oper.yml"
+# args <- NULL; args[1] <- "recipes/conf_user.yml"; args[2] <- "recipes/seasonal_complex.yml"
+
+conf <- read_yaml(args[1])
+recipe <- read_yaml(args[2])
 
 # Load required libraries
 source("tools/libs.R")
 
-check_conf(conf, file = args[1]) # check_conf(conf, file = "conf_user.yml")
-verifications <- check_recipe(recipe, file = args[2], conf) 
-# verifications <- check_recipe(recipe, file = "recipe.yml", conf)
-
-for (indep in verifications$independent) {  
-  variable <- indep$name
-  fcst.sdate <- recipe$Analysis$Time$forecast_sdate
-  fcst.freq <- indep$freq
-  hcst.inityear <- recipe$Analysis$Time$hindcast_start
-  hcst.endyear <- recipe$Analysis$Time$hindcast_end
-  fcst.name <- recipe$Analysis$Datasets$System[[1]]$name
-  ref.name <- recipe$Analysis$Datasets$Reference[[1]]$name
-  ltmin <- recipe$Analysis$Time$leadtimemin 
-  ltmax <- recipe$Analysis$Time$leadtimemax
-  lats.min <- recipe$Analysis$Region$Regional[[1]]$latmin
-  lats.max <- recipe$Analysis$Region$Regional[[1]]$latmax
-  lons.min <- recipe$Analysis$Region$Regional[[1]]$lonmin
-  lons.max <- recipe$Analysis$Region$Regional[[1]]$lonmax
-  stream <- "fcst"
-  #source("modules/data_load/load.R")
-}
+# Create output folder and log:
+logger <- prepare_outputs(recipe = recipe, file = args[2], conf = conf) 
+
+# Checks:
+check_conf(conf, file = args[1], logger) 
+verifications <- check_recipe(recipe, file = args[2], conf, logger) 
+
+# Go to verification code: 
+source("modules/verifications.R")
 
diff --git a/conf/archive.yml b/conf/archive.yml
index 0981e0d919ddbe2ffa90695c8d2232a6620872a5..f55a2e90e91e5c5a3b9117c9d2fd08e2d5c58678 100644
--- a/conf/archive.yml
+++ b/conf/archive.yml
@@ -4,10 +4,10 @@ archive:
     esarchive: "/esarchive/oper/S2S-downscale/" 
     system5c3s: 
         src: "exp/ecmwf/system5c3s/"
-        daily_mean: {"tas":"tas_f6h/","rsds":"rsds_s0-24h/",
-                     "prlr":"prlr_s0-24h/","sfcWind":"sfcWind_f6h/"}
-        monthly_mean: {"tas":"tas_f6h/","rsds":"rsds_s0-24h/",
-                     "prlr":"prlr_s0-24h/","sfcWind":"sfcWind_f6h/"}
+        daily_mean: {"tas":"_f6h/","rsds":"_s0-24h/",
+                     "prlr":"_s0-24h/","sfcWind":"_f6h/"}
+        monthly_mean: {"tas":"_f6h/","rsds":"_s0-24h/",
+                     "prlr":"_s0-24h/","sfcWind":"_f6h/"}
         lat_decreasing_sort: TRUE
         lon_circular_sort:
           ini: 0
@@ -15,17 +15,26 @@ archive:
         nmember:
           fcst: 51
           hcst: 25
+        regrid: "/esarchive/exp/ecmwf/system5c3s/monthly_mean/tas_f6h/tas_20180501.nc"
 
     era5: 
         src: "recon/ecmwf/era5/"
-        daily_mean: {"tas":"tas_f1h-r1440x721cds/","rsds":"rsds_f1h-r1440x721cds/",
-                     "prlr":"prlr_f1h-r1440x721cds/","sfcWind":"sfcWind_f1h-r1440x721cds/"}
+        daily_mean: {"tas":"_f1h-r1440x721cds/","rsds":"_f1h-r1440x721cds/",
+                     "prlr":"_f1h-r1440x721cds/","sfcWind":"_f1h-r1440x721cds/"}
+        monthly_mean: {"tas":"_f1h-r1440x721cds/"}
+        lat_decreasing_sort: TRUE
+        lon_circular_sort:
+          ini: 0
+          end: 361
+        regrid: "/esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h-r1440x721cds/tas_201805.nc"
     era5land: 
         src: "recon/ecmwf/era5land/"
-        daily_mean: {"tas":"tas_f1h/","rsds":"rsds_f1h/",
-                     "prlr":"prlr_f1h/","sfcWind":"sfcWind_f1h/"}
+        daily_mean: {"tas":"_f1h/","rsds":"_f1h/",
+                     "prlr":"_f1h/","sfcWind":"_f1h/"}
+        regrid: "/esarchive/recon/ecmwf/era5land/daily_mean/tas_f1h/tas_201805.nc"
     uerra: 
         src: "recon/ecmwf/uerra_mescan/"
-        daily_mean: {"tas":"tas_f6h/"}
+        daily_mean: {"tas":"_f6h/"}
+        grid: "/esarchive/recon/ecmwf/uerra_mescan/daily_mean/tas_f6h/tas_201805.nc"
 
 
diff --git a/conf/indicators_table.yml b/conf/indicators_table.yml
new file mode 100644
index 0000000000000000000000000000000000000000..8e8ce0fdbc391ca49d164783bd1ef8ce6fbd0e75
--- /dev/null
+++ b/conf/indicators_table.yml
@@ -0,0 +1,19 @@
+
+gdd:
+  longname: 'Growing Degree Days'
+  ECVs: tas
+  freq: daily_mean
+  fun: AccumulationExceedingThreshold
+
+gst:
+  longname: 'Growing Season Temperature'
+  ECVs: tas
+  freq: daily_mean
+  fun: PeriodMean
+
+spr32:
+  longname: 'Spring Heat Stress Days - 32°C'
+  ECVs: tasmax
+  freq: daily_mean
+  fun: TotalTimeExceedingThreshold
+
diff --git a/main.R b/main.R
deleted file mode 100644
index 648b943a528d23d5243c80c2aced4efbefb67ed7..0000000000000000000000000000000000000000
--- a/main.R
+++ /dev/null
@@ -1,40 +0,0 @@
-#!/usr/bin/env Rscript
-args = commandArgs(trailingOnly = TRUE)
-
-# Author: N.Pérez-Zanón 26th July 2021
-# execution: Rscript OperationalCS.R conf_user.yml recipe.yml
-# This code should process and check the recipe and build the workflow
-print(args)
-library(yaml)
-conf <- read_yaml(args[1]) # conf <- read_yaml("recipes/conf_user.yml")
-recipe <- read_yaml(args[2]) # recipe <- read_yaml("recipes/seasonal_oper.yml")
-
-# Load required libraries
-source("tools/libs.R")
-
-#check_conf(conf, file = args[1]) 
-#verifications <- check_recipe(recipe, file = args[2], conf) 
-# check_conf(conf, file = "conf_user.yml")
-# verifications <- check_recipe(recipe, file = "recipe.yml", conf)
-
-recipe <- recipe$Analysis
-
-variable <- recipe$Variables$ECVs[[1]]$name
-fcst.sdate <- recipe$Time$sdates
-fcst.freq <- recipe$Variables$ECVs[[1]]$freq
-hcst.inityear <- recipe$Time$hcst_start
-hcst.endyear <- recipe$Time$hcst_end
-
-fcst.name <- recipe$Datasets$System[[1]]$name
-ref.name <- recipe$Datasets$Reference[[1]]$name
-ltmin <- recipe$Time$leadtimemin 
-ltmax <- recipe$Time$leadtimemax
-lats.min <- recipe$Region$Regional[[1]]$latmin
-lats.max <- recipe$Region$Regional[[1]]$latmax
-lons.min <- recipe$Region$Regional[[1]]$lonmin
-lons.max <- recipe$Region$Regional[[1]]$lonmax
-stream <- "fcst"
-
-source("modules/data_load/load.R")
-
-
diff --git a/modules/data_load/load.R b/modules/data_load/load.R
deleted file mode 100644
index fac981c67f4052d1754612da56725645a32512c6..0000000000000000000000000000000000000000
--- a/modules/data_load/load.R
+++ /dev/null
@@ -1,210 +0,0 @@
-  
-  
-fcst.month <- substr(fcst.sdate,5,6) 
-fcst.year <- substr(fcst.sdate,1,4) 
-
-file_dates.fcst <- paste(fcst.year, fcst.month, sep = "")
-file_dates <- paste(strtoi(hcst.inityear):strtoi(hcst.endyear),
-                       fcst.month, sep = "")
-                       
-add_dims <- function(data){
-  default_dims <- c(sweek = 1, sday = 1, syear = length(data))
-  default_dims[names(dim(data))] <- dim(data)
-  dim(data) <- default_dims
-  return(data)
-}
-
-file_dates <- add_dims(file_dates)
-file_dates.fcst <- add_dims(file_dates.fcst)
-
-table <- read_yaml(paste0(conf$code_dir, "conf/archive.yml"))
-dataset_descrip <- table$archive[[fcst.name]]
-fcst.nmember <- dataset_descrip$nmember$fcst
-hcst.nmember <- dataset_descrip$nmember$hcst
-fcst.dir <- dataset_descrip$src 
-hcst.dir <- dataset_descrip$src
-
-store_freq <- dataset_descrip[[fcst.freq]][[variable]]
-freq.hcst <- unlist(dataset_descrip$store_freq[store_freq])
-
-refere_descrip <- table[which(names(table) == ref.name)][[1]]
-freq.obs <- unlist(refere_descrip$store_freq[store_freq])
-obs.dir <- refere_descrip$path
-
-if ("accum" %in% names(refere_descrip)) {
-  accum <- unlist(refere_descrip$accum[store_freq][[1]])
-} else {
-  accum <- FALSE
-}
-# -----------
-  obs.path <- paste0(
-              obs.dir,fcst.freq,"/$var$_",
-              freq.obs,"/$var$_$file_date$.nc")
-
-  hcst.path <- paste0(
-               hcst.dir,fcst.freq,"/", "/$var$_",
-               freq.hcst,"/$var$_$file_date$01",
-               ".nc")
-
-  fcst.path <- paste0(
-               fcst.dir,fcst.freq,"/", "/$var$_",
-               freq.hcst,"/$var$_$file_date$01",
-               ".nc")
-  #-------------------------------------------------------------------
-  # Timeseries load 
-  #-------------------------------------------------------------------
-  
-  if (tolower(stream) == "fcst"){
-    
-    fcst <- Start(dat = fcst.path,
-                  var = variable,
-                  file_date = file_dates.fcst,
-                  time = indices(ltmin:ltmax),
-                  latitude = values(list(lats.min, lats.max)),
-                  latitude_reorder = Sort(decreasing = TRUE),
-                  longitude = values(list(lons.min, lons.max)),
-                  longitude_reorder = CircularSort(0, 361),
-                  synonims = list(latitude=c('lat','latitude'),
-                                  longitude=c('lon','longitude'),
-                                  member=c('ensemble')), 
-                  member = indices(1:fcst.nmember),
-                  return_vars = list(latitude = 'dat',
-                                     longitude = 'dat',
-                                     time = 'file_date'),
-                  split_multiselected_dims = TRUE,
-                  retrieve = T)
-  }
-
-  hcst<-Start(dat = hcst.path,
-              var = variable,
-              file_date = file_dates,
-              time = indices(ltmin:ltmax),
-              latitude = values(list(lats.min, lats.max)),
-              latitude_reorder = Sort(decreasing = TRUE),
-              longitude = values(list(lons.min, lons.max)),
-              longitude_reorder = CircularSort(0, 361),
-              synonims = list(latitude=c('lat','latitude'),
-                              longitude=c('lon','longitude'),
-                              member=c('ensemble')), 
-              member = indices(1:hcst.nmember),
-              return_vars = list(latitude = 'dat',
-                                 longitude = 'dat',
-                                 time = 'file_date'),
-              split_multiselected_dims = TRUE,
-              retrieve = T)
-
-  hcst.NA_files <- c(attributes(hcst)$NotFoundFiles)
-  hcst.NA_files <- hcst.NA_files[!is.na(hcst.NA_files)]
-  try(hcst.NA_files <- hcst.NA_files[order(hcst.NA_files)])
-
-  dates <- attr(hcst, 'Variables')$common$time
-  dates_file <- sapply(dates, format, '%Y%m%d')
-  dates_file <- format(as.Date(dates_file, '%Y%m%d'), "%Y%m")
-  dim(dates_file) <- dim(Subset(hcst, 
-                                along=c('dat','var',
-                                'latitude', 'longitude', 'member'), 
-                                list(1,1,1,1,1), drop="selected")) 
-
-  obs <- Start(dat = obs.path,        
-               var = variable,
-               file_date = dates_file,
-               latitude = values(list(lats.min, lats.max)),
-               latitude_reorder = Sort(decreasing = TRUE),
-               longitude = values(list(lons.min, lons.max)),
-               longitude_reorder = CircularSort(0, 361),
-               synonims = list(latitude=c('lat','latitude'),
-                               longitude=c('lon','longitude')), 
-               transform = CDORemapper,
-               transform_params = list(
-               grid = paste0(
-                  "/esarchive/exp/ecmwf/system5c3s/",fcst.freq,
-                  "/",variable,"_",freq.hcst,
-                  "/",variable,"_","201805","01.nc"
-               ),
-               method = 'con',
-               crop = c(lons.min,lons.max,lats.min,lats.max)), 
-               transform_vars = c('latitude', 'longitude'),
-               return_vars = list(latitude = 'dat',
-                                  longitude = 'dat',
-                                  time = 'file_date'),
-               split_multiselected_dims = TRUE,
-               retrieve = TRUE)
-
-  dates_file <- paste0(dates_file, '01')
-  dim(dates_file) <- dim(Subset(hcst,
-                                along=c('dat','var',
-                                'latitude', 'longitude', 'member'),
-                                list(1,1,1,1,1), drop="selected"))
-
-  file_dates <- paste0(file_dates, '01')
-  dim(file_dates) <- dim(Subset(hcst,
-                                along=c('dat','var',
-                                'latitude', 'longitude', 'member', 'time'),
-                                list(1,1,1,1,1,1), drop="selected"))
-
-
-
-  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,fcst.freq,"/",variable,"_",
-                         freq.obs,"obs.grid","/",variable,"_",obs.NA_dates,".nc")
-
-  if (any(is.na(hcst))){
-    print(" ERROR: MISSING HCST VALUES FOUND DURING LOADING # " )
-    print(" ################################################# " )
-    print(" ######           MISSING FILES               #### " )
-    print(" ################################################# " )
-    print("hcst files:")
-    print(hcst.NA_files)
-    print(" ################################################# " )
-    print(" ################################################# " )
-    quit(status=1)
-  }
-
-  if (any(is.na(obs)) && !identical(obs.NA_dates,character(0))){
-    print(" ERROR: MISSING OBS VALUES FOUND DURING  LOADING # " )
-    print(" ################################################# " )
-    print(" ######           MISSING FILES               #### " )
-    print(" ################################################# " )
-    print("obs files:")
-    print(obs.NA_files)
-    print(" ################################################# " )
-    print(" ################################################# " )
-    quit(status=1)
-  }
-
-  print("######### DATA LOADING COMPLETED SUCCESFULLY ##############")
-
-  default_dims <- c(dat=1, var=1, sweek = 1, 
-                    sday = 1, syear = 1, time=1, 
-                    latitude=1, longitude = 1, member = 1)
-
-  default_dims[names(dim(obs))] <- dim(obs)
-  dim(obs) <- default_dims
-
-  lon <- attr(obs, 'Variables')$dat1$longitude
-  lat <- attr(obs, 'Variables')$dat1$latitude
-  hcst.times <- attr(hcst, 'Variables')$common$time
-  hcst.times <- sort(unique(sapply(as.character(hcst.times),
-                     substr, 1, 10)))
-
-  #obs<-Subset(obs,c('dat','var'),list(1,1),drop='selected')
-  #hcst<-Subset(hcst,c('dat','var'),list(1,1),drop='selected')
-  #if (stream == "fcst"){ 
-  #  fcst<-Subset(fcst,c('dat','var'),list(1,1),drop='selected')
-  #}
-
-  #filters negative values in accum vars
-  if (accum){
-     obs[obs < 0 ] <- 0
-     hcst[hcst < 0 ] <- 0
-     if (stream == "fcst"){ fcst[fcst < 0 ] <- 0 }
-  }
-
-  sdates.hcst <- file_dates
-  sdates.fcst <- file_dates.fcst
-  leadtimes.hcst <- dates_file
-
diff --git a/modules/data_load/seas5.load.R b/modules/data_load/seas5.load.R
new file mode 100644
index 0000000000000000000000000000000000000000..4ebd922bb15001b3546d56ef3aa20d50ac7f18fc
--- /dev/null
+++ b/modules/data_load/seas5.load.R
@@ -0,0 +1,249 @@
+  
+  
+  fcst.month <- substr(fcst.sdate,5,6) 
+  fcst.year <- substr(fcst.sdate,1,4) 
+
+  file_dates.fcst <- paste(fcst.year, fcst.month, sep = "")
+  file_dates <- paste(strtoi(hcst.inityear):strtoi(hcst.endyear),
+                         fcst.month, sep = "")
+                         
+
+  file_dates <- add_dims(file_dates)
+  file_dates.fcst <- add_dims(file_dates.fcst)
+# Take parameters from conf/archive for datasets:
+table <- read_yaml(paste0(conf$code_dir, "conf/archive.yml"))
+dataset_descrip <- table$archive[which(names(table$archive) == fcst.name)][[1]]
+freq.hcst <- unlist(dataset_descrip[store.freq][[1]][variable])
+reference_descrip <- table$archive[which(names(table$archive) == 
+                                tolower(ref.name))][[1]]
+freq.obs <- unlist(reference_descrip[store.freq][[1]][variable])
+obs.dir <- reference_descrip$src
+fcst.dir <- dataset_descrip$src 
+hcst.dir <- dataset_descrip$src
+fcst.nmember <- dataset_descrip$nmember$fcst
+hcst.nmember <- dataset_descrip$nmember$hcst
+
+if ("accum" %in% names(reference_descrip)) {
+  accum <- unlist(reference_descrip$accum[store.freq][[1]])
+} else {
+  accum <- FALSE
+}
+# -----------
+  obs.path <- paste0("/esarchive/",
+              obs.dir, store.freq, "/$var$",
+              freq.obs,"$var$_$file_date$.nc")
+
+  hcst.path <- paste0("/esarchive/",
+               hcst.dir, store.freq, "/$var$",
+               freq.hcst,"$var$_$file_date$01",
+               ".nc")
+
+  fcst.path <- paste0("/esarchive/",
+               fcst.dir, store.freq, "/$var$",
+               freq.hcst,"$var$_$file_date$01",
+               ".nc")
+#-------------------------------------------------------------------
+# Regrid:
+if (tolower(recipe$Analysis$Regrid$type) == 'reference') {
+  fcst.gridtype <- reference_descrip$regrid 
+  fcst.gridmethod <- recipe$Analysis$Regrid$method
+  fcst.tranform <- CDORemapper
+  obs.gridtype <- NULL
+  obs.gridmethod <- NULL
+  obs.tranform <- NULL
+} else if (tolower(recipe$Analysis$Regrid$type) == 'system') {
+  fcst.gridtype <- NULL
+  fcst.gridmethod <- NULL
+  fcst.transform <- NULL
+  obs.gridtype <- dataset_descrip$regrid
+  obs.gridmethod <- recipe$Analysis$Regrid$method
+  obs.transform <- CDORemapper 
+} else {
+  fcst.gridtype <- recipe$Analysis$Regrid$type
+  fcst.gridmethod <- recipe$Analysis$Regrid$method
+  fcst.transform <- CDORemapper
+  obs.gridtype <- recipe$Analysis$Regrid$type
+  obs.gridmethod <- recipe$Analysis$Regrid$method
+  obs.transform <- CDORemapper
+}
+
+
+  # Timeseries load 
+  #-------------------------------------------------------------------
+  
+  if (tolower(stream) == "fcst"){
+    
+    fcst <- Start(dat = fcst.path,
+                  var = variable,
+                  file_date = file_dates.fcst,
+                  time = indices(ltmin:ltmax),
+  	              latitude = values(list(lats.min, lats.max)),
+                  latitude_reorder = Sort(decreasing = 
+                                          dataset_descrip$lat_decreasing_sort),
+                  longitude = values(list(lons.min, lons.max)),
+                  longitude_reorder = CircularSort(
+                            dataset_descrip$lon_circular_sort$ini,
+                            dataset_descrip$lon_circular_sort$end),
+                  transform = fcst.transform,
+                  transform_params = list(grid = fcst.gridtype,
+                                          method = fcst.gridmethod,
+                       		          crop = c(lons.min, lons.max,
+                                                   lats.min, lats.max)),
+                  transform_vars = c('latitude', 'longitude'),
+                  synonims = list(latitude = c('lat','latitude'),
+                                  longitude = c('lon','longitude'),
+                                  member = c('ensemble')), 
+                  member = indices(1:fcst.nmember),
+                  return_vars = list(latitude = 'dat',
+                                     longitude = 'dat',
+                                     time = 'file_date'),
+                  split_multiselected_dims = TRUE,
+                  retrieve = TRUE)
+  }
+
+  hcst <- Start(dat = hcst.path,
+                var = variable,
+                file_date = file_dates,
+                time = indices(ltmin:ltmax),
+                latitude = values(list(lats.min, lats.max)),
+                latitude_reorder = Sort(decreasing =  
+                                        dataset_descrip$lat_decreasing_sort),
+                longitude = values(list(lons.min, lons.max)),
+                longitude_reorder = CircularSort(
+                            dataset_descrip$lon_circular_sort$ini,
+                            dataset_descrip$lon_circular_sort$end),
+                transform = fcst.transform,
+                transform_params = list(grid = fcst.gridtype,
+                                        method = fcst.gridmethod,
+                                        crop = c(lons.min, lons.max,
+                                                 lats.min, lats.max)),
+                transform_vars = c('latitude', 'longitude'),
+                synonims = list(latitude = c('lat','latitude'),
+                                longitude = c('lon','longitude'),
+                                member = c('ensemble')), 
+                member = indices(1:hcst.nmember),
+                return_vars = list(latitude = 'dat',
+                                   longitude = 'dat',
+                                   time = 'file_date'),
+                split_multiselected_dims = TRUE,
+                retrieve = TRUE)
+
+  hcst.NA_files <- c(attributes(hcst)$NotFoundFiles)
+  hcst.NA_files <- hcst.NA_files[!is.na(hcst.NA_files)]
+  try(hcst.NA_files <- hcst.NA_files[order(hcst.NA_files)])
+
+  dates <- attr(hcst, 'Variables')$common$time
+  dates_file <- sapply(dates, format, '%Y%m%d')
+  dates_file <- format(as.Date(dates_file, '%Y%m%d'), "%Y%m")
+  dim(dates_file) <- dim(Subset(hcst, 
+                                along=c('dat','var',
+                                'latitude', 'longitude', 'member'), 
+                                list(1,1,1,1,1), drop="selected")) 
+
+  obs <- Start(dat = obs.path,		  
+               var = variable,
+  	       file_date = dates_file,
+               latitude = values(list(lats.min, lats.max)),
+               latitude_reorder = Sort(decreasing =  
+                                       reference_descrip$lat_decreasing_sort),
+               longitude = values(list(lons.min, lons.max)),
+               longitude_reorder = CircularSort(
+                            reference_descrip$lon_circular_sort$ini,
+                            reference_descrip$lon_circular_sort$end),
+               synonims = list(latitude = c('lat','latitude'),
+                               longitude = c('lon','longitude')), 
+               transform = obs.transform,
+	       transform_params = list(grid = obs.gridtype,
+		                       method = obs.gridmethod,
+		                       crop = c(lons.min, lons.max,
+                                                lats.min, lats.max)), 
+	       transform_vars = c('latitude', 'longitude'),
+  	       return_vars = list(latitude = 'dat',
+                                  longitude = 'dat',
+                                  time = 'file_date'),
+  	       split_multiselected_dims = TRUE,
+  	       retrieve = TRUE)
+
+  dates_file <- paste0(dates_file, '01')
+  dim(dates_file) <- dim(Subset(hcst,
+                                along=c('dat','var',
+                                'latitude', 'longitude', 'member'),
+                                list(1,1,1,1,1), drop="selected"))
+
+  file_dates <- paste0(file_dates, '01')
+  dim(file_dates) <- dim(Subset(hcst,
+                                along=c('dat','var',
+                                'latitude', 'longitude', 'member', 'time'),
+                                list(1,1,1,1,1,1), drop="selected"))
+
+
+
+  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(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(logger,
+      paste(" ERROR: MISSING OBS VALUES FOUND DURING  LOADING # ",
+            " ################################################# ",
+            " ######           MISSING FILES               #### ",
+            " ################################################# ",
+            "obs files:",
+            obs.NA_files,
+            " ################################################# ",
+            " ################################################# ",
+            sep="\n"))
+    quit(status=1)
+  }
+
+  info(logger,
+       "######### DATA LOADING COMPLETED SUCCESFULLY ##############")
+
+  default_dims <- c(dat = 1, var = 1, sweek = 1, 
+                    sday = 1, syear = 1, time = 1, 
+                    latitude = 1, longitude = 1, member = 1)
+
+  default_dims[names(dim(obs))] <- dim(obs)
+  dim(obs) <- default_dims
+
+  lon <- attr(obs, 'Variables')$dat1$longitude
+  lat <- attr(obs, 'Variables')$dat1$latitude
+  hcst.times <- attr(hcst, 'Variables')$common$time
+  hcst.times <- sort(unique(sapply(as.character(hcst.times),
+                     substr, 1, 10)))
+
+  #obs<-Subset(obs,c('dat','var'),list(1,1),drop='selected')
+  #hcst<-Subset(hcst,c('dat','var'),list(1,1),drop='selected')
+  #if (stream == "fcst"){ 
+  #  fcst<-Subset(fcst,c('dat','var'),list(1,1),drop='selected')
+  #}
+
+  #filters negative values in accum vars
+  if (accum){
+     obs[obs < 0 ] <- 0
+     hcst[hcst < 0 ] <- 0
+     if (stream == "fcst"){ fcst[fcst < 0 ] <- 0 }
+  }
+
+  sdates.hcst <- file_dates
+  sdates.fcst <- file_dates.fcst
+  leadtimes.hcst <- dates_file
+
diff --git a/modules/verifications.R b/modules/verifications.R
new file mode 100644
index 0000000000000000000000000000000000000000..2749201a888c6b0d9d12a9432d6f685cf1e3217e
--- /dev/null
+++ b/modules/verifications.R
@@ -0,0 +1,49 @@
+for (indep in verifications$independent) {
+print(indep)
+  if (length(indep) == 1) {
+    info(logger,
+         paste(" #*****************************************#",
+               " # Starting Independent verification of Indicator ",
+               " indep[[1]]$name",
+               sep = "\n"))
+    ind_table <- read_yaml(paste0(conf$code_dir,
+                           "conf/indicators_table.yml"))
+    variable <- ind_table[indep[[1]]][[1]]$ECVs
+    store.freq <- ind_table[indep[[1]]][[1]]$freq
+    ind.fun <- ind_table[indep[[1]]][[1]]$fun
+  } else {
+    info(logger,
+         paste(" #*****************************************#",
+               " # Starting Independent verification of ECV ",
+               " indep[[1]]$name, indep[[1]]$freq",
+               sep = "\n"))
+    variable <- indep$name
+    store.freq <- indep$freq
+    ind.fun <- NULL
+  }
+  hcst.inityear <- recipe$Analysis$Time$hcst_start
+  hcst.endyear <- recipe$Analysis$Time$hcst_end
+  ltmin <- recipe$Analysis$Time$leadtimemin
+  ltmax <- recipe$Analysis$Time$leadtimemax
+  lats.min <- recipe$Analysis$Region$Regional[[1]]$latmin
+  lats.max <- recipe$Analysis$Region$Regional[[1]]$latmax
+  lons.min <- recipe$Analysis$Region$Regional[[1]]$lonmin
+  lons.max <- recipe$Analysis$Region$Regional[[1]]$lonmax
+  stream <- verifications$stream
+  # Is multimodel FALSE
+  if (recipe$Analysis$Datasets$Multimodel) {
+    # is there a multimodel load step in S2S4E backend?
+  } else {
+    for (sys in 1:length(recipe$Analysis$Datasets$System)) {
+      fcst.name <- recipe$Analysis$Datasets$System[[sys]]$name
+      for (ref in 1:length(recipe$Analysis$Datasets$Reference)) {
+        ref.name <- recipe$Analysis$Datasets$Reference[[ref]]$name
+        for (sdate in verifications$fcst.sdate) {
+          fcst.sdate <- sdate
+          source("modules/data_load/seas5.load.R")
+        }
+      }
+    }
+  }
+}
+
diff --git a/recipes/conf_user.yml b/recipes/conf_user.yml
index 2860105dd973aa6c8d9fc623175812d1e84943e7..490da5e25486b063e4aad9107caa6f4d1da8c654 100644
--- a/recipes/conf_user.yml
+++ b/recipes/conf_user.yml
@@ -1,4 +1,4 @@
 # Users configuration to point correct directories
 
-output_dir: /esarchive/scratch/lpalma/tmp/
-code_dir: /esarchive/scratch/lpalma/git/auto-s2s/
+output_dir: /esarchive/scratch/nperez/git/auto-s2s/out-logs/
+code_dir: /esarchive/scratch/nperez/git/auto-s2s/
diff --git a/recipes/seasonal_complex.yml b/recipes/seasonal_complex.yml
new file mode 100644
index 0000000000000000000000000000000000000000..dd0137cd86d324373d1a8ce4682fafde702ebd75
--- /dev/null
+++ b/recipes/seasonal_complex.yml
@@ -0,0 +1,48 @@
+Description:
+  Author: N.Pérez-Zanón
+  Info: This is a test to transform s2s4e data-analysis for SEAS5
+
+Analysis: 
+  Horizon: Seasonal
+  Variables:
+    ECVs: 
+      - {name: tas, freq: monthly_mean}
+      - {name: tas, freq: daily_mean}
+    Indicators:
+      - {name: gdd}
+  Datasets:
+    System: 
+      - name: system5c3s
+    Multimodel: False
+    Reference: 
+      - {name: ERA5}
+  Time:
+    sdate: 
+      fcst_year: 2021
+      fcst_month: [07, 08]
+      fcst_day: 01
+    hcst_start: 2000
+    hcst_end: 2002
+    leadtimemin: 2
+    leadtimemax: 4
+  Region:
+    Global: TRUE
+    Aggregation: False
+    Regional: 
+      - {latmin: -90, latmax: 90, lonmin: 0, lonmax: 360}
+      - {latmin: -10, latmax: 10, lonmin: 0, lonmax: 20} 
+  Regrid:
+    method: bicubic
+    type: system 
+  Workflow:
+    Calibration:
+      method: SBC
+    Skill:
+      metric: RPSS
+    Indicators:
+      index: False 
+  Output_format: S2S4E
+
+Run:
+  Loglevel: INFO
+  Terminal: True
diff --git a/recipes/seasonal_oper.yml b/recipes/seasonal_oper.yml
index b9115974bf75ddc2608d914c4f1c5cf830774a44..9598f4495ff4cc95219895bea6e315cc55df984c 100644
--- a/recipes/seasonal_oper.yml
+++ b/recipes/seasonal_oper.yml
@@ -14,9 +14,12 @@ Analysis:
       - name: system5c3s
     Multimodel: False
     Reference: 
-      - {name: ERA5}
+      - {name: era5}
   Time:
-    sdates: 20210701
+    sdate: 
+      fcst_year: 2021
+      fcst_month: 07
+      fcst_day: 01
     hcst_start: 2000
     hcst_end: 2002
     leadtimemin: 2
@@ -28,13 +31,17 @@ Analysis:
       - {latmin: -90, latmax: 90, lonmin: 0, lonmax: 360}
       - {latmin: -10, latmax: 10, lonmin: 0, lonmax: 20} 
   Regrid:
-    method: ycon
-    type: 180x90
+    method: bilinear
+    type: system
   Workflow:
     Calibration:
       method: SBC
     Skill:
       metric: RPSS
     Indicators:
-      index: False 
+      index: FALSE 
   Output_format: S2S4E
+
+Run:
+  Loglevel: INFO
+  Terminal: TRUE
diff --git a/tools/add_dims.R b/tools/add_dims.R
new file mode 100644
index 0000000000000000000000000000000000000000..aa770eb81825b94e57c3fe04c53c6443ef8fa6bd
--- /dev/null
+++ b/tools/add_dims.R
@@ -0,0 +1,9 @@
+# Function created in S2S$E-backend. Copied from:
+# https://earth.bsc.es/gitlab/es/S2S4E-backend-BSC/-/blob/master/data-analysis/SEAS5/seas5.load.R
+add_dims <- function(data){
+  default_dims <- c(sweek = 1, sday = 1, syear = length(data))
+  default_dims[names(dim(data))] <- dim(data)
+  dim(data) <- default_dims
+  return(data)
+ }
+
diff --git a/tools/check_conf.R b/tools/check_conf.R
index 38ea52366b59806a437b060b2bdbc2251dd9b320..70ffa9f0dcf5f06355e3b9d1030dde97273a7bd4 100644
--- a/tools/check_conf.R
+++ b/tools/check_conf.R
@@ -1,19 +1,21 @@
 
 
 
-check_conf <- function(conf, file) {
+check_conf <- function(conf, file, logger) {
   # conf: yaml user config already read it
   # file: name of the file to clarify messages
-  message(paste("Checking user configuration", file))
+  info(logger, paste("Checking user configuration", file))
   if (!('output_dir' %in% names(conf))) {
-    stop(paste("The user configuration file should have element 'output_dir'."))
+    error(logger,
+         paste("The user configuration file should have element 'output_dir'."))
   }
   if (!is.character(conf$output_dir)) {
-    stop (paste("The element 'output_dir' in", file,"file ",
-          "should be a character string indicating the path ",
-          "where to save the outputs."))
+    error(logger,
+          paste("The element 'output_dir' in", file,"file ",
+                "should be a character string indicating the path ",
+                "where to save the outputs."))
   }
-  message("User configuration checked succsessfully.")
+  info(logger, "User configuration checked succsessfully.")
 }
 
 
diff --git a/tools/check_recipe.R b/tools/check_recipe.R
index 767592335bbfe1313eba1d7b9a5094b00da4b7af..3fea5594c65f3d1cb691c1cfc45bdf745d00319b 100644
--- a/tools/check_recipe.R
+++ b/tools/check_recipe.R
@@ -1,35 +1,95 @@
-check_recipe <- function(recipe, file, conf) {
+check_recipe <- function(recipe, file, conf, logger) {
   # recipe: yaml recipe already read it
   # file: name of the file to clarify messages
   # conf: configuration file read to load the tables of definitions
   # output: errors or the total number of workflow (vars x regions) to compute
-  message(paste("Checking recipe", file))
+  info(logger, paste("Checking recipe", file))
+  # create output dirs:
   if (!any(names(recipe) %in% "Analysis")) {
-    stop("The recipe should has a element 'Analysis'.")
+    error(logger, "The recipe should has a element 'Analysis'.")
   }
-  params <- c('Horizon', 'Time', 'Variables', 'Region', 'Regrid', 'Workflow', 'Datasets')
+  params <- c('Horizon', 'Time', 'Variables', 'Region', 
+              'Regrid', 'Workflow', 'Datasets')
   if (!all(params %in% names(recipe$Analysis))) {
-    stop(paste("The element 'Analysis' in the recipe should contain these",
+    error(logger, 
+          paste("The element 'Analysis' in the recipe should contain these",
          "elements:", paste(params, collapse = " ")))
   }
   horizons <- c('Subseasonal', 'Seasonal', 'Decadal')
   if (!any(horizons %in% recipe$Analysis$Horizon)) {
-    stop("The element 'Horizon' in the recipe should be one of the followings:",
+    error(logger,
+         "The element 'Horizon' in the recipe should be one of the followings:",
          paste(horizons, collapse = " "))
   }
-  # more checks
+  # Check temporal settings and
+  # count the number of verifications
+  time_settings <- c('sdate', 'leadtimemin', 'leadtimemax', 'hcst_start', 'hcst_end')
+  if (!all(time_settings %in% names(recipe$Analysis$Time))) {
+    error(logger, 
+       paste("The element 'Time' in the recipe should contain these elements:",
+             paste(time_settings, collapse = " ")))
+  }
+  if (is.null(recipe$Analysis$Time$sdate$fcst_year) ||
+    recipe$Analysis$Time$sdate$fcst_year == 'None') {
+    stream <- "hindcast"
+    recipe$Analysis$Time$sdate$fcst_year <- 'YYYY'
+  } else {
+    stream <- "fcst"
+  }
+  if (length(recipe$Analysis$Time$sdate$fcst_day) > 1 &&
+      tolower(recipe$Analysis$Horizon) != "subseasonal") {
+    warn(logger,
+         paste("Only subseasonal verification allows multiple forecast days."),
+               "Element fcst_day in recipe set as 1.")
+    recipe$Analysis$Time$sdate$fcst_day <- '01'
+  }
+  if (is.null(recipe$Analysis$Time$sdate$fcst_month)) {
+    error(logger,
+          paste("The element 'fcst_month' in the recipe should be defined."))
+  }
+  fcst.sdate <- NULL
+  for (year in recipe$Analysis$Time$sdate$fcst_year) {
+    for (month in recipe$Analysis$Time$sdate$fcst_month) {
+      for (day in recipe$Analysis$Time$sdate$fcst_day) {
+        fcst.sdate <- c(fcst.sdate, 
+                             paste0(year,
+                               sprintf("%02d", as.numeric(month)),
+                               sprintf("%02d", as.numeric(day))))
+      }
+    }
+  }
+  fcst.sdate <- list(stream = stream, fcst.sdate = fcst.sdate)   
+  # Regrid checks:
+  if (length(recipe$Analysis$Regrid) != 2) {
+    error(logger,
+          "The 'Regrid' element should specified the 'method' and 'type'.")
+  }
+# more checks
   # ...
   # calculate number of workflows to create for each variable and 
   if (length(recipe$Analysis$Horizon) > 1) {
-     stop("Only 1 Horizon can be specified in the recipe")
+     error(logger, "Only 1 Horizon can be specified in the recipe")
   }
   nvar <- length(recipe$Analysis$Variables)
+  if (nvar > 2) {
+    error(logger,
+          "Only two type of Variables can be listed: ECVs and Indicators.")
+  }
+  # remove NULL or None Indicators or ECVs from the recipe:
+  if (!is.list(recipe$Analysis$Variables$Indicators)) {
+    recipe$Analysis$Variables <- recipe$Analysis$Variables[
+           -which(names(recipe$Analysis$Variables) == 'Indicators')]
+  }
+  if (!is.list(recipe$Analysis$Variables$ECVs)) {
+    recipe$Analysis$Variables <- recipe$Analysis$Variables[
+           -which(names(recipe$Analysis$Variables) == 'ECVs')]
+  }
   if (length(recipe$Analysis$Region) > 3) {
-    stop("The section 'Region' in the recipe should be checked.",
+    error(logger, "The section 'Region' in the recipe should be checked.",
          "Global, Countries and Regional elements are expected.")
   }
   if (!is.logical(recipe$Analysis$Region$Global)) {
-    stop("The element 'Global' in Region should be logical.")
+    error(logger, "The element 'Global' in Region should be logical.")
   }
   nregions <- 0
   if (recipe$Analysis$Region$Global == TRUE) {
@@ -43,8 +103,9 @@ check_recipe <- function(recipe, file, conf) {
     limits <- c('latmin', 'latmax', 'lonmin', 'lonmax')
     for (i in 1:length(recipe$Analysis$Region$Regional)) {
       if (!all(limits %in% names(recipe$Analysis$Region$Regional[[i]]))) {
-        stop("Each region defined in element 'Regional' should have 4 elements:",
-             paste(limits, collapse = " "))
+        error(logger,
+           "Each region defined in element 'Regional' should have 4 elements:",
+              paste(limits, collapse = " "))
       }
       # are numeric? class list mode list 
     }
@@ -52,8 +113,9 @@ check_recipe <- function(recipe, file, conf) {
   # Check workflow: need to define restrictions? 
   # e.g. only one calibration method
   nverifications <- check_number_of_dependent_verifications(recipe, conf)
-  message("Recipe checked succsessfully.")
-  return(nverifications)
+  info(logger, paste("Start Dates", paste(fcst.sdate, collapse = " ")))
+  info(logger, "Recipe checked succsessfully.")
+  return(append(nverifications, fcst.sdate))
 }
 
 check_number_of_dependent_verifications <- function(recipe, conf) {
@@ -64,79 +126,116 @@ check_number_of_dependent_verifications <- function(recipe, conf) {
   independent_verifications <- NULL
   dependent_verifications <- NULL
   dep <- 1
-  if (is.null(recipe$Analysis$Variables$Indicators) || 
-      (length(recipe$Analysis$Variables$Indicators) == 1 && 
-       is.null(recipe$Analysis$Variables$ECVs))) {
-  # look for duplications
-    vars <- unlist(recipe$Analysis$Variables)
-    if (any(duplicated(vars))) {
-      stop("There are duplicated variables.")
+  # check workflow order:
+  if (all(c('Calibration', 'Indicators') %in% names(recipe$Analysis$Workflow))) {
+    cal_pos <- which(names(recipe$Analysis$Workflow) == 'Calibration')
+    ind_pos <- which(names(recipe$Analysis$Workflow) == 'Indicators')
+    if (cal_pos < ind_pos) {
+      workflow_independent <- FALSE
+    } else {
+      workflow_independent <- TRUE
     }
-    independent_verifications <- as.vector(vars)
-  } else if (length(names(recipe$Analysis$Variables)) == 2 || 
-             length(recipe$Analysis$Variables$Indicators) > 1) {
-    ecvs <- recipe$Analysi$Variables$ECVs$name
-    inds <- recipe$Analysi$Variables$Indicators$name
-    ind_table <- read_yaml(paste0(conf$code_dir,
-                           "conf/vars-dict.yml"))
-    ind_table <- ind_table[unlist(inds)]
-    vars_for_ind <- NULL
-    # look for ecvs required to compute indicators
-    for (ind in unlist(inds)) {
-      if (is.null(ind_table[ind][[1]])) {
-        stop("Check the name of the indicators.")
-      }
-      vars_for_ind <- c(vars_for_ind, ind_table[ind][[1]]$ECVs)
-    }
-    if (any(vars_for_ind %in% unlist(ecvs)) || 
-        any(duplicated(vars_for_ind))) {
-      # check workflow order to know the number of verifications to calculate
-      if (all(c('Calibration', 'Indicators') %in% names(recipe$Analysis$Workflow))) {
-        cal_pos <- which(names(recipe$Analysis$Workflow) == 'Calibration')
-        ind_pos <- which(names(recipe$Analysis$Workflow) == 'Indicators')
-        if (cal_pos < ind_pos) {
-          workflow_independent <- FALSE
-        } else {
-          workflow_independent <- TRUE
+  }
+  if (workflow_independent) {
+    independent_verifications <- append(recipe$Analysis$Variables$ECVs,
+                                        recipe$Analysis$Variables$Indicators)
+  } else { 
+    if (is.null(recipe$Analysis$Variables$Indicators) || 
+       (length(recipe$Analysis$Variables$Indicators) == 1 && 
+        is.null(recipe$Analysis$Variables$ECVs))) {
+      independent_verifications <- append(recipe$Analysis$Variables$ECVs,
+                                          recipe$Analysis$Variables$Indicators)
+    } else { 
+      ecvs <- recipe$Analysi$Variables$ECVs
+      inds <- recipe$Analysi$Variables$Indicators
+      ind_table <- read_yaml(paste0(conf$code_dir,
+                             "conf/indicators_table.yml"))
+      # first, loop on ecvs if any and compare to indicators
+      done <- NULL # to gather the indicators reviewed
+      if (!is.null(ecvs)) {
+        for (i in 1:length(ecvs)) {
+          dependent <- list(ecvs[[i]])
+          for (j in 1:length(inds)) {
+            if (ind_table[inds[[j]]$name][[1]]$ECVs == ecvs[[i]]$name) {
+              if (ind_table[inds[[j]]$name][[1]]$freq == ecvs[[i]]$freq) {
+              # they are dependent
+                dependent <- append(dependent, inds[[j]])
+                done <- append(done, inds[[j]])
+              }
+            }
+          }
+          if (length(dependent) == 1) {
+            dependent <- NULL
+            independent_verifications <- append(independent_verifications,
+                                              list(ecvs[[i]]))
+          } else {
+            dependent_verifications <- append(dependent_verifications,
+                                              list(dependent))
+          }  
         }
-      }
-    }
-    vars <- unique(c(unlist(ecvs), unlist(vars_for_ind)))
-    for (var in unlist(vars)) {
-      if (var %in% vars_for_ind) {
-        if (workflow_independent) {
-          if (var %in% unlist(ecvs)) {
-            varname <- var
+        # There are indicators not reviewed yet?
+        if (length(done) < length(inds)) {
+          if (length(inds) == 1) {
+            independent_verifications <- append(independent_verifications,
+                                                inds)
           } else {
-            varname <- NULL # avoid writting a ecv not requested by needed by the ind
+            done <- NULL
+            for (i in 1:(length(inds) - 1)) {
+              dependent <- list(inds[[i]]$name)
+              if (is.na(match(unlist(dependent), unlist(done)))) {
+                for (j in (i+1):length(inds)) {
+                  if (ind_table[inds[[i]]$name][[1]]$ECVs == 
+                      ind_table[inds[[j]]$name][[1]]$ECVs) {
+                    if (ind_table[inds[[i]]$name][[1]]$freq ==
+                        ind_table[inds[[j]]$name][[1]]$freq) {
+                      dependent <- append(dependent, inds[[j]]$name)
+                      done <- dependent
+                    }
+                  }
+                }
+              } 
+              if (length(dependent) == 1) {
+                independent_verifications <- dependent
+                dependent <- NULL     
+              } else {
+                dependent_verifications <- dependent
+              }
+            } 
+          }
+        }
+      } else { # there are only Indicators:
+        done <- NULL
+        for (i in 1:(length(inds) - 1)) {
+          dependent <- list(inds[[i]]$name)
+          if (is.na(match(unlist(dependent), unlist(done)))) {
+            for (j in (i+1):length(inds)) {
+              if (ind_table[inds[[i]]$name][[1]]$ECVs ==
+                  ind_table[inds[[j]]$name][[1]]$ECVs) {
+                if (ind_table[inds[[i]]$name][[1]]$freq ==
+                    ind_table[inds[[j]]$name][[1]]$freq) {
+                  dependent <- append(dependent, inds[[j]]$name)
+                  done <- dependent
+                }
+              }
+            }
           }
-          independent_verifications <- c(independent_verifications, varname,
-                         as.vector(unlist(inds[which(unlist(vars_for_ind %in% var))])))
-        } else {
-          if (var %in% unlist(ecvs)) {
-            varname <- var
+          if (length(dependent) == 1) {
+            independent_verifications <- dependent
+            dependent <- NULL
           } else {
-            varname <- NULL
+            dependent_verifications <- dependent
           }
-          dependent_verifications[[dep]] <- c(varname,
-                         as.vector(unlist(inds[which(unlist(vars_for_ind %in% var))])))
-          dep = dep + 1
         }
-      } else {
-        independent_verifications <- var
-      }
+      } 
     }
-  } else {
-   stop("Only two types of 'Variables', named 'ECVs' and 'Indicators',",
-         " are allawed.")
   }
   if (!is.null(independent_verifications)) {   
-    message("The variables for independent verification are ", 
-            paste(independent_verifications, collapse = " "))
+    info(logger, paste("The variables for independent verification are ", 
+         paste(independent_verifications, collapse = " ")))
   }
   if (!is.null(dependent_verifications)) {
-    message("The variables for dependent verification are: ",
-            paste(dependent_verifications, collapse = " "))
+    info(logger, paste("The variables for dependent verification are: ",
+         paste(dependent_verifications, collapse = " ")))
   }
   # remove unnecessary names in objects to be removed
   return(list(independent = independent_verifications, 
diff --git a/tools/libs.R b/tools/libs.R
index 9ac904f17fb7af1a457af42ab29d2b45b2ec9ac2..8d9856d7697de1c89031bcd98a1dd8afb2a32163 100644
--- a/tools/libs.R
+++ b/tools/libs.R
@@ -1,10 +1,11 @@
-
+library(log4r)
+library(startR)
+library(ClimProjDiags)
+library(multiApply)
 # library(s2dverification)
 #  library(ncdf4)
 #  library(abind)
 #  library(easyVerification)
-#  library(multiApply)
-#  library(startR)
 #  library(easyNCDF)
 #  library(CSTools)
 #
@@ -21,6 +22,7 @@
 #  source("R_CST_MergeDims.R")
 #setwd(conf$code_dir)
 #  # To be removed when new package is done by library(CSOperational)
-  source("/esarchive/scratch/lpalma/git/auto-s2s/tools/check_recipe.R")
-  source("/esarchive/scratch/lpalma/git/auto-s2s/tools/check_conf.R")
-
+  source(paste0(conf$code_dir, "tools/check_recipe.R"))
+  source(paste0(conf$code_dir, "tools/check_conf.R"))
+  source(paste0(conf$code_dir, "tools/prepare_outputs.R"))
+  source(paste0(conf$code_dir, "tools/add_dims.R")) # Not sure if necessary yet
diff --git a/tools/prepare_outputs.R b/tools/prepare_outputs.R
new file mode 100644
index 0000000000000000000000000000000000000000..08f0063e00e6602ce6489c0e0909f27eccd8284a
--- /dev/null
+++ b/tools/prepare_outputs.R
@@ -0,0 +1,34 @@
+prepare_outputs <- function(recipe, file, conf) {
+# recipe: the content of the readed recipe
+# file: the recipe file name
+# conf: the content of the readed config user file
+  # Create output folders:
+  folder_name <- paste0(gsub(".yml", "", gsub("/", "_", file)), "_",
+    gsub(" ", "", gsub(":", "", gsub("-", "", Sys.time()))))
+  dir.create(file.path(conf$output_dir, folder_name, 'plots'), recursive = TRUE)
+  dir.create(file.path(conf$output_dir, folder_name, 'outputs'))
+  dir.create(file.path(conf$output_dir, folder_name, 'logs'))
+  file.copy(file, file.path(conf$output_dir, folder_name, 'logs'))
+  logfile <- file.path(conf$output_dir, folder_name, 'logs', 'log.txt')
+  # Set default behaviour of log output file:
+  if (is.null(recipe$Run)) {
+    recipe$Run <- list(Loglevel = 'INFO', Terminal = TRUE)
+  }
+  if (is.null(recipe$Run$Loglevel)) { 
+    recipe$Run$Loglevel <- 'INFO'
+  }
+  if (!is.logical(recipe$Run$Terminal)) {
+    recipe$Run$Terminal <- TRUE
+  }
+  if (recipe$Run$Terminal) {
+    logger <- logger(threshold = recipe$Run$Loglevel,  
+              appenders = list(console_appender(layout = default_log_layout()),
+                               file_appender(logfile, append = TRUE,   
+                               layout = default_log_layout())))
+  } else {
+    logger <- logger(threshold = recipe$Run$Loglevel,
+                appenders = list(file_appender(logfile, append = TRUE,             
+                               layout = default_log_layout())))
+  }
+  return(logger)
+}
diff --git a/tools/test_check_number_of_independent_verifications.R b/tools/test_check_number_of_independent_verifications.R
new file mode 100644
index 0000000000000000000000000000000000000000..846dc5be7697653b3b722157de288becfbc7329e
--- /dev/null
+++ b/tools/test_check_number_of_independent_verifications.R
@@ -0,0 +1,161 @@
+library(testthat)
+test_that("A few combinations", {
+ source("/esarchive/scratch/nperez/git/auto-s2s/tools/check_recipe.R")
+ recipe <- list(Analysis = 
+                list(Variables = list(ECVs = list(
+                                  list(name = 'tas', freq = 'daily_mean')), 
+                     Indicators = list(list(name = 'gdd'))),
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+ conf <- list(code_dir = "/esarchive/scratch/nperez/git/auto-s2s/")
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = NULL,
+                   dependent = list(list(list(name = 'tas', freq = 'daily_mean'),
+                                    name = 'gdd'))))
+ recipe <- list(Analysis = 
+                list(Variables = list(ECVs = list(
+                                 list(name = 'tas', freq = 'daily_mean')),
+                                      Indicators = list(name = 'gdd')),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC'))) 
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tas', freq = 'daily_mean'), 
+                               name = 'gdd'), dependent = NULL))
+
+ recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                          list(name = 'tas', freq = 'daily_mean'))),
+                     Workflow = list(Indicators = FALSE, Calibration = 'SBC')))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tas', freq = 'daily_mean')),
+                   dependent = NULL))
+
+ recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                         list(name = 'tas', freq = 'daily_mean'), 
+                         list(name = 'tasmax', freq = 'daily_mean'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tas', freq = 'daily_mean'),
+                                    list(name = 'tasmax', freq = 'daily_mean')),
+                   dependent = NULL))
+
+ recipe <- list(Analysis =
+                list(Variables = list(Indicators = list(list(name = 'gdd'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'gdd')), dependent = NULL))
+ 
+ recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                                   list(name = 'tas', freq = 'daily_mean')), 
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'))), 
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tas', freq = 'daily_mean'),
+                       list(name =  'gdd'), list(name = 'gst')), dependent = NULL))
+
+ recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                                    list(name = 'tas', freq = 'daily_mean')), 
+                                      Indicators = list(list(name = 'gdd'),
+                                                        list(name = 'gst'))),  
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = NULL,
+                   dependent = list(list(list(name = 'tas', freq = 'daily_mean'),
+                                         name = 'gdd', name = 'gst'))))
+# Dependent workflow: cal & ind
+  recipe <- list(Analysis =
+                list(Variables = list(
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'))),
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = NULL, dependent = list('gdd', 'gst')))
+
+  recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                                    list(name = 'tasmax', freq = 'daily_mean')),
+                                      Indicators = list(list(name = 'gdd'),
+                                                        list(name = 'gst'))),
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tasmax', freq = 'daily_mean')),
+                   dependent = list('gdd', 'gst')))
+
+  recipe <- list(Analysis =
+                 list(Variables = list(ECVs = list(
+                                    list(name = 'tasmax', freq = 'daily_mean'), 
+                                    list(name = 'tas', freq = 'daily_mean')),
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'))),
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tasmax', freq = 'daily_mean')),
+                   dependent = list(list(list(name = 'tas', freq = 'daily_mean'),
+                                     name = 'gdd', name = 'gst'))))
+
+  recipe <- list(Analysis =
+                 list(Variables = list(ECVs = list(
+                                  list(name = 'tasmax', freq = 'daily_mean'),
+                                  list(name = 'tas', freq = 'daily_mean')),
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'), 
+                                                        list(name = 'spr32'))),
+                     Workflow = list(Calibration = 'SBC', Indicators = TRUE)))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = NULL, 
+                   dependent = list(list(list(name = 'tasmax', freq = 'daily_mean'),
+                                         name = 'spr32'), 
+                                    list(list(name = 'tas', freq = 'daily_mean'),
+                                         name = 'gdd', name ='gst'))))
+# Independent workflow: ind & cal
+  recipe <- list(Analysis =
+                list(Variables = list(
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+ expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'gdd'), list(name = 'gst')),
+                   dependent = NULL))
+
+  recipe <- list(Analysis =
+                list(Variables = list(ECVs = list(
+                       list(name ='tasmax', freq = 'daily_freq')),
+                                      Indicators = list(list(name = 'gdd'), 
+                                                        list(name = 'gst'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tasmax', freq = 'daily_freq'),
+                                   list(name = 'gdd'), list(name = 'gst')),
+                   dependent = NULL))
+
+  recipe <- list(Analysis =
+                 list(Variables = list(ECVs = list(
+                           list(name = 'tasmax', freq = 'daily_mean'),
+                           list(name = 'tas', freq = 'daily_mean')),
+                                      Indicators = list(list(name = 'gdd'),
+                                                        list(name = 'gst'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tasmax', freq = 'daily_mean'),
+                                      list(name = 'tas', freq = 'daily_mean'),
+                                      list(name = 'gdd'), 
+                                      list(name = 'gst')), dependent = NULL))
+
+  recipe <- list(Analysis =
+                 list(Variables = list(ECVs = list(
+                                      list(name = 'tasmax', freq = 'daily_mean'),
+                                       list(name = 'tas', freq = 'daily_mean')),
+                                      Indicators = list(list(name = 'gdd'),
+                                                        list(name = 'gst'),
+                                                        list(name = 'spr32'))),
+                     Workflow = list(Indicators = TRUE, Calibration = 'SBC')))
+  expect_equal(check_number_of_dependent_verifications(recipe, conf),
+              list(independent = list(list(name = 'tasmax', freq = 'daily_mean'),
+                                      list(name = 'tas', freq = 'daily_mean'),
+                                      list(name = 'gdd'), list(name = 'gst'), 
+                                      list(name = 'spr32')),
+                   dependent = NULL))
+})
+