diff --git a/NAMESPACE b/NAMESPACE
index 627bf51ce445554a9aa98ee491eb9163808fc2dc..aab3df2ebc5e6f91664caaafddca495b1d08b86c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(AnimateMap)
+export(Ano)
 export(Clim)
 export(ColorBar)
 export(ConfigAddEntry)
@@ -21,6 +22,7 @@ export(InsertDim)
 export(LeapYear)
 export(Load)
 export(MeanDims)
+export(PlotAno)
 export(PlotClim)
 export(PlotEquiMap)
 export(PlotLayout)
@@ -32,6 +34,7 @@ export(RMSSS)
 export(Regression)
 export(Reorder)
 export(Season)
+export(Smoothing)
 export(ToyModel)
 export(Trend)
 export(clim.colors)
@@ -46,6 +49,7 @@ import(methods)
 import(multiApply)
 import(ncdf4)
 import(parallel)
+import(plyr)
 importFrom(ClimProjDiags,Subset)
 importFrom(abind,abind)
 importFrom(abind,adrop)
diff --git a/NEWS.md b/NEWS.md
index 2dd1570b127f4cf597633b8adbf810e47e82d485..d89c1dd7311175acd32c214396dcf7d32cbf27b2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,4 +1,6 @@
 # s2dv 0.0.2 (Release date: 2020-)
+- New functions: Ano(), Composite(), PlotAno(), and Smoothing().
+- Change the default value of PlotClim() parameter 'fileout' to NULL.
 - Change Regression() parameter 'time_dim' to 'reg_dim', and enable the inputs to be vectors.
 - Change Trend() parameter 'time_dim' default value from 'sdate' to 'ftime'.
 - Change the default of Season() parameter 'time_dim' from 'sdate' to 'ftime'.
