diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b1cf9976e0a017b4076e024242f9c1a9a7e61487..39a31ecb9326dc181d17a8f0f7ad1da5a68bce8f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,7 +4,7 @@ build:
   stage: build
   script:
     - module load R
-#    - module load CDO
+    - module load CDO
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual --run-donttest s2dverification_*.tar.gz
     - R -e 'covr::package_coverage()'
diff --git a/NAMESPACE b/NAMESPACE
index 4e8557898444268933dfffceb4f5d7d3c8460183..8561a4f02ff2318d1146b35d48cd4102ada054ce 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -93,6 +93,7 @@ import(methods)
 import(ncdf4)
 import(parallel)
 import(plyr)
+importFrom(abind,abind)
 importFrom(grDevices,bmp)
 importFrom(grDevices,col2rgb)
 importFrom(grDevices,colorRampPalette)
diff --git a/R/Load.R b/R/Load.R
index cded9565fe77ae89528bd4c701d171982647c795..b5600af2a0e2b7c5383b71bf7a20fffd426e5081 100644
--- a/R/Load.R
+++ b/R/Load.R
@@ -840,6 +840,7 @@
 #'  } 
 #'@import parallel bigmemory methods
 #'@importFrom stats ts window na.omit
+#'@importFrom abind abind
 #'@export
 Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL, 
                  nmemberobs = NULL, nleadtime = NULL, leadtimemin = 1, 
