diff --git a/.Rbuildignore b/.Rbuildignore
index b2d8e5fcebca62bff5e5380a881580283874cd54..fa596e707601da63df8c53cf4f087a70a953dbea 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -5,4 +5,4 @@
 ./.nc$
 .*^(?!data)\.RData$
 .*\.gitlab-ci.yml$
-#^tests$
+^tests$
diff --git a/DESCRIPTION b/DESCRIPTION
index 318b382f4467a01310bc67ef2ad5c6e475de24f2..52999f2a9e329e9db763e3e2f5e482329bd4dc7b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -85,4 +85,4 @@ VignetteBuilder: knitr
 License: Apache License 2.0
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 7.0.1
+RoxygenNote: 7.0.2
diff --git a/NEWS.md b/NEWS.md
index 154c630724e79b37f607edd687f8fba9538c4457..6d66f0a29c3e7b22226d91d6ca505b8e9a462961 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,11 +1,13 @@
 ### CSTools 4.0.1
-<<<<<<< HEAD
 **Submission date to CRAN: XX-06-2021**
 
 - New features:
     + Dynamical Bias Correction method: `CST_ProxiesAttractors` and `CST_DynBiasCorrection` 
       (optionally `Predictability`)
     + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'.
+    + Use cases
+    + CST_SaveExp includes parameter extra_string
+    + PlotCombinedMap includes parameter cex_bar_titles
 
 - Fixes:
     + Calibration retains correlation absolute value 
diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R
index 9c689ff7156d825d22cdb1bde8d4a6785c854dbc..f9290b18d739636b66fe61feeabb86e774e9b22e 100644
--- a/R/CST_SaveExp.R
+++ b/R/CST_SaveExp.R
@@ -13,6 +13,7 @@
 #'folder tree: destination/experiment/variable/. By default the function 
 #'creates and saves the data into the folder "CST_Data" in the working 
 #'directory.
+#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters.
 #'
 #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}}
 #'
@@ -29,7 +30,7 @@
 #'}
 #'
 #'@export
-CST_SaveExp <- function(data, destination = "./CST_Data") {
+CST_SaveExp <- function(data, destination = "./CST_Data", extra_string = NULL) {
   if (!is.character(destination) & length(destination) > 1) {
     stop("Parameter 'destination' must be a character string of one element ",
          "indicating the name of the file (including the folder if needed) ",
@@ -55,7 +56,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
   SaveExp(data = data$data, lon = data$lon, lat = data$lat,
           Dataset = names(data$Datasets), var_name = var_name,
           units = units, cdo_grid_name = cdo_grid_name, projection = projection,
-          startdates = sdates, Dates = time_values, destination) 
+          startdates = sdates, Dates = time_values, destination, 
+          extra_string = extra_string) 
 } 
 #'Save an experiment in a format compatible with CST_Load
 #'@description This function is created for compatibility with CST_Load/Load for saving post-processed datasets such as those calibrated of downscaled with CSTools functions
@@ -73,8 +75,9 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
 #'@param cdo_grid_name a character string indicating the name of the grid e.g.: 'r360x181'
 #'@param projection a character string indicating the projection name
 #'@param destination a character string indicating the path where to store the NetCDF files
+#'@param extra_string a character string to be include as part of the file name, for instance, to identify member or realization.
 #'
-#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members
+#'@return the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters.
 #' The path will be created with the name of the variable and each Datasets.
 #' 
 #'@import ncdf4
@@ -102,7 +105,8 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
 #'}
 #'@export
 SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
-                    cdo_grid_name, projection, destination) {  
+                    cdo_grid_name, projection, destination, 
+                    extra_string = NULL) {  
   dimname <- names(dim(data))
   if (any(dimname == "ftime")) {
     dimname[which(dimname == "ftime")] <- "time"
@@ -136,7 +140,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
     stop("Element 'data' in parameter 'data' has more than one 'sdate'",
          " dimension.")
   }
-
+  if (!is.null(extra_string)) {
+    if (!is.character(extra_string)) {
+      stop("Parameter 'extra_string' must be a character string.")
+    }
+  }
   dataset_pos <- which(dimname == "dataset" | dimname == "dat")
   dims <- dim(data)
   if (length(dataset_pos) == 0) {
@@ -227,7 +235,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
                              NULL, 'time'),
           fun = .saveExp, var_name = var_name, units = units,
           dims_var = dims_var, cdo_grid_name = cdo_grid_name, projection = projection,
-          destination = path)
+          destination = path, extra_string = extra_string)
   }
 }
 