diff --git a/R/Ano.R b/R/Ano.R
new file mode 100644
index 0000000000000000000000000000000000000000..75a3edfef98ec1876c42541286532a444b3bd4d6
--- /dev/null
+++ b/R/Ano.R
@@ -0,0 +1,99 @@
+#'Compute forecast or observation anomalies
+#'
+#'This function computes anomalies from a multidimensional data array and a 
+#'climatology array. 
+#'
+#'@param data  A numeric array with named dimensions, representing the model or
+#'  observational data to be calculated the anomalies. It should involve all
+#'  the dimensions in parameter 'clim', and it can have more other dimensions.
+#'@param clim A numeric array with named dimensions, representing the 
+#'  climatologies to be deducted from parameter 'data'. It can be generated by 
+#'  Clim(). The dimensions should all be involved in parameter 'data' with the
+#'  same length.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return An array with same dimensions as parameter 'data' but with different 
+#'  dimension order. The dimensions in parameter 'clim' are ordered first.
+#'
+#'@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)
+#'ano_exp <- Reorder(ano_exp, c(1, 2, 4, 3))
+#'ano_obs <- Reorder(ano_obs, c(1, 2, 4, 3))
+#'\donttest{
+#'PlotAno(ano_exp, ano_obs, startDates, 
+#'        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
+#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+#'}
+#'@import multiApply
+#'@export
+  Ano <- function(data, clim, ncores = NULL) {
+
+  # Check inputs 
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be a numeric array.")
+  }
+  if (is.null(dim(data))) {  #is vector
+    dim(data) <- c(length(data))
+    names(dim(data)) <- 'tmp_name'
+  }
+  if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
+    stop("Parameter 'data' must have dimension names.")
+  }
+  ## clim
+  if (is.null(clim)) {
+    stop("Parameter 'clim' cannot be NULL.")
+  }
+  if (!is.numeric(clim)) {
+    stop("Parameter 'clim' must be a numeric array.")
+  }
+  if (is.null(dim(clim))) {  #is vector
+    dim(clim) <- c(length(clim))
+    names(dim(clim)) <- 'tmp_name'
+  }
+  if (any(is.null(names(dim(clim))))| any(nchar(names(dim(clim))) == 0)) {
+    stop("Parameter 'clim' must have dimension names.")
+  }
+  for (i in 1:length(dim(clim))) {
+    if (!(names(dim(clim))[i] %in% names(dim(data)))) {
+      stop("Parameter 'data' must have all the dimensions of parameter 'clim'.")
+    } else {
+      pos <- names(dim(data))[which(names(dim(clim))[i] == names(dim(data)))]
+      if ((dim(clim)[i] != dim(data)[pos])) {
+        stop("Some dimensions of parameter 'clim' have different length from parameter 'data'.")
+      }
+    }
+  }
+  ## 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 Ano
+
+  res <- Apply(list(data),
+               target_dims = names(dim(clim)),
+               output_dims = names(dim(clim)),
+               fun = .Ano,
+               clim = clim,
+               ncores = ncores)$output1
+
+    return(invisible(res))
+  }
+
+  .Ano <- function(data, clim) {
+    ano <- data - clim
+    return(ano)
+  }
diff --git a/R/PlotAno.R b/R/PlotAno.R
new file mode 100644
index 0000000000000000000000000000000000000000..08a9ab88ceff4bfbcd3130203f42e031f0a3ceed
--- /dev/null
+++ b/R/PlotAno.R
@@ -0,0 +1,297 @@
+#'Plot Anomaly time series
+#'
+#'Plots time series of raw or smoothed anomalies of any variable output from 
+#'\code{Load()} or \code{Ano()} or or \code{Ano_CrossValid()} or 
+#'\code{Smoothing()}.
+#'
+#'@param exp_ano A numerical array containing the experimental data:\cr
+#'  c(nmod/nexp, nmemb/nparam, nsdates, nltime).
+#'@param obs_ano A numerical array containing the observational data:\cr
+#'  c(nobs, nmemb, nsdates, nltime)
+#'@param sdates A character vector of start dates in the format of 
+#'  c('YYYYMMDD','YYYYMMDD').
+#'@param toptitle Main title for each experiment: c('',''), optional.
+#'@param ytitle Title of Y-axis for each experiment: c('',''), optional.
+#'@param limits c(lower limit, upper limit): limits of the Y-axis, optional.
+#'@param legends List of observational dataset names, optional.
+#'@param freq 1 = yearly, 12 = monthly, 4 = seasonal, ... Default: 12.
+#'@param biglab TRUE/FALSE for presentation/paper plot. Default = FALSE.
+#'@param fill TRUE/FALSE if the spread between members should be filled. 
+#'  Default = TRUE.
+#'@param memb TRUE/FALSE if all members/only the ensemble-mean should be 
+#'  plotted.\cr
+#'  Default = TRUE.
+#'@param ensmean TRUE/FALSE if the ensemble-mean should be plotted. 
+#'  Default = TRUE.
+#'@param linezero TRUE/FALSE if a line at y=0 should be added. 
+#'  Default = FALSE.
+#'@param points TRUE/FALSE if points instead of lines should be shown. 
+#'  Default = FALSE.
+#'@param vlines List of x location where to add vertical black lines, optional.
+#'@param sizetit Multiplicative factor to scale title size, optional.
+#'@param fileout Name of the output file for each experiment: c('',''). 
+#'  Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff. If filenames 
+#'  with different extensions are passed, it will be considered only the first 
+#'  one and it will be extended to the rest. The default value is NULL, which
+#'  the pop-up window shows.
+#'@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 \dots 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 lend 
+#'  lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page 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`.
+#'
+#'@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)
+#'smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
+#'smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
+#'smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
+#'smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
+#'PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
+#'        toptitle = paste('smoothed anomalies'), ytitle = c('K', 'K', 'K'), 
+#'        legends = 'ERSST', biglab = FALSE)
+#'
+#'@importFrom grDevices dev.cur dev.new dev.off 
+#'@importFrom stats ts
+#'@export
+PlotAno <- function(exp_ano, obs_ano = NULL, sdates, toptitle = rep('', 15),
+                    ytitle = rep('', 15), limits = NULL, legends = NULL, 
+                    freq = 12, biglab = FALSE, fill = TRUE, memb = TRUE, 
+                    ensmean = TRUE, linezero = FALSE, points = FALSE, 
+                    vlines = NULL, sizetit = 1, 
+                    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", "lty", "lwd", "mai", "mgp", "new", "pch", "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(exp_ano)) != 4 ) {
+    stop("4 dim needed : c(nexp/nobs, nmemb, nsdates, nltime)")
+  }
+  nexp <- dim(exp_ano)[1]
+  nmemb <- dim(exp_ano)[2]
+  nleadtime <- dim(exp_ano)[4]
+  nsdates <- dim(exp_ano)[3]
+  if (is.null(obs_ano) == FALSE) { 
+    nobs <- dim(obs_ano)[1] 
+    if (length(dim(obs_ano)) != 4 ) {
+      stop("4 dim needed : c(nexp/nobs, nmemb, nsdates, nltime)") 
+    }
+    if (dim(obs_ano)[3] != nsdates | dim(obs_ano)[4] != nleadtime ) {
+      stop("obs and exp must have same number of sdates & ltimes") 
+    }
+  } else {
+    nobs <- 0
+  }
+  if (is.null(limits) == TRUE) {
+    if (memb) {
+      ll <- min(min(exp_ano, na.rm = TRUE), min(obs_ano, na.rm = TRUE), na.rm = TRUE)
+      ul <- max(max(exp_ano, na.rm = TRUE), max(obs_ano, na.rm = TRUE), na.rm = TRUE)
+    }
+    else{
+      ll <- min(min(MeanDims(exp_ano, 2), na.rm = TRUE), min(obs_ano, na.rm = TRUE), 
+                na.rm = TRUE)
+      ul <- max(max(MeanDims(exp_ano, 2), na.rm = TRUE), max(obs_ano, na.rm = TRUE),
+                na.rm = TRUE)
+    }
+    if (nobs > 0) {
+      if (biglab) {
+        ul <- ul + 0.3 * (ul - ll)
+      } else {
+        ul <- ul + 0.2 * (ul - ll)
+      }
+    }
+  } else {
+    ll <- limits[1]
+    ul <- limits[2]
+  }
+  #
+  #  Define some plot parameters
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  #
+  yearinit <- as.integer(substr(sdates[1], 1, 4))
+  moninit <- as.integer(substr(sdates[1], 5, 6))
+  lastyear <- as.integer(substr(sdates[nsdates], 1, 4)) + (moninit + (
+                         nleadtime - 1) * 12 / freq - 1) %/% 12
+  lastmonth <- (moninit + (nleadtime - 1) * (12 / freq) - 1) %% 12 + 1
+  empty_ts <- ts(start = c(yearinit, (moninit - 1) %/% (12 / freq) + 1), 
+                 end = c(lastyear, (lastmonth - 1) %/% (12 / freq) + 1), 
+                 frequency = freq)
+  color <- c("red4", "orange4", "lightgoldenrod4", "olivedrab4", "green4",
+             "lightblue4", "dodgerblue4", "mediumpurple4", "mediumorchid4",
+             "deeppink4")
+  color <- c(color, color, color, color, color, color, color, color, color, 
+             color, color)
+  colorblock <- c("red1", "orange1", "lightgoldenrod1", "olivedrab1", "green1",
+                  "lightblue1", "dodgerblue1", "mediumpurple1", "mediumorchid1",
+                   "deeppink1")
+  colorblock <- c(colorblock, colorblock, colorblock, colorblock, colorblock,
+                  colorblock, colorblock, colorblock, colorblock, colorblock)
+  type <- c(1, 3, 2, 4)
+  thickness <- c(1, 3, 2, 2)
+  #
+  #  Loop on the experiments : one plot for each
+  # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  # 
+  for (jexp in 1:nexp) {
+    #
+    #  Define plot layout
+    # ~~~~~~~~~~~~~~~~~~~~
+    #
+
+    # Open connection to graphical device
+    if (!is.null(fileout)) {
+      saveToFile(fileout[jexp])
+    } 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, 1.1, 0.5, 0), mgp = c(2.8, 0.9, 0))
+      par(cex = 1.3, cex.lab = 2, cex.axis = 1.8)
+      cexmain <- 2.2
+      legsize <- 1.5
+    } else {
+      par(mai = c(0.8, 0.8, 0.5, 0.3), mgp = c(2, 0.5, 0))
+      par(cex = 1.3, cex.lab = 1.5, cex.axis = 1.1)
+      cexmain <- 1.5
+      legsize <- 1
+    }
+    plot(empty_ts, ylim = c(ll, ul), xlab = "Time (years)", ylab = ytitle[jexp],
+          main = toptitle[jexp], cex.main = cexmain * sizetit)
+    # 
+    #  Plot experimental data + all observational datasets sdate by sdate
+    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    #
+    for (jdate in 1:nsdates) {
+      year0 <- as.integer(substr(sdates[jdate], 1, 4))
+      mon0 <- as.integer(substr(sdates[jdate], 5, 6))
+      start <- (year0 - yearinit) * freq + 1
+      end <- start + nleadtime - 1
+      var <- array(dim = c(nmemb, length(empty_ts)))
+      var[, start:end] <- exp_ano[jexp, , jdate, ]
+      #
+      #  Compute parameters for filling max-min over members 
+      # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~    
+      #
+      if (fill) { 
+        par(new = TRUE)
+        bordup <- array(dim = nleadtime)
+        borddown <- array(dim = nleadtime)
+        for (jt in 1:nleadtime) {
+          bordup[jt] <- max(exp_ano[jexp, , jdate, jt], na.rm = TRUE)
+          borddown[jt] <- min(exp_ano[jexp, , jdate, jt], na.rm = TRUE)
+        }
+        tmp <- c(start:end)
+        xout <- is.na(bordup + borddown) 
+        tmp <- tmp[which(xout == FALSE)]
+        xx <- c(tmp, rev(tmp))
+        bordup <- bordup[which(xout == FALSE)]
+        borddown <- borddown[which(xout == FALSE)]
+        yy <- c(bordup, rev(borddown))
+        #
+        #  Plotting 
+        # ~~~~~~~~~~
+        # 
+        if (jdate == 1) {
+          matplot(t(var), type = "l", lty = 1, lwd = 1, ylim = c(ll, ul),
+                  col = color[jdate], xlab = "", ylab = "", axes = FALSE)
+        } 
+        # Max-min member range
+        polygon(xx, yy, col = colorblock[jdate], border = NA) 
+      }
+      if (ensmean) {  # Ensemble-mean
+        par(new = TRUE)
+        if (points) {
+          plot(MeanDims(t(var), 2), type = "p", lty = 1, lwd = 4, 
+               ylim = c(ll, ul), col = color[jdate], xlab = "", ylab = "",
+               axes = FALSE)
+        } else {
+          plot(MeanDims(t(var), 2), type = "l", lty = 1, lwd = 4, 
+               ylim = c(ll, ul), col = color[jdate], xlab = "", ylab = "",
+               axes = FALSE)
+        }
+      }
+      if (memb) {
+        par(new = TRUE)  # All members
+        if (points) { 
+          matpoints(t(var), type = "p", lty = 1, lwd = 1, pch = 20, 
+                    ylim = c(ll, ul), col = color[jdate], xlab = "", ylab = "",
+                    axes = FALSE)
+        } else {
+          matplot(t(var), type = "l", lty = 1, lwd = 1, ylim = c(ll, ul),
+                  col = color[jdate], xlab = "", ylab = "", axes = FALSE)
+        }
+      }
+      if (nobs > 0) {   
+        for (jobs in 1:nobs) { 
+          for (jmemb in 1:dim(obs_ano)[2]) {
+            var <- array(dim = length(empty_ts))
+            var[start:end] <- obs_ano[jobs, jmemb, jdate, ]
+            par(new = TRUE)  # Observational datasets
+            if (points) {
+              plot(var, type = "p", lty = 1, lwd = 4, pch = 20, 
+                   ylim = c(ll, ul), col = 1, ylab = "", xlab = "", 
+                   axes = FALSE)
+            } else {
+              plot(var, lty = type[jobs], lwd = thickness[jobs], type = "l",
+                   ylim = c(ll, ul), col = 1, ylab = "", xlab = "", 
+                   axes = FALSE)
+            }
+          }
+        }
+      }
+    }
+    if (linezero) {
+      abline(h = 0, col = "black")
+    }
+    if (is.null(vlines) == FALSE) {
+      for (x in vlines) {
+        abline(v = x, col = "black")
+      }
+    }
+    if (is.null(legends) == FALSE) {
+      if (points) {
+        legend('topleft', legends[1:nobs], lty = 3, lwd = 10, col = 1, 
+               cex = legsize)
+      } else {
+        legend('topleft', ul, legends[1:nobs], lty = type[1:nobs], 
+               lwd = thickness[1:nobs], col = 1, cex = legsize)
+      }
+    }
+
+    # If the graphic was saved to file, close the connection with the device
+    if(!is.null(fileout)) dev.off()
+  }
+}
diff --git a/R/PlotClim.R b/R/PlotClim.R
index 45eb2ce3e81823e0f1c23c7200001714d18d7fe9..aff96e40d1a6addd91f0573be83e456ae44144c6 100644
--- a/R/PlotClim.R
+++ b/R/PlotClim.R
@@ -31,8 +31,8 @@
 #'@param res Resolution of the device (file or window) to plot in. See 
 #'  ?Devices and the creator function of the corresponding device.
 #'@param fileout Name of output file. Extensions allowed: eps/ps, jpeg, png, 
