diff --git a/DESCRIPTION b/DESCRIPTION
index 7fcf115d0c7fc117938b624e467bb914e7388602..90e0e83155393c8d19d68ea8c92b9c0121f6153f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: CSIndicators
 Title: Climate Services' Indicators Based on Sub-Seasonal to Decadal Predictions
-Version: 1.0.0
+Version: 1.0.1
 Authors@R: c(
   person("Eva", "Rifà", , "eva.rifarovira@bsc.es", role = c("cre")),
   person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut"), comment = c(ORCID = "0000-0001-8568-3071")),
@@ -27,9 +27,7 @@ Depends:
     R (>= 3.6.0)
 Imports:
     multiApply (>= 2.1.1),
-    s2dv,
-    stats,
-    ClimProjDiags
+    stats
 Suggests:
     testthat,
     CSTools,
diff --git a/NAMESPACE b/NAMESPACE
index 133942a9fa0a8c645c323c2b3755abdc386efe8d..d80accbbe3b6781acb27cb8c953a9c6b00a1824a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,9 +26,6 @@ export(TotalTimeExceedingThreshold)
 export(WindCapacityFactor)
 export(WindPowerDensity)
 import(multiApply)
-importFrom(ClimProjDiags,Subset)
-importFrom(s2dv,InsertDim)
-importFrom(s2dv,Reorder)
 importFrom(stats,approxfun)
 importFrom(stats,ecdf)
 importFrom(stats,quantile)
diff --git a/NEWS.md b/NEWS.md
index e28aab8829eab58cde5a114901d086598396db80..44285d2f5c49d8a63ca79ed8ce235d6ee3c1dacd 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,16 +1,22 @@
-# CSIndicators 1.0.0 (Release date: 2023-04-05)
+# CSIndicators 1.0.1 (Release date: 2023-05-18)  
 **Fixes**  
-- Correct vignettes figures links.
+- Add EnergyIndicators vignette figures  
+- Remove ClimProjDiags dependency  
+- Remove s2dv dependency  
+
+# CSIndicators 1.0.0 (Release date: 2023-04-05)  
+**Fixes**  
+- Correct vignettes figures links.  
 
 **New features**  
-- Exceeding Threshold functions to allow between thresholds or equal threshold options.
-- New s2dv_cube object development for all the functions, unit tests, examples and vignettes.
+- Exceeding Threshold functions to allow between thresholds or equal threshold options.  
+- New s2dv_cube object development for all the functions, unit tests, examples and vignettes.  
 
-# CSIndicators 0.0.2 (Release date: 2022-10-21)
+# CSIndicators 0.0.2 (Release date: 2022-10-21)  
 **Fixes**  
-- Correct figures of EnergyIndicators vignette.
-- Sanity check correction in functions CST_PeriodAccumulation, CST_AbsToProbs, CST_AccumulationExceedingThreshold, CST_MergeRefToExp, CST_PeriodMean, CST_QThreshold, CST_SelectPeriodOnData, CST_Threshold, TotalSpellTimeExceedingThreshold, CST_TotalTimeExceedingThreshold, CST_WindCapacityFactor and CST_WindPowerDensity.
-- Revise examples using s2dv::InsertDim in MergeRefToExp().
+- Correct figures of EnergyIndicators vignette.  
+- Sanity check correction in functions CST_PeriodAccumulation, CST_AbsToProbs, CST_AccumulationExceedingThreshold, CST_MergeRefToExp, CST_PeriodMean, CST_QThreshold, CST_SelectPeriodOnData, CST_Threshold, TotalSpellTimeExceedingThreshold, CST_TotalTimeExceedingThreshold, CST_WindCapacityFactor and CST_WindPowerDensity.  
+- Revise examples using s2dv::InsertDim in MergeRefToExp().  
 
-# CSIndicators 0.0.1 (Release date: 2021-05-07)
-- This package is intended for sub-seasonal, seasonal and decadal climate predictions, but its methods are also applicable to other time-scales. Additionally, the outputs of the functions in this package are compatible with 'CSTools'. This package was developed in the context of H2020 MED-GOLD (776467) and S2S4E (776787) projects. Lledó et al. (2019) <doi:10.1016/j.renene.2019.04.135>.
+# CSIndicators 0.0.1 (Release date: 2021-05-07)  
+- This package is intended for sub-seasonal, seasonal and decadal climate predictions, but its methods are also applicable to other time-scales. Additionally, the outputs of the functions in this package are compatible with 'CSTools'. This package was developed in the context of H2020 MED-GOLD (776467) and S2S4E (776787) projects. Lledó et al. (2019) <doi:10.1016/j.renene.2019.04.135>.  
diff --git a/R/AccumulationExceedingThreshold.R b/R/AccumulationExceedingThreshold.R
index 7fd78f49f9de6d3fad3a65043ed7a46f2a19e143..e346b5310504a9e0a9b0fbbf289d332d4cbeb8cb 100644
--- a/R/AccumulationExceedingThreshold.R
+++ b/R/AccumulationExceedingThreshold.R
@@ -170,7 +170,6 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
 #'GDD <- AccumulationExceedingThreshold(data, threshold = 0, start = list(1, 4),
 #'                                      end = list(31, 10))
 #'@import multiApply
-#'@importFrom s2dv Reorder
 #'@export
 AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FALSE,
                                            dates = NULL, start = NULL, end = NULL,
diff --git a/R/MergeRefToExp.R b/R/MergeRefToExp.R
index fa3dcaf8415b8a36c659b36f0247bd74bb87b308..434cae35c0db9d2d0312c02473a0b5a8ef490078 100644
--- a/R/MergeRefToExp.R
+++ b/R/MergeRefToExp.R
@@ -60,8 +60,6 @@
 #'                              start2 = list(1, 7), end2 = list(21, 9))
 #' 
 #'@import multiApply
-#'@importFrom ClimProjDiags Subset
-#'@importFrom s2dv InsertDim
 #'@export
 CST_MergeRefToExp <- function(data1, data2, start1, end1, start2, end2,
                               time_dim = 'ftime', sdate_dim = 'sdate',
@@ -201,8 +199,6 @@ CST_MergeRefToExp <- function(data1, data2, start1, end1, start2, end2,
 #'                          time_dim = 'time')
 #'
 #'@import multiApply
-#'@importFrom ClimProjDiags Subset
-#'@importFrom s2dv InsertDim
 #'@export
 MergeRefToExp <- function(data1, dates1, start1, end1, data2, dates2, start2, 
                           end2, time_dim = 'ftime', sdate_dim = 'sdate',
@@ -252,8 +248,8 @@ MergeRefToExp <- function(data1, dates1, start1, end1, data2, dates2, start2,
     dif_dims <- which(names(dim(data2)) %in% names(dim(data1)) == FALSE)
     if (length(dif_dims) > 0) {
       for (i in dif_dims) {
-        data1 <- InsertDim(data1, posdim = i, lendim = dim(data2)[i], 
-                           name = names(dim(data2))[i])
+        data1 <- .insertdim(data1, posdim = i, lendim = dim(data2)[i], 
+                            name = names(dim(data2))[i])
       }
     }
   }
@@ -262,8 +258,8 @@ MergeRefToExp <- function(data1, dates1, start1, end1, data2, dates2, start2,
     dif_dims <- which(names(dim(data1)) %in% names(dim(data2)) == FALSE)
     if (length(dif_dims) > 0) {
       for (i in dif_dims) {
-        data2 <- InsertDim(data2, posdim = i, lendim = dim(data1)[i],
-                           name = names(dim(data1))[i])
+        data2 <- .insertdim(data2, posdim = i, lendim = dim(data1)[i],
+                            name = names(dim(data1))[i])
       }
     }
   }
diff --git a/R/QThreshold.R b/R/QThreshold.R
index e86b95a0dcb634326356a6283db8d837e1b3c185..49217dd20e8986ecb786d8f4232212743132a989 100644
--- a/R/QThreshold.R
+++ b/R/QThreshold.R
@@ -68,7 +68,6 @@
 #'exp_probs <- CST_QThreshold(exp, threshold)
 #' 
 #'@import multiApply
-#'@importFrom ClimProjDiags Subset
 #'@export
 CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
                            time_dim = 'ftime', memb_dim = 'member', 
@@ -162,11 +161,12 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
 #'thres_q <- QThreshold(data, threshold)
 #' 
 #'@import multiApply
-#'@importFrom ClimProjDiags Subset
 #'@export
 QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
                        time_dim = 'time', memb_dim = 'member', sdate_dim = 'sdate',
                        ncores = NULL) {
+  # Initial checks
+  ## data
   if (is.null(data)) {
     stop("Parameter 'data' cannot be NULL.")
   }
@@ -177,6 +177,10 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
     dim(data) <- c(length(data), 1)
     names(dim(data)) <- c(memb_dim, sdate_dim)
   }
+  if (is.null(names(dim(data)))) {
+    stop("Parameter 'data' must have named dimensions.")
+  }
+  ## threshold
   if (is.null(threshold)) {
       stop("Parameter 'threshold' cannot be NULL.")
   }
@@ -189,8 +193,8 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
   } else if (length(threshold) == 1) {
     dim(threshold) <- NULL
   }
-  if (is.null(names(dim(data)))) {
-    stop("Parameter 'data' must have named dimensions.")
+  if (sdate_dim %in% names(dim(threshold))) {
+    stop("Parameter threshold cannot have dimension 'sdate_dim'.") 
   }
   if (is.null(names(dim(threshold))) && length(threshold) > 1) {
     stop("Parameter 'threshold' must have named dimensions.")
@@ -206,9 +210,9 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
       }
       if (time_dim %in% names(dim(threshold))) {
         if (dim(threshold)[time_dim] == dim(data)[time_dim]) {
-          if (!is.null(dim(dates)) && sdate_dim %in% dim(dates)) {
-            dates_thres <- Subset(dates, along = sdate_dim, indices = 1) 
-            threshold <- SelectPeriodOnData(threshold, dates_thres, start, end,
+          if (!is.null(dim(dates)) && sdate_dim %in% names(dim(dates))) {
+            dates_thres <- .arraysubset(dates, dim = sdate_dim, value = 1)
+            threshold <- SelectPeriodOnData(data = threshold, dates = dates_thres, start, end,
                                             time_dim = time_dim, ncores = ncores)
           } else {
             threshold <- SelectPeriodOnData(threshold, dates, start, end,
@@ -231,9 +235,7 @@ QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
     }
   } else {
     target_thres <- NULL
-    if (sdate_dim %in% names(dim(threshold))) {
-      stop("Parameter threshold cannot have dimension 'sdate_dim'.") 
-    }
+
     if (memb_dim %in% names(dim(data))) {
       if (memb_dim %in% names(dim(threshold))) {
       # comparison member by member
diff --git a/R/SelectPeriodOnData.R b/R/SelectPeriodOnData.R
index b9cf8ac3966bca04dee0c806b0c2d3421de259f2..94bcfe9f3daa66e63e8a3375369c279d2279a2fc 100644
--- a/R/SelectPeriodOnData.R
+++ b/R/SelectPeriodOnData.R
@@ -100,7 +100,6 @@ CST_SelectPeriodOnData <- function(data, start, end, time_dim = 'ftime',
 #'dim(Dates) <- c(ftime = 214, sdate = 3)
 #'Period <- SelectPeriodOnData(data, Dates, start = list(21, 6), end = list(21, 9))
 #'@import multiApply
-#'@importFrom ClimProjDiags Subset
 #'@export
 SelectPeriodOnData <- function(data, dates, start, end, 
                                time_dim = 'ftime', ncores = NULL) {
@@ -150,12 +149,11 @@ SelectPeriodOnData <- function(data, dates, start, end,
   names_data <- sort(names(dim(data)))
   if (!all(names_res %in% names_data)) {
     dim_remove <- names_res[-which(names_res %in% names_data)]
-    indices <- as.list(rep(1, length(dim_remove)))
-    res <- Subset(res, along = dim_remove, indices, drop = 'selected')
+    res <- .arraysubset(res, dim = dim_remove, value = 1)
+    dim(res) <- dim(res)[-which(names(dim(res)) %in% dim_remove)]
   }
 
   pos <- match(names(dim(data)), names(dim(res)))
   res <- aperm(res, pos)
   return(res)
-}
-
+}
\ No newline at end of file
diff --git a/R/SelectPeriodOnDates.R b/R/SelectPeriodOnDates.R
index 09633dd605100f8fb964ad30acb1792d3f227cac..fcb1a4cb752deb02461ea33499e9e67d2efc1e7b 100644
--- a/R/SelectPeriodOnDates.R
+++ b/R/SelectPeriodOnDates.R
@@ -20,7 +20,6 @@
 #'the vector dates during the period requested from \code{start} to \code{end}.
 #'
 #'@import multiApply
-#'@importFrom s2dv Reorder
 #'
 #'@examples
 #'Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
@@ -66,7 +65,8 @@ SelectPeriodOnDates <- function(dates, start, end,
     res <- as.POSIXct(res, origin = '1970-01-01', tz = 'UTC')  
   } else {
     if (!all(names(dim(res)) == names(dim(dates)))) {
-      res <- s2dv::Reorder(res, names(dim(dates)))
+      pos <- match(names(dim(dates)), names(dim(res)))
+      res <- aperm(res, pos)
     }  
     res <- dates[res]
     dim(res) <- dims
diff --git a/R/zzz.R b/R/zzz.R
index 9b0c6488fb9a91a5c8b5e86bffaa454a7be48379..cf9163970f76c2da74b3d811c1fa0b71445beeab 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -24,6 +24,34 @@
   return(position)
 }
 
+# Function to subset dimension indices of an array
+.arraysubset <- function(x, dim, value, drop = FALSE) { 
+  indices <- rep(list(bquote()), length(dim(x)))
+  if (is.character(dim)) {
+    dim <- which(names(dim(x)) %in% dim)
+  }
+  indices[dim] <- value
+  call <- as.call(c(list(as.name("["), quote(x)), indices, drop = drop))
+  eval(call)
+}
+
+# Function to insert a dimension in an array
+.insertdim <- function(data, posdim, lendim, name = NULL) {
+  names(lendim) <- name
+  data <- array(data, dim = c(dim(data), lendim))
+  ## Reorder dimension
+  if (posdim == 1) {
+    order <- c(length(dim(data)), 1:(length(dim(data)) - 1))
+    data <- aperm(data, order)
+  } else if (posdim == length(dim(data))) {  # last dim
+
+  } else { # middle dim
+    order <- c(1:(posdim - 1), length(dim(data)), posdim:(length(dim(data)) - 1))
+    data <- aperm(data, order)
+  }
+  return(data)
+}
+
 
 #=======================
 # Read a powercurve file
@@ -60,4 +88,4 @@ wind2CF <- function(wind, pc) {
   power <- wind2power(wind, pc)
 	CF <- power / pc$attr$RatedPower
 	return(CF)
-}
+}
\ No newline at end of file
diff --git a/inst/doc/paper-figure-PlotForecastPDF.R b/inst/doc/paper-figure-PlotForecastPDF.R
index b8867a61c882f4cb6453a64226c0539aae74aa8a..bebc49725a02dec4c256a398974c532705d3833a 100644
--- a/inst/doc/paper-figure-PlotForecastPDF.R
+++ b/inst/doc/paper-figure-PlotForecastPDF.R
@@ -34,9 +34,9 @@ c(hcst, hcst_ref) %<-% CST_Load(var = 'prlr',
                                 leadtimemin = 1, leadtimemax = 214, 
                                 nmember = 25, output = "lonlat") 
 hcst$data <- hcst$data * 3600 * 24 * 1000
-attributes(hcst$attrs$Variable)$units <- 'mm'
+hcst$attrs$Variable$metadata$prlr$units <- 'mm'
 hcst_ref$data <- hcst_ref$data * 3600 * 24 * 1000
-attributes(hcst_ref$attrs$Variable)$units <- 'mm'
+hcst_ref$attrs$Variable$metadata$prlr$units <- 'mm'
 
 
 c(fcst, obs) %<-% CST_Load(var = 'prlr',
@@ -49,9 +49,9 @@ c(fcst, obs) %<-% CST_Load(var = 'prlr',
                            leadtimemin = 1, leadtimemax = 214,
                            nmember = 50, output = "lonlat")
 fcst$data <- fcst$data * 1000 * 3600 * 24
-attributes(fcst$attrs$Variable)$units <- 'mm'
+fcst$attrs$Variable$metadata$prlr$units <- 'mm'
 obs$data <- obs$data  * 1000 * 3600 * 24
-attributes(obs$attrs$Variable)$units <- 'mm'
+obs$attrs$Variable$metadata$prlr$units <- 'mm'
 
 
 fcst_QM <- CST_QuantileMapping(exp = hcst,
diff --git a/tests/testthat/test-QThreshold.R b/tests/testthat/test-QThreshold.R
index 7572bd065e95930e3e427db62fde411e519362d1..41cc3e5312d1e7d56d698bd0c1f6432772e06608 100644
--- a/tests/testthat/test-QThreshold.R
+++ b/tests/testthat/test-QThreshold.R
@@ -87,6 +87,25 @@ test_that("Sanity checks", {
     c(sdate = 2, time = 5, member = 3, lat = 2)
   )
 
+  # test different common dimensions
+  
+  exp <- array(1:61, dim = c(ftime = 61, sdate = 3))
+  threshold <- array(1:61, dim = c(ftime = 61))
+  Dates <- c(seq(as.Date("01-05-2000", format = "%d-%m-%Y"), 
+                as.Date("30-06-2000", format = "%d-%m-%Y"), by = 'day'),
+            seq(as.Date("01-05-2001", format = "%d-%m-%Y"), 
+                as.Date("30-06-2001", format = "%d-%m-%Y"), by = 'day'),
+            seq(as.Date("01-05-2002", format = "%d-%m-%Y"), 
+                as.Date("30-06-2002", format = "%d-%m-%Y"), by = 'day'))
+  dim(Dates) <- c(ftime = 61, sdate = 3, syear = 1)
+  res <- QThreshold(data = exp, dates = Dates, 
+                    start = list(21, 4), end = list(21, 6), threshold = threshold, 
+                    time_dim = 'ftime', sdate_dim = 'sdate')
+  expect_equal(
+    dim(res),
+    c(sdate = 3, ftime = 52)
+  )
+
 })
 
 ##############################################
diff --git a/vignettes/EnergyIndicators.Rmd b/vignettes/EnergyIndicators.Rmd
index f4a1a04b722a6b250efadc015808dc5691811969..5a65c4759ee4932b76b791e1f3787fa6dcafaf25 100644
--- a/vignettes/EnergyIndicators.Rmd
+++ b/vignettes/EnergyIndicators.Rmd
@@ -30,24 +30,38 @@ Although wind turbines cannot extract all of the kinetic energy in the wind, and
 As an example, we simulate a time series of 1000 wind speed values from a Weibull distribution with scale factor of 6 and a shape factor of 2, which represent a sample of wind speed values obtained at a single location. The Weibull distribution is often assumed to fit observed wind speed values to a probability distribution function. Then, each instantaneous wind speed value is converted to its equivalent WPD.
 The `mean` and `sd` of the WPD can be employed to summarize the wind resource in that location. Otherwise, we can plot the histograms to see the full distribution of values:
 
-```{r, fig.width=7}
+```
 library(CSIndicators)
 set.seed(1)
-oldpar <- par(no.readonly = TRUE) 
 wind <- rweibull(n = 1000, shape = 2, scale = 6)
 WPD <- WindPowerDensity(wind)
 mean(WPD)
+```
+
+```
+## [1] 170.6205
+```
+
+```
 sd(WPD)
+```
+
+```
+## [1] 251.1349
+```
+
+```
 par(mfrow = c(1, 2))
 hist(wind, breaks = seq(0, 20))
 hist(WPD, breaks = seq(0, 4000, 200))
 ```
+![WPD](./Figures/WPD_histogram.png)
 
 As you can see the histogram of the WPD is highly skewed, even if the wind speed was only a little skewed!
 
 If not specified, an air density of 1.225 kg/m^3 is assumed. Otherwise, the parameter `ro` can be set to a fixed value (for instance the mean air density at the site elevation could be used), or a timeseries of density values measured at each time stamp can be used to obtain more accurate results.
 
-```{r}
+```
 WPD <- WindPowerDensity(wind, ro = 1.15)
 ```
 
@@ -61,16 +75,16 @@ Notice that power curves are intended to be used with 10-minutal steady wind spe
 
 Following on the previous example, we will compute now the CF that would be obtained from our sample of 1000 wind speed values when using a turbine of class IEC I, and compare it to the CF values for a class III:
 
-```{r, fig.width=7}
+```
 WCFI <- WindCapacityFactor(wind, IEC_class = "I")
 WCFIII <- WindCapacityFactor(wind, IEC_class = "III")
 par(mfrow = c(1, 3))
 hist(wind, breaks = seq(0, 20))
 hist(WCFI, breaks = seq(0, 1, 0.05), ylim = c(0, 500))
 hist(WCFIII, breaks = seq(0, 1, 0.05), ylim = c(0, 500))
-par(oldpar) 
 ```
 
+![WCF](./Figures/WCF_histogram.png)
 
 From the CF histograms we can see that, for this particular wind speed distribution, the IEC I turbine (designed for high winds) producess less energy than the IEC III turbine, which is more suitable for this range of wind speed values.