diff --git a/R/Start.R b/R/Start.R
index 702b77633cfda496eedf865bce1fb5b75b849c5a..92eb16dc107628e8e4d04e666c74878f9c373d54 100644
--- a/R/Start.R
+++ b/R/Start.R
@@ -4108,7 +4108,90 @@ Start <- function(..., # dim = indices/selectors,
         }
       }
     }
+    # Retrieve variable metadata
+    # Compare array_of_metadata_flags with array_of_files_to_load to know which files to take for metadata
+    if (!is.null(metadata_dims)) {
+      array_of_metadata_flags <- array(FALSE, dim = dim(array_of_files_to_load))
+      metadata_indices_to_load <- as.list(rep(1, length(dim(array_of_files_to_load))))
+      names(metadata_indices_to_load) <- names(dim(array_of_files_to_load))
+      metadata_indices_to_load[metadata_dims] <- as.list(rep(TRUE, length(metadata_dims)))
+      array_of_metadata_flags <- do.call('[<-', c(list(array_of_metadata_flags),  metadata_indices_to_load,
+                                                  list(value = rep(TRUE, prod(dim(array_of_files_to_load)[metadata_dims])))))
+
+      if (tail(names(dim(array_of_files_to_load)), 1) != found_pattern_dim) {
+        tmp1 <- s2dv::Reorder(array_of_files_to_load, c(2:length(dim(array_of_files_to_load)), 1))
+        tmp2 <- s2dv::Reorder(array_of_metadata_flags, c(2:length(dim(array_of_metadata_flags)), 1))
+        files_for_metadata <- tmp1[tmp2]
+      } else {
+        files_for_metadata <- array_of_files_to_load[array_of_metadata_flags]
+      }
 
+      # Get variable name
+      #NOTE: This part probably will fail when one netCDF file has more than one variable.
+      if (found_pattern_dim %in% metadata_dims) {  # metadata_dims has "dat"
+        if (any(metadata_dims %in% c('var', 'variable'))) { # metadata_dim is c('dat', 'var')
+          how_many_vars <-  length(dat[[1]][['selectors']]$var[[1]])
+        } else if (length(metadata_dims) > 1) {  # metadata_dims is c('dat', xxx)
+          how_many_vars <- length(dat[[1]][['selectors']][[metadata_dims[which(found_pattern_dim != metadata_dims)]]][[1]])
+        } else {  # metadata_dims is 'dat'
+          how_many_vars <- 1
+        }
+        tmp_var <- matrix(NA, how_many_vars, length(dat))
+        for (i_dat in 1:dim(array_of_metadata_flags)[found_pattern_dim]) {
+          if (any(metadata_dims %in% c('var', 'variable'))) { # metadata_dims has "var"
+            tmp_var[, i_dat] <- dat[[i_dat]][['selectors']]$var[[1]]
+          } else if (length(metadata_dims) > 1) {  # metadata_dims is c('dat', xxx)
+            tmp_var[, i_dat] <- rep(dat[[i_dat]][['selectors']]$var[[1]][1],
+                                    length(dat[[1]][['selectors']][[metadata_dims[which(found_pattern_dim != metadata_dims)]]][[1]]))
+          } else {  # metadata_dims is 'dat'
+            tmp_var[, i_dat] <- dat[[i_dat]][['selectors']]$var[[1]][1]
+          }
+        }
+
+        # if metadat_dims = c('dat', 'var') and [dat = 2, var = 2], tmp_var has length 4, like c('tas', 'tos', 'tas', 'tos').
+        # if metadata_dims = 'dat' and [dat = 2], tmp_var has length 2 like c('tas', 'tos').
+        tmp_var <- c(tmp_var)
+
+      } else {  # metadata_dims doesn't have "dat"
+        if (any(metadata_dims %in% c('var', 'variable'))) { # metadata_dims has "var"
+          tmp_var <- dat[[1]][['selectors']]$var[[1]]
+        } else {
+          tmp_var <- rep(dat[[1]][['selectors']]$var[[1]][1], length(dat[[1]][['selectors']][[metadata_dims]][[1]]))
+        }
+        # if metadata_dims = 'var' and [var = 2], tmp_var has length 2 like c('tas', 'tos')
+        # if metadata_dims = 'table' and [table = 2], tmp_var has length 1 like 'tas'
+      }
+
+      loaded_metadata <- vector('list', length = length(files_for_metadata))
+      for (i_file in 1:length(files_for_metadata)) {
+        #NOTE: Not use ncatt_get() because it only gets the attr shown with ncdump -h
+        tmp <- file_opener(files_for_metadata[i_file])
+        if (!is.null(tmp)) {  # if file exists
+          loaded_metadata[[i_file]][[1]] <- tmp$var[[tmp_var[i_file]]]
+          names(loaded_metadata[[i_file]]) <- tmp_var[i_file]
+          file_closer(tmp)
+        }
+      }
+      # Find loaded_metadata_files identical as "retrieve = T" case. If dataset_has_files is F, deduct that dataset from counting
+      ind_loaded_metadata_has_values <- which(!sapply(loaded_metadata, is.null)) # c(1, 2, 4)
+      if (!all(dataset_has_files)) {  # If dataset_has_files has F, deduct that dataset from counting
+        if (found_pattern_dim %in% metadata_dims) {  # metadata_dims has "dat" 
+          dataset_has_files_expand <- rep(dataset_has_files, each = how_many_vars)
+          i_ind <- 1
+          while (i_ind <= length(ind_loaded_metadata_has_values)) { # 3, 4, 8
+            if (ind_loaded_metadata_has_values[i_ind] > i_ind) {
+              ind_loaded_metadata_has_values[i_ind] <- ind_loaded_metadata_has_values[i_ind] - length(which(!dataset_has_files_expand[1:dataset_has_files_expand[i_ind]]))
+            }
+            i_ind <- i_ind + 1
+          }
+        }
+      }
+      loaded_metadata_files <- as.character(ind_loaded_metadata_has_values)
+      loaded_metadata <- loaded_metadata[which(!sapply(loaded_metadata, is.null))]
+      return_metadata <- create_metadata_list(array_of_metadata_flags, metadata_dims, pattern_dims,
+                                              loaded_metadata_files, loaded_metadata, dat_names,
+                                              dataset_has_files)
+    }
   }
   # Print the warnings from transform
   if (!is.null(c(warnings1, warnings2, warnings3))) {
@@ -4150,27 +4233,28 @@ Start <- function(..., # dim = indices/selectors,
     file_selectors[[dat[[i]][['name']]]] <- dat[[i]][['selectors']][which(names(dat[[i]][['selectors']]) %in% found_file_dims[[i]])]
   }
 
+  # Prepare attr Variables
+  if (all(sapply(return_metadata, is.null))) {
+    # We don't have metadata of the variable (e.g., tas). The returned metadata list only 
+    # contains those are specified in argument "return_vars".
+    Variables_list <- c(list(common = picked_common_vars), picked_vars)
+    .warning(paste0("Metadata cannot be retrieved. The reason may be the ",
+                    "non-existence of the first file. Use parameter 'metadata_dims'",
+                    " to assign to file dimensions along which to return metadata, ",
+                    "or check the existence of the first file."))
+  } else {
+    # Add the metadata of the variable (e.g., tas) into the list of picked_vars or
+    # picked_common_vars.
+    Variables_list <- combine_metadata_picked_vars(
+                        return_metadata, picked_vars, picked_common_vars,
+                        metadata_dims, pattern_dims, length(dat))
+  }
+
   if (retrieve) {
     if (!silent) {
       .message("Successfully retrieved data.")
     }
 
-    if (all(sapply(return_metadata, is.null))) {
-      # We don't have metadata of the variable (e.g., tas). The returned metadata list only 
-      # contains those are specified in argument "return_vars".
-      Variables_list <- c(list(common = picked_common_vars), picked_vars)
-      .warning(paste0("Metadata cannot be retrieved. The reason may be the ",
-                      "non-existence of the first file. Use parameter 'metadata_dims'",
-                      " to assign to file dimensions along which to return metadata, ",
-                      "or check the existence of the first file."))
-    } else {
-      # Add the metadata of the variable (e.g., tas) into the list of picked_vars or
-      # picked_common_vars.
-      Variables_list <- combine_metadata_picked_vars(
-                          return_metadata, picked_vars, picked_common_vars,
-                          metadata_dims, pattern_dims, length(dat))
-    }
-
     attributes(data_array) <- c(attributes(data_array), 
                                 list(Variables = Variables_list,
                                      Files = array_of_files_to_load, 
@@ -4200,7 +4284,7 @@ Start <- function(..., # dim = indices/selectors,
     start_call[['retrieve']] <- TRUE
     attributes(start_call) <- c(attributes(start_call),
                                 list(Dimensions = final_dims_fake,
-                                     Variables = c(list(common = picked_common_vars), picked_vars),
+                                     Variables = Variables_list,
                                      ExpectedFiles = array_of_files_to_load,
                                      FileSelectors = file_selectors,
                                      PatternDim = found_pattern_dim,
@@ -4237,6 +4321,7 @@ Start <- function(..., # dim = indices/selectors,
   second_round_indices <- work_piece[['second_round_indices']]
   #print("2")
   file_to_open <- work_piece[['file_path']]
+  # Get data and metadata
   sub_array <- file_data_reader(file_to_open, NULL, 
                                 work_piece[['file_selectors']],
                                 first_round_indices, synonims)
diff --git a/R/zzz.R b/R/zzz.R
index 22805ff6e46751fbe4c991a08f5f5bb6a38220b1..1e56e291fa6286d0f15df806c634dc962eb29ac2 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -1363,7 +1363,10 @@ combine_metadata_picked_vars <- function(return_metadata, picked_vars, picked_co
       sublist_names <- lapply(return_metadata, names)[[kk]]
       if (!is.null(sublist_names)) {
         for (jj in 1:length(sublist_names)) {
-          picked_vars[[kk]][[sublist_names[jj]]] <- return_metadata[[kk]][[jj]]
+          if (!is.null(return_metadata[[kk]][[jj]])) {
+            picked_vars[[kk]] <- c(picked_vars[[kk]], list(return_metadata[[kk]][[jj]]))
+            names(picked_vars[[kk]])[length(picked_vars[[kk]])] <- names(return_metadata[[kk]][jj])
+          }
         }
       }
     }
diff --git a/inst/doc/usecase.md b/inst/doc/usecase.md
index 47ee89e526d4a21718f7f9651ad436db76271d0f..80614d01c8420629dad7c930eadd1406f32fd698 100644
--- a/inst/doc/usecase.md
+++ b/inst/doc/usecase.md
@@ -78,7 +78,7 @@ The problem may occur when the dimension number of the splitted selector is more
 
    4. [Use two functions in workflow](inst/doc/usecase/ex2_4_two_func.R)  
 
-   5. [RainFARM precipitation downscaling](https://earth.bsc.es/gitlab/es/startR/-/blob/develop-RainFARMCase/inst/doc/usecase/ex2_5_rainFARM.R)  
+   5. [RainFARM precipitation downscaling](/inst/doc/usecase/ex2_5_rainFARM.R)  
       This example shows how to apply a statistical downscaling function with startR and simultaneously (considering the memory size if unnecessary dimensions are included) saves the data by chunks (e.g., chunking those dimensions which are not required for downscaling) in the esarchive format. It is not recommended to save big outputs. Consider to perform some analysis and then retrieve the result instead of saving data. This is a simplified example of RainFARM for more information visit: https://www.medscope-project.eu/products/data/. 
 Find more explanation of this use case in FAQ [How-to-27](inst/doc/faq.md#27-utilize-chunk-number-in-the-function).
 
diff --git a/tests/testthat/test-Start-metadata_dims.R b/tests/testthat/test-Start-metadata_dims.R
index ce30eec3100f661df7c99f3706af8a314597f366..4251c71d040ec4bdd8b4f19743ce8eb017a87a0a 100644
--- a/tests/testthat/test-Start-metadata_dims.R
+++ b/tests/testthat/test-Start-metadata_dims.R
@@ -19,10 +19,6 @@ suppressWarnings(
                 retrieve = T
                )
 )
-  expect_equal(
-  length(attr(data, 'Variables')),
-  2
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "system5_m1")
@@ -45,6 +41,41 @@ suppressWarnings(
   tolerance = 0.0001
   )
 
+
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system5_m1', path = repos)),
+                var = 'tas',
+                sdate = '20170101',
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:2),
+                lon = indices(2:3),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = 'dat',  # it can be omitted since it is automatically specified as 'dat'
+                retrieve = F
+               )
+)
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  c('time', 'tas')
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$common$tas),
+  22
+  )
+
 })
 
 
@@ -69,10 +100,6 @@ suppressWarnings(
                 retrieve = T
                )
 )
-  expect_equal(
-  length(attr(data, 'Variables')),
-  3
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "system4_m1", "system5_m1")
@@ -107,6 +134,52 @@ suppressWarnings(
   c(247.2570, 248.5016),
   tolerance = 0.0001
   )
+
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = 'tas',
+                sdate = '20170101',
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:2),
+                lon = indices(2:3),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = 'dat',  # it can be omitted since it is automatically specified as 'dat'
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  'time'
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system4_m1),
+  c("lon", "lat", "tas")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat", "tas")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system5_m1$tas),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system4_m1$tas),
+  22
+  )
+
+
 })
 
 test_that("3. One data set, two vars", {
@@ -127,10 +200,6 @@ suppressWarnings(
                 retrieve = TRUE
                )
 )
-  expect_equal(
-  length(attr(data, 'Variables')),
-  2
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "dat1")
@@ -161,6 +230,43 @@ suppressWarnings(
   c(250.00110, 25.04345),
   tolerance = 0.0001
   )
+
+
+suppressWarnings(
+  dataF <- Start(dat = repos,
+                var = var,
+                time = indices(1),
+                lat = indices(9:10),
+                lon = indices(10:11),
+                return_vars = list(lat = NULL, lon = NULL),
+                metadata_dims = 'var',
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "dat1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  c('lat', 'lon', 'tas', 'clt')
+  )
+  expect_equal(
+  is.null(attr(dataF, 'Variables')$dat1),
+  TRUE
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$common$tas),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$common$clt),
+  22
+  )
+
 })
 
 test_that("4. Two data sets, two vars", {
@@ -184,10 +290,6 @@ suppressWarnings(
                 retrieve = T
                )
 )
-  expect_equal(
-  length(attr(data, 'Variables')),
-  3
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "system4_m1", "system5_m1")
@@ -213,6 +315,14 @@ suppressWarnings(
   11
   )
   expect_equal(
+  attr(data, 'Variables')$system4_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 512)
+  )
+  expect_equal(
+  attr(data, 'Variables')$system5_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 1296)
+  )
+  expect_equal(
   data[1, , 1, 1, 1, 2, 2],
   c(247.227219, 6.370782),
   tolerance = 0.0001
@@ -223,6 +333,57 @@ suppressWarnings(
   tolerance = 0.0001
   )
 
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = c('tas', 'sfcWind'),
+                sdate = '20170101',
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:2),
+                lon = indices(1:2),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = 'dat',
+                retrieve = F
+               )
+)
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  'time'
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system4_m1),
+  c("lon", "lat", "tas")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat", "tas")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system5_m1$tas),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system4_m1$tas),
+  22
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system4_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 512)
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system5_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 1296)
+  )
+
 #-------------------------------------------------------------
 suppressWarnings(
   data <- Start(dat = list(list(name = 'system4_m1', path = repos2),
@@ -243,10 +404,6 @@ suppressWarnings(
                )
 )
 
-  expect_equal(
-  length(attr(data, 'Variables')),
-  3
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "system4_m1", "system5_m1")
@@ -289,6 +446,83 @@ suppressWarnings(
   c(248.781540, 5.794801),
   tolerance = 0.0001
   )
+
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = c('tas', 'sfcWind'),
+                sdate = '20170101',
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:2),
+                lon = indices(1:2),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = c('dat', 'var'),
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  'time'
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system4_m1),
+  c("lon", "lat", "tas", "sfcWind")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat", "tas", "sfcWind")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system5_m1$tas),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system5_m1$sfcWind),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system4_m1$tas),
+  22
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system4_m1$sfcWind),
+  22
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system4_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 512)
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system5_m1$tas$dim[[1]][1:2],
+  list(name = 'longitude', len = 1296)
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system4_m1$sfcWind$dim[[1]][1:2],
+  list(name = 'longitude', len = 512)
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system5_m1$sfcWind$dim[[1]][1:2],
+  list(name = 'lon', len = 1296)
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system4_m1$tas$name,
+  'tas'
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system5_m1$sfcWind$name,
+  'sfcWind'
+  )
+
 })
 
 test_that("5. Specify metadata_dims with another file dimension", {
@@ -313,10 +547,6 @@ suppressWarnings(
                )
 )
 
-  expect_equal(
-  length(attr(data, 'Variables')),
-  3
-  )
   expect_equal(
   names(attr(data, 'Variables')),
   c("common", "system4_m1", "system5_m1")
@@ -337,7 +567,167 @@ suppressWarnings(
   length(attr(data, 'Variables')$common$tas),
   12
   )
+  expect_equal(
+  attr(data, 'Variables')$common[[3]]$dim[[4]]$units,
+  "hours since 2016-01-01 00:00:00"
+  )
+
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = c('tas', 'sfcWind'),
+                sdate = c('20150101', '20160101', '20170101'),
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:10),
+                lon = indices(1:10),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = 'sdate',
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  c('time', 'tas', 'tas', 'tas')
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system4_m1),
+  c("lon", "lat")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$common[[2]]),
+  22
+  )
+  expect_equal(
+  attr(data, 'Variables')$common$time,
+  attr(dataF, 'Variables')$common$time
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$common[[3]]$dim[[4]]$units,
+  "hours since 2016-01-01 00:00:00"
+  )
+
+#------------------------------------------------------------------
+suppressWarnings(
+  data <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = c('tas', 'sfcWind'),
+                sdate = c('20150101', '20160101', '20170101'),
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:10),
+                lon = indices(1:10),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = c('dat', 'sdate'),
+                retrieve = T
+               )
+)
+
+  expect_equal(
+  names(attr(data, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(data, 'Variables')$common),
+  c('time')
+  )
+  expect_equal(
+  names(attr(data, 'Variables')$system4_m1),
+  c("lon", "lat", "tas", "tas", "tas")
+  )
+  expect_equal(
+  names(attr(data, 'Variables')$system5_m1),
+  c("lon", "lat", "tas", "tas", "tas")
+  )
+  expect_equal(
+  length(attr(data, 'Variables')$system4_m1[[3]]),
+  12
+  )
+  expect_equal(
+  length(attr(data, 'Variables')$system4_m1[[4]]),
+  12
+  )
+  expect_equal(
+  length(attr(data, 'Variables')$system5_m1[[5]]),
+  12
+  )
+  expect_equal(
+  attr(data, 'Variables')$system4_m1[[3]]$dim[[4]]$units,
+  "hours since 2015-01-01 00:00:00"
+  )
+  expect_equal(
+  attr(data, 'Variables')$system5_m1[[4]]$dim[[4]]$units,
+  "hours since 2016-01-01 00:00:00"
+  )
 
+suppressWarnings(
+  dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = c('tas', 'sfcWind'),
+                sdate = c('20150101', '20160101', '20170101'),
+                ensemble = indices(1),
+                time = indices(1),
+                lat = indices(1:10),
+                lon = indices(1:10),
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = c('dat', 'sdate'),
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  names(attr(dataF, 'Variables')),
+  c("common", "system4_m1", "system5_m1")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$common),
+  c('time')
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system4_m1),
+  c("lon", "lat", "tas", "tas", "tas")
+  )
+  expect_equal(
+  names(attr(dataF, 'Variables')$system5_m1),
+  c("lon", "lat", "tas", "tas", "tas")
+  )
+  expect_equal(
+  length(attr(dataF, 'Variables')$system4_m1[[3]]),
+  22
+  )
+  expect_equal(
+  attr(data, 'Variables')$common$time,
+  attr(dataF, 'Variables')$common$time
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system4_m1[[3]]$dim[[4]]$units,
+  "hours since 2015-01-01 00:00:00"
+  )
+  expect_equal(
+  attr(dataF, 'Variables')$system5_m1[[4]]$dim[[4]]$units,
+  "hours since 2016-01-01 00:00:00"
+  )
 })
 
 
