diff --git a/R/AbsToProbs.R b/R/AbsToProbs.R
index 247ba38d218d2a1fd9d438dd47145cc1dc5ac864..db8891ea5b9ded36537da1686af149b452872223 100644
--- a/R/AbsToProbs.R
+++ b/R/AbsToProbs.R
@@ -5,8 +5,8 @@
 #'(start dates) are provided, the function will create the Cumulative 
 #'Distribution Function excluding the corresponding initialization. 
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to define the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/AccumulationExceedingThreshold.R b/R/AccumulationExceedingThreshold.R
index 0b312e9dfc6878585f177369b0058af132a3e3dd..feee031aef0cee3387ad08d38ab6c70ae118c1d1 100644
--- a/R/AccumulationExceedingThreshold.R
+++ b/R/AccumulationExceedingThreshold.R
@@ -12,8 +12,8 @@
 #'                  temperatures and 10°C between April 1st and October 31st}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in
-#'  package CSTools. 
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
 #'  object or a multidimensional array with named dimensions. It must be in the 
 #'  same units and with the common dimensions of the same length as parameter 
diff --git a/R/PeriodAccumulation.R b/R/PeriodAccumulation.R
index 41c9b6783395cdffc3450c37294f3fd4e61fb7fd..89a4d37f5162e57a23bcb699457651eaa337ac04 100644
--- a/R/PeriodAccumulation.R
+++ b/R/PeriodAccumulation.R
@@ -10,8 +10,8 @@
 #'                   August 21st to October 21st}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/PeriodMax.R b/R/PeriodMax.R
index d0004e0831119d3ff6166a7c20a0810987e3c0fc..8c56156a73d2f063b90bb63f4e91be1ead6ac5a2 100644
--- a/R/PeriodMax.R
+++ b/R/PeriodMax.R
@@ -12,8 +12,8 @@
 #'                    that prevails during the wettest month.}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/PeriodMean.R b/R/PeriodMean.R
index 69ffc8c7cb75047d54ed074ef999b3d672f1cdd2..9106927edf7a3c446a31f7f5da3ef3374582312a 100644
--- a/R/PeriodMean.R
+++ b/R/PeriodMean.R
@@ -10,8 +10,8 @@
 #'                    maximum temperature from April 1st to May 31st}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/PeriodMin.R b/R/PeriodMin.R
index afe5eb83f71db31dc24ddf62ca5a0de340102c20..74547cc7e609d4f46d227048c2aebf57faa2d4b5 100644
--- a/R/PeriodMin.R
+++ b/R/PeriodMin.R
@@ -12,8 +12,8 @@
 #'                    that prevails during the driest month.}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/PeriodVariance.R b/R/PeriodVariance.R
index 77bf68dc70667e126b82c5009f478868b5e65439..f95b8e4e5033e39704934a1670eca8ca2c94f313 100644
--- a/R/PeriodVariance.R
+++ b/R/PeriodVariance.R
@@ -16,8 +16,8 @@
 #'                    and is expressed as a percentage}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start An optional parameter to defined the initial date of the period 
 #'  to select from the data by providing a list of two elements: the initial 
 #'  date of the period and the initial month of the period. By default it is set
diff --git a/R/QThreshold.R b/R/QThreshold.R
index 19cda3c64cca5003abbb294bdceda1a72cf0d9fd..6bd8d83eeae4e936a9790f67bca1d015adc6a9d8 100644
--- a/R/QThreshold.R
+++ b/R/QThreshold.R
@@ -22,8 +22,8 @@
 #'        the sample used must be especified in sdate_dim parameter.}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param threshold An 's2dv_cube' object as output of a 'CST_' function in the 
 #'  same units as parameter 'data' and with the common dimensions of the element
 #'  'data' of the same length. A single scalar is also possible.
diff --git a/R/SelectPeriodOnData.R b/R/SelectPeriodOnData.R
index bef70be2399f9f9256e91635056a7e63fe2864d3..a4cc07c1e7493e44ae6c944a740e9bd1eda6de3a 100644
--- a/R/SelectPeriodOnData.R
+++ b/R/SelectPeriodOnData.R
@@ -2,8 +2,8 @@
 #'
 #' Auxiliary function to subset data for a specific period.
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param start A parameter to defined the initial date of the period to select 
 #'  from the data by providing a list of two elements: the initial date of the 
 #'  period and the initial month of the period.
diff --git a/R/Threshold.R b/R/Threshold.R
index ee4fa38110acb87d5f2b7786d98a553f57c34b0d..314bc704e1e67cb7356804f89825e1ff3b9e729a 100644
--- a/R/Threshold.R
+++ b/R/Threshold.R
@@ -5,8 +5,8 @@
 #'day (time). This function calculates the corresponding value of a percentile 
 #'given a dataset.
 #'
-#'@param data An 's2dv_cube' object as provided function \code{CST_Load} in
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param threshold A single scalar or vector indicating the relative 
 #'  threshold(s). It must contain values between 0 and 1.
 #'@param start An optional parameter to defined the initial date of the period 
