diff --git a/modules/Loading/Loading_decadal.R b/modules/Loading/Loading_decadal.R
index 43b3c54fef12287133915203af83bb036c5356d5..74c97e299f7e7e232fbf2fb3e261bec9673a394e 100644
--- a/modules/Loading/Loading_decadal.R
+++ b/modules/Loading/Loading_decadal.R
@@ -133,7 +133,8 @@ load_datasets <- function(recipe) {
   if (length(variable) > 1) {
     Start_default_arg_list <- c(Start_default_arg_list, 
       list(table = table, grid = grid, version = version,
-           table_depends = 'var', grid_depends = 'var', version_depends = 'var'))
+           table_depends = 'var', grid_depends = 'var', version_depends = 'var',
+           metadata_dims = 'var'))
   }
 
   if (!multi_path) {
@@ -277,9 +278,13 @@ load_datasets <- function(recipe) {
 #-------------------------------------------
 # Step 3. Load the reference
 #-------------------------------------------
-  obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq,
-                        paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]]))
-  obs.files <- paste0('$var$_$file_date$.nc')
+  obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src,
+                        store.freq, "$var$$var_dir$", "$var$_$file_date$.nc")
+  var_dir_obs <- archive$Reference[[ref.name]][[store.freq]][variable] # list(tas = "_f1h-r1440x721cds", tos = "_f1h-r1440x721cds")
+
+#  obs.path <- file.path(archive$src, archive$Reference[[ref.name]]$src, store.freq,
+#                        paste0(variable, archive$Reference[[ref.name]][[store.freq]][[variable]]))
+#  obs.files <- paste0('$var$_$file_date$.nc')
 
   # Get from startR_cube
 #  dates <- attr(hcst, 'Variables')$common$time
@@ -299,8 +304,10 @@ load_datasets <- function(recipe) {
   # Restore correct dimensions
   dim(dates) <- dim(dates_file)
 
-  obs <- Start(dat = file.path(obs.path, obs.files),
+  obs <- Start(dat = obs.path,
                var = variable,
+               var_dir = var_dir_obs,
+               var_dir_depends = 'var',
                file_date = unique(format(dates, '%Y%m')),
                time = dates, # [sday, sweek, syear, time]
                time_across = 'file_date',
@@ -327,8 +334,10 @@ load_datasets <- function(recipe) {
 # Method 2: reshape hcst time attr's date into an array with time dim then as obs date selector
 #////////////////
 
-  obs <- Start(dat = file.path(obs.path, obs.files),
+  obs <- Start(dat = obs.path,
                var = variable,
+               var_dir = var_dir_obs,
+               var_dir_depends = 'var',
                file_date = dates_file, #dates_arr, # [sday, sweek, syear, time]
                split_multiselected_dims = TRUE,
                latitude = values(list(lats.min, lats.max)),
@@ -344,6 +353,7 @@ load_datasets <- function(recipe) {
                                longitude = c('lon','longitude')),
                return_vars = list(latitude = NULL, longitude = NULL,
                                   time = 'file_date'),
+               metadata_dims = 'var',
                silent = !DEBUG,
                retrieve = TRUE)
   }
@@ -353,6 +363,9 @@ load_datasets <- function(recipe) {
 # sday sweek syear  time 
 #    1     1     2    14 
 
+  # Remove var_dir dimension
+  obs <- Subset(obs, along = "var_dir", indices = 1, drop = "selected")
+
   # Only ensemble dim could be different
   if (!identical(dim(obs), dim(hcst$data)[-9])) {
     error(recipe$Run$logger,
@@ -480,33 +493,42 @@ load_datasets <- function(recipe) {
 #-------------------------------------------
   # Remove negative values in accumulative variables
   dictionary <- read_yaml("conf/variable-dictionary.yml")
-  if (dictionary$vars[[variable]]$accum) {
-    info(recipe$Run$logger, 
-	 " Accumulated variable: setting negative values to zero.")
-    obs$data[obs$data < 0] <- 0
-    hcst$data[hcst$data < 0] <- 0
-    if (!is.null(fcst)) {
-      fcst$data[fcst$data < 0] <- 0
-    }
-  }
-  # Convert precipitation to mm/day
-  ## TODO: Make a function?
-  if (variable == "prlr") {
-    # Verify that the units are m/s and the same in obs and hcst
-    if (((attr(obs$Variable, "variable")$units == "m s-1") ||
-	 (attr(obs$Variable, "variable")$units == "m s**-1")) &&
-	((attr(hcst$Variable, "variable")$units == "m s-1") ||
-	 (attr(hcst$Variable, "variable")$units == "m s**-1"))) {
-      
+  for (var_idx in 1:length(variable)) {
+    var_name <- variable[var_idx]
+    if (dictionary$vars[[var_name]]$accum) {
       info(recipe$Run$logger,
-	   "Converting precipitation from m/s to mm/day.")
-      obs$data <- obs$data*86400*1000
-      attr(obs$Variable, "variable")$units <- "mm/day"
-      hcst$data <- hcst$data*86400*1000
-      attr(hcst$Variable, "variable")$units <- "mm/day"
+           paste0("Accumulated variable ", var_name,
+                  ": setting negative values to zero."))
+      # obs$data[, var_idx, , , , , , , ] <- pmax(Subset(obs$data,
+      #                                                along = "var",
+      #                                                indices = var_idx, F), 0)
+      obs$data[, var_idx, , , , , , , ][obs$data[, var_idx, , , , , , , ] < 0] <- 0
+      hcst$data[, var_idx, , , , , , , ][hcst$data[, var_idx, , , , , , , ] < 0] <- 0
       if (!is.null(fcst)) {
-        fcst$data <- fcst$data*86400*1000
-        attr(fcst$Variable, "variable")$units <- "mm/day"
+        fcst$data[, var_idx, , , , , , , ][fcst$data[, var_idx, , , , , , , ] < 0] <- 0
+      }
+    }
+    
+    # Convert prlr from m/s to mm/day
+    ## TODO: Make a unit conversion function
+    if (variable[[var_idx]] == "prlr") {
+      # Verify that the units are m/s and the same in obs and hcst
+      if (((obs$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+           (obs$attrs$Variable$metadata[[var_name]]$units == "m s**-1")) &&
+          ((hcst$attrs$Variable$metadata[[var_name]]$units == "m s-1") ||
+           (hcst$attrs$Variable$metadata[[var_name]]$units == "m s**-1"))) {
+        info(recipe$Run$logger, "Converting precipitation from m/s to mm/day.")
+        obs$data[, var_idx, , , , , , , ] <-
+          obs$data[, var_idx, , , , , , , ]*86400*1000
+        obs$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+        hcst$data[, var_idx, , , , , , , ] <-
+          hcst$data[, var_idx, , , , , , , ]*86400*1000
+        hcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+        if (!is.null(fcst)) {
+          fcst$data[, var_idx, , , , , , , ] <- 
+            fcst$data[, var_idx, , , , , , , ]*86400*1000
+          fcst$attrs$Variable$metadata[[var_name]]$units <- "mm/day"
+        }
       }
     }
   }