@@ -374,6 +764,18 @@ data <- Start(repos = mask_path,
   length(attr(data, 'Variables')$common$nav_lat),
   8
   )
+#NOTE: The following code doesn't work, and since this netCDF file doesn't follow the convention (one var per file), we leave this development to the future.
+#suppressWarnings(
+#dataF <- Start(repos = mask_path,
+#              var = c('nav_lon', 'nav_lat'),
+#              t = 'first',
+#              z = 'first',
+#              x = 'all',
+#              y = 'all',
+#              return_vars = list(var_names = NULL),
+#              var_var = 'var_names',
+#              retrieve = F)
+#)
 
 })
 
@@ -396,7 +798,7 @@ suppressWarnings(
                 return_vars = list(time = 'sdate',
                                    lon = 'dat',
                                    lat = 'dat'),
-                metadata_dims = 'dat',  # it can be omitted since it is automatically specified as 'dat'
+                metadata_dims = 'dat', 
                 retrieve = T
                )
 )
@@ -409,6 +811,10 @@ suppressWarnings(
   list(lon = NULL, lat = NULL)
   )
   expect_equal(
+  names(attr(data, "Variables")$system5_m1),
+  c('lon', 'lat', 'tas')
+  )
+  expect_equal(
   length(attr(data, "Variables")$system5_m1$lon),
   1296
   )