-#'  pdf, bmp and tiff. \cr
-#'  Default = 'output_plotclim.eps'.
+#'  pdf, bmp and tiff. The default value is NULL, which the figure is shown
+#'  in a pop-up window.
 #'@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 
@@ -45,11 +45,9 @@
 #'# Load sample data as in Load() example:
 #'example(Load)
 #'clim <- Clim(sampleData$mod, sampleData$obs)
-#'  \donttest{
 #'PlotClim(clim$clim_exp, clim$clim_obs, toptitle = paste('climatologies'), 
 #'         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-#'         listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
-#'  }
+#'         listobs = c('ERSST'), biglab = FALSE, fileout = NULL)
 #'
 #'@importFrom grDevices dev.cur dev.new dev.off 
 #'@importFrom stats ts
@@ -58,7 +56,7 @@ PlotClim <- function(exp_clim, obs_clim = NULL, toptitle = '', ytitle = '',
                      monini = 1, freq = 12, limits = NULL, 
                      listexp = c('exp1', 'exp2', 'exp3'), 
                      listobs = c('obs1', 'obs2', 'obs3'), biglab = FALSE, 
-                     leg = TRUE, sizetit = 1, fileout = 'output_plotclim.eps',
+                     leg = TRUE, sizetit = 1, 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
diff --git a/R/Smoothing.R b/R/Smoothing.R
new file mode 100644
index 0000000000000000000000000000000000000000..b2b11c7539b49d26b0cc889e252055914f1ac0df
--- /dev/null
+++ b/R/Smoothing.R
@@ -0,0 +1,125 @@
+#'Smooth an array along one dimension
+#'
+#'Smooth an array of any number of dimensions along one dimension.
+#'
+#'@param data A numerical array to be smoothed along one of its dimension 
+#'  (typically the forecast time dimension).
+#'@param runmeanlen An integer indicating the running mean length of sampling 
+#'  units (typically months). The default value is 12.
+#'@param time_dim A character string indicating the name of the dimension to be 
+#'  smoothed along. The default value is 'ftime'.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return A numerical array with the same dimensions as parameter 'data' but 
+#'  the 'time_dim' dimension is moved to the first. The head and tail part which
+#'  do not have enough neighboring data for smoothing is assigned as NA. 
+#'
+#'@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)
+#'smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
+#'smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
+#'smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
+#'smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
+#'  \donttest{
+#'PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
+#'        toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
+#'        fileout = "tos_smoothed_ano.png")
+#'  }
+#'@import plyr multiApply
+#'@export
+Smoothing <- function(data, time_dim = 'ftime', runmeanlen = 12, ncores = NULL) {
+
+  # Check data
+  ## data
+  if (is.null(data)) {
+    stop("Parameter 'data' cannot be NULL.")
+  }
+  if (!is.numeric(data)) {
+    stop("Parameter 'data' must be a numeric array.")
+  }
+  if (is.null(dim(data))) {  #is vector
+    dim(data) <- c(length(data))
+    names(dim(data)) <- time_dim
+  }
+  if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
+    stop("Parameter 'data' must have 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(data))) {
+    stop("Parameter 'time_dim' is not found in 'data' dimension.")
+  }
+  ## runmeanlen
+  if (any(!is.numeric(runmeanlen) | runmeanlen %% 1 != 0 | 
+          runmeanlen <= 0 | length(runmeanlen) > 1)) {
+    stop("Parameter 'runmeanlen' must be a positive integer.")
+  }
+  time_dim_length <- dim(data)[which(names(dim(data)) == time_dim)]
+  if (runmeanlen >= time_dim_length & time_dim_length %% 2 == 0) {
+    stop(paste0("Parameter 'runmeanlen' must be within [1, ", time_dim_length - 1, 
+                "]."))
+  }
+  if (runmeanlen > time_dim_length & time_dim_length %% 2 != 0) {
+    stop(paste0("Parameter 'runmeanlen' must be within [1, ", time_dim_length,
+                "]."))
+  }
+  ## 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 Smoothing
+
+  if (runmeanlen == 1) {
+    output <- Reorder(data, c(time_dim, 
+                              names(dim(data))[-which(names(dim(data)) == time_dim)]))
+  } else {
+
+    output <- Apply(list(data),
+                    target_dims = time_dim,
+                    fun = .Smoothing,
+                    output_dims = time_dim,
+                    time_dim = time_dim, runmeanlen = runmeanlen, 
+                    ncores = ncores)$output1
+  }
+  return(output)
+}
+
+.Smoothing <- function(data, runmeanlen = 12, time_dim = 'ftime') {
+  # data: [time_dim]
+
+  nmr1 <- floor(runmeanlen / 2)
+  nltime <- length(data)
+  smoothed <- rep(NA, length(data))
+
+  # We do a loop for all values which have the complete window. 
+  # Other values are left to NA.
+  # If the window length is even, need to weight the half to the values
+  # at both ends.
+  if ((runmeanlen %% 2) == 0) {
+    for (jtime in (1 + nmr1):(nltime - nmr1)) {
+      # calculate the two edges
+      edge <- (data[jtime - nmr1] + data[jtime + nmr1]) / (runmeanlen * 2)
+      # calculate the complete window
+      smoothed[jtime] <- (sum(data[(jtime - nmr1 + 1):(jtime + nmr1 - 1)]) / runmeanlen) + edge
+    }
+  } else {
+    for (jtime in (1 + nmr1):(nltime - nmr1)) {
+      # calculate the complete window
+      smoothed[jtime] <- sum(data[(jtime - nmr1):(jtime + nmr1)]) / runmeanlen
+    }
+  }
+
+  return(smoothed)
+}
diff --git a/man/Ano.Rd b/man/Ano.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2dd4dead6f14828bf18cc108b050e53c6f3fa98c
--- /dev/null
+++ b/man/Ano.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Ano.R
+\name{Ano}
+\alias{Ano}
+\title{Compute forecast or observation anomalies}
+\usage{
+Ano(data, clim, ncores = NULL)
+}
+\arguments{
+\item{data}{A numeric array with named dimensions, representing the model or
+observational data to be calculated the anomalies. It should involve all
+the dimensions in parameter 'clim', and it can have more other dimensions.}
+
+\item{clim}{A numeric array with named dimensions, representing the 
+climatologies to be deducted from parameter 'data'. It can be generated by 
+Clim(). The dimensions should all be involved in parameter 'data' with the
+same length.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+An array with same dimensions as parameter 'data' but with different 
+ dimension order. The dimensions in parameter 'clim' are ordered first.
+}
+\description{
+This function computes anomalies from a multidimensional data array and a 
+climatology array.
+}
+\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)
+ano_exp <- Reorder(ano_exp, c(1, 2, 4, 3))
+ano_obs <- Reorder(ano_obs, c(1, 2, 4, 3))
+\donttest{
+PlotAno(ano_exp, ano_obs, startDates, 
+       toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
+       legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
+}
+}
+
diff --git a/man/PlotAno.Rd b/man/PlotAno.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..18bf0c903db72df4d7ff929b6f8360be5adf19ef
--- /dev/null
+++ b/man/PlotAno.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotAno.R
+\name{PlotAno}
+\alias{PlotAno}
+\title{Plot Anomaly time series}
+\usage{
+PlotAno(exp_ano, obs_ano = NULL, sdates, toptitle = rep("", 15),
+  ytitle = rep("", 15), limits = NULL, legends = NULL, freq = 12,
+  biglab = FALSE, fill = TRUE, memb = TRUE, ensmean = TRUE,
+  linezero = FALSE, points = FALSE, vlines = NULL, sizetit = 1,
+  fileout = NULL, width = 8, height = 5, size_units = "in", res = 100,
+  ...)
+}
+\arguments{
+\item{exp_ano}{A numerical array containing the experimental data:\cr
+c(nmod/nexp, nmemb/nparam, nsdates, nltime).}
+
+\item{obs_ano}{A numerical array containing the observational data:\cr
+c(nobs, nmemb, nsdates, nltime)}
+
+\item{sdates}{A character vector of start dates in the format of 
+c('YYYYMMDD','YYYYMMDD').}
+
+\item{toptitle}{Main title for each experiment: c('',''), optional.}
+
+\item{ytitle}{Title of Y-axis for each experiment: c('',''), optional.}
+
+\item{limits}{c(lower limit, upper limit): limits of the Y-axis, optional.}
+
+\item{legends}{List of observational dataset names, optional.}
+
+\item{freq}{1 = yearly, 12 = monthly, 4 = seasonal, ... Default: 12.}
+
+\item{biglab}{TRUE/FALSE for presentation/paper plot. Default = FALSE.}
+
+\item{fill}{TRUE/FALSE if the spread between members should be filled. 
+Default = TRUE.}
+
+\item{memb}{TRUE/FALSE if all members/only the ensemble-mean should be 
+plotted.\cr
+Default = TRUE.}
+
+\item{ensmean}{TRUE/FALSE if the ensemble-mean should be plotted. 
+Default = TRUE.}
+
+\item{linezero}{TRUE/FALSE if a line at y=0 should be added. 
+Default = FALSE.}
+
+\item{points}{TRUE/FALSE if points instead of lines should be shown. 
+Default = FALSE.}
+
+\item{vlines}{List of x location where to add vertical black lines, optional.}
+
+\item{sizetit}{Multiplicative factor to scale title size, optional.}
+
+\item{fileout}{Name of the output file for each experiment: c('',''). 
+Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff. If filenames 
+with different extensions are passed, it will be considered only the first 
+one and it will be extended to the rest. The default value is NULL, which
+the pop-up window shows.}
+
+\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{\dots}{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 lend 
+lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page 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{
+Plots time series of raw or smoothed anomalies of any variable output from 
+\code{Load()} or \code{Ano()} or or \code{Ano_CrossValid()} or 
+\code{Smoothing()}.
+}
+\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)
+smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
+smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
+smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
+smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
+PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
+       toptitle = paste('smoothed anomalies'), ytitle = c('K', 'K', 'K'), 
+       legends = 'ERSST', biglab = FALSE)
+
+}
+
diff --git a/man/PlotClim.Rd b/man/PlotClim.Rd
index 35ab17d1028933b98f4178e6773b26d61bee751f..b62ff442b7cb6e82a4f1994a037b8f66518f70d0 100644
--- a/man/PlotClim.Rd
+++ b/man/PlotClim.Rd
@@ -7,8 +7,8 @@
 PlotClim(exp_clim, obs_clim = NULL, toptitle = "", ytitle = "",
   monini = 1, freq = 12, limits = NULL, listexp = c("exp1", "exp2",
   "exp3"), listobs = c("obs1", "obs2", "obs3"), biglab = FALSE,
-  leg = TRUE, sizetit = 1, fileout = "output_plotclim.eps", width = 8,
-  height = 5, size_units = "in", res = 100, ...)
+  leg = TRUE, sizetit = 1, fileout = NULL, width = 8, height = 5,
+  size_units = "in", res = 100, ...)
 }
 \arguments{
 \item{exp_clim}{Matrix containing the experimental data with dimensions:\cr
@@ -39,8 +39,8 @@ c(nobs, nmemb, nltime) or c(nobs, nltime)}
 \item{sizetit}{Multiplicative factor to scale title size, optional.}
 
 \item{fileout}{Name of output file. Extensions allowed: eps/ps, jpeg, png, 
-pdf, bmp and tiff. \cr
-Default = 'output_plotclim.eps'.}
+pdf, bmp and tiff. The default value is NULL, which the figure is shown
+in a pop-up window.}
 
 \item{width}{File width, in the units specified in the parameter size_units 
 (inches by default). Takes 8 by default.}
