diff --git a/conf/archive_subseasonal.yml b/conf/archive_subseasonal.yml
new file mode 100644
index 0000000000000000000000000000000000000000..ed7027da15e020d6585da4f08cc2aeb4f5ac7309
--- /dev/null
+++ b/conf/archive_subseasonal.yml
@@ -0,0 +1,112 @@
+# TO BE REVIEWED
+
+esarchive:
+    src: "/esarchive/"
+    System:
+        NCEP-CFSv2:
+            name: "NCEP CFSv2"
+            institution: "NOAA NCEP" #?
+            src: "exp/ncep/cfs-v2/"
+            monthly_mean: {"tas":"monthly_mean/tas_f6h/", "prlr":"monthly_mean/prlr_f6h/",
+                           "tasmax":"monthly_mean/tasmax_f6h/", "tasmin":"monthly_mean/tasmin_f6h/"}
+            nmember:
+              fcst: 48
+              hcst: 12
+            calendar: "gregorian"
+            time_stamp_lag: "0"
+            reference_grid: "conf/grid_description/griddes_ncep-cfsv2.txt"
+    Reference:
+        ERA5:
+            name: "ERA5"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "recon/ecmwf/era5/"
+            daily_mean: {"tas":"daily_mean/tas_f1h-r1440x721cds/",
+                         "rsds":"daily_mean/rsds_f1h-r1440x721cds/",
+                         "prlr":"daily_mean/prlr_f1h-r1440x721cds/",
+                         "g300":"daily_mean/g300_f1h-r1440x721cds/",
+                         "g500":"daily_mean/g500_f1h-r1440x721cds/",
+                         "g850":"daily_mean/g850_f1h-r1440x721cds/",
+                         "sfcWind":"daily_mean/sfcWind_f1h-r1440x721cds/",
+                         "tasmax":"daily/tasmax_f1h-r1440x721cds/",
+                         "tasmin":"daily/tasmin_f1h-r1440x721cds/",
+                         "ta300":"daily_mean/ta300_f1h-r1440x721cds/",
+                         "ta500":"daily_mean/ta500_f1h-r1440x721cds/", 
+                         "ta850":"daily_mean/ta850_f1h-r1440x721cds/",
+                         "hurs":"daily_mean/hurs_f1h-r1440x721cds/"}
+            monthly_mean: {"tas":"monthly_mean/tas_f1h-r1440x721cds/",
+                           "psl":"monthly_mean/psl_f1h-r1440x721cds/",
+                           "prlr":"monthly_mean/prlr_f1h-r1440x721cds/",
+                           "rsds":"monthly_mean/rsds_f1h-r1440x721cds/",
+                           "g300":"monthly_mean/g300_f1h-r1440x721cds/",
+                           "g500":"monthly_mean/g500_f1h-r1440x721cds/",
+                           "g850":"monthly_mean/g850_f1h-r1440x721cds/",
+                           "sfcWind":"monthly_mean/sfcWind_f1h-r1440x721cds/",
+                           "tasmax":"monthly_mean/tasmax_f1h-r1440x721cds/",
+                           "tasmin":"monthly_mean/tasmin_f1h-r1440x721cds/",
+                           "ta300":"montly_mean/ta300_f1h-r1440x721cds/",
+                           "ta500":"monthly_mean/ta500_f1h-r1440x721cds/", 
+                           "ta850":"monthly_mean/ta850_f1h-r1440x721cds/",
+                           "tos":"monthly_mean/tos_f1h-r1440x721cds/",
+                           "sic":"monthly_mean/sic_f1h-r1440x721cds/"}
+            calendar: "standard"
+            reference_grid: "/esarchive/recon/ecmwf/era5/monthly_mean/tas_f1h-r1440x721cds/tas_201805.nc"
+            land_sea_mask: "/esarchive/recon/ecmwf/era5/constant/lsm-r1440x721cds/sftof.nc"
+        ERA5-Land:
+            name: "ERA5-Land"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "recon/ecmwf/era5land/"
+            daily_mean: {"tas":"daily_mean/tas_f1h/", "rsds":"daily_mean/rsds_f1h/",
+                         "prlr":"daily_mean/prlr_f1h/", "sfcWind":"daily_mean/sfcWind_f1h/",
+                         "tasmin":"daily/tasmin/", "tasmax":"daily/tasmax/"}
+            monthly_mean: {"tas":"monthly_mean/tas_f1h/","tasmin":"monthly_mean/tasmin_f24h/",
+                           "tasmax":"monthly_mean/tasmax_f24h/", "prlr":"monthly_mean/prlr_f1h/",
+                           "sfcWind":"monthly_mean/sfcWind_f1h/", "rsds":"monthly_mean/rsds_f1h/",
+                           "tdps":"monthly_mean/tdps_f1h/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/recon/ecmwf/era5land/daily_mean/tas_f1h/tas_201805.nc"
+        UERRA:
+            name: "ECMWF UERRA"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "recon/ecmwf/uerra_mescan/"
+            daily_mean: {"tas":"daily_mean/tas_f6h/"}
+            monthly_mean: {"tas":"monthly_mean/tas_f6h/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/recon/ecmwf/uerra_mescan/daily_mean/tas_f6h/tas_201805.nc"
+        CERRA:
+            name: "ECMWF CERRA"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "recon/ecmwf/cerra/"
+            daily_mean: {"hurs":"daily_mean/hurs_f3h-r2631x1113/", "ps":"daily_mean/ps_f3h-r2631x1113/",
+                         "sfcWind":"daily_mean/sfcWind_f3h-r2631x1113/",
+                         "tas":"daily_mean/tas_f3h-r2631x1113/", "winddir":"daily_mean/tas_f3h-r2631x1113/"}
+            monthly_mean: {"hurs":"monthly_mean/hurs_f3h-r2631x1113/", "ps":"monthly_mean/ps_f3h-r2631x1113/",
+                           "sfcWind":"monthly_mean/sfcWind_f3h-r2631x1113/",
+                           "tas":"monthly_mean/tas_f3h-r2631x1113/",
+                           "winddir":"monthly_mean/winddir_f3h-r2631x1113/",
+                           "tasmin":"monthly_mean/tasmin_f24h-r2631x1113/",
+                           "tasmax":"monthly_mean/tasmax_f24h-r2631x1113/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/recon/ecmwf/cerra/monthly_mean/tas_f3h-r2631x1113/tas_200506.nc"
+        CERRA-Land:
+            name: "ECMWF CERRA-Land"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "recon/ecmwf/cerraland/"
+            daily_mean: {"prlr":"daily_mean/prlr_f6h-r2631x1113/"}
+            monthly_mean: {"prlr":"monthly_mean/prlr_f6h-r2631x1113/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/recon/ecmwf/cerraland/monthly_mean/prlr_f6h-r2631x1113/prlr_200412.nc"
+        HadCRUT5:
+            name: "HadCRUT5"
+            institution: "Met Office"
+            src: "obs/ukmo/hadcrut_v5.0_analysis/"
+            monthly_mean: {"tasanomaly":"monthly_mean/tasanomaly/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/obs/ukmo/hadcrut_v5.0_analysis/monthly_mean/tasanomaly/tasanomaly_202001.nc"
+        BEST:
+            name: "BEST"
+            institution: "European Centre for Medium-Range Weather Forecasts"
+            src: "obs/berkeleyearth/berkeleyearth/"
+            daily_mean: {"tas":"daily_mean/tas/"}
+            monthly_mean: {"tas":"monthly_mean/tas/"}
+            calendar: "proleptic_gregorian"
+            reference_grid: "/esarchive/obs/berkeleyearth/berkeleyearth/monthly_mean/tas/tas_201805.nc"
diff --git a/modules/Visualization/R/plot_ensemble_mean.R b/modules/Visualization/R/plot_ensemble_mean.R
index ff91adf4cef2d5b453f5a6fa61ca06d7310883e8..91e5cc2be89f515f3ab6f89f7471c77f81fb1dcc 100644
--- a/modules/Visualization/R/plot_ensemble_mean.R
+++ b/modules/Visualization/R/plot_ensemble_mean.R
@@ -7,7 +7,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
   latitude <- fcst$coords$lat
   longitude <- fcst$coords$lon
   archive <- get_archive(recipe)
-  if (recipe$Analysis$Datasets$System$name == 'Multimodel'){
+  if (recipe$Analysis$Datasets$System$name == 'Multimodel') {
     system_name <- paste0('Multimodel-',
                           recipe$Analysis$Datasets$Multimodel$approach)
   } else {
@@ -18,9 +18,12 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
   if (tolower(recipe$Analysis$Horizon) == "seasonal") {
     init_month <- as.numeric(substr(recipe$Analysis$Time$sdate,
                                     start = 1, stop = 2))
+  } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+    init_week <- recipe$Analysis$Time$sdate
   } else {
     ## TODO: Sort out decadal initial month (is it always January?)
-    init_month <- 1
+    init_month <- 1 
+    init_week <- 1
   }
   if (!is.null(recipe$Analysis$Workflow$Visualization$projection)) {
     projection <- tolower(recipe$Analysis$Workflow$Visualization$projection)
@@ -40,7 +43,6 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                                                     "sday", "sweek"),
                                           indices = list(1, var, 1, 1),
                                           drop = 'selected')
-
     var_ens_mean <- Reorder(var_ens_mean, c("syear",
                                             "time",
                                             "longitude",
@@ -79,21 +81,13 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         cols <- col_fun(length(brks) - 1)
       }
     }
-
+    
     for (i_syear in start_date) {
-      if (length(start_date) == 1) {
-        i_var_ens_mean <- ClimProjDiags::Subset(var_ens_mean,
+      i_var_ens_mean <- ClimProjDiags::Subset(var_ens_mean,
                                                 along = c("syear"),
                                                 indices = which(start_date == i_syear),
                                                 drop = 'selected')
-        outfile <- paste0(outdir[[var]], "forecast_ensemble_mean-", start_date)
-      } else {
-        i_var_ens_mean <- ClimProjDiags::Subset(var_ens_mean,
-                                                along = c("syear"),
-                                                indices = which(start_date == i_syear),
-                                                drop = 'selected')
-        outfile <- paste0(outdir[[var]], "forecast_ensemble_mean-", i_syear)
-      }
+      outfile <- paste0(outdir[[var]], "forecast_ensemble_mean-", i_syear)
       # Mask
       if (!is.null(mask)) {
         outfile <- paste0(outfile, "_enscormask")
@@ -116,56 +110,83 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
         var_dots <- as.numeric(var_dots <= 0)
         dim(var_dots) <- dim_dots
       }
-      toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
-                         "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear)
+      if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                           "\n", "Forecast Ensemble Mean / ", "Issued on ", 
+                           format(ymd(i_syear), "%d-%m-%Y"))
+      } else {
+        toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                           "\n", "Forecast Ensemble Mean / ", "Init.: ", i_syear)
+      }
       if (is.null(attributes(fcst$attrs$time_bounds))) {
-        months <- lubridate::month(
-                    fcst$attrs$Dates[1, 1, which(start_date == i_syear), ],
-                    label = T, abb = F)
+        if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+          time_labels <- lubridate::month( ## is time_labels an appropriate name? 
+            fcst$attrs$Dates[1, 1, which(start_date == i_syear), ],
+            label = T, abb = F)
+        } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+          time_labels <- fcst$attrs$Dates[1, 1, which(start_date == i_syear), ] 
+          monday <- ymd_hms(time_labels) - days(wday(ymd_hms(time_labels), week_start = 1) - 1)
+          sunday <- monday + days(6)
+          time_labels <- paste0("Valid from ", format(monday,"%d-%m"), " to ", format(sunday, "%d-%m"))
+        }
       } else {
         if (length(attributes(fcst$attrs$time_bounds)$plotting_attr) > 1) {
-          months <- unlist(
-          lapply(1:length(fcst$attrs$time_bounds$start), function(i) {
-            ftime_ini <- attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i]
-            ftime_end <- attributes(fcst$attrs$time_bounds)$plotting_attr$end[i]
-            # labels for file name:
-            ftime <- paste0(ftime_ini, "-", ftime_end)
-            # title names:
-            ftime_ini <- init_month + ftime_ini - 1
-            ftime_ini <- ifelse(ftime_ini > 12, ftime_ini - 12, ftime_ini)
-            ftime_ini <- month.name[ftime_ini]
-            ftime_end <- init_month + ftime_end - 1
-            ftime_end <- ifelse(ftime_end > 12, ftime_end - 12, ftime_end)
-            ftime_end <- month.name[ftime_end]
-            toptitle <- paste(ftime_ini, "to", ftime_end)}))
+          time_labels <- unlist(
+            lapply(1:length(fcst$attrs$time_bounds$start), function(i) {
+              ftime_ini <- attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i]
+              ftime_end <- attributes(fcst$attrs$time_bounds)$plotting_attr$end[i]
+              # labels for file name:
+              ftime <- paste0(ftime_ini, "-", ftime_end)
+              # title names:
+              if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+                ftime_ini <- init_month + ftime_ini - 1
+                ftime_ini <- ifelse(ftime_ini > 12, ftime_ini - 12, ftime_ini)
+                ftime_ini <- month.name[ftime_ini]
+                ftime_end <- init_month + ftime_end - 1
+                ftime_end <- ifelse(ftime_end > 12, ftime_end - 12, ftime_end)
+                ftime_end <- month.name[ftime_end]
+                toptitle <- paste(ftime_ini, "to", ftime_end)
+              } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+                ## set to make ftime_ini be a Monday
+                ## if init_week always Thursday, the 2nd addend can be simplified to +4
+                ftime_ini <- ymd(init_week) + 
+                             ((8 - wday(ymd(init_week), week_start = 1)) %% 7) + 
+                             weeks(ftime_ini)
+                ## set to make ftime_end be a Sunday
+                ftime_end <- ymd(init_week) + 
+                             ((8 - wday(ymd(init_week), week_start = 1)) %% 7) + 
+                             weeks(ftime_end) + 6
+                toptitle <- paste("Valid from", ftime_ini, "to", ftime_end)
+                }            
+              })) 
         } else {
-          months <- attributes(fcst$attrs$time_bounds)$plotting_attr[[1]]
+          time_labels <- attributes(fcst$attrs$time_bounds)$plotting_attr[[1]]
         }
       }
       years <- lubridate::year(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ])
