diff --git a/R/PlotLayout.R b/R/PlotLayout.R
index 00ba6747e98d4213ecb0be30e4b6b6e69235fa2c..c442bf77311df2a1f6b5aa5149b2ffbb3f66ca25 100644
--- a/R/PlotLayout.R
+++ b/R/PlotLayout.R
@@ -548,26 +548,45 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
   title_margin <- 0.5 * title_cex * title_margin_scale
   subtitle_cex <- 1.5 * subtitle_scale
   subtitle_margin <- 0.5 * sqrt(nrow * ncol) * subtitle_cex * subtitle_margin_scale
-  mat_layout <- 1:(nrow * ncol) + ifelse(drawleg != FALSE, 1, 0)
+  mat_layout <- 1:(nrow * ncol)
+  if (drawleg != FALSE) {
+    if (all(fun %in% 'PlotMostLikelyQuantileMap')) {  #multi_colorbar
+      multi_colorbar <- TRUE
+      cat_dim <- list(...)$cat_dim
+      if (is.null(cat_dim)) cat_dim <- 'bin'  # default
+      nmap <- as.numeric(dim(var[[1]])[cat_dim])
+      minimum_value <- ceiling(1 / nmap * 10 * 1.1) * 10
+      display_range = c(minimum_value, 100)
+      mat_layout <- mat_layout + nmap
+    } else {
+      multi_colorbar <- FALSE
+      mat_layout <- mat_layout + 1
+    }
+  }
   mat_layout <- matrix(mat_layout, nrow, ncol, byrow = layout_by_rows)
   fsu <- figure_size_units <- 10  # unitless
   widths <- rep(fsu, ncol)
   heights <- rep(fsu, nrow)
-  n_figures <- nrow * ncol
-  if (length(row_titles) > 0) {
-    mat_layout <- cbind(rep(0, dim(mat_layout)[1]), mat_layout)
-    widths <- c(((subtitle_cex + subtitle_margin / 2) * cs / device_size[1]) * ncol * fsu, widths)
-  }
-  if (length(col_titles) > 0) {
-    mat_layout <- rbind(rep(0, dim(mat_layout)[2]), mat_layout)
-    heights <- c(((subtitle_cex + subtitle_margin) * cs / device_size[2]) * nrow * fsu, heights)
-  }
+   # Useless
+#  n_figures <- nrow * ncol
+
   if (drawleg != FALSE) {
     if (drawleg == 'N') {
       mat_layout <- rbind(rep(1, dim(mat_layout)[2]), mat_layout)
       heights <- c(round(bar_scale * 2 * nrow), heights)
     } else if (drawleg == 'S') {
-      mat_layout <- rbind(mat_layout, rep(1, dim(mat_layout)[2]))
+      if (multi_colorbar) {
+        new_mat_layout <- c()
+        for (i_col in 1:ncol) {
+          new_mat_layout <- c(new_mat_layout, rep(mat_layout[, i_col], nmap))
+        }
+        new_mat_layout <- matrix(new_mat_layout, nrow, nmap * ncol)
+        colorbar_row <- rep(1:nmap, each = ncol)
+        mat_layout <- rbind(new_mat_layout, as.numeric(colorbar_row))
+        widths <- rep(widths, nmap)
+      } else {
+        mat_layout <- rbind(mat_layout, rep(1, dim(mat_layout)[2]))
+      }
       heights <- c(heights, round(bar_scale * 2 * nrow))
     } else if (drawleg == 'W') {
       mat_layout <- cbind(rep(1, dim(mat_layout)[1]), mat_layout)
@@ -576,8 +595,20 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
       mat_layout <- cbind(mat_layout, rep(1, dim(mat_layout)[1]))
       widths <- c(widths, round(bar_scale * 3 * ncol))
     }
-    n_figures <- n_figures + 1
+   # Useless
+#    n_figures <- n_figures + 1
   }
+
+  # row and col titles
+  if (length(row_titles) > 0) {
+    mat_layout <- cbind(rep(0, dim(mat_layout)[1]), mat_layout)
+    widths <- c(((subtitle_cex + subtitle_margin / 2) * cs / device_size[1]) * ncol * fsu, widths)
+  }
+  if (length(col_titles) > 0) {
+    mat_layout <- rbind(rep(0, dim(mat_layout)[2]), mat_layout)
+    heights <- c(((subtitle_cex + subtitle_margin) * cs / device_size[2]) * nrow * fsu, heights)
+  }
+  # toptitle
   if (toptitle != '') {
     mat_layout <- rbind(rep(0, dim(mat_layout)[2]), mat_layout)
     heights <-  c(((title_cex + title_margin) * cs / device_size[2]) * nrow * fsu, heights)
@@ -590,13 +621,32 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
       bar_extra_margin[2] <- bar_extra_margin[2] + (subtitle_cex + subtitle_margin / 2) * 
                                                    bar_left_shift_scale
     }