@@ -74,11 +74,9 @@ c(nobs, nmemb, nltime) or c(nobs, nltime) for the observational data
 # Load sample data as in Load() example:
 example(Load)
 clim <- Clim(sampleData$mod, sampleData$obs)
- \donttest{
 PlotClim(clim$clim_exp, clim$clim_obs, toptitle = paste('climatologies'), 
         ytitle = 'K', monini = 11, listexp = c('CMIP5 IC3'), 
-        listobs = c('ERSST'), biglab = FALSE, fileout = 'tos_clim.eps')
- }
+        listobs = c('ERSST'), biglab = FALSE, fileout = NULL)
 
 }
 
diff --git a/man/Smoothing.Rd b/man/Smoothing.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..ba62ca17eab6aa023119e504b37abdae19157ea7
--- /dev/null
+++ b/man/Smoothing.Rd
@@ -0,0 +1,46 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Smoothing.R
+\name{Smoothing}
+\alias{Smoothing}
+\title{Smooth an array along one dimension}
+\usage{
+Smoothing(data, time_dim = "ftime", runmeanlen = 12, ncores = NULL)
+}
+\arguments{
+\item{data}{A numerical array to be smoothed along one of its dimension 
+(typically the forecast time dimension).}
+
+\item{time_dim}{A character string indicating the name of the dimension to be 
+smoothed along. The default value is 'ftime'.}
+
+\item{runmeanlen}{An integer indicating the running mean length of sampling 
+units (typically months). The default value is 12.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A numerical array with the same dimensions as parameter 'data' but 
+ the 'time_dim' dimension is moved to the first. The head and tail part which
+ do not have enough neighboring data for smoothing is assigned as NA.
+}
+\description{
+Smooth an array of any number of dimensions along one dimension.
+}
+\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)
+smooth_ano_exp <- Smoothing(ano_exp, time_dim = 'ftime', runmeanlen = 12)
+smooth_ano_obs <- Smoothing(ano_obs, time_dim = 'ftime', runmeanlen = 12)
+smooth_ano_exp <- Reorder(smooth_ano_exp, c(2, 3, 4, 1))
+smooth_ano_obs <- Reorder(smooth_ano_obs, c(2, 3, 4, 1))
+ \donttest{
+PlotAno(smooth_ano_exp, smooth_ano_obs, startDates, 
+       toptitle = "Smoothed Mediterranean mean SST", ytitle = "K",
+       fileout = "tos_smoothed_ano.png")
+ }
+}
+
diff --git a/tests/testthat/test-Ano.R b/tests/testthat/test-Ano.R
new file mode 100644
index 0000000000000000000000000000000000000000..4d64ce6cd6e4e37655dd28cfccb529e760f155b5
--- /dev/null
+++ b/tests/testthat/test-Ano.R
@@ -0,0 +1,79 @@
+context("s2dv::Ano test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  dat1 <- array(rnorm(72), c(dat = 1,member = 3, sdate = 4, ftime = 6))
+  set.seed(2)
+  clim1 <- array(rnorm(12), c(dat = 1,member = 3, sdate = 4))
+
+  #dat2
+
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+    Ano(c()),
+    "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    Ano(c(NA, NA)),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    Ano(array(1:10, dim = c(2, 5))),
+    "Parameter 'data' must have dimension names."
+  )
+  expect_error(
+    Ano(dat1, c()),
+    "Parameter 'clim' cannot be NULL."
+  )
+  expect_error(
+    Ano(dat1, c(NA, NA)),
+    "Parameter 'clim' must be a numeric array."
+  )
+  expect_error(
+    Ano(dat1, array(1:10, dim = c(2, 5))),
+    "Parameter 'clim' must have dimension names."
+  )
+  expect_error(
+    Ano(dat1, array(1:10, dim = c(dat = 1, member = 3, sdate = 4, a = 2))),
+    "Parameter 'data' must have all the dimensions of parameter 'clim'."
+  )
+  expect_error(
+    Ano(dat1, array(1:10, dim = c(dat = 1, member = 3, sdate = 2))),
+    "Some dimensions of parameter 'clim' have different length from parameter 'data'."
+  )
+  expect_error(
+    Ano(dat1, clim1, ncore = 3.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+  expect_equal(
+    dim(Ano(dat1, clim1)),
+    c(dat = 1, member = 3, sdate = 4, ftime = 6)
+  )
+  expect_equal(
+    mean(Ano(dat1, clim1)),
+    -0.1434844,
+    tolerance = 0.0001
+  )
+  expect_equal(
+    min(Ano(dat1, clim1)),
+    -2.971104,
+    tolerance = 0.0001
+  )
+  expect_equal(
+    Ano(dat1, clim1)[1, 2, , 4],
+    c(-0.24416258, -0.08427184, 0.79636122, -0.05306879),
+    tolerance = 0.0001
+  )
+
+})
diff --git a/tests/testthat/test-Smoothing.R b/tests/testthat/test-Smoothing.R
new file mode 100644
index 0000000000000000000000000000000000000000..ad3dcb9599a82bc32ffa5e0ddc7caba460adbd15
--- /dev/null
+++ b/tests/testthat/test-Smoothing.R
@@ -0,0 +1,118 @@
+context("s2dv::Smoothing test")
+
+##############################################
+  # dat1
+  set.seed(1)
+  dat1 <- array(rnorm(720), c(dat = 1, sdate = 3, ftime = 120, lon = 2)) 
+
+  # dat2
+  set.seed(2)
+  dat2 <- array(rnorm(363), c(member = 3, time = 121)) 
+
+##############################################
+
+test_that("1. Input checks", {
+
+  expect_error(
+  Smoothing(c()),
+  "Parameter 'data' cannot be NULL."
+  )
+  expect_error(
+    Smoothing(c(NA, NA)),
+    "Parameter 'data' must be a numeric array."
+  )
+  expect_error(
+    Smoothing(array(1:10, dim = c(2, 5))),
+    "Parameter 'data' must have dimension names."
+  )
+  expect_error(
+    Smoothing(dat1, time_dim = 'a'),
+    "Parameter 'time_dim' is not found in 'data' dimension."
+  )
+  expect_error(
+    Smoothing(dat1, time_dim = 3),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    Smoothing(dat1, runmeanlen = 0),
+    "Parameter 'runmeanlen' must be a positive integer."
+  )
+  expect_error(
+    Smoothing(dat1, runmeanlen = 0.5),
+    "Parameter 'runmeanlen' must be a positive integer."
+  )
+  expect_error(
+    Smoothing(dat1, runmeanlen = 120),
+    "Parameter 'runmeanlen' must be within [1, 119].",
+    fixed = TRUE
+  )
+  expect_error(
+    Smoothing(dat2, time_dim = 'time', runmeanlen = 122),
+    "Parameter 'runmeanlen' must be within [1, 121].",
+    fixed = TRUE
+  )
+  expect_error(
+    Smoothing(dat1, ncore = 3.5),
+    "Parameter 'ncores' must be a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+  expect_equal(
+    dim(Smoothing(dat1)),
+    c(ftime = 120, dat = 1, sdate = 3, lon = 2)
+  )
+  expect_equal(
+    mean(Smoothing(dat1), na.rm = TRUE),
+    -0.01902183,
+    tolerance = 0.0001
+  )
+  expect_equal(
+    length(Smoothing(dat1)[which(is.na(Smoothing(dat1)))]),
+    72
+  )
+  expect_equal(
+    Smoothing(dat1)[1:7, 1, 2, 2],
+    c(NA, NA, NA, NA, NA, NA, 0.01114821),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    length(Smoothing(dat1, runmeanlen = 21)[which(is.na(Smoothing(dat1, runmeanlen = 21)))]),
+    120
+  )
+  expect_equal(
+    min(Smoothing(dat1, runmeanlen = 21), na.rm = TRUE),
+    -0.6057508,
+    tolerance = 0.0001
+  )
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+  expect_equal(
+    dim(Smoothing(dat2, time_dim = 'time')),
+    c(time = 121, member = 3)
+  )
+  expect_equal(
+    mean(Smoothing(dat2, time_dim = 'time'),  na.rm = TRUE),
+    0.07041814,
+    tolerance = 0.0001
+  )
+  expect_equal(
+    length(Smoothing(dat2, time_dim = 'time')[which(is.na(Smoothing(dat2, time_dim = 'time')))]),
+    36
+  )
+  expect_equal(
+    min(Smoothing(dat2, time_dim = 'time', runmeanlen = 120), na.rm = TRUE),
+    -0.0246432,
+    tolerance = 0.0001
+  )
+
+})
+
+