diff --git a/modules/Loading/R/dates2load.R b/modules/Loading/R/dates2load.R
index 7cf80b8f080a097eca647c6d349cdfb15cedf50d..1d3fcbaa4b82d2693deba8c10b76353a276a4034 100644
--- a/modules/Loading/R/dates2load.R
+++ b/modules/Loading/R/dates2load.R
@@ -26,6 +26,7 @@ dates2load <- function(recipe, logger) {
    # hcst dates
     file_dates <- paste0(strtoi(recipe$hcst_start):strtoi(recipe$hcst_end),
                          recipe$sdate)
+    ## TODO: Add dims?
     file_dates <- .add_dims(file_dates)
   } else if (temp_freq == "weekly_mean") {
     sday <- recipe$sday_window
diff --git a/modules/Loading/R/get_timeidx.R b/modules/Loading/R/get_timeidx.R
index dfd23ee083cee7325175990077b2e9ef17436b9c..2796b2021d3ce341508a6ec11ab4799cb2a53188 100644
--- a/modules/Loading/R/get_timeidx.R
+++ b/modules/Loading/R/get_timeidx.R
@@ -35,7 +35,6 @@ get_timeidx <- function(sdates, ltmin, ltmax,
       day_seq <- seq(idx_min[sdate], idx_max[sdate], by = 'days')
       indxs[sdate,] <- day_seq[!(format(day_seq, "%m%d") == "0229")]
     }
-    # browser()
     indxs <- as.POSIXct(indxs*86400, 
                         tz = 'UTC', origin = '1970-01-01')
     lubridate::hour(indxs) <- 12 
@@ -44,10 +43,27 @@ get_timeidx <- function(sdates, ltmin, ltmax,
                     time = time_length)
 
   } else if (time_freq == "monthly_mean") {
-  
-    idx_min <- ltmin
-    idx_max <- ltmax
-    indxs <- indices(idx_min:idx_max)
+    # idx_min <- ltmin
+    # idx_max <- ltmax
+    # indxs <- indices(idx_min:idx_max)
+    indxs <- Apply(sdates,
+                   margins = list(x = "syear"),
+                   fun = function(x, idx_min, idx_max) {
+                     x <- as.POSIXct(x, format = "%Y%m%d", tz = "UTC") %m+% months(idx_min:idx_max)
+                     return(as.array(x))
+                   },
+                   output_dims = c("time"),
+                   idx_min = ltmin - 1,
+                   idx_max = ltmax - 1,
+                   ncores = recipe$Analysis$ncores
+                   )[[1]]
+    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")
   }
 
   # TODO: 6 hourly case
diff --git a/modules/Loading/R/load_seasonal.R b/modules/Loading/R/load_seasonal.R
index 919d54246d41ee082e3b26b9839c5b9cde9f0d55..1e6c9367be006160210168efb9d524ccacccd874 100644
--- a/modules/Loading/R/load_seasonal.R
+++ b/modules/Loading/R/load_seasonal.R
@@ -39,13 +39,6 @@ load_seasonal <- function(recipe) {
                              time_freq=store.freq)
   }
 
-  ## TODO: Examine this verifications part, verify if it's necessary
-  # stream <- verifications$stream
-  # sdates <- verifications$fcst.sdate
-  
-  ## TODO: define fcst.name
-  ##fcst.name <- recipe$Analysis$Datasets$System[[sys]]$name
-  
   # get datasets dict:
   archive <- get_archive(recipe)
   exp_descrip <- archive$System[[exp.name]]
@@ -92,15 +85,10 @@ load_seasonal <- function(recipe) {
   #-------------------------------------------------------------------
   circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max)
   
-  if (recipe$Analysis$Variables$freq == "monthly_mean"){
-    split_multiselected_dims = TRUE
-  } else {
-    split_multiselected_dims = FALSE
-  }
-
   # Load hindcast 
   #-------------------------------------------------------------------