@@ -421,6 +827,49 @@ suppressWarnings(
   array(c(NA, "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20170101.nc"),
         dim = c(dat = 2, var = 1, sdate = 1))
   )
+
+suppressWarnings(
+ dataF <- Start(dat = list(list(name = 'system4_m1', path = repos2),
+                           list(name = 'system5_m1', path = repos)),
+                var = var,
+                sdate = '20170101',
+                ensemble = indices(1),
+                time = indices(1),
+                lat = 'all',
+                lon = 'all',
+                synonims = list(lat = c('lat', 'latitude'),
+                                lon = c('lon', 'longitude')),
+                return_vars = list(time = 'sdate',
+                                   lon = 'dat',
+                                   lat = 'dat'),
+                metadata_dims = 'dat', 
+                retrieve = F
+               )
+)
+
+  expect_equal(
+  attr(dataF, "Variables")$system4_m1,
+  list(lon = NULL, lat = NULL)
+  )
+  expect_equal(
+  length(attr(dataF, "Variables")$system5_m1$lon),
+  1296
+  )
+  expect_equal(
+  names(attr(dataF, "Variables")$system5_m1),
+  c('lon', 'lat', 'tas')
+  )
+  expect_equal(
+  length(attr(dataF, "Variables")$system5_m1$tas),
+  22
+  )
+  expect_equal(
+  attr(dataF, 'ExpectedFiles'),
+  array(c("/esarchive/exp/ecmwf/system4_m1/monthly_mean/tas_f2h/tas_20170101.nc", "/esarchive/exp/ecmwf/system5_m1/monthly_mean/tas_f6h/tas_20170101.nc"),
+        dim = c(dat = 2, var = 1, sdate = 1))
+  )
+
+
 })
 
 test_that("8. Two data sets, both have files but the first file is missing", {
@@ -457,10 +906,6 @@ data <- Start(dataset = path_list,
   5500
   )
   expect_equal(
-  length(attr(data, "Variables")$MPI_ESM),
-  3
-  )
-  expect_equal(
   length(attr(data, "Variables")$MPI_ESM$lon),
   30
   )
@@ -489,4 +934,43 @@ data <- Start(dataset = path_list,
   )
 
 
+suppressWarnings(
+dataF <- Start(dataset = path_list,
+              var = 'tasmin',
+              member = list(c('r1i1p1f1', 'r2i1p1f2')),
+              sdate = paste0(2018),
+              chunk = list(c('20181101-20281231', '20181101-20181230')),
+              time = indices(1), #'all',
+              lat = values(list(0, 14)),
+              lon = values(list(0, 28)),
+              synonims = list(time = c('fmonth','time'),
+                              lon = c('lon', 'longitude'),
+                              lat = c('lat', 'latitude')),
+              return_vars = list(lat = 'dataset', lon = 'dataset'),
+              lat_reorder = Sort(),
+              num_procs = 1,
+              retrieve = F)
+)
+
+  expect_equal(
+  length(attr(dataF, "Variables")$MPI_ESM$lon),
+  30
+  )
+  expect_equal(
+  names(attr(dataF, "Variables")$MPI_ESM),
+  c('lat', 'lon', 'tasmin')
+  )
+  expect_equal(
+  length(attr(dataF, "Variables")$MPI_ESM$tasmin),
+  22
+  )
+  expect_equal(
+  names(attr(dataF, "Variables")$HadGEM3),
+  c('lat', 'lon')
+  )
+  expect_equal(
+  length(attr(dataF, "Variables")$HadGEM3$lon),
+  34
+  )
+
 })
diff --git a/tests/testthat/test-Start-metadata_reshaping.R b/tests/testthat/test-Start-metadata_reshaping.R
index ac332cdc95e0608fa320706bf26e8ad4bddfcc5d..7e9c2801dcf7e659fd06741504122574c80d16a3 100644
--- a/tests/testthat/test-Start-metadata_reshaping.R
+++ b/tests/testthat/test-Start-metadata_reshaping.R
@@ -651,7 +651,7 @@ data <- Start(dat = paste0('/esarchive/recon/ecmwf/erainterim/6hourly/',
 )
 expect_equal(
 names(attributes(data)$Variables$common),
-NULL
+"tas"
 )
 
 })