diff --git a/conf/archive.yml b/conf/archive.yml
index 88872cb85e4bda71134af98417c64720df05a28c..a0caf86aa4f19b0ec1f1cf9938a6b2657d758c76 100644
--- a/conf/archive.yml
+++ b/conf/archive.yml
@@ -16,7 +16,7 @@ esarchive:
                            "tasmin":"_f24h/", "tasmax":"_f24h/",
                            "ta300":"_f12h/", "ta500":"_f12h/", "ta850":"_f12h/",
                            "g300":"_f12h/", "g500":"_f12h/", "g850":"_f12h/",
-                           "tdps":"_f6h/", "tos":"_f6h/"}
+                           "tdps":"_f6h/", "psl":"_f6h/", "tos":"_f6h/"}
             nmember:
               fcst: 51
               hcst: 25
@@ -135,7 +135,7 @@ esarchive:
             daily_mean: {"tas":"_f1h-r1440x721cds/",
                          "rsds":"_f1h-r1440x721cds/",
                          "prlr":"_f1h-r1440x721cds/",
-                          "g300":"_f1h-r1440x721cds/",
+                         "g300":"_f1h-r1440x721cds/",
                          "g500":"_f1h-r1440x721cds/",
                          "g850":"_f1h-r1440x721cds/",
                          "sfcWind":"_f1h-r1440x721cds/",
@@ -146,6 +146,7 @@ esarchive:
                          "ta850":"_f1h-r1440x721cds/",
                          "hurs":"_f1h-r1440x721cds/"}
             monthly_mean: {"tas":"_f1h-r1440x721cds/",
+                           "psl":"_f1h-r1440x721cds/",
                            "prlr":"_f1h-r1440x721cds/",
                            "rsds":"_f1h-r1440x721cds/",
                            "g300":"_f1h-r1440x721cds/",
@@ -157,8 +158,7 @@ esarchive:
                            "ta300":"_f1h-r1440x721cds/",
                            "ta500":"_f1h-r1440x721cds/", 
                            "ta850":"_f1h-r1440x721cds/",
-                           "tos":"_f1h-r1440x721cds"
-            }
+                           "tos":"_f1h-r1440x721cds/"}
             calendar: "standard"
             reference_grid: "/esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h-r1440x721cds/tas_201805.nc"
         ERA5-Land:
