diff --git a/R/QThreshold.R b/R/QThreshold.R
index 0f20858fc8dd989797aa0759bf3205502148201d..19cda3c64cca5003abbb294bdceda1a72cf0d9fd 100644
--- a/R/QThreshold.R
+++ b/R/QThreshold.R
@@ -177,6 +177,7 @@ CST_QThreshold <- function(data, threshold, start = NULL, end = NULL,
 #'                      start = list(21, 4), end = list(21, 6))
 #' 
 #'@import multiApply
+#'@importFrom ClimProjDiags Subset
 #'@export
 QThreshold <- function(data, threshold, dates = NULL, start = NULL, end = NULL,
                        time_dim = 'time', memb_dim = 'member', 
@@ -231,7 +232,7 @@ 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% names(dim(dates))) {
-            dates_thres <- .arraysubset(dates, dim = sdate_dim, value = 1)
+            dates_thres <- Subset(dates, along = sdate_dim, indices = 1)
             threshold <- SelectPeriodOnData(data = threshold, dates = dates_thres, start, end,
                                             time_dim = time_dim, ncores = ncores)
           } else {
diff --git a/R/SelectPeriodOnData.R b/R/SelectPeriodOnData.R
index 016133be313cfae097ddfd2f63abd9386d34129d..bef70be2399f9f9256e91635056a7e63fe2864d3 100644
--- a/R/SelectPeriodOnData.R
+++ b/R/SelectPeriodOnData.R
@@ -100,6 +100,7 @@ CST_SelectPeriodOnData <- function(data, start, end, time_dim = 'time',
 #'dim(Dates) <- c(time = 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 = 'time', ncores = NULL) {
@@ -149,11 +150,10 @@ 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)]
-    res <- .arraysubset(res, dim = dim_remove, value = 1)
-    dim(res) <- dim(res)[-which(names(dim(res)) %in% dim_remove)]
+    indices <- as.list(rep(1, length(dim_remove)))
+    res <- Subset(res, along = dim_remove, indices, drop = 'selected')
   }
-
   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 cf9163970f76c2da74b3d811c1fa0b71445beeab..0724f066e6c2ce674c779159f96ccee4ff38576c 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -24,35 +24,6 @@
   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
 # Create the approximation function