diff --git a/modules/Loading/R/dates2load.R b/modules/Loading/R/dates2load.R
index 16e24647338c60eaf93d62ef41b58b34b84e387f..66aba97b685d6502f515635b1d9cb39fa43089c1 100644
--- a/modules/Loading/R/dates2load.R
+++ b/modules/Loading/R/dates2load.R
@@ -28,6 +28,7 @@ dates2load <- function(recipe, logger) {
   # if (nchar(recipe$sdate) == 8 &&) {
   #   recipe$sdate <- substr(recipe$sdate, 5, 8)
   # }
+  # Case 1: Seasonal monthly
   if (temp_freq == "monthly_mean") {
     # hcst dates
     file_dates <- paste0(strtoi(recipe$hcst_start):strtoi(recipe$hcst_end),
@@ -42,7 +43,9 @@ dates2load <- function(recipe, logger) {
            paste("January 1993 start date is not available for", system_name,
                  "and has been removed from the list of start dates."))
     }
-    file_dates <- .add_dims(file_dates)
+    
+    dim(file_dates) <- c(syear = length(file_dates))
+  # Case 2: subseasonal weekly
   } else if (temp_freq == "weekly_mean") {
     sday <- recipe$sday_window
     if (is.null(sday)) {
@@ -73,16 +76,19 @@ dates2load <- function(recipe, logger) {
                                      hcst.end = as.numeric(recipe$hcst_end),
                                      ftime_min = recipe$ftime_min,
                                      ftime_max = recipe$ftime_max, out = 'hcst')
+  # Case 3: daily
+  ## TODO: distinguish between seasonal daily and subseasonal daily?
   } else {
     file_dates <- paste0(strtoi(recipe$hcst_start):strtoi(recipe$hcst_end),
                          recipe$sdate)
-    file_dates <- .add_dims(file_dates)
+    dim(file_dates) <- c(syear = length(file_dates))
   }
+
   # fcst dates (if fcst_year empty it creates an empty object)
   if (!is.null(recipe$fcst_year)) {
     if (temp_freq == "monthly_mean" || temp_freq == "daily_mean") {
       file_dates.fcst <- paste0(recipe$fcst_year, recipe$sdate)
-      file_dates.fcst <- .add_dims(file_dates.fcst)
+      dim(file_dates.fcst) <- c(syear = length(file_dates.fcst))    
     } else if (temp_freq == "weekly_mean") {
       file_dates.fcst <- subseas_file_dates(startdate = fcst.sdate_to_start,
                                             n.skill.weeks = n.skill.weeks,
diff --git a/modules/Loading/R/get_forecast_times.R b/modules/Loading/R/get_forecast_times.R
index 621c7d2778596255a3dc7b5f883dcddf20fbccd3..e5cd4fe454fc539762ef3a9409af3eacf065c276 100644
--- a/modules/Loading/R/get_forecast_times.R
+++ b/modules/Loading/R/get_forecast_times.R
@@ -41,7 +41,9 @@ get_forecast_times <- function(sdates, ltmin, ltmax,
                         tz = 'UTC', origin = '1970-01-01')
     lubridate::hour(indxs) <- 12 
     lubridate::minute(indxs) <- 01
-    dim(indxs) <- c(file_date = length(sdates),
+    dim(indxs) <- c(sday = 1,
+                    sweek = 1,
+                    syear = length(sdates),
                     time = time_length)
 
   } else if (time_freq == "monthly_mean") {
@@ -59,13 +61,16 @@ get_forecast_times <- function(sdates, ltmin, ltmax,
                    idx_max = ltmax - 1,
                    ncores = ncores
                    )[[1]]
+    indxs <- s2dv::Reorder(indxs, c("syear", "time"))
     indxs_dims <- dim(indxs)
     indxs <- as.POSIXct(indxs, origin = "1970-01-01", tz = "UTC")
     lubridate::day(indxs) <- round(days_in_month(indxs)/2 + 1)
     lubridate::hour(indxs) <- 12
     lubridate::minute(indxs) <- 01
-    dim(indxs) <- indxs_dims
-    names(dim(indxs)) <- c("time", "file_date")
+    dim(indxs) <- c(sday = 1,
+                    sweek = 1,
+                    syear = indxs_dims[['syear']],
+                    time = indxs_dims[['time']])
   }
 
   # TODO: 6 hourly case
diff --git a/modules/Loading/R/load_seasonal.R b/modules/Loading/R/load_seasonal.R
index 1824b7d5b9408bca1be3b2374e702ca6a3cf55ba..dd40ecec5b854b426f9e7eccc4ba021a708bf106 100644
--- a/modules/Loading/R/load_seasonal.R
+++ b/modules/Loading/R/load_seasonal.R
@@ -69,13 +69,13 @@ load_seasonal <- function(recipe) {
   
   # -----------
   obs.path <- paste0(archive$src_ref, obs.dir, "$var_dir$",
-                     "$var$_$file_date$.nc")
+                     "$var$_$syear$.nc")
   
   hcst.path <- paste0(archive$src_sys, hcst.dir, "$var_dir$",
-                      "$var$_$file_date$.nc")
+                      "$var$_$syear$.nc")
   
   fcst.path <- paste0(archive$src_sys, hcst.dir, "$var_dir$",
-                      "$var$_$file_date$.nc")
+                      "$var$_$syear$.nc")
   
   # Define regrid parameters:
   #-------------------------------------------------------------------
@@ -87,13 +87,13 @@ load_seasonal <- function(recipe) {
   
   # Load hindcast 
   #-------------------------------------------------------------------
-
-  dim(sdates$hcst) <- dim(sdates$hcst)['syear']
+  
   hcst <- Start(dat = hcst.path,
                 var = variable,
                 var_dir = var_dir_exp,
-                file_date = sdates$hcst,
+                syear = sdates$hcst,
                 time = idxs$hcst,
+                time_across = 'syear',
                 var_dir_depends = 'var',
                 latitude = values(list(lats.min, lats.max)),
                 latitude_reorder = Sort(),
@@ -108,37 +108,19 @@ load_seasonal <- function(recipe) {
                                 longitude = c('lon', 'longitude'),
                                 ensemble = c('member', 'ensemble', 'lev')),
                 ensemble = indices(1:hcst.nmember),
-                metadata_dims = 'var', # change to just 'var'?
+                metadata_dims = 'var',
                 return_vars = list(latitude = 'dat',
                                    longitude = 'dat',
-                                   time = 'file_date'),
-                split_multiselected_dims = FALSE,
+                                   time = 'syear'),
+                split_multiselected_dims = TRUE,
+                merge_across_dims = TRUE,
                 retrieve = TRUE)
- 
+
   # Remove var_dir dimension
   if ("var_dir" %in% names(dim(hcst))) {
     hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected")
   }
   
-  # Adjusts dims for daily case, could be removed if startR allows
-  # multidim split
-  names(dim(hcst))[which(names(dim(hcst)) == 'file_date')] <- "syear"
-  default_dims <- c(dat = 1, var = 1, sday = 1, 
-                    sweek = 1, syear = 1, time = 1, 
-                    latitude = 1, longitude = 1, ensemble = 1)
-  default_dims[names(dim(hcst))] <- dim(hcst)
-  dim(hcst) <- default_dims
-  # Change time attribute dimensions
-  default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-  names(dim(attr(hcst, "Variables")$common$time))[which(names(
-    dim(attr(hcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-  default_time_dims[names(dim(attr(hcst, "Variables")$common$time))] <- 
-    dim(attr(hcst, "Variables")$common$time)
-  dim(attr(hcst, "Variables")$common$time) <- default_time_dims
-  
-  # Convert hcst to s2dv_cube object
-  ## TODO: Give correct dimensions to $Dates
-  ## (sday, sweek, syear instead of file_date)
   hcst <- as.s2dv_cube(hcst)
   # Adjust dates for models where the time stamp goes into the next month
   if (recipe$Analysis$Variables$freq == "monthly_mean") {
@@ -151,13 +133,14 @@ load_seasonal <- function(recipe) {
     # the call uses file_date instead of fcst_syear so that it can work
     # with the daily case and the current version of startR not allowing
     # multiple dims split
-    
+ 
     fcst <- Start(dat = fcst.path,
                   var = variable,
                   var_dir = var_dir_exp,
-                  var_dir_depends = 'var',
-                  file_date = sdates$fcst,
+                  syear = sdates$fcst,
                   time = idxs$fcst,
+                  time_across = 'syear',
+                  var_dir_depends = 'var',
                   latitude = values(list(lats.min, lats.max)),
                   latitude_reorder = Sort(),
                   longitude = values(list(lons.min, lons.max)),
@@ -174,30 +157,15 @@ load_seasonal <- function(recipe) {
                   metadata_dims = 'var',
                   return_vars = list(latitude = 'dat',
                                      longitude = 'dat',
-                                     time = 'file_date'),
-                  split_multiselected_dims = FALSE,
+                                     time = 'syear'),
+                  split_multiselected_dims = TRUE,
+                  merge_across_dims = TRUE,
                   retrieve = TRUE)
-    
+
     if ("var_dir" %in% names(dim(fcst))) {
       fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected")
     }
     
-    # Adjusts dims for daily case, could be removed if startR allows
-    # multidim split
-    names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear"
-    default_dims <- c(dat = 1, var = 1, sday = 1,
-                      sweek = 1, syear = 1, time = 1,
-                      latitude = 1, longitude = 1, ensemble = 1)
-    default_dims[names(dim(fcst))] <- dim(fcst)
-    dim(fcst) <- default_dims
-    # Change time attribute dimensions
-    default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-    names(dim(attr(fcst, "Variables")$common$time))[which(names(
-      dim(attr(fcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-    default_time_dims[names(dim(attr(fcst, "Variables")$common$time))] <- 
-      dim(attr(fcst, "Variables")$common$time)
-    dim(attr(fcst, "Variables")$common$time) <- default_time_dims
-    
     # Convert fcst to s2dv_cube
     fcst <- as.s2dv_cube(fcst)
     # Adjust dates for models where the time stamp goes into the next month
@@ -223,13 +191,12 @@ load_seasonal <- function(recipe) {
     
     dates_file <- format(as.Date(idxs$hcst, '%Y%m%d'), "%Y%m")
     dim(dates_file) <- dim(idxs$hcst)
-    dates_file <- Reorder(dates_file, order = c("file_date", "time"))
-    
+    # dates_file dims: [sday, sweek, syear, time]
     obs <- Start(dat = obs.path,
                  var = variable,
                  var_dir = var_dir_obs,
                  var_dir_depends = 'var',
-                 file_date = dates_file,
+                 syear = dates_file,
                  latitude = values(list(lats.min, lats.max)),
                  latitude_reorder = Sort(),
                  longitude = values(list(lons.min, lons.max)),
@@ -244,19 +211,16 @@ load_seasonal <- function(recipe) {
                  metadata_dims = 'var',
                  return_vars = list(latitude = 'dat',
                                     longitude = 'dat',
-                                    time = 'file_date'),
+                                    time = 'syear'),
                  split_multiselected_dims = TRUE,
                  retrieve = TRUE)
     
   } else if (store.freq %in% c("daily_mean", "daily")) {
-    
+
+    # $syear$ must be replaced in obs.path due to startR issue
+    obs.path <- gsub("$syear$", "$file_date$", obs.path, fixed = TRUE)
     # Get year and month for file_date
-    dates_file <- sapply(dates, format, '%Y%m')
-    dim(dates_file) <- dim(dates)
-    # Set hour to 12:00 to ensure correct date retrieval for daily data
-    lubridate::hour(dates) <- 12 
-    lubridate::minute(dates) <- 00
-    dim(idxs$hcst) <- dim(dates_file)
+    dates_file <- sapply(idxs$hcst, format, '%Y%m')
 
     obs <- Start(dat = obs.path,
                  var = variable,
@@ -291,19 +255,13 @@ load_seasonal <- function(recipe) {
   if ("var_dir" %in% names(dim(obs))) {
     obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected")
   }
-  # Adds ensemble dim to obs (for consistency with hcst/fcst)
-  names(dim(obs))[which(names(dim(obs)) == 'file_date')] <- "syear"
-  default_dims <- c(dat = 1, var = 1, sday = 1, 
-                    sweek = 1, syear = 1, time = 1, 
+  # Add ensemble dim to obs (for consistency with hcst/fcst)
+  default_dims <- c(dat = 1, var = 1, sday = 1,
+                    sweek = 1, syear = 1, time = 1,
                     latitude = 1, longitude = 1, ensemble = 1)
   default_dims[names(dim(obs))] <- dim(obs)
   dim(obs) <- default_dims
 
-  names(dim(attr(obs, "Variables")$common$time))[which(names(
-    dim(attr(obs, "Variables")$common$time)) == 'file_date')] <- "syear"
-  default_time_dims[names(dim(attr(obs, "Variables")$common$time))] <- 
-    dim(attr(obs, "Variables")$common$time)
-  dim(attr(obs, "Variables")$common$time) <- default_time_dims
   # Convert obs to s2dv_cube
   obs <- as.s2dv_cube(obs)
   
diff --git a/modules/Loading/R/load_tas_tos.R b/modules/Loading/R/load_tas_tos.R
index c19b79923bdc5a676b5b071bcfbd8cb095bb5a25..a8c68986194fefc44fc55bc0f1f6fd5eda70d3eb 100644
--- a/modules/Loading/R/load_tas_tos.R
+++ b/modules/Loading/R/load_tas_tos.R
@@ -71,13 +71,13 @@ load_tas_tos <- function(recipe) {
   
   # -----------
   obs.path <- paste0(archive$src_ref, obs.dir, "$var_dir$",
-                     "$var$_$file_date$.nc")
+                     "$var$_$syear$.nc")
   
   hcst.path <- paste0(archive$src_sys, hcst.dir, "$var_dir$",
-                      "$var$_$file_date$.nc")
+                      "$var$_$syear$.nc")
   
   fcst.path <- paste0(archive$src_sys, hcst.dir, "$var_dir$",
-                      "/$var$_$file_date$.nc")
+                      "/$var$_$syear$.nc")
   
   # Define regrid parameters:
   #-------------------------------------------------------------------
@@ -89,12 +89,13 @@ load_tas_tos <- function(recipe) {
   
   # Load hindcast data without regrid
   #-------------------------------------------------------------------
-  dim(sdates$hcst) <- dim(sdates$hcst)['syear']
+
   hcst <- Start(dat = hcst.path,
                 var = variable,
                 var_dir = var_dir_exp,
-                file_date = sdates$hcst,
+                syear = sdates$hcst,
                 time = idxs$hcst,
+                time_across = 'syear',
                 var_dir_depends = 'var',
                 latitude = values(list(lats.min, lats.max)),
                 latitude_reorder = Sort(decreasing = TRUE),
@@ -104,43 +105,25 @@ load_tas_tos <- function(recipe) {
                                 longitude = c('lon', 'longitude'),
                                 ensemble = c('member', 'ensemble', 'lev')),
                 ensemble = indices(1:hcst.nmember),
-                metadata_dims = 'var', # change to just 'var'?
+                metadata_dims = 'var',
                 return_vars = list(latitude = 'dat',
                                    longitude = 'dat',
-                                   time = 'file_date'),
-                split_multiselected_dims = FALSE,
+                                   time = 'syear'),
+                split_multiselected_dims = TRUE,
+                merge_across_dims = TRUE,
                 retrieve = TRUE)
-  
 
   # Remove var_dir dimension
   if ("var_dir" %in% names(dim(hcst))) {
     hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected")
   }
   
-  # Adjusts dims for daily case, could be removed if startR allows
-  # multidim split
-  names(dim(hcst))[which(names(dim(hcst)) == 'file_date')] <- "syear"
-  default_dims <- c(dat = 1, var = 1, sday = 1, 
-                    sweek = 1, syear = 1, time = 1, 
-                    latitude = 1, longitude = 1, ensemble = 1)
-  default_dims[names(dim(hcst))] <- dim(hcst)
-  browser()
-  dim(hcst) <- default_dims
-  # Change time attribute dimensions
-  default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-  names(dim(attr(hcst, "Variables")$common$time))[which(names(
-    dim(attr(hcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-  default_time_dims[names(dim(attr(hcst, "Variables")$common$time))] <- 
-    dim(attr(hcst, "Variables")$common$time)
-  dim(attr(hcst, "Variables")$common$time) <- default_time_dims
-  
   # Define sea-ice grid points based of sea-ice concentration threshold
   ice_hcst <- hcst[,3,,,,,,,] >= sic.threshold
   
   # Replace Tos with Tas for data points with sea ice
   hcst[,2,,,,,,,][ice_hcst] <- hcst[,1,,,,,,,][ice_hcst] 
   
-  
   # Convert hcst to s2dv_cube object
   ## TODO: Give correct dimensions to $Dates
   ## (sday, sweek, syear instead of file_date)
@@ -166,16 +149,14 @@ load_tas_tos <- function(recipe) {
   # Load forecast data without regrid
   #-------------------------------------------------------------------
   if (!is.null(recipe$Analysis$Time$fcst_year)) {
-    # the call uses file_date instead of fcst_syear so that it can work
-    # with the daily case and the current version of startR not allowing
-    # multiple dims split
-    
+
     fcst <- Start(dat = fcst.path,
                   var = variable,
                   var_dir = var_dir_exp,
-                  var_dir_depends = 'var',
-                  file_date = sdates$fcst,
+                  syear = sdates$fcst,
                   time = idxs$fcst,
+                  time_across = 'syear',
+                  var_dir_depends = 'var',
                   latitude = values(list(lats.min, lats.max)),
                   latitude_reorder = Sort(decreasing = TRUE),
                   longitude = values(list(lons.min, lons.max)),
@@ -187,31 +168,15 @@ load_tas_tos <- function(recipe) {
                   metadata_dims = 'var',
                   return_vars = list(latitude = 'dat',
                                      longitude = 'dat',
-                                     time = 'file_date'),
-                  split_multiselected_dims = FALSE,
+                                     time = 'syear'),
+                  split_multiselected_dims = TRUE,
+                  merge_across_dims = TRUE,
                   retrieve = TRUE)
- 
+
     if ("var_dir" %in% names(dim(fcst))) {
       fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected")
     }
     
-    # Adjusts dims for daily case, could be removed if startR allows
-    # multidim split
-    names(dim(fcst))[which(names(dim(fcst)) == 'file_date')] <- "syear"
-    default_dims <- c(dat = 1, var = 1, sday = 1,
-                      sweek = 1, syear = 1, time = 1,
-                      latitude = 1, longitude = 1, ensemble = 1)
-    default_dims[names(dim(fcst))] <- dim(fcst)
-    dim(fcst) <- default_dims
-    # Change time attribute dimensions
-    default_time_dims <- c(sday = 1, sweek = 1, syear = 1, time = 1)
-    names(dim(attr(fcst, "Variables")$common$time))[which(names(
-      dim(attr(fcst, "Variables")$common$time)) == 'file_date')] <- "syear"
-    default_time_dims[names(dim(attr(fcst, "Variables")$common$time))] <- 
-      dim(attr(fcst, "Variables")$common$time)
-    dim(attr(fcst, "Variables")$common$time) <- default_time_dims
-    
-    
     # Define sea-ice grid points based of sea-ice concentration threshold
     ice_fcst <- fcst[,3,,,,,,,] >= sic.threshold 
     
@@ -254,53 +219,42 @@ load_tas_tos <- function(recipe) {
   
   # Separate Start() call for monthly vs daily data
   if (frequency == "monthly_mean") {
-    
+   
     dates_file <- format(as.Date(idxs$hcst, '%Y%m%d'), "%Y%m")
     dim(dates_file) <- dim(idxs$hcst)
-    dates_file <- Reorder(dates_file, order = c("file_date", "time"))
-    
-    # Define variables for blended tas-tos datasets
-    if (recipe$Analysis$Datasets$Reference$name == 'BEST'){
-      variable <- 'tas' 
-      var_dir_obs <- reference_descrip[[frequency]][variable]
-    }
-    
+    # dates_file dims: [sday, sweek, syear, time]
     obs <- Start(dat = obs.path,
                  var = variable,
                  var_dir = var_dir_obs,
                  var_dir_depends = 'var',
-                 file_date = dates_file,
+                 syear = dates_file,
                  latitude = values(list(lats.min, lats.max)),
                  latitude_reorder = Sort(decreasing = TRUE),
                  longitude = values(list(lons.min, lons.max)),
                  longitude_reorder = circularsort,
+                 transform_vars = c('latitude', 'longitude'),
                  synonims = list(latitude = c('lat','latitude'),
                                  longitude = c('lon','longitude')),
                  metadata_dims = 'var',
                  return_vars = list(latitude = 'dat',
                                     longitude = 'dat',
-                                    time = 'file_date'),
+                                    time = 'syear'),
                  split_multiselected_dims = TRUE,
                  retrieve = TRUE)
-    
-    
+
   } else if (frequency == "daily_mean") {
-    
+
+    # $syear$ must be replaced in obs.path due to startR issue
+    obs.path <- gsub("$syear$", "$file_date$", obs.path, fixed = TRUE)
     # Get year and month for file_date
-    dates_file <- sapply(dates, format, '%Y%m')
-    dim(dates_file) <- dim(dates)
-    # Set hour to 12:00 to ensure correct date retrieval for daily data
-    lubridate::hour(dates) <- 12 
-    lubridate::minute(dates) <- 00
-    # Restore correct dimensions
-    dim(dates) <- dim(dates_file)
-    
+    dates_file <- sapply(idxs$hcst, format, '%Y%m')
+
     obs <- Start(dat = obs.path,
                  var = variable,
                  var_dir = var_dir_obs,
                  var_dir_depends = 'var',
                  file_date = sort(unique(dates_file)),
-                 time = dates,
+                 time = idxs$hcst,
                  time_var = 'time',
                  time_across = 'file_date',
                  merge_across_dims = TRUE,
@@ -308,7 +262,6 @@ load_tas_tos <- function(recipe) {
                  latitude = values(list(lats.min, lats.max)),
                  latitude_reorder = Sort(decreasing = TRUE),
                  longitude = values(list(lons.min, lons.max)),
-                 longitude_reorder = circularsort,
                  synonims = list(latitude = c('lat','latitude'),
                                  longitude = c('lon','longitude')),
                  metadata_dims = 'var',
@@ -317,7 +270,7 @@ load_tas_tos <- function(recipe) {
                                     time = 'file_date'),
                  split_multiselected_dims = TRUE,
                  retrieve = TRUE)
-    
+   
   }
   
   # Remove var_dir dimension
@@ -325,7 +278,6 @@ load_tas_tos <- function(recipe) {
     obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected")
   }
   # Adds ensemble dim to obs (for consistency with hcst/fcst)
-  names(dim(obs))[which(names(dim(obs)) == 'file_date')] <- "syear"
   default_dims <- c(dat = 1, var = 1, sday = 1, 
                     sweek = 1, syear = 1, time = 1, 
                     latitude = 1, longitude = 1, ensemble = 1)
diff --git a/modules/Loading/R/mask_tas_tos.R b/modules/Loading/R/mask_tas_tos.R
index 35f5832c74501cb4692b51d949e403d3a10e4bb0..a358c76b88ad381b87c7780dbde972827d2f90b2 100644
--- a/modules/Loading/R/mask_tas_tos.R
+++ b/modules/Loading/R/mask_tas_tos.R
@@ -1,7 +1,3 @@
-library(multiApply)
-library(startR)
-library(s2dv)
-
 mask_tas_tos <- function(input_data, mask_path, lsm_var_name = lsm_var_name, lon, lat, region = region,
                          lon_dim = 'lon', lat_dim = 'lat', 
                          ncores = NULL){
diff --git a/modules/Saving/R/save_forecast.R b/modules/Saving/R/save_forecast.R
index 807a6bcb2b663922a0872ea41ddefe678182c07e..5198d278f0df69fe40c9f7db754d97a9f0080a39 100644
--- a/modules/Saving/R/save_forecast.R
+++ b/modules/Saving/R/save_forecast.R
@@ -125,7 +125,8 @@ save_forecast <- function(recipe,
         fcst.sdate <- format(fcst.sdate, '%Y%m%d')
   
       } else {
-        fcst.sdate <- data_cube$attrs$load_parameters$dat1$file_date[[1]][i]
+        load_param <- which(names(data_cube$attrs$load_parameters$dat1) %in% c("file_date", "syear"))
+        fcst.sdate <- data_cube$attrs$load_parameters$dat1[[load_param]][[1]][i]
       }
   
       # Get time dimension values and metadata
diff --git a/modules/Saving/R/save_observations.R b/modules/Saving/R/save_observations.R
index b357d5c3cd56c2ae00d41080b1ba94535159cf32..22ba37c2f8164cacddfea7586ba580e660f00d20 100644
--- a/modules/Saving/R/save_observations.R
+++ b/modules/Saving/R/save_observations.R
@@ -112,7 +112,8 @@ save_observations <- function(recipe,
         fcst.sdate <- init_date + lubridate::years(syears_val[i] - lubridate::year(init_date))
       } else {
         if (store.freq == "monthly_mean") {
-          fcst.sdate <- data_cube$attrs$load_parameters$dat1$file_date[[1]][i]
+          load_param <- which(names(data_cube$attrs$load_parameters$dat1) %in% c("file_date", "syear"))
+          fcst.sdate <- data_cube$attrs$load_parameters$dat1[[load_param]][[1]][i]
           fcst.sdate <- as.Date(paste0(fcst.sdate, "01"), '%Y%m%d')
         } else {
           fcst.sdate <- as.Date(data_cube$attrs$Dates[i])
diff --git a/modules/Saving/R/save_probabilities.R b/modules/Saving/R/save_probabilities.R
index 1f6a597a9c624316582d3d950bebe27898b37023..070deef9b24f3b8f4de9947e4084c088564ba2b8 100644
--- a/modules/Saving/R/save_probabilities.R
+++ b/modules/Saving/R/save_probabilities.R
@@ -110,7 +110,8 @@ save_probabilities <- function(recipe,
         fcst.sdate <- init_date + lubridate::years(syears_val[i] - lubridate::year(init_date))
         fcst.sdate <- format(fcst.sdate, '%Y%m%d')
       } else {
-        fcst.sdate <- data_cube$attrs$load_parameters$dat1$file_date[[1]][i]
+        load_param <- which(names(data_cube$attrs$load_parameters$dat1) %in% c("file_date", "syear"))
+        fcst.sdate <- data_cube$attrs$load_parameters$dat1[[load_param]][[1]][i]
       }
 
       # Get time dimension values and metadata