@@ -252,7 +260,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
 #Dates <- as.Date(c("1900-11-01", "1900-12-01", "1901-01-01", "1901-02-01", "1901-03-01"))
 #.saveExp(data, sdate, Dates, var_name, units, dims_var, cdo_grid_name = 'r360x181', projection = 'none', destination)
 .saveExp <- function(data, sdate, Dates, var_name, units, dims_var, 
-                     cdo_grid_name, projection, destination) {
+                     cdo_grid_name, projection, destination, extra_string) {
   dim_names <- names(dim(data))
   if (any(dim_names != c('longitude', 'latitude', 'member', 'time'))) {
     data <- Reorder(data, c('longitude', 'latitude', 'member', 'time'))
@@ -266,7 +274,11 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
   datanc  <- ncvar_def(name = var_name,
                        units = units,
                        dim = dims_var, missval = -99999)
-  file_name <- paste0(var_name, "_", sdate, ".nc")
+  if (is.null(extra_string)) {
+    file_name <- paste0(var_name, "_", sdate, ".nc")
+  } else {
+    file_name <- paste0(var_name, "_", extra_string, "_", sdate, ".nc")
+  }
   full_filename <- file.path(destination, file_name)
   file_nc <- nc_create(full_filename, datanc)
   ncvar_put(file_nc, datanc, data)
diff --git a/R/PlotCombinedMap.R b/R/PlotCombinedMap.R
index 7169e4c43a8bab1e9a30accca7433b0702760d12..8af4a14db9f69e258e1eb5bf28961cace1be9def 100644
--- a/R/PlotCombinedMap.R
+++ b/R/PlotCombinedMap.R
@@ -22,6 +22,7 @@
 #'  layers via the parameter 'dot_symbol'.
 #'@param bar_titles Optional vector of character strings providing the titles to be shown on top of each of the colour bars.
 #'@param legend_scale Scale factor for the size of the colour bar labels. Takes 1 by default.
+#'@param cex_bar_titles Scale factor for the sizes of the bar titles. Takes 1.5 by default.
 #'@param fileout File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff
 #'@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.
@@ -69,6 +70,7 @@ PlotCombinedMap <- function(maps, lon, lat,
                             mask = NULL, col_mask = 'grey',
                             dots = NULL,
                             bar_titles = NULL, legend_scale = 1,
+                            cex_bar_titles = 1.5,
                             fileout = NULL, width = 8, height = 5, 
                             size_units = 'in', res = 100, 
                             ...) {
@@ -422,7 +424,7 @@ PlotCombinedMap <- function(maps, lon, lat,
              draw_separators = TRUE, extra_margin = c(2, 0, 2, 0), 
              label_scale = legend_scale * 1.5)
     if (!is.null(bar_titles)) {
-      mtext(bar_titles[[k]], 3, line = -3, cex = 1.5)
+      mtext(bar_titles[[k]], 3, line = -3, cex = cex_bar_titles)
     }
   }
   
diff --git a/inst/doc/UseCase1_WindEvent_March2018.R b/inst/doc/UseCase1_WindEvent_March2018.R
index 2baa597bdfd933b59abc7c1dacb3ff854d1a36a7..d3bbb9363ba0c9058efa573bae6f0a9512163a80 100644
--- a/inst/doc/UseCase1_WindEvent_March2018.R
+++ b/inst/doc/UseCase1_WindEvent_March2018.R
@@ -1,7 +1,22 @@
+#!/usr/bin/env Rscript
+rm(list=ls()); gc();
 
-library(CSTools)
+### Creation date: 3rd June 2021
+# Author: N. Pérez-Zanón
+# Refer CSTools package and manuscript when using it.
+# ----------------------------------------
+# Wind speed bias adjustment for assessment of an extreme event
+# The System5-ECMWF downloaded from C3S seasonal forecasts initialized 
+# in December 2017, January 2018 and February 2018 
+# This code includes the bias adjustent and the results visualization 
+# ----------------------------------------
+
+#library(CSTools)
 library(s2dv)
 library(ragg)
+library(multiApply)
+output_dir <- "/esarchive/scratch/nperez/CSTools_manuscript/v20210603/"
+
 
 exp_path <- list(name = "ECMWFS5", 
                  path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$START_DATE$.nc")
@@ -15,6 +30,7 @@ months_in_advance <- c('02', '01', '12')
 wind_fsct_BC <- list()
 wind_ref_terciles <- NULL
 wind_ref_extremes <- NULL
+wind_thres_latlon <- NULL
   # Observations March 2018
   wind_obs <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path),
                        sdates = '20180301', nmember = 1,
@@ -63,12 +79,13 @@ for (mm in 1:3) {
                         grid = 'r360x181')
   str(wind_ref$Dates)
   dim(wind_ref$data)
+  print(wind_ref$Dates$start)
+                                    
   wind_ref_terciles <- rbind(wind_ref_terciles,
                              quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.3, 0.6)))
   wind_ref_extremes <- rbind(wind_ref_extremes,
                              quantile(MeanDims(wind_ref$data, c('lat', 'lon')), c(0.1, 0.9))) 
-  #source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R")
-  library(multiApply) 
+  source("/esarchive/scratch/nperez/git/cstools/R/CST_BiasCorrection.R")
   wind_fsct <- CST_BiasCorrection(exp = wind_hcst,
                                   obs = wind_ref, 
                                   exp_cor = wind_fcst)
@@ -118,17 +135,20 @@ for (mm in 1:3) {
                 } 
                 return(y) })$output1
 
-
+  wind_thres_latlon <- abind::abind(wind_thres_latlon, thres, along = 4)
   source("/esarchive/scratch/nperez/git/cstools/R/PlotCombinedMap.R")
   source("/esarchive/scratch/nperez/git/cstools/R/PlotMostLikelyQuantileMap.R")
-  agg_png(paste0("/esarchive/scratch/nperez/CSTools_manuscript/Wind/MostLikely_", mm, "_obstercile.png"),
-          width = 1000, height = 1000, units = 'px', res = 144) 
-  PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon, lat = wind_fsct$lat,
+  agg_png(paste0(output_dir, "Wind_MostLikely_", mm, "_obstercile.png"),
+          width = 1050, height = 1000, units = 'px', res = 144) 
+  PlotMostLikelyQuantileMap(probs = Mean_PB, lon = wind_fsct$lon,
+                          lat = wind_fsct$lat, sizetit = 1.5,
                           intylat = 2, intxlon = 2, 
-                          coast_width = 1.5, legend_scale = 0.8, cat_dim = 'bin', 
+                          coast_width = 1.5, legend_scale = 0.8, 
+                          cat_dim = 'bin', dot_size = 2.5,
+                          axes_label_scale = 1.6, degree_sym = TRUE, 
                           dots = filtered_obs_terciles[,,1,1,1,1],
                           toptitle = c(paste0("Initialized on ", 
-                                        month.name[as.numeric(months_in_advance[mm])]))) 
+                              month.name[as.numeric(months_in_advance[mm])]))) 
   dev.off() 
 }
 