-
+      
       if (recipe$Analysis$Workflow$Visualization$multi_panel) { 
         # Define name of output file and titles
-        titles <- as.vector(months)
+        if (tolower(recipe$Analysis$Horizon) %in% c("seasonal", "subseasonal")) {
+          titles <- as.vector(time_labels)
+        } else {
+          titles <- NULL
+        }
         # Plots
-        PlotLayout(PlotEquiMap, c('longitude', 'latitude'),
-                   i_var_ens_mean, longitude, latitude,
-                   mask = mask,
-                   dots = dots,
-                   filled.continents = F,
-                   toptitle = toptitle,
-                   title_scale = 0.7,
-                   subtitle_scale = 0.8,
-                   subtitle_margin_scale = 2,
-                   titles = titles,
-                   units = units,
-                   cols = cols,
-                   brks = brks,
-                   fileout = paste0(outfile, ".pdf"),
-                   bar_label_digits = 4,
-                   extra_margin = rep(1, 4),
-                   bar_label_scale = 1.5,
-                   axes_label_scale = 1.1)
+        output_configuration <- output_conf$Multipanel$forecast_ensemble_mean
+        base_args <- list(fun = "PlotEquiMap",
+                          plot_dims = c('longitude', 'latitude'),
+                          var = i_var_ens_mean, lon = longitude,
+                          lat = latitude, mask = mask, dots = dots,
+                          filled.continents = FALSE, toptitle = toptitle,
+                          title_scale = 0.7, subtitle_scale = 0.8,
+                          subtitle_margin_scale = 2, titles = titles,
+                          units = units, cols = cols, brks = brks,
+                          fileout = paste0(outfile, ".pdf"),
+                          bar_label_digits = 4, extra_margin = rep(1, 4),
+                          bar_label_scale = 1.5, axes_label_scale = 1.1)
+        base_args[names(output_configuration)] <- output_configuration
+        do.call(PlotLayout, base_args)
       } else {
         # Define function and parameters depending on projection
         if (projection == 'cylindrical_equidistant') {
@@ -194,17 +215,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                             style = 'point', brks = brks, cols = cols)
         }
         # Loop over forecast times
