diff --git a/R/CDORemap.R b/R/CDORemap.R
index 11aa2a0670dcbacf177329ea4931502ae98a215a..efd9bf216c27e02795f6864938360c5f23e7d4f4 100644
--- a/R/CDORemap.R
+++ b/R/CDORemap.R
@@ -60,7 +60,7 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
       }
     }
     names(array_dims) <- c(new_lat_dim_name, new_lon_dim_name)
-    data_array <- array(NA, array_dims)
+    data_array <- array(as.numeric(NA), array_dims)
   }
   if (!(is.logical(data_array) || is.numeric(data_array)) || !is.array(data_array)) {
     stop("Parameter 'data_array' must be a numeric array.")
@@ -500,6 +500,12 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
     }
     names(attr(lons, 'variables')) <- lon_var_name
     names(attr(lats, 'variables')) <- lat_var_name
+    if (!is.null(attr(lons, 'variables')[[1]][['dim']])) {
+      attr(lons, 'variables')[[1]][['dim']] <- NULL
+    }
+    if (!is.null(attr(lats, 'variables')[[1]][['dim']])) {
+      attr(lats, 'variables')[[1]][['dim']] <- NULL
+    }
     lons_lats_taken <- FALSE
     for (i in 1:total_slices) {
       tmp_file <- tempfile('R_CDORemap_', write_dir, fileext = '.nc')
@@ -535,8 +541,33 @@ CDORemap <- function(data_array = NULL, lons, lats, grid, method,
         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)
+        found_var_names <- names(ncdf_remapped$var)
+        found_lon_var_name <- which(found_var_names %in% .KnownLonNames())
+        found_lat_var_name <- which(found_var_names %in% .KnownLatNames())
+        if (length(found_lon_var_name) > 0) {
+          found_lon_var_name <- found_var_names[found_lon_var_name[1]]
+        } else {
+          found_lon_var_name <- NULL
+        }
+        if (length(found_lat_var_name) > 0) {
+          found_lat_var_name <- found_var_names[found_lat_var_name[1]]
+        } else {
+          found_lat_var_name <- NULL
+        }
+        if (length(found_lon_var_name) > 0) {
+          found_lons <- ncvar_get(ncdf_remapped, found_lon_var_name, 
+                                  collapse_degen = FALSE)
+        } else {
+          found_lons <- ncdf_remapped$dim[[found_lon_dim]]$vals
+          dim(found_lons) <- found_lon_dim_size
+        }
+        if (length(found_lat_var_name) > 0) {
+          found_lats <- ncvar_get(ncdf_remapped, found_lat_var_name, 
+                                  collapse_degen = FALSE)
+        } else {
+          found_lats <- ncdf_remapped$dim[[found_lat_dim]]$vals
+          dim(found_lats) <- found_lat_dim_size
+        }
         if (length(dim(lons)) == length(dim(found_lons))) {
           new_lon_name <- lon_dim
         } else {
diff --git a/R/Clim.R b/R/Clim.R
index 615f845a898c00a72e1a7c6ab2c7191327818b73..4fdfa038f9ef4b5396870173f9d35346a914a9eb 100644
--- a/R/Clim.R
+++ b/R/Clim.R
@@ -47,8 +47,8 @@ Clim <- function(var_exp, var_obs, memb = TRUE, kharin = FALSE, NDV = FALSE) {
     trend_obs <- array(dim = dim(var_exp))
     trend_exp <- array(dim = dim(var_exp))
     for (jdate in 1:dimexp[3]) {
-      trend_exp[, , jdate, , , , ] <- tmp_exp$trend[, , 4, , , , ] 
-                                      + jdate * tmp_exp$trend[, , 2, , , , ]
+      trend_exp[, , jdate, , , , ] <- tmp_exp$trend[, , 4, , , , ] +
+                                       jdate * tmp_exp$trend[, , 2, , , , ]
       tmp_obs2 <- MeanListDim(tmp_obs$trend,c(2, 1))
       trend_obs[, , jdate, , , , ] <- InsertDim(InsertDim(tmp_obs2[4, , , , ] + 
                                       jdate * tmp_obs2[2, , , , ], 1, dimexp[1]),
@@ -66,13 +66,13 @@ Clim <- function(var_exp, var_obs, memb = TRUE, kharin = FALSE, NDV = FALSE) {
     reg_obs <- array(dim = dim(var_exp))
     reg_exp <- array(dim = dim(var_exp))
     for (jdate in 1:dimexp[3]) {
-      reg_exp[, , jdate, , , , ] <- tmp_exp$regression[, , 4, , , , ] 
-                                    + iniexp[, , jdate, , , , ] * 
+      reg_exp[, , jdate, , , , ] <- tmp_exp$regression[, , 4, , , , ] +
+                                     iniexp[, , jdate, , , , ] * 
                                     tmp_exp$regression[, , 2, , , , ]
       tmp_obs2 <- MeanListDim(tmp_obs$regression,c(2, 1))
-      reg_obs[, , jdate, , , , ] <- InsertDim(InsertDim(tmp_obs2[4, , , , ] 
-                                    + MeanListDim(iniobs,c(2, 1))[jdate, , , , ] 
-                                    * tmp_obs2[2, , , , ], 1, dimexp[1]), 
+      reg_obs[, , jdate, , , , ] <- InsertDim(InsertDim(tmp_obs2[4, , , , ] +
+                                     MeanListDim(iniobs,c(2, 1))[jdate, , , , ] *
+                                     tmp_obs2[2, , , , ], 1, dimexp[1]), 
                                     2, dimexp[2])
     }
     out_clim_exp <- reg_exp - reg_obs + InsertDim(InsertDim(InsertDim(
diff --git a/R/ColorBar.R b/R/ColorBar.R
index 4f92a2a159e9ba504a14a955270f240e338a95be..d47c3103a25f16c8054d649e72c74b0c6109fc24 100644
--- a/R/ColorBar.R
+++ b/R/ColorBar.R
@@ -223,6 +223,10 @@ ColorBar <- function(brks = NULL, cols = NULL, vertical = TRUE,
     stop("Parameter 'subsampleg' must be numeric.")
   }
   subsampleg <- round(subsampleg)
+  draw_labels <- TRUE
+  if ((subsampleg) < 1) {
+    draw_labels <- FALSE
+  }
 
   # Check plot
   if (!is.logical(plot)) {
@@ -433,7 +437,11 @@ ColorBar <- function(brks = NULL, cols = NULL, vertical = TRUE,
     labels <- c(labels, extra_labels)
     tick_reorder <- sort(at, index.return = TRUE)
     at <- tick_reorder$x
-    labels <- labels[tick_reorder$ix]
+    if (draw_labels) {
+      labels <- labels[tick_reorder$ix]
+    } else {
+      labels <- FALSE
+    }
     axis(d, at = at, tick = draw_ticks, labels = labels, cex.axis = cex_labels, tcl = cex_ticks)
     par(saved_pars)
   }
diff --git a/R/Load.R b/R/Load.R
index d11e809d22449803f00bcf4a121579b2347aff1a..5cebec71106778ac3e1205aa6f432f332e552853 100644
--- a/R/Load.R
+++ b/R/Load.R
@@ -18,9 +18,12 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   }
   load_parameters <- lapply(parameter_names, get, envir = environment())
   names(load_parameters) <- parameter_names
-  parameters_to_show <- c('var', 'exp', 'obs', 'sdates', 'grid', 'output', 'storefreq')
+  parameters_to_show <- c('var', 'exp', 'obs', 'sdates', 'nmember', 'leadtimemin',
+                          'leadtimemax', 'latmin', 'latmax', 'lonmin', 'lonmax',
+                          'output', 'grid', 'storefreq')
   load_parameters <- c(load_parameters[parameters_to_show], load_parameters[-match(parameters_to_show, names(load_parameters))])
-  message(paste("* The load call you issued is:\n*   Load(", 
+  if (!silent) {
+    message(paste("* The load call you issued is:\n*   Load(", 
             paste(strwrap(
               paste(unlist(lapply(names(load_parameters[1:length(parameters_to_show)]), 
               function(x) paste(x, '=', 
@@ -46,23 +49,7 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
               collapse = ', '), width = getOption('width') - 9, indent = 0, exdent = 8), collapse = '\n*'),
             ", ...)\n* See the full call in '$load_parameters' after Load() finishes.", sep = ''))
 
-  # .message("* The load call you issued is:")
-  # .message("*   Load(")
-  # .message( 
-  #   strwrap(
-  #     paste(
-  #       unlist(lapply(names(load_parameters[1:length(parameters_to_show)]), 
-  #             function(x) paste(x, '=', 
-  #               if (x == 'sdates' && length(load_parameters[[x]]) > 4) {
-  #                 paste0("c('", load_parameters[[x]][1], "', '", load_parameters[[x]][2], 
-  #                        "', ..., '", tail(load_parameters[[x]], 1), "')")
-  #               } else {
-  #                 paste(deparse(load_parameters[[x]]), collapse = '')
-  #               }))), 
-  #       collapse = ', '), 
-  #     width = getOption('width') - 9, indent = 0, exdent = 8))
-  # .message("*        , ...)")
-  # .message("* See the full call in '$load_parameters' after Load() finishes.")
+  }
 
   # Run Load() error-aware, so that it always returns something
   errors <- try({
@@ -392,8 +379,10 @@ Load <- function(var, exp = NULL, obs = NULL, sdates, nmember = NULL,
   replace_globs <- path_glob_permissive %in% c('no', 'partial')
 
   # If not all data has been provided in 'exp' and 'obs', configuration file is read.
-  if (length(exps_to_fetch) > 0 || length(obs_to_fetch) > 0) {
-    .message("Some 'path's not explicitly provided in 'exp' and 'obs', so will now proceed to open the configuration file.")
+  if ((length(exps_to_fetch) > 0 || length(obs_to_fetch) > 0)) {
+    if (!silent) {
+      .message("Some 'path's not explicitly provided in 'exp' and 'obs', so will now proceed to open the configuration file.")
+    }
     data_info <- ConfigFileOpen(configfile, silent, TRUE)
 
     # Check that the var, exp and obs parameters are right and keep the entries
diff --git a/R/Subset.R b/R/Subset.R
index a084a1c2ff3526c460bdcc87e8741cfa355d0386..d10629c59dea260cc4d5b9e37dcb14706aae6cc3 100644
--- a/R/Subset.R
+++ b/R/Subset.R
@@ -80,6 +80,10 @@ Subset <- function(x, along, indices, drop = FALSE) {
       names(metadata[['dim']]) <- dim_names[-dims_to_drop]
       metadata[['dimensions']] <- dim_names[-dims_to_drop]
     }
+    if (length(metadata[['dim']]) == 0) {
+      metadata['dim'] <- list(NULL)
+      metadata['dimensions'] <- list(NULL)
+    }
   } else if (is.character(dim_names)) {
     names(metadata[['dim']]) <- dim_names
     metadata[['dimensions']] <- dim_names
diff --git a/R/Utils.R b/R/Utils.R
index 603d47163f48771ae6310b7e27e13666e98939dd..31f9ae071fbb78f4d398cf447cdf907f9aaaa645 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -1565,33 +1565,41 @@
 #    '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)
+  if (!(is.null(array1) || is.null(array2))) {
+    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)
+            }
           }
         }
       }
     }
+    if (!(along %in% names(dim(array2)))) {
+      stop("The dimension specified in 'along' is not present in the ",
+           "provided arrays.")
+    }
+    array1 <- abind(array1, array2, along = which(names(dim(array1)) == along))
+    names(dim(array1)) <- names(dim(array2))
+  } else if (is.null(array1)) {
+    array1 <- array2
   }
-  array1 <- abind(array1, array2, along = which(names(dim(array1)) == along))
-  names(dim(array1)) <- names(dim(array2))
   array1
 }
diff --git a/Subset.R b/Subset.R
new file mode 100644
index 0000000000000000000000000000000000000000..c82b0c67f6804ab2a3ce9127ede530a6fbe4cbb2
--- /dev/null
+++ b/Subset.R
@@ -0,0 +1,93 @@
+Subset <- function(x, along, indices, drop = FALSE) {
+  # Check x
+  if (!is.array(x)) {
+    stop("Input array 'x' must be a numeric array.")
+  }
+
+  # Take the input array dimension names
+  dim_names <- attr(x, 'dimensions')
+  if (!is.character(dim_names)) {
+    dim_names <- names(dim(x))
+  }
+  if (!is.character(dim_names)) {
+    if (any(sapply(along, is.character))) {
+      stop("The input array 'x' doesn't have labels for the dimensions but the parameter 'along' contains dimension names.")
+    }
+  }
+
+  # Check along
+  if (any(sapply(along, function(x) !is.numeric(x) && !is.character(x)))) {
+    stop("All provided dimension indices in 'along' must be integers or character strings.")
+  }
+  if (any(sapply(along, is.character))) {
+    req_dimnames <- along[which(sapply(along, is.character))]
+    if (length(unique(req_dimnames)) < length(req_dimnames)) {
+      stop("The parameter 'along' must not contain repeated dimension names.")
+    }
+    along[which(sapply(along, is.character))] <- match(req_dimnames, dim_names)
+    if (any(is.na(along))) {
+      stop("Could not match all dimension names in 'indices' with dimension names in input array 'x'.")
+    }
+    along <- as.numeric(along)
+  }
+
+  # Check indices
+  if (!is.list(indices)) {
+    indices <- list(indices)
+  }
+
+  # Check parameter drop
+  dims_to_drop <- c()
+  if (is.character(drop)) {
+    if (drop == 'all') {
+      drop <- TRUE
+    } else if (any(drop %in% c('selected', 'non-selected', 'none'))) {
+      if (drop == 'selected') {
+        dims_to_drop <- along[which(sapply(indices, length) == 1)]
+      } else if (drop == 'non-selected') {
+        dims_to_drop <- dim(x) == 1
+        dims_to_drop[along] <- FALSE
+        dims_to_drop <- which(dims_to_drop)
+      }
+      drop <- FALSE
+    } else {
+      stop("Parameter 'drop' must be one of TRUE, FALSE, 'all', 'selected', 'non-selected', 'none'.")
+    }
+  } else if (!is.logical(drop)) {
+    stop("Parameter 'drop' must be one of TRUE, FALSE, 'all', 'selected', 'non-selected', 'none'.")
+  }
+
+  # Take the subset
+  nd <- length(dim(x))
+  index <- as.list(rep(TRUE, nd))
+  index[along] <- indices
+  subset <- eval(as.call(c(as.name("["), as.name("x"), index, drop = drop)))
+  # If dropped all dimensions, need to drop dimnames too
+  if (is.character(dim_names) && drop == TRUE) {
+    dim_names_to_remove <- unique(c(along[which(sapply(indices, length) == 1)],
+                                    which(dim(x) == 1)))
+    if (length(dim_names_to_remove) > 0) {
+      dim_names <- dim_names[-dim_names_to_remove]
+    }
+  }
+
+  # Amend the final dimensions and put dimnames and attributes
+  metadata <- attributes(x)
+  metadata[['dim']] <- dim(subset)
+  if (length(dims_to_drop) > 0) {
+    metadata[['dim']] <- metadata[['dim']][-dims_to_drop]
+    if (is.character(dim_names)) {
+      names(metadata[['dim']]) <- dim_names[-dims_to_drop]
+      if ('dimensions' %in% names(attributes(x))) {
+        metadata[['dimensions']] <- dim_names[-dims_to_drop]
+      }
+    }
+  } else if (is.character(dim_names)) {
+    names(metadata[['dim']]) <- dim_names
+    if ('dimensions' %in% names(attributes(x))) {
+      metadata[['dimensions']] <- dim_names
+    }
+  }
+  attributes(subset) <- metadata
+  subset
+}
diff --git a/man/ColorBar.Rd b/man/ColorBar.Rd
index fb7cc7e508b23f6d87773077c3313436a1f03827..51ca567029edd49ca38e0f0daac7816acee1ba29 100644
--- a/man/ColorBar.Rd
+++ b/man/ColorBar.Rd
@@ -33,7 +33,7 @@ TRUE/FALSE for vertical/horizontal colour bar (disregarded if plot = FALSE).
   }
   \item{subsampleg}{
 The first of each subsampleg breaks will be ticked on the colorbar.\cr
-Takes by default an approximation of a value that yields a readable tick arrangement (extreme breaks always ticked). See the code of the function for details or use 'extra_labels' for customized tick arrangements.
+Takes by default an approximation of a value that yields a readable tick arrangement (extreme breaks always ticked). If set to 0 or lower, no labels are drawn. See the code of the function for details or use 'extra_labels' for customized tick arrangements.
   }
   \item{bar_limits}{
 Vector of two numeric values with the extremes of the range of values represented in the colour bar. If 'var_limits' go beyond this interval, the drawing of triangle extremes is triggered at the corresponding sides, painted in 'col_inf' and 'col_sup'. Either of them can be set as NA and will then take as value the corresponding extreme in 'var_limits' (hence a triangle end won't be triggered for these sides). Takes as default the extremes of 'brks' if available, else the same values as 'var_limits'.
diff --git a/man/PlotLayout.Rd b/man/PlotLayout.Rd
index 9e0c2cb7937db03122e02a89a83365f2cc0bb54d..34331a3776bb2eddc5179e7c06abaa3d7f85afe0 100644
--- a/man/PlotLayout.Rd
+++ b/man/PlotLayout.Rd
@@ -41,7 +41,7 @@ Multi-dimensional array with at least the dimensions expected by the specified p
 Parameters to be sent to the plotting function 'fun'. If multiple arrays are provided in 'var' and multiple functions are provided in 'fun', the parameters provided through \dots will be sent to all the plot functions, as common parameters. To specify concrete arguments for each of the plot functions see parameter 'special_args'.
   }
   \item{special_args}{
-List of sub-lists, each sub-list having extra arguments for each of the plot functions provided in 'fun'.
+List of sub-lists, each sub-list having specific extra arguments for each of the plot functions provided in 'fun'. If you want to fix a different value for each plot in the layout you can do so by a) splitting your array into a list of sub-arrays (each with the data for one plot) and providing it as parameter 'var', b) providing a list of named sub-lists in 'special_args', where the names of each sub-list match the names of the parameters to be adjusted, and each value in a sub-list contains the value of the corresponding parameter. 
   }
   \item{nrow}{
 Numeric value to force the number of rows in the automatically generated layout. If higher than the required, this will yield blank cells in the layout (which can then be populated). If lower than the required the function will stop. By default it is configured to arrange the layout in a shape as square as possible. Blank cells can be manually populated after with customized plots (see SwitchTofigure).