diff --git a/R/CDORemap.R b/R/CDORemap.R
index f3909baefe83817c6e2a2779f854ccbe82a10638..8a210ca11948e994214eff7c053b19eec2330ce6 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -491,15 +491,18 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
     }
     names(attr(lons, 'variables')) <- lon_var_name
     names(attr(lats, 'variables')) <- lat_var_name
+    lons_lats_taken <- FALSE
     for (i in 1:total_slices) {
       tmp_file <- tempfile('R_CDORemap_', write_dir, fileext = '.nc')
       tmp_file2 <- tempfile('R_CDORemap_', write_dir, fileext = '.nc')
       if (!is.null(dims_to_iterate)) {
         slice_indices <- which(slices_to_iterate == i, arr.ind = TRUE)
         subset <- Subset(data_array, dims_to_iterate, as.list(slice_indices), drop = 'selected')
+#        dims_before_crop <- dim(subset)
         # Make sure subset goes along with metadata
         ArrayToNetCDF(setNames(list(subset, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file)
       } else {
+#        dims_before_crop <- dim(data_array)
         ArrayToNetCDF(setNames(list(data_array, lons, lats), c('var', lon_var_name, lat_var_name)), tmp_file)
       }
       sellonlatbox <- ''
@@ -508,7 +511,6 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
                                            ',', lat_extremes[1], ',', lat_extremes[2], ' -')
       }
       err <- try({
-## TODO: Here add sellonlatbox. Also check constantin's issue, may contain hint. Also search if possible to crop without
         system(paste0("cdo -s ", sellonlatbox, "remap", method, ",", grid, " ", tmp_file, " ", tmp_file2))
       })
       file.remove(tmp_file)
@@ -516,30 +518,45 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
         stop("CDO remap failed.")
       }
       ncdf_remapped <- nc_open(tmp_file2)
