PlotMultiMap <- function(var, fun, ..., nrow = NULL, ncol = NULL, 
                         drawleg = 'S', titles = NULL, nbrks = 21, 
                         digits = 1, color_fun = clim.colors, sizeleg = 1,
                         subsampleg = 1) {
  if (is.array(var)) {
    var <- list(var)
  } else if (is.list(var)) {
    if (!all(sapply(var, is.numeric) || !all(sapply(var, is.array)))) {
      stop("Parameter 'var' must be an array or a list of arrays.")
    }
  } else {
    stop("Parameter 'var' must be an array or a list of arrays.")
  }

  if (!is.function(fun)) {
    stop("Parameter 'fun' must be a function.")
  }

  if (!is.null(nrow)) {
    if (!is.numeric(nrow)) {
      stop("Parameter 'nrow' must be numeric or NULL.")
    }
    nrow <- round(nrow)
  }

  if (!is.null(ncol)) {
    if (!is.numeric(ncol)) {
      stop("Parameter 'ncol' must be numeric or NULL.")
    }
    ncol <- round(ncol)
  }

  if (is.character(drawleg)) {
    if (drawleg %in% c('up', 'u', 'U', 'top', 't', 'T', 'north', 'n', 'N')) {
      drawleg <- 'N'
    } else if (drawleg %in% c('down', 'd', 'D', 'bottom', 'b', 'B', 'south', 's', 'S')) {
      drawleg <- 'S'
    } else if (drawleg %in% c('right', 'r', 'R', 'east', 'e', 'E')) {
      drawleg <- 'E'
    } else if (drawleg %in% c('left', 'l', 'L', 'west', 'w', 'W')) {
      drawleg <- 'W'
    }
  } else if (!is.logical(drawleg)) {
    stop("Parameter 'drawleg' must be either TRUE, FALSE or one of the positions in ?MultiMap.")
  }

  if (!all(sapply(titles, is.character))) {
    stop("Parameter 'titles' must be a vector of character strings.")
  }

  if (!is.numeric(nbrks)) {
    stop("Parameter 'nbrks' must be numeric.")
  }
  nbrks <- round(nbrks)

  if (!is.numeric(digits)) {
    stop("Parameter 'digits' must be numeric.")
  }
  digits <- round(digits)

  if (!is.function(color_fun)) {
    stop("Parameter 'color_fun' must be a function.")
  }

  if (!is.numeric(sizeleg)) {
    stop("Parameter 'sizeleg' must be numeric.")
  }
  sizeleg <- round(sizeleg)

  # Count the total number of maps and reorder each array of maps to have the lat and lon dimensions at the end.
  n_maps <- 0
  nlats <- NULL
  nlons <- NULL
  map_array_i <- 1
  for (map_array in var) {
    dimnames <- NULL
    if (!is.null(names(map_array))) {
      dimnames <- names(map_array)
    } else if (!is.null(attr(map_array, 'dimensions'))) {
      dimnames <- attr(map_array, 'dimensions')
    }
    if (!is.null(dimnames)) {
      if (!('lon' %in% dimnames || 'lat' %in% dimnames)) {
        stop("Parameter 'var' contains arrays without the dimensions 'lat' and/or 'lon'.")
      }
      lat_dim <- which(dimnames == 'lat')[1]
      lon_dim <- which(dimnames == 'lon')[1]
      map_array_nlats <- dim(map_array)[lat_dim]
      map_array_nlons <- dim(map_array)[lon_dim]
      n_maps <- n_maps + prod(dim(map_array)[-c(lat_dim, lon_dim)])
      var[[map_array_i]] <- aperm(var[[map_array_i]], c((1:length(dim(map_array)))[-c(lat_dim, lon_dim)], lat_dim, lon_dim))
    } else {
      .warning(paste0("Assuming the ", map_array_i, "th array provided in 'var' has 'lat' and 'lon' as last dimensions (right-most)."))
      dims <- tail(c(1, 1, dim(map_array)), 2)
      map_array_nlats <- dims[1]
      map_array_nlons <- dims[2]
      if (length(dim(var[[map_array_i]])) < 2) {
        dim(var[[map_array_i]]) <- dims
      }
      n_maps <- n_maps + prod(head(c(1, 1, dim(map_array)), length(dim(map_array))))
    }
    if (length(dim(var[[map_array_i]])) == 2) {
      dim(var[[map_array_i]]) <- c(1, dim(var[[map_array_i]]))
    }
    if (is.null(nlats)) {
      nlats <- map_array_nlats
      nlons <- map_array_nlons
    } else if (map_array_nlons != nlons || map_array_nlats != nlats) {
      stop(paste0("All the provided arrays in 'var' must have the same number of latitudes and longitudes (the ", map_array_i, "th array was the first unmatching array)."))
    }
    map_array_i <- map_array_i + 1
  }

  if (is.null(nrow) && is.null(ncol)) {
    ncol <- ceiling(sqrt(n_maps))
    nrow <- ceiling(n_maps/ncol)
  } else if (is.null(ncol)) {
    ncol <- ceiling(n_maps/nrow)
  } else if (is.null(nrow)) {
    nrow <- ceiling(n_maps/ncol)
  } else if (nrow * ncol < n_maps) {
    stop("There are more maps in 'var' than cells defined by 'nrow' x 'ncol'.")
  }

  if (is.logical(drawleg) && drawleg) {
    if (nrow > ncol) {
      drawleg <- 'S'
    } else {
      drawleg <- 'E'
    }
  }

  bc <- GenBrksCols(unlist(var), nbrks = nbrks, digits = digits, 
                    color_fun = color_fun)
  mat_layout <- rep(0, nrow*ncol)
  mat_layout[1:n_maps] <- 1:n_maps
  mat_layout <- matrix(mat_layout, nrow, ncol, byrow = TRUE)
  widths <- rep(1, ncol)
  heights <- rep(1, nrow)
  if (drawleg == 'N') {
    mat_layout <- rbind(rep(n_maps + 1, ncol), mat_layout)
    heights <- c(round(1 + sizeleg*nrow), heights * 10)
  } else if (drawleg == 'S') {
    mat_layout <- rbind(mat_layout, rep(n_maps + 1, ncol))
    heights <- c(heights * 10, round(1 + sizeleg*nrow))
  } else if (drawleg == 'W') {
    mat_layout <- cbind(rep(n_maps + 1, nrow), mat_layout)
    widths <- c(round(2 + sizeleg*ncol), widths * 10)
  } else if (drawleg == 'E') {
    mat_layout <- cbind(mat_layout, rep(n_maps + 1, nrow))
    widths <- c(widths * 10, round(2 + sizeleg*ncol))
  }
  layout(mat_layout, widths, heights)

  map_number <- 1
  lapply(var, function(x) {
    apply(x, (1:length(dim(x)))[1:(length(dim(x)) - 2)], function(y) {
      fun(y, toptitle = titles[map_number], brks = bc$brks, cols = bc$cols, drawleg = FALSE, ...)
      map_number <<- map_number + 1
    })
  })

  if (drawleg != FALSE) {
    ColorBar(bc$brks, bc$cols, vert = (drawleg %in% c('E', 'W')),
             subsampleg = subsampleg)
  }
}