-  
+
+  dim(sdates$hcst) <- dim(sdates$hcst)['syear']
   hcst <- Start(dat = hcst.path,
                 var = variable,
                 var_dir = var_dir_exp,
@@ -124,31 +112,29 @@ load_seasonal <- function(recipe) {
                 return_vars = list(latitude = 'dat',
                                    longitude = 'dat',
                                    time = 'file_date'),
-                split_multiselected_dims = split_multiselected_dims,
+                split_multiselected_dims = FALSE,
                 retrieve = TRUE)
-  
+ 
   # Remove var_dir dimension
   if ("var_dir" %in% names(dim(hcst))) {
     hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected")
   }
   
-  if (store.freq %in% c("daily_mean", "daily")) {
-    # 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
-  }
+  # 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
@@ -189,30 +175,28 @@ load_seasonal <- function(recipe) {
                   return_vars = list(latitude = 'dat',
                                      longitude = 'dat',
                                      time = 'file_date'),
-                  split_multiselected_dims = split_multiselected_dims,
+                  split_multiselected_dims = FALSE,
                   retrieve = TRUE)
     
     if ("var_dir" %in% names(dim(fcst))) {
       fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected")
     }
     
-    if (store.freq %in% c("daily_mean", "daily")) {
-      # 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
-    }
+    # 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)
@@ -237,8 +221,9 @@ load_seasonal <- function(recipe) {
   # Separate Start() call for monthly vs daily data
   if (store.freq == "monthly_mean") {
     
-    dates_file <- format(as.Date(dates, '%Y%m%d'), "%Y%m")
-    dim(dates_file) <- dim(dates)
+    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"))
     
     obs <- Start(dat = obs.path,
                  var = variable,
@@ -302,18 +287,23 @@ load_seasonal <- function(recipe) {
                  retrieve = TRUE)
   }
   
-  
   # Remove var_dir dimension
   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, 
                     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 f866407a707612f6617213fdfcfde67581f268b0..d035f05d04a89afbf6f950a98559227d7293ff9f 100644
--- a/modules/Loading/R/load_tas_tos.R
+++ b/modules/Loading/R/load_tas_tos.R
@@ -87,14 +87,9 @@ load_tas_tos <- function(recipe) {
   #-------------------------------------------------------------------
   circularsort <- check_latlon(lats.min, lats.max, lons.min, lons.max)
   
-  if (recipe$Analysis$Variables$freq == "monthly_mean"){
-    split_multiselected_dims = TRUE
-  } else {
-    split_multiselected_dims = FALSE
-  }
-  
   # Load hindcast data without regrid
   #-------------------------------------------------------------------
+  dim(sdates$hcst) <- dim(sdates$hcst)['syear']
   hcst <- Start(dat = hcst.path,
                 var = variable,
                 var_dir = var_dir_exp,
@@ -113,7 +108,7 @@ load_tas_tos <- function(recipe) {
                 return_vars = list(latitude = 'dat',
                                    longitude = 'dat',
                                    time = 'file_date'),
-                split_multiselected_dims = split_multiselected_dims,
+                split_multiselected_dims = FALSE,
                 retrieve = TRUE)
   
 
@@ -122,23 +117,22 @@ load_tas_tos <- function(recipe) {
     hcst <- Subset(hcst, along = "var_dir", indices = 1, drop = "selected")
   }
   
-  if (recipe$Analysis$Variables$freq == "daily_mean") {
-    # 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
-  }
+  # 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
@@ -194,30 +188,29 @@ load_tas_tos <- function(recipe) {
                   return_vars = list(latitude = 'dat',
                                      longitude = 'dat',
                                      time = 'file_date'),
-                  split_multiselected_dims = split_multiselected_dims,
+                  split_multiselected_dims = FALSE,
                   retrieve = TRUE)
  
     if ("var_dir" %in% names(dim(fcst))) {
       fcst <- Subset(fcst, along = "var_dir", indices = 1, drop = "selected")
     }
     
-    if (recipe$Analysis$Variables$freq == "daily_mean") {
-      # 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
-    }
+    # 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 
@@ -262,9 +255,9 @@ load_tas_tos <- function(recipe) {
   # Separate Start() call for monthly vs daily data
   if (frequency == "monthly_mean") {
     
-    dates_file <- format(as.Date(dates, '%Y%m%d'), "%Y%m")
-    dim(dates_file) <- dim(dates)
-    
+    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'){
@@ -332,20 +325,18 @@ 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)
   default_dims[names(dim(obs))] <- dim(obs)
   dim(obs) <- default_dims
   
-  if(!recipe$Analysis$Datasets$Reference$name %in% c('HadCRUT4','HadCRUT5','BEST','GISTEMPv4')){
-    
+  if (!recipe$Analysis$Datasets$Reference$name %in% c('HadCRUT4', 'HadCRUT5', 'BEST', 'GISTEMPv4')) {
     # Define sea-ice grid points based of sea-ice concentration threshold
-    ice_obs <- (obs[,3,,,,,,,]) >= sic.threshold
-    
+    ice_obs <- (obs[,3,,,,,,,]) >= sic.threshold    
     # Replace NA values with False
     ice_obs[is.na(ice_obs)] <- FALSE
-
     # Replace Tos with Tas for datapoints with sea ice
     obs[,2,,,,,,,][ice_obs] <- obs[,1,,,,,,,][ice_obs] 
   }
@@ -376,7 +367,7 @@ load_tas_tos <- function(recipe) {
   #------------------------------------------------------------------- 
   
   # Regrid reference to system grid:
-  if(recipe$Analysis$Regrid$type == 'to_system'){
+  if (recipe$Analysis$Regrid$type == 'to_system') {
     
     aux <- CDORemap(data_array = obs$data, ## Not regridding to desired grid when latitudes are ordered descending 
                     lons = obs$coords$longitude,
@@ -397,7 +388,7 @@ load_tas_tos <- function(recipe) {
   }
   
   # Regrid system to reference grid:
-  if(recipe$Analysis$Regrid$type == 'to_reference'){
+  if (recipe$Analysis$Regrid$type == 'to_reference') {
     
     aux <- CDORemap(data_array = hcst$data, 
                     lons = hcst$coords$longitude, lats = hcst$coords$latitude, 
@@ -433,7 +424,7 @@ load_tas_tos <- function(recipe) {
   }
   
   # Regrid all data to user defined grid:
-  if(!recipe$Analysis$Regrid$type %in% c('to_system','to_reference')){
+  if (!recipe$Analysis$Regrid$type %in% c('to_system', 'to_reference')) {
   
     aux <- CDORemap(data_array = hcst$data, 
                     lons = hcst$coords$longitude, lats = hcst$coords$latitude,