-      found_dim_names <- sapply(ncdf_remapped$var$var$dim, '[[', 'name')
-      found_lon_dim <- found_dim_names[which(found_dim_names %in% .KnownLonNames())[1]]
-      found_lat_dim <- found_dim_names[which(found_dim_names %in% .KnownLatNames())[1]]
-      found_lon_dim_size <- length(ncdf_remapped$dim[[found_lon_dim]]$vals)
-      found_lat_dim_size <- length(ncdf_remapped$dim[[found_lat_dim]]$vals)
-      found_lons <- ncvar_get(ncdf_remapped, 'lon', collapse_degen = FALSE)
-      found_lats <- ncvar_get(ncdf_remapped, 'lat', collapse_degen = FALSE)
-      if (length(dim(found_lons)) > 1) {
-        if (found_lon_dim < found_lat_dim) {
-          names(dim(found_lons)) <- c(found_lon_dim, found_lat_dim)
+      if (!lons_lats_taken) {
+        found_dim_names <- sapply(ncdf_remapped$var$var$dim, '[[', 'name')
+        found_lon_dim <- found_dim_names[which(found_dim_names %in% .KnownLonNames())[1]]
+        found_lat_dim <- found_dim_names[which(found_dim_names %in% .KnownLatNames())[1]]
+        found_lon_dim_size <- length(ncdf_remapped$dim[[found_lon_dim]]$vals)
+        found_lat_dim_size <- length(ncdf_remapped$dim[[found_lat_dim]]$vals)
+        found_lons <- ncvar_get(ncdf_remapped, 'lon', collapse_degen = FALSE)
+        found_lats <- ncvar_get(ncdf_remapped, 'lat', collapse_degen = FALSE)
+        if (length(dim(lons)) == length(dim(found_lons))) {
+          new_lon_name <- lon_dim
         } else {
-          names(dim(found_lons)) <- c(found_lat_dim, found_lon_dim)
+          new_lon_name <- found_lon_dim
         }
-      } else {
-        names(dim(found_lons)) <- found_lon_dim
-      }
-      if (length(dim(found_lats)) > 1) {
-        if (found_lon_dim < found_lat_dim) {
-          names(dim(found_lats)) <- c(found_lon_dim, found_lat_dim)
+        if (length(dim(lats)) == length(dim(found_lats))) {
+          new_lat_name <- lat_dim
         } else {
-          names(dim(found_lats)) <- c(found_lat_dim, found_lon_dim)
+          new_lat_name <- found_lat_dim
         }
-      } else {
-        names(dim(found_lats)) <- found_lat_dim
+        if (length(dim(found_lons)) > 1) {
+          if (which(sapply(ncdf_remapped$var$lon$dim, '[[', 'name') == found_lon_dim) < 
+              which(sapply(ncdf_remapped$var$lon$dim, '[[', 'name') == found_lat_dim)) {
+            names(dim(found_lons)) <- c(new_lon_name, new_lat_name)
+          } else {
+            names(dim(found_lons)) <- c(new_lat_name, new_lon_name)
+          }
+        } else {
+          names(dim(found_lons)) <- new_lon_name
+        }
+        if (length(dim(found_lats)) > 1) {
+          if (which(sapply(ncdf_remapped$var$lat$dim, '[[', 'name') == found_lon_dim) < 
+              which(sapply(ncdf_remapped$var$lat$dim, '[[', 'name') == found_lat_dim)) {
+            names(dim(found_lats)) <- c(new_lon_name, new_lat_name)
+          } else {
+            names(dim(found_lats)) <- c(new_lat_name, new_lon_name)
+          }
+        } else {
+          names(dim(found_lats)) <- new_lat_name
+        }
+        lons_lats_taken <- TRUE
       }
       if (!is.null(dims_to_iterate)) {
         if (is.null(result_array)) {
@@ -559,7 +576,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
         new_dims <- dim(data_array)
         new_dims[c(lon_dim, lat_dim)] <- c(found_lon_dim_size, found_lat_dim_size)
         result_array <- ncvar_get(ncdf_remapped, 'var', collapse_degen = FALSE)
-        names(dim(result_array)) <- names(new_dims)
+        dim(result_array) <- new_dims
       }
       nc_close(ncdf_remapped)
       file.remove(tmp_file2)
@@ -576,8 +593,6 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
     }
     attribute_backup[['dim']][which(names(dim(result_array)) == lon_dim)] <- dim(result_array)[lon_dim]
     attribute_backup[['dim']][which(names(dim(result_array)) == lat_dim)] <- dim(result_array)[lat_dim]
-    new_lon_name <- names(dim(found_lons))[which(names(dim(found_lons)) %in% .KnownLonNames())]
-    new_lat_name <- names(dim(found_lats))[which(names(dim(found_lats)) %in% .KnownLatNames())]
     names(attribute_backup[['dim']])[which(names(dim(result_array)) == lon_dim)] <- new_lon_name
     names(attribute_backup[['dim']])[which(names(dim(result_array)) == lat_dim)] <- new_lat_name
     if (!is.null(attribute_backup[['variables']]) && (length(attribute_backup[['variables']]) > 0)) {
diff --git a/R/Utils.R b/R/Utils.R
index 210cdfafe5b36af8391cadf7532a18f6499c1ce3..c5f034c35b33e13e5c318a6be59281f1071baaed 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -230,12 +230,12 @@
         grids_types <- unlist(lapply(grids_info, function (x) x[grep('gridtype', x) + 1]))
         grids_matches <- unlist(lapply(grids_info, function (x) {
           nlons <- if (length(grep('xsize', x)) > 0) {
-                     as.integer(x[grep('xsize', x) + 1])
+                     as.numeric(x[grep('xsize', x) + 1])
                    } else {
                      NA
                    }
           nlats <- if (length(grep('ysize', x)) > 0) {
-                    as.integer(x[grep('ysize', x) + 1])
+                    as.numeric(x[grep('ysize', x) + 1])
                   } else {
                     NA
                   }
@@ -279,14 +279,14 @@
         # Now we calculate the common grid type and its lons and lats
         if (length(grep('^t\\d{1,+}grid$', work_piece[['grid']])) > 0) {
           common_grid_type <- 'gaussian'
-          common_grid_res <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2])
+          common_grid_res <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2])
           nlonlat <- .t2nlatlon(common_grid_res)
           common_grid_lats <- nlonlat[1]
           common_grid_lons <- nlonlat[2]
         } else if (length(grep('^r\\d{1,+}x\\d{1,+}$', work_piece[['grid']])) > 0) {
           common_grid_type <- 'lonlat'
-          common_grid_lons <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2])
-          common_grid_lats <- as.integer(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][3])
+          common_grid_lons <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][2])
+          common_grid_lats <- as.numeric(strsplit(work_piece[['grid']], '[^0-9]{1,+}')[[1]][3])
         } else {
           stop("Error: Only supported grid types in parameter 'grid' are t<RES>grid and r<NX>x<NY>")
         }