diff --git a/example_scripts/execute_NAO.R b/example_scripts/execute_NAO.R
new file mode 100644
index 0000000000000000000000000000000000000000..39d435c543409b5d2b5b1fd4a7e8fe33cbdc4f5d
--- /dev/null
+++ b/example_scripts/execute_NAO.R
@@ -0,0 +1,24 @@
+
+recipe_file <- "recipes/examples/NAO_recipe.yml"
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("tools/prepare_outputs.R")
+source("modules/Anomalies/Anomalies.R")
+recipe <- prepare_outputs(recipe_file)
+
+#for (smonth in 1:12) {
+  data <- load_datasets(recipe)
+  gc()
+  data <- compute_anomalies(recipe, data)
+#  data <- readRDS("../Test_NAOmodule.RDS")
+source("modules/Indices/Indices.R")
+  nao_s2dv <- Indices(data = data, recipe = recipe)
+source("modules/Skill/Skill.R")
+  # todo parameter agg to get it from the nao_s2dv?
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nao_s2dv,
+                                         agg = 'region')
+
+
diff --git a/example_scripts/execute_NAO_decadal.R b/example_scripts/execute_NAO_decadal.R
new file mode 100644
index 0000000000000000000000000000000000000000..41570aa475c4d61fd7fdb98b2eaa12e57c9a9129
--- /dev/null
+++ b/example_scripts/execute_NAO_decadal.R
@@ -0,0 +1,24 @@
+recipe_file <- "recipes/examples/recipe_model_decadal_NAO.yml"
+
+
+source("modules/Loading/Loading_decadal.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("tools/prepare_outputs.R")
+source("modules/Anomalies/Anomalies.R")
+recipe <- prepare_outputs(recipe_file)
+
+#for (smonth in 1:12) {
+  data <- load_datasets(recipe)
+  gc()
+  data <- compute_anomalies(recipe, data)
+# data <- readRDS("../data_decadal_nao.RDS")
+source("modules/Indices/Indices.R")
+  nao_s2dv <- Indices(data = data, recipe = recipe)
+source("modules/Skill/Skill.R")
+  # todo parameter agg to get it from the nao_s2dv?
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nao_s2dv,
+                                         agg = 'region')
+
+
diff --git a/example_scripts/execute_Nino.R b/example_scripts/execute_Nino.R
new file mode 100644
index 0000000000000000000000000000000000000000..d64513ad29fd3045284e2590765b81ccded3865d
--- /dev/null
+++ b/example_scripts/execute_Nino.R
@@ -0,0 +1,34 @@
+recipe_file <- "recipes/examples/Nino_recipe.yml"
+
+
+source("modules/Loading/Loading.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("tools/prepare_outputs.R")
+source("modules/Anomalies/Anomalies.R")
+recipe <- prepare_outputs(recipe_file)
+
+#for (smonth in 1:12) {
+  data <- load_datasets(recipe)
+  gc()
+  data <- compute_anomalies(recipe, data)
+# saveRDS(data, file = "Test_Ninomodule.RDS")
+#  data <- readRDS("Test_Ninomodule.RDS")
+source("modules/Indices/Indices.R")
+  nino_s2dv <- Indices(data = data, recipe = recipe)
+source("modules/Skill/Skill.R")
+  # todo parameter agg to get it from the nao_s2dv?
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nino_s2dv$Nino1,
+                                         agg = 'region')
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nino_s2dv$Nino3,
+                                         agg = 'region')
+
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nino_s2dv$Nino4,
+                                         agg = 'region')
+
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nino_s2dv$Nino3.4,
+                                         agg = 'region')
+
+
+
diff --git a/example_scripts/execute_Nino_decadal.R b/example_scripts/execute_Nino_decadal.R
new file mode 100644
index 0000000000000000000000000000000000000000..1a7c3184a625a690fa3e98d98429547a9c250f32
--- /dev/null
+++ b/example_scripts/execute_Nino_decadal.R
@@ -0,0 +1,24 @@
+recipe_file <- "recipes/examples/recipe_Nino_decadal.yml"
+
+
+source("modules/Loading/Loading_decadal.R")
+source("modules/Calibration/Calibration.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("tools/prepare_outputs.R")
+source("modules/Anomalies/Anomalies.R")
+recipe <- prepare_outputs(recipe_file)
+
+#for (smonth in 1:12) {
+  data <- load_datasets(recipe)
+  gc()
+  data <- compute_anomalies(recipe, data)
+# data <- readRDS("../data_decadal_nao.RDS")
+source("modules/Indices/Indices.R")
+  nino_s2dv <- Indices(data = data, recipe = recipe)
+source("modules/Skill/Skill.R")
+  # todo parameter agg to get it from the nao_s2dv?
+  skill_metrics <- compute_skill_metrics(recipe = recipe, data = nino_s2dv[[1]],
+                                         agg = 'region')
+
+
diff --git a/modules/Indices/Indices.R b/modules/Indices/Indices.R
new file mode 100644
index 0000000000000000000000000000000000000000..fb9a7277587f5fb63f4f8151a090a6245f571a5b
--- /dev/null
+++ b/modules/Indices/Indices.R
@@ -0,0 +1,111 @@
+source("modules/Indices/R/compute_nao.R")
+source("modules/Indices/R/compute_nino.R")
+source("modules/Indices/R/drop_indices_dims.R")
+source("modules/Saving/Saving.R")
+Indices <- function(recipe, data) {
+  # Define parameters
+  nao <- NULL
+  if ("nao" %in% tolower(names(recipe$Analysis$Workflow$Indices))) {
+    if (!is.null(recipe$Analysis$Workflow$Indices$NAO$obsproj)) {
+      obsproj <- recipe$Analysis$Workflow$Indices$NAO$obsproj
+    } else {
+      obsproj <- TRUE
+    }
+    if (!is.null(recipe$Analysis$Workflow$Indices$NAO$plot_ts)) {
+      plot_ts <- recipe$Analysis$Workflow$Indices$NAO$plot_ts
+    } else {
+      plot_ts <- TRUE
+    }
+    if (!is.null(recipe$Analysis$Workflow$Indices$NAO$plot_sp)) {
+      plot_sp <- recipe$Analysis$Workflow$Indices$NAO$plot_sp
+    } else {
+      plot_sp <- TRUE
+    }
+    if (!is.null(recipe$Analysis$Workflow$Indices$NAO$alpha)) {
+      alpha <- recipe$Analysis$Workflow$Indices$NAO$alpha
+    } else {
+      alpha <- 0.05
+    }
+    nao <- compute_nao(data = data, recipe = recipe,
+                       obsproj = obsproj,
+                       plot_ts = plot_ts, plot_sp = plot_sp,
+                       alpha = alpha)
+  }
+  ninos <- list()
+  num_ninos <- sum(tolower(substr(names(recipe$Analysis$Workflow$Indices),
+                                  1, 4)) == "nino")
+  if (num_ninos > 0) { 
+    for (nins in 1:num_ninos) {
+      if ("nino1+2" %in% 
+          tolower(names(recipe$Analysis$Workflow$Indices)[nins])) {
+        region <- c(lonmin = -90, lonmax = -80, latmin = -10, latmax = 0)
+      } else if ("nino3" %in%
+           tolower(names(recipe$Analysis$Workflow$Indices)[nins])) {
+        region <- c(lonmin = -150, lonmax = -90, latmin = -5, latmax = 5)
+      } else if ("nino3.4" %in%
+           tolower(names(recipe$Analysis$Workflow$Indices)[nins])) {
+        region <- c(lonmin = -170, lonmax = -120, latmin = -5, latmax = 5)
+      } else if ("nino4" %in%
+           tolower(names(recipe$Analysis$Workflow$Indices)[nins])) {
+        region <- c(lonmin = 160, lonmax = -150, latmin = -5, latmax = 5)
+      }
+      if (is.null(recipe$Analysis$ncores)) {
+        ncores <- 1
+      } else {
+        ncores <- recipe$Analysis$ncores
+      }
+      if (is.null(recipe$Analysis$remove_NAs)) {
+        na.rm <- TRUE
+      } else {
+        na.rm <- recipe$Analysis$remove_NAs
+      }
+      if (is.null(recipe$Analysis$Workflow$Indices[nins]$standardised)) {
+        standardised <- TRUE
+      } else {
+        standardised <-recipe$Analysis$Workflow$Indices[nins]$standardised
+      }
+      if (is.null(recipe$Analysis$Workflow$Indices[nins]$running_mean)) {
+        running_mean <- NULL
+      } else {
+        running_mean <-recipe$Analysis$Workflow$Indices[nins]$running_mean
+      }
+      if (!is.null(recipe$Analysis$Workflow$Indices[nins]$plot_ts)) {
+        plot_ts <- recipe$Analysis$Workflow$Indices[nins]$plot_ts
+      } else {
+        plot_ts <- TRUE
+      }
+      if (!is.null(recipe$Analysis$Workflow$Indices[nins]$plot_sp)) {
+        plot_sp <- recipe$Analysis$Workflow$Indices[nins]$plot_sp
+      } else {
+        plot_sp <- TRUE
+      }
+      if (!is.null(recipe$Analysis$Workflow$Indices[nins]$alpha)) {
+        alpha <- recipe$Analysis$Workflow$Indices[nins]$alpha
+      } else {
+        alpha <- 0.05
+      }
+      if (!is.null(recipe$Analysis$Workflow$Indices[nins]$save)) {
+        save <- recipe$Analysis$Workflow$Indices[nins]$save
+      } else {
+        save <- 'all'
+      }         
+      nino <- compute_nino(data = data, recipe = recipe, 
+                           standardised = standardised,
+                           region = region, save = save, plot_ts = plot_ts,
+                           plot_sp = plot_sp, alpha = alpha,
+                           running_mean = running_mean, ncores = ncores,
+                           na.rm = na.rm)
+      ninos[[nins]] <- nino
+      names(ninos)[nins] <- names(recipe$Analysis$Workflow$Indices)[nins]    
+    } 
+  }
+  if (is.null(nao)) {
+    return(ninos)
+    info(recipe$Run$logger,
+         "##### EL NINO INDICES COMPUTED SUCCESSFULLY #####")
+  } else {
+    return(nao)
+    info(recipe$Run$logger,
+         "##### NAO INDEX COMPUTED SUCCESSFULLY #####")
+  }
+}
diff --git a/modules/Indices/R/compute_nao.R b/modules/Indices/R/compute_nao.R
new file mode 100644
index 0000000000000000000000000000000000000000..4be27b8932415c9407a1f2df5833e2b45cf60554
--- /dev/null
+++ b/modules/Indices/R/compute_nao.R
@@ -0,0 +1,444 @@
+
+compute_nao <- function(data, recipe, obsproj, plot_ts, plot_sp,
+                        alpha, logo = NULL) {
+  ## TODO: if fcst object in data compute the nao too
+  if (!is.null(data$fcst)) {
+    warning("NAO computed only for hindcast data.")
+  }
+  # Check if subsetting the data for the region is required
+  lons <- data$hcst$coords$longitude
+  lats <- data$hcst$coords$latitude
+  subset_needed <- FALSE
+  nao_region <- c(lonmin = -80, lonmax = 40,
+                  latmin = 20, latmax = 80)
+  if (any(lons < 0)) {  #[-180, 180] 
+    if (!(min(lons) > -90 & min(lons) < -70 & 
+          max(lons) < 50 & max(lons) > 30)) {
+      subset_needed <- TRUE
+    }
+  } else {  #[0, 360]
+    if (any(lons >= 50 & lons <= 270)) {
+      susbset_needed <- TRUE
+    } else {
+      lon_E <- lons[which(lons < 50)]
+      lon_W <- lons[-which(lons < 50)]
+      if (max(lon_E) < 30 | min(lon_W) > 290) {
+        subset_needed <- TRUE
+      }
+    }
+  }
+  if (any(max(lats) > 80 & min(lats) < 20)) {
+    subset_needed <- TRUE
+  }
+  if (subset_needed) {
+    warning("The data is being subsetted for 20N-80N and 80W-40E region.")
+
+      hcst1 <- ClimProjDiags::SelBox(data$hcst$data,
+                                     lon = as.vector(lons),
+                                     lat = as.vector(lats),
+                                     region = nao_region,
+                                     londim = "longitude",
+                                     latdim = "latitude")
+      obs1 <- ClimProjDiags::SelBox(data$obs$data,
+                                    lon = as.vector(lons),
+                                    lat = as.vector(lats),
+                                    region = nao_region,
+                                    londim = "longitude",
+                                    latdim = "latitude")
+      hcst <- s2dv_cube(data = hcst1$data, lat = hcst1$lat, lon = hcst1$lon,
+                        Variable = c(data$hcst$Variable[1], level = 'surface'),
+                        data$hcst$Va, Datasets = data$hcst$Datasets,
+                        time_dims = c('syear', 'time'),
+                        Dates = data$hcst$Dates)
+      obs <- s2dv_cube(data = obs1$data, lat = obs1$lat, lon = obs1$lon,
+                       Variable = c(data$obs$Variable[1], level = 'surface'),
+                       Datasets = data$obs$Datasets, time_dims = c('syear', 'time'))
+# TODO check and create data object for the next steps
+     data <- list(hcst = hcst, obs = obs)
+     lons <- data$hcst$coords$longitude
+     lats <- data$hcst$coords$latitude
+     obs1 <- hcst1 <- NULL
+     gc()
+  }
+  nao <- NAO(exp = data$hcst$data, obs = data$obs$data,
+             lat = data$hcst$coords$latitude,
+             lon = data$hcst$coords$longitude,
+             time_dim = 'syear',
+             memb_dim = 'ensemble', space_dim = c('latitude', 'longitude'),
+             ftime_dim = 'time', ftime_avg = NULL,
+             obsproj = obsproj, ncores = recipe$Analysis$ncores)
+  ## Standardisation:
+  nao$exp <- Apply(list(nao$exp), target_dims = c('syear', 'ensemble'),
+                   fun = function(x) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(x, c(1,2), function(x) {(x-means)/sd})},
+                   ncores = recipe$Analysis$ncores)$output1
+  nao$obs <- Apply(list(nao$obs), target_dims = c('syear', 'ensemble'),
+                   fun = function(x) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(x, c(1,2), function(x) {(x-means)/sd})},
+                   ncores = recipe$Analysis$ncores)$output1
+  nao$exp <- InsertDim(nao$exp, posdim = 1, lendim = 1, name = 'region')
+  nao$obs <- InsertDim(nao$obs, posdim = 1, lendim = 1, name = 'region')
+  hcst_dates <- data$hcst$attrs$Dates
+  hcst_dates <- drop(data$hcst$attrs$Dates)
+
+  if (!("time" %in% names(dim(hcst_dates)))) {
+     if (is.null(dim(hcst_dates))) {
+       hcst_dates <- array(hcst_dates, c(syear = length(hcst_dates)))
+     }
+     hcst_dates <- InsertDim(hcst_dates, pos = 1, len = 1, name = 'time')
+     hcst_dates <- as.POSIXct(hcst_dates, origin = '1970-01-01', tz = 'UTC')
+  }
+  nao <- list(hcst = s2dv_cube(
+                       data = nao$exp, 
+                       varName = "nao",
+                       metadata = list(
+                         region = list(name = "NAO region",
+                                      lats_range = paste0(range(lats)),
+                                      lons_range = paste0(range(lons))),
+                         time = data$hcst$attrs$Variable$metadata$time,
+                       nao = list(units = 'adim', 
+                                  longname = 'North Atlantic Oscillation')),
+                       Dates = hcst_dates,
+                       Dataset = recipe$Analysis$Datasets$System$name),
+               obs = s2dv_cube(
+                       data = nao$obs, 
+                       varName = "nao",
+                       metadata = list(
+                         region = list(name = "NAO region",
+                                       lats_range = paste0(range(lats)),
+                                       lons_range = paste0(range(lons))),
+                         time = data$obs$attrs$Variable$metadata$time,
+                         nao = list(units = 'adim', 
+                                    longname = 'North Atlantic Oscillation')),
+                         Dates = data$obs$attrs$Dates,
+                         Dataset = recipe$Analysis$Datasets$Reference$name))
+  if (recipe$Analysis$Workflow$Indices$NAO$save == 'all') {
+    file_dest <- paste0(recipe$Run$output_dir, "/outputs/Indices/")
+    if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+      # Use startdates param from SaveExp to correctly name the files:
+      if (length(data$hcst$attrs$source_files) == dim(data$hcst$data)['syear']) {
+        file_dates <- Apply(data$hcst$attrs$source_files, target_dim = NULL,
+                            fun = function(x) {
+                              pos <- which(strsplit(x, "")[[1]] == ".")
+                              res <- substr(x, pos-8, pos-1)
+                            })$output1
+      }
+    } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+      file_dates <- paste0('s',
+             recipe$Analysis$Time$hcst_start : recipe$Analysis$Time$hcst_end)
+    }
+    # need to recover original dimensions after saving to make Skill module work
+    nao_original_dims_hcst <- nao$hcst$data
+    nao$hcst$data <- .drop_indices_dims(nao$hcst$data)
+    CST_SaveExp(data = nao$hcst,
+          destination = file_dest, 
+          startdates = as.vector(file_dates),
+          dat_dim = NULL, sdate_dim = 'syear', 
+          ftime_dim = 'time', var_dim = NULL, 
+          memb_dim = 'ensemble')
+    nao_original_dims_obs <- nao$obs$data
+    nao$obs$data <- .drop_indices_dims(nao$obs$data)
+    CST_SaveExp(data = nao$obs, #res,
+          destination = file_dest,
+          startdates = as.vector(file_dates),
+          dat_dim = NULL, sdate_dim = 'syear',
+          ftime_dim = 'time', var_dim = NULL,
+          memb_dim = NULL)
+    nao$hcst$data <- nao_original_dims_hcst
+    nao$obs$data <- nao_original_dims_obs
+    nao_original_dims_hcst <- nao_original_dims_obs <- NULL
+    gc()
+  }
+  # Read variable long_name to plot it
+  conf <- yaml::read_yaml("conf/variable-dictionary.yml")
+  var_name <- conf$vars[[which(names(conf$vars) ==
+                               recipe$Analysis$Variables$name)]]$long
+
+  if (plot_ts) {
+    dir.create(paste0(recipe$Run$output_dir, "/plots/Indices/"),
+               showWarnings = F, recursive = T)
+    source("modules/Indices/R/plot_deterministic_forecast.R")
+    for (tstep in 1:dim(nao$hcst$data)['time']) {
+        mes <- as.numeric(substr(recipe$Analysis$Time$sdate, 1,2)) +
+               (tstep - 1) + (recipe$Analysis$Time$ftime_min - 1)
+        mes <- ifelse(mes > 12, mes - 12, mes)
+      fmonth <- sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min)
+      obs <- Subset(nao$obs$data, along = 'time', ind = tstep)
+      exp <- Subset(nao$hcst$data, along = 'time', ind = tstep) 
+      if (gsub(".", "", recipe$Analysis$Datasets$System$name) == "") {
+        system <- recipe$Analysis$Datasets$System$name
+      } else {
+        system <-gsub(".", "", recipe$Analysis$Datasets$System$name)
+      }
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        toptitle <- paste("NAO Index\n",  
+                          month.abb[mes],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end, "\nMethod:",
+                          ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                          "Pobs", "Pmod"), "(Doblas-Reyes et al., 2003)")
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_",
+                           system, "_", recipe$Analysis$Datasets$Reference$name,
+                           "_s", recipe$Analysis$Time$sdate, "_ftime",
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png") 
+        caption <- paste0("NAO method: ",
+                          ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                 "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                          "Nominal start date: 1st of ",
+                          month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                          "\n",
+                          "Forecast month: ", fmonth, "\n")
+        xlabs <- as.numeric(substr(file_dates, 1, 4))
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+        toptitle <- paste("NAO Index\n",  
+                          "Lead time", fmonth,
+                          " / Start dates", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end, "\nMethod:",
+                          ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                 "Pobs", "Pmod"), "(Doblas-Reyes et al., 2003)")
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_",
+                           system, "_", recipe$Analysis$Datasets$Reference$name,
+                           "_ftime",
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png") 
+        caption <- paste0("NAO method: ",
+                          ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                 "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                          "Start date month: ",
+                          month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                          "\n",
+                          "Lead time: ", fmonth, "\n")
+        xlabs <- file_dates
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      plot_deterministic_forecast(obs, exp, 
+                                  time_dim = 'syear',
+                                  member_dim = 'ensemble', style = 'boxplot',
+                                  xlabs = xlabs,
+                                  title = toptitle, fileout = plotfile,
+                                  caption = caption, caption_line = 6.5,
+                                  legend_text = c(
+                                  recipe$Analysis$Datasets$Reference$name,
+                                  recipe$Analysis$Datasets$System$name))
+    }
+  }
+  if (plot_sp) {
+    ## TODO: To be removed when s2dv is released:
+    source("modules/Visualization/R/tmp/PlotRobinson.R")
+    source("modules/Indices/R/correlation_eno.R")
+    source("modules/Visualization/R/get_proj_code.R")
+    dir.create(paste0(recipe$Run$output_dir, "/plots/Indices/"),
+               showWarnings = F, recursive = T)
+    # Get correct code for stereographic projection
+    projection_code <- get_proj_code(proj_name = "stereographic")
+    correl_obs <- Apply(list(data$obs$data, nao$obs$data), 
+                        target_dims = 'syear', fun = .correlation_eno, 
+                        time_dim = 'syear', method = 'pearson', alpha = alpha, 
+                        test.type = 'two-sided', pval = FALSE, 
+                        ncores = recipe$Analysis$ncores)
+    correl_hcst <- Apply(list(data$hcst$data, nao$hcst$data), 
+                         target_dims = c('syear', 'ensemble'),
+                         fun = function(x, y) {
+                         x <- apply(x, 1, mean, na.rm = TRUE)
+                         y <- apply(y, 1, mean, na.rm = TRUE)
+                         dim(y) <- c(syear = length(y))
+                         dim(x) <- c(syear = length(x))
+                         res <- .correlation_eno(x, y,
+                           time_dim = 'syear', method = 'pearson', alpha = alpha, 
+                           test.type = 'two-sided', pval = FALSE)}, 
+                           ncores = recipe$Analysis$ncores)
+    correl_hcst_full <- Apply(list(data$hcst$data, nao$hcst$data), 
+                              target_dims = c('syear', 'ensemble'), 
+                              fun = function(x,y) {
+                              dim(y) <- c(syear = length(y))
+                              dim(x) <- c(syear = length(x))
+                           res <- .correlation_eno(x, y,
+                              time_dim = 'syear', method = 'pearson', alpha = alpha, 
+                              test.type = 'two-sided', pval = FALSE)}, 
+                                ncores = recipe$Analysis$ncores)
+    
+    for (tstep in 1:dim(nao$obs$data)['time']) {
+      fmonth <- sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min)
+      ## Observations
+      map <- drop(Subset(correl_obs$r, along = 'time', ind = tstep))
+      sig <- drop(Subset(correl_obs$sign, along = 'time', ind = tstep))
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        mes <- as.numeric(substr(recipe$Analysis$Time$sdate, 1,2)) +
+               (tstep - 1) + (recipe$Analysis$Time$ftime_min - 1)
+        mes <- ifelse(mes > 12, mes - 12, mes)
+        toptitle <- paste(recipe$Analysis$Datasets$Reference$name, "\n",
+                          "NAO Index -",var_name, "\n",
+                          " Correlation /",  month.abb[mes], 
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_",
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Nominal start date: 1st of ",
+                           month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
+        toptitle <- paste(recipe$Analysis$Datasets$Reference$name, "\n",
+                          "NAO Index -",var_name, "\n",
+                          " Correlation / Start dates ", 
+                          recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_",
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Start date: month ",
+                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      if (gsub(".", "", recipe$Analysis$Datasets$System$name) == "") {
+        system <- recipe$Analysis$Datasets$System$name
+      } else {
+        system <- gsub(".", "", recipe$Analysis$Datasets$System$name)
+      }
+
+      PlotRobinson(map, lon = lons, lat = lats, target_proj = projection_code,
+                   lat_dim = 'latitude', lon_dim = 'longitude',
+                   legend = 's2dv', style = 'polygon',
+                   toptitle = toptitle, crop_coastlines = nao_region,
+                   caption = caption, mask = sig, bar_extra_margin = c(4,0,4,0),
+                   fileout = plotfile, width = 8, height = 6,
+                   brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+      ## Ensemble-mean
+      map <- drop(Subset(correl_hcst$r, along = 'time', ind = tstep))
+      sig <- drop(Subset(correl_hcst$sign, along = 'time', ind = tstep))
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        toptitle <- paste(recipe$Analysis$Datasets$System$name, "\n",
+                          "NAO Index -", var_name, "\n",
+                          " Correlation /", month.abb[mes],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_ensmean_",
+                           system,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Correlation ensemble mean\n",
+                           "Nominal start date: 1st of ",
+                           month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+        toptitle <- paste(recipe$Analysis$Datasets$System$name,"\n",
+                          "NAO Index -",var_name, "\n",
+                          " Correlation / Start dates ", 
+                          recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_ensmean_",
+                           system, 
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Correlation ensemble mean\n",
+                           "Start date month: ",
+                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+
+      PlotRobinson(map, lon = lons, lat = lats, target_proj = projection_code,
+                   lat_dim = 'latitude', lon_dim = 'longitude',
+                   legend = 's2dv', bar_extra_margin = c(4,0,4,0),
+                   toptitle = toptitle, style = 'polygon',
+                   caption = caption, mask = sig, crop_coastline = nao_region,
+                   fileout = plotfile, width = 8, height = 6,
+                   brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+      
+      # Full hcst corr
+      map <- drop(Subset(correl_hcst_full$r, along = 'time', ind = tstep))
+      sig <- drop(Subset(correl_hcst_full$sign, along = 'time', ind = tstep))
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        toptitle <- paste(recipe$Analysis$Datasets$System$name,"\n",
+                          "NAO Index -",var_name, "\n",
+                          " Correlation /", month.abb[mes],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_member_",
+                           system,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Correlation all members\n",
+                           "Nominal start date: 1st of ",
+                           month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+        toptitle <- paste(recipe$Analysis$Datasets$System$name,"\n",
+                          "NAO Index -",var_name, "\n",
+                          " Correlation / Start dates ", 
+                          recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/NAO_correlation_",
+                           recipe$Analysis$Variable$name, "_member_",
+                           system, 
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("NAO method: ",
+                           ifelse(recipe$Analysis$Workflow$Indices$NAO$obsproj,
+                                  "Pobs", "Pmod"), " (Doblas-Reyes et al., 2003)\n",
+                           "Correlation all members\n",
+                           "Start date month: ",
+                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      PlotRobinson(map, lon = lons, lat = lats, target_proj = projection_code,
+                   lat_dim = 'latitude', lon_dim = 'longitude',
+                   legend = 's2dv', bar_extra_margin = c(4,0,4,0),
+                   toptitle = toptitle, style = 'polygon',
+                   caption = caption, mask = sig, crop_coastline = nao_region,
+                   fileout = plotfile, width = 8, height = 6,
+                   brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+    } # end tstep loop
+  }
+  return(nao)
+}
diff --git a/modules/Indices/R/compute_nino.R b/modules/Indices/R/compute_nino.R
new file mode 100644
index 0000000000000000000000000000000000000000..7ed50c65c73d1e7eb3f5f1e3394c9b01a6f049e4
--- /dev/null
+++ b/modules/Indices/R/compute_nino.R
@@ -0,0 +1,438 @@
+compute_nino <- function(data, recipe, region, standardised = TRUE,
+                         running_mean = NULL, ftime_avg = NULL,
+                         plot_ts = TRUE, plot_sp = TRUE, alpha = 0.5,
+                         save = TRUE, na.rm = TRUE, logo = NULL,
+                         ncores = NULL) {
+  #ftime_avg = 1:2
+  if (!is.null(data$fcst)) {
+    warn(recipe$Run$logger, "Nino computed only for hindcast data.")
+  }
+  var <- recipe$Analysis$Variables$name
+  if (!(var %in% c('tos', 'sst'))) {
+    warn(recipe$Run$logger, "Variable name is not one of the expected sst or tos")
+  }
+  var_units <- data$hcst$attrs$Variable$metadata[[var]]$units
+  nino_hcst <- WeightedMean(data = data$hcst$data, 
+                       lon = as.vector(data$hcst$coords$longitude),
+                       lat = as.vector(data$hcst$coords$latitude),
+                       region = region,
+                       londim = 'longitude',
+                       latdim = 'latitude',
+                       na.rm = na.rm,
+                       ncores = ncores)
+  nino_obs <- WeightedMean(data = data$obs$data,
+                       lon = as.vector(data$hcst$coords$longitude),
+                       lat = as.vector(data$hcst$coords$latitude),
+                       region = region,
+                       londim = 'longitude',
+                       latdim = 'latitude',
+                       na.rm = na.rm,
+                       ncores = ncores)
+  if (standardised) {
+  nino_hcst <- Apply(list(nino_hcst), target_dims = c('syear', 'ensemble'),
+                   fun = function(x) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(x, c(1,2), function(x) {(x-means)/sd})},
+                   ncores = recipe$Analysis$ncores)$output1
+  nino_obs <- Apply(list(nino_obs), target_dims = c('syear', 'ensemble'),
+                   fun = function(x) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(x, c(1,2), function(x) {(x-means)/sd})},
+                   ncores = recipe$Analysis$ncores)$output1
+    var_units <- 'adim'
+  }
+  if (!is.null(running_mean)) {
+    nino_hcst <- Smoothing(nino_hcst, 
+                           runmeanlen = running_mean,
+                           time_dim = 'time',
+                           ncores = ncores)
+    nino_obs <- Smoothing(nino_obs,
+                          runmeanlen = running_mean,
+                          time_dim = 'time',
+                          ncores = ncores)
+  }
+  if (!is.null(ftime_avg)) {
+    if (recipe$Analysis$Variable$freq == 'monthly_mean') {
+      monini <- as.numeric(substr(recipe$Analysis$Time$sdate, 1, 2))
+      moninf <- ftime_avg[1]
+      monsup <- ftime_avg[length(ftime_avg)]
+      nino_hcst <- Season(nino_hcst, posdim = 'time', 
+                          monini = monini,
+                          moninf = moninf, monsup = monsup)
+      nino_obs <- Season(nino_obs, posdim = 'time',
+                          monini = monini,
+                          moninf = moninf,
+                          monsup = monsup)
+    } else {
+      warning("'ftime_avg' only available for monthly data.")
+    }
+  }
+  if (all(region == c(-90, -80, -10, 0))) {
+    region_name <- "1+2"
+    nino_name <- "nino12"
+  } else if (all(region == c(-150, -90, -5, 5))) {
+    region_name <- "3"
+    nino_name <- "nino3"
+  } else if (all(region == c(-170, -120, -5, 5))) {
+    region_name <- "3.4"
+    nino_name <- "nino34"
+  } else if (all(region == c(160, -150, -5, 5))) {
+    region_name <- "4"
+    nino_name <- "nino4"
+  } else {
+    stop("Unknown nino region")
+  }
+  nino_hcst <- InsertDim(nino_hcst, posdim = 1, lendim = 1, name = 'region')
+  nino_obs <- InsertDim(nino_obs, posdim = 1, lendim = 1, name = 'region')
+  dims_dates_not_null <- dim(data$hcst$attrs$Dates)[which(dim(data$hcst$attrs$Dates) > 1)]
+  hcst_dates <- Subset(data$hcst$attrs$Dates, along = names(dims_dates_not_null),
+                       indices = lapply(dims_dates_not_null, function(x){1:x}),
+                       drop = "non-selected")	
+   if (!("time" %in% names(dim(hcst_dates)))) {
+       hcst_dates <- InsertDim(hcst_dates, pos = 1, len = 1, name = 'time')
+       ## TODO: recover dates format
+       hcst_dates <- as.POSIXct(hcst_dates, origin = '1970-01-01', tz = 'UTC')
+    }
+  nino <- list(hcst = s2dv_cube(
+                      data = nino_hcst,
+                      varName = nino_name,
+                      metadata = list(
+                      region = list(name = paste("Nino", region_name, "region"),
+                                    lats_range = paste(region[3:4]),
+                                    lons_range = paste(region[1:2])),
+                      time = data$hcst$attrs$Variable$metadata$time,
+                      nino = list(units = var_units,
+                             longname = paste("El Niño", region_name, "Index"))),
+                      Dates = hcst_dates,
+                      Dataset = recipe$Analysis$Datasets$System$name),
+               obs = s2dv_cube(
+                     data = nino_obs,
+                     varName = nino_name,
+                     metadata = list(
+                     region = list(name = paste("Nino", region_name, "region"),
+                                   lats_range = paste(region[3:4]),
+                                   lons_range = paste(region[1:2])),
+                     time = data$obs$attrs$Variable$metadata$time,
+                     nino = list(units = var_units,
+                            longname = paste("El Niño", region_name, "Index"))),
+                     Dates = data$obs$attrs$Dates,
+                     Dataset = recipe$Analysis$Datasets$Reference$name))
+
+  if (save) {
+    file_dest <- paste0(recipe$Run$output_dir, "/outputs/Indices/")
+    # Use startdates param from SaveExp to correctly name the files:
+    if (tolower(recipe$Analysis$Horizon) == "seasonal") { 
+      if (length(data$hcst$attrs$source_files) == dim(data$hcst$data)['syear']) {
+        file_dates <- Apply(data$hcst$attrs$source_files, target_dim = NULL,
+                            fun = function(x) {
+                            pos <- which(strsplit(x, "")[[1]] == ".")
+                            res <- substr(x, pos-8, pos-1)
+                            })$output1
+      }
+    } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+      file_dates <- paste0('s',recipe$Analysis$Time$hcst_start : recipe$Analysis$Time$hcst_end)
+    }
+    nino$hcst$data <- .drop_indices_dims(nino_hcst)
+    CST_SaveExp(data = nino$hcst,
+          destination = file_dest,
+          startdates = as.vector(file_dates),
+          dat_dim = NULL, sdate_dim = 'syear',
+          ftime_dim = 'time', var_dim = NULL,
+          memb_dim = 'ensemble')
+    res <- .drop_indices_dims(nino_obs)
+    if (!("time" %in% names(dim(res)))) {
+      res <- InsertDim(res, pos = 1, len = 1, name = 'time')
+    }
+    nino$obs$data <- res
+    CST_SaveExp(data = nino$obs,
+          destination = file_dest,
+          startdates = as.vector(file_dates),
+          dat_dim = NULL, sdate_dim = 'syear',
+          ftime_dim = 'time', var_dim = NULL,
+          memb_dim = NULL)
+    nino$hcst$data <- nino_hcst
+    nino$obs$data <- nino_obs
+    res <- NULL
+    gc()
+  }
+ # Read variable long_name to plot it
+  conf <- yaml::read_yaml("conf/variable-dictionary.yml")
+  var_name <- conf$vars[[which(names(conf$vars) ==
+                               recipe$Analysis$Variables$name)]]$long
+  if (plot_ts) {
+    dir.create(paste0(recipe$Run$output_dir, "/plots/Indices/"),
+               showWarnings = F, recursive = T)
+    source("modules/Indices/R/plot_deterministic_forecast.R")
+    for (tstep in 1:dim(nino$hcst$data)['time']) {
+      mes <- as.numeric(substr(recipe$Analysis$Time$sdate, 1,2)) +
+               (tstep - 1) + (recipe$Analysis$Time$ftime_min - 1)
+      mes <- ifelse(mes > 12, mes - 12, mes)
+      fmonth <- sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min)
+      obs <- Subset(nino$obs$data, along = 'time', ind = tstep)
+      exp <- Subset(nino$hcst$data, along = 'time', ind = tstep)
+      if (gsub(".", "", recipe$Analysis$Datasets$System$name) == "") {
+        system <- recipe$Analysis$Datasets$System$name
+      } else {
+        system <- gsub(".", "", recipe$Analysis$Datasets$System$name)
+      }
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        toptitle <- paste("Ni\u00F1o", region_name, "SST Index\n",
+                          month.abb[mes],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                           nino_name, "_",
+                           system, "_", recipe$Analysis$Datasets$Reference$name,
+                           "_s", recipe$Analysis$Time$sdate, "_ftime",
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), 
+                           ".png")
+        caption <-  paste0("Nominal start date: 1st of ",
+                           month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                           "\n",
+                           "Forecast month: ", fmonth)
+        xlabs = as.numeric(substr(file_dates, 1, 4))
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
+        toptitle <- paste("Ni\u00F1o", region_name, "SST Index\n",
+                          "Lead time", fmonth,
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/", nino_name, "_",
+                           system, "_", recipe$Analysis$Datasets$Reference$name,
+                           "_ftime",
+                           sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min), ".png")
+       caption <- paste0("Start date month: ",
+                          month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                          "\n",
+                          "Lead time: ", fmonth, "\n")
+        xlabs <- substr(file_dates, 2,5)
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      plot_deterministic_forecast(obs, exp,
+                                  time_dim = 'syear',
+                                  member_dim = 'ensemble', style = 'boxplot',
+                                  xlabs = xlabs,
+                                  ylab = var_units, 
+                                  title = toptitle, fileout = plotfile,
+                                  caption =  caption,
+                                  legend_text = c(
+                                  recipe$Analysis$Datasets$Reference$name,
+                                  recipe$Analysis$Datasets$System$name))
+    }
+  }
+  if (plot_sp) {
+    source("modules/Visualization/R/tmp/PlotRobinson.R")
+    source("modules/Indices/R/correlation_eno.R")
+    source("modules/Visualization/R/get_proj_code.R")
+    lons <- data$hcst$coords$longitude
+    lats <- data$hcst$coords$latitude
+    # Get code for Robinson projection depending on GEOS/GDAL/PROJ version
+    projection_code <- get_proj_code("robinson")
+    # Avoid rewriten longitudes a shift is neeced:
+    dir.create(paste0(recipe$Run$output_dir, "/plots/Indices/"),
+               showWarnings = F, recursive = T)
+    correl_obs <- Apply(list(data$obs$data, nino$obs$data),
+                        target_dims = 'syear',
+                        fun = .correlation_eno,
+                        time_dim = 'syear', method = 'pearson', alpha = alpha,
+                        test.type = 'two-sided', pval = FALSE,
+                        ncores = recipe$Analysis$ncores)
+    correl_hcst <- Apply(list(data$hcst$data, nino$hcst$data),
+                         target_dims = c('syear', 'ensemble'),
+                         fun = function(x, y) {
+                         x <- apply(x, 1, mean, na.rm = TRUE)
+                         y <- apply(y, 1, mean, na.rm = TRUE)
+                         dim(y) <- c(syear = length(y))
+                         dim(x) <- c(syear = length(x))
+                         res <- .correlation_eno(x, y,
+                           time_dim = 'syear', method = 'pearson', alpha = alpha,
+                           test.type = 'two-sided', pval = FALSE)},
+                           ncores = recipe$Analysis$ncores)
+    correl_hcst_full <- Apply(list(data$hcst$data, nino$hcst$data),
+                              target_dims = c('syear', 'ensemble'),
+                              fun = function(x,y) {
+                              dim(y) <- c(syear = length(y))
+                              dim(x) <- c(syear = length(x))
+                           res <- .correlation_eno(x, y,
+                              time_dim = 'syear', method = 'pearson', alpha = alpha,
+                              test.type = 'two-sided', pval = FALSE)},
+                              ncores = recipe$Analysis$ncores)
+    for (tstep in 1:dim(nino$obs$data)['time']) {
+      map <- Subset(correl_obs$r, along = 'time', ind = tstep, drop = T)
+      sig <- Subset(correl_obs$sig, along = 'time', ind = tstep, drop = T)
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        mes <- as.numeric(substr(recipe$Analysis$Time$sdate, 1,2)) +
+               (tstep - 1) + (recipe$Analysis$Time$ftime_min - 1)
+        mes <- ifelse(mes > 12, mes - 12, mes)
+        fmonth <- sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min)
+        toptitle <- paste(recipe$Analysis$Datasets$Reference$name, "\n",
+                          "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                          " Correlation /",
+                          month.abb[mes],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/", nino_name,
+                            "_correlation_",
+                           recipe$Analysis$Variable$name, "_",
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("Nominal start date: 1st of ",
+                month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                         "\n",
+                         "Forecast month: ", fmonth, "\n",
+                         "Pearson correlation ; alpha = ", alpha)
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
+        fmonth <- sprintf("%02d", tstep - 1 + recipe$Analysis$Time$ftime_min)
+        toptitle <- paste(recipe$Analysis$Datasets$Reference$name, "\n",
+                          "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                          "Correlation / Start dates",
+                          recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                           nino_name, "_correlation_",
+                           recipe$Analysis$Variable$name, "_",
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("Start date: month ",
+                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)     
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      
+      PlotRobinson(map, lon = lons, lat = lats, 
+                   target_proj = projection_code, #"ESRI:54030",
+                   lat_dim = 'latitude', lon_dim = 'longitude',
+                   legend = 's2dv', style = 'point',
+                   toptitle = toptitle, bar_extra_margin = c(4,0,4,0), 
+                   caption = caption, mask = sig,
+                   fileout = plotfile, width = 8, height = 6,
+                   brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+ 
+      ## Ensemble-mean
+      map <- Subset(correl_hcst$r, along = 'time', ind = tstep)
+      sig <- Subset(correl_hcst$sig, along = 'time', ind = tstep)
+      if (gsub(".", "", recipe$Analysis$Datasets$System$name) == "") {
+        system <- recipe$Analysis$Datasets$System$name
+      } else {
+        system <-gsub(".", "", recipe$Analysis$Datasets$System$name)
+      }
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+       toptitle <- paste(recipe$Analysis$Datasets$System$name, "\n",
+                          "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                          "Correlation /",
+                          month.abb[as.numeric(fmonth)],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                           nino_name, "_correlation_",
+                           recipe$Analysis$Variable$name, "_ensmean_",
+                           system,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("Ensemble mean\n",
+                           "Nominal start date: 1st of ",
+                month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                         "Pearson correlation ; alpha = ", alpha)
+      } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+        toptitle <- paste(recipe$Analysis$Datasets$System$name, "\n",
+                          "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                          "Correlation / Start dates",
+                          recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                           nino_name, "_correlation_",
+                           recipe$Analysis$Variable$name, "_ensmean_",
+                           system,
+                           recipe$Analysis$Datasets$Reference$name,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("Correlation ensemble mean\n",
+                           "Start date month: ",
+                           month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                           "\n",
+                           "Forecast month: ", fmonth, "\n",
+                           "Pearson correlation; alpha = ", alpha)
+      } else {
+        toptitle <- NULL
+        warning("The plot title is not defined for this time horizon. ",
+                "The plots will not have a title.")
+      }
+      PlotRobinson(map, lon = lons, lat = lats,
+                   target_proj = projection_code, #"ESRI:54030",
+                   lat_dim = 'latitude', lon_dim = 'longitude',
+                   legend = 's2dv', style = 'point',
+                   toptitle = toptitle, bar_extra_margin = c(4,0,4,0),
+                   caption = caption, mask = sig,
+                   fileout = plotfile, width = 8, height = 6,
+                   brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+
+     # Full hcst corr
+      map <- Subset(correl_hcst_full$r, along = 'time', ind = tstep)
+      sig <- Subset(correl_hcst_full$sig, along = 'time', ind = tstep)
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+        toptitle <- paste(recipe$Analysis$Datasets$System$name, "\n",
+                          "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                          " Correlation /",
+                          month.abb[as.numeric(fmonth)],
+                          "/", recipe$Analysis$Time$hcst_start, "-",
+                          recipe$Analysis$Time$hcst_end)
+        plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                           nino_name, "_correlation_",
+                           recipe$Analysis$Variable$name, "_member_",
+                           system,
+                           "_s", recipe$Analysis$Time$sdate,
+                           "_ftime", fmonth, ".png")
+        caption <-  paste0("Individual members\n",
+                           "Nominal start date: 1st of ",
+                month.name[as.numeric(substr(recipe$Analysis$Time$sdate, 1,2))],
+                         "\n",
+                         "Forecast month: ", fmonth, "\n",
+                         "Pearson correlation ; alpha = ", alpha)
+    } else if (tolower(recipe$Analysis$Horizon) == "decadal"){
+      toptitle <- paste(recipe$Analysis$Datasets$System$name, "\n",
+                        "Ni\u00F1o", region_name, "SST Index -",var_name, "\n",
+                        "Correlation / Start dates",
+                        recipe$Analysis$Time$hcst_start, "-",
+                        recipe$Analysis$Time$hcst_end)
+      plotfile <- paste0(recipe$Run$output_dir, "/plots/Indices/",
+                         nino_name, "_correlation_",
+                         recipe$Analysis$Variable$name, "_ensmean_",
+                         system,
+                         recipe$Analysis$Datasets$Reference$name,
+                         "_ftime", fmonth, ".png")
+      caption <-  paste0("Correlation ensemble mean\n",
+                         "Start date month: ",
+                         month.name[get_archive(recipe)$System[[recipe$Analysis$Datasets$System$name]]$initial_month],
+                         "\n",
+                         "Forecast month: ", fmonth, "\n",
+                         "Pearson correlation; alpha = ", alpha)
+    } else {
+      toptitle <- NULL
+      warning("The plot title is not defined for this time horizon. ",
+              "The plots will not have a title.")
+    }
+    PlotRobinson(map, lon = lons, lat = lats,
+                 target_proj = projection_code, #"ESRI:54030",
+                 lat_dim = 'latitude', lon_dim = 'longitude',
+                 legend = 's2dv', style = 'point',
+                 toptitle = toptitle, bar_extra_margin = c(4,0,4,0),
+                 caption = caption, mask = sig,
+                 fileout = plotfile, width = 8, height = 6,
+                 brks = seq(-1, 1, 0.2), cols = brewer.pal(10, 'PuOr'))
+
+ 
+    }
+  }  
+  return(nino)
+} 
diff --git a/modules/Indices/R/correlation_eno.R b/modules/Indices/R/correlation_eno.R
new file mode 100644
index 0000000000000000000000000000000000000000..62a3d92b9b275decd4589474c5bf67a3a5b8ccc7
--- /dev/null
+++ b/modules/Indices/R/correlation_eno.R
@@ -0,0 +1,45 @@
+.correlation_eno <- function(exp, obs, time_dim, method, alpha, test.type, pval){
+  
+  cor <- NULL
+  cor$r = cor(x = exp, y = obs, method = method) # Correlation coefficient
+  
+  n_eff = s2dv::Eno(data = obs, time_dim = time_dim, na.action = na.pass, ncores = 1)
+  
+  if (test.type == 'one-sided'){
+    
+    t_alpha_n2 = qt(p=alpha, df = n_eff-2, lower.tail = FALSE)
+    t = cor$r * sqrt(n_eff-2) / sqrt(1-cor$r^2)
+    
+    if (anyNA(c(t,t_alpha_n2)) == FALSE & t >= t_alpha_n2 & cor$r > 0){
+      cor$sign = TRUE
+    } else {
+      cor$sign = FALSE
+    }
+    
+    # cor$n_eff <- n_eff
+    
+    if (isTRUE(pval)){
+      cor$pval <- pt(q = t, df = n_eff-2, lower.tail = FALSE)
+    }
+    
+  } else if (test.type == 'two-sided'){
+    
+    t_alpha2_n2 = qt(p=alpha/2, df = n_eff-2, lower.tail = FALSE)
+    t = abs(cor$r) * sqrt(n_eff-2) / sqrt(1-cor$r^2)
+    
+    if (anyNA(c(t,t_alpha2_n2)) == FALSE & t >= t_alpha2_n2){
+      cor$sign = TRUE
+    } else {
+      cor$sign = FALSE
+    }
+    
+    cor$n_eff <- n_eff
+    
+    if (isTRUE(pval)){
+      cor$pval <- 2 * pt(q = t, df = n_eff-2, lower.tail = FALSE)
+    }
+    
+  } else {stop('test.type not supported')}
+  
+  return(cor)
+}
diff --git a/modules/Indices/R/drop_indices_dims.R b/modules/Indices/R/drop_indices_dims.R
new file mode 100644
index 0000000000000000000000000000000000000000..aab97a170b7ebc96221bda9d3f80d3b3d5682f79
--- /dev/null
+++ b/modules/Indices/R/drop_indices_dims.R
@@ -0,0 +1,87 @@
+# version victoria https://earth.bsc.es/gitlab/es/auto-s2s/-/blob/dev-Loading-multivar/modules/Skill/Skill.R
+# metric_array is in this case the index
+.drop_indices_dims <- function(metric_array) {
+  # Define dimensions that are not essential for saving
+  droppable_dims <- c("var", "dat", "sday", "sweek", "ensemble", "nobs",
+                      "nexp", "exp_memb", "obs_memb", "bin")
+  # Select non-essential dimensions of length 1
+  dims_to_drop <- intersect(names(which(dim(metric_array) == 1)),
+                            droppable_dims)
+  drop_indices <- grep(paste(dims_to_drop, collapse = "|"),
+                       names(dim(metric_array)))
+  # Drop selected dimensions
+  metric_array <- abind::adrop(metric_array, drop = drop_indices)
+  # If array has memb dim (Corr case), change name to 'ensemble'
+  if ("exp_memb" %in% names(dim(metric_array))) {
+    names(dim(metric_array))[which(names(dim(metric_array)) ==
+                                   "exp_memb")] <- "ensemble"
+  }
+  return(metric_array)
+}
+
+
+
+
+
+## TODO: Replace with ClimProjDiags::Subset and add var and dat dimensions
+#.drop_dims <- function(metric_array) {
+#  # Drop all singleton dimensions
+#  dims <- dim(metric_array)
+#  metric_array <- drop(metric_array)
+#  if ("region" %in% names(dims)) {
+#    if (!is.array(metric_array)) {
+#      dim(metric_array) <- length(metric_array)
+#      if (!all(dims > 1)) {
+#        if ("syear" %in% names(dims)) {
+#          names(dim(metric_array)) <- "syear"
+#        } else {
+#          names(dim(metric_array)) <- "time"
+#        }
+#      } else {
+#        names(dim(metric_array)) <- names(dims[which(dims > 1)])
+#      }
+#    }
+#    if (!("time" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c("time" = 1, dim(metric_array))
+#    }
+#    # If latitude was singleton, add it back
+#    if (!("region" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "region" = 1)
+#    }
+#    # If array has memb dim (Corr case), change name to 'ensemble'
+#    if ("exp_memb" %in% names(dim(metric_array))) {
+#      names(dim(metric_array))[which(names(dim(metric_array)) ==
+#                                   "exp_memb")] <- "ensemble"
+#    }
+#  } else if (all(c("latitude", "longiguted") %in% names(dims))) {
+#    # If the array becomes a vector as a result, restore dimensions.
+#    ## This applies to the case of 'corr' with one lon and one lat
+#    if (!is.array(metric_array)) {
+#      dim(metric_array) <- length(metric_array)
+#      names(dim(metric_array)) <- names(dims[which(dims > 1)])
+#    }
+#    # If the array becomes a vector as a result, restore dimensions.
+#    ## This applies to the case of 'corr' with one lon and one lat
+#    # If time happened to be a singleton dimension, add it back in the array
+#   if (!("time" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c("time" = 1, dim(metric_array))
+#    }
+#    # If latitude was singleton, add it back
+#    if (!("latitude" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "latitude" = 1)
+#    }
+#    # If longitude was singleton, add it back
+#    if (!("longitude" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "longitude" = 1)
+#    }
+#    # If array has memb dim (Corr case), change name to 'ensemble'
+#    if ("exp_memb" %in% names(dim(metric_array))) {
+#      names(dim(metric_array))[which(names(dim(metric_array)) ==
+#                                   "exp_memb")] <- "ensemble"
+#    } else {
+#      stop("what dimensions")
+#    }
+#  }
+#  return(metric_array)
+#}
+
diff --git a/modules/Indices/R/plot_deterministic_forecast.R b/modules/Indices/R/plot_deterministic_forecast.R
new file mode 100644
index 0000000000000000000000000000000000000000..da107e4732ec027f60729b7520a8d4264787a1a0
--- /dev/null
+++ b/modules/Indices/R/plot_deterministic_forecast.R
@@ -0,0 +1,176 @@
+
+plot_deterministic_forecast <- function(obs, fcst, title = NULL,
+                                        xlabs = NULL, ylab = '', ylims = NULL,
+                                        xlabs_pos = NULL, ylabs_pos = NULL,
+                                        style = 'shading', col_obs = 'black',
+                                        col_hcst = 'blue', col_fcst = 'red', 
+                                        member_dim = 'member',
+                                        time_dim = 'year',
+                                        logo = NULL, caption = NULL,
+                                        caption_line = NULL,
+                                        legend_text = c('Observations', 'Hindcast',
+                                                        'Forecast'),
+                                        width = 1010, height = 510, res = NA,
+                                        fileout = NULL) {
+  n_obs <- as.numeric(dim(obs)[time_dim])
+  n_fcst <- as.numeric(dim(fcst)[time_dim])
+  
+  if (is.null(ylims)) {
+    ylims <- c(-max(abs(fcst) + 0.1, abs(obs) + 0.1, na.rm = TRUE),
+               max(abs(fcst) + 0.1, abs(obs) + 0.1, na.rm = TRUE))
+  }
+  if (is.null(xlabs)) {
+    xlabs <- 1:n_fcst
+  }
+  if (is.null(xlabs_pos)) {
+    xlabs_pos <- -2 + seq(from = 1, to = n_fcst, by = 5)
+  }
+  if (is.null(ylabs_pos)) {
+    ylabs_pos <- ylims[1] - 0.045 + par('usr')[3]
+  }
+  ## Opening the output file
+  if (!is.null(fileout)) {
+    png(filename = fileout, width = width, height = height, res = res)
+  }
+  if (identical(style, 'boxplot')) {
+    if (!is.null(legend_text)) {
+      right_mar <- max(c(nchar(legend_text), 8)) + 1 
+    } else {
+      right_mar <- 8
+    }
+    par(mar = c(7, 4.1, 4.1, right_mar))
+  }
+  ## Empty figure with the labs
+  plot(x = 1:n_fcst, y = rep(0,n_fcst), col = col_obs, type = 'l', lwd = 1,
+       xlab = '', ylab = ylab, main = title, xlim = c(1,n_fcst), ylim = ylims,
+       cex.lab = 1.2, cex.main = 1.2, cex.axis = 1.2, frame.plot = TRUE,
+       xaxt = 'n')
+  axis(1, at = 1:n_fcst, labels = xlabs, cex = 1.5)
+  #text(cex = 1.2, x = xlabs_pos, y = ylabs_pos, labels = xlabs[seq(from = 1, to = n_fcst, by = 5)], srt = 45, pos = 1, xpd = TRUE)
+  
+  ## Auxiliary lines
+  abline(h = 0, lw = 1)
+  abline(v = seq(from = 1, to = n_fcst, by = 5), lty = 2, lwd = 1, col = 'grey')
+  #abline(h = seq(from = floor(ylims[1]), to = ceiling(ylims[2]), by = 1), lty = 2, col = 'grey')
+  ## Ensemble spread
+  if (identical(style, 'boxplot')) {
+    ## Hindcast
+    for (y in 1:n_obs) {
+      boxplot(x = ClimProjDiags::Subset(x = fcst,
+                                        along = time_dim,
+                                        indices = y,
+                                        drop = 'selected'),
+              at = y, notch = FALSE, range = 0, add = TRUE,
+              col = adjustcolor(col_hcst, alpha.f = 0.60),
+              boxwex = 1, axes = F)
+    }
+    ## Forecast
+    if (n_fcst > n_obs) {
+      for (y in (n_obs+1):n_fcst) {
+        boxplot(x = ClimProjDiags::Subset(x = fcst,
+                                          along = time_dim,
+                                          indices = y,
+                                          drop = 'selected'),
+                at = y, notch = FALSE, range = 0, add = TRUE,
+                col = adjustcolor(col_fcst, alpha.f = 0.60),
+                boxwex = 1, axes = F)
+      }
+    }
+  } else if (identical(style, 'shading')) {
+    fcst_quantiles <- multiApply::Apply(data = fcst, target_dims = member_dim,
+                                        fun = function(x) {quantile(x = as.vector(x), type = 8, na.rm = FALSE)},
+                                        ncores = 1)$output1
+    ## Hindcast
+    polygon(x = c(1:n_obs, rev(1:n_obs)),
+            y = c(fcst_quantiles[4,1:n_obs], rev(fcst_quantiles[5,1:n_obs])),
+            col = adjustcolor(col_hcst, alpha.f = 0.10),
+            border = NA)
+    polygon(x = c(1:n_obs, rev(1:n_obs)),
+            y = c(fcst_quantiles[2,1:n_obs], rev(fcst_quantiles[4,1:n_obs])),
+            col = adjustcolor(col_hcst, alpha.f = 0.30),
+            border = NA)
+    polygon(x = c(1:n_obs, rev(1:n_obs)),
+            y = c(fcst_quantiles[1,1:n_obs],
+                  rev(fcst_quantiles[2,1:n_obs])),
+            col = adjustcolor(col_hcst, alpha.f = 0.10),
+            border = NA)
+    ## Forecast
+    if (n_fcst > n_obs){
+      polygon(x = c((n_obs):n_fcst, rev((n_obs):n_fcst)),
+              y = c(fcst_quantiles[4,(n_obs):n_fcst],
+                    rev(fcst_quantiles[5,(n_obs):n_fcst])),
+              col = adjustcolor(col_fcst, alpha.f = 0.10),
+              border = NA)
+      polygon(x = c((n_obs):n_fcst, rev((n_obs):n_fcst)),
+              y = c(fcst_quantiles[2,(n_obs):n_fcst],
+                    rev(fcst_quantiles[4,(n_obs):n_fcst])),
+              col = adjustcolor(col_fcst, alpha.f = 0.30),
+              border = NA)
+      polygon(x = c((n_obs):n_fcst, rev((n_obs):n_fcst)),
+              y = c(fcst_quantiles[1,(n_obs):n_fcst],
+                    rev(fcst_quantiles[2,(n_obs):n_fcst])),
+              col = adjustcolor(col_fcst, alpha.f = 0.10), 
+              border = NA)
+    }
+    ## NOTEV: What is 'fsct_stype'?
+  } else {stop('fcst_stype must be either boxplot or shading')}
+  
+  ## Observations
+  lines(x = 1:n_obs, y = obs, col = col_obs, lty = 'solid', lwd = 5)
+  
+  ## Ensemble mean
+  ## Hindcast
+  lines(x = 1:n_obs,
+        y = multiApply::Apply(data = fcst,
+                              target_dims = member_dim, 
+                              fun = mean,
+                              na.rm = FALSE,
+                              ncores = 1)$output1[1:n_obs],
+        col = col_hcst, type = 'l', lwd = 5)
+  ## Forecast
+  if (n_fcst > n_obs) {
+    lines(x = (n_obs):n_fcst,
+          y = multiApply::Apply(data = fcst,
+                                target_dims = member_dim,
+                                fun = mean,
+                                na.rm = FALSE,
+                                ncores = 1)$output1[(n_obs):n_fcst],
+          col = col_fcst, type = 'l', lwd = 5)
+  }
+  
+  ## Legend with the skill scores
+  if (!is.null(legend)) {
+    par(xpd = TRUE)
+    legend(par('usr')[2], par('usr')[4],#x = 1, y = ylims[2], 
+           text.col = c(col_obs, col_hcst), cex = 1.2, 
+           legend = legend_text[1:2], lty = 1, lwd = 3, 
+           col = c(col_obs, col_hcst, bg="transparent"), bty = 'n')
+  }
+  if (!is.null(caption)) {
+    if (is.null(caption_line)) {
+      caption_line <- 4.5
+     }
+     mtext(caption, side = 1, line = caption_line, cex = 1.2, adj = 0)
+  }
+  if (identical(style, 'boxplot')){
+    if (!is.null(logo)) {
+      rasterImage(logo, ybottom = par('usr')[3] - 2, ytop = par('usr')[3] - 1.1, #ylims[2] + 2, 
+                  xleft = par('usr')[2] - 1.5 , 
+                  xright = par('usr')[2] + dim(logo)[2]/dim(logo)[1] * width/height, xpd = T)
+    }
+    data_samp <- 1:100
+    data_samp <- (data_samp - mean(data_samp))/ sqrt(var(data_samp))
+    par(xpd = T)
+    box_params <- boxplot(data_samp,  add = T, at = par('usr')[2] + n_fcst * 0.1)
+    text(x = par('usr')[2] + n_fcst * 0.1 + .5, y = box_params$stats[1,], "Minimum", adj = c(0,0))
+    text(x = par('usr')[2] + n_fcst * 0.1 + .5, y = box_params$stats[2,], "Q1", adj = c(0,0))
+    text(x = par('usr')[2] + n_fcst * 0.1 + .5, y = box_params$stats[3,], "Median", adj = c(0,0))
+    text(x = par('usr')[2] + n_fcst * 0.1 + .5, y = box_params$stats[4,], "Q3", adj = c(0,0))
+    text(x = par('usr')[2] + n_fcst * 0.1 + .5, y = box_params$stats[5,], "Maximum", adj = c(0,0))
+  }
+  
+  ## Closing the output file
+  if (!is.null(fileout)){dev.off()}
+  
+}
+
diff --git a/modules/Saving/R/drop_dims.R b/modules/Saving/R/drop_dims.R
new file mode 100644
index 0000000000000000000000000000000000000000..7361faa8ea1056b414b5e2ed67740e4ceeb38bae
--- /dev/null
+++ b/modules/Saving/R/drop_dims.R
@@ -0,0 +1,86 @@
+# version victoria https://earth.bsc.es/gitlab/es/auto-s2s/-/blob/dev-Loading-multivar/modules/Skill/Skill.R
+.drop_dims <- function(metric_array) {
+  # Define dimensions that are not essential for saving
+  droppable_dims <- c("dat", "sday", "sweek", "ensemble", "nobs",
+                      "nexp", "exp_memb", "obs_memb", "bin")
+  # Select non-essential dimensions of length 1
+  dims_to_drop <- intersect(names(which(dim(metric_array) == 1)),
+                            droppable_dims)
+  drop_indices <- grep(paste(dims_to_drop, collapse = "|"),
+                       names(dim(metric_array)))
+  # Drop selected dimensions
+  metric_array <- abind::adrop(metric_array, drop = drop_indices)
+  # If array has memb dim (Corr case), change name to 'ensemble'
+  if ("exp_memb" %in% names(dim(metric_array))) {
+    names(dim(metric_array))[which(names(dim(metric_array)) ==
+                                   "exp_memb")] <- "ensemble"
+  }
+  return(metric_array)
+}
+
+
+
+
+
+## TODO: Replace with ClimProjDiags::Subset and add var and dat dimensions
+#.drop_dims <- function(metric_array) {
+#  # Drop all singleton dimensions
+#  dims <- dim(metric_array)
+#  metric_array <- drop(metric_array)
+#  if ("region" %in% names(dims)) {
+#    if (!is.array(metric_array)) {
+#      dim(metric_array) <- length(metric_array)
+#      if (!all(dims > 1)) {
+#        if ("syear" %in% names(dims)) {
+#          names(dim(metric_array)) <- "syear"
+#        } else {
+#          names(dim(metric_array)) <- "time"
+#        }
+#      } else {
+#        names(dim(metric_array)) <- names(dims[which(dims > 1)])
+#      }
+#    }
+#    if (!("time" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c("time" = 1, dim(metric_array))
+#    }
+#    # If latitude was singleton, add it back
+#    if (!("region" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "region" = 1)
+#    }
+#    # If array has memb dim (Corr case), change name to 'ensemble'
+#    if ("exp_memb" %in% names(dim(metric_array))) {
+#      names(dim(metric_array))[which(names(dim(metric_array)) ==
+#                                   "exp_memb")] <- "ensemble"
+#    }
+#  } else if (all(c("latitude", "longiguted") %in% names(dims))) {
+#    # If the array becomes a vector as a result, restore dimensions.
+#    ## This applies to the case of 'corr' with one lon and one lat
+#    if (!is.array(metric_array)) {
+#      dim(metric_array) <- length(metric_array)
+#      names(dim(metric_array)) <- names(dims[which(dims > 1)])
+#    }
+#    # If the array becomes a vector as a result, restore dimensions.
+#    ## This applies to the case of 'corr' with one lon and one lat
+#    # If time happened to be a singleton dimension, add it back in the array
+#   if (!("time" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c("time" = 1, dim(metric_array))
+#    }
+#    # If latitude was singleton, add it back
+#    if (!("latitude" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "latitude" = 1)
+#    }
+#    # If longitude was singleton, add it back
+#    if (!("longitude" %in% names(dim(metric_array)))) {
+#      dim(metric_array) <- c(dim(metric_array), "longitude" = 1)
+#    }
+#    # If array has memb dim (Corr case), change name to 'ensemble'
+#    if ("exp_memb" %in% names(dim(metric_array))) {
+#      names(dim(metric_array))[which(names(dim(metric_array)) ==
+#                                   "exp_memb")] <- "ensemble"
+#    } else {
+#      stop("what dimensions")
+#    }
+#  }
+#  return(metric_array)
+#}
+
diff --git a/modules/Saving/R/get_filename.R b/modules/Saving/R/get_filename.R
index 1c92565179ab74ee0968323258b3efa18e396f77..f18a6fb48e1c893467fcae1160aba392833e9393 100644
--- a/modules/Saving/R/get_filename.R
+++ b/modules/Saving/R/get_filename.R
@@ -17,7 +17,8 @@ get_filename <- function(dir, recipe, var, date, agg, file.type) {
 
   switch(tolower(agg),
 	 "country" = {gg <- "-country"},
-         "global" = {gg <- ""})
+         "global" = {gg <- ""},
+         "region" = {gg <- "-region"})
 
   system <- gsub('.','', recipe$Analysis$Datasets$System$name, fixed = T)
   reference <- gsub('.','', recipe$Analysis$Datasets$Reference$name, fixed = T)
diff --git a/modules/Saving/R/get_global_attributes.R b/modules/Saving/R/get_global_attributes.R
new file mode 100644
index 0000000000000000000000000000000000000000..25c47279c19911fbf47dd9de8d9be635913958ba
--- /dev/null
+++ b/modules/Saving/R/get_global_attributes.R
@@ -0,0 +1,20 @@
+get_global_attributes <- function(recipe, archive) {
+  # Generates metadata of interest to add to the global attributes of the
+  # netCDF files.
+  parameters <- recipe$Analysis
+  hcst_period <- paste0(parameters$Time$hcst_start, " to ",
+                        parameters$Time$hcst_end)
+  current_time <- paste0(as.character(Sys.time()), " ", Sys.timezone())
+  system_name <- parameters$Datasets$System$name
+  reference_name <- parameters$Datasets$Reference$name
+
+  attrs <- list(reference_period = hcst_period,
+                institution_system = archive$System[[system_name]]$institution,
+                institution_reference = archive$Reference[[reference_name]]$institution,
+                system = system_name,
+                reference = reference_name,
+                calibration_method = parameters$Workflow$Calibration$method,
+                computed_on = current_time)
+
+  return(attrs)
+}
diff --git a/modules/Saving/R/get_latlon.R b/modules/Saving/R/get_latlon.R
new file mode 100644
index 0000000000000000000000000000000000000000..25d536380bd67300fb0ca57a839f467470b6f9ea
--- /dev/null
+++ b/modules/Saving/R/get_latlon.R
@@ -0,0 +1,17 @@
+get_latlon <- function(latitude, longitude) {
+  # Adds dimensions and metadata to lat and lon
+  # latitude: array containing the latitude values
+  # longitude: array containing the longitude values
+  dim(longitude) <- length(longitude)
+  metadata <- list(longitude = list(units = 'degrees_east'))
+  attr(longitude, 'variables') <- metadata
+  names(dim(longitude)) <- 'longitude'
+
+  dim(latitude) <- length(latitude)
+  metadata <- list(latitude = list(units = 'degrees_north'))
+  attr(latitude, 'variables') <- metadata
+  names(dim(latitude)) <- 'latitude'
+
+  return(list(lat=latitude, lon=longitude))
+
+}
diff --git a/modules/Saving/R/get_time.R b/modules/Saving/R/get_time.R
new file mode 100644
index 0000000000000000000000000000000000000000..5426c177f1bd31e14b27ca9e831eb880da183f04
--- /dev/null
+++ b/modules/Saving/R/get_time.R
@@ -0,0 +1,30 @@
+get_times <- function(store.freq, fcst.horizon, leadtimes, sdate, calendar) {
+  # Generates time dimensions and the corresponding metadata.
+  ## TODO: Subseasonal
+
+  switch(fcst.horizon,
+         "seasonal" = {time <- leadtimes; ref <- 'hours since ';
+                       stdname <- paste(strtoi(leadtimes), collapse=", ")},
+         "subseasonal" = {len <- 4; ref <- 'hours since ';
+                          stdname <- ''},
+         "decadal" = {time <- leadtimes; ref <- 'hours since ';
+                      stdname <- paste(strtoi(leadtimes), collapse=", ")})
+
+  dim(time) <- length(time)
+  sdate <- as.Date(sdate, format = '%Y%m%d') # reformatting
+  metadata <- list(time = list(units = paste0(ref, sdate, 'T00:00:00'),
+                   calendar = calendar))
+  attr(time, 'variables') <- metadata
+  names(dim(time)) <- 'time'
+
+  sdate <- 1:length(sdate)
+  dim(sdate) <- length(sdate)
+  metadata <- list(sdate = list(standard_name = paste(strtoi(sdate),
+                                                      collapse=", "),
+                                units = paste0('Init date')))
+  attr(sdate, 'variables') <- metadata
+  names(dim(sdate)) <- 'sdate'
+
+  return(list(time=time))
+}
+
diff --git a/modules/Saving/R/save_metrics.R b/modules/Saving/R/save_metrics.R
index 609537c0226b39c1ac4058eec7029fb9098720f6..a232d89069f6e7298b88c124e0aea869ca4c5367 100644
--- a/modules/Saving/R/save_metrics.R
+++ b/modules/Saving/R/save_metrics.R
@@ -1,17 +1,16 @@
 save_metrics <- function(recipe,
                          skill,
+                         dictionary = NULL,
                          data_cube,
                          agg = "global",
                          outdir = NULL) {
   # This function adds metadata to the skill metrics in 'skill'
   # and exports them to a netCDF file inside 'outdir'.
-  
   # Define grid dimensions and names
   lalo <- c('longitude', 'latitude')
   archive <- get_archive(recipe)
   dictionary <- read_yaml("conf/variable-dictionary.yml")
 
-  # Add global and variable attributes
   global_attributes <- .get_global_attributes(recipe, archive)
   ## TODO: Sort out the logic once default behavior is decided
   if ((!is.null(recipe$Analysis$Workflow$Anomalies$compute)) && 
@@ -22,7 +21,6 @@ save_metrics <- function(recipe,
     global_attributes <- c(list(from_anomalies = "No"),
                            global_attributes)
   }
-
   # Time indices and metadata
   fcst.horizon <- tolower(recipe$Analysis$Horizon)
   store.freq <- recipe$Analysis$Variables$freq
@@ -67,7 +65,6 @@ save_metrics <- function(recipe,
 
   times <- .get_times(store.freq, fcst.horizon, leadtimes, fcst.sdate, calendar)
   time <- times$time
-
   # Loop over variable dimension
   for (var in 1:data_cube$dims[['var']]) {
     # Subset skill arrays
@@ -98,6 +95,9 @@ save_metrics <- function(recipe,
       if (tolower(agg) == "country") {
         sdname <- paste0(metric, " region-aggregated metric")
         dims <- c('Country', 'time')
+      } else if (tolower(agg) == "region") {
+        sdname <- paste0(metric, " region-aggregated metric")
+        dims <- c('region', 'time')
       } else {
         sdname <- paste0(metric)
         dims <- c(lalo, 'time')
@@ -113,6 +113,14 @@ save_metrics <- function(recipe,
     if (tolower(agg) == "country") { 
       country <- get_countries(grid)
       ArrayToNc(append(country, time, subset_skill), outfile)
+    } else if (tolower(agg) == "region") {
+      region <- array(1:dim(skill[[1]])['region'], c(dim(skill[[1]])['region']))
+      # TODO: check metadata when more than 1 region is store in the data array
+      metadata <- list(region = list(long_name = data_cube$attrs$Variable$metadata$region$name))
+      attr(region, 'variables') <-  metadata
+      vars <- list(region, time)
+      vars <- c(vars, subset_skill)
+      ArrayToNc(vars, outfile)
     } else {
       latitude <- data_cube$coords$lat[1:length(data_cube$coords$lat)]
       longitude <- data_cube$coords$lon[1:length(data_cube$coords$lon)]
@@ -125,3 +133,4 @@ save_metrics <- function(recipe,
     info(recipe$Run$logger, "##### SKILL METRICS SAVED TO NETCDF FILE #####")
   }
 }
+
diff --git a/modules/Saving/Saving.R b/modules/Saving/Saving.R
index c52991d38337c59a40a422b4b5fd526a5365e761..2d067ba394fdc4a3bc0ea47f6033fb0326c8c6a2 100644
--- a/modules/Saving/Saving.R
+++ b/modules/Saving/Saving.R
@@ -10,10 +10,15 @@ source("modules/Saving/R/save_metrics.R")
 source("modules/Saving/R/save_corr.R")
 source("modules/Saving/R/save_probabilities.R")
 source("modules/Saving/R/save_percentiles.R")
-
+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) {
+		      skill_metrics = NULL,
+		      probabilities = NULL,
+                      agg = 'global',
+		      archive = NULL) {
   # Wrapper for the saving functions.
   # recipe: The auto-s2s recipe
   # archive: The auto-s2s archive
@@ -68,7 +73,7 @@ save_data <- function(recipe, data,
   if (!is.null(skill_metrics)) {
       save_metrics(recipe = recipe,
                    skill = skill_metrics,
-                   data_cube = data$hcst,
+                   data_cube = data$hcst, agg = agg,
                    outdir = outdir[var])
   }
   if (!is.null(corr_metrics)) {
diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R
index f22d6b9b8bc75bd6639471705b90a5deb03a86f0..ecb89224149ebc3783a229b3b2532109e9712d70 100644
--- a/modules/Skill/Skill.R
+++ b/modules/Skill/Skill.R
@@ -10,6 +10,7 @@
 source("modules/Skill/R/compute_quants.R")
 source("modules/Skill/R/compute_probs.R")
 source("modules/Skill/R/s2s.metrics.R")
+source("modules/Saving/R/drop_dims.R")
 ## TODO: Remove with the next release of s2dv
 source("modules/Skill/R/tmp/RPS.R")
 source("modules/Skill/R/tmp/RPSS.R")
@@ -54,9 +55,9 @@ source("modules/Skill/R/CRPS_clim.R")
 #                    " 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) {
+# 				  clim_data$hcst = NULL,
+# 				  clim_obs = NULL) {
+compute_skill_metrics <- function(recipe, data, agg = 'global') {
 
   # data$hcst: s2dv_cube containing the hindcast
   # obs: s2dv_cube containing the observations
@@ -316,6 +317,7 @@ compute_skill_metrics <- function(recipe, data) {
   }
   info(recipe$Run$logger, "##### SKILL METRIC COMPUTATION COMPLETE #####")
   # Save outputs
+#NURIA: I think change the output_dir is a problem for future savings
   recipe$Run$output_dir <- paste0(recipe$Run$output_dir,
                                   "/outputs/Skill/")
   # Separate 'corr' from the rest of the metrics because of extra 'ensemble' dim
@@ -323,7 +325,7 @@ compute_skill_metrics <- function(recipe, data) {
     corr_metric_names <- grep("^corr", names(skill_metrics))
     if (length(corr_metric_names) == 0) {
       save_metrics(recipe = recipe, skill = skill_metrics,
-                   data_cube = data$hcst)
+                   data_cube = data$hcst, agg = agg)
     } else {
       # Save corr
       if (length(skill_metrics[corr_metric_names]) > 0) {
@@ -333,7 +335,7 @@ compute_skill_metrics <- function(recipe, data) {
       # Save other skill metrics
       if (length(skill_metrics[-corr_metric_names]) > 0) {
         save_metrics(recipe = recipe, skill = skill_metrics[-corr_metric_names],
-                     data_cube = data$hcst)
+                     data_cube = data$hcst, agg = agg)
       }
     }
   }
@@ -419,7 +421,6 @@ compute_probabilities <- function(recipe, data) {
         }
       }
     }
-
     # Rearrange dimensions and return probabilities
     named_probs <- lapply(named_probs, function(x) {.drop_dims(x)})
     named_quantiles <- lapply(named_quantiles, function(x) {.drop_dims(x)})
@@ -464,21 +465,3 @@ compute_probabilities <- function(recipe, data) {
   }
 }
 
-.drop_dims <- function(metric_array) {
-  # Define dimensions that are not essential for saving
-  droppable_dims <- c("dat", "sday", "sweek", "ensemble", "nobs",
-                      "nexp", "exp_memb", "obs_memb", "bin")
-  # Select non-essential dimensions of length 1
-  dims_to_drop <- intersect(names(which(dim(metric_array) == 1)),
-                            droppable_dims)
-  drop_indices <- grep(paste(dims_to_drop, collapse = "|"),
-                       names(dim(metric_array)))
-  # Drop selected dimensions
-  metric_array <- abind::adrop(metric_array, drop = drop_indices)
-  # If array has memb dim (Corr case), change name to 'ensemble'
-  if ("exp_memb" %in% names(dim(metric_array))) {
-    names(dim(metric_array))[which(names(dim(metric_array)) ==
-                                   "exp_memb")] <- "ensemble"
-  }
-  return(metric_array)
-}
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index 3263076d5977932169f27e5c94a8230b944acd13..07e35706d617d16b6e9b33cb2d08418597692224 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -56,7 +56,6 @@ plot_ensemble_mean <- function(recipe, fcst, outdir) {
       brks <- pretty(range(var_ens_mean, na.rm = T), n = 15, min.n = 8)
     }
     cols <- grDevices::hcl.colors(length(brks) - 1, palette, rev = rev)
-    options(bitmapType = "cairo")
 
     for (i_syear in start_date) {
       if (length(start_date) == 1) {
diff --git a/modules/Visualization/R/plot_skill_metrics.R b/modules/Visualization/R/plot_skill_metrics.R
index 834cc8cc05578df530afadf799db71383f464cb0..59f487c2c038bc4e00155da32f4d9931a10d61d3 100644
--- a/modules/Visualization/R/plot_skill_metrics.R
+++ b/modules/Visualization/R/plot_skill_metrics.R
@@ -111,7 +111,6 @@ plot_skill_metrics <- function(recipe, data_cube, skill_metrics,
           col_inf <- colorbar[1]
           col_sup <- colorbar[length(colorbar)]
         }
-        options(bitmapType = "cairo")
         # Reorder dimensions
         skill <- Reorder(skill, c("time", "longitude", "latitude"))
         # If the significance has been requested and the variable has it,
diff --git a/recipes/examples/NAO_recipe.yml b/recipes/examples/NAO_recipe.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6db51af3bdff0e01484c088e99972448f43fd0c1
--- /dev/null
+++ b/recipes/examples/NAO_recipe.yml
@@ -0,0 +1,60 @@
+Description:
+  Author: nperez
+  Info: ECMWF SEAS5 Seasonal Forecast Example recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: psl
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: ECMWF-SEAS5 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0301' ## MMDD
+   # fcst_year:  # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 2 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 2 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 20 # Mandatory, int: minimum latitude
+    latmax: 80 # Mandatory, int: maximum latitude
+    lonmin: -80 # Mandatory, int: minimum longitude
+    lonmax: 40 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: yes
+      cross_validation: no
+      save: none
+    Indices:
+      NAO: {ftime_avg: no, obsproj: TRUE, save: 'all', plot_ts: TRUE, plot_sp: yes}
+    Calibration:
+      method: raw # Mandatory, str: Calibration method. See docu.
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rps rpss crps crpss EnsSprErr
+      save: 'all'
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: none
+    Indicators:
+      index: no
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/vagudets/auto-s2s-outputs/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
+
+
diff --git a/recipes/examples/Nino_recipe.yml b/recipes/examples/Nino_recipe.yml
new file mode 100644
index 0000000000000000000000000000000000000000..257b6f532d0575a1500d0115b0806467f5c110e2
--- /dev/null
+++ b/recipes/examples/Nino_recipe.yml
@@ -0,0 +1,61 @@
+Description:
+  Author: nperez
+  Info: ECMWF SEAS5 Seasonal Forecast Example recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: tos
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: ECMWF-SEAS5 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0101' ## MMDD
+   # fcst_year:  # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2016' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 2 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 2 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: -90 # Mandatory, int: minimum latitude
+    latmax: 90 # Mandatory, int: maximum latitude
+    lonmin: 0 # Mandatory, int: minimum longitude
+    lonmax: 359.9 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "to_system"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: yes
+      cross_validation: no
+      save: none
+    Indices:
+      Nino1+2: {save: 'all', plot_ts: TRUE, plot_sp: TRUE}
+      Nino3: {save: 'all', plot_ts: TRUE, plot_sp: TRUE}
+      Nino3.4: {save: 'all', plot_ts: TRUE, plot_sp: TRUE}
+      Nino4: {save: 'all', plot_ts: TRUE, plot_sp: TRUE}
+    Calibration:
+      method: raw # Mandatory, str: Calibration method. See docu.
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rps rpss crps crpss EnsSprErr
+      save: 'all'
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: none
+    Indicators:
+      index: no
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: TRUE
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/nperez/scorecards_data/input_test/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/examples/recipe_Nino_decadal.yml b/recipes/examples/recipe_Nino_decadal.yml
new file mode 100644
index 0000000000000000000000000000000000000000..6eb9619f86958361244241d1d6e0f7a77aef31f1
--- /dev/null
+++ b/recipes/examples/recipe_Nino_decadal.yml
@@ -0,0 +1,55 @@
+Description:
+  Author: Carlos Delgado
+  Info: Test for NAO indices in decadal prediction
+Analysis:
+  Horizon: Decadal
+  Variables:
+    name: psl
+    freq: monthly_mean
+  Datasets:
+    System: 
+      name: EC-Earth3-i4
+      member: 'all'
+    Multimodel: no
+    Reference:
+      name: JRA-55
+  Time:
+    fcst_year: 
+    hcst_start: 1996
+    hcst_end: 2016
+    ftime_min: 3
+    ftime_max: 4 
+  Region:
+    latmin: -90
+    latmax: 90
+    lonmin: 0
+    lonmax:  359.9
+  Regrid:
+    method: bilinear
+    type: "r180x90"#to_system
+  Workflow:
+    Anomalies: 
+      compute: yes
+      cross_validation: no
+      save: none
+    Indices:
+      Nino3.4: {save: 'all', plot_ts: TRUE, plot_sp: yes}
+    Calibration:
+      method: raw
+      save: none
+    Skill:
+      metric: EnsCorr RPSS
+      save: "all"
+    Probabilities:
+      percentiles: [[1/3, 2/3]]
+      save: none
+    Indicators:
+      index: FALSE
+  ncores: 8 # Optional, int: number of cores, defaults to 1
+  remove_NAs: FALSE # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/nperez/scorecards_data/input_test/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/recipes/examples/recipe_model_decadal_NAO.yml b/recipes/examples/recipe_model_decadal_NAO.yml
new file mode 100644
index 0000000000000000000000000000000000000000..cd1b2bdd1b7d797dbf62f5c9235362de4c2632b6
--- /dev/null
+++ b/recipes/examples/recipe_model_decadal_NAO.yml
@@ -0,0 +1,55 @@
+Description:
+  Author: Carlos Delgado
+  Info: Test for NAO indices in decadal prediction
+Analysis:
+  Horizon: Decadal
+  Variables:
+    name: psl
+    freq: monthly_mean
+  Datasets:
+    System: 
+      name: EC-Earth3-i4
+      member: 'all'
+    Multimodel: no
+    Reference:
+      name: JRA-55
+  Time:
+    fcst_year: 
+    hcst_start: 1995
+    hcst_end: 2016
+    ftime_min: 3
+    ftime_max: 5
+  Region:
+    latmin: 20
+    latmax: 80
+    lonmin: -80
+    lonmax:  40
+  Regrid:
+    method: bilinear
+    type: "r180x90"#to_system
+  Workflow:
+    Anomalies: 
+      compute: yes
+      cross_validation: no
+      save: none
+    Indices:
+      NAO: {obsproj: TRUE, save: 'all', plot_ts: TRUE, plot_sp: yes}
+    Calibration:
+      method: raw
+      save: none
+    Skill:
+      metric: EnsCorr RPSS
+      save: "all"
+    Probabilities:
+      percentiles: [[1/3, 2/3]]
+      save: none
+    Indicators:
+      index: FALSE
+  ncores: 8 # Optional, int: number of cores, defaults to 1
+  remove_NAs: FALSE # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: /esarchive/scratch/nperez/scorecards_data/input_test/
+  code_dir: /esarchive/scratch/nperez/git/s2s-suite/
diff --git a/tests/recipes/recipe-seasonal_NAO.yml b/tests/recipes/recipe-seasonal_NAO.yml
new file mode 100644
index 0000000000000000000000000000000000000000..95d5a86e64be204c57681672a5c2859c10933987
--- /dev/null
+++ b/tests/recipes/recipe-seasonal_NAO.yml
@@ -0,0 +1,62 @@
+Description:
+  Author: nperez
+  Info: ECMWF SEAS5 Seasonal Forecast Example recipe (monthly mean, tas)
+
+Analysis:
+  Horizon: seasonal # Mandatory, str: either subseasonal, seasonal, or decadal
+  Variables:
+    name: psl
+    freq: monthly_mean
+  Datasets:
+    System:
+      name: ECMWF-SEAS5 # Mandatory, str: system5c3s system21_m1 system35c3s
+    Multimodel: no # Mandatory, bool: Either yes/true or no/false
+    Reference:
+      name: ERA5 # Mandatory, str: Reference codename. See docu.
+  Time:
+    sdate: '0301' ## MMDD
+   # fcst_year:  # Optional, int: Forecast year 'YYYY'
+    hcst_start: '1993' # Mandatory, int: Hindcast start year 'YYYY'
+    hcst_end: '2000' # Mandatory, int: Hindcast end year 'YYYY'
+    ftime_min: 2 # Mandatory, int: First leadtime time step in months 
+    ftime_max: 2 # Mandatory, int: Last leadtime time step in months
+  Region:
+    latmin: 20 # Mandatory, int: minimum latitude
+    latmax: 80 # Mandatory, int: maximum latitude
+    lonmin: -80 # Mandatory, int: minimum longitude
+    lonmax: 40 # Mandatory, int: maximum longitude
+  Regrid:
+    method: bilinear # Mandatory, str: Interpolation method. See docu.
+    type: "r180x90"
+    #type: /esarchive/scratch/nmilders/gitlab/git_clones/auto-s2s/conf/grid_description.txt #'r360x180' # Mandatory, str: to_system, to_reference, or CDO-accepted grid.
+  Workflow:
+    Anomalies:
+      compute: yes
+      cross_validation: no
+      save: none
+    Indices:
+      NAO: {ftime_avg: no, obsproj: TRUE, save: 'all', plot_ts: TRUE, plot_sp: yes}
+    Calibration:
+      method: raw # Mandatory, str: Calibration method. See docu.
+      save: none
+    Skill:
+      metric: mean_bias EnsCorr rps rpss crps crpss EnsSprErr
+      save: 'all'
+    Probabilities:
+      percentiles: [[1/3, 2/3], [1/10, 9/10]] # frac: Quantile thresholds.
+      save: none
+    Indicators:
+      index: no
+  ncores: 4 # Optional, int: number of cores, defaults to 1
+  remove_NAs: # Optional, bool: Whether NAs are removed, defaults to FALSE
+  Output_format: scorecards
+  logo: yes
+Run:
+  Loglevel: INFO
+  Terminal: yes
+  output_dir: ./tests/out-logs/
+  code_dir: /esarchive/scratch/vagudets/repos/auto-s2s/
+  #output_dir: /esarchive/scratch/nperez/scorecards_data/input_test/
+  #code_dir: /esarchive/scratch/nperez/git/s2s-suite/
+
+
diff --git a/tests/testthat/test-seasonal_NAO.R b/tests/testthat/test-seasonal_NAO.R
new file mode 100644
index 0000000000000000000000000000000000000000..7a03ecbff94e6915540eb1518b0a6af4868196cf
--- /dev/null
+++ b/tests/testthat/test-seasonal_NAO.R
@@ -0,0 +1,261 @@
+context("Seasonal NAO monthly data")
+
+source("modules/Loading/Loading.R")
+source("modules/Anomalies/Anomalies.R")
+source("modules/Skill/Skill.R")
+source("modules/Saving/Saving.R")
+source("modules/Indices/Indices.R")
+
+recipe_file <- "tests/recipes/recipe-seasonal_NAO.yml"
+recipe <- prepare_outputs(recipe_file, disable_checks = F)
+archive <- read_yaml(paste0(recipe$Run$code_dir, "conf/archive.yml"))$archive
+
+# Load datasets
+suppressWarnings({invisible(capture.output(
+data <- load_datasets(recipe)
+))})
+
+# Compute anomalies
+suppressWarnings({invisible(capture.output(
+ano_data <- compute_anomalies(recipe, data)
+))})
+
+# Compute index NAO
+suppressWarnings({invisible(capture.output(
+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')
+))})
+
+
+outdir <- get_dir(recipe = recipe, variable = data$hcst$attrs$Variable$varName)
+
+# ------- TESTS --------
+
+test_that("1. Loading", {
+
+expect_equal(
+is.list(data),
+TRUE
+)
+expect_equal(
+names(data),
+c("hcst", "fcst", "obs")
+)
+expect_equal(
+class(data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(data$fcst),
+"NULL"
+)
+expect_equal(
+class(data$obs),
+"s2dv_cube"
+)
+expect_equal(
+names(data$hcst),
+c("data", "dims", "coords", "attrs")
+)
+expect_equal(
+names(data$hcst),
+names(data$obs)
+)
+expect_equal(
+dim(data$hcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 8, time = 1, latitude = 30, longitude = 61, ensemble = 25)
+)
+expect_equal(
+dim(data$hcst$attrs$Dates),
+c(sday = 1, sweek = 1, syear = 8, time = 1)
+)
+expect_equal(
+as.vector(drop(data$hcst$data)[1:2,1,2,3]),
+c(101607.2, 101523.7),
+tolerance = 0.0001
+)
+expect_equal(
+mean(data$hcst$data),
+101450,
+tolerance = 0.0001
+)
+expect_equal(
+range(data$hcst$data),
+c(98909.3, 103299.8),
+tolerance = 0.0001
+)
+expect_equal(
+(data$hcst$attrs$Dates)[1],
+as.POSIXct("1993-04-30 18:00:00", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$attrs$Dates)[2],
+as.POSIXct("1994-04-30 18:00:00", tz = 'UTC')
+)
+expect_equal(
+(data$hcst$attrs$Dates)[5],
+as.POSIXct("1997-04-30 18:00:00", tz = 'UTC')
+)
+expect_equal(
+(data$obs$attrs$Dates)[8],
+as.POSIXct("2000-04-16", tz = 'UTC')
+)
+
+})
+
+test_that("2. Anomalies", {
+
+expect_equal(
+is.list(ano_data),
+TRUE
+)
+expect_equal(
+names(ano_data),
+c("hcst", "obs", "fcst", "hcst.full_val", "obs.full_val")
+)
+expect_equal(
+class(ano_data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(ano_data$fcst),
+"NULL"
+)
+expect_equal(
+dim(ano_data$hcst$data),
+c(dat = 1, var = 1, sday = 1, sweek = 1, syear = 8, time = 1, latitude = 30, longitude = 61, ensemble = 25)
+)
+expect_equal(
+mean(ano_data$hcst$data),
+-1.044053e-13,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(ano_data$hcst$data)[1:3, 3, 4, 1]),
+c(-165.688477, -233.336914, 3.358398),
+tolerance = 0.0001
+)
+expect_equal(
+range(ano_data$hcst$data),
+c(-2045.577, 2027.675),
+tolerance = 0.0001
+)
+
+})
+
+#=====================================
+test_that("3. Indices", {
+
+expect_equal(
+is.list(nao_data),
+TRUE
+)
+expect_equal(
+names(nao_data),
+c("hcst", "obs") 
+)
+expect_equal(
+class(nao_data$hcst),
+"s2dv_cube"
+)
+expect_equal(
+class(nao_data$fcst),
+"NULL"
+)
+expect_equal(
+dim(nao_data$hcst$data),
+c(region = 1, syear = 8, ensemble = 25, dat = 1, var = 1, sday = 1, sweek = 1, time = 1)
+)
+expect_equal(
+mean(nao_data$hcst$data),
+8.239937e-18,
+tolerance = 0.0001
+)
+expect_equal(
+as.vector(drop(nao_data$hcst$data)[1:3, 1]),
+c(0.6778843, -1.6014398, 0.3849058),
+tolerance = 0.0001
+)
+expect_equal(
+range(nao_data$hcst$data),
+c(-3.038780, 2.193189),
+tolerance = 0.0001
+)
+
+outputs <- paste0(recipe$Run$output_dir, "/plots/")
+expect_equal(
+all(list.files(outputs, recursive = T) %in%
+paste0("Indices/",
+  c("NAO_ECMWF-SEAS5_ERA5_s0301_ftime02.png",
+    "NAO_correlation_psl_ERA5_s0301_ftime02.png",
+    "NAO_correlation_psl_ensmean_ECMWF-SEAS5_s0301_ftime02.png",
+    "NAO_correlation_psl_member_ECMWF-SEAS5_s0301_ftime02.png"))),
+TRUE)
+
+expect_equal(
+length(list.files(outputs, recursive = T)),
+4
+)
+
+
+})
+
+#======================================
+test_that("3. Metrics", {
+
+expect_equal(
+is.list(skill_metrics),
+TRUE
+)
+expect_equal(
+names(skill_metrics),
+c("mean_bias", "enscorr",
+  "enscorr_significance", "rps", "rps_clim", "rpss", "rpss_significance",
+  "crps", "crps_clim", "crpss", "crpss_significance",  "enssprerr")
+)
+expect_equal(
+class(skill_metrics$rpss),
+"array"
+)
+expect_equal(
+dim(skill_metrics$rpss),
+c(region = 1, var = 1, time = 1)
+)
+expect_equal(
+dim(skill_metrics$rpss_significance),
+dim(skill_metrics$rpss)
+)
+expect_equal(
+skill_metrics$rpss[1],
+-0.06814118,
+tolerance = 0.0001
+)
+expect_equal(
+skill_metrics$rpss_significance[1],
+FALSE
+)
+
+})
+
+test_that("4. Saving", {
+outputs <- paste0(recipe$Run$output_dir, "/outputs/")
+expect_equal(
+all(list.files(outputs, recursive = T) %in%
+c(paste0("Indices/ECMWF-SEAS5/nao/", paste0("nao_", 1993:2000, "0301.nc")),
+  paste0("Indices/ERA5/nao/", paste0("nao_", 1993:2000, "0301.nc")),
+  "Skill/ECMWF-SEAS5/nao/scorecards_ECMWF-SEAS5_ERA5_nao-skill_1993-2000_s03.nc")
+), TRUE)
+
+expect_equal(
+length(list.files(outputs, recursive = T)),
+17
+)
+
+})
+
+# Delete files
+unlink(recipe$Run$output_dir, recursive = T)
diff --git a/tools/check_recipe.R b/tools/check_recipe.R
index abd0c8e2a24a9b311526a502103444557096f3d7..c105353455c984dd9d3e412b0dd0275c979f1624 100644
--- a/tools/check_recipe.R
+++ b/tools/check_recipe.R
@@ -296,6 +296,49 @@ check_recipe <- function(recipe) {
       }
     }
   }
+  # Indices
+  if ("Indices" %in% names(recipe$Analysis$Workflow)) {
+    nino_indices <- paste0("nino", c("1+2", "3", "3.4", "4"))
+    indices <- c("nao", nino_indices)
+    if (!("anomalies" %in% tolower(names(recipe$Analysis$Workflow)))) {
+      error(recipe$Run$logger,
+            paste0("Indices uses Anomalies as input, but Anomalies are missing",
+                   "in the recipe."))
+      error_status <- T
+    } else if (!(recipe$Analysis$Workflow$Anomalies$compute)) {
+        error(recipe$Run$logger,
+              paste0("Indices uses Anomalies as input, but the parameter",
+                     "'Anomalies:compute' is set as no/False."))
+      error_status <- T
+    }
+    recipe_indices <- tolower(names(recipe$Analysis$Workflow$Indices))
+    if (!all(recipe_indices %in% indices)) {
+      error(recipe$Run$logger,
+            paste0("Some of the indices under 'Indices' are not available.",
+                   "The available Indices are: 'NAO', 'Nino1+2', 'Nino3', ",
+                   "'Nino3.4' and 'Nino4'."))
+      error_status <- T
+    }
+    # Check that variables correspond with indices requested
+    recipe_variables <- strsplit(recipe$Analysis$Variables$name, ", | |,")[[1]]
+    if (("nao" %in% recipe_indices) &&
+        (!all(recipe_variables %in% c("psl", "z500")))) {
+      error(recipe$Run$logger,
+            paste0("It is not possible to compute the NAO with some of the ",
+                   "variables requested. To compute the NAO, please make sure",
+                   "your recipe requests only psl and/or z500."))
+      error_status <- T
+    }
+    if ((any(nino_indices %in% recipe_indices)) &&
+        (!all(recipe_variables %in% c("tos", "sst")))) {
+      error(recipe$Run$logger,
+            paste0("It is not possible to compute El Nino indices with some ",
+                   "of the variables requested. To compute El Nino, please ",
+                   "make sure your recipe requests only tos."))
+      error_status <- T
+    }
+  }
+
   # Skill
   if (("Skill" %in% names(recipe$Analysis$Workflow)) &&
       (is.null(recipe$Analysis$Workflow$Skill$metric))) {
diff --git a/tools/libs.R b/tools/libs.R
index e1179e8dc948372f18ba95a6909f48340f509df3..a67f9549fe2236d8c014ebf8f4af0bca111b6500 100644
--- a/tools/libs.R
+++ b/tools/libs.R
@@ -1,3 +1,4 @@
+# Libraries
 library(log4r)
 library(docopt)
 library(startR)
@@ -19,6 +20,7 @@ library(cowplot)
 library(stringr)
 library(pryr) # To check mem usage.
 
+# Functions
 ## To be removed when new package is done by library(CSOperational)
 source("tools/check_recipe.R")
 source("tools/prepare_outputs.R")
@@ -28,3 +30,7 @@ source("tools/read_atomic_recipe.R")
 source("tools/write_autosubmit_conf.R")
 source("tools/get_archive.R")
 # source("tools/add_dims.R") # Not sure if necessary yet
+
+# Settings
+options(bitmapType = 'cairo')
+pdf(NULL)