Start.R

                    }
                    
                    # Sorting the transformed variable and working out the indices again after transform.
                    if (!is.null(dim_reorder_params[[inner_dim]])) {
                      transformed_subset_var_reorder <- dim_reorder_params[[inner_dim]](transformed_subset_var)
                      transformed_subset_var <- transformed_subset_var_reorder$x
                      #NOTE: The fix here solves the mis-ordered lon when across_meridian. 
                      transformed_subset_var_unorder <- transformed_subset_var_reorder$ix
                      #                      transformed_subset_var_unorder <- sort(transformed_subset_var_reorder$ix, index.return = TRUE)$ix
                    } else {
                      transformed_subset_var_unorder <- 1:length(transformed_subset_var)
                    }
                    if (!selectors_are_indices) {  # selectors are values
                      sub_array_of_sri <- selector_checker(
                        sub_array_of_selectors, transformed_subset_var,
                        tolerance = if (aiat) {
                          tolerance_params[[inner_dim]]
                          } else {
                            NULL
                          })
                      if (!is.list(sub_array_of_sri)) {
                        sub_array_of_sri <- unique(sub_array_of_sri)
                      }
                    } else {  # selectors are indices
                      # Need to transfer to values first, then use the values to get the new
                      # indices in transformed_subset_var.
                      if (is.list(sub_array_of_selectors)) {
                        ori_values <- vars_to_transform[[var_with_selectors_name]][sub_array_of_selectors[[1]]:sub_array_of_selectors[[2]]]
                      } else {
                        ori_values <- vars_to_transform[[var_with_selectors_name]][sub_array_of_selectors]
                      }
                      sub_array_of_sri <- selector_checker(
                        ori_values, transformed_subset_var,
                        tolerance = if (aiat) {
                          tolerance_params[[inner_dim]]
                          } else {
                            NULL
                          })
                      # Here may need to further modify considering aiat. If aiat = FALSE,
                      # (i.e., indices are taken before transform), unique() is needed.
                      sub_array_of_sri <- unique(sub_array_of_sri)
                    }
                    
                    # Check if selectors fall out of the range of the transform grid
                    # It may happen when original lon is [-180, 180] while want to regrid to
                    # [0, 360], and lon selector = [-20, -10].
                    if (any(is.na(sub_array_of_sri))) {
                      stop(paste0("The selectors of ",
                                  inner_dim, " are out of range of transform grid '",
                                  transform_params$grid, "'. Use parameter '",
                                  inner_dim, "_reorder' or change ", inner_dim,
                                  " selectors."))
                    }
                    
                    if (goes_across_prime_meridian) {
                      
                      if (sub_array_of_sri[[1]] == sub_array_of_sri[[2]]) {
                        # global longitude
                        sub_array_of_sri <- c(1:length(transformed_subset_var))
                      } else {
                        # the common case, i.e., non-global
#                        # NOTE: Because sub_array_of_sri order is exchanged due to 
#                        # previous development, here [[1]] and [[2]] should exchange
#                        sub_array_of_sri <- c(1:sub_array_of_sri[[1]],
#                                              sub_array_of_sri[[2]]:length(transformed_subset_var))
                        #NOTE: the old code above is not suitable for all the possible cases.
                        #      If sub_array_of_selectors is not exactly the value in transformed_subset_var, sub_array_of_sri[[1]] will be larger than sub_array_of_sri[[2]].
                        #      Though here is not global case, we already have transformed_subset_var cropped as the desired region, so it is okay to use the whole length. Not sure if it will cause other problems...
                        sub_array_of_sri <- 1:length(transformed_subset_var)
                      }
                      
                    } else if (is.list(sub_array_of_sri)) {
                      sub_array_of_sri <- sub_array_of_sri[[1]]:sub_array_of_sri[[2]]
                    }

#========================================================

# Instead of using values to find sri, directly use the destination grid to count.
#NOTE: sub_array_of_sri seems to start at 1 always (because crop = c(lonmin, lonmax, latmin, latmax) already?)
                    if (chunks[[inner_dim]]["n_chunks"] > 1) {
                      sub_array_of_sri <- sub_array_of_sri[get_chunk_indices(
                        length(sub_array_of_sri),
                        chunks[[inner_dim]]["chunk"],
                        chunks[[inner_dim]]["n_chunks"],
                        inner_dim)]
                    }
#========================================================

                    ordered_sri <- sub_array_of_sri
                    sub_array_of_sri <- transformed_subset_var_unorder[sub_array_of_sri]

