diff --git a/NAMESPACE b/NAMESPACE
index 3e35710e2d54618b5a2311a1a03f839de76e59f3..4e8557898444268933dfffceb4f5d7d3c8460183 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -56,6 +56,7 @@ export(PlotBoxWhisker)
 export(PlotClim)
 export(PlotEquiMap)
 export(PlotLayout)
+export(PlotMatrix)
 export(PlotSection)
 export(PlotStereoMap)
 export(PlotVsLTime)
diff --git a/R/PlotMatrix.R b/R/PlotMatrix.R
new file mode 100644
index 0000000000000000000000000000000000000000..32a40c99693519596007710da03ff9163d3f528d
--- /dev/null
+++ b/R/PlotMatrix.R
@@ -0,0 +1,201 @@
+#'Function to convert any numerical table to a grid of coloured squares.
+#'
+#'This function converts a numerical data matrix into a coloured 
+#'grid. It is useful for a slide or article to present tabular results as 
+#'colors instead of numbers.
+#'
+#'@param var A numerical matrix containing the values to be displayed in a 
+#'  colored image.
+#'@param brks A vector of the color bar intervals. The length must be one more 
+#'  than the parameter 'cols'. Use ColorBar() to generate default values.
+#'@param cols A vector of valid color identifiers for color bar. The length
+#'  must be one less than the parameter 'brks'. Use ColorBar() to generate 
+#'  default values.
+#'@param toptitle A string of the title of the grid. Set NULL as default. 
+#'@param title.color A string of valid color identifier to decide the title 
+#'  color. Set "royalblue4" as default.
+#'@param xtitle A string of title of the x-axis. Set NULL as default.
+#'@param ytitle A string of title of the y-axis. Set NULL as default.
+#'@param xlabels A vector of labels of the x-axis. The length must be 
+#'  length of the column of parameter 'var'. Set the sequence from 1 to the 
+#'  length of the column of parameter 'var' as default.
+#'@param xvert A logical value to decide whether to place x-axis labels 
+#'  vertically. Set FALSE as default, which keeps the labels horizontally. 
+#'@param ylabels A vector of labels of the y-axis The length must be 
+#'  length of the row of parameter 'var'. Set the sequence from 1 to the 
+#'  length of the row of parameter 'var' as default.
+#'@param line An integer specifying the distance between the title of the 
+#'  x-axis and the x-axis. Set 3 as default. Adjust if the x-axis labels 
+#'  are long.
+#'@param figure.width A positive number as a ratio adjusting the width of the 
+#'  grids. Set 1 as default.
+#'@param legend A logical value to decide to draw the grid color legend or not. 
+#'  Set TRUE as default.
+#'@param legend.width A number between 0 and 0.5 to adjust the legend width.
+#'  Set 0.15 as default.
+#'@param fileout A string of full directory path and file name indicating where 
+#'  to save the plot. If not specified (default), a graphics device will pop up. 
+#'@param size_units A string indicating the units of the size of the device 
+#'  (file or window) to plot in. Set 'px' as default. See ?Devices and the 
+#'  creator function of the corresponding device.
+#'@param res A positive number indicating resolution of the device (file or window) 
+#'  to plot in. See ?Devices and the creator function of the corresponding device.
+#'@param ... The additional parameters to be passed to function ColorBar() in 
+#'  s2dverification for color legend creation.
+#'@return A figure in popup window by default, or saved to the specified path.
+#'
+#'@examples 
+#'#Example with random data
+#' PlotMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
+#'            cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59',
+#'                       '#e34a33','#b30000', '#7f0000'),
+#'            brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1),
+#'            toptitle = "Mean Absolute Error", 
+#'            xtitle = "Forecast time (month)", ytitle = "Start date",
+#'            xlabels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
+#'                        "Aug", "Sep", "Oct", "Nov", "Dec"))
+#'@importFrom grDevices dev.new dev.off dev.cur 
+#'@export
+PlotMatrix <- function(var, brks = NULL, cols = NULL, 
+                       toptitle = NULL, title.color = "royalblue4", 
+                       xtitle = NULL, ytitle = NULL, xlabels = NULL, xvert = FALSE, 
+                       ylabels = NULL, line = 3, figure.width = 1, legend = TRUE, 
+                       legend.width = 0.15, fileout = NULL, 
+                       size_units = 'px', res = 100, ...) {
+
+  # Check variables:  
+  if (!is.matrix(var))
+    stop("Input values are not a matrix")
+  if (!is.numeric(var))
+    stop("Input values are not always numbers")
+
+  # Build: brks, cols 
+  colorbar <- ColorBar(brks = brks, cols = cols, vertical = FALSE, 
+                       plot = FALSE, ...) 
+  brks <- colorbar$brks
+  cols <- colorbar$cols
+
+  n.cols <- length(cols) ## number of colours
+  n.brks <- length(brks) ## number of intervals
+
+  if (n.brks != n.cols + 1)
+    stop("There must be one break more than the number of colors")
+  ncols <- ncol(var) ## number of columns of the image
+  nrows <- nrow(var) ## number of rows of the image
+  if (ncols < 2)
+    stop("Matrix must have at least two columns")
+  if (nrows < 2)
+    stop("Matrix must have at least two rows")
+  if (!is.null(xlabels) && length(xlabels) != ncols)
+      stop(paste0("The number of x labels must be equal to the number of ",
+                  "columns of the data matrix"))
+  if (!is.null(ylabels) && length(ylabels) != nrows)
+      stop(paste0("The number of y labels must be equal to the number of ",
+                  "rows of the data matrix"))
+  if (!is.numeric(figure.width) || figure.width < 0)
+      stop("figure.width must be a positive number")
+  if (!is.numeric(legend.width) || legend.width < 0 || legend.width > 0.5)
+      stop("legend.width must be a number from 0 to 0.5")
+
+    
+  # If there is any filenames to store the graphics, process them
+  # to select the right device 
+  if (!is.null(fileout)) {
+    deviceInfo <- .SelectDevice(fileout = fileout,
+                                width = 80 * ncols * figure.width, 
+                                height = 80 * nrows,
+                                units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  # Open connection to graphical device
+  if (!is.null(fileout)) {
+    saveToFile(fileout)
+  } else if (names(dev.cur()) == 'null device') {
+    dev.new(units = size_units, res = res,
+            width = 8 * figure.width, height = 5)
+  }
+
+  if (!is.null(fileout)) {
+
+  # Draw empty plot:
+    par(mar = c(4, 4, 1, 0), fig = c(0.1, 1 - legend.width, 0.1, 0.9))
+    plot(1, 1, type = "n", xaxt = "n", yaxt = "n", ylim = c(0.5, nrows + 0.5), 
+         xlim = c(-0.5, ncols - 1 + 0.5), ann = F,  bty = "n")
+    
+  # Add axes titles:
+    label.size <- 1.2 * (max(ncols, nrows) / 10) ^ 0.5  
+    mtext(side = 1, text = xtitle, line = line, cex = label.size, font = 3)
+    mtext(side = 2, text = ytitle, line = 3, cex = label.size, font = 3)
+
+  # Add plot title:
+    if (is.null(title.color)) title.color <- "royalblue4"
+    mtext(side = 3, text = toptitle, cex = 1.75 * (nrows / 10) ^ 0.7, 
+          col = title.color)
+
+  # Add axis labels:
+    axis.size <- (max(ncols, nrows) / 10) ^ 0.3
+    if (is.null(xlabels)) xlabels = 1:ncols
+    if (is.null(ylabels)) ylabels = 1:nrows
+    axis(1, at = seq(0, ncols - 1), las = ifelse(xvert, 2, 1), labels = xlabels, 
+         cex.axis = axis.size, tck = 0, lwd = 0, line = - 1 - (nrows / 10 - 1))  ## Add x-axis labels
+    axis(2, at = seq(1, nrows), las = 1, labels = rev(ylabels), 
+         cex.axis = axis.size, tck = 0, lwd = 0, line = 0.5 - ncols / 10)  ## Add y-axis labels
+
+  } else {
+
+  # Draw empty plot:
+    par(mar = c(4, 4, 1, 0), fig = c(0.1, 1 - legend.width, 0.1, 0.9))
+    plot(1, 1, type = "n", xaxt = "n", yaxt = "n", ylim = c(0.5, nrows + 0.5),
+         xlim = c(-0.5, ncols - 1 + 0.5), ann = F,  bty = "n")
+
+  # Add axes titles:
+    label.size <- 1.2 # * (max(ncols, nrows) / 10) ^ 0.5
+    mtext(side = 1, text = xtitle, line = line, cex = label.size, font = 3)
+    mtext(side = 2, text = ytitle, line = 3, cex = label.size, font = 3)
+
+  # Add plot title:
+    if (is.null(title.color)) title.color <- "royalblue4"
+    mtext(side = 3, text = toptitle, cex = 1.5,  #* (nrows / 10) ^ 0.7, 
+          col = title.color)
+
+  # Add axis labels:
+    axis.size <- 1 #(max(ncols, nrows) / 10) ^ 0.3
+    if (is.null(xlabels)) xlabels = 1:ncols
+    if (is.null(ylabels)) ylabels = 1:nrows
+    axis(1, at = seq(0, ncols - 1), las = ifelse(xvert, 2, 1), labels = xlabels, 
+         cex.axis = axis.size, tck = 0, lwd = 0, line = - 1 - (nrows / 10 - 1))  ## Add x-axis labels
+    axis(2, at = seq(1, nrows), las = 1, labels = rev(ylabels),
+         cex.axis = axis.size, tck = 0, lwd = 0, line = 0.5 - ncols / 10)  ## Add y-axis labels
+
+  }
+
+  # Create an array of colors instead of numbers (it starts all gray):                
+  array.colors <- array("gray", c(nrows, ncols))
+  for (int in n.cols:1) array.colors[var <= brks[int + 1]] <- cols[int]
+
+  # fill with colors the cells in the figure:
+  for (p in 1:nrows) {
+    for (l in 0:(ncols - 1)) {
+      polygon(c(0.5 + l - 1, 0.5 + l - 1, 1.5 + l - 1, 1.5 + l - 1),
+              c(-0.5 + nrows + 1 - p, 0.5 + nrows + 1 - p, 
+                 0.5 + nrows + 1 - p, -0.5 + nrows + 1 - p), 
+                 col = array.colors[p, 1 + l], border = "black")
+    }
+  }
+
+  # Draw color legend:
+  if (legend) {
+    par(fig = c(1 - legend.width - 0.01, 
+                1 - legend.width + legend.width * min(1, 10 / ncols), 
+                0.3, 0.8), new = TRUE)
+    #legend.label.size <- (max(ncols, nrows) / 10) ^ 0.4
+    ColorBar(brks = brks, cols = cols, vertical = TRUE, ...)
+  }
+    
+  # If the graphic was saved to file, close the connection with the device
+  if (!is.null(fileout)) dev.off()
+   invisible(list(brks = brks, cols = cols))
+ 
+}
diff --git a/R/Trend.R b/R/Trend.R
index 15ec7ebe91338be57dbc456994e6942ffbfc8ac0..56819bb7858d8eed13254818303394c627c33906 100644
--- a/R/Trend.R
+++ b/R/Trend.R
@@ -1,37 +1,41 @@
 #'Computes the Trend of the Ensemble Mean
 #'
 #'Computes the trend along the forecast time of the ensemble mean by least 
-#'square fitting, and the associated error interval.\cr
-#'Trend() also provides the time series of the detrended ensemble mean 
-#'forecasts.\cr
-#'The confidence interval relies on a student-T distribution.\cr\cr
-#'.Trend provides the same functionality but taking a matrix ensemble members 
+#'square fitting, and the associated error interval. The trend would be 
+#'provided in number of units per month or year. This function also returns the 
+#'time series of the detrended ensemble mean forecasts. The confidence interval 
+#'relies on a student-T distribution.\cr
+#'.Trend() provides the same functionality but taking a matrix ensemble members 
 #'as input (exp).
 #'
-#'@param var Array of any number of dimensions up to 10.
-#'@param exp M by N matrix of M forecasts from N ensemble members.
-#'@param interval Number of months/years between 2 points along posTR 
-#'  dimension.\cr
-#'  Default = 1.\cr
-#'  The trend would be provided in number of units per month or year.
-#'@param siglev Confidence level for the computation of confidence interval. 
-#'  0.95 by default.
-#'@param conf Whether to compute the confidence levels or not. TRUE by default.
-#'@param posTR Position along which to compute the trend.
+#'@param var An array of any number of dimensions up to 10.
+#'@param exp An M by N matrix representing M forecasts from N ensemble members.
+#'@param interval A number of months/years between 2 points along posTR 
+#'  dimension. Set 1 as default.
+#'@param siglev A numeric value indicating the confidence level for the 
+#'  computation of confidence interval. Set 0.95 as default.
+#'@param conf A logical value indicating whether to compute the confidence 
+#'  levels or not. Set TRUE as default.
+#'@param posTR An integer indicating the position along which to compute the 
+#'  trend.
 #'
 #'@return 
 #'\item{$trend}{
-#'  The intercept and slope coefficients for the least squares fitting of the 
-#'  trend. 
+#'  An array with same dimensions as parameter 'var' except along the posTR 
+#'  dimension, which is replaced by a length 4 (or length 2 if conf = FALSE) 
+#'  dimension, corresponding to the lower limit of the confidence interval 
+#'  (only present if conf = TRUE), the slope, the upper limit of the confidence 
+#'  interval (only present if conf = TRUE), and the intercept.
 #'}
-#'\item{$conf.int}{
-#'  Corresponding to the limits of the \code{siglev}\% confidence interval 
-#'  (only present if \code{conf = TRUE}) for the slope coefficient.
-#'} 
 #'\item{$detrended}{
 #'  Same dimensions as var with linearly detrended var along the posTR 
 #'  dimension.
 #'}
+#'Only in .Trend:
+#'\item{$conf.int}{
+#'  Corresponding to the limits of the \code{siglev}\% confidence interval 
+#'  (only present if \code{conf = TRUE}) for the slope coefficient.
+#'} 
 #'
 #'@keywords datagen
 #'@author History:\cr
diff --git a/man/PlotMatrix.Rd b/man/PlotMatrix.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..18d15c1199b45924e89e4dcb23200f5c77e61960
--- /dev/null
+++ b/man/PlotMatrix.Rd
@@ -0,0 +1,89 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotMatrix.R
+\name{PlotMatrix}
+\alias{PlotMatrix}
+\title{Function to convert any numerical table to a grid of coloured squares.}
+\usage{
+PlotMatrix(var, brks = NULL, cols = NULL, toptitle = NULL,
+  title.color = "royalblue4", xtitle = NULL, ytitle = NULL,
+  xlabels = NULL, xvert = FALSE, ylabels = NULL, line = 3,
+  figure.width = 1, legend = TRUE, legend.width = 0.15, fileout = NULL,
+  size_units = "px", res = 100, ...)
+}
+\arguments{
+\item{var}{A numerical matrix containing the values to be displayed in a 
+colored image.}
+
+\item{brks}{A vector of the color bar intervals. The length must be one more 
+than the parameter 'cols'. Use ColorBar() to generate default values.}
+
+\item{cols}{A vector of valid color identifiers for color bar. The length
+must be one less than the parameter 'brks'. Use ColorBar() to generate 
+default values.}
+
+\item{toptitle}{A string of the title of the grid. Set NULL as default.}
+
+\item{title.color}{A string of valid color identifier to decide the title 
+color. Set "royalblue4" as default.}
+
+\item{xtitle}{A string of title of the x-axis. Set NULL as default.}
+
+\item{ytitle}{A string of title of the y-axis. Set NULL as default.}
+
+\item{xlabels}{A vector of labels of the x-axis. The length must be 
+length of the column of parameter 'var'. Set the sequence from 1 to the 
+length of the column of parameter 'var' as default.}
+
+\item{xvert}{A logical value to decide whether to place x-axis labels 
+vertically. Set FALSE as default, which keeps the labels horizontally.}
+
+\item{ylabels}{A vector of labels of the y-axis The length must be 
+length of the row of parameter 'var'. Set the sequence from 1 to the 
+length of the row of parameter 'var' as default.}
+
+\item{line}{An integer specifying the distance between the title of the 
+x-axis and the x-axis. Set 3 as default. Adjust if the x-axis labels 
+are long.}
+
+\item{figure.width}{A positive number as a ratio adjusting the width of the 
+grids. Set 1 as default.}
+
+\item{legend}{A logical value to decide to draw the grid color legend or not. 
+Set TRUE as default.}
+
+\item{legend.width}{A number between 0 and 0.5 to adjust the legend width.
+Set 0.15 as default.}
+
+\item{fileout}{A string of full directory path and file name indicating where 
+to save the plot. If not specified (default), a graphics device will pop up.}
+
+\item{size_units}{A string indicating the units of the size of the device 
+(file or window) to plot in. Set 'px' as default. See ?Devices and the 
+creator function of the corresponding device.}
+
+\item{res}{A positive number indicating resolution of the device (file or window) 
+to plot in. See ?Devices and the creator function of the corresponding device.}
+
+\item{...}{The additional parameters to be passed to function ColorBar() in 
+s2dverification for color legend creation.}
+}
+\value{
+A figure in popup window by default, or saved to the specified path.
+}
+\description{
+This function converts a numerical data matrix into a coloured 
+grid. It is useful for a slide or article to present tabular results as 
+colors instead of numbers.
+}
+\examples{
+#Example with random data
+PlotMatrix(var = matrix(rnorm(n = 120, mean = 0.3), 10, 12),
+           cols = c('white','#fef0d9','#fdd49e','#fdbb84','#fc8d59',
+                      '#e34a33','#b30000', '#7f0000'),
+           brks = c(-1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1),
+           toptitle = "Mean Absolute Error", 
+           xtitle = "Forecast time (month)", ytitle = "Start date",
+           xlabels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
+                       "Aug", "Sep", "Oct", "Nov", "Dec"))
+}
+
diff --git a/man/Trend.Rd b/man/Trend.Rd
index da1fe9560cb63ebde5aceea25cde13a3db4b98b8..e7035e3dce50dbf9b97d2ce6c8bab1c40ef9b06e 100644
--- a/man/Trend.Rd
+++ b/man/Trend.Rd
@@ -10,43 +10,47 @@ Trend(var, posTR = 2, interval = 1, siglev = 0.95, conf = TRUE)
 .Trend(exp, interval = 1, siglev = 0.95, conf = TRUE)
 }
 \arguments{
-\item{var}{Array of any number of dimensions up to 10.}
+\item{var}{An array of any number of dimensions up to 10.}
 
-\item{posTR}{Position along which to compute the trend.}
+\item{posTR}{An integer indicating the position along which to compute the 
+trend.}
 
-\item{interval}{Number of months/years between 2 points along posTR 
-dimension.\cr
-Default = 1.\cr
-The trend would be provided in number of units per month or year.}
+\item{interval}{A number of months/years between 2 points along posTR 
+dimension. Set 1 as default.}
 
-\item{siglev}{Confidence level for the computation of confidence interval. 
-0.95 by default.}
+\item{siglev}{A numeric value indicating the confidence level for the 
+computation of confidence interval. Set 0.95 as default.}
 
-\item{conf}{Whether to compute the confidence levels or not. TRUE by default.}
+\item{conf}{A logical value indicating whether to compute the confidence 
+levels or not. Set TRUE as default.}
 
-\item{exp}{M by N matrix of M forecasts from N ensemble members.}
+\item{exp}{An M by N matrix representing M forecasts from N ensemble members.}
 }
 \value{
 \item{$trend}{
- The intercept and slope coefficients for the least squares fitting of the 
- trend. 
+ An array with same dimensions as parameter 'var' except along the posTR 
+ dimension, which is replaced by a length 4 (or length 2 if conf = FALSE) 
+ dimension, corresponding to the lower limit of the confidence interval 
+ (only present if conf = TRUE), the slope, the upper limit of the confidence 
+ interval (only present if conf = TRUE), and the intercept.
 }
-\item{$conf.int}{
- Corresponding to the limits of the \code{siglev}\% confidence interval 
- (only present if \code{conf = TRUE}) for the slope coefficient.
-} 
 \item{$detrended}{
  Same dimensions as var with linearly detrended var along the posTR 
  dimension.
 }
+Only in .Trend:
+\item{$conf.int}{
+ Corresponding to the limits of the \code{siglev}\% confidence interval 
+ (only present if \code{conf = TRUE}) for the slope coefficient.
+}
 }
 \description{
 Computes the trend along the forecast time of the ensemble mean by least 
-square fitting, and the associated error interval.\cr
-Trend() also provides the time series of the detrended ensemble mean 
-forecasts.\cr
-The confidence interval relies on a student-T distribution.\cr\cr
-.Trend provides the same functionality but taking a matrix ensemble members 
+square fitting, and the associated error interval. The trend would be 
+provided in number of units per month or year. This function also returns the 
+time series of the detrended ensemble mean forecasts. The confidence interval 
+relies on a student-T distribution.\cr
+.Trend() provides the same functionality but taking a matrix ensemble members 
 as input (exp).
 }
 \examples{