@@ -136,7 +156,7 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', '
                      jan = as.vector(MeanDims(wind_fsct_BC[[2]]$data, c('lat', 'lon'))),
                      feb = as.vector(MeanDims(wind_fsct_BC[[1]]$data, c('lat', 'lon'))))
 
-  wind_obs <- CST_Load(var = 'windagl100', obs = list(obs_path),
+  wind_obs_areave <- CSTools::CST_Load(var = 'windagl100', obs = list(obs_path),
                        sdates = '20180301', nmember = 1,
                         leadtimemin = 1, leadtimemax = 1,
                         storefreq = "monthly", sampleperiod = 1,
@@ -147,13 +167,45 @@ visual <- data.frame(dec = as.vector(MeanDims(wind_fsct_BC[[3]]$data, c('lat', '
 print("IS DATA LOADED")
 
 print("Wait")
-agg_png("/esarchive/scratch/nperez/CSTools_manuscript/Wind/PlotForecast_IP.png",
-    width = 1000, height = 1000, units = 'px',res = 144)
-PlotForecastPDF(visual, tercile.limits = wind_ref_terciles,
+agg_png(paste0(output_dir, "Wind_PlotForecast_IP.png"),
+    width = 1000, height = 1000, units = 'px',res = 150)
+CSTools::PlotForecastPDF(visual, tercile.limits = wind_ref_terciles,
                 extreme.limits = wind_ref_extremes,
-                var.name = "Wind Speed 100 m (m/s)", title = "Bias Corrected forecasts at IP", 
-    fcst.names = c("December", "January", "February"), obs = as.vector(wind_obs$data))
+                var.name = "Wind Speed 100 m (m/s)",
+                title = "Bias Corrected forecasts at IP", 
+                fcst.names = c("December", "January", "February"), 
+                obs = as.vector(wind_obs_areave$data))
 dev.off()
 
+# Plotting observed terciles:
+names(dim(wind_thres_latlon)) <- c('thres', 'lat', 'lon', 'sdate')
+wind_thres_latlon <- ClimProjDiags::Subset(wind_thres_latlon, indices = 1, along = 'sdate')
+wind_obs_obstercile <- Apply(list(wind_obs$data, wind_thres_latlon), 
+                             target_dims = list(NULL, 'thres'),
+                             fun = function(x, tercile) {
+                               if (x <= tercile[1]) {
+                                 res <- 1
+                               } else if (x > tercile[2]) {
+                                 res <- 3
+                               } else {
+                                 res <- 2
+                               }
+                               return(res)
+                             })$output1
+                              
+agg_png(paste0(output_dir, "MostLikely_Observed_obstercile.png"),
+        width = 1000, height = 1000, units = 'px', res = 144)
+
+s2dv::PlotEquiMap(wind_obs_obstercile, 
+                  lon = wind_obs$lon, lat = wind_obs$lat,
+                  brks = c(0,1,2,3),
+                  cols = c("#6BAED6FF", "#FFEDA0FF", "#FC4E2AFF"),
+                          intylat = 2, intxlon = 2,
+                          coast_width = 1.5, filled.continents = FALSE,
+                          toptitle = c("Observed terciles March 2018"))
+dev.off()
+
+# All gridpoints are above normal observed tercile.
+
 print("DONE")
  
diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R
new file mode 100644
index 0000000000000000000000000000000000000000..146382dc8983caa3873a8b13665214b6d0dfb7d0
--- /dev/null
+++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF100.R
@@ -0,0 +1,160 @@
+#!/usr/bin/env Rscript
+rm(list=ls()); gc();
+
+### Creation date: 9th June 2021
+# Author: N. Pérez-Zanón
+# Adapted from the original version Terzago et al. 2020
+# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view
+# Refer CSTools package and manuscript when using it.
+# ----------------------------------------
+# This code is using divided in 4 steps plus visualization of the results
+# Taking advantage of 
+# STEP 1: Compute Slope values for RainFARM downscaling method
+# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast
+# STEP 3: Compute Weights for RainFARM downscaling method
+# STEP 4: Apply RainFARM downscaling method
+# ----------------------------------------
+# Note: STEP 3 requires a machine with internet connexion. 
+#       In this file, the lines are commented since they have been run and the 
+#       result saved on disk, then the result is loaded.
+# ----------------------------------------
+#
+# Load required libraries and setup output directory:
+library(CSTools)
+library(ClimProjDiags)
+library(zeallot)
+library(ragg)
+dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end
+
+# --------------------------------------------
+# STEP 1:
+# --------------------------------------------
+era5 <- list(name = 'era5',
+             path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+years <- unlist(lapply(1993:2018, function(x){
+                       paste0(x, sprintf("%02d",1:12), '01')}))
+era5 <- CST_Load(var = 'prlr',
+                 exp = list(era5),
+                 sdates = years, nmember = 1,
+                 storefreq = "daily", sampleperiod = 1,
+                 latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25,
+                 output = 'lonlat', nprocs = 1)
+era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result
+era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26))
+slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5)
+save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2)
+slope_plot <- slope
+
+# --------------------------------------------
+# STEP 2:
+# --------------------------------------------
+StartDates <- paste0(1993:2018, '1101')
+exp <- list(name = 'ecmwfS5',
+            path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")
+obs <- list(name = 'era5',
+            path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+
+c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs),
+                          sdates = StartDates, nmember = 25,
+                          storefreq = "daily", sampleperiod = 1,
+                          latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+                          output = 'lonlat', nprocs = 1)
+exp <- CST_SplitDim(exp, split_dim = c('ftime'))
+obs <- CST_SplitDim(obs, split_dim = c('ftime'))
+exp$data <- exp$data * 24 * 3600 * 1000
+obs$data <- obs$data * 24 * 3600 * 1000
+exp$data[which(exp$data < 0)] <- 0
+exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT",
+                              wet.day = FALSE,
+                              sample_dims = c('member', 'sdate', 'ftime'),
+                              ncores = 4)
+save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2)
+
+# --------------------------------------------
+# STEP 3:
+# --------------------------------------------
+#library(raster); library(s2dv); library(CSTools)
+#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45)
+#wc_month <- lapply(1:12, FUN = function(x) {
+#                   res <- crop(worldclim[[x]],
+#                               extent(3.5, 11.5, 41.5, 49.5))
+#                   res <- as.array(res)
+#                   names(dim(res)) <- c('lat', 'lon', 'month')
+#                   return(res)
+#                   })
+#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)),
+#                 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5))))
+#lons <- unique(xy[,1])
+#lats <- unique(xy[,2])
+#wc_month <- unlist(wc_month)
+#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12)
+#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly'))
+#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats,
+#                       Datasets = 'WorldClim')
+#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 100)
+#save(weight, file = paste0(dir_output, 'weightsRF100.RDS'))
+
+load(paste0(dir_output, 'weightsRF100.RDS'))
+
+# --------------------------------------------
+# Visualization
+# --------------------------------------------
+agg_png(paste0(dir_output, "RF100_WeightsDec.png"),
+        width = 1000, height = 1100, units = 'px',res = 144)
+PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat,
+            filled.continents = FALSE, title_scale = 1,
+            intylat = 2, intxlon = 2,
+            toptitle = 'December Weights RF 100')
+dev.off()
+
+# --------------------------------------------
+# STEP 4:
+# --------------------------------------------
+
+weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6))
+slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), 
+                drop = 'non-selected')
+k = 1 # To create the member ID
+for (realizations in 1:10) {
+  for (member in 1:25) {
+    result <- exp.qm # to store the data
+    result$data <- NULL
+    for (month in 1:8) {
+      data <- exp.qm # to take the correct piece of data
+      data$data <- data$data[1, member, , , , , month]
+      fs <- CST_RainFARM(data, nf = 100,
+                        weights = weights, slope = slope[month],
+                        kmin = 1, nens = 1, verbose = TRUE,
+                        nprocs = 8,
+                        drop_realization = TRUE)
+      result$data <- abind::abind(result$data, fs$data, along = 5)
+      if (month == 2 & member == 1 & realization == 1) {
+        # ----------------------------
+        # Visualization:
+        # ----------------------------
+        agg_png(paste0(dir_output, "RF100_Down_11dec.png"),
+                width = 1000, height = 1100, units = 'px',res = 144)
+        PlotEquiMap(fs$data[1,11,,],lon = fs$lon, lat = fs$lat,
+                    filled.continents = FALSE, bar_limits = c(0,40),
+                    intylat = 2, intxlon = 2, title_scale = 1,
+                    triangle_ends = c(TRUE, FALSE),
+                    toptitle = 'Downsacaled RF 100', units = 'precipitation (mm)')
+        dev.off()
+    }
+    result$lon <- fs$lon
+    result$lat <- fs$lat
+    result <- CST_MergeDims(result, merge_dims = c("ftime", "monthly"),
+                         na.rm = TRUE)
+    result$Dataset <- paste0('RF100_ECMWFC3S_QM_member_', member, '_real_', 
+                            realizations)
+    result$Dates[[1]] <- exp$Dates[[1]]
+    CST_SaveExp(result, destination = dir_output, 
+                extra_string = paste0('member', k))
+    gc() 
+    k = k + 1
+    rm(list= list('fs', 'result', 'data'))
+  }
+} 
+
+
diff --git a/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R
new file mode 100644
index 0000000000000000000000000000000000000000..ee2df1943eef5547d21401f55c8fe10d2853132d
--- /dev/null
+++ b/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R
@@ -0,0 +1,215 @@
+#!/usr/bin/env Rscript
+rm(list=ls()); gc();
+
+### Creation date: 3rd June 2021
+# Author: N. Pérez-Zanón
+# Adapted from the original version Terzago et al. 2020
+# https://drive.google.com/file/d/1qp2gbtKdBl4XmsyOeaEhFENwpeUuJwkf/view
+# Refer CSTools package and manuscript when using it.
+# ----------------------------------------
+# This code is divided in 4 steps plus visualization of the results
+# STEP 1: Compute Slope values for RainFARM downscaling method
+# STEP 2: Apply Quantile Mapping Correction to Seasonal Precipitation Forecast
+# STEP 3: Compute Weights for RainFARM downscaling method
+# STEP 4: Apply RainFARM downscaling method
+# ----------------------------------------
+# Note: STEP 3 requires a machine with internet connexion. 
+#       In this file, the lines are commented since they have been run and the 
+#       result saved on disk, then the result is loaded.
+# ----------------------------------------
+#
+# Load required libraries and setup output directory:
+library(CSTools)
+library(ClimProjDiags)
+library(zeallot)
+library(ragg)
+dir_output <- '/esarchive/scratch/nperez/CSTools_manuscript/v20210603/' #slash end
+# --------------------------------------------
+# STEP 1:
+# --------------------------------------------
+era5 <- list(name = 'era5',
+             path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+years <- unlist(lapply(1993:2018, function(x){
+                       paste0(x, sprintf("%02d",1:12), '01')}))
+era5 <- CST_Load(var = 'prlr',
+                 exp = list(era5),
+                 sdates = years, nmember = 1,
+                 storefreq = "daily", sampleperiod = 1,
+                 latmin = 37.5, latmax = 53.25, lonmin = 2.5, lonmax = 18.25,
+                 output = 'lonlat', nprocs = 1)
+era5$data <- era5$data * 24 * 3600 * 1000 # with or without this line -> same result
+era5 <- CST_SplitDim(era5, split_dim = 'sdate', indices = rep(1:12, 26))
+slope <- CST_RFSlope(era5, time_dim = c('sdate', 'ftime'), kmin = 5)
+save(slope, file = paste0(dir_output, 'Slope.RDS'), version = 2)
+slope_plot <- slope
+
+# --------------------------------------------
+# STEP 2:
+# --------------------------------------------
+StartDates <- paste0(1993:2018, '1101')
+exp <- list(name = 'ecmwfS5',
+            path = "/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc")
+obs <- list(name = 'era5',
+            path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+#obs <- list(name = 'era5', path = '/esarchive/scratch/nperez/ERA5/daily_total_prlr_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
+
+c(exp, obs) %<-% CST_Load(var = 'prlr', exp = list(exp), obs = list(obs),
+                          sdates = StartDates, nmember = 25,
+                          storefreq = "daily", sampleperiod = 1,
+                          latmin = 42, latmax = 49, lonmin = 4, lonmax = 11,
+                          output = 'lonlat', nprocs = 1)
+exp <- CST_SplitDim(exp, split_dim = c('ftime'))
+obs <- CST_SplitDim(obs, split_dim = c('ftime'))
+exp$data <- exp$data * 24 * 3600 * 1000
+obs$data <- obs$data * 24 * 3600 * 1000
+exp$data[which(exp$data < 0)] <- 0
+exp.qm <- CST_QuantileMapping(exp, obs, method = "QUANT",
+                              wet.day = FALSE,
+                              sample_dims = c('member', 'sdate', 'ftime'),
+                              ncores = 4)
+save(exp.qm, file = paste0(dir_output, 'ExpQM.RDS'), version = 2)
+save(exp, file = paste0(dir_output, 'Exp.RDS'), version = 2)
+# --------------------------------------------
+# STEP 3:
+# --------------------------------------------
+#library(raster); library(s2dv); library(CSTools)
+#worldclim <- getData("worldclim", var = "prec", res = 0.5, lon = 5, lat = 45)
+#wc_month <- lapply(1:12, FUN = function(x) {
+#                   res <- crop(worldclim[[x]],
+#                               extent(3.5, 11.5, 41.5, 49.5))
+#                   res <- as.array(res)
+#                   names(dim(res)) <- c('lat', 'lon', 'month')
+#                   return(res)
+#                   })
+#xy <- xyFromCell(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5)),
+#                 1:length(crop(worldclim[[1]], extent(3.5, 11.5, 41.5, 49.5))))
+#lons <- unique(xy[,1])
+#lats <- unique(xy[,2])
+#wc_month <- unlist(wc_month)
+#dim(wc_month) <- c(lat = length(lats), lon = length(lons), monthly = 12)
+#wc_month <- Reorder(wc_month, c('lon', 'lat', 'monthly'))
+#wc_month <- s2dv_cube(data = wc_month, lon = lons, lat = lats,
+#                       Datasets = 'WorldClim')
+#weight <- CST_RFWeights(wc_month, lon = 4:11, lat = 49:42, nf = 4)
+#save(weight, file = paste0(dir_output, 'weightsRF4.RDS'))
+
+load(paste0(dir_output, 'weightsRF4.RDS'))
+
+# --------------------------------------------
+# STEP 4:
+# --------------------------------------------
+Rprof()
+weights <- Subset(weight$data, along = 'monthly', indices = c(11, 12, 1:6))
+slope <- Subset(slope, along = 'monthly', indices = c(11, 12, 1:6), 
+                drop = 'non-selected')
+fs <- CST_RainFARM(exp.qm, nf = 4,
+                   weights = weights, slope = slope,
+                   kmin = 1, nens = 10, verbose = TRUE,
+                   time_dim = c("member", "sdate", "ftime"), nprocs = 8,
+                   drop_realization = TRUE)
+newfs <- CST_MergeDims(fs, merge_dims = c("ftime", "monthly"),
+                         na.rm = TRUE)
+
+newfs$Dates[[1]] <- exp$Dates[[1]]
+CST_SaveExp(newfs, destination = paste0(dir_output, 'RF4/'))
+
+Rprof(NULL)
+profile.info <- summaryRprof(paste0(dir_output, "Rprof.out"))
+# --------------------------------------------
+# Visualization
+# --------------------------------------------
+library(s2dv)
+agg_png(paste0(dir_output, "EXP_11dec.png"),
+        width = 800, height = 900, units = 'px',res = 144)
+PlotEquiMap(exp$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat,
+            filled.continents = FALSE, bar_limits = c(0,40),
+           intylat = 2, intxlon = 2, title_scale = 0.8,
+           bar_label_scale = 1.3, axes_label_scale = 1.2,
+           triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5,
+           toptitle = 'SEAS5', units = 'precipitation (mm)')
+dev.off()
+agg_png(paste0(dir_output, "EXPQM_11dec.png"),
+        width = 800, height = 900, units = 'px',res = 144) 
+PlotEquiMap(exp.qm$data[1,1,1,11,,,2],lon = exp$lon, lat = exp$lat,
+            filled.continents = FALSE, bar_limits = c(0,40),
+           intylat = 2, intxlon = 2, title_scale = 0.8,
+           bar_label_scale = 1.3, axes_label_scale = 1.2,
+           triangle_ends = c(TRUE, FALSE), degree_sym = TRUE, units_scale = 1.5,
+           toptitle = 'Bias Adjusted', units = 'precipitation (mm)')
+dev.off()
+agg_png(paste0(dir_output, "RF4_Down_11dec.png"),
+        width = 800, height = 900, units = 'px',res = 144)
+PlotEquiMap(fs$data[1,1,1,11,,,2],lon = fs$lon, lat = fs$lat,
+            filled.continents = FALSE, bar_limits = c(0,40), 
+           intylat = 2, intxlon = 2, title_scale = 0.8,
+           bar_label_scale = 1.3, axes_label_scale = 1.2,
+           triangle_ends = c(TRUE, TRUE), degree_sym = TRUE, units_scale = 1.5,
+           toptitle = 'Downscaled nf 4', units = 'precipitation (mm)')
+dev.off()
+agg_png(paste0(dir_output, "RF4_WeightsDec.png"),
+        width = 800, height = 900, units = 'px',res = 144)
+PlotEquiMap(weight$data[,,12], lon = weight$lon, lat = weight$lat,
+            filled.continents = FALSE, title_scale = 0.8,
+           bar_label_scale = 1.3, axes_label_scale = 1.2,
+            intylat = 2, intxlon = 2, degree_sym = TRUE,
+            toptitle = 'December Weights nf 4')
+dev.off()
+agg_png(paste0(dir_output, "Slope.png"),
+        width = 700, height = 700, units = 'px',res = 144)
+plot(1:12, slope_plot, type = 'b', main = 'Slope', pch = 16, xlab = 'month',
+      ylab = 'Slope', bty = 'n', cex.main = 1.5, cex.lab = 1.3, cex = 1.5)
+lines(12, slope_plot[12], type = 'p', col = 'red', pch = 16, cex = 1.3)
+dev.off()
+
+# Plot ForecastPDF
+library(abind)
+#obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T)
+#print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T))
+#expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T)
+#exp.qmts <- MeanDims(exp.qm$data,  c('lat', 'lon'), na.rm = T)
+print("Quantiles gridpoint")
+print(quantile(as.vector(exp.qm$data[1,,,11,5,4,2]), c(0.1,0.3,0.6,0.9)))
+data <- rbind(exp$data[1, ,1,11,5,4,2], exp.qm$data[1,,1,11,5,4,2])
+empty <- array(NA, c(dataset = 2, members = 225))
+names(dim(data)) <- names(dim(empty))
+data <- abind(data, empty, along = 2)
+names(dim(data)) <- names(dim(empty))
+data <- abind(data, fs$data[1,,1,11,15,18 ,2], along = 1)
+names(dim(data)) <- names(dim(empty))
+print(dim(data))
+print(quantile(obs$data[1,,,11,5,4,2]))
+agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_GridPoint.png"),
+        width = 1000, height = 1000, units = 'px', res = 144)
+PlotForecastPDF(data, tercile.limits = c(0.02, 1.17), 
+                extreme.limits = c(0.0155555, 7.75), color.set = 'hydro', 
+                var.name = "Precipitation (mm)", add.ensmemb = 'no',
+               title = "Forecasts issued on Nov 1993 for 11th December 1993",
+               fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4"))
+dev.off()
+
+obsts <- MeanDims(obs$data, c('lat', 'lon'), na.rm = T)
+print(quantile(obsts, c(0.1, 0.3, 0.6, 0.9), na.rm = T))
+expts <- MeanDims(exp$data, c('lat', 'lon'), na.rm = T)
+exp.qmts <- MeanDims(exp.qm$data,  c('lat', 'lon'), na.rm = T)
+empty <- array(NA, c(dataset = 2, member = 225, sdate = 26, ftime = 31, monthly = 8))
+data <- abind(expts, exp.qmts, along = 1)
+names(dim(data)) <- names(dim(expts))
+data <- abind(data, empty, along = 2)
+names(dim(data)) <- names(dim(expts))
+fsts <- MeanDims(fs$data, c('lat', 'lon'), na.rm = T)
+data <- abind(data, fsts, along = 1)
+names(dim(data)) <- c('dataset', 'members', 'sdate', 'ftime', 'monthly')
+agg_png(paste0(dir_output, "/FiguresPDF_RF4_11December_Areave.png"),
+        width = 1400, height = 800, units = 'px', res = 144)
+PlotForecastPDF(data[,,1,11,2], tercile.limits = c(0.67, 2.5), 
+                extreme.limits = c(0.09, 7.3), color.set = 'hydro', 
+                var.name = "Precipitation (mm)", 
+               title = "Forecasts issued on Nov 1993 for 11th December 1993",
+               fcst.names = c("SEAS5", "Bias Adjusted", "Downscaled nf 4"))
+dev.off()
+
+
+
+
+
+
diff --git a/inst/doc/UseCase3_data_preparation_SCHEME_model.R b/inst/doc/UseCase3_data_preparation_SCHEME_model.R
index 0ff104dbd5457f6a99369237b24607479587eb01..ada24ef2c762ca78eef563a2fca2c97a75a44903 100644
--- a/inst/doc/UseCase3_data_preparation_SCHEME_model.R
+++ b/inst/doc/UseCase3_data_preparation_SCHEME_model.R
@@ -1,3 +1,6 @@
+# Author: Bert Van Schaeybroeck
+# Use Case 3: Seasonal forecasts for a river flow
+# ----------------------------------------------- 
 rm(list = ls())
 library(CSTools)
 library(s2dverification)
@@ -18,6 +21,7 @@ domain.high.res <- "greece_high_res"
 domain.low.res <- "greece_low_res"
 sdate.mon.to.use <- 5 #Month of startdate for ECMWF Sys5 possibilities are 5 (May) or 11 (November)
 sdate.day.to.use <- 1 #Day of startdate for ECMWF Sys5 only possibility is 1
+make.plot.msk <- F #mask to indicate if figures need to be made based on the output of the Analog function.
 
 #LOCAL PARAMETERS (to be adjusted for each working system)
 #---------------------------------------------------------
@@ -210,6 +214,10 @@ exp.msl.eur.merge <- CST_MergeDims(
 	merge_dims = c("sdate", "ftime"),
 	rename_dim = "sdate")
 
+obs.msl.eur.merge.an <- CST_Anomaly(exp = obs.msl.eur.merge, dim_anom = 3)
+exp.msl.eur.merge.an <- CST_Anomaly(exp = exp.msl.eur.merge, dim_anom = 3)
+
+
 #Load observational and forecast set of variable that needs to be calibrated and downscaled:
 file.to.load <- paste0(dir.rdata, "data_all.RData")
 if(file.exists(file.to.load)){
@@ -333,8 +341,8 @@ lon.low.res <- as.vector(cal.merge$lon)
 lat.low.res <- as.vector(cal.merge$lat)
 lon.high.res <- as.vector(obs.high.res$lon)
 lat.high.res <- as.vector(obs.high.res$lat)
-lon.eur <- as.vector(obs.msl.eur.merge$lon)
-lat.eur <- as.vector(obs.msl.eur.merge$lat)
+lon.eur <- as.vector(obs.msl.eur.merge.an$lon)
+lat.eur <- as.vector(obs.msl.eur.merge.an$lat)
 
 #amount of lead times in months. For ECMWF Sys5 it is 7:
 amt.lead.mon <- as.numeric(dim(cal.merge$data)["monthly"])
@@ -387,9 +395,9 @@ i.mbr.obs <- 1
 for(i.mbr in seq(1, amt.mbr)){
 	for(i.mon in seq(1, amt.lead.mon)){
 		for(i.sdate in seq(1, amt.sdate)){
-			i.mbr <- 1
-			i.mon = 1
-			i.sdate = 24
+			#i.mbr <- 1
+			#i.mon = 1
+			#i.sdate = 24
 			date.to.use <- sub.time[ i.sdate, i.mon]
 			date.an.lst <- sub.time[ , i.mon]
 			cat("i.mbr = ", i.mbr, ", i.mon =", i.mon, ", i.sdate = ", 
@@ -400,10 +408,10 @@ for(i.mbr in seq(1, amt.mbr)){
 			exp.low.res.tmp <- exp.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon]
 			cal.low.res.tmp <- cal.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon]
 			#Extract the large-scale pressure field of that day
-			exp.msl.eur.tmp <- exp.msl.eur.merge$data[i.dataset, i.mbr, i.sdate, , , i.mon]
+			exp.msl.eur.tmp <- exp.msl.eur.merge.an$data[i.dataset, i.mbr, i.sdate, , , i.mon]
 			
 			#Extract all observations that will be used to find analogs 
-			obs.msl.eur.tmp <- obs.msl.eur.merge$data[i.dataset, i.mbr.obs, , , , i.mon]#-i.sdate
+			obs.msl.eur.tmp <- obs.msl.eur.merge.an$data[i.dataset, i.mbr.obs, , , , i.mon]#-i.sdate
 			obs.low.res.tmp <- obs.low.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate
 			obs.high.res.tmp <- obs.high.res.merge$data[i.dataset, i.mbr.obs, , , , i.mon] #-i.sdate
 			names(dim(obs.high.res.tmp)) <- c("time", "lat", "lon")
@@ -428,9 +436,77 @@ for(i.mbr in seq(1, amt.mbr)){
 					excludeTime = date.to.use,
 					AnalogsInfo = T, 
 					nAnalogs = 1000)
+				
+				
+				if(make.plot.msk){
+					corr.date <- as.character(res$dates[1]) #select the date of the most
+					corr.dex <- which(date.an.lst == corr.date)
+					
+					#The following figure shows the uncalibrated raw model field (analogous to Fig. 9a)
+					file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, 
+						"_sdate_", date.to.use, "_exp.low.res.pdf")
+					pdf(file = file.fig)
+					PlotEquiMap(
+						exp.low.res.tmp[ , ], 
+						lon = obs.low.res.merge$lon, 
+						lat = obs.low.res.merge$lat,
+						filled.continents = F, 
+						intylat = 2, 
+						intxlon = 2, 
+						title_scale = 0.7, #bar_limits = c(0, 60),
+						units = "precipitation (mm)") 
+					dev.off()
+					
+					#The following figure includes the calibrated model field (analogous to Fig. 9b)
+					file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, 
+						"_sdate_", date.to.use, "_cal.low.res.pdf")
+					pdf(file = file.fig)
+					PlotEquiMap(
+						cal.low.res.tmp, 
+						lon = obs.low.res.merge$lon, 
+						lat = obs.low.res.merge$lat,
+						filled.continents = F,
+						intylat = 2, 
+						intxlon = 2, 
+						title_scale = 0.7, #bar_limits = c(0, 60),
+						units = "precipitation (mm)")
+					
+					#The following figure includes the analog upscaled field (analogous to Fig. 9c)
+					file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, 
+						"_sdate_", date.to.use, "_obs.low.res.pdf")
+					pdf(file = file.fig)
+					PlotEquiMap(
+						obs.low.res.tmp[corr.dex, , ], 
+						lon = obs.low.res.merge$lon, 
+						lat = obs.low.res.merge$lat,
+						filled.continents = F, 
+						intylat = 2, 
+						intxlon = 2, 
+						title_scale = 0.7, #bar_limits = c(0, 60),
+						units = "precipitation (mm)") 
+					dev.off()
+					
+					#The following figure includes the analog field (analogous to Fig. 9d)
+					file.fig <- paste0("mbr_", i.mbr, "_mon_", i.mon, 
+						"_sdate_", date.to.use, "_obs.high.res.pdf")
+					pdf(file = file.fig)
+					PlotEquiMap(
+						obs.high.res.tmp[corr.dex, , ], 
+						lon = obs.high.res.merge$lon, 
+						lat = obs.high.res.merge$lat,
+						filled.continents = F,
+						intylat = 2, 
+						intxlon = 2, 
+						title_scale = 0.7, #bar_limits = c(0, 60),
+						units = "precipitation (mm)")
+					dev.off()
+					
+					
+				}
 			}
 		}
 	}
 }
 
 