-        for (i in 1:length(months)) {
-          # Get forecast time label
-          if (is.null(attributes(fcst$attrs$time_bounds))) {
-            forecast_time <- match(months[i], month.name) - init_month + 1
-            if (forecast_time < 1) {
-              forecast_time <- forecast_time + 12
-            }
-            forecast_time <- sprintf("%02d", forecast_time)
-          } else {
-            forecast_time <- months
-          }
+        for (i in 1:length(time_labels)) {  
           # Get mask subset
           if (!is.null(mask)) {
             mask_i <- Subset(var_mask, along = 'time', indices = i, drop = TRUE)
@@ -217,20 +228,33 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
           } else {
             dots_i <- NULL
           }
-         # Define plot title
+          # Define plot title
           if (tolower(recipe$Analysis$Horizon) == 'seasonal') {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
                                "\n", "Ensemble Mean / ",
-                               months[i], " ", years[i],
+                               time_labels[i], " ", years[i],
                                " / Start date: ",
                                format(as.Date(i_syear, format="%Y%m%d"),
-                                              "%d-%m-%Y"))
+                                      "%d-%m-%Y"))
+          } else if (tolower(recipe$Analysis$Horizon) == 'subseasonal') {
+            toptitle <- paste0(system_name, " / ",
+                               str_to_title(var_long_name),
+                               "\n", "Ensemble Mean / ",
+                               "Issued on ", 
+                               format(as.Date(i_syear, format="%Y%m%d"), 
+                                      "%d-%m-%Y"),
+                               "\n", time_labels[i], " of ",  years[i])
+                               ## "Week: ", time_labels[i], 
+                               ## " of ", years[i],
+                               ## " / Start date: ",
+                               ## format(as.Date(i_syear, format="%Y%m%d"),
+                               ##        "%d-%m-%Y"))
           } else {
             toptitle <- paste0(system_name, " / ",
                                str_to_title(var_long_name),
                                "\n", "Ensemble Mean / ",
-                               months[i], 
+                               months[i],
                                " / Start date: ",
                                i_syear)
           }