@@ -612,15 +612,15 @@
                          ' 2>/dev/null'), intern = TRUE), split = ' ')
         years <- strsplit(system(paste('cdo showyear ', filein, 
                           ' 2>/dev/null'), intern = TRUE), split = ' ')
-        mons <- as.integer(mons[[1]][which(mons[[1]] != "")])
-        years <- as.integer(years[[1]][which(years[[1]] != "")])
+        mons <- as.numeric(mons[[1]][which(mons[[1]] != "")])
+        years <- as.numeric(years[[1]][which(years[[1]] != "")])
         time_indices <- ts(time_indices, start = c(years[1], mons[1]), 
                            end = c(years[length(years)], mons[length(mons)]),
                            frequency = 12)
         ltimes_list <- list()
         for (sdate in work_piece[['startdates']]) {
-          selected_time_indices <- window(time_indices, start = c(as.integer(
-                                   substr(sdate, 1, 4)), as.integer(substr(sdate, 5, 6))), 
+          selected_time_indices <- window(time_indices, start = c(as.numeric(
+                                   substr(sdate, 1, 4)), as.numeric(substr(sdate, 5, 6))), 
                                    end = c(3000, 12), frequency = 12, extend = TRUE)
           selected_time_indices <- selected_time_indices[work_piece[['leadtimes']]]
           ltimes_list <- c(ltimes_list, list(selected_time_indices))
@@ -1158,7 +1158,7 @@
 
   if (give_warning) {
     .warning(paste0("Too complex path pattern specified for ", dataset_name, 
-                    ". Double check carefully the 'source_files' fetched for this dataset or specify a simpler path pattern."))
+                    ". Double check carefully the '$Files' fetched for this dataset or specify a simpler path pattern."))
   }
 
   if (permissive) {
@@ -1168,6 +1168,60 @@
   }
 }
 
+.FindTagValue <- function(path_with_globs_and_tag, actual_path, tag) {
+  tag <- paste0('\\$', tag, '\\$')
+  path_with_globs_and_tag <- paste0('^', path_with_globs_and_tag, '$')
+  parts <- strsplit(path_with_globs_and_tag, '*', fixed = TRUE)[[1]]
+  parts <- as.list(parts[grep(tag, parts)])
+  longest_couples <- c()
+  pos_longest_couples <- c()
+  found_value <- NULL
+  for (i in 1:length(parts)) {
+    parts[[i]] <- strsplit(parts[[i]], tag)[[1]]
+    if (length(parts[[i]]) == 1) {
+      parts[[i]] <- c(parts[[i]], '')
+    }
+    len_parts <- sapply(parts[[i]], nchar)
+    len_couples <- len_parts[-length(len_parts)] + len_parts[2:length(len_parts)]
+    pos_longest_couples <- c(pos_longest_couples, which.max(len_couples))
+    longest_couples <- c(longest_couples, max(len_couples))
+  }
+  chosen_part <- which.max(longest_couples)
+  parts[[chosen_part]] <- parts[[chosen_part]][pos_longest_couples[chosen_part]:(pos_longest_couples[chosen_part] + 1)]
+  if (nchar(parts[[chosen_part]][1]) >= nchar(parts[[chosen_part]][2])) {
+    if (nchar(parts[[chosen_part]][1]) > 0) {
+      matches <- gregexpr(parts[[chosen_part]][1], actual_path)[[1]]
+      if (length(matches) == 1) {
+        match_left <- matches
+        actual_path <- substr(actual_path, match_left + attr(match_left, 'match.length'), nchar(actual_path))
+      }
+    }
+    if (nchar(parts[[chosen_part]][2]) > 0) {
+      matches <- gregexpr(parts[[chosen_part]][2], actual_path)[[1]]
+      if (length(matches) == 1) {
+        match_right <- matches
+        found_value <- substr(actual_path, 0, match_right - 1)
+      }
+    }
+  } else {
+    if (nchar(parts[[chosen_part]][2]) > 0) {
+      matches <- gregexpr(parts[[chosen_part]][2], actual_path)[[1]]
+      if (length(matches) == 1) {
+        match_right <- matches
+        actual_path <- substr(actual_path, 0, match_right - 1)
+      }
+    }
+    if (nchar(parts[[chosen_part]][1]) > 0) {
+      matches <- gregexpr(parts[[chosen_part]][1], actual_path)[[1]]
+      if (length(matches) == 1) {
+        match_left <- matches
+        found_value <- substr(actual_path, match_left + attr(match_left, 'match.length'), nchar(actual_path))
+      }
+    }
+  }
+  found_value
+}
+
 .FilterUserGraphicArgs <- function(excludedArgs, ...) {
   # This function filter the extra graphical parameters passed by the user in
   # a plot function, excluding the ones that the plot function uses by default.
@@ -1440,3 +1494,90 @@
   plot(0, type = 'n', axes = FALSE, ann = FALSE)
   par(mfg = next_attempt)
 }