-    ColorBar(colorbar$brks, colorbar$cols, vertical, subsampleg, 
-             bar_limits, var_limits,
-             triangle_ends = triangle_ends, col_inf = colorbar$col_inf,
-             col_sup = colorbar$col_sup, color_fun, plot = TRUE, draw_bar_ticks, 
-             draw_separators, triangle_ends_scale, bar_extra_labels,
-             units, units_scale, bar_label_scale, bar_tick_scale, 
-             bar_extra_margin, bar_label_digits)
+
+    if (multi_colorbar) {  # multiple colorbar
+      if (!is.null(list(...)$bar_titles)) {
+        bar_titles <- list(...)$bar_titles
+      } else {
+        bar_titles <- NULL
+      }
+      GradientCatsColorBar(nmap = nmap, 
+                           brks = brks, cols = cols, vertical = vertical, subsampleg = subsampleg,
+                           bar_limits = display_range, var_limits = var_limits,
+                           triangle_ends = triangle_ends, plot = TRUE,
+                           draw_separators = draw_separators,
+                           bar_titles = bar_titles, title_scale = units_scale,
+                           label_scale = bar_label_scale, extra_margin = bar_extra_margin)
+
+    } else {  # one colorbar
+      ColorBar(brks = colorbar$brks, cols = colorbar$cols, vertical = vertical, subsampleg = subsampleg,
+               bar_limits = bar_limits, var_limits = var_limits,
+               triangle_ends = triangle_ends, col_inf = colorbar$col_inf,
+               col_sup = colorbar$col_sup, color_fun = color_fun, plot = TRUE, draw_ticks = draw_bar_ticks,
+               draw_separators = draw_separators, triangle_ends_scale = triangle_ends_scale,
+               extra_labels = bar_extra_labels,
+               title = units, title_scale = units_scale, label_scale = bar_label_scale, tick_scale = bar_tick_scale,
+               extra_margin = bar_extra_margin, label_digits = bar_label_digits)
+
+    }
   }
 
   # Draw titles
@@ -665,24 +715,26 @@ PlotLayout <- function(fun, plot_dims, var, ..., special_args = NULL,
       # For each of the arrays provided in that array
       apply(x, (1:length(dim(x)))[-plot_dim_indices], 
             function(y) {
-        # Do the plot
-        fun_args <- c(list(y, toptitle = titles[plot_number]), list(...), 
+        # Do the plot. colorbar is not drew.
+        fun_args <- c(list(y, toptitle = titles[plot_number], drawleg = FALSE), list(...),
                       special_args[[array_number]])
-        funct <- fun[[array_number]]
+#        funct <- fun[[array_number]]
         if (fun[[array_number]] %in% c('PlotEquiMap', 'PlotStereoMap')) {
           fun_args <- c(fun_args, list(brks = colorbar$brks, cols = colorbar$cols, 
                                        col_inf = colorbar$col_inf, 
                                        col_sup = colorbar$col_sup, 
-                                       drawleg = FALSE,
-                                       title_scale = subplot_titles_scale))
+                                       title_scale = subplot_titles_scale  # when all the functions have this argument, put it above in fun_args
+                                       ))
         } else if (fun[[array_number]] == c('PlotSection')) {
-          fun_args <- c(fun_args, list(brks = colorbar$brks, cols = colorbar$cols, 
-                                       drawleg = FALSE))
+          fun_args <- c(fun_args, list(brks = colorbar$brks, cols = colorbar$cols))
 
+        } else if (fun[[array_number]] %in% 'PlotMostLikelyQuantileMap') {
+          #TODO: pre-generate colorbar params? like above
+          fun_args <- c(fun_args, list(brks = brks, cols = cols))
         }
         do.call(fun[[array_number]], fun_args)
         plot_number <<- plot_number + 1
-            })
+      })
     }
     array_number <<- array_number + 1
   })
