diff --git a/DESCRIPTION b/DESCRIPTION
index 7fcf115d0c7fc117938b624e467bb914e7388602..9d6c9bd992b78770f4aef6fdc552e90a447686f6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -28,8 +28,7 @@ Depends:
 Imports:
     multiApply (>= 2.1.1),
     s2dv,
-    stats,
-    ClimProjDiags
+    stats
 Suggests:
     testthat,
     CSTools,
diff --git a/NAMESPACE b/NAMESPACE
index 133942a9fa0a8c645c323c2b3755abdc386efe8d..0fa3f245d7c6bb9bea3ec06c4d158233e2f7f4d8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,7 +26,6 @@ export(TotalTimeExceedingThreshold)
 export(WindCapacityFactor)
 export(WindPowerDensity)
 import(multiApply)
-importFrom(ClimProjDiags,Subset)
 importFrom(s2dv,InsertDim)
 importFrom(s2dv,Reorder)
 importFrom(stats,approxfun)
diff --git a/R/MergeRefToExp.R b/R/MergeRefToExp.R
index fa3dcaf8415b8a36c659b36f0247bd74bb87b308..17868263bb8c6bffa40315462337ccd88f0a14d8 100644
--- a/R/MergeRefToExp.R
+++ b/R/MergeRefToExp.R
@@ -60,7 +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,
@@ -201,7 +200,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, 
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/zzz.R b/R/zzz.R
index 9b0c6488fb9a91a5c8b5e86bffaa454a7be48379..5e20873cf754a33700167f34801e98444427d854 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -24,6 +24,17 @@
   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)
+}
+
 
 #=======================
 # Read a powercurve file
@@ -60,4 +71,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/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)
+  )
+
 })
 
 ##############################################