@@ -238,9 +262,10 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
           if (identical(fun, PlotRobinson)) {
             ## TODO: Customize technical details
             base_args[['caption']] <-
-              paste0("Nominal start date: ", start_date, "\n",
-                     "Forecast month: ", sprintf("%02d", i), "\n",
-                     "Reference: ", recipe$Analysis$Datasets$Reference)
+              paste0("Nominal start date: ", ymd(start_date), "\n",
+                     "Forecast week: ", sprintf("%02d", i), "\n",
+                     "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
+                     "Units: ", units)
           }
           # Modify base arguments
           base_args[[1]] <- i_var_ens_mean[i, , ]
@@ -249,8 +274,8 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
           } else {
             if (length(attributes(fcst$attrs$time_bounds)$plotting_attr) > 1) {
               fileout <- paste0(outfile, "_ft",
-                attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i], "-",
-                attributes(fcst$attrs$time_bounds)$plotting_attr$end[i], ".pdf")
+                                attributes(fcst$attrs$time_bounds)$plotting_attr$ini[i], "-",
+                                attributes(fcst$attrs$time_bounds)$plotting_attr$end[i], ".pdf")
             } else {
               fileout <- paste0(outfile, "_ft", months[i], ".pdf")
             }
@@ -262,7 +287,7 @@ plot_ensemble_mean <- function(recipe, fcst, mask = NULL, dots = NULL, outdir, o
                   args = c(base_args,
                            list(toptitle = toptitle,
                                 fileout = fileout)))
-        }
+        }        
       }
     }
   }
diff --git a/modules/Visualization/R/plot_metrics.R b/modules/Visualization/R/plot_metrics.R
index 2df06e8bdf87027e3a42da24b78821c5bdc1de82..5b71523bdf285edbae0ad670b4454ad077f2249a 100644
--- a/modules/Visualization/R/plot_metrics.R
+++ b/modules/Visualization/R/plot_metrics.R
@@ -1,4 +1,5 @@
 library(stringr)