###########################old##################################
#                    if (chunks[[inner_dim]]["n_chunks"] > 1) {
#                      tmp <- which(transformed_subset_var >= min(sub_sub_array_of_values) &
#                                   transformed_subset_var <= max(sub_sub_array_of_values))
#                      sub_array_of_sri <- sub_array_of_sri[tmp]
#                    }
################################################################

                    # In this case, the tvi are not defined and the 'transformed_subset_var'
                    # will be taken instead of the var transformed before in the code.
                    if (debug) {
                      if (inner_dim %in% dims_to_check) {
                        print("-> FIRST INDEX:")
#                        print(first_index)
                        print("NOTE: Check function generate_sub_array_of_fri() in zzz.R")
                        print("-> LAST INDEX:")
#                        print(last_index)
                        print("NOTE: Check function generate_sub_array_of_fri() in zzz.R")
                        print("-> STRUCTURE OF FIRST ROUND INDICES:")
                        print(str(sub_array_of_fri))
                        print("-> STRUCTURE OF SECOND ROUND INDICES:")
                        print(str(sub_array_of_sri))
                        print("-> STRUCTURE OF TRANSFORMED VARIABLE INDICES:")
                        print(str(tvi))
                      }
                    }
                    ###                    # If the selectors are expressed after transformation
                    ###                    } else {
                    ###if (debug) {
                    ###if (inner_dim %in% dims_to_check) {
                    ###print("-> SELECTORS REQUESTED AFTER TRANSFORM.")
                    ###}
                    ###}
                    ###                      if (goes_across_prime_meridian) {
                    ###                        sub_array_of_indices <- c(sub_array_of_indices[[1]]:m,
                    ###                                                    1:sub_array_of_indices[[2]])
                    ###                      }
                    ###                      first_index <- min(unlist(sub_array_of_indices))
                    ###                      last_index <- max(unlist(sub_array_of_indices))
                    ###                      first_index_before_transform <- max(transform_indices(first_index, m, n) - beta, 1)
                    ###                      last_index_before_transform <- min(transform_indices(last_index, m, n) + beta, n)
                    ###                      sub_array_of_fri <- first_index_before_transform:last_index_before_transform
                    ###                      n_of_extra_cells <- round(beta / n * m)
                    ###                      if (is.list(sub_array_of_indices) && (length(sub_array_of_indices) > 1)) {
                    ###                        sub_array_of_sri <- 1:(last_index - first_index + 1) 
                    ###                        if (is.null(tvi)) {
                    ###                          tvi <- sub_array_of_sri + first_index - 1
                    ###                        }
                    ###                      } else {
                    ###                        sub_array_of_sri <- sub_array_of_indices - first_index + 1
                    ###                        if (is.null(tvi)) {
                    ###                          tvi <- sub_array_of_indices
                    ###                        }
                    ###                      }
                    ###                      sub_array_of_sri <- sub_array_of_sri + n_of_extra_cells
                    sri <- do.call('[[<-', c(list(x = sri), as.list(selector_store_position),
                                             list(value = sub_array_of_sri)))

                  } else {  # !with_transform
                    sub_array_of_fri <- generate_sub_array_of_fri(
                      with_transform, goes_across_prime_meridian, sub_array_of_indices, n, beta,
                      is_circular_dim)
                  }

                  # Reorder sub_array_of_fri if reordering function is used.
                  # It was index in the assigned order (e.g., in [-180, 180] if CircularSort(-180, 180)), and here is changed to the index in the original order.
                  if (!is.null(var_unorder_indices)) {
                    if (is.null(ordered_fri)) {
                      ordered_fri <- sub_array_of_fri
                    }
                    sub_array_of_fri <- var_unorder_indices[sub_array_of_fri]
                  }
                  fri <- do.call('[[<-', c(list(x = fri), as.list(selector_store_position),
                                           list(value = sub_array_of_fri)))

                  #NOTE: This part existed always but never was used. taken_chunks
                  #      is related to merge_across_dims, but I don't know how it is
                  #      used (maybe for higher efficiency?)
#                  if (!is.null(crossed_file_dim)) {
#                    taken_chunks[selector_store_position[[crossed_file_dim]]] <- TRUE
#                  } else {
                  taken_chunks <- TRUE
