diff --git a/NAMESPACE b/NAMESPACE
index 0e2aec97d063b1b2138af5c1b637b5f9f54a1b13..393534f4adc87919cc28481f0247247e26567b96 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -23,6 +23,7 @@ export(ConfigRemoveEntry)
 export(ConfigShowDefinitions)
 export(ConfigShowSimilarEntries)
 export(ConfigShowTable)
+export(Consist_Trend)
 export(Corr)
 export(EOF)
 export(Eno)
@@ -46,6 +47,7 @@ export(PlotLayout)
 export(PlotMatrix)
 export(PlotSection)
 export(PlotStereoMap)
+export(PlotVsLTime)
 export(ProbBins)
 export(ProjectField)
 export(REOF)
diff --git a/R/Consist_Trend.R b/R/Consist_Trend.R
new file mode 100644
index 0000000000000000000000000000000000000000..b02aa5fe10cf64b6c574fb2314dd3ad91ab38029
--- /dev/null
+++ b/R/Consist_Trend.R
@@ -0,0 +1,201 @@
+#'Compute trend using only model data for which observations are available
+#'
+#'Compute the linear trend for a time series by least square fitting together
+#'with the associated error interval for both the observational and model data.
+#'The 95\% confidence interval and detrended observational and model data are 
+#'also provided.\cr
+#'The function doesn't do the ensemble mean, so if the input data have the 
+#'member dimension, ensemble mean needs to be computed beforehand.
+#'
+#'@param exp A named numeric array of experimental data, with at least two 
+#'  dimensions 'time_dim' and 'dat_dim'.
+#'@param obs A named numeric array of observational data, same dimensions as  
+#'  parameter 'exp' except along 'dat_dim'.
+#'@param dat_dim A character string indicating the name of the dataset 
+#'  dimensions. If data at some point of 'time_dim' are not complete along 
+#'  'dat_dim' in both 'exp' and 'obs', this point in all 'dat_dim' will be 
+#'  discarded. The default value is 'dataset'.
+#'@param time_dim A character string indicating the name of dimension along  
+#'  which the trend is computed. The default value is 'sdate'.
+#'@param interval A positive numeric indicating the unit length between two 
+#' points along 'time_dim' dimension. The default value is 1.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return
+#'A list containing:
+#'\item{$trend}{
+#'  A numeric array of the trend coefficients of model and observational data 
+#'  with dimensions c(stats = 2, nexp + nobs, the rest dimensions of 'exp' and 
+#'  'obs' except time_dim), where 'nexp' is the length of 'dat_dim' in 'exp' 
+#'  and 'nobs' is the length of 'dat_dim' in 'obs. The 'stats' dimension 
+#'  contains the intercept and the slope.
+#'}
+#'\item{$conf.lower}{
+#'  A numeric array of the lower limit of 95\% confidence interval with 
+#'  dimensions same as $trend. The 'stats' dimension contains the lower 
+#'  confidence level of the intercept and the slope.
+#'}
+#'\item{$conf.upper}{
+#'  A numeric array of the upper limit of 95\% confidence interval with 
+#'  dimensions same as $trend. The 'stats' dimension contains the upper 
+#'  confidence level of the intercept and the slope.
+#'}
+#'\item{$detrended_exp}{
+#'  A numeric array of the detrended model data with the same dimensions as 
+#'  'exp'.
+#'}
+#'\item{$detrended_obs}{
+#'  A numeric array of the detrended observational data with the same 
+#'  dimensions as 'obs'.
+#'}
+#'
+#'@examples
+#'#'# Load sample data as in Load() example:
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'runmean_months <- 12
+#'smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
+#'smooth_ano_obs <- Smoothing(ano_obs, runmeanlen = runmean_months)
+#'dim_to_mean <- 'member'  # average along members
+#'years_between_startdates <- 5
+#'trend <- Consist_Trend(MeanDims(smooth_ano_exp, dim_to_mean, na.rm = TRUE), 
+#'                       MeanDims(smooth_ano_obs, dim_to_mean, na.rm = TRUE), 
+#'                       interval = years_between_startdates)
+#'#Bind data for plotting
+#'trend_bind <- abind::abind(trend$conf.lower[2, , ], trend$trend[2, , ], 
+#'                           trend$conf.upper[2, , ], trend$trend[1, , ], along = 0)
+#'trend_bind <- Reorder(trend_bind, c(2, 1, 3))
+#'\donttest{
+#'PlotVsLTime(trend_bind, toptitle = "trend", ytitle = "K/(5 years)", 
+#'            monini = 11, limits = c(-0.8, 0.8), listexp = c('CMIP5 IC3'), 
+#'            listobs = c('ERSST'), biglab = FALSE, hlines = c(0)) 
+#'PlotAno(InsertDim(trend$detrended_exp, 2, 1), InsertDim(trend$detrended_obs, 2, 1), 
+#'        startDates, "Detrended tos anomalies", ytitle = 'K', 
+#'        legends = 'ERSST', biglab = FALSE)
+#'}
+#'
+#'@import multiApply
+#'@export
+Consist_Trend <- function(exp, obs, dat_dim = 'dataset', time_dim = 'sdate', interval = 1,
+                          ncores = NULL) {
+  # Check inputs 
+  ## exp and obs
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a numeric array.")
+  }
+  if (is.null(dim(exp)) | is.null(dim(obs))) {
+    stop(paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
+                "containing time_dim and dat_dim."))
+  }
+  if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
+     any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  if(!all(names(dim(exp)) %in% names(dim(obs))) | 
+     !all(names(dim(obs)) %in% names(dim(exp)))) {
+    stop("Parameter 'exp' and 'obs' must have the same dimension names.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) > 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
+    stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## dat_dim
+  if (!is.character(dat_dim)) {
+    stop("Parameter 'dat_dim' must be a character string.")
+  }
+  if (!all(dat_dim %in% names(dim(exp))) | !all(dat_dim %in% names(dim(obs)))) {
+    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  for (i in 1:length(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim[i])]
+    name_obs <- name_obs[-which(name_obs == dat_dim[i])]
+  }
+  if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimension expect 'dat_dim'."))
+  }
+  ## interval
+  if (!is.numeric(interval) | interval <= 0 | length(interval) > 1) {
+    stop("Parameter 'interval' must be a positive number.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+      length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+
+  ###############################
+  # Calculate Consist_Trend
+
+  output <- Apply(list(exp, obs),
+                  target_dims = list(c(dat_dim, time_dim),
+                                     c(dat_dim, time_dim)),
+                  fun = .Consist_Trend,
+                  output_dims = list(trend = c('stats', dat_dim),
+                                     conf.lower = c('stats', dat_dim),
+                                     conf.upper = c('stats', dat_dim),
+                                     detrended_exp = c(dat_dim, time_dim),
+                                     detrended_obs = c(dat_dim, time_dim)),
+                  interval = interval, 
+                  ncores = ncores)
+
+  return(output)
+}
+
+.Consist_Trend <- function(exp, obs, interval = 1) {
+  # exp: [nexp, sdate]
+  # obs: [nobs, sdate]
+  
+  #  Find common points  
+  nan <- apply(exp, 2, mean, na.rm = FALSE) + apply(obs, 2, mean, na.rm = FALSE)  # [sdate]
+  exp[, is.na(nan)] <- NA
+  obs[, is.na(nan)] <- NA  
+
+  #  Compute trends 
+  res_exp <- apply(exp, 1, .Trend, interval = interval, polydeg = 1)
+  res_obs <- apply(obs, 1, .Trend, interval = interval, polydeg = 1)
+  exp_trend <- lapply(res_exp, '[[', 'trend')
+  exp_trend <- do.call(abind::abind, c(exp_trend, along = 2)) # [stats = 2, dat]
+  obs_trend <- lapply(res_obs, '[[', 'trend')
+  obs_trend <- do.call(abind::abind, c(obs_trend, along = 2))
+  # bind along 'dat'
+  res_trend <- abind::abind(exp_trend, obs_trend, along = 2)  # [stats = 2, dat = (nexp + nobs)]
+
+  # Compute conf.lower
+  exp_conf.lower <- lapply(res_exp, '[[', 'conf.lower')
+  exp_conf.lower <- do.call(abind::abind, c(exp_conf.lower, along = 2)) # [stats = 2, dat]
+  obs_conf.lower <- lapply(res_obs, '[[', 'conf.lower')
+  obs_conf.lower <- do.call(abind::abind, c(obs_conf.lower, along = 2))
+  res_conf.lower <- abind::abind(exp_conf.lower, obs_conf.lower, along = 2)
+
+  # Compute conf.upper
+  exp_conf.upper <- lapply(res_exp, '[[', 'conf.upper')
+  exp_conf.upper <- do.call(abind::abind, c(exp_conf.upper, along = 2)) # [stats = 2, dat]
+  obs_conf.upper <- lapply(res_obs, '[[', 'conf.upper')
+  obs_conf.upper <- do.call(abind::abind, c(obs_conf.upper, along = 2))
+  res_conf.upper <- abind::abind(exp_conf.upper, obs_conf.upper, along = 2)
+
+  # Compute detrended
+  exp_detrended <- lapply(res_exp, '[[', 'detrended')
+  exp_detrended <- do.call(abind::abind, c(exp_detrended, along = 0))
+  obs_detrended <- lapply(res_obs, '[[', 'detrended')
+  obs_detrended <- do.call(abind::abind, c(obs_detrended, along = 0))
+
+  return(invisible(list(trend = res_trend,
+                        conf.lower = res_conf.lower, conf.upper = res_conf.upper,
+                        detrended_exp = exp_detrended, detrended_obs = obs_detrended)))
+}
diff --git a/R/PlotVsLTime.R b/R/PlotVsLTime.R
new file mode 100644
index 0000000000000000000000000000000000000000..94c82e0735c61c4ecf0f8186e09fbe681560a8cb
--- /dev/null
+++ b/R/PlotVsLTime.R
@@ -0,0 +1,265 @@
+#'Plot a score along the forecast time with its confidence interval
+#'
+#'Plot the correlation (\code{Corr()}), the root mean square error 
+#'(\code{RMS()}) between the forecast values and their observational 
+#'counterpart, the slope of their trend (\code{Trend()}), the 
+#'InterQuartile range, maximum-mininum, standard deviation or median absolute 
+#'Deviation of the ensemble members (\code{Spread()}), or the ratio between 
+#'the ensemble spread and the RMSE of the ensemble mean (\code{RatioSDRMS()}) 
+#'along the forecast time for all the input experiments on the same figure 
+#'with their confidence intervals.
+#'
+#'@param var Matrix containing any Prediction Score with dimensions:\cr
+#'  (nexp/nmod, 3/4 ,nltime)\cr
+#'  or (nexp/nmod, nobs, 3/4 ,nltime).
+#'@param toptitle Main title, optional.
+#'@param ytitle Title of Y-axis, optional.
+#'@param monini Starting month between 1 and 12. Default = 1.
+#'@param freq 1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.
+#'@param nticks Number of ticks and labels on the x-axis, optional.
+#'@param limits c(lower limit, upper limit): limits of the Y-axis, optional.
+#'@param listexp List of experiment names, optional.
+#'@param listobs List of observation names, optional.
+#'@param biglab TRUE/FALSE for presentation/paper plot. Default = FALSE.
+#'@param hlines c(a,b, ..) Add horizontal black lines at Y-positions a,b, ...\cr
+#'  Default = NULL.
+#'@param leg TRUE/FALSE if legend should be added or not to the plot. 
+#'  Default = TRUE.
+#'@param siglev TRUE/FALSE if significance level should replace confidence 
+#'  interval.\cr
+#'  Default = FALSE.
+#'@param sizetit Multiplicative factor to change title size, optional.
+#'@param show_conf TRUE/FALSE to show/not confidence intervals for input 
+#'  variables.
+#'@param fileout Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+#'  pdf, bmp and tiff. The default value is NULL.
+#'@param width File width, in the units specified in the parameter size_units 
+#'  (inches by default). Takes 8 by default.
+#'@param height File height, in the units specified in the parameter 
+#'  size_units (inches by default). Takes 5 by default.
+#'@param size_units Units of the size of the device (file or window) to plot 
+#'  in. Inches ('in') by default. See ?Devices and the creator function of the 
+#'  corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. See 
+#'  ?Devices and the creator function of the corresponding device.
+#'@param ... Arguments to be passed to the method. Only accepts the following
+#'  graphical parameters:\cr
+#'  adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+#'  csi cxy err family fg fig font font.axis font.lab font.main font.sub 
+#'  lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+#'  smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+#'  For more information about the parameters see `par`.
+#'
+#'@details
+#'Examples of input:\cr
+#'Model and observed output from \code{Load()} then \code{Clim()} then 
+#'\code{Ano()} then \code{Smoothing()}:\cr
+#'(nmod, nmemb, nsdate, nltime) and (nobs, nmemb, nsdate, nltime)\cr
+#'then averaged over the members\cr
+#'\code{Mean1Dim(var_exp/var_obs, posdim = 2)}:\cr
+#'(nmod, nsdate, nltime) and (nobs, nsdate, nltime)\cr
+#'then passed through\cr
+#'  \code{Corr(exp, obs, posloop = 1, poscor = 2)} or\cr
+#'  \code{RMS(exp, obs, posloop = 1, posRMS = 2)}:\cr
+#'  (nmod, nobs, 3, nltime)\cr
+#'would plot the correlations or RMS between each exp & each obs as a function 
+#'of the forecast time.
+#'
+#'@examples
+#'# Load sample data as in Load() example:
+#'example(Load)
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'runmean_months <- 12
+#'smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+#'smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+#'dim_to_mean <- 'member'  # mean along members
+#'required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+#'leadtimes_per_startdate <- 60
+#'corr <- Corr(MeanDims(smooth_ano_exp, dim_to_mean), 
+#'             MeanDims(smooth_ano_obs, dim_to_mean), 
+#'             comp_dim = required_complete_row, 
+#'             limits = c(ceiling((runmean_months + 1) / 2), 
+#'                        leadtimes_per_startdate - floor(runmean_months / 2)))
+#'# Combine corr results for plotting
+#'corr_combine <- abind::abind(corr$conf.lower, corr$corr, corr$conf.upper, corr$p.val, along = 0)
+#'corr_combine <- Reorder(corr_combine, c(2, 3, 1, 4))
+#'\donttest{
+#'PlotVsLTime(corr_combine, toptitle = "correlations", ytitle = "correlation", 
+#'            monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), 
+#'            listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1))
+#'  }
+#'
+#'@importFrom grDevices dev.cur dev.new dev.off
+#'@importFrom stats ts
+#'@export
+PlotVsLTime <- function(var, toptitle = '', ytitle = '', monini = 1, freq = 12, 
+               nticks = NULL, limits = NULL, 
+               listexp = c('exp1', 'exp2', 'exp3'), 
+               listobs = c('obs1', 'obs2', 'obs3'), biglab = FALSE, hlines = NULL, 
+               leg = TRUE, siglev = FALSE, sizetit = 1, show_conf = TRUE,
+               fileout = NULL, 
+               width = 8, height = 5, size_units = 'in', res = 100, ...) {
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "cex.lab", "cex.main", "col", "fin", "lab", "las", "lend", "lty", "lwd", "mai", "mgp", "new", "pin", "ps", "pty")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  # 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 = width, height = height, units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  #
+  #  Get some arguments
+  # ~~~~~~~~~~~~~~~~~~~~
+  #
+  if (length(dim(var)) == 3) {
+    var <- InsertDim(var, posdim = 2, lendim = 1)
+  } else if (length(dim(var)) != 4) {
+    stop("Parameter 'var' should have 3 or 4 dimensions: c(n. exp[, n. obs], 3/4, n. lead-times)")
+  }
+  nleadtime <- dim(var)[4]
+  nexp <- dim(var)[1]
+  nobs <- dim(var)[2]
+  if (is.null(limits) == TRUE) {
+	if (all(is.na(var > 0))) {
+      ll <- ul <- 0
+	} else {
+      ll <- min(var, na.rm = TRUE)
+      ul <- max(var, na.rm = TRUE)
+    }
+    if (biglab) {
+      ul <- ul + 0.4 * (ul - ll)
+    } else {
+      ul <- ul + 0.3 * (ul - ll)
+    }
+  } else {
+    ll <- limits[1]
+    ul <- limits[2]
+  }
+  lastyear <- (monini + (nleadtime - 1) * 12 / freq - 1) %/% 12
+  lastmonth <- (monini + (nleadtime - 1) * 12 / freq - 1) %% 12 + 1
+  empty_ts <- ts(start = c(0000, (monini - 1) %/% (12 / freq) + 1), 
+                 end = c(lastyear, (lastmonth - 1) %/% (12 / freq) + 1),
+                 frequency = freq)
+  empty <- array(dim = length(empty_ts))
+  #
+  #  Define some plot parameters
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #
+  if (is.null(nticks)) {
+    if (biglab) {
+      nticks <- 5
+    } else {
+      nticks <- 10
+    }
+  }
+  labind <- seq(1, nleadtime, max(nleadtime %/% nticks, 1))
+  months <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep",
+              "Oct", "Nov", "Dec")
+  labyear <- ((labind - 1) * 12 / freq + monini - 1) %/% 12
+  labmonth <- months[((labind - 1) * 12 / freq + monini -1 ) %% 12 + 1]
+  for (jx in 1:length(labmonth)) {
+    y2o3dig <- paste("0", as.character(labyear[jx]), sep = "")
+    labmonth[jx] <- paste(labmonth[jx], "\nYr ", substr(y2o3dig, nchar(y2o3dig)
+                    - 1, nchar(y2o3dig)), sep = "")
+  }
+  color <- c("red1", "dodgerblue1", "green1", "orange1", "lightblue1", 
+             "deeppink1", "mediumpurple1", "lightgoldenrod1", "olivedrab1", 
+             "mediumorchid1")
+  type <- c(1, 3, 2, 4)
+  thickness <- array(dim = c(4, 4))
+  thickness[, 1] <- c(1, 2, 1, 1.5)
+  thickness[, 2] <- c(8, 12, 8, 10)
+  thickness[, 3] <- thickness[, 1]
+  thickness[, 4] <- c(4, 6, 4, 5)
+  if (siglev == TRUE) {
+    lines <- c("n", "l", "n", "l")
+  } else {
+    lines <- c("l", "l", "l", "n")
+  }
+  #
+  #  Define plot layout
+  # ~~~~~~~~~~~~~~~~~~~~
+  #
+
+  # 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 = width, height = height)
+  }
+
+  # Load the user parameters
+  par(userArgs)
+  
+  if (biglab) {
+    par(mai = c(1.25, 1.4, 0.5, 1), mgp = c(4, 2.5, 0))
+    par(cex = 1.3, cex.lab = 2, cex.axis = 1.8)
+    cexmain <- 2.2
+    legsize <- 1.5
+  } else {
+    par(mai = c(1, 1.1, 0.5, 0), mgp = c(3, 1.8, 0))
+    par(cex = 1.3, cex.lab = 1.5, cex.axis = 1.1)
+    cexmain <- 1.5
+    legsize <- 1
+  }
+  plot(empty, ylim = c(ll, ul), xlab = "Time (months)", ylab = ytitle, 
+       main = toptitle, cex.main = cexmain*sizetit, axes = FALSE)
+  axis(1, at = labind, labels = labmonth)
+  axis(2)
+  box()
+  if (is.null(hlines) != TRUE) { 
+    for (jy in 1:length(hlines)) {
+      par(new = TRUE) 
+      abline(h = hlines[jy])
+    }
+  }
+  #
+  #  Loop on experimental & observational data
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #                
+  legendnames <- array(dim = nobs * nexp)
+  legendthick <- array(dim = nobs * nexp)
+  legendsty <- array(dim = nobs * nexp)
+  legendcol <- array(dim = nobs * nexp)
+  ind <- 1
+  if (show_conf == TRUE) {
+    start_line <- dim(var)[3]
+    end_line <- 1
+  } else {
+    start_line <- 2
+    end_line <- 2
+  }
+  for (jt in seq(start_line, end_line, -1)) {
+    ind <- 1
+    for (jexp in 1:nexp) {
+      for (jobs in 1:nobs) {
+        par(new = TRUE)
+        plot(var[jexp, jobs, jt, ], type = lines[jt], ylim = c(ll, ul), 
+             col = color[jexp], lty = type[jobs], lwd = thickness[jobs, jt], 
+             ylab = "", xlab = "", axes = FALSE)
+        legendnames[ind] <- paste(listexp[jexp], 'vs', listobs[jobs])
+        legendthick[ind] <- thickness[jobs, 1] * 3
+        legendsty[ind] <- type[jobs]
+        legendcol[ind] <- color[jexp]
+        ind <- ind + 1
+      }
+    }
+  }
+  if (leg) {
+    if (nobs == 1) {
+      legendnames <- listexp[1:nexp]
+    }  
+    legend(1, ul, legendnames, lty = legendsty, lwd = legendthick,
+           col = legendcol, cex = legsize)
+  }
+  
+  # If the graphic was saved to file, close the connection with the device
+  if(!is.null(fileout)) dev.off()
+}
diff --git a/man/Consist_Trend.Rd b/man/Consist_Trend.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2ac7d42655c91de8713d5a2ff7a865b6e5da5aa9
--- /dev/null
+++ b/man/Consist_Trend.Rd
@@ -0,0 +1,100 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Consist_Trend.R
+\name{Consist_Trend}
+\alias{Consist_Trend}
+\title{Compute trend using only model data for which observations are available}
+\usage{
+Consist_Trend(
+  exp,
+  obs,
+  dat_dim = "dataset",
+  time_dim = "sdate",
+  interval = 1,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numeric array of experimental data, with at least two 
+dimensions 'time_dim' and 'dat_dim'.}
+
+\item{obs}{A named numeric array of observational data, same dimensions as
+parameter 'exp' except along 'dat_dim'.}
+
+\item{dat_dim}{A character string indicating the name of the dataset 
+dimensions. If data at some point of 'time_dim' are not complete along 
+'dat_dim' in both 'exp' and 'obs', this point in all 'dat_dim' will be 
+discarded. The default value is 'dataset'.}
+
+\item{time_dim}{A character string indicating the name of dimension along
+which the trend is computed. The default value is 'sdate'.}
+
+\item{interval}{A positive numeric indicating the unit length between two 
+points along 'time_dim' dimension. The default value is 1.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A list containing:
+\item{$trend}{
+ A numeric array of the trend coefficients of model and observational data 
+ with dimensions c(stats = 2, nexp + nobs, the rest dimensions of 'exp' and 
+ 'obs' except time_dim), where 'nexp' is the length of 'dat_dim' in 'exp' 
+ and 'nobs' is the length of 'dat_dim' in 'obs. The 'stats' dimension 
+ contains the intercept and the slope.
+}
+\item{$conf.lower}{
+ A numeric array of the lower limit of 95\% confidence interval with 
+ dimensions same as $trend. The 'stats' dimension contains the lower 
+ confidence level of the intercept and the slope.
+}
+\item{$conf.upper}{
+ A numeric array of the upper limit of 95\% confidence interval with 
+ dimensions same as $trend. The 'stats' dimension contains the upper 
+ confidence level of the intercept and the slope.
+}
+\item{$detrended_exp}{
+ A numeric array of the detrended model data with the same dimensions as 
+ 'exp'.
+}
+\item{$detrended_obs}{
+ A numeric array of the detrended observational data with the same 
+ dimensions as 'obs'.
+}
+}
+\description{
+Compute the linear trend for a time series by least square fitting together
+with the associated error interval for both the observational and model data.
+The 95\% confidence interval and detrended observational and model data are 
+also provided.\cr
+The function doesn't do the ensemble mean, so if the input data have the 
+member dimension, ensemble mean needs to be computed beforehand.
+}
+\examples{
+#'# Load sample data as in Load() example:
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+runmean_months <- 12
+smooth_ano_exp <- Smoothing(ano_exp, runmeanlen = runmean_months)
+smooth_ano_obs <- Smoothing(ano_obs, runmeanlen = runmean_months)
+dim_to_mean <- 'member'  # average along members
+years_between_startdates <- 5
+trend <- Consist_Trend(MeanDims(smooth_ano_exp, dim_to_mean, na.rm = TRUE), 
+                      MeanDims(smooth_ano_obs, dim_to_mean, na.rm = TRUE), 
+                      interval = years_between_startdates)
+#Bind data for plotting
+trend_bind <- abind::abind(trend$conf.lower[2, , ], trend$trend[2, , ], 
+                          trend$conf.upper[2, , ], trend$trend[1, , ], along = 0)
+trend_bind <- Reorder(trend_bind, c(2, 1, 3))
+\donttest{
+PlotVsLTime(trend_bind, toptitle = "trend", ytitle = "K/(5 years)", 
+           monini = 11, limits = c(-0.8, 0.8), listexp = c('CMIP5 IC3'), 
+           listobs = c('ERSST'), biglab = FALSE, hlines = c(0)) 
+PlotAno(InsertDim(trend$detrended_exp, 2, 1), InsertDim(trend$detrended_obs, 2, 1), 
+       startDates, "Detrended tos anomalies", ytitle = 'K', 
+       legends = 'ERSST', biglab = FALSE)
+}
+
+}
diff --git a/man/PlotVsLTime.Rd b/man/PlotVsLTime.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..05e2b422189793b9c359177ae3ef77ae34cc6ad4
--- /dev/null
+++ b/man/PlotVsLTime.Rd
@@ -0,0 +1,144 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotVsLTime.R
+\name{PlotVsLTime}
+\alias{PlotVsLTime}
+\title{Plot a score along the forecast time with its confidence interval}
+\usage{
+PlotVsLTime(
+  var,
+  toptitle = "",
+  ytitle = "",
+  monini = 1,
+  freq = 12,
+  nticks = NULL,
+  limits = NULL,
+  listexp = c("exp1", "exp2", "exp3"),
+  listobs = c("obs1", "obs2", "obs3"),
+  biglab = FALSE,
+  hlines = NULL,
+  leg = TRUE,
+  siglev = FALSE,
+  sizetit = 1,
+  show_conf = TRUE,
+  fileout = NULL,
+  width = 8,
+  height = 5,
+  size_units = "in",
+  res = 100,
+  ...
+)
+}
+\arguments{
+\item{var}{Matrix containing any Prediction Score with dimensions:\cr
+(nexp/nmod, 3/4 ,nltime)\cr
+or (nexp/nmod, nobs, 3/4 ,nltime).}
+
+\item{toptitle}{Main title, optional.}
+
+\item{ytitle}{Title of Y-axis, optional.}
+
+\item{monini}{Starting month between 1 and 12. Default = 1.}
+
+\item{freq}{1 = yearly, 12 = monthly, 4 = seasonal, ... Default = 12.}
+
+\item{nticks}{Number of ticks and labels on the x-axis, optional.}
+
+\item{limits}{c(lower limit, upper limit): limits of the Y-axis, optional.}
+
+\item{listexp}{List of experiment names, optional.}
+
+\item{listobs}{List of observation names, optional.}
+
+\item{biglab}{TRUE/FALSE for presentation/paper plot. Default = FALSE.}
+
+\item{hlines}{c(a,b, ..) Add horizontal black lines at Y-positions a,b, ...\cr
+Default = NULL.}
+
+\item{leg}{TRUE/FALSE if legend should be added or not to the plot. 
+Default = TRUE.}
+
+\item{siglev}{TRUE/FALSE if significance level should replace confidence 
+interval.\cr
+Default = FALSE.}
+
+\item{sizetit}{Multiplicative factor to change title size, optional.}
+
+\item{show_conf}{TRUE/FALSE to show/not confidence intervals for input 
+variables.}
+
+\item{fileout}{Name of output file. Extensions allowed: eps/ps, jpeg, png, 
+pdf, bmp and tiff. The default value is NULL.}
+
+\item{width}{File width, in the units specified in the parameter size_units 
+(inches by default). Takes 8 by default.}
+
+\item{height}{File height, in the units specified in the parameter 
+size_units (inches by default). Takes 5 by default.}
+
+\item{size_units}{Units of the size of the device (file or window) to plot 
+in. Inches ('in') by default. See ?Devices and the creator function of the 
+corresponding device.}
+
+\item{res}{Resolution of the device (file or window) to plot in. See 
+?Devices and the creator function of the corresponding device.}
+
+\item{...}{Arguments to be passed to the method. Only accepts the following
+graphical parameters:\cr
+adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+csi cxy err family fg fig font font.axis font.lab font.main font.sub 
+lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page pch plt 
+smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog \cr
+For more information about the parameters see `par`.}
+}
+\description{
+Plot the correlation (\code{Corr()}), the root mean square error 
+(\code{RMS()}) between the forecast values and their observational 
+counterpart, the slope of their trend (\code{Trend()}), the 
+InterQuartile range, maximum-mininum, standard deviation or median absolute 
+Deviation of the ensemble members (\code{Spread()}), or the ratio between 
+the ensemble spread and the RMSE of the ensemble mean (\code{RatioSDRMS()}) 
+along the forecast time for all the input experiments on the same figure 
+with their confidence intervals.
+}
+\details{
+Examples of input:\cr
+Model and observed output from \code{Load()} then \code{Clim()} then 
+\code{Ano()} then \code{Smoothing()}:\cr
+(nmod, nmemb, nsdate, nltime) and (nobs, nmemb, nsdate, nltime)\cr
+then averaged over the members\cr
+\code{Mean1Dim(var_exp/var_obs, posdim = 2)}:\cr
+(nmod, nsdate, nltime) and (nobs, nsdate, nltime)\cr
+then passed through\cr
+ \code{Corr(exp, obs, posloop = 1, poscor = 2)} or\cr
+ \code{RMS(exp, obs, posloop = 1, posRMS = 2)}:\cr
+ (nmod, nobs, 3, nltime)\cr
+would plot the correlations or RMS between each exp & each obs as a function 
+of the forecast time.
+}
+\examples{
+# Load sample data as in Load() example:
+example(Load)
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+runmean_months <- 12
+smooth_ano_exp <- Smoothing(data = ano_exp, runmeanlen = runmean_months)
+smooth_ano_obs <- Smoothing(data = ano_obs, runmeanlen = runmean_months)
+dim_to_mean <- 'member'  # mean along members
+required_complete_row <- 'ftime'  # discard startdates for which there are NA leadtimes
+leadtimes_per_startdate <- 60
+corr <- Corr(MeanDims(smooth_ano_exp, dim_to_mean), 
+            MeanDims(smooth_ano_obs, dim_to_mean), 
+            comp_dim = required_complete_row, 
+            limits = c(ceiling((runmean_months + 1) / 2), 
+                       leadtimes_per_startdate - floor(runmean_months / 2)))
+# Combine corr results for plotting
+corr_combine <- abind::abind(corr$conf.lower, corr$corr, corr$conf.upper, corr$p.val, along = 0)
+corr_combine <- Reorder(corr_combine, c(2, 3, 1, 4))
+\donttest{
+PlotVsLTime(corr_combine, toptitle = "correlations", ytitle = "correlation", 
+           monini = 11, limits = c(-1, 2), listexp = c('CMIP5 IC3'), 
+           listobs = c('ERSST'), biglab = FALSE, hlines = c(-1, 0, 1))
+ }
+
+}
diff --git a/tests/testthat/test-Consist_Trend.R b/tests/testthat/test-Consist_Trend.R
new file mode 100644
index 0000000000000000000000000000000000000000..2dd2214a8ea31e851fae8564921c66266d54342c
--- /dev/null
+++ b/tests/testthat/test-Consist_Trend.R
@@ -0,0 +1,166 @@
+context("s2dv::Consist_Trend tests")
+
+##############################################
+  # dat1
+  set.seed(1)
+  exp1 <- array(rnorm(30), dim = c(dataset = 2, sdate = 5, ftime = 3))
+  set.seed(2)
+  obs1 <- array(rnorm(15), dim = c(dataset = 1, sdate = 5, ftime = 3))
+  # dat2
+  exp2 <- exp1
+  set.seed(1)
+  exp2[1, 1, 1] <- NA
+  obs2 <- obs1
+  obs2[1, 2, 3] <- NA
+
+
+##############################################
+test_that("1. Input checks", {
+
+  # exp and obs (1)
+  expect_error(
+  Consist_Trend(c(), c()),
+  "Parameter 'exp' and 'obs' cannot be NULL."
+  )
+  expect_error(
+  Consist_Trend(c('b'), c('a')),
+  "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+  Consist_Trend(c(1:10), c(2:4)),
+  paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
+         "containing time_dim and dat_dim.")
+  )
+  expect_error(
+  Consist_Trend(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
+  "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+  Consist_Trend(array(1:10, dim = c(a = 2, c = 5)), array(1:4, dim = c(a = 2, b = 2))),
+  "Parameter 'exp' and 'obs' must have the same dimension names."
+  )
+  # time_dim
+  expect_error(
+  Consist_Trend(exp1, obs1, time_dim = c('sdate', 'ftime')),
+  "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+  Consist_Trend(exp1, obs1, time_dim = 'asd'),
+  "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # dat_dim
+  expect_error(
+  Consist_Trend(exp1, obs1, dat_dim = c(1, 2)),
+  "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+  Consist_Trend(exp1, obs1, dat_dim = c('member')),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # exp and obs (2)
+  expect_error(
+  Consist_Trend(array(1:10, dim = c(dataset = 2, member = 5, sdate = 4, ftime = 3)),
+       array(1:4, dim = c(dataset = 2, member = 2, sdate = 5, ftime = 3))),
+  paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimension expect 'dat_dim'.")
+  )
+  # interval
+  expect_error(
+  Consist_Trend(exp1, obs1, interval = 0),
+  "Parameter 'interval' must be a positive number."
+  )
+  # ncores
+  expect_error(
+  Consist_Trend(exp1, obs1, ncores = 0),
+  "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+expect_equal(
+names(Consist_Trend(exp1, obs1)),
+c('trend', 'conf.lower', 'conf.upper', 'detrended_exp', 'detrended_obs')
+)
+expect_equal(
+dim(Consist_Trend(exp1, obs1)$trend),
+c(stats = 2, dataset = 3, ftime = 3)
+)
+expect_equal(
+dim(Consist_Trend(exp1, obs1)$conf.lower),
+c(stats = 2, dataset = 3, ftime = 3)
+)
+expect_equal(
+dim(Consist_Trend(exp1, obs1)$conf.upper),
+c(stats = 2, dataset = 3, ftime = 3)
+)
+expect_equal(
+dim(Consist_Trend(exp1, obs1)$detrended_exp),
+c(dataset = 2, sdate = 5, ftime = 3)
+)
+expect_equal(
+dim(Consist_Trend(exp1, obs1)$detrended_obs),
+c(dataset = 1, sdate = 5, ftime = 3)
+)
+expect_equal(
+Consist_Trend(exp1, obs1)$trend[, 2, 1],
+c(0.8287843, -0.1835020),
+tolerance = 0.0001
+)
+expect_equal(
+Consist_Trend(exp1, obs1)$conf.lower[, 2, 2],
+c(-5.449028, -0.943639),
+tolerance = 0.0001 
+)
+expect_equal(
+Consist_Trend(exp1, obs1)$conf.upper[, 2, 2],
+c(3.176215, 1.656969),
+tolerance = 0.0001
+)
+expect_equal(
+Consist_Trend(exp1, obs1)$detrended_exp[, 2, 1],
+c(-0.449003, 1.133500),
+tolerance = 0.0001
+)
+expect_equal(
+Consist_Trend(exp1, obs1)$detrended_obs[, 2, 1],
+c(0.2836287),
+tolerance = 0.0001
+)
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+expect_equal(
+Consist_Trend(exp2, obs2)$trend[, 2, 1],
+c(1.7520623, -0.4143214),
+tolerance = 0.0001
+)
+expect_equal(
+Consist_Trend(exp2, obs2)$detrended_exp[1, , 1],
+c(NA, -0.3160783, 0.4098429, 0.1285491, -0.2223137),
+tolerance = 0.0001
+)
+expect_equal(
+Consist_Trend(exp2, obs2)$detrended_obs[1, , 1],
+c(NA, -0.4826962, 1.2716524, -1.0952163, 0.3062600),
+tolerance = 0.0001
+)
+expect_equal(
+mean(Consist_Trend(exp2, obs2)$detrended_obs, na.rm = TRUE)*10^18,
+2.118364,
+tolerance = 0.0001
+)
+expect_equal(
+mean(Consist_Trend(exp2, obs2)$trend, na.rm = TRUE),
+0.1662461,
+tolerance = 0.0001
+)
+
+})
+
+