+library(lubridate)
 
 plot_metrics <- function(recipe, data_cube, metrics,
                           outdir, significance = F, output_conf) {
@@ -8,8 +9,7 @@ plot_metrics <- function(recipe, data_cube, metrics,
   # metrics: list of named metric arrays with named dimensions
   # outdir: output directory
   # significance: T/F, whether to display the significance dots in the plots
-
-  # Abort if frequency is daily
+   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
     error(recipe$Run$logger, "Visualization functions not yet implemented
           for daily data.")
@@ -20,8 +20,8 @@ plot_metrics <- function(recipe, data_cube, metrics,
     stop("The element 'metrics' must be a list of named arrays.")
   }
   
-  latitude <- data_cube$coords$lat
-  longitude <- data_cube$coords$lon
+  latitude <- data_cube$coords$latitude
+  longitude <- data_cube$coords$longitude
   archive <- get_archive(recipe)
   if (recipe$Analysis$Datasets$System$name == 'Multimodel'){
     system_name <- paste0('Multimodel-',
@@ -31,19 +31,35 @@ plot_metrics <- function(recipe, data_cube, metrics,
   }
   hcst_period <- paste0(recipe$Analysis$Time$hcst_start, "-",
                         recipe$Analysis$Time$hcst_end)
+  start_date <- recipe$Analysis$Time$sdate
   if (tolower(recipe$Analysis$Horizon) == "seasonal") {
     init_month <- as.numeric(substr(recipe$Analysis$Time$sdate,
                                     start = 1, stop = 2))
-  } else {
-    ## TODO: Sort out decadal initial month (is it always January?)
-    init_month <- 1
-  }
-  month_label <- tolower(month.name[init_month])
-  month_abbreviation <- month.abb[init_month]
-  # Get months
-  months <- lubridate::month(Subset(data_cube$attrs$Dates,
+    months <- lubridate::month(Subset(data_cube$attrs$Dates,
                                     "syear", indices = 1),
                              label = T, abb = F,locale = "en_GB")
+    month_label <- tolower(month.name[init_month])
+    month_abbreviation <- month.abb[init_month]
+  } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+    init_week <- recipe$Analysis$Time$sdate
+    weeks <- Subset(data_cube$attr$Dates, along = 'sweek', 
+                    indices = (recipe$Analysis$Workflow$Skill$sweek_window +1)/2)    
+    weeks <- Subset(weeks, along = 'syear', indices = 1)
+    ## to revise:
+    #if (dim(data_cube$attrs$Dates)['sweek'] > 1) {
+    #  weeks <- Subset(data_cube$attrs$Dates, along = 'sweek',
+    #                  indices = (recipe$Analysis$Workflow$Skill$sweek_window +1)/2)                      
+    #} else {
+    #  weeks <- Subset(data_cube$attrs$Dates, along = 'syear', indices = 1)
+    #}     
+    week_label <- init_week 
+    # week_abbreviation <- week.abb[init_week]
+    } else {
+ 
+    init_month <- 1 
+    init_week <- 1
+  }
+   
   if (!is.null(recipe$Analysis$Workflow$Visualization$projection)) {
     projection <- tolower(recipe$Analysis$Workflow$Visualization$projection)
   } else {
@@ -191,61 +207,73 @@ plot_metrics <- function(recipe, data_cube, metrics,
         } else {
           metric_significance <- NULL
         }
+        
         # Define output file name and titles
         if (tolower(recipe$Analysis$Horizon) == "seasonal") {
           outfile <- paste0(outdir[var], name, "-", month_label)
+        } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+          outfile <- paste0(outdir[var], name, "-", week_label)
         } else {
           outfile <- paste0(outdir[var], name)
         }
         # Get variable name and long name
+
+        
         var_name <- data_cube$attrs$Variable$varName[[var]]
         var_long_name <- data_cube$attrs$Variable$metadata[[var_name]]$long_name
         # Multi-panel or single-panel plots
         if (recipe$Analysis$Workflow$Visualization$multi_panel) {
           # Define titles
-          toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+          if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+            toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
                              "\n", display_name, " / ", hcst_period)
-          ## time_bounds in data_cube to know if Time_aggregation was applied
-          if (!is.null(attributes(data_cube$attrs$time_bounds))) {
-            info(recipe$Run$logger, "Using plotting attrs from time_bounds.")
-            if (length(attributes(data_cube$attrs$time_bounds)$plotting_attr) > 1) {
-              titles <- unlist(
+            ## time_bounds in data_cube to know if Time_aggregation was applied
+            if (!is.null(attributes(data_cube$attrs$time_bounds))) {
+              info(recipe$Run$logger, "Using plotting attrs from time_bounds.")
+              if (length(attributes(data_cube$attrs$time_bounds)$plotting_attr) > 1) {
+                titles <- unlist(
                         lapply(1:length(attributes(data_cube$attrs$time_bounds)$plotting_attr$ini_ftime),
                                function(x) {
                                paste("Forecast time",
                                      attributes(data_cube$attrs$time_bounds)$plotting_attr$ini_ftime[x],
                                      "to",
                                      attributes(data_cube$attrs$time_bounds)$plotting_attr$end_ftime[x])}))
+              } else {
+                titles <- attributes(data_cube$attrs$time_bounds)$plotting_attr[[1]]
+              }
             } else {
-              titles <- attributes(data_cube$attrs$time_bounds)$plotting_attr[[1]]
+              titles <- as.vector(months)
             }
-          } else {
-            titles <- as.vector(months)
+          } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+              toptitle <- paste0(system_name, " / ", str_to_title(var_long_name),
+                             "\n", display_name, " / ", "Issued on ",
+                             format(as.Date(recipe$Analysis$Time$sdate, "%Y%m%d"), "%Y-%m-%d"), 
+                             " / ", hcst_period)
+              monday <- ymd_hms(weeks) - days(wday(ymd_hms(weeks), week_start = 1) - 1) # monday of the week
+              sunday <- monday + days(6) # sunday of the week               
+              titles <- paste("Valid from", format(monday,"%d-%m"), "to", format(sunday, "%d-%m"))
+	  } else {
+            toptitle <- "Unknown"
+            titles <- "Unknown"
           }
           ## TODO: Combine PlotLayout with PlotRobinson?
-          suppressWarnings(
-          PlotLayout(PlotEquiMap, c('longitude', 'latitude'),
-                     asplit(metric, MARGIN=1), # Splitting array into a list
-                     longitude, latitude,
-                     special_args = metric_significance, 
-                     dot_symbol = 20,
-                     toptitle = toptitle,
-                     title_scale = 0.6,
-                     titles = titles,
-                     filled.continents = F,
-                     brks = brks,
-                     cols = cols,
-                     col_inf = col_inf,
-                     col_sup = col_sup,
-                     fileout = paste0(outfile, ".pdf"),
-                     bar_label_digits = 3,
-                     bar_extra_margin = rep(0.9, 4),
-                     extra_margin = rep(1, 4),
-                     bar_label_scale = 1.5,
-                     axes_label_scale = 1.3,
-                     width = 11,#default i
-                     height = 11)
-          )
+          output_configuration <- output_conf$Multipanel$plot_metrics
+          base_args <- list(fun = "PlotEquiMap",
+                            plot_dims = c('longitude', 'latitude'),
+                            var = asplit(metric, MARGIN = 1),
+                            lon = longitude, lat = latitude,
+                            special_args = metric_significance,
+                            dot_symbol = 20, toptitle = toptitle,
+                            title_scale = 0.6, titles = titles,
+                            filled.continents = FALSE, brks = brks,
+                            cols = cols, col_inf = col_inf, col_sup = col_sup,
+                            fileout = paste0(outfile, ".pdf"),
+                            bar_label_digits = 3,
+                            bar_extra_margin = rep(0.9, 4),
+                            extra_margin = rep(1, 4), bar_label_scale = 1.5,
+                            axes_label_scale = 1.3, width = 11, height = 11)
+          base_args[names(output_configuration)] <- output_configuration
+          do.call(PlotLayout, base_args)
         } else {
           # Define function and parameters depending on projection
           if (projection == 'cylindrical_equidistant') {
@@ -281,45 +309,62 @@ plot_metrics <- function(recipe, data_cube, metrics,
           # Loop over forecast times
           for (i in 1:dim(metric)[['time']]) {
             # Get forecast time label
-            # Case without time aggregation:
-            if (is.null(attributes(data_cube$attrs$time_bounds))) {
-              forecast_time <- match(months[i], month.name) - init_month + 1
-                        
-              if (forecast_time < 1) {
-                forecast_time <- forecast_time + 12
-              }
-              forecast_time <- sprintf("%02d", forecast_time)
-              # Define plot title
-              toptitle <- paste(system_name, "/",
-                                str_to_title(var_long_name),
-                                "\n", display_name, "/", months[i], "/",
-                                hcst_period)
-            } else {
-              if (length(attributes(data_cube$attrs$time_bounds)$plotting_attr) > 1) {
-                forecast_time_ini <-attributes(data_cube$attrs$time_bounds)$plotting_attr$ini[i]
-                forecast_time_end <- attributes(data_cube$attrs$time_bounds)$plotting_attr$end[i]
-                # labels for file name:
-                forecast_time <- paste0(forecast_time_ini, "-", forecast_time_end)
-                # title names:
-                forecast_time_ini <- init_month + forecast_time_ini - 1
-                forecat_time_ini <- ifelse(forecast_time_ini > 12, forecast_time_ini - 12, forecast_time_ini) 
-                forecast_time_ini <- month.name[forecast_time_ini]
-                forecast_time_end <- init_month + forecast_time_end - 1
-                forecat_time_end <- ifelse(forecast_time_end > 12, forecast_time_end - 12, forecast_time_end)
-                forecast_time_end <- month.name[forecast_time_end]
+            if (recipe$Analysis$Horizon == "Seasonal") {
+              # Case without time aggregation:
+              if (is.null(attributes(data_cube$attrs$time_bounds))) {
+                forecast_time <- match(months[i], month.name) - init_month + 1                  
+                if (forecast_time < 1) {
+                  forecast_time <- forecast_time + 12
+                }
+                forecast_time <- sprintf("%02d", forecast_time)
+                # Define plot title
                 toptitle <- paste(system_name, "/",
                                   str_to_title(var_long_name),
-                                  "\n", display_name, "/", forecast_time_ini, "to", 
-                                  forecast_time_end, "/",
+                                  "\n", display_name, "/", months[i], "/",
                                   hcst_period)
               } else {
-                forecast_time <- attributes(data_cube$attrs$time_bounds)$plotting_attr[[1]][i]
-                toptitle <- paste(system_name, "/",
-                                  str_to_title(var_long_name),
-                                  "\n", display_name, "/", 
-                                  forecast_time, "/",
-                                  hcst_period)
+                if (length(attributes(data_cube$attrs$time_bounds)$plotting_attr) > 1) {
+                  forecast_time_ini <-attributes(data_cube$attrs$time_bounds)$plotting_attr$ini[i]
+                  forecast_time_end <- attributes(data_cube$attrs$time_bounds)$plotting_attr$end[i]
+                  # labels for file name:
+                  forecast_time <- paste0(forecast_time_ini, "-", forecast_time_end)
+                  # title names:
+                  forecast_time_ini <- init_month + forecast_time_ini - 1
+                  forecat_time_ini <- ifelse(forecast_time_ini > 12, forecast_time_ini - 12, forecast_time_ini) 
+                  forecast_time_ini <- month.name[forecast_time_ini]
+                  forecast_time_end <- init_month + forecast_time_end - 1
+                  forecat_time_end <- ifelse(forecast_time_end > 12, forecast_time_end - 12, forecast_time_end)
+                  forecast_time_end <- month.name[forecast_time_end]
+                  toptitle <- paste(system_name, "/",
+                                    str_to_title(var_long_name),
+                                    "\n", display_name, "/", forecast_time_ini, "to", 
+                                    forecast_time_end, "/",
+                                    hcst_period)
+                } else {
+                  forecast_time <- attributes(data_cube$attrs$time_bounds)$plotting_attr[[1]][i]
+                  toptitle <- paste(system_name, "/",
+                                    str_to_title(var_long_name),
+                                    "\n", display_name, "/", 
+                                    forecast_time, "/",
+                                    hcst_period)
+                }
               }
+            } else if (tolower(recipe$Analysis$Horizon == "subseasonal")) {
+              # Get forecast time label
+              forecast_time <- weeks[i]
+              monday <- ymd_hms(weeks[i]) - days(wday(ymd_hms(weeks[i]), week_start = 1) - 1)
+              sunday <- monday + days(6)
+              toptitle <- paste(system_name, "/",
+                                str_to_title(var_long_name),
+                                "\n", display_name, " / ", "Issued on ",
+                                format(as.Date(recipe$Analysis$Time$sdate, "%Y%m%d"), "%Y-%m-%d"),
+                                " / ", hcst_period, "\n",
+                                paste("Valid from", format(monday,"%d-%m"), "to", 
+                                      format(sunday, "%d-%m"), "of", year(monday)))
+                                # "/ valid week", format(weeks[i],
+                                #  "%Y-%m-%d"), "/", hcst_period)
+            } else {
+              warning("Plotting decadal?")
             }
             # Modify base arguments
             base_args[[1]] <- metric[i, , ]
@@ -332,12 +377,15 @@ plot_metrics <- function(recipe, data_cube, metrics,
             if (identical(fun, PlotRobinson)) {
               ## TODO: Customize alpha and other technical details depending on the metric
               base_args[['caption']] <- 
-                paste0("Nominal start date: ",
-                       "1st of ", str_to_title(month_label), "\n",
-                       "Forecast month: ", forecast_time, "\n",
+                paste0("Nominal start date: ", ymd(start_date), "\n",
+                       "Forecast week: ",  sprintf("%02d", i), "\n", ## This is specific for subseasonal, would need a loop to specify time horizon
                        "Reference: ", recipe$Analysis$Datasets$Reference, "\n",
+                       "Units: ", data_cube$attrs$Variable$metadata[[var_name]]$units, "\n",
                        significance_caption)
             }
+            if (tolower(recipe$Analysis$Horizon == 'subseasonal')) {
+	     forecast_time <- format(forecast_time, "%Y%m%d")
+            } 
             fileout <- paste0(outfile, "_ft", forecast_time, ".pdf")
             # Plot
             info(recipe$Run$logger,
@@ -352,6 +400,6 @@ plot_metrics <- function(recipe, data_cube, metrics,
       }
     }
   }
-  info(recipe$Run$logger,
+  info(recipe$Run$logger, 
        "##### SKILL METRIC PLOTS SAVED TO OUTPUT DIRECTORY #####")
 }
diff --git a/modules/Visualization/R/plot_most_likely_terciles_map.R b/modules/Visualization/R/plot_most_likely_terciles_map.R
index 739f2d10535191d39175abe11c78016526e1bf6e..094aa7a0de9b93558c9b86c55554f5e8299e4450 100644
--- a/modules/Visualization/R/plot_most_likely_terciles_map.R
+++ b/modules/Visualization/R/plot_most_likely_terciles_map.R
@@ -17,13 +17,13 @@ plot_most_likely_terciles <- function(recipe,
                                       dots, 
                                       outdir,
                                       output_conf) {
-
+  
   ## TODO: Add 'anomaly' to plot title
   # Abort if frequency is daily
   if (recipe$Analysis$Variables$freq %in% c("daily", "daily_mean")) {
     stop("Visualization functions not yet implemented for daily data.")
   }
-
+  
   latitude <- fcst$coords$lat
   longitude <- fcst$coords$lon
   archive <- get_archive(recipe)
@@ -38,9 +38,14 @@ plot_most_likely_terciles <- function(recipe,
   if (tolower(recipe$Analysis$Horizon) == "seasonal") {
     init_month <- as.numeric(substr(recipe$Analysis$Time$sdate,
                                     start = 1, stop = 2))
+    
+  } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+    init_week <- as.numeric(substr(recipe$Analysis$Time$sdate,
+                                   start = 1, stop = 2))
   } else {
     ## TODO: Sort out decadal initial month (is it always January?)
     init_month <- 1
+    init_week <- 1   
   }
   # Retrieve and rearrange probability bins for the forecast
   if (is.null(probabilities$probs_fcst$prob_b33) ||
@@ -49,14 +54,14 @@ plot_most_likely_terciles <- function(recipe,
     stop("The forecast tercile probability bins are not present inside ",
          "'probabilities', the most likely tercile map cannot be plotted.")
   }
-
+  
   probs_fcst <- abind(probabilities$probs_fcst$prob_b33,
                       probabilities$probs_fcst$prob_33_to_66,
                       probabilities$probs_fcst$prob_a66,
                       along = 0)
   names(dim(probs_fcst)) <- c("bin",
                               names(dim(probabilities$probs_fcst$prob_b33)))
-
+  
   ## TODO: Improve this section
   # Drop extra dims, add time dim if missing:
   for (var in 1:fcst$dims[['var']]) {
@@ -140,7 +145,13 @@ plot_most_likely_terciles <- function(recipe,
       years <- lubridate::year(fcst$attrs$Dates[1, 1, which(start_date == i_syear), ])
       if (recipe$Analysis$Workflow$Visualization$multi_panel) {
         ## TODO: Ensure this works for daily and sub-daily cases
-        titles <- as.vector(months)
+        if (tolower(recipe$Analysis$Horizon) == "seasonal") {
+          titles <- as.vector(months)
+        } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+          titles <- as.vector(weeks)
+        } else {
+          titles <- as.vector(months)
+        }
         # Plots
         ## NOTE: PlotLayout() and PlotMostLikelyQuantileMap() are still being worked
         ## on. 
@@ -190,6 +201,8 @@ plot_most_likely_terciles <- function(recipe,
                           return_leg = T,
                           triangle_ends = c(F, T) , width = 10, height = 8)
         base_args[names(output_configuration)] <- output_configuration
+      }
+      if (tolower(recipe$Analysis$Horizon) == "seasonal") {
         for (i in 1:length(months)) {
           # Get forecast time label
           if (is.null(attributes(fcst$attrs$time_bounds))) {
@@ -202,6 +215,70 @@ plot_most_likely_terciles <- function(recipe,
             forecast_time <- months
           }
 
+          # Get mask subset
+          if (!is.null(mask)) {
+            mask_i <- Subset(var_mask, along = 'time',
+                             indices = i, drop = TRUE)
+          } else {
+            mask_i <- NULL
+          }
+          # Get dots subset
+          if (!is.null(dots)) {
+            dots_i <- Subset(var_dots, along = 'time',
+                             indices = i, drop = TRUE)
+          } else {
+            dots_i <- NULL
+          }
+          # Define plot title
+          toptitle <- paste0(system_name, " / ",
+                             str_to_title(var_long_name),
+                             "\n", "Most Likely Tercile / ",
+                             months[i], " ", years[i],
+                             " / Start date: ", 
+                             format(as.Date(i_syear, format="%Y%m%d"), 
+                                    "%d-%m-%Y"))
+          # Plot
+          fileout <- paste0(outfile, "_ft", forecast_time, ".pdf")
+          base_args$probs <- i_var_probs[i, , , ]
+          base_args$mask <- mask_i
+          base_args$dots <- dots_i
+          cb_info <- do.call(PlotMostLikelyQuantileMap,
+                             args = c(base_args,
+                                      list(toptitle = toptitle,
+                                           fileout = fileout)))
+          # Add color bars with 1 range for normal category:
+          for (i_bar in 1:cb_info$nmap) {
+            tmp <- cb_info
+            tmp$brks <- tmp$brks[[i_bar]]
+            tmp$cols <- tmp$cols[[i_bar]]
+            tmp$bar_limits <- tmp$bar_limits[[i_bar]]
+            tmp$col_sup <- tmp$col_sup[[i_bar]]
+            tmp$title <- tmp$bar_titles[i_bar]
+            tmp$bar_titles <- NULL
+            tmp$nmap <- NULL
+            tmp$var_limits <- NULL
+            if (length(cb_info$brks[[i_bar]]) > 2) {
+              # plot colorbar as normal
+              do.call(ColorBar, tmp)
+            } else {
+              # plot a square
+              tmp$brks <- 4 
+              tmp$draw_ticks <- F
+              tmp$box_label <- "> 40"
+              tmp$triangle_ends <- c(F, F) 
+              tmp$draw_separators <- FALSE
+              do.call(ColorBar_onebox, tmp)
+            }
+          }
+        }
+      } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") { 
+        for (i in 1:length(weeks)) {
+          # Get forecast time label
+          forecast_time <- match(weeks[i], week.name) - init_week + 1
+          if (forecast_time < 1) {
+            forecast_time <- forecast_time + 12
+          }
+          forecast_time <- sprintf("%02d", forecast_time)
           # Get mask subset
           if (!is.null(mask)) {
             mask_i <- Subset(var_mask, along = 'time', indices = i, drop = TRUE)
@@ -263,10 +340,10 @@ plot_most_likely_terciles <- function(recipe,
             tmp$nmap <- NULL
             tmp$var_limits <- NULL
             if (length(cb_info$brks[[i_bar]]) > 2) {
-            # plot colorbar as normal
+              # plot colorbar as normal
               do.call(ColorBar, tmp)
             } else {
-            # plot a square
+              # plot a square
               tmp$brks <- 4 
               tmp$draw_ticks <- F
               tmp$box_label <- "> 40"
@@ -275,11 +352,11 @@ plot_most_likely_terciles <- function(recipe,
               do.call(ColorBar_onebox, tmp)
             }
           }
-          dev.off()
         }
       }
     }
+    dev.off()
   }
-  info(recipe$Run$logger,
-       "##### MOST LIKELY TERCILE PLOTS SAVED TO OUTPUT DIRECTORY #####")
+info(recipe$Run$logger,
+     "##### MOST LIKELY TERCILE PLOTS SAVED TO OUTPUT DIRECTORY #####")
 }
diff --git a/modules/Visualization/Visualization.R b/modules/Visualization/Visualization.R
index 3750baead9284492691714ef50d356f7a5faa8ba..52ce506875362ab4e7f04bdd5e96aa3e493a8810 100644
--- a/modules/Visualization/Visualization.R
+++ b/modules/Visualization/Visualization.R
@@ -77,8 +77,9 @@ Visualization <- function(recipe,
   # Plot skill metrics
   if ("skill_metrics" %in% plots) {
     if (!is.null(skill_metrics)) {
-      plot_metrics(recipe, data$hcst, skill_metrics, outdir,
-                   significance, output_conf = output_conf)
+      plot_metrics(recipe = recipe, data_cube = data$hcst, 
+                         metrics = skill_metrics, outdir = outdir,
+                         significance = significance, output_conf = output_conf)
     } else {
       error(recipe$Run$logger,
             paste0("The skill metric plots have been requested, but the ",
@@ -158,7 +159,6 @@ Visualization <- function(recipe,
                    "there is no fcst element in 'data'"))
     }
   }
-
   # Plot Most Likely Terciles
   if ("most_likely_terciles" %in% plots) {
     if ((!is.null(probabilities)) && (!is.null(data$fcst))) {
diff --git a/modules/Visualization/output_size.yml b/modules/Visualization/output_size.yml
index 0cd945be4b2d9cf578df057e2f2dc63c8ff241a2..f5259fb4302683daa2cf472bc1569ccbf898aa92 100644
--- a/modules/Visualization/output_size.yml
+++ b/modules/Visualization/output_size.yml
@@ -10,6 +10,7 @@ region: #units inches
         dot_size: 1.7
         dot_symbol: 4
         font.main: 1
+        colNA: "white"
       forecast_ensemble_mean:
         width: 8.5
         height: 8.5
@@ -19,15 +20,67 @@ region: #units inches
         dot_symbol: 4
         dot_size: 1.7
         font.main: 1
+        colNA: "white"
       most_likely_terciles: 
         width: 8.5
         height: 8.5
         dot_size: 2
         plot_margin: !expr c(0, 4.1, 4.1, 2.1)
+        colNA: "white"
     Multipanel:
+      forecast_ensemble_mean: 
+        width: 8.5
+        height: 8.5
     Robinson:
       skill_metrics: {width: 8, height: 5}
+      colNA: "white"
   NA-EU: #Norht Atlantic European region
+  Iberia: #latmin: 34, latmax: 46, lonmin: -10, lonmax: 5
+    PlotEquiMap:
+      skill_metrics:
+        width: 8
+        height: 7.5
+        intylat: 2
+        intxlon: 2 
+        axes_label_scale: 0.8
+        bar_label_scale: 1.2
+        bar_extra_margin: !expr c(2,1,0.5,1)
+        dot_size: 1.7
+        dot_symbol: 4
+        font.main: 1
+        colNA: "white"
+      forecast_ensemble_mean:
+        width: 8
+        height: 7.5
+        intylat: 2
+        intxlon: 2
+        axes_label_scale: 0.8
+        bar_label_scale: 1.2
+        bar_extra_margin: !expr c(2,1,0.5,1)
+        dot_symbol: 4
+        dot_size: 1.7
+        font.main: 1
+        colNA: "white"
+      most_likely_terciles:
+        width: 8
+        height: 7.5       
+        intylat: 2
+        intxlon: 2
+        dot_size: 2
+        plot_margin: !expr c(0, 4.1, 4.1, 2.1)
+        colNA: "white"
+    Multipanel:
+      forecast_ensemble_mean: 
+        width: 8.5
+        height: 9.5
+        title_margin_scale: 4
+      skill_metrics:
+        width: 8.5
+        height: 9.5
+        title_margin_scale: 4
+    Robinson:
+      skill_metrics: {width: 8, height: 5}
+      colNA: "white"
   Mediterranean:
   Global:
 # Add other regions
diff --git a/tools/get_archive.R b/tools/get_archive.R
index 11602698b6a07ad4d82acde880291906263fc201..56309ccb1db11ca6018116b382bd52131e30aacc 100644
--- a/tools/get_archive.R
+++ b/tools/get_archive.R
@@ -2,7 +2,10 @@ get_archive <- function(recipe) {
   if (tolower(recipe$Analysis$Horizon) == "seasonal") {
     archive <- 
       read_yaml(paste0("conf/archive.yml"))[[recipe$Run$filesystem]]
-  } else if (tolower(recipe$Analysis$Horizon) == "decadal") {
+  } else if (tolower(recipe$Analysis$Horizon) == "subseasonal") {
+    archive <-
+      read_yaml(paste0("conf/archive_subseasonal.yml"))[[recipe$Run$filesystem]]
+  } else {
     archive <- 
       read_yaml(paste0("conf/archive_decadal.yml"))[[recipe$Run$filesystem]]
   }