@@ -850,6 +851,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
                  configfile = NULL, varmin = NULL, varmax = NULL, 
                  silent = FALSE, nprocs = NULL, dimnames = NULL, 
                  remapcells = 2, path_glob_permissive = 'partial') {
+
   #library(parallel)
   #library(bigmemory)
 
@@ -1052,7 +1054,6 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   if (!is.character(storefreq) || !(storefreq %in% c('monthly', 'daily'))) {
     stop("Error: parameter 'storefreq' is wrong, can take value 'daily' or 'monthly'.")
   }
-
   # sampleperiod
   if (is.null(sampleperiod) || !is.numeric(sampleperiod)) {
     stop("Error: parameter 'sampleperiod' is wrong. It should be numeric.")
@@ -1382,6 +1383,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   dims2define <- TRUE
   is_file_per_member_exp <- rep(nmod, FALSE)
   exp_work_pieces <- list()
+  first_time_step_list <- NULL
   jmod <- 1
   while (jmod <= nmod) {
     first_dataset_file_found <- FALSE
@@ -1427,6 +1429,21 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       if (is_file_per_member_exp[jmod]) {
         replace_values[["MEMBER_NUMBER"]] <- '*'
       }
+      if (jsdate == 1) {
+        work_piecetime <- list(dataset_type = dataset_type,
+                           filename = .ConfigReplaceVariablesInString(quasi_final_path,
+                                                                      replace_values),
+                           namevar = namevar, grid = grid, remap = remap, 
+                           remapcells = remapcells,
+                           is_file_per_member = is_file_per_member_exp[jmod],
+                           is_file_per_dataset = FALSE,
+                           lon_limits = c(lonmin, lonmax),
+                           lat_limits = c(latmin, latmax), dimnames = exp[[jmod]][['dimnames']],
+                           single_dataset = single_dataset)
+        looking_time <- .LoadDataFile(work_piecetime, explore_dims = TRUE, 
+                                      silent = silent)
+        first_time_step_list <- c(first_time_step_list, list(looking_time$time_dim))
+      }
       # If the dimensions of the output matrices are still to define, we try to read
       # the metadata of the data file that corresponds to the current iteration
       if (dims2define) {
@@ -1531,6 +1548,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       jsdate <- jsdate + 1
     }
     replace_values[extra_vars] <- NULL
+    #first_dataset_file_found <- FALSE
     jmod <- jmod + 1
   }
   if (dims2define && length(exp) > 0) {
@@ -1553,7 +1571,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   # the current date.
   if (is.null(exp) || dims == 0) {
     if (is.null(leadtimemax)) {
-      diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d')
+      diff <- Sys.time() - as.POSIXct(sdates[1], format = '%Y%m%d', tz = "UTC")
       if (diff > 0) {
         .warning("Loading observations only and no 'leadtimemax' specified. Data will be loaded from each starting date to current time.")
         if (storefreq == 'monthly') { 
@@ -1571,7 +1589,52 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     }
     leadtimes <- seq(leadtimemin, leadtimemax, sampleperiod)
   }
-  
+  # If there are differences in the first time stamp in exp files:
+  if (!is.null(exp)) {
+    in_date <- lapply(first_time_step_list, function(x) {
+        origin <- as.POSIXct(
+          paste(strsplit(x$time_units, " ")[[1]][c(3,4)],
+                collapse = " "), tz = 'UTC')
+        units <- strsplit(x$time_units, " ")[[1]][1]
+        if (units == 'hours') {
+          exp_first_time_step <- as.POSIXct(
+            x$first_time_step_in_file *
+            3600, origin = origin, tz = 'UTC')
+        } else if (units == 'days') {
+          exp_first_time_step <- as.POSIXct(
+            x$first_time_step_in_file *
+            86400, origin = origin, tz = 'UTC')
+        }
+        day <- as.numeric(format(exp_first_time_step, "%d"))
+        return(day)
+    })
+    exp_first_time_step <- min(unlist(in_date))
+    if (max(unlist(in_date)) > 1) {
+      leadtimes <- seq(exp_first_time_step, leadtimemax + max(unlist(in_date)) - 1, 
+                       sampleperiod)
+    }
+    if (leadtimemin > 1 & length(in_date) > 1) {
+      lags <- lapply(in_date, function(x) {x - in_date[[1]]})
+      new_leadtimemin <- lapply(lags, function(x) {leadtimemin - x})
+      new_leadtimemax <- lapply(lags, function(x) {leadtimemax - x})
+      jmod <- 2
+      npieces <- length(exp_work_pieces)/nmod
+      while (jmod <= nmod) {
+        jpiece <- 1
+        while (jpiece <= npieces) {
+          exp_work_pieces[[npieces * (jmod - 1) + jpiece]]$leadtimes <- 
+            seq(new_leadtimemin[[jmod]], new_leadtimemax[[jmod]], sampleperiod)
+          jpiece <- jpiece + 1
+        }
+        jmod <- jmod + 1
+      }
+    }
+    lag <- 1 - in_date[[1]]
+    leadtimes <- seq(leadtimemin - lag, leadtimemax #+ max(unlist(in_date)) + lag,
+                     - lag,
+                     sampleperiod)
+    exp_first_time_step <- leadtimemin - lag
+  }
   # Now we start iterating over observations. We try to find the output matrix
   # dimensions and we build anyway the work pieces corresponding to the observational
   # data that time-corresponds the experimental data or the time-steps until the
@@ -1635,6 +1698,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
         found_data <- .LoadDataFile(work_piece, explore_dims = TRUE, silent = silent)
         found_dims <- found_data$dims
         var_long_name <- found_data$var_long_name
+        first_time_step_list <- c(first_time_step_list, list(found_data$time_dim))
         units <- found_data$units
         if (!is.null(found_dims)) {
           is_2d_var <- found_data$is_2d_var
@@ -1709,7 +1773,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
           }
           day <- as.integer(day)
           startdate <- as.POSIXct(paste(substr(sdate, 1, 4), '-',
-                       substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = '')) + 
+                       substr(sdate, 5, 6), '-', day, ' 12:00:00', sep = ''), tz = "UTC") + 
                        (leadtimemin - 1) * 86400
           year <- as.integer(substr(startdate, 1, 4))
           month <- as.integer(substr(startdate, 6, 7))
@@ -1732,7 +1796,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
             ## This condition must be fulfilled to put all the month time steps
             ## in the dimension of length nleadtimes. Otherwise it must be cut:
             #(length(leadtimes) - 1)*sampleperiod + 1 - (jleadtime - 1)*sampleperiod >= days_in_month - day + 1
-            obs_file_indices <- seq(day, min(days_in_month, (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod)
+            
+            ## The first time step in exp could be different from sdate:
+            if (jleadtime == 1 & !is.null(exp)) {
+             if (is.null(first_time_step_list[[1]])) {
+                stop("Check 'time' variable in the experimental files ",
+                     "since not units or first time step have been found.")
+              } else {
+                  day <- leadtimes[1]
+              }
+            }
+          obs_file_indices <- seq(day, min(days_in_month,
+                  (length(leadtimes) - jleadtime) * sampleperiod + day), sampleperiod)
           } else {
             obs_file_indices <- 1
           }
@@ -1828,7 +1903,8 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
           }
           
           if (storefreq == 'daily') {
-            startdate <- startdate + 86400 * sampleperiod * length(obs_file_indices)
+            startdate <- startdate + 86400 * sampleperiod * 
+                         max(obs_file_indices)
             year <- as.integer(substr(startdate, 1, 4))
             month <- as.integer(substr(startdate, 6, 7))
             day <- as.integer(substr(startdate, 9, 10))
@@ -2142,49 +2218,54 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     attr(variable, 'daily_agg_cellfun') <- 'none'
     attr(variable, 'monthly_agg_cellfun') <- 'none'
     attr(variable, 'verification_time') <- 'none'
-
-    number_ftime <- NULL 
+    
+    number_ftime <- NULL
     if (is.null(var_exp)) {
       mod_data <- NULL
-    } else { 
+    } else {
       dim_reorder <- length(dim_exp):1
       dim_reorder[2:3] <- dim_reorder[3:2]
       old_dims <- dim_exp
       dim_exp <- dim_exp[dim_reorder]
-      mod_data <- aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder)
+      mod_data <-
+        aperm(array(bigmemory::as.matrix(var_exp), dim = old_dims), dim_reorder)
       attr(mod_data, 'dimensions') <- names(dim_exp)
       names(dim(mod_data)) <- names(dim_exp)
       number_ftime <- dim_exp[["ftime"]]
     }
-
+    
     if (is.null(var_obs)) {
       obs_data <- NULL
-    } else {                     
+    } else {
       dim_reorder <- length(dim_obs):1
       dim_reorder[2:3] <- dim_reorder[3:2]
       old_dims <- dim_obs
       dim_obs <- dim_obs[dim_reorder]
-      obs_data <- aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder)
+      obs_data <-
+        aperm(array(bigmemory::as.matrix(var_obs), dim = old_dims), dim_reorder)
       attr(obs_data, 'dimensions') <- names(dim_obs)
       names(dim(obs_data)) <- names(dim_obs)
       if (is.null(number_ftime)) {
         number_ftime <- dim_obs[["ftime"]]
       }
     }
-
     if (is.null(latitudes)) {
       lat <- 0
       attr(lat, 'cdo_grid_name') <- 'none'
       attr(lat, 'first_lat') <- 'none'
-      attr(lat, 'last_lat') <- 'none' 
+      attr(lat, 'last_lat') <- 'none'
     } else {
       lat <- latitudes
-      attr(lat, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid
+      attr(lat, 'cdo_grid_name') <-
+        if (is.null(grid))
+          'none'
+      else
+        grid
       attr(lat, 'first_lat') <- tail(lat, 1)
       attr(lat, 'last_lat') <- head(lat, 1)
     }
     attr(lat, 'projection') <- 'none'
-
+    
     if (is.null(longitudes)) {
       lon <- 0
       attr(lon, 'cdo_grid_name') <- 'none'
@@ -2194,14 +2275,18 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       attr(lon, 'last_lon') <- 'none'
     } else {
       lon <- longitudes
-      attr(lon, 'cdo_grid_name') <- if (is.null(grid)) 'none' else grid
+      attr(lon, 'cdo_grid_name') <-
+        if (is.null(grid))
+          'none'
+      else
+        grid
       attr(lon, 'data_across_gw') <- data_across_gw
       attr(lon, 'array_across_gw') <- array_across_gw
       attr(lon, 'first_lon') <- lon[which.min(abs(lon - lonmin))]
       attr(lon, 'last_lon') <- lon[which.min(abs(lon - lonmax))]
     }
     attr(lon, 'projection') <- 'none'
-
+    
     dates <- list()
     ## we must put a start and end time for each prediction c(start date, forecast time)
     if (storefreq == 'minutely') {
@@ -2213,22 +2298,33 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
     } else if (storefreq == 'monthly') {
       store_period <- 'month'
     }
-
+    
     addTime <- function(date, period, n = 1) {
       seq(date, by = paste(n, period), length = 2)[2]
     }
-
+    
     # We build dates, a list with components start and end.
     # Start is a list with as many components as start dates.
     # Each component is a vector of the initial POSIXct date of each
     # forecast time step
+    if (!is.null(exp)) {
+      if (storefreq == 'daily' & leadtimes[[1]] > 1) {
+        origin <- leadtimes[[1]] - 1
+        leadtimemin <- 1
+      } else {
+        origin <- 0
+      }
+    } else {
+      origin <- 0
+    }
     dates[["start"]] <- do.call(c, lapply(sdates,
       function(x) {
-        do.call(c, lapply((0:(number_ftime - 1)) * sampleperiod,
+        do.call(c, lapply((origin:(origin + number_ftime - 1)) * sampleperiod,
           function(y) {
-            addTime(as.POSIXct(x, format = "%Y%m%d"), store_period, y + leadtimemin - 1)
+            addTime(as.POSIXct(x, format = "%Y%m%d", tz = "UTC"), store_period, y + leadtimemin - 1)
           }))
       }))
+    attr(dates[["start"]], "tzone") <- "UTC"
     # end is similar to start, but contains the end dates of each forecast 
     # time step
     dates[["end"]] <- do.call(c, lapply(dates[["start"]],
@@ -2238,6 +2334,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
             addTime(y, store_period)
           }))
       }))
+    attr(dates[["end"]], "tzone") <- "UTC"
 
     tags_to_find <- c('START_DATE', 'MEMBER_NUMBER', 'YEAR', 'MONTH', 'DAY')
     position_of_tags <- na.omit(match(tags_to_find, names(replace_values)))
@@ -2249,14 +2346,19 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
       models <- list()
       for (jmod in 1:length(exp)) {
         member_names <- paste0("Member_", 1:nmember[jmod])
-        models[[exp[[jmod]][["name"]]]] <- list(
-          InitializationDates = lapply(member_names,
+
+        InitializationDates <- do.call(c, lapply(member_names,
             function(x) {
               do.call(c, lapply(sdates, function(y) {
-                as.POSIXct(y, format = "%Y%m%d")
+                as.POSIXct(y, format = "%Y%m%d", tz = "UTC")
               }))
-            }),
+            }))
+        attr(InitializationDates, "tzone") <- "UTC"
+
+        models[[exp[[jmod]][["name"]]]] <- list(
+          InitializationDates = InitializationDates,
           Members = member_names)
+
         names(models[[exp[[jmod]][["name"]]]]$InitializationDates) <- member_names
         attr(models[[exp[[jmod]][["name"]]]], 'dataset') <- exp[[jmod]][["name"]]
         attr(models[[exp[[jmod]][["name"]]]], 'source') <- {
diff --git a/R/Utils.R b/R/Utils.R
index 4460018139d3bf55bfd8cc3f9b1b2a4f6f17f850..a7415c512985a01be2853b7a04a2532c6679615c 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -591,6 +591,9 @@
         nltime <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$len
         expected_dims <- c(expected_dims, time_dimname)
         dim_matches <- match(expected_dims, var_dimnames)
+        first_time_step_in_file <- fnc$var[[namevar]][['dim']][[match(time_dimname,
+                                        var_dimnames)]]$vals[1]
+        time_units <- fnc$var[[namevar]][['dim']][[match(time_dimname, var_dimnames)]]$units
       } else {
         if (!is.null(old_members_dimname)) {
           expected_dims[which(expected_dims == 'lev')] <- old_members_dimname
@@ -631,7 +634,7 @@
         time_indices <- ts(time_indices, start = c(years[1], mons[1]), 
                            end = c(years[length(years)], mons[length(mons)]),
                            frequency = 12)
-        ltimes_list <- list()
+                ltimes_list <- list()
         for (sdate in work_piece[['startdates']]) {
           selected_time_indices <- window(time_indices, start = c(as.numeric(
                                    substr(sdate, 1, 4)), as.numeric(substr(sdate, 5, 6))), 
@@ -938,7 +941,9 @@
   if (explore_dims) {
     list(dims = dims, is_2d_var = is_2d_var, grid = grid_name, 
          units = units, var_long_name = var_long_name, 
-         data_across_gw = data_across_gw, array_across_gw = array_across_gw)
+         data_across_gw = data_across_gw, array_across_gw = array_across_gw,
+         time_dim = list(first_time_step_in_file = first_time_step_in_file,
+                        time_units = time_units))
   } else {
     ###if (!silent && !is.null(progress_connection) && !is.null(work_piece[['progress_amount']])) {
     ###  foobar <- writeBin(work_piece[['progress_amount']], progress_connection)
diff --git a/man/Composite.Rd b/man/Composite.Rd
index 3d9a3fd8eba3c149af352eb0027b8688dce3e8c4..b957d9fdb0cd71cab5db8da10f12313bdae30c00 100644
--- a/man/Composite.Rd
+++ b/man/Composite.Rd
@@ -84,7 +84,6 @@ data <- 1 : (4 * 5 * 6)
 dim(data) <- c(lon = 4, lat = 5, case = 6)
 occ <- c(1, 1, 2, 2, 3, 3) 
 res <- Composite(data, occ, K = 4)
-
 }
 \author{
 History:
diff --git a/tests/testthat/test-Load.R b/tests/testthat/test-Load.R
new file mode 100644
index 0000000000000000000000000000000000000000..645ceaae964b81e6bf30591598ab273c0cf7881a
--- /dev/null
+++ b/tests/testthat/test-Load.R
@@ -0,0 +1,274 @@
+context("Testing Load with ECMWF System5c3s daily data")
+test_that("First time step is correctly interpreted:", {
+  system5 <- Load(var = 'prlr',
+                  exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+                  sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+                  nmember = 5,leadtimemin = 1,
+                  leadtimemax = NULL,
+                  storefreq = "daily", sampleperiod = 1,
+                  latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+                  output = 'lonlat', nprocs = 1,
+                  method = 'conservative', grid = 'r360x181',
+                  maskmod = vector("list", 15), maskobs = vector("list", 15),
+                  configfile = NULL, varmin = NULL, varmax = NULL,
+                  silent = FALSE, dimnames = NULL,
+                  remapcells = 2, path_glob_permissive = 'partial',
+                  nmemberobs = NULL, nleadtime = NULL)
+
+  
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  
+  expect_equal(dim(data$mod), c(dataset = 1, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'), 
+                                       as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+                                   seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                                       as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$end, dates)
+  
+  expect_equal(data$mod, system5$mod)
+  
+  era5 <- Load(var = 'prlr', exp = NULL,
+              obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'),
+               nmember = 1, leadtimemin = 1,
+               leadtimemax = 215,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = 1, nleadtime = NULL)
+  # The values of observation when loading exp and obs simultaneously corresponds to the second day of the observations in the file:
+  expect_equal(data$obs[1,1,1,,,], era5$obs[1,1,1,2:215,,])
+  expect_equal(data$obs[1,1,2,,,], era5$obs[1,1,2,2:215,,]) 
+
+ # Test when 2 observational datasets are requested:
+  data <- Load(var = 'prlr',
+             exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+             obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc'), list(name = 'erainterim')),
+             sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+             nmember = 5, leadtimemin = 1,
+             leadtimemax = NULL,
+             storefreq = "daily", sampleperiod = 1,
+             latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+             output = 'lonlat', nprocs = 1,
+             method = 'conservative', grid = 'r360x181',
+             maskmod = vector("list", 15), maskobs = vector("list", 15),
+             configfile = NULL, varmin = NULL, varmax = NULL,
+             silent = FALSE, dimnames = NULL,
+             remapcells = 2, path_glob_permissive = 'partial',
+             nmemberobs = NULL, nleadtime = NULL)  
+  expect_equal(dim(data$obs), c(dataset = 2, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 1, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(data$obs[1,1,,,,], era5$obs[1,1,,2:215,,])
+  expect_equal(system5$mod, data$mod)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(dates, data$Dates$start)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'  
+  expect_equal(dates, data$Dates$end)
+  
+  eraint <- Load(var = 'prlr', exp = NULL,
+              obs = list(list(name = 'erainterim')),
+              sdates = c('19931101', '19941101'),
+              nmember = 1, leadtimemin = 1,
+              leadtimemax = 215,
+              storefreq = "daily", sampleperiod = 1,
+              latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+              output = 'lonlat', nprocs = 1,
+              method = 'conservative', grid = 'r360x181',
+              maskmod = vector("list", 15), maskobs = vector("list", 15),
+              configfile = NULL, varmin = NULL, varmax = NULL,
+              silent = FALSE, dimnames = NULL,
+              remapcells = 2, path_glob_permissive = 'partial',
+              nmemberobs = 1, nleadtime = NULL)
+  expect_equal(data$obs[2,1,,,,], eraint$obs[1,1,,2:215,,])
+
+  # Test for 2 experimental datasets
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(sum(is.na(data$mod)), 0)
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                                       as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+                                   seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                                       as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-04", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-04", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$end, dates)
+  expect_equal(data$obs[,,,,,], era5$obs[,,,2:215,,])
+  expect_equal(system5$exp, data$exp[1,,,,,])
+  
+  system4 <- Load(var = 'prlr',
+                  exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 1, leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+
+  expect_equal(system4$exp, data$exp[2,,,,,])
+  expect_equal(system4$Dates$start[c(-1,-216)], data$Dates$start)
+  
+  
+  # Test for 2 experimental datasets oposite order
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"), 
+                          list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 1,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 215,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 215,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod[1,,,,,])), 0)
+  expect_equal(sum(is.na(data$mod[2,,,-215,,])),0)
+  expect_equal(sum(is.na(data$mod[2,,,215,,])), 640)
+  expect_equal(data$Dates$start, system4$Dates$start)
+  expect_equal(data$obs, era5$obs)
+  # Test leadtimemin > 1
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 2,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 213,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2,ftime = 213,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(system4$mod[,1,,3:215,,], data$mod[2,1,,,,])
+  expect_equal(system5$mod[,,,2:214,,], data$mod[1,,,,,]) 
+  expect_equal(era5$obs[,,,3:215,,], data$obs[1,,,,,], tolerance = 0.001)
+  dates <- c(seq(as.POSIXct("1993-11-03", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-03", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+  
+  data <- Load(var = 'prlr',
+               exp = list(list(name = 'ecmwfS4', path = "/esarchive/exp/ecmwf/system4_m1/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$START_DATE$.nc"),
+                          list(name = 'ecmwfS5', path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")),
+               obs = list(list(name = 'era5', path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')),
+               sdates = c('19931101', '19941101'), #paste0(1993:2018, '1101')
+               nmember = 5,leadtimemin = 2,
+               leadtimemax = NULL,
+               storefreq = "daily", sampleperiod = 1,
+               latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+               output = 'lonlat', nprocs = 1,
+               method = 'conservative', grid = 'r360x181',
+               maskmod = vector("list", 15), maskobs = vector("list", 15),
+               configfile = NULL, varmin = NULL, varmax = NULL,
+               silent = FALSE, dimnames = NULL,
+               remapcells = 2, path_glob_permissive = 'partial',
+               nmemberobs = NULL, nleadtime = NULL)
+  expect_equal(system4$mod[1,1,,2:215,,], data$mod[1,1,,1:214,,])
+  expect_equal(system5$mod[1,1,,1:214,,], data$mod[2,1,,1:214,,])
+  expect_equal(dim(data$obs), c(dataset = 1, member = 1, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(dim(data$mod), c(dataset = 2, member = 5, sdate = 2, ftime = 214,
+                                lat = 8, lon = 8))
+  expect_equal(sum(is.na(data$obs)), 0)
+  expect_equal(sum(is.na(data$mod)), 0)
+  expect_equal(data$obs[,,,1,,], era5$obs[,,,2,,])
+  dates <- c(seq(as.POSIXct("1993-11-02", tz = 'UTC'),
+                 as.POSIXct("1994-06-03", tz = 'UTC'), "d"),
+             seq(as.POSIXct("1994-11-02", tz = 'UTC'),
+                 as.POSIXct("1995-06-03", tz = 'UTC'), "d"))
+  attributes(dates)$tzone <- 'UTC'
+  expect_equal(data$Dates$start, dates)
+})