diff --git a/modules/Crossval/Crossval_NAO.R b/modules/Crossval/Crossval_NAO.R
new file mode 100644
index 0000000000000000000000000000000000000000..d83448b1263911b969d4b566a733e528c93497d0
--- /dev/null
+++ b/modules/Crossval/Crossval_NAO.R
@@ -0,0 +1,385 @@
+# Full-cross-val workflow
+## This code should be valid for individual months and temporal averages
+source("modules/Crossval/R/tmp/GetProbs.R")
+source("modules/Crossval/R/tmp/NAO.R")  
+source("modules/Crossval/R/tmp/ProjectField.R")
+source("modules/Crossval/R/tmp/EOF.R")
+source("modules/Crossval/R/tmp/Utils.R")
+
+Crossval_NAO <- function(recipe, data) {
+  cross.method <- recipe$Analysis$cross.method
+  # TODO move check
+  if (is.null(cross.method)) {
+    cross.method <- 'leave-one-out'
+  }
+  categories <- recipe$Analysis$Workflow$Probabilities$percentiles
+  categories <- lapply(categories, function (x) { 
+                       sapply(x, function(y) {
+                              eval(parse(text = y))})})
+  obsproj <- recipe$Analysis$Workflow$Indices$NAO$obsproj
+  if (is.null(obsproj)) {
+    obsproj <- TRUE
+  }
+  ncores <- recipe$Analysis$ncores
+  na.rm <- recipe$Analysis$remove_NAs 
+  ## data dimensions 
+  sdate_dim <- dim(data$hcst$data)['syear']
+  orig_dims <- names(dim(data$hcst$data))
+  # spatial dims
+  nlats <- dim(data$hcst$data)['latitude']
+  nlons <- dim(data$hcst$data)['longitude']
+  agg = 'region'
+
+  # TODO fix it to use new version https://earth.bsc.es/gitlab/external/cstools/-/blob/dev-cross-indices/R/CST_Calibration.R#L570
+  cross <- CSTools:::.make.eval.train.dexes(eval.method = cross.method,
+                                            amt.points = sdate_dim,
+                                            amt.points_cor = NULL) # k = ?
+  .cross_nao <- function(indices, exp, obs, cross, categories, recipe) {
+    ## Objects to return:
+    # lims_nao_hcst_tr: hcst percentiles (per category)
+    # lims_nao_obs_tr: obs percentiles (per category)
+    # ano_nao_ev: hcst cross-validated anomalies
+    # nao_obs_ev: obs cross-validated anomalies
+    # nao_obs_tr: obs training anomalies 
+ 
+    cross <- cross[[indices]]
+    print(paste("Crossval Index:", cross$eval.dexes))
+ 
+    # subset years: Subset works at BSC not at Athos
+    ## training indices
+    obs_tr <- Subset(data$obs$data, along = 'syear',
+   		     indices = cross$train.dexes)
+    hcst_tr <- Subset(data$hcst$data, along = 'syear', 
+  		      indices = cross$train.dexes)
+    ## evaluation indices
+    hcst_ev <- Subset(data$hcst$data, along = 'syear',
+  		      indices = cross$eval.dexes)
+    obs_ev <- Subset(data$obs$data, along = 'syear',
+  		     indices = cross$eval.dexes)
+    # compute climatology:
+    clim_obs_tr <- MeanDims(obs_tr, 'syear')
+    clim_hcst_tr <- MeanDims(hcst_tr, c('syear', 'ensemble'))
+    # compute anomalies:
+    ano_obs_tr <- s2dv::Ano(obs_tr, clim_obs_tr,
+                            ncores = ncores)
+    ano_hcst_tr <- s2dv::Ano(hcst_tr, clim_hcst_tr,
+                             ncores = ncores)
+    ano_hcst_ev <- s2dv::Ano(hcst_ev, clim_hcst_tr,
+                             ncores = ncores)
+    ano_obs_ev <- s2dv::Ano(obs_ev, clim_obs_tr,
+                            ncores = ncores)
+    # compute NAO:
+    nao <- NAO(exp = ano_hcst_tr, obs = ano_obs_tr, exp_cor = ano_hcst_ev,
+               ftime_avg = NULL, time_dim = 'syear',
+               memb_dim = 'ensemble',
+               space_dim = c('latitude', 'longitude'),
+               ftime_dim = 'time', obsproj = obsproj,
+               lat = data$obs$coords$latitude,
+               lon = data$obs$coords$longitude, 
+               ncores = recipe$Analysis$ncores)
+
+    nao_obs_ev <- NAO(exp = ano_hcst_tr, obs = ano_obs_tr, exp_cor = ano_obs_ev,
+                      ftime_avg = NULL, time_dim = 'syear',
+                      memb_dim = 'ensemble',
+                      space_dim = c('latitude', 'longitude'),
+                      ftime_dim = 'time', obsproj = obsproj,
+                      lat = data$obs$coords$latitude,
+                      lon = data$obs$coords$longitude,
+                      ncores = recipe$Analysis$ncores)$exp_cor
+    #Standarisation:
+    # Need the nao_hcst (for the train.dexes) to standarize the eval.dexes?
+    nao_hcst_ev <- Apply(list(nao$exp, nao$exp_cor),
+                         target_dims = c('syear', 'ensemble'),
+                         fun = function(x, y) {
+                           sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                           means <- mean(as.vector(x), na.rm = TRUE)
+                           res <- apply(y, c(1,2), function(z) {(z-means)/sd})},
+                         ncores = recipe$Analysis$ncores)$output1
+    nao_obs_ev <- Apply(list(nao$obs, nao_obs_ev),
+                        target_dims = c('syear','ensemble'),
+                        fun = function(x, y) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(y, c(1,2), function(z) {(z-means)/sd})},
+                        ncores = recipe$Analysis$ncores)$output1
+    nao_obs_tr <- Apply(list(nao$obs), target_dims = 'syear',
+                        fun = function(x) {
+                          sd <- sqrt(var(as.vector(x), na.rm = TRUE))
+                          means <- mean(as.vector(x), na.rm = TRUE)
+                          res <- apply(x, 1, function(z) {(z-means)/sd})},
+                        ncores = recipe$Analysis$ncores,
+                        output_dims = 'syear')$output1
+    nao_hcst_tr <- 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(z) {(z-means)/sd})},
+                          ncores = recipe$Analysis$ncores)$output1
+    # compute category limits
+    lims_nao_hcst_tr <- Apply(nao_hcst_tr, target_dims = c('syear', 'ensemble'),
+                              fun = function(x, prob_lims) {
+                                lapply(prob_lims, function(ps) {
+                                  as.array(quantile(as.vector(x), ps,
+                                                    na.rm = TRUE))})},
+                              prob_lims = list(unlist(categories)),
+			      output_dims = "bin",
+                              ncores = ncores)$output1
+    lims_nao_obs_tr <- Apply(nao_obs_tr, target_dims = c('syear', 'ensemble'),
+                             fun = function(x, prob_lims) {
+                               lapply(prob_lims, function(ps) {
+                                 as.array(quantile(as.vector(x), ps,
+                                                   na.rm = TRUE))})},
+                             prob_lims = list(unlist(categories)),
+			     output_dims = "bin",
+                             ncores = ncores)$output1
+    # 
+    no_dims <- which(names(dim(exp)) %in%
+			 c("latitude", "longitude", "syear"))
+    nao_hcst_ev <- MergeDims(nao_hcst_ev, c('syear', 'ensemble'), 'ensemble')
+    nao_hcst_ev <- Reorder(nao_hcst_ev, names(dim(exp))[-no_dims])
+    nao_obs_ev <- MergeDims(nao_obs_ev, c('syear', 'ensemble'), 'ensemble')
+    nao_obs_ev <- Reorder(nao_obs_ev, names(dim(obs))[-no_dims])
+    nao_obs_tr <- MergeDims(nao_obs_tr, c('syear', 'ensemble'), 'ensemble')
+    nao_obs_tr <- Reorder(nao_obs_tr, names(dim(obs))[-no_dims])
+    return(list(nao_hcst_ev = nao_hcst_ev,
+                nao_obs_ev = nao_obs_ev,
+                nao_obs_tr = nao_obs_tr,
+                lims_nao_hcst_tr = lims_nao_hcst_tr,
+                lims_nao_obs_tr = lims_nao_obs_tr))
+  }
+  indices <- array(1:length(cross), dim = c(syear = length(cross)))
+  res <- Apply(data = list(indices = indices,
+                           exp = data$hcst$data,
+                           obs = data$obs$data),
+               target_dims = list(indices = NULL,
+                                  obs = names(dim(data$obs$data)),
+                                  exp = names(dim(data$hcst$data))),
+               fun = .cross_nao,
+               ncores = 1, ## TODO: Explore options for core distribution
+               cross = cross,
+               categories = categories,
+               recipe = recipe)
+
+  info(recipe$Run$logger,
+       "#### NAO Cross-validation loop ended #####")
+  # Rearrange limits into a list:
+  lims_tmp_hcst <- list()
+  lims_tmp_obs <- list()
+  last_index <- 0
+  for (ps in 1:length(categories)) {
+    indices <- last_index + 1:length(categories[[ps]])
+    lims_category_hcst <- Subset(res$lims_nao_hcst_tr, along = "bin",
+                                 indices = list(indices),
+                                 drop = FALSE)
+    lims_tmp_hcst[[ps]] <- lims_category_hcst
+    lims_category_obs <- Subset(res$lims_nao_obs_tr, along = "bin",
+                                indices = list(indices),
+                                drop = FALSE)
+    lims_tmp_obs[[ps]] <- lims_category_obs
+    last_index <- indices[length(indices)]
+  }
+  res$lims_nao_hcst_tr <- lims_tmp_hcst
+  res$lims_nao_obs_tr <- lims_tmp_obs
+  rm(lims_tmp_hcst, lims_tmp_obs)
+  gc()
+
+  # Forecast anomalies:
+  nao_fcst <- NULL
+  if (!is.null(data$fcst)) {
+    clim_hcst <- Apply(data$hcst$data,
+                       target_dims = c('syear', 'ensemble'),
+                       mean,
+                       na.rm = na.rm,
+                       ncores = ncores)$output1
+    clim_obs <- Apply(data$obs$data,
+                       target_dims = c('syear', 'ensemble'),
+                       mean,
+                       na.rm = na.rm,
+                       ncores = ncores)$output1
+    ano_fcst <- Ano(data = data$fcst$data, clim = clim_hcst)
+    ano_hcst <- Ano(data = data$hcst$data, clim = clim_hcst)
+    ano_obs <- Ano(data= data$obs$data, clim = clim_obs)
+    nao_fcst <- NAO(exp = ano_hcst, obs = ano_obs, exp_cor = ano_fcst,
+                    ftime_avg = NULL, time_dim = 'syear',
+                    memb_dim = 'ensemble',
+                    space_dim = c('latitude', 'longitude'),
+                    ftime_dim = 'time', obsproj = obsproj,
+                    lat = data$obs$coords$latitude,
+                    lon = data$obs$coords$longitude,
+                    ncores = recipe$Analysis$ncores)
+    # Terciles limits using the whole hindcast period:
+    lims_fcst <- Apply(nao_fcst$exp, target_dims = c('syear', 'ensemble'),
+                       fun = function(x, prob_lims) {
+                         lapply(prob_lims, function(ps) {
+                         quantile(as.vector(x), ps, na.rm = na.rm)})},
+                       output_dims = lapply(categories, function(x) {'cat'}),
+                       prob_lims = categories,
+                       ncores = ncores)
+  }
+
+  fcst_probs <- lapply(categories, function(x){NULL})
+  hcst_probs_ev <- lapply(categories, function(x){NULL})
+  obs_probs_ev <- lapply(categories, function(x){NULL})
+
+  # Compute Probabilities
+  for (ps in 1:length(categories)) {
+    hcst_probs_ev[[ps]] <- GetProbs(res$nao_hcst_ev, time_dim = 'syear',
+                              prob_thresholds = NULL,
+                              bin_dim_abs = 'bin',
+                              indices_for_quantiles = NULL,
+                              memb_dim = 'ensemble',
+                              abs_thresholds = res$lims_nao_hcst_tr[[ps]],
+                              ncores = ncores)
+    obs_probs_ev[[ps]] <- GetProbs(res$nao_obs_ev, time_dim = 'syear',
+                                   prob_thresholds = NULL,
+                                   bin_dim_abs = 'bin',
+                                   indices_for_quantiles = NULL,
+                                   memb_dim = 'ensemble',
+                                   abs_thresholds = res$lims_nao_obs_tr[[ps]],
+                                   ncores = ncores)
+    if (!is.null(data$fcst)) {
+      fcst_probs[[ps]] <- GetProbs(nao_fcst$exp_cor, time_dim = 'syear',
+                                   prob_thresholds = NULL,
+                                   bin_dim_abs = 'bin',
+                                   indices_for_quantiles = NULL,
+                                   memb_dim = 'ensemble',
+                                   abs_thresholds = lims_fcst[[ps]],
+                                   ncores = ncores)
+    }
+  }
+  # Convert to s2dv_cubes
+  nao_hcst <- data$hcst
+  nao_hcst$data <- res$nao_hcst_ev
+  nao_obs <- data$obs
+  nao_obs$data <- res$nao_obs_ev
+
+  info(recipe$Run$logger, 
+       "#### NAO  and Probabilities Done #####") 
+  if (recipe$Analysis$Workflow$Indices$NAO$save != 'none') {
+    info(recipe$Run$logger, "##### START SAVING NAO #####")
+    recipe$Run$output_dir <- paste0(recipe$Run$output_dir,
+                                    "/outputs/Anomalies/")
+    # Save forecast
+    if ((recipe$Analysis$Workflow$Indices$NAO$save %in% 
+           c('all', 'exp_only', 'fcst_only')) && !is.null(data$fcst)) {
+      save_forecast(recipe = recipe, data_cube = nao_fcst, type = 'fcst')
+    }
+    # Save hindcast
+    if (recipe$Analysis$Workflow$Indices$save %in%
+        c('all', 'exp_only')) {
+      save_forecast(recipe = recipe, data_cube = nao_hcst, type = 'hcst')
+    }
+    # Save observation
+    if (recipe$Analysis$Workflow$Indices$NAO$save == 'all') {
+      save_observations(recipe = recipe, data_cube = nao_obs)
+    }
+  }
+  # Save probability bins
+  probs_hcst <- list()
+  probs_fcst <- list()
+  probs_obs <- list()
+  all_names <- NULL
+  # Make categories rounded number to use as names:
+  categories <- recipe$Analysis$Workflow$Probabilities$percentiles
+  categories <- lapply(categories, function (x) {
+                       sapply(x, function(y) {
+                              round(eval(parse(text = y)),2)})})
+
+  for (ps in 1:length(categories)) {
+    for (perc in 1:(length(categories[[ps]]) + 1)) {
+      if (perc == 1) {
+        name_elem <- paste0("below_", categories[[ps]][perc])
+      } else if (perc == length(categories[[ps]]) + 1) {
+        name_elem <- paste0("above_", categories[[ps]][perc-1])
+      } else {
+        name_elem <- paste0("from_", categories[[ps]][perc-1], 
+                            "_to_", categories[[ps]][perc])
+      }
+      probs_hcst <- append(list(Subset(hcst_probs_ev[[ps]],
+                             along = 'bin', indices = perc, drop = 'all')),
+                           probs_hcst)
+      probs_obs <- append(list(Subset(obs_probs_ev[[ps]],
+                              along = 'bin', indices = perc, drop = 'all')),
+                          probs_obs)
+      if (!is.null(data$fcst)) {
+        probs_fcst <- append(list(Subset(fcst_probs[[ps]],
+                              along = 'bin', indices = perc, drop = 'all')),
+                             probs_fcst)
+      }
+      all_names <- c(all_names, name_elem)
+    }
+  }
+  names(probs_hcst) <- all_names
+  if (!('var' %in% names(dim(probs_hcst[[1]])))) {
+    probs_hcst <- lapply(probs_hcst, function(x) {
+                     dim(x) <- c(var = 1, dim(x))
+                     return(x)})
+  }
+  names(probs_obs) <- all_names
+  if (!('var' %in% names(dim(probs_obs[[1]])))) {
+    probs_obs <- lapply(probs_obs, function(x) {
+                     dim(x) <- c(var = 1, dim(x))
+                     return(x)})
+  }
+
+  if (!is.null(data$fcst)) {
+    names(probs_fcst) <- all_names
+    if (!('var' %in% names(dim(probs_fcst[[1]])))) {
+      probs_fcst <- lapply(probs_fcst, function(x) {
+                     dim(x) <- c(var = 1, dim(x))
+                     return(x)})
+    }
+    if (!('syear' %in% names(dim(probs_fcst[[1]])))) {
+      probs_fcst <- lapply(probs_fcst, function(x) {
+                     dim(x) <- c(syear = 1, dim(x))
+                     return(x)})
+    }
+  }
+  if (recipe$Analysis$Workflow$Probabilities$save %in%
+    c('all', 'bins_only')) {
+    save_probabilities(recipe = recipe, probs = probs_hcst,
+                       data_cube = data$hcst, agg = agg,
+                       type = "hcst")
+    save_probabilities(recipe = recipe, probs = probs_obs,
+                       data_cube = data$hcst, agg = agg,
+                       type = "obs")
+    # TODO Forecast
+    if (!is.null(probs_fcst)) {
+      save_probabilities(recipe = recipe, probs = probs_fcst,
+                         data_cube = data$fcst, agg = agg, 
+                         type = "fcst")
+    }
+  }
+  return(list(hcst = nao_hcst, obs = nao_obs, fcst = nao_fcst,
+              # Mean Bias for the NAO doesn't make sense, so
+              # not returning fullvals
+              #hcst.full_val = data$hcst, obs.full_val = data$obs,
+              cat_lims = list(hcst_tr = res$lims_nao_hcst_tr,
+                              obs_tr = res$lims_nao_obs_tr),
+              probs = list(hcst_ev = hcst_probs_ev,
+                           obs_ev = obs_probs_ev),
+              ref_obs_tr = res$nao_obs_tr))
+}
+
+
+## The result contains the inputs for Skill_full_crossval.
+## this is a list with the required elements:
+  ## probs is a list with
+    ## probs$hcst_ev and probs$obs_ev 
+    ## probs$hcst_ev will have as many elements in the $Probabilities$percentiles
+    ## each element will be an array with 'cat' dimension
+    ## the same for probs$obs_ev
+  ## hcst is a s2dv_cube for the post-processed hindcast for the evalutaion indices
+    ## in this case cross validated anomalies
+  ## obs is a s2dv_cube for the post-processed obs 
+    ## in this case cross validated anomalies
+  ## fcst is a s2dv_cube for the post-processed fcst
+    ## in this case cross anomalies with the full hindcast period
+    ## this object is not required for skill assessment
+  ## hcst.full_val and obs.full_val are the original data to compute mean bias
+  ## cat_lims used to compute the probabilities
+    ## this object is not required for skill assessment
+  ## ref_obs_tr is an array with the cross-validate observed anomalies
+    ## to be used as reference forecast in the CRPSS and CRPS_clim
+    ## it is computed from the training indices
+
diff --git a/modules/Visualization/R/plot_deterministic_forecast.R b/modules/Visualization/R/plot_deterministic_forecast.R
new file mode 100644
index 0000000000000000000000000000000000000000..30f13b78db1654bcb35fa7a86b6f4ad93d7c6393
--- /dev/null
+++ b/modules/Visualization/R/plot_deterministic_forecast.R
@@ -0,0 +1,181 @@
+
+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)) {
+    if (all(is.na(fcst)) && all(is.na(obs))) {
+      # No data, return random ylims
+      ylims <- c(-1, 1)
+    } else {
+      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/Visualization/plot_NAO.R b/modules/Visualization/plot_NAO.R
new file mode 100644
index 0000000000000000000000000000000000000000..9eb6286984308cc05fe1cfb4c9609b4bd715fc7f
--- /dev/null
+++ b/modules/Visualization/plot_NAO.R
@@ -0,0 +1,307 @@
+plot_NAO <- function(recipe, nao, data) {
+
+# 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
+  plot_ts <- TRUE
+  plot_sp <- TRUE
+  alpha <- 0.05
+  lons <- data$hcst$coords$longitude
+  lats <- data$hcst$coords$latitude
+   nao_region <- c(lonmin = -80, lonmax = 40,
+                  latmin = 20, latmax = 80)
+
+  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)
+  }
+
+  if (plot_ts) {
+    dir.create(paste0(recipe$Run$output_dir, "/plots/Indices/"),
+               showWarnings = F, recursive = T)
+    source("modules/Visualization/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)
+        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), ".pdf") 
+        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)
+        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), ".pdf") 
+        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, ".pdf")
+        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, ".pdf")
+        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, ".pdf")
+        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, ".pdf")
+        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, ".pdf")
+        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, ".pdf")
+        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
+  }
+}