+
+# Function to permute arrays of non-atomic elements (e.g. POSIXct)
+.aperm2 <- function(x, new_order) {
+  y <- array(1:length(x), dim = dim(x))
+  y <- aperm(y, new_order)
+  old_dims <- dim(x)
+  x <- x[as.vector(y)]
+  dim(x) <- old_dims[new_order]
+  x
+}
+
+# This function is a helper for the function .MergeArrays.
+# It expects as inputs two named numeric vectors, and it extends them
+# with dimensions of length 1 until an ordered common dimension
+# format is reached.
+.MergeArrayDims <- function(dims1, dims2) {
+  new_dims1 <- c()
+  new_dims2 <- c()
+  while (length(dims1) > 0) {
+    if (names(dims1)[1] %in% names(dims2)) {
+      pos <- which(names(dims2) == names(dims1)[1])
+      dims_to_add <- rep(1, pos - 1)
+      if (length(dims_to_add) > 0) {
+        names(dims_to_add) <- names(dims2[1:(pos - 1)])
+      }
+      new_dims1 <- c(new_dims1, dims_to_add, dims1[1])
+      new_dims2 <- c(new_dims2, dims2[1:pos])
+      dims1 <- dims1[-1]
+      dims2 <- dims2[-c(1:pos)]
+    } else {
+      new_dims1 <- c(new_dims1, dims1[1])
+      new_dims2 <- c(new_dims2, 1)
+      names(new_dims2)[length(new_dims2)] <- names(dims1)[1]
+      dims1 <- dims1[-1]
+    }
+  }
+  if (length(dims2) > 0) {
+    dims_to_add <- rep(1, length(dims2))
+    names(dims_to_add) <- names(dims2)
+    new_dims1 <- c(new_dims1, dims_to_add)
+    new_dims2 <- c(new_dims2, dims2)
+  }
+  list(new_dims1, new_dims2)
+}
+
+# This function takes two named arrays and merges them, filling with
+# NA where needed.
+# dim(array1)
+#          'b'   'c'         'e'   'f'
+#           1     3           7     9
+# dim(array2)
+#    'a'   'b'         'd'         'f'   'g'
+#     2     3           5           9     11
+# dim(.MergeArrays(array1, array2, 'b'))
+#    'a'   'b'   'c'   'e'   'd'   'f'   'g'
+#     2     4     3     7     5     9     11
+.MergeArrays <- function(array1, array2, along) {
+  if (!(identical(names(dim(array1)), names(dim(array2))) &&
+      identical(dim(array1)[-which(names(dim(array1)) == along)],
+                dim(array2)[-which(names(dim(array2)) == along)]))) {
+    new_dims <- .MergeArrayDims(dim(array1), dim(array2))
+    dim(array1) <- new_dims[[1]]
+    dim(array2) <- new_dims[[2]]
+    for (j in 1:length(dim(array1))) {
+      if (names(dim(array1))[j] != along) {
+        if (dim(array1)[j] != dim(array2)[j]) {
+          if (which.max(c(dim(array1)[j], dim(array2)[j])) == 1) {
+            na_array_dims <- dim(array2)
+            na_array_dims[j] <- dim(array1)[j] - dim(array2)[j]
+            na_array <- array(dim = na_array_dims)
+            array2 <- abind(array2, na_array, along = j)
+            names(dim(array2)) <- names(na_array_dims)
+          } else {
+            na_array_dims <- dim(array1)
+            na_array_dims[j] <- dim(array2)[j] - dim(array1)[j]
+            na_array <- array(dim = na_array_dims)
+            array1 <- abind(array1, na_array, along = j)
+            names(dim(array1)) <- names(na_array_dims)
+          }
+        }
+      }
+    }
+  }
+  array1 <- abind(array1, array2, along = which(names(dim(array1)) == along))
+  names(dim(array1)) <- names(dim(array2))
+  array1
+}