repos <- paste0('/esarchive/exp/ecmwf/system4_m1/6hourly/',
                  '$var$/$var$_$sdate$.nc')
  
  system4 <- Start(dat = repos,
                   var = 'tas',
                   sdate = sapply(1994, function(x) paste0(x, sprintf('%02d', 5:12), '01')),
                   time = indices(seq(1, 124, 4)),  #first time step per day
                   ensemble = 'all',
                   latitude = indices(1:10),
                   longitude = indices(1:10),
                   return_vars = list(latitude = NULL,
                                      longitude = NULL,
                                      time = c('sdate')))

#-------- Check exp data -----------
 attr(system4, 'Dimensions')
#      dat       var     sdate      time  ensemble  latitude longitude 
#        1         1         8        31        51        10        10 
#-----------------------------------

# ------- retrieve the attributes for obs load-in ----------
  dates <- attr(system4, 'Variables')$common$time
#> dim(dates)
#sdate  time 
#    8    31 
  dates_file <- sort(unique(gsub('-', '', sapply(as.character(dates),
  substr, 1, 7))))
#> dates_file
#[1] "199405" "199406" "199407" "199408" "199409" "199410" "199411" "199412"
# -----------------------------------------------------------
  
  repos_obs <- paste0('/esarchive/recon/ecmwf/erainterim/6hourly/',
                      '$var$/$var$_$file_date$.nc')
  erai <-    Start(dat = repos_obs,
                   var = 'tas',
                   file_date = dates_file,
                   time = values(dates),  #use dates from exp
                   latitude = indices(1:10),
                   longitude = indices(1:10),
                   time_var = 'time',
                   #because time is assigned by 'values', set the tolerance to avoid too distinct match
                   time_tolerance = as.difftime(1, units = 'hours'), 
                   #time is defined by dates, which dimension is [sdate = 8, time = 31]
                   time_across = 'file_date',  
                   return_vars = list(latitude = NULL,
                                      longitude = NULL,
                                      time = 'file_date'),
                   #combine time and file_date
                   merge_across_dims = TRUE,
                   #exclude the additional NAs generated by merge_across_dims
                   merge_across_dims_narm = TRUE,
                   #split time, because it is two-dimensional
                   split_multiselected_dims = TRUE,
                   retrieve = TRUE)

#------- Check erai -----------
dim(erai)
#      dat       var     sdate      time  latitude longitude 
#        1         1         8        31        10        10 

erai[1, 1, 1, , 1, 1]
# [1] 255.0120 256.8095 254.3654 254.6059 257.0551 255.5087 256.8167 257.9717
# [9] 258.7491 259.2942 259.6682 260.7215 260.0988 261.2605 263.3590 265.6683
#[17] 262.4813 262.6136 263.0591 262.8377 261.7276 263.9910 264.7755 266.0213
#[25] 268.5927 267.8699 268.9210 269.4702 267.6735 267.9255 268.2216

erai[1, 1, , 1, 2, 2]
#[1] 254.5793 269.6221 274.5021 274.0269 269.5855 253.7458 243.9750 244.2415

#------------------------------