#                  }
                }
              } else {
              # The inner dim goes across a file dim (e.g., time_across = 'sdate')
                if (debug) {
                  if (inner_dim %in% dims_to_check) {
                    print("-> THE INNER DIMENSION GOES ACROSS A FILE DIMENSION.")
                  }
                }
                # If "<inner_dim>_across = <crossed_file_dim> + merge_across_dims = FALSE + chunk over <inner_dim>", return error because this instance is not logically correct.
                if (chunks[[inner_dim]]["n_chunks"] > 1 & inner_dim %in% inner_dims_across_files) {
                  stop("Chunk over dimension '", inner_dim, "' is not allowed because '",
                       inner_dim, "' is across '",
                       names(inner_dims_across_files)[which(inner_dim %in% inner_dims_across_files)], "'.")
                }

                if (inner_dim %in% names(dim(sub_array_of_selectors))) {
                  if (is.null(var_with_selectors_name)) {
                    if (!largest_dims_length | (largest_dims_length & length(unique(inner_dim_lengths)) <= 1)) {  #old code
                      maximal_indice <- data_dims[inner_dim] * chunk_amount
                    } else { # files have different length of inner dim
                      maximal_indice <- sum(inner_dim_lengths)
                    }

                    if (any(na.omit(unlist(sub_array_of_selectors)) < 1) ||
                        any(na.omit(unlist(sub_array_of_selectors)) > maximal_indice)) {
                      stop("Provided indices out of range for dimension '", inner_dim, "' ",
                           "for dataset '", dat[[i]][['name']], "' (accepted range: 1 to ",
                           maximal_indice, ").")
                    }
                  } else {
                    if (inner_dim %in% names(dim(sub_array_of_values))) {
                      # NOTE: Put across-inner-dim at the 1st position.
                      # POSSIBLE PROB!! Only organize inner dim, the rest dims may not in the same order as sub_array_of_selectors below.
                      inner_dim_pos_in_sub_array <- which(names(dim(sub_array_of_values)) == inner_dim)
                      if (inner_dim_pos_in_sub_array != 1) {
                        new_sub_array_order <- (1:length(dim(sub_array_of_values)))[-inner_dim_pos_in_sub_array]
                        new_sub_array_order <- c(inner_dim_pos_in_sub_array, new_sub_array_order)
                        sub_array_of_values <- .aperm2(sub_array_of_values, new_sub_array_order)
                      }
                    }
                  }
                  
                  # NOTE: Put across-inner-dim at the 1st position.
                  # POSSIBLE PROB!! Only organize inner dim, the rest dims may not in the same order as sub_array_of_values above.
                  inner_dim_pos_in_sub_array <- which(names(dim(sub_array_of_selectors)) == inner_dim)
                  if (inner_dim_pos_in_sub_array != 1) {
                    new_sub_array_order <- (1:length(dim(sub_array_of_selectors)))[-inner_dim_pos_in_sub_array]
                    new_sub_array_order <- c(inner_dim_pos_in_sub_array, new_sub_array_order)
                    sub_array_of_selectors <- .aperm2(sub_array_of_selectors, new_sub_array_order)
                  }
                  sub_array_of_indices <- selector_checker(sub_array_of_selectors, sub_array_of_values,
                                                           tolerance = tolerance_params[[inner_dim]])
                  # It is needed to expand the indices here, otherwise for 
                  # values(list(date1, date2)) only 2 values are picked.
                  if (is.list(sub_array_of_indices)) {
                    sub_array_of_indices <- sub_array_of_indices[[1]]:sub_array_of_indices[[2]]
                  }
                  sub_array_of_indices <- sub_array_of_indices[get_chunk_indices(length(sub_array_of_indices),
                                                                             chunks[[inner_dim]]['chunk'],
                                                                             chunks[[inner_dim]]['n_chunks'],
                                                                             inner_dim)]
                  sub_array_is_list <- FALSE
                  if (is.list(sub_array_of_indices)) {
                    sub_array_is_list <- TRUE
                    sub_array_of_indices <- unlist(sub_array_of_indices)
                  }

                  # "indices_chunk" refers to in which file the 
                  # sub_array_of_indices is; "transformed_indices" 
                  # refers to the indices of sub_array_of_indices in each file.
                  if (!largest_dims_length | 
                      (largest_dims_length & length(unique(inner_dim_lengths)) <= 1)) {
                    # old code; all the files have the same length of inner_dim
                    if (is.null(var_with_selectors_name)) {
                      indices_chunk <- floor((sub_array_of_indices - 1) / data_dims[inner_dim]) + 1
                      transformed_indices <- ((sub_array_of_indices - 1) %% data_dims[inner_dim]) + 1
                    } else {
                      indices_chunk <- floor((sub_array_of_indices - 1) / var_full_dims[inner_dim]) + 1
                      transformed_indices <- ((sub_array_of_indices - 1) %% var_full_dims[inner_dim]) + 1
                    }
                  } else {  # files have different inner dim length
                    indices_chunk <- c()
                    for (item in 1:length(inner_dim_lengths)) {
                      tmp <- which(sub_array_of_indices <= cumsum(inner_dim_lengths)[item])
                      indices_chunk <- c(indices_chunk, rep(item, length(tmp) - length(indices_chunk)))
                    }
                    sub_array_of_indices_by_file <- split(sub_array_of_indices, indices_chunk)
                    for (item in 2:length(sub_array_of_indices_by_file)) {
                      sub_array_of_indices_by_file[[item]] <- sub_array_of_indices_by_file[[item]] - cumsum(inner_dim_lengths)[item - 1]
                    }
                    transformed_indices <- unlist(sub_array_of_indices_by_file, use.names = FALSE)
                  }

                  if (sub_array_is_list) {
                    sub_array_of_indices <- as.list(sub_array_of_indices)
                  }
                  if (debug) {
                    if (inner_dim %in% dims_to_check) {
                      print("-> GOING TO ITERATE ALONG CHUNKS.")
                    }
                  }

                  for (chunk in 1:chunk_amount) {
                    if (!is.null(names(selector_store_position))) {
                      selector_store_position[crossed_file_dim] <- chunk
                    } else {
                      selector_store_position <- chunk
                    }
                    sub_array_of_indices <- transformed_indices[which(indices_chunk == chunk)]

                    #NOTE: This 'with_transform' part is probably not tested because 
                    #      here is for the inner dim that goes across a file dim, which
                    #      is normally not lat and lon dimension. If in the future, we 
                    #      can interpolate time, this part needs to be examined.
                    if (with_transform) {
                      # If the provided selectors are expressed in the world
                      # before transformation
                      if (!aiat) {
                        first_index <- min(unlist(sub_array_of_indices))
                        last_index <- max(unlist(sub_array_of_indices))
                        sub_array_of_fri <- max(c(first_index - beta, 1)):min(c(last_index + beta, n))
                        sub_array_of_sri <- transform_indices(unlist(sub_array_of_indices) - first_index + 1, n, m)
                        if (is.list(sub_array_of_indices)) {
                          if (length(sub_array_of_sri) > 1) {
                            sub_array_of_sri <- sub_array_of_sri[[1]]:sub_array_of_sri[[2]]
                          }
                        }
                        ##TODO: TRANSFORM SUBSET VARIABLE AS ABOVE, TO COMPUTE SRI
                        # If the selectors are expressed after transformation
                      } else {
                        first_index <- min(unlist(sub_array_of_indices))
                        last_index <- max(unlist(sub_array_of_indices))
                        first_index_before_transform <- max(transform_indices(first_index, m, n) - beta, 1)
                        last_index_before_transform <- min(transform_indices(last_index, m, n) + beta, n)
                        sub_array_of_fri <- first_index_before_transform:last_index_before_transform
                        if (is.list(sub_array_of_indices) && (length(sub_array_of_indices) > 1)) {
                          sub_array_of_sri <- 1:(last_index - first_index + 1) + 
                            round(beta / n * m) 
                        } else {
                          sub_array_of_sri <- sub_array_of_indices - first_index + 1 +
                            round(beta / n * m)
                        }
                        ##TODO: FILL IN TVI
                      }
                      sri <- do.call('[[<-', c(list(x = sri), as.list(selector_store_position),
                                               list(value = sub_array_of_sri)))
                      if (length(sub_array_of_sri) > 0) {
                        taken_chunks[chunk] <- TRUE
                      }
                    } else {
                      sub_array_of_fri <- sub_array_of_indices
                      if (length(sub_array_of_fri) > 0) {
                        taken_chunks[chunk] <- TRUE
                      }
                    }

                    if (!is.null(var_unorder_indices)) {
                      ordered_fri <- sub_array_of_fri
                      sub_array_of_fri <- var_unorder_indices[sub_array_of_fri]
                    }
                    fri <- do.call('[[<-', c(list(x = fri), as.list(selector_store_position),
                                             list(value = sub_array_of_fri)))
                  }
                  if (debug) {
                    if (inner_dim %in% dims_to_check) {
                      print("-> FINISHED ITERATING ALONG CHUNKS")
                    }
                  }
                } else {
                  stop("Provided array of indices for dimension '", inner_dim, "', ",
                       "which goes across the file dimension '", crossed_file_dim, "', but ",
                       "the provided array does not have the dimension '", inner_dim, 
                       "', which is mandatory.")
                }
              }
            }
          }
          if (debug) {
            if (inner_dim %in% dims_to_check) {
              print("-> PROCEEDING TO CROP VARIABLES")
            }
          }
          #if ((length(selector_array) == 1) && (selector_array %in% c('all', 'first', 'last'))) {
          #if (!is.null(var_with_selectors_name) || (is.null(var_with_selectors_name) && is.character(selector_array) &&
          #    (length(selector_array) == 1) && (selector_array %in% c('all', 'first', 'last')))) {
          empty_chunks <- which(!taken_chunks)
          if (length(empty_chunks) >= length(taken_chunks)) {
            stop("Selectors do not match any of the possible values for the dimension '", inner_dim, "'.")
          }
          if (length(empty_chunks) > 0) {
            #                # Get the first group of chunks to remove, and remove them. 
            #                # E.g., from c(1, 2, 4, 5, 6, 8, 9) remove only 1 and 2
            #                dist <- abs(rev(empty_chunks) - c(rev(empty_chunks)[1] - 1, head(rev(empty_chunks), length(rev(empty_chunks)) - 1)))
            #                if (all(dist == 1)) {
            #                  start_chunks_to_remove <- NULL
            #                } else {
            #                  first_chunk_to_remove <- tail(which(dist > 1), 1)
            #                  start_chunks_to_remove <- rev(rev(empty_chunks)[first_chunk_to_remove:length(empty_chunks)])
            #                }
            #                # Get the last group of chunks to remove, and remove them. 
            #                # E.g., from c(1, 2, 4, 5, 6, 8, 9) remove only 8 and 9
            #                dist <- abs(empty_chunks - c(empty_chunks[1] - 1, head(empty_chunks, length(empty_chunks) - 1)))
            #                if (all(dist == 1)) {
            #                  first_chunk_to_remove <- 1
            #                } else {
            #                  first_chunk_to_remove <- tail(which(dist > 1), 1)
            #                }
            #                end_chunks_to_remove <- empty_chunks[first_chunk_to_remove:length(empty_chunks)]
            #                chunks_to_keep <- which(!((1:length(taken_chunks)) %in% c(start_chunks_to_remove, end_chunks_to_remove)))
            chunks_to_keep <- which(taken_chunks)
            dims_to_crop[[crossed_file_dim]] <- c(dims_to_crop[[crossed_file_dim]], list(chunks_to_keep))
            #                found_indices <- Subset(found_indices, crossed_file_dim, chunks_to_keep)
            #                # Crop dataset variables file dims.
            #                for (picked_var in names(picked_vars[[i]])) {
            #                  if (crossed_file_dim %in% names(dim(picked_vars[[i]][[picked_var]]))) {
            #                    picked_vars[[i]][[picked_var]] <- Subset(picked_vars[[i]][[picked_var]], crossed_file_dim, chunks_to_keep)
            #                  }
            #                }
          }
          #}
          dat[[i]][['selectors']][[inner_dim]] <- list(fri = fri, sri = sri)
          # Crop dataset variables inner dims.
          # Crop common variables inner dims.
          types_of_var_to_crop <- 'picked'
          if (with_transform) {
            types_of_var_to_crop <- c(types_of_var_to_crop, 'transformed')
          }
          if (!is.null(dim_reorder_params[[inner_dim]])) {
            types_of_var_to_crop <- c(types_of_var_to_crop, 'reordered')
          }
          for (type_of_var_to_crop in types_of_var_to_crop) {
            if (type_of_var_to_crop == 'transformed') {
              if (is.null(tvi)) {
                if (!is.null(dim_reorder_params[[inner_dim]])) {
                  crop_indices <- unique(unlist(ordered_sri))
                } else {
                  crop_indices <- unique(unlist(sri))
                }
              } else {
                crop_indices <- unique(unlist(tvi))
              }
              vars_to_crop <- transformed_vars[[i]]
              common_vars_to_crop <- transformed_common_vars
            } else if (type_of_var_to_crop == 'reordered') {
              crop_indices <- unique(unlist(ordered_fri))
              vars_to_crop <- picked_vars_ordered[[i]]
              common_vars_to_crop <- picked_common_vars_ordered
            } else {
              #TODO: If fri has different indices in each list, the crop_indices should be 
              #      separated for each list. Otherwise, picked_common_vars later will be wrong.
              crop_indices <- unique(unlist(fri))
              vars_to_crop <- picked_vars[[i]]
              common_vars_to_crop <- picked_common_vars
            }
            for (var_to_crop in names(vars_to_crop)) {
              if (inner_dim %in% names(dim(vars_to_crop[[var_to_crop]]))) {
                if (!is.null(crop_indices)) {
                  if (type_of_var_to_crop == 'transformed') {
                    if (!aiat) {
                      if (!(length(selector_array) == 1 &
                            all(selector_array %in% c('all', 'first', 'last')))) {
                        vars_to_crop[[var_to_crop]] <- Subset(transformed_subset_var, inner_dim, crop_indices)
                      } else {
                        vars_to_crop[[var_to_crop]] <-
                          Subset(transformed_var_with_selectors, inner_dim, crop_indices)
                      }
                    } else {
                      vars_to_crop[[var_to_crop]] <- Subset(vars_to_crop[[var_to_crop]], inner_dim, crop_indices)
                    }
                  } else {
                    vars_to_crop[[var_to_crop]] <- Subset(vars_to_crop[[var_to_crop]], inner_dim, crop_indices)
                  }
                }
              }
            }
            if (i == length(dat)) {
              for (common_var_to_crop in names(common_vars_to_crop)) {
                if (inner_dim %in% names(dim(common_vars_to_crop[[common_var_to_crop]]))) {

                  if (type_of_var_to_crop == 'transformed' & !aiat) {
                    if (!(length(selector_array) == 1 &
                          all(selector_array %in% c('all', 'first', 'last')))) {
                      common_vars_to_crop[[common_var_to_crop]] <- 
                        Subset(transformed_subset_var, inner_dim, crop_indices)
                    } else {
                      common_vars_to_crop[[common_var_to_crop]] <-
                        Subset(transformed_var_with_selectors, inner_dim, crop_indices)
                    }
                  } else {
                    if (!is.null(crossed_file_dim)) {  #merge_across_dims, crossed_file_dim is the depended file dim
                      #NOTE: When is not this case??? Maybe this condition is not needed
                      if (any(crossed_file_dim %in% names(dim(common_vars_to_crop[[common_var_to_crop]])))) {
                        tmp <- common_vars_to_crop[[common_var_to_crop]]
                        tmp_attributes <- attributes(common_vars_to_crop[[common_var_to_crop]])
                        dim_extra_ind <- which(!names(dim(tmp)) %in% c(crossed_file_dim, inner_dim))
                        if (!identical(dim_extra_ind, integer(0))) {
                          tmp_list <- asplit(tmp, dim_extra_ind)
                          dim_file_ind <- which(names(dim(tmp_list[[1]])) %in% crossed_file_dim)
                          tmp_list <- lapply(tmp_list, asplit, dim_file_ind)
                        } else {  # only crossed_file_dim and inner_dim
                          dim_file_ind <- which(names(dim(tmp)) %in% crossed_file_dim)
                          tmp_list <- asplit(tmp, dim_file_ind)
                          # Add another layer to be consistent with the first case above
                          tmp_list <- list(tmp_list)
                        }
                        max_fri_length <- max(sapply(fri, length))
                        for (i_extra_dim in 1:length(tmp_list)) {
                          for (i_fri in 1:length(fri)) {
                            tmp_list[[i_extra_dim]][[i_fri]] <-
                              tmp_list[[i_extra_dim]][[i_fri]][fri[[i_fri]]]
  
                            if (length(tmp_list[[i_extra_dim]][[i_fri]]) != max_fri_length) {
                              tmp_list[[i_extra_dim]][[i_fri]] <-
                                c(tmp_list[[i_extra_dim]][[i_fri]], rep(NA, max_fri_length - length(tmp_list[[i_extra_dim]][[i_fri]])))
                            }
                          }
                        }
                        # Change list back to array
                        tmp_new_dim <- c(max_fri_length, dim(tmp)[crossed_file_dim], dim(tmp)[dim_extra_ind])
                        names(tmp_new_dim) <- c(inner_dim, crossed_file_dim, names(dim(tmp))[dim_extra_ind])
                        common_vars_to_crop[[common_var_to_crop]] <-
                          array(unlist(tmp_list), dim = tmp_new_dim)
                        # Reorder back
                        common_vars_to_crop[[common_var_to_crop]] <-
                          aperm(common_vars_to_crop[[common_var_to_crop]], match(names(dim(tmp)), names(tmp_new_dim)))
                        # Put attributes back
                        tmp <- which(!names(tmp_attributes) %in% names(attributes(common_vars_to_crop[[common_var_to_crop]])))
                        attributes(common_vars_to_crop[[common_var_to_crop]]) <- 
                          c(attributes(common_vars_to_crop[[common_var_to_crop]]), 
                            tmp_attributes[tmp])

                        if ('time' %in% synonims[[common_var_to_crop]]) {
                          # Convert number back to time
                          common_vars_to_crop[[common_var_to_crop]] <-
                            as.POSIXct(common_vars_to_crop[[common_var_to_crop]],
                                       origin = "1970-01-01", tz = 'UTC')
                        }
                      }
                    } else {  # old code

                      common_vars_to_crop[[common_var_to_crop]] <- Subset(common_vars_to_crop[[common_var_to_crop]], inner_dim, crop_indices)
                    }

                  }

                }
              }
            }
            if (type_of_var_to_crop == 'transformed') {
              if (!is.null(vars_to_crop)) {
                transformed_vars[[i]] <- vars_to_crop
              }
              if (i == length(dat)) {
                transformed_common_vars <- common_vars_to_crop
              }
            } else if (type_of_var_to_crop == 'reordered') {
              if (!is.null(vars_to_crop)) {
                picked_vars_ordered[[i]] <- vars_to_crop
              }
              if (i == length(dat)) {
                picked_common_vars_ordered <- common_vars_to_crop
              }
            } else {
              if (!is.null(vars_to_crop)) {
                picked_vars[[i]] <- vars_to_crop
              }
              if (i == length(dat)) {
                #NOTE: To avoid redundant run
                if (inner_dim %in% names(common_vars_to_crop)) {
                  picked_common_vars <- common_vars_to_crop
                }
              }
            }
          }
          #}
        }
        # After the selectors have been picked (using the original variables), 
        # the variables are transformed. At that point, the original selectors
        # for the transformed variables are also kept in the variable original_selectors.
        #print("L")
      }
    }
  }
  #  if (!is.null(transformed_common_vars)) {
  #    picked_common_vars[names(transformed_common_vars)] <- transformed_common_vars
  #  }
  # Remove the trailing chunks, if any.
  for (file_dim in names(dims_to_crop)) {
    #    indices_to_keep <- min(sapply(dims_to_crop[[file_dim]], min)):max(sapply(dims_to_crop[[file_dim]], max))
    ## TODO: Merge indices in dims_to_crop with some advanced mechanism?
    indices_to_keep <- unique(unlist(dims_to_crop[[file_dim]]))
    array_of_files_to_load <- Subset(array_of_files_to_load, file_dim, indices_to_keep)
    array_of_not_found_files <- Subset(array_of_not_found_files, file_dim, indices_to_keep)
    for (i in 1:length(dat)) {
      # Crop selectors
      for (selector_dim in names(dat[[i]][['selectors']])) {
        if (selector_dim == file_dim) {
          for (j in 1:length(dat[[i]][['selectors']][[selector_dim]][['fri']])) {
            dat[[i]][['selectors']][[selector_dim]][['fri']][[j]] <- dat[[i]][['selectors']][[selector_dim]][['fri']][[j]][indices_to_keep]
          }
          for (j in 1:length(dat[[i]][['selectors']][[selector_dim]][['sri']])) {
            dat[[i]][['selectors']][[selector_dim]][['sri']][[j]] <- dat[[i]][['selectors']][[selector_dim]][['sri']][[j]][indices_to_keep]
          }
        }
        if (file_dim %in% names(dim(dat[[i]][['selectors']][[selector_dim]][['fri']]))) {
          dat[[i]][['selectors']][[selector_dim]][['fri']] <- Subset(dat[[i]][['selectors']][[selector_dim]][['fri']], file_dim, indices_to_keep)
          dat[[i]][['selectors']][[selector_dim]][['sri']] <- Subset(dat[[i]][['selectors']][[selector_dim]][['sri']], file_dim, indices_to_keep)
        }
      }
      # Crop dataset variables file dims.
      for (picked_var in names(picked_vars[[i]])) {
        if (file_dim %in% names(dim(picked_vars[[i]][[picked_var]]))) {
          picked_vars[[i]][[picked_var]] <- Subset(picked_vars[[i]][[picked_var]], file_dim, indices_to_keep)
        }
      }
      for (transformed_var in names(transformed_vars[[i]])) {
        if (file_dim %in% names(dim(transformed_vars[[i]][[transformed_var]]))) {
          transformed_vars[[i]][[transformed_var]] <- Subset(transformed_vars[[i]][[transformed_var]], file_dim, indices_to_keep)
        }
      }
    }
    # Crop common variables file dims.
    for (picked_common_var in names(picked_common_vars)) {
      if (file_dim %in% names(dim(picked_common_vars[[picked_common_var]]))) {
        picked_common_vars[[picked_common_var]] <- Subset(picked_common_vars[[picked_common_var]], file_dim, indices_to_keep)
      }
    }
    for (transformed_common_var in names(transformed_common_vars)) {
      if (file_dim %in% names(dim(transformed_common_vars[[transformed_common_var]]))) {
        transformed_common_vars[[transformed_common_var]] <- Subset(transformed_common_vars[[transformed_common_var]], file_dim, indices_to_keep)
      }
    }
  }
  # Calculate the size of the final array.
  total_inner_dims <- NULL
  for (i in 1:length(dat)) {
    if (dataset_has_files[i]) {
      inner_dims <- expected_inner_dims[[i]]
      inner_dims <- sapply(inner_dims, 
                           function(x) {
                             if (!all(sapply(dat[[i]][['selectors']][[x]][['sri']], is.null))) {
                               max(sapply(dat[[i]][['selectors']][[x]][['sri']], length))
                             } else {
                               if (length(var_params[[x]]) > 0) {
                                 if (var_params[[x]] %in% names(transformed_vars[[i]])) {
                                   length(transformed_vars[[i]][[var_params[[x]]]])
                                 } else if (var_params[[x]] %in% names(transformed_common_vars)) {
                                   length(transformed_common_vars[[var_params[[x]]]])
                                 } else {
                                   max(sapply(dat[[i]][['selectors']][[x]][['fri']], length))
                                 }
                               } else {
                                 max(sapply(dat[[i]][['selectors']][[x]][['fri']], length))
                               }
                             }
                           })
      names(inner_dims) <- expected_inner_dims[[i]]
      if (is.null(total_inner_dims)) {
        total_inner_dims <- inner_dims
      } else {
        new_dims <- .MergeArrayDims(total_inner_dims, inner_dims)
        total_inner_dims <- new_dims[[3]]
      }
    }
  }
  new_dims <- .MergeArrayDims(dim(array_of_files_to_load), total_inner_dims)
  final_dims <- new_dims[[3]][dim_names]
  # final_dims_fake is the vector of final dimensions after having merged the 
  # 'across' file dimensions with the respective 'across' inner dimensions, and
  # after having broken into multiple dimensions those dimensions for which 
  # multidimensional selectors have been provided.
  # final_dims will be used for collocation of data, whereas final_dims_fake 
  # will be used for shaping the final array to be returned to the user.
  final_dims_fake <- final_dims
  if (merge_across_dims) {
    final_dims_fake <- dims_merge(inner_dims_across_files, final_dims_fake)
  }
  #=========================================================================
  # Find the dimension to split if split_multiselected_dims = TRUE.
  # If there is no dimension able to be split, change split_multiselected_dims to FALSE.
  all_split_dims <- NULL
  inner_dim_has_split_dim <- NULL
  if (split_multiselected_dims) {
    tmp <- dims_split(dim_params, final_dims_fake)
    final_dims_fake <- tmp[[1]]
    # all_split_dims is a list containing all the split dims
    all_split_dims <- tmp[[2]]

   if (is.null(all_split_dims)) {
     split_multiselected_dims <- FALSE
     .warning(paste0("Not found any dimensions able to be split. The parameter ",
                     "'split_multiselected_dims' is changed to FALSE."))
   } else {
    tmp_fun <- function (x, y) {
      any(names(dim(x)) %in% y)
    }
    if (!is.null(picked_common_vars)) {
      inner_dim_has_split_dim <- names(which(unlist(lapply(
                                   picked_common_vars, tmp_fun, names(all_split_dims)))))
      if (!identical(inner_dim_has_split_dim, character(0))) {
        # If merge_across_dims also, it will be replaced later
        saved_reshaped_attr <- attr(picked_common_vars[[inner_dim_has_split_dim]], 'variables')
      }
    }
   }
  }
  #======================================================================
  # If only merge_across_dims and merge_across_dims_narm and no split_multiselected_dims,
  # the length of inner across dim (e.g., time) needs to be adjusted. Sum up the actual length
  # without potential NAs.
  if (merge_across_dims) {
    # Prepare the arguments for later use
    across_inner_dim <- inner_dims_across_files[[1]]  #TODO: more than one?
    # Get the length of each inner_dim ('time') along each file_dim ('file_date')  
    length_inner_across_dim <- lapply(dat[[i]][['selectors']][[across_inner_dim]][['fri']], length)
    dims_of_merge_dim <- dim(picked_common_vars[[across_inner_dim]])
    # Save attributes for later use. If split_multiselected_dims, this variable has been created above but is replaced here
    saved_reshaped_attr <- attr(picked_common_vars[[across_inner_dim]], 'variables')

    if (merge_across_dims_narm & !split_multiselected_dims) {
      final_dims_fake <- merge_narm_dims(final_dims_fake, across_inner_dim, length_inner_across_dim)
    }
  }
  
  if (!silent) {
    .message("Detected dimension sizes:")
    longest_dim_len <- max(sapply(names(final_dims_fake), nchar))
    longest_size_len <- max(sapply(paste0(final_dims_fake, ''), nchar))
    sapply(names(final_dims_fake), 
           function(x) {
             message(paste0("*   ", paste(rep(' ', longest_dim_len - nchar(x)), collapse = ''), 
                            x, ": ", paste(rep(' ', longest_size_len - nchar(paste0(final_dims_fake[x], ''))), collapse = ''), 
                            final_dims_fake[x]))
           })
    bytes <- prod(c(final_dims_fake, 8))
    dim_sizes <- paste(final_dims_fake, collapse = ' x ')
    if (retrieve) {
      .message(paste("Total size of requested data:"))
    } else {
      .message(paste("Total size of involved data:"))
    }
    .message(paste(dim_sizes, " x 8 bytes =", 
                   format(structure(bytes, class = "object_size"), units = "auto")), 
             indent = 2)
  }
  
  # NOTE: If split_multiselected_dims + merge_across_dims, the dim order may need to be changed.
  #       The inner_dim needs to be the first dim among split dims.
  # TODO: Cannot control the rest dims are in the same order or not...
  #       Suppose users put the same order of across inner and file dims.
  if (split_multiselected_dims & merge_across_dims) {
    # TODO: More than one split?
    inner_dim_pos_in_split_dims <- which(names(all_split_dims[[1]]) == inner_dims_across_files)  

    # if inner_dim is not the first, change!
    if (inner_dim_pos_in_split_dims != 1) {
      # Save the current final_dims_fake for reordering it back later
      final_dims_fake_output <- final_dims_fake
      tmp <- reorder_split_dims(all_split_dims[[1]], inner_dim_pos_in_split_dims, final_dims_fake)
      final_dims_fake <- tmp[[1]]
      all_split_dims[[1]] <- tmp[[2]]
    }
  }
  if (merge_across_dims | split_multiselected_dims) {
    if (!merge_across_dims & split_multiselected_dims & identical(inner_dim_has_split_dim, character(0))) {
      final_dims_fake_metadata <- NULL
    } else {
      if (!merge_across_dims & split_multiselected_dims & !is.null(picked_common_vars)) {
        if (any(names(all_split_dims[[1]]) %in% names(dim(picked_common_vars[[inner_dim_has_split_dim]]))) &
            names(all_split_dims)[1] != inner_dim_has_split_dim) {
          if (inner_dim_has_split_dim %in% names(final_dims)) {
            stop("Detect inner dimension in the split array, but merge_across_dims is not used. The output dimensions will be repeated. Check if the dimensions and parameters are correctly defined.")
          } else {
            # Only split no merge, time dim is not explicitly defined because the 
            # length is 1, the sdate dim to be split having 'time' as one dimension.
            # --> Take 'time' dim off from picked_common_vars.
            dim(picked_common_vars[[inner_dim_has_split_dim]]) <- dim(picked_common_vars[[inner_dim_has_split_dim]])[-which(names(dim(picked_common_vars[[inner_dim_has_split_dim]])) == inner_dim_has_split_dim)]
          }
        }
      }
      final_dims_fake_metadata <- find_final_dims_fake_metadata(
        merge_across_dims, split_multiselected_dims, picked_common_vars = picked_common_vars[[inner_dim_has_split_dim]], across_inner_dim,
        final_dims_fake, dims_of_merge_dim, all_split_dims)
    }
  }
  
  # The following several lines will only run if retrieve = TRUE
  if (retrieve) {
    
    ########## CREATING THE SHARED MATRIX AND DISPATCHING WORK PIECES ###########
    # TODO: try performance of storing all in cols instead of rows
    # Create the shared memory array, and a pointer to it, to be sent
    # to the work pieces.
    if (is.null(ObjectBigmemory)) {
        data_array <- bigmemory::big.matrix(nrow = prod(final_dims), ncol = 1)
    } else {
        data_array <- bigmemory::big.matrix(nrow = prod(final_dims), ncol = 1,
                                            backingfile = ObjectBigmemory,
                                            init = NA)
    }
    shared_matrix_pointer <- bigmemory::describe(data_array)
    if (is.null(ObjectBigmemory)) {
        name_bigmemory_obj <- attr(shared_matrix_pointer, 'description')$sharedName
    } else {
        name_bigmemory_obj <- attr(shared_matrix_pointer, 'description')$filename
    }

    #warning(paste("SharedName:", attr(shared_matrix_pointer, 'description')$sharedName))
    #warning(paste("Filename:", attr(shared_matrix_pointer, 'description')$filename))
    #if (!is.null(ObjectBigmemory)) {
    #  attr(shared_matrix_pointer, 'description')$sharedName <- ObjectBigmemory
    #}
    if (is.null(num_procs)) {
      num_procs <- future::availableCores()
    }
    # Creating a shared tmp folder to store metadata from each chunk
    array_of_metadata_flags <- array(FALSE, dim = dim(array_of_files_to_load))
    if (!is.null(metadata_dims)) {
      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])))))
    }
    metadata_file_counter <- 0
    metadata_folder <- tempfile('metadata')
    dir.create(metadata_folder)
    # Build the work pieces, each with:
    # - file path
    # - total size (dims) of store array
    # - start position in store array
    # - file selectors (to provide extra info. useful e.g. to select variable)
    # - indices to take from file
    work_pieces <- list()
    for (i in 1:length(dat)) {
      if (dataset_has_files[i]) {
        # metadata_file_counter may be changed by the following function
        work_pieces <- build_work_pieces(
                         work_pieces = work_pieces, i = i, selectors = dat[[i]][['selectors']], 
                         file_dims = found_file_dims[[i]],
                         inner_dims = expected_inner_dims[[i]], final_dims = final_dims,
                         found_pattern_dim = found_pattern_dim, 
                         inner_dims_across_files = inner_dims_across_files,
                         array_of_files_to_load = array_of_files_to_load,
                         array_of_not_found_files = array_of_not_found_files,
                         array_of_metadata_flags = array_of_metadata_flags,
                         metadata_file_counter = metadata_file_counter,
                         depending_file_dims = depending_file_dims, transform = transform,
                         transform_vars = transform_vars, picked_vars = picked_vars[[i]],
                         picked_vars_ordered = picked_vars_ordered[[i]],
                         picked_common_vars = picked_common_vars,
                         picked_common_vars_ordered = picked_common_vars_ordered, 
                         metadata_folder = metadata_folder, debug = debug)
      }
    }
    #print("N")
    if (debug) {
      print("-> WORK PIECES BUILT")
    }
    
    # Calculate the progress %s that will be displayed and assign them to 
    # the appropriate work pieces.
    work_pieces <- retrieve_progress_message(work_pieces, num_procs, silent)

    # NOTE: In .LoadDataFile(), metadata is saved in metadata_folder/1 or /2 etc. Before here,
    #       the path name is created in work_pieces but the path hasn't been built yet.
    if (num_procs == 1) {
      res3 <- .withWarnings(
        lapply(work_pieces, .LoadDataFile, 
                            shared_matrix_pointer = shared_matrix_pointer,
                            file_data_reader = file_data_reader, 
                            synonims = synonims,
                            transform = transform, 
                            transform_params = transform_params,
                            transform_crop_domain = transform_crop_domain,
                            silent = silent, debug = debug)
      )
      found_files <- res3$value

      if (!is.null(res3$warnings)) {
        if (is.null(warnings3)) {
          warnings3 <- res3$warnings
        } else {
          warnings3 <- c(warnings3, res3$warnings)
        }
      }
      
    } else {
      cluster <- parallel::makeCluster(num_procs, outfile = "")
      # Send the heavy work to the workers
      work_errors <- try({
        res3 <- .withWarnings(
        found_files <- parallel::clusterApplyLB(cluster, work_pieces, .LoadDataFile, 
                                                shared_matrix_pointer = shared_matrix_pointer,
                                                file_data_reader = file_data_reader,
                                                synonims = synonims,
                                                transform = transform, 
                                                transform_params = transform_params,
                                                transform_crop_domain = transform_crop_domain,
                                                silent = silent, debug = debug)
        )
        found_files <- res3$value

        if (!is.null(res3$warnings)) {
          if (is.null(warnings3)) {
            warnings3 <- res3$warnings
          } else {
            warnings3 <- c(warnings3, res3$warnings)
          }
        }

      })
      parallel::stopCluster(cluster)
    }
    
    if (!silent) {
      # if (progress_message != '')
      if (length(work_pieces) / num_procs >= 2 && !silent) {
        .message("\n", tag = '')
      }
    }
    #print("P")

    # If merge_across_dims = TRUE, there might be additional NAs due to unequal
    # inner_dim ('time') length across file_dim ('file_date').
    # If merge_across_dims_narm = TRUE, add additional lines to remove these NAs.
    # TODO: Now it assumes that only one '_across'. Add a for loop for more-than-one case. 
    if (merge_across_dims & (split_multiselected_dims | merge_across_dims_narm)) {
      if (!merge_across_dims_narm) {
        data_array_tmp <- array(bigmemory::as.matrix(data_array), dim = final_dims)
        tmp <- match(names(final_dims), names(dims_of_merge_dim))
        if (any(diff(tmp[!is.na(tmp)]) < 0)) { #need to reorder
          picked_common_vars[[across_inner_dim]] <- .aperm2(picked_common_vars[[across_inner_dim]], tmp[!is.na(tmp)])
        }
        metadata_tmp <- picked_common_vars[[across_inner_dim]]
      } else {
        tmp <- remove_additional_na_from_merge(
                 data_array = bigmemory::as.matrix(data_array),
                 merge_dim_metadata = picked_common_vars[[across_inner_dim]],
                 inner_dims_across_files, final_dims,
                 length_inner_across_dim)
        data_array_tmp <- tmp$data_array
        metadata_tmp <- tmp$merge_dim_metadata
      }

      if (length(data_array_tmp) != prod(final_dims_fake)) {
        stop(paste0("After reshaping, the data do not fit into the expected output dimension. ",
                    "Check if the reshaping parameters are used correctly."))
      }
      if (length(metadata_tmp) != prod(final_dims_fake_metadata)) {
        stop(paste0("After reshaping, the metadata do not fit into the expected output dimension. ",
                    "Check if the reshaping parameters are used correctly or contact support."))
      }

      #NOTE: When one file contains values for dicrete dimensions, rearrange the 
      #      chunks (i.e., work_piece) is necessary.
      if (split_multiselected_dims) {
        tmp <- rebuild_array_merge_split(
                 data_array = data_array_tmp, metadata = metadata_tmp, indices_chunk,
                 all_split_dims, final_dims_fake, across_inner_dim, length_inner_across_dim)
        data_array_tmp <- tmp$data_array
        metadata_tmp <- tmp$metadata
      }

      data_array <- array(data_array_tmp, dim = final_dims_fake)
      metadata_tmp <- array(metadata_tmp, dim = final_dims_fake_metadata)

      # If split_multiselected_dims + merge_across_dims, the dimension order may change above.
      # To get the user-required dim order, we need to reorder the array again.
      if (split_multiselected_dims) {
        if (inner_dim_pos_in_split_dims != 1) {
          correct_order <- match(names(final_dims_fake_output), names(final_dims_fake))
          data_array <- .aperm2(data_array, correct_order)
          correct_order_metadata <- match(names(final_dims_fake_output), names(all_split_dims[[across_inner_dim]]))
          metadata_tmp <- .aperm2(metadata_tmp, correct_order_metadata[!is.na(correct_order_metadata)])
        }
      }
      # Convert numeric back to dates
      if ('time' %in% synonims[[across_inner_dim]]) {
        metadata_tmp <- as.POSIXct(metadata_tmp, origin = "1970-01-01", tz = 'UTC')
      }

      picked_common_vars[[across_inner_dim]] <- metadata_tmp
      attr(picked_common_vars[[across_inner_dim]], 'variables') <- saved_reshaped_attr

    } else {  # ! (merge_across_dims + split_multiselected_dims) (old version)
      data_array <- array(bigmemory::as.matrix(data_array), dim = final_dims_fake)
      if (merge_across_dims) {
        # merge_across_dims = TRUE but (merge_across_dims_narm = F & split_multiselected_dims = F)

        inner_dim_pos <- which(names(dims_of_merge_dim) == inner_dims_across_files)
        file_dim_pos <- which(names(dims_of_merge_dim) == names(inner_dims_across_files))
        if (file_dim_pos < inner_dim_pos) {  #need to reorder
          tmp <- seq(1, length(dims_of_merge_dim))
          tmp[inner_dim_pos] <- file_dim_pos
          tmp[file_dim_pos] <- inner_dim_pos
          picked_common_vars[[across_inner_dim]] <- .aperm2(picked_common_vars[[across_inner_dim]], tmp)
        }
        metadata_tmp <- array(picked_common_vars[[across_inner_dim]], dim = final_dims_fake_metadata)
        # Convert numeric back to dates
        if ('time' %in% synonims[[across_inner_dim]]) {
          metadata_tmp <- as.POSIXct(metadata_tmp, origin = "1970-01-01", tz = 'UTC')
        }
        picked_common_vars[[across_inner_dim]] <- metadata_tmp
        attr(picked_common_vars[[across_inner_dim]], 'variables') <- saved_reshaped_attr
      }
      if (split_multiselected_dims & !is.null(picked_common_vars)) {
        if (!identical(inner_dim_has_split_dim, character(0))) {
          metadata_tmp <- array(picked_common_vars[[inner_dim_has_split_dim]], dim = final_dims_fake_metadata)
          # Convert numeric back to dates
          if (is(picked_common_vars[[inner_dim_has_split_dim]], 'POSIXct')) {
            metadata_tmp <- as.POSIXct(metadata_tmp, origin = "1970-01-01", tz = 'UTC')
          }
          picked_common_vars[[inner_dim_has_split_dim]] <- metadata_tmp
          attr(picked_common_vars[[inner_dim_has_split_dim]], 'variables') <- saved_reshaped_attr
        }
      }
    }
   
    gc()
    
    # Load metadata and remove the metadata folder
    if (!is.null(metadata_dims)) {
      loaded_metadata_files <- list.files(metadata_folder)

      if (!identical(loaded_metadata_files, character(0))) {  # old code