+
diff --git a/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh b/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh
new file mode 100644
index 0000000000000000000000000000000000000000..2e1dd93945d0794a4451584920c45e8791d46c38
--- /dev/null
+++ b/inst/doc/launch_UseCase2_PrecipitationDownscaling_RF4.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+#BSUB -W 6:00
+#BSUB -n 16 
+#BSUB -M 7000
+#BSUB -J RainFARM_Downsc
+#BSUB -o /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.out
+#BSUB -e /esarchive/scratch/nperez/CSTools_manuscript/v20210603/lsf-%J.err
+
+module load R
+Rscript /esarchive/scratch/nperez/git/cstools/inst/doc/UseCase2_PrecipitationDownscaling_RainFARM_RF4.R
diff --git a/man/CST_SaveExp.Rd b/man/CST_SaveExp.Rd
index ddd9164e1c9d88e284d4653d1b05941699e2adf4..92c972826018a806b67c84aa82644350a01130c6 100644
--- a/man/CST_SaveExp.Rd
+++ b/man/CST_SaveExp.Rd
@@ -5,7 +5,7 @@
 \title{Save CSTools objects of class 's2dv_cube' containing experiments or observed 
 data in NetCDF format}
 \usage{
-CST_SaveExp(data, destination = "./CST_Data")
+CST_SaveExp(data, destination = "./CST_Data", extra_string = NULL)
 }
 \arguments{
 \item{data}{an object of class \code{s2dv_cube}.}
@@ -15,6 +15,8 @@ to save the data. NetCDF file for each starting date are saved into the
 folder tree: destination/experiment/variable/. By default the function 
 creates and saves the data into the folder "CST_Data" in the working 
 directory.}