diff --git a/R/TotalSpellTimeExceedingThreshold.R b/R/TotalSpellTimeExceedingThreshold.R
index 1bc1bb3862aea135dcbe74d3fd08169b91aa1627..cd2676e0f1d71585de7cb71ff1867b7adfa29cbe 100644
--- a/R/TotalSpellTimeExceedingThreshold.R
+++ b/R/TotalSpellTimeExceedingThreshold.R
@@ -16,8 +16,8 @@
 #'by using function \code{AbsToProbs}. See section @examples.
 #'@seealso [Threshold()] and [AbsToProbs()].
 #'
-#'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
 #'  object or a multidimensional array with named dimensions. It must be in the 
 #'  same units and with the common dimensions of the same length as parameter 
diff --git a/R/TotalTimeExceedingThreshold.R b/R/TotalTimeExceedingThreshold.R
index 5c6bb628bb0c4b35845cb5e6d095253d694eec71..fa57208fe4309955c682bc2c0fc2ccd144981218 100644
--- a/R/TotalTimeExceedingThreshold.R
+++ b/R/TotalTimeExceedingThreshold.R
@@ -20,8 +20,8 @@
 #'                    32 between April 21st and June 21st}
 #'}
 #'
-#'@param data An 's2dv_cube' object as provided by function \code{CST_Load} in 
-#'  package CSTools.
+#'@param data An 's2dv_cube' object as provided function \code{CST_Start} or 
+#'  \code{CST_Load} in package CSTools.
 #'@param threshold If only one threshold is used, it can be an 's2dv_cube' 
 #'  object or a multidimensional array with named dimensions. It must be in the 
 #'  same units and with the common dimensions of the same length as parameter 
diff --git a/R/WindPowerDensity.R b/R/WindPowerDensity.R
index 3dc835a2208d2f998928760357da7bf2d86b5e12..4a2e51e94171d3ef81e07a9da2f80fdfe1b05981 100644
--- a/R/WindPowerDensity.R
+++ b/R/WindPowerDensity.R
@@ -7,7 +7,7 @@
 #'it will give inaccurate results if used with period means.
 #'
 #'@param wind An 's2dv_cube' object with instantaneous wind speeds expressed in 
-#'  m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge.
+#'  m/s obtained from CST_Start or s2dv_cube functions from CSTools pacakge.
 #'@param ro A scalar, or alternatively a multidimensional array with the same
 #'  dimensions as wind, with the air density expressed in kg/m^3. By default it
 #'  takes the value 1.225, the standard density of air at 15ºC and 1013.25 hPa.