diff --git a/R/Utils.R b/R/Utils.R
index c13691499e7db59f4247411b0ad79498f9c0a65b..c2c17eb2698c1378b7a4c2147c1002da9b3f9f3a 100644
--- a/R/Utils.R
+++ b/R/Utils.R
@@ -1720,3 +1720,92 @@
   
   return(anom)
 }
+
+#TODO: Remove from s2dv when PlotLayout can get colorbar info from plotting function directly. 
+#      The function is temporarily here because PlotLayout() needs to draw the colorbars of
+#      PlotMostLikelyQuantileMap().
+#Draws Color Bars for Categories
+#A wrapper of s2dv::ColorBar to generate multiple color bars for different 
+#categories, and each category has different color set.
+GradientCatsColorBar <- function(nmap, brks = NULL, cols = NULL, vertical = TRUE, subsampleg = NULL,
+                                 bar_limits, var_limits = NULL,
+                                 triangle_ends = NULL, plot = TRUE,
+                                 draw_separators = FALSE,
+                                 bar_titles = NULL, title_scale = 1, label_scale = 1, extra_margin = rep(0, 4),
+                                 ...) {
+  # bar_limits
+  if (!is.numeric(bar_limits) || length(bar_limits) != 2) {
+    stop("Parameter 'bar_limits' must be a numeric vector of length 2.")
+  }
+
+  # Check brks
+  if (is.null(brks) || (is.numeric(brks) && length(brks) == 1)) {
+    num_brks <- 5
+    if (is.numeric(brks)) {
+      num_brks <- brks
+    }
+    brks <- seq(from = bar_limits[1], to = bar_limits[2], length.out = num_brks)
+  }
+  if (!is.numeric(brks)) {
+    stop("Parameter 'brks' must be a numeric vector.")
+  }
+  # Check cols
+  col_sets <- list(c("#A1D99B", "#74C476", "#41AB5D", "#238B45"),
+                   c("#6BAED6FF", "#4292C6FF", "#2171B5FF", "#08519CFF"),
+                   c("#FFEDA0FF", "#FED976FF", "#FEB24CFF", "#FD8D3CFF"),
+                   c("#FC4E2AFF", "#E31A1CFF", "#BD0026FF", "#800026FF"),
+                   c("#FCC5C0", "#FA9FB5", "#F768A1", "#DD3497"))
+  if (is.null(cols)) {
+    if (length(col_sets) >= nmap) {
+      chosen_sets <- 1:nmap
+      chosen_sets <- chosen_sets + floor((length(col_sets) - length(chosen_sets)) / 2)
+    } else {
+      chosen_sets <- array(1:length(col_sets), nmap)
+    }
+    cols <- col_sets[chosen_sets]
+  } else {
+    if (!is.list(cols)) {
+      stop("Parameter 'cols' must be a list of character vectors.")
+    }
+    if (!all(sapply(cols, is.character))) {
+      stop("Parameter 'cols' must be a list of character vectors.")
+    }
+    if (length(cols) != nmap) {
+      stop("Parameter 'cols' must be a list of the same length as the number of ",
+           "maps in 'maps'.")
+    }
+  }
+  for (i in 1:length(cols)) {
+    if (length(cols[[i]]) != (length(brks) - 1)) {
+      cols[[i]] <- grDevices::colorRampPalette(cols[[i]])(length(brks) - 1)
+    }
+  }
+
+  # Check bar_titles
+  if (is.null(bar_titles)) {
+    if (nmap == 3) {
+      bar_titles <- c("Below normal (%)", "Normal (%)", "Above normal (%)")
+    } else if (nmap == 5) {
+      bar_titles <- c("Low (%)", "Below normal (%)",
+                         "Normal (%)", "Above normal (%)", "High (%)")
+    } else {
+      bar_titles <- paste0("Cat. ", 1:nmap, " (%)")
+    }
+  }
+
+  if (plot) {
+    for (k in 1:nmap) {
+      s2dv::ColorBar(brks = brks, cols = cols[[k]], vertical = FALSE, subsampleg = subsampleg,
+#                     bar_limits = bar_limits, var_limits = var_limits,
+                     triangle_ends = triangle_ends, plot = TRUE,
+                     draw_separators = draw_separators,
+                     title = bar_titles[[k]], title_scale = title_scale,
+                     label_scale = label_scale, extra_margin = extra_margin)
+    }
+  } else {
+    #TODO: col_inf and col_sup
+    return(list(brks = brks, cols = cols))
+  }
+
+}
+