+
+\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization. It would be added to the file name between underscore characters.}
 }
 \description{
 This function allows to divide and save a object of class 
diff --git a/man/PlotCombinedMap.Rd b/man/PlotCombinedMap.Rd
index e631761e8188e591b045a905852a0be76c0904fa..3d6661e14c45a3f9d5906125bce7d5a187e73d5a 100644
--- a/man/PlotCombinedMap.Rd
+++ b/man/PlotCombinedMap.Rd
@@ -19,6 +19,7 @@ PlotCombinedMap(
   dots = NULL,
   bar_titles = NULL,
   legend_scale = 1,
+  cex_bar_titles = 1.5,
   fileout = NULL,
   width = 8,
   height = 5,
@@ -61,6 +62,8 @@ layers via the parameter 'dot_symbol'.}
 
 \item{legend_scale}{Scale factor for the size of the colour bar labels. Takes 1 by default.}
 
+\item{cex_bar_titles}{Scale factor for the sizes of the bar titles. Takes 1.5 by default.}
+
 \item{fileout}{File where to save the plot. If not specified (default) a graphics device will pop up. Extensions allowed: eps/ps, jpeg, png, pdf, bmp and tiff}
 
 \item{width}{File width, in the units specified in the parameter size_units (inches by default). Takes 8 by default.}
diff --git a/man/SaveExp.Rd b/man/SaveExp.Rd
index 40ace2dbab6cc9bb1c84292ce7387f6a3e991dc0..345974b22a80c52f5cce724fd4f5c560ea876834 100644
--- a/man/SaveExp.Rd
+++ b/man/SaveExp.Rd
@@ -15,7 +15,8 @@ SaveExp(
   Dates,
   cdo_grid_name,
   projection,
-  destination
+  destination,
+  extra_string = NULL
 )
 }
 \arguments{
@@ -40,9 +41,11 @@ SaveExp(
 \item{projection}{a character string indicating the projection name}
 
 \item{destination}{a character string indicating the path where to store the NetCDF files}
+
+\item{extra_string}{a character string to be include as part of the file name, for instance, to identify member or realization.}
 }
 \value{
-the function creates as many files as sdates per dataset. Each file could contain multiple members
+the function creates as many files as sdates per dataset. Each file could contain multiple members. It would be added to the file name between underscore characters.
 The path will be created with the name of the variable and each Datasets.
 }
 \description{