diff --git a/man/CST_AbsToProbs.Rd b/man/CST_AbsToProbs.Rd
index ef8f42dc9a36098904cbf51c568762722cb20e27..45a733d6b47efbe9ba2868fe64b38f2c2961ac5d 100644
--- a/man/CST_AbsToProbs.Rd
+++ b/man/CST_AbsToProbs.Rd
@@ -15,8 +15,8 @@ CST_AbsToProbs(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to define the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_AccumulationExceedingThreshold.Rd b/man/CST_AccumulationExceedingThreshold.Rd
index 2bdee762b40ce882de9fd3eacd78ec47f95812b8..0bea015823d0b08ecfa65db23385c6bf8c190aec 100644
--- a/man/CST_AccumulationExceedingThreshold.Rd
+++ b/man/CST_AccumulationExceedingThreshold.Rd
@@ -17,8 +17,8 @@ CST_AccumulationExceedingThreshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
 object or a multidimensional array with named dimensions. It must be in the 
diff --git a/man/CST_PeriodAccumulation.Rd b/man/CST_PeriodAccumulation.Rd
index 2c8fb5a16ceb168837496da30360329aab142306..0820a0dc57848cd6cf257636cec0a12d3568a7da 100644
--- a/man/CST_PeriodAccumulation.Rd
+++ b/man/CST_PeriodAccumulation.Rd
@@ -14,8 +14,8 @@ CST_PeriodAccumulation(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_PeriodMax.Rd b/man/CST_PeriodMax.Rd
index 7b01d145482423cf55a375ba9bc9a31fb7f6d4eb..a372a1d8bb52ee97153157f75398b771ec0a2e9e 100644
--- a/man/CST_PeriodMax.Rd
+++ b/man/CST_PeriodMax.Rd
@@ -14,8 +14,8 @@ CST_PeriodMax(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_PeriodMean.Rd b/man/CST_PeriodMean.Rd
index ab420665dabdb4fa8f94a40f8c0ed2d65552bda0..323adae91e962df2f10fe09496d7ac76b10d5bb8 100644
--- a/man/CST_PeriodMean.Rd
+++ b/man/CST_PeriodMean.Rd
@@ -14,8 +14,8 @@ CST_PeriodMean(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_PeriodMin.Rd b/man/CST_PeriodMin.Rd
index 939569986386893c495d8789bdbaca904393f1d4..9d3833f4937576a39e8b0b97dd6205511b538e4c 100644
--- a/man/CST_PeriodMin.Rd
+++ b/man/CST_PeriodMin.Rd
@@ -14,8 +14,8 @@ CST_PeriodMin(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_PeriodVariance.Rd b/man/CST_PeriodVariance.Rd
index e28bf9551727659ffb53e257ec2c29fe7a8829b2..b04b4ed03fd1298e16163fe185d2326c76e52e08 100644
--- a/man/CST_PeriodVariance.Rd
+++ b/man/CST_PeriodVariance.Rd
@@ -14,8 +14,8 @@ CST_PeriodVariance(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{An optional parameter to defined the initial date of the period 
 to select from the data by providing a list of two elements: the initial 
diff --git a/man/CST_QThreshold.Rd b/man/CST_QThreshold.Rd
index b1683756da67d4d64403e876226b2d941af59d25..a07d130ef89a2fe4e3a8491072dddef497a91027 100644
--- a/man/CST_QThreshold.Rd
+++ b/man/CST_QThreshold.Rd
@@ -16,8 +16,8 @@ CST_QThreshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{threshold}{An 's2dv_cube' object as output of a 'CST_' function in the 
 same units as parameter 'data' and with the common dimensions of the element
diff --git a/man/CST_SelectPeriodOnData.Rd b/man/CST_SelectPeriodOnData.Rd
index 5f12633efb1d7c2d364ff4a8c24a6b644936cd95..5e4eff434bec51f97c5334e0c5120b05330a15f0 100644
--- a/man/CST_SelectPeriodOnData.Rd
+++ b/man/CST_SelectPeriodOnData.Rd
@@ -7,8 +7,8 @@
 CST_SelectPeriodOnData(data, start, end, time_dim = "time", ncores = NULL)
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{start}{A parameter to defined the initial date of the period to select 
 from the data by providing a list of two elements: the initial date of the 
diff --git a/man/CST_Threshold.Rd b/man/CST_Threshold.Rd
index 2b4ea2faec47b364e745eba504b28ea74f7a7fd4..07571b54ca01270fe3c57a4047fb711c6d70a29f 100644
--- a/man/CST_Threshold.Rd
+++ b/man/CST_Threshold.Rd
@@ -17,8 +17,8 @@ CST_Threshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided function \code{CST_Load} in
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{threshold}{A single scalar or vector indicating the relative 
 threshold(s). It must contain values between 0 and 1.}
diff --git a/man/CST_TotalSpellTimeExceedingThreshold.Rd b/man/CST_TotalSpellTimeExceedingThreshold.Rd
index 940478f685741c1d16343bad9c7330981110a9e5..831bf2cf9f3af2f549fff3f64734bd065ac6b52a 100644
--- a/man/CST_TotalSpellTimeExceedingThreshold.Rd
+++ b/man/CST_TotalSpellTimeExceedingThreshold.Rd
@@ -16,8 +16,8 @@ CST_TotalSpellTimeExceedingThreshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
 object or a multidimensional array with named dimensions. It must be in the 
diff --git a/man/CST_TotalTimeExceedingThreshold.Rd b/man/CST_TotalTimeExceedingThreshold.Rd
index cd9916398ed31952339db3f0072b86fcf8506adb..db9a3d6cabff166e688ce90ab2db93ab9ebb4834 100644
--- a/man/CST_TotalTimeExceedingThreshold.Rd
+++ b/man/CST_TotalTimeExceedingThreshold.Rd
@@ -16,8 +16,8 @@ CST_TotalTimeExceedingThreshold(
 )
 }
 \arguments{
-\item{data}{An 's2dv_cube' object as provided by function \code{CST_Load} in 
-package CSTools.}
+\item{data}{An 's2dv_cube' object as provided function \code{CST_Start} or 
+\code{CST_Load} in package CSTools.}
 
 \item{threshold}{If only one threshold is used, it can be an 's2dv_cube' 
 object or a multidimensional array with named dimensions. It must be in the 
diff --git a/man/CST_WindPowerDensity.Rd b/man/CST_WindPowerDensity.Rd
index d37fd9dea74355f771e1f4edd8245ba65774dce5..ec82d8a0ef0e5d8b696960fc6c632aa46d484e92 100644
--- a/man/CST_WindPowerDensity.Rd
+++ b/man/CST_WindPowerDensity.Rd
@@ -15,7 +15,7 @@ CST_WindPowerDensity(
 }
 \arguments{
 \item{wind}{An 's2dv_cube' object with instantaneous wind speeds expressed in 
-m/s obtained from CST_Load or s2dv_cube functions from CSTools pacakge.}
+m/s obtained from CST_Start or s2dv_cube functions from CSTools pacakge.}
 
 \item{ro}{A scalar, or alternatively a multidimensional array with the same
 dimensions as wind, with the air density expressed in kg/m^3. By default it
diff --git a/vignettes/AgriculturalIndicators.Rmd b/vignettes/AgriculturalIndicators.Rmd
index 3c9cf7d23a9b4fd68b8f819638c85f328bd7c761..6345d6c9cf8fc8f8c5693c43a7c4414623aa5f94 100644
--- a/vignettes/AgriculturalIndicators.Rmd
+++ b/vignettes/AgriculturalIndicators.Rmd
@@ -2,6 +2,8 @@
 title: "Agricultural Indicators"
 author: "Earth Sciences department, Barcelona Supercomputing Center (BSC)"
 date: "`r Sys.Date()`"
+revisor: "Eva Rifà"
+revision date: "October 2023"
 output: rmarkdown::html_vignette
 vignette: >
   %\VignetteEngine{knitr::knitr}
@@ -58,47 +60,88 @@ library(s2dv)
 
 To obtain the precipitation forecast and observation, we load the daily precipitation (**prlr** given in `var`) data sets of ECMWF SEAS5 seasonal forecast and ERA5 reanalysis for the four starting dates 20130401-20160401 (provided in `sdates`) with the entire 7-month forecast time, April-October (214 days in total given in parameter `leadtimemax`).
 
-The pathways of SEAS5 and ERA5 are given in the lists with some **whitecards (inside two dollar signs)** used to replace the variable name and iterative items such as year and month. See details of requirements in Section 4 in vignette [Data retrieval and storage](https://cran.r-project.org/package=CSTools/vignettes/Data_Considerations.html) from CSTools package.
+The pathways of SEAS5 and ERA5 are given in the lists with some **whitecards (inside two dollar signs)** used to replace the variable name and iterative items such as year and month. See details of requirements in Section 5 in vignette [Data retrieval and storage](https://cran.r-project.org/package=CSTools/vignettes/Data_Considerations.html) from CSTools package.
 
 The spatial domain covers part of Douro Valley of Northern Portugal lon=[352.25, 353], lat=[41, 41.75]. These four values are provided in `lonmin`, `lonmax`, `latmin` and `latmax`.
 
 With `grid` set to **r1440x721**, the SEAS5 forecast would be interpolated to the 0.25-degree ERA5 grid by using the **bicubic** method given in `method`.
 
 
-```
-S5path_prlr <- list(path = '/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_s0-24h/$VAR_NAME$_$START_DATE$.nc')
-
-path_ERA5prlr_CDS <- list(path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
-
+```r
 sdates <- paste0(2013:2016, '04', '01')
 
-c(prlr_exp, prlr_obs) %<-% CST_Load(var = 'prlr', 
-                                    exp = list(S5path_prlr), 
-                                    obs = list(path_ERA5prlr_CDS), 
-                                    sdates = sdates, 
-                                    lonmax = 353, lonmin = 352.25,
-                                    latmax = 41.75, latmin = 41, 
-                                    storefreq = 'daily', 
-                                    leadtimemin = 1, leadtimemax = 214, 
-                                    nmember = 3, output = "lonlat", 
-                                    grid = "r1440x721", method = 'bicubic')
+lat_min = 41
+lat_max = 41.75
+lon_min = 352.25
+lon_max = 353
+
+S5path_prlr <- paste0("/esarchive/exp/ecmwf/system5c3s/daily_mean/$var$_s0-24h/$var$_$sdate$.nc")
+prlr_exp <- CST_Start(dataset = S5path_prlr,
+                      var = "prlr",
+                      member = startR::indices(1:3),
+                      sdate = sdates,
+                      ftime = startR::indices(1:214),
+                      lat = startR::values(list(lat_min, lat_max)),
+                      lat_reorder = startR::Sort(decreasing = TRUE), 
+                      lon = startR::values(list(lon_min, lon_max)),
+                      lon_reorder = startR::CircularSort(0, 360),
+                      synonims = list(lon = c('lon', 'longitude'),
+                                      lat = c('lat', 'latitude'),
+                                      member = c('member', 'ensemble'),
+                                      ftime = c('ftime', 'time')),
+                      transform = startR::CDORemapper,
+                      transform_extra_cells = 2,
+                      transform_params = list(grid = "r1440x721",
+                                              method = "bicubic"),
+                      transform_vars = c('lat', 'lon'),
+                      return_vars = list(lat = NULL, 
+                                         lon = NULL, ftime = 'sdate'),
+                                         retrieve = TRUE)
+
+dates_exp <- prlr_exp$attrs$Dates
+
+path_ERA5prlr_CDS <- paste0("/esarchive/recon/ecmwf/era5/daily_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc")
+prlr_obs <- CST_Start(dataset = path_ERA5prlr_CDS,
+                      var = "prlr",
+                      date = unique(format(dates_exp, '%Y%m')),
+                      ftime = startR::values(dates_exp), 
+                      ftime_across = 'date',
+                      ftime_var = 'ftime',
+                      merge_across_dims = TRUE,
+                      split_multiselected_dims = TRUE,
+                      lat = startR::values(list(lat_min, lat_max)),
+                      lat_reorder = startR::Sort(decreasing = TRUE),
+                      lon = startR::values(list(lon_min, lon_max)),
+                      lon_reorder = startR::CircularSort(0, 360),
+                      synonims = list(lon = c('lon', 'longitude'),
+                                      lat = c('lat', 'latitude'),
+                                      ftime = c('ftime', 'time')),
+                      transform = startR::CDORemapper,
+                      transform_extra_cells = 2,
+                      transform_params = list(grid = "r1440x721",
+                                              method = "bicubic"),
+                      transform_vars = c('lat', 'lon'),
+                      return_vars = list(lon = NULL,
+                                         lat = NULL,
+                                         ftime = 'date'),
+                      retrieve = TRUE)
 ```
 
 The output contains data and metadata for the experiment and the observations. The elements `prlr_exp$data` and `prlr_obs$data` have dimensions:
 
 
-```
+```r
 dim(prlr_exp$data)
-#dataset  member   sdate   ftime     lat     lon
-#      1       3       4     214       4       4
+# dataset     var  member   sdate   ftime     lat     lon 
+#       1       1       3       4     214       4       4 
 dim(prlr_obs$data)
-#dataset  member   sdate   ftime     lat     lon
-#      1       1       4     214       4       4
+# dataset     var   sdate   ftime     lat     lon 
+#       1       1       4     214       4       4 
 ```
 
 To compute **SprR** of forecast and observation, we can run:
 
-```
+```r
 SprR_exp <- CST_PeriodAccumulation(prlr_exp, start = list(21, 4), end = list(21, 6))
 SprR_obs <- CST_PeriodAccumulation(prlr_obs, start = list(21, 4), end = list(21, 6))
 ```
@@ -109,41 +152,41 @@ As mentioned, these parameters are optional, the function will take the entire t
 
 The dimensions of SprR forecasts and observations are:
 
-```
+```r
 dim(SprR_exp$data)
-#dataset  member   sdate     lat     lon
-#      1       3       4       4       4
+# dataset     var  member   sdate     lat     lon 
+#       1       1       3       4       4       4 
 dim(SprR_obs$data)
-#dataset  member   sdate     lat     lon 
-#      1       1       4       4       4 
+# dataset     var   sdate     lat     lon 
+#       1       1       4       4       4  
 ```
 
 The forecast SprR for the 1st member from 2013-2016 of the 1st grid point in mm are:
 
-```
-SprR_exp$data[1, 1, , 1, 1] * 86400 * 1000
-#[1]  93.23205 230.41904 194.01412 226.52614
+```r
+SprR_exp$data[1, 1, 1, , 1, 1] * 86400 * 1000
+# [1]  93.23236 230.41754 194.01401 226.52564
 ```
 
 Dry springs will delay vegetative growth and reduce vigour and leaf area total surface. Fungal disease pressure will be lower and therefore there will be less need for protective and / or curative treatments, translating as less costs. Wet springs will promote higher vigour, increase the risk of fungal disease and disrupt vineyard operations as it may prevent machinery from getting in the vineyard due to mud. They are usually associated with higher costs.
 
 On the other hand, another moisture-related indicators, **HarvestR**, can be computed by using `PeriodAccumulation` as well, with the defined period as the following lines.
 
-```
+```r
 HarvestR_exp <- CST_PeriodAccumulation(prlr_exp, start = list(21, 8), end = list(21, 10))
 HarvestR_obs <- CST_PeriodAccumulation(prlr_obs, start = list(21, 8), end = list(21, 10))
 ```
 
 The forecast HarvestR for the 1st member from 2013-2016 of the 1st grid point in mm are:
 
-```
-HarvestR_exp$data[1, 1, , 1, 1] * 86400 * 1000
-#[1]  52.30026  42.88068 156.87961  32.18579
+```r
+HarvestR_exp$data[1, 1, 1, , 1, 1] * 86400 * 1000
+# [1]  52.30058  42.88070 156.87922  32.18567
 ```
 
 To compute the 2013-2016 ensemble-mean bias of forecast HarvestR, run
 
-```
+```r
 fcst <- drop(HarvestR_exp$data) * 86400 * 1000
 obs <- drop(HarvestR_obs$data) * 86400 * 1000
 
@@ -152,7 +195,7 @@ Bias <- MeanDims((fcst - InsertDim(obs, 1, dim(fcst)['member'])), 'member')
 
 To plot the map of ensemble-mean bias of HarvestR forecast, run
 
-```
+```r
 cols <- c('#b2182b', '#d6604d', '#f4a582', '#fddbc7', '#d1e5f0', 
           '#92c5de', '#4393c3', '#2166ac')
 
@@ -178,42 +221,81 @@ For the function `PeriodMean`, we use Growing Season Temperature (**GST**) as an
 
 Firstly, we prepare a sample data of daily mean temperature of SEAS5 and ERA5 data sets with the same starting dates, spatial domain, interpolation grid and method by running
 
-```
-S5path <- list(path = '/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$_mean/$VAR_NAME$_f6h/$VAR_NAME$_$YEAR$$MONTH$01.nc')
-ERA5path <- list(path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$_mean/$VAR_NAME$_f1h-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
-
-c(tas_exp, tas_obs) %<-% CST_Load(var = 'tas', exp = list(S5path), obs = list(ERA5path),
-                                  sdates = sdates, lonmax = 353, lonmin = 352.25,
-                                  latmax = 41.75, latmin = 41,
-                                  storefreq = 'daily',
-                                  leadtimemin = 1, leadtimemax = 214,
-                                  nmember = 3, output = "lonlat",
-                                  grid = "r1440x721", method = 'bicubic')
-```
-
-The output contains observations `tas_dv$obs$data` and forecast `tas_dv$exp$data`, and their dimensions and summaries are like
-
-```
+```r
+S5path <- paste0("/esarchive/exp/ecmwf/system5c3s/daily_mean/$var$_f6h/$var$_$sdate$.nc")
+tas_exp <- CST_Start(dataset = S5path,
+                     var = "tas",
+                     member = startR::indices(1:3),
+                     sdate = sdates,
+                     ftime = startR::indices(1:214),
+                     lat = startR::values(list(lat_min, lat_max)),
+                     lat_reorder = startR::Sort(decreasing = TRUE), 
+                     lon = startR::values(list(lon_min, lon_max)),
+                     lon_reorder = startR::CircularSort(0, 360),
+                     synonims = list(lon = c('lon', 'longitude'),
+                                     lat = c('lat', 'latitude'),
+                                     member = c('member', 'ensemble'),
+                                     ftime = c('ftime', 'time')),
+                     transform = startR::CDORemapper,
+                     transform_extra_cells = 2,
+                     transform_params = list(grid = "r1440x721",
+                                             method = "bicubic"),
+                     transform_vars = c('lat', 'lon'),
+                     return_vars = list(lat = NULL, 
+                                        lon = NULL, ftime = 'sdate'),
+                                        retrieve = TRUE)
+dates_exp <- tas_exp$attrs$Dates
+
+ERA5path <- paste0("/esarchive/recon/ecmwf/era5/daily_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc")
+tas_obs <- CST_Start(dataset = ERA5path,
+                     var = "tas",
+                     date = unique(format(dates_exp, '%Y%m')),
+                     ftime = startR::values(dates_exp), 
+                     ftime_across = 'date',
+                     ftime_var = 'ftime',
+                     merge_across_dims = TRUE,
+                     split_multiselected_dims = TRUE,
+                     lat = startR::values(list(lat_min, lat_max)),
+                     lat_reorder = startR::Sort(decreasing = TRUE),
+                     lon = startR::values(list(lon_min, lon_max)),
+                     lon_reorder = startR::CircularSort(0, 360),
+                     synonims = list(lon = c('lon', 'longitude'),
+                                     lat = c('lat', 'latitude'),
+                                     ftime = c('ftime', 'time')),
+                     transform = startR::CDORemapper,
+                     transform_extra_cells = 2,
+                     transform_params = list(grid = "r1440x721",
+                                             method = "bicubic"),
+                     transform_vars = c('lat', 'lon'),
+                     return_vars = list(lon = NULL,
+                                        lat = NULL,
+                                        ftime = 'date'),
+                     retrieve = TRUE)
+```
+
+The output contains observations `tas_obs$data` and forecast `tas_exp$data`, and their dimensions and summaries are like
+
+```r
 dim(tas_obs$data)
-#dataset  member   sdate   ftime     lat     lon 
-#      1       1       4     214       4       4 
+# dataset     var   sdate   ftime     lat     lon 
+#       1       1       4     214       4       4 
 
 dim(tas_exp$data)
-#dataset  member   sdate   ftime     lat     lon 
-#      1       3       4     214       4       4 
+# dataset     var  member   sdate   ftime     lat     lon 
+#       1       1       3       4     214       4       4 
 
 summary(tas_obs$data - 273.15)
 #  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-#  3.63   13.97   17.25   17.29   20.75   30.21
+# 3.627  13.974  17.248  17.294  20.752  30.206 
 
 summary(tas_exp$data - 273.15)
-# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-# 0.54   11.65   16.56   16.50   21.25   31.41 
+#  Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+# 0.5363 11.6517 16.5610 16.4961 21.2531 31.4063
 
 ```
 To compute the GST for both observation and forecast, run the following lines
 
-```
+```r
 # change the unit of temperature from °C to K
 
 tas_exp$data <- tas_exp$data - 273.15
@@ -230,7 +312,7 @@ Since the period considered for GST is the entire period for starting month of A
  
 The summaries and dimensions of the output are as follows:
 
-```
+```r
 summary(GST_exp$data)
 # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 # 14.23   15.78   16.50   16.50   17.17   18.70 
@@ -240,17 +322,17 @@ summary(GST_obs$data)
 # 15.34   16.77   17.22   17.29   18.00   18.75
 
 dim(GST_exp$data)
-#dataset  member   sdate     lat     lon 
-#      1       3       4       4       4 
+# dataset     var  member   sdate     lat     lon 
+#       1       1       3       4       4       4 
 
 dim(GST_obs$data)
-#dataset  member   sdate     lat     lon 
-#      1       1       4       4       4 
+# dataset     var   sdate     lat     lon 
+#       1       1       4       4       4 
 ```
 
 Here, we plot the 2013-2016 mean climatology of ERA5 GST by running 
 
-```
+```r
 # compute ERA5 GST climatology
 GST_Clim <- MeanDims(drop(GST_obs$data), 'sdate')
 
@@ -288,29 +370,68 @@ Here, we take SU35 as example, therefore the daily temperature maximum of the en
 
 Load SEAS5 and ERA5 daily temperature maximum by running
 
-```
-S5path <- list(path = '/esarchive/exp/ecmwf/system5c3s/$STORE_FREQ$/$VAR_NAME$/$VAR_NAME$_$YEAR$$MONTH$01.nc')
-ERA5path <- list(path = '/esarchive/recon/ecmwf/era5/$STORE_FREQ$/$VAR_NAME$-r1440x721cds/$VAR_NAME$_$YEAR$$MONTH$.nc')
-
-c(tasmax_exp, tasmax_obs) %<-% CST_Load(var = 'tasmax', exp = list(S5path),
-                                        obs = list(ERA5path), sdates = sdates,
-                                        lonmax = 353, lonmin = 352.25,
-                                        latmax = 41.75, latmin = 41, storefreq = 'daily',
-                                        leadtimemin = 1, leadtimemax = 214, nmember = 3,
-                                        output = "lonlat", grid = "r1440x721", 
-                                        method = 'bicubic', nprocs = 1)
+```r
+S5path <- paste0("/esarchive/exp/ecmwf/system5c3s/daily/$var$/$var$_$sdate$.nc")
+tasmax_exp <- CST_Start(dataset = S5path,
+                        var = "tasmax",
+                        member = startR::indices(1:3),
+                        sdate = sdates,
+                        ftime = startR::indices(1:214),
+                        lat = startR::values(list(lat_min, lat_max)),
+                        lat_reorder = startR::Sort(decreasing = TRUE), 
+                        lon = startR::values(list(lon_min, lon_max)),
+                        lon_reorder = startR::CircularSort(0, 360),
+                        synonims = list(lon = c('lon', 'longitude'),
+                                          lat = c('lat', 'latitude'),
+                                          member = c('member', 'ensemble'),
+                                          ftime = c('ftime', 'time')),
+                        transform = startR::CDORemapper,
+                        transform_extra_cells = 2,
+                        transform_params = list(grid = "r1440x721",
+                                                method = "bicubic"),
+                        transform_vars = c('lat', 'lon'),
+                        return_vars = list(lat = NULL, 
+                                          lon = NULL, ftime = 'sdate'),
+                                          retrieve = TRUE)
+dates_exp <- tasmax_exp$attrs$Dates
+
+ERA5path <- paste0("/esarchive/recon/ecmwf/era5/daily/$var$-r1440x721cds/$var$_$date$.nc")
+tasmax_obs <- CST_Start(dataset = ERA5path,
+                        var = "tasmax",
+                        date = unique(format(dates_exp, '%Y%m')),
+                        ftime = startR::values(dates_exp), 
+                        ftime_across = 'date',
+                        ftime_var = 'ftime',
+                        merge_across_dims = TRUE,
+                        split_multiselected_dims = TRUE,
+                        lat = startR::values(list(lat_min, lat_max)),
+                        lat_reorder = startR::Sort(decreasing = TRUE),
+                        lon = startR::values(list(lon_min, lon_max)),
+                        lon_reorder = startR::CircularSort(0, 360),
+                        synonims = list(lon = c('lon', 'longitude'),
+                                          lat = c('lat', 'latitude'),
+                                          ftime = c('ftime', 'time')),
+                        transform = startR::CDORemapper,
+                        transform_extra_cells = 2,
+                        transform_params = list(grid = "r1440x721",
+                                                method = "bicubic"),
+                        transform_vars = c('lat', 'lon'),
+                        return_vars = list(lon = NULL,
+                                          lat = NULL,
+                                          ftime = 'date'),
+                        retrieve = TRUE)
 ```
 
 Check the unit of temperature to from °C to K for the comparison with the threshold defined (for example 35°C here). 
 
-```
+```r
 tasmax_exp$data <- tasmax_exp$data - 273.15
 tasmax_obs$data <- tasmax_obs$data - 273.15
 ```
 
 Computing SU35 for forecast and observation by running
 
-```
+```r
 threshold <- 35
 SU35_exp <- CST_TotalTimeExceedingThreshold(tasmax_exp, threshold = threshold,
                                             start = list(1, 4), end = list(31, 10))
@@ -320,40 +441,40 @@ SU35_obs <- CST_TotalTimeExceedingThreshold(tasmax_obs, threshold = threshold,
 
 The summaries of SU35 forecasts and observations are given below.
 
-```
+```r
 summary(SU35_exp$data)
 # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-# 0.00    2.00    5.00    7.12   12.00   26.00 
+# 0.000   2.000   5.000   7.135  12.000  26.000 
 
 summary(SU35_obs$data)
 # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-# 0.000   0.000   1.000   2.609   5.000  10.000 
+# 0.000   0.000   1.000   2.609   5.000  10.000
 
 ```
 As shown in the summaries, SEAS5 SU35 forecasts are overestimated by 5 days in terms of mean value.
 
 Therefore, `CST_BiasCorrection` is used to bias adjust the SU35 forecasts.
 
-```
+```r
 res <- CST_BiasCorrection(obs = SU35_obs, exp = SU35_exp)
 SU35_exp_BC <- drop(res$data)
 summary(SU35_exp_BC)
 #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-# -1.419   0.000   1.613   2.831   4.756  17.768 
+# -1.523   0.000   1.613   2.830   4.756  17.768 
 ```
 
 Since there are negative values after bias adjustment, all negative data is converted to zero.
 
-```
+```r
 SU35_exp_BC[SU35_exp_BC < 0] <- 0
 summary(SU35_exp_BC)
 #   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-#  0.000   0.000   1.613   2.943   4.756  17.768 
+#  0.000   0.000   1.613   2.941   4.756  17.768 
 ```
 
 Plot the bias-adjusted SU35 forecast in 2016 by running
 
-```
+```r
 SU35_obs_Y2016 <- drop(SU35_obs$data)[4, , ]
 SU35_exp_Y2016 <- MeanDims(drop(SU35_exp$data)[, 4, , ], 'member')
 SU35_exp_BC_Y2016 <- MeanDims(SU35_exp_BC[, 4, , ], 'member')
@@ -414,34 +535,34 @@ The revised definition of SU35 is to reduce the potential influence induced by t
 
 As mentioned, the forecast is translated to its percentile by using the function `ABsToProbs` by running
 
-```
+```r
 exp_percentile <- AbsToProbs(tasmax_exp$data)
 S5txP <- aperm(drop(exp_percentile), c(2, 1, 3, 4, 5))
 ```
 
 After that, based on 35 of threshold, the percentile corresponding to each observational value can be calculated as follows.
 
-```
+```r
 obs_percentile <- QThreshold(tasmax_obs$data, threshold = 35)
 obs_percentile <- drop(obs_percentile)
 ```
 
 After translating both forecasts and observations into probabilities, the comparison can then be done by running
 
-```
+```r
 SU35_exp_Percentile <- TotalTimeExceedingThreshold(S5txP, threshold = obs_percentile, 
                                                    time_dim = 'ftime')
 ```
 
 Compute the same ensemble-mean SU35 **with percentile adjustment** in 2016 by running
 
-```
+```r
 SU35_exp_per_Y2016 <- MeanDims(SU35_exp_Percentile[4, , , ], 'member')
 ```
 
 Plot the same map for comparison
 
-```
+```r
 toptitle <- 'SU35 forecast with percentile adjustment in 2016'
 PlotEquiMap(SU35_exp_per_Y2016, 
             lon = tasmax_obs$coords$lon, lat = tasmax_obs$coords$lat,
@@ -468,14 +589,14 @@ The definition of GDD is the summation of daily differences between daily averag
 
 *Note: The data is in degrees Celsiusi at this point*
 
-```
+```r
 GDD_exp <- CST_AccumulationExceedingThreshold(tas_exp, threshold = 10, diff = TRUE)
 GDD_obs <- CST_AccumulationExceedingThreshold(tas_obs, threshold = 10, diff = TRUE)
 ```
 
 The summaries of GDD are 
 
-```
+```r
 summary(GDD_exp$data)
 # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 # 1021    1331    1480    1469    1596    1873 
@@ -487,7 +608,7 @@ summary(GDD_obs$data)
 
 To compute the correlation coefficient for the period from 2013-2016, run the following lines
 
-```
+```r
 # reorder the dimension
 fcst <- Reorder(drop(GDD_exp$data), c(4, 3, 2, 1))
 obs <- Reorder(drop(GDD_obs$data), c(3, 2, 1))
@@ -498,7 +619,7 @@ GDD_Corr <- Reorder(EnsCorr, c(2, 1))
 
 To plot the map of correlation coefficient of GDD for the 2013-2016 period.
 
-```
+```r
 cols <- c("#f7fcf5", "#e5f5e0", "#c7e9c0", "#a1d99b", "#74c476")
 toptitle <- '2013-2016 correlation coefficient of GDD'
 PlotEquiMap(GDD_Corr, lon = tas_obs$coords$lon, lat = tas_obs$coords$lat, 
@@ -522,32 +643,32 @@ One of the critical agricultural indicators related to dry spell is the **Warm S
 
 The maximum temperature data used in Section 3. Since the daily maximum temperature needs to compare to its 90th percentile, the function `Threshold` in the `CSIndicators` package is required to compute the percentile of observations used for each day. Here the same period (2013-2016) is considered. 
 
-```
-tx_p <- CST_Threshold(tasmax_obs, threshold = 0.9)
+```r
+tx_p <- CST_Threshold(tasmax_obs, threshold = 0.9, memb_dim = NULL)
 ```
 
 The output will be the 90th percentile of each day of each grid point derived by using all the years in the data.See the dimension and summary as below.
 
-```
+```r
 dim(tx_p$data)
-#dataset   ftime     lat     lon 
-#      1     214       4       4 
+# dataset     var   ftime     lat     lon 
+#       1       1     214       4       4 
 
 summary(tx_p$data)
-# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 # 13.83   22.08   26.08   26.22   30.72   36.72 
 ```
 
 With the prepared threshold (90th percentile), the WSDI can be computed by running
 
-```
+```r
 WSDI_exp <- CST_TotalSpellTimeExceedingThreshold(tasmax_exp, threshold = tx_p, spell = 6)
 WSDI_obs <- CST_TotalSpellTimeExceedingThreshold(tasmax_obs, threshold = tx_p, spell = 6)
 ```
 
 After checking the summaries, compute the Fair Ranked Probability Skill Score (FRPSS) of WSDI by running the following lines
 
-```
+```r
 # Reorder the data
 fcst <- Reorder(drop(WSDI_exp$data), c(4, 3, 2, 1))
 obs <- Reorder(drop(WSDI_obs$data), c(3, 2, 1))
@@ -559,7 +680,7 @@ summary(fcst)
 
 summary(obs)
 # Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-# 9.00   19.00   22.50   23.16   26.00   33.00 
+# 9.00   19.00   22.50   23.25   26.00   33.00 
  
 # compute FRPSS
 f <- veriApply('FairRpss', fcst = fcst, obs = obs, ensdim = 4, tdim = 3, prob = 1:2/3)$skillscore
@@ -568,7 +689,7 @@ WSDI_FRPSS <- Reorder(f, c(2,1))
 
 Plot the map of WSDI FRPSS for the period from 2013-2016
 
-```
+```r
 cols <- c("#edf8fb", "#ccece6", "#99d8c9", "#66c2a4")
 toptitle <- 'SEAS5 WSDI FRPSS (2013-2016)'