diff --git a/R/AccumulationExceedingThreshold.R b/R/AccumulationExceedingThreshold.R
index b8ae9ae798909c22d8a9ee4b4804ed09af63ac0c..054a8b71a7aca73b3ceb75a94eb7da924aa5673b 100644
--- a/R/AccumulationExceedingThreshold.R
+++ b/R/AccumulationExceedingThreshold.R
@@ -113,7 +113,7 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
                                           na.rm = na.rm, ncores = ncores)
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -128,10 +128,10 @@ CST_AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff =
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -436,6 +436,7 @@ AccumulationExceedingThreshold <- function(data, threshold, op = '>', diff = FAL
                      ncores = ncores)$output1
     }
   }
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/PeriodAccumulation.R b/R/PeriodAccumulation.R
index 8986e3ff97e070c87550fb8c4c840617f1d2b55f..8b0129fe9bcc07ca2f983f02a1183040edde09a9 100644
--- a/R/PeriodAccumulation.R
+++ b/R/PeriodAccumulation.R
@@ -129,12 +129,12 @@ CST_PeriodAccumulation <- function(data, start = NULL, end = NULL,
     data$attrs$Dates <- Dates 
   }
   if (is.null(rollwidth)) {
-    data$coords[[time_dim]] <- NULL
+    data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
     if (!is.null(dim(Dates))) {
       # Create time_bounds
       time_bounds <- NULL
-      time_bounds$start <- Subset(Dates, time_dim, 1, drop = 'selected')
-      time_bounds$end <- Subset(Dates, time_dim, dim(Dates)[time_dim], drop = 'selected')
+      time_bounds$start <- Subset(Dates, time_dim, 1, drop = FALSE)
+      time_bounds$end <- Subset(Dates, time_dim, dim(Dates)[time_dim], drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -321,7 +321,7 @@ PeriodAccumulation <- function(data, dates = NULL, start = NULL, end = NULL,
     pos <- match(dimnames, names(dim(total)))
     total <- aperm(total, pos)
   }
-
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/PeriodMax.R b/R/PeriodMax.R
index 3ae23ec4d519572de99059183cc980fb5e40df57..f038afffe317f2db9d28a1a4d50b1b106dc11656 100644
--- a/R/PeriodMax.R
+++ b/R/PeriodMax.R
@@ -83,8 +83,8 @@ CST_PeriodMax <- function(data, start = NULL, end = NULL,
   
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
-
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
+  
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
       Dates <- SelectPeriodOnDates(dates = Dates, start = start, end = end,
@@ -98,10 +98,10 @@ CST_PeriodMax <- function(data, start = NULL, end = NULL,
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -207,6 +207,8 @@ PeriodMax <- function(data, dates = NULL, start = NULL, end = NULL,
   }
   total <- Apply(list(data), target_dims = time_dim, fun = max,
                  na.rm = na.rm, ncores = ncores)$output1
+  
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/PeriodMean.R b/R/PeriodMean.R
index f58bbeb27c9fc1e12efd8de06b058efd3be737be..bebbedd1047328ba9f27b9ac2a72796b01fddd76 100644
--- a/R/PeriodMean.R
+++ b/R/PeriodMean.R
@@ -83,7 +83,7 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
   
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -98,10 +98,10 @@ CST_PeriodMean <- function(data, start = NULL, end = NULL,
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -206,6 +206,8 @@ PeriodMean <- function(data, dates = NULL, start = NULL, end = NULL,
   }
   total <- Apply(list(data), target_dims = time_dim, fun = mean,
                  na.rm = na.rm, ncores = ncores)$output1
+  
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/PeriodMin.R b/R/PeriodMin.R
index 99ce8016a76fca9d80920abb2b44cf89094ad79d..91f1cd294395621b598dc91b3a229eba577df997 100644
--- a/R/PeriodMin.R
+++ b/R/PeriodMin.R
@@ -83,7 +83,7 @@ CST_PeriodMin <- function(data, start = NULL, end = NULL,
   
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -98,10 +98,10 @@ CST_PeriodMin <- function(data, start = NULL, end = NULL,
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -207,6 +207,7 @@ PeriodMin <- function(data, dates = NULL, start = NULL, end = NULL,
   }
   total <- Apply(list(data), target_dims = time_dim, fun = min,
                  na.rm = na.rm, ncores = ncores)$output1
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/PeriodVariance.R b/R/PeriodVariance.R
index 4e1e93e89181fe40b8aa903ec4df8d7fd6a371f4..5c52466ba07d6b4e2abde7d6f18d05d3c7256988 100644
--- a/R/PeriodVariance.R
+++ b/R/PeriodVariance.R
@@ -87,7 +87,7 @@ CST_PeriodVariance <- function(data, start = NULL, end = NULL,
   
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -102,10 +102,10 @@ CST_PeriodVariance <- function(data, start = NULL, end = NULL,
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -216,6 +216,7 @@ PeriodVariance <- function(data, dates = NULL, start = NULL, end = NULL,
   total <- Apply(list(data), target_dims = time_dim, 
                  fun = .periodvariance,
                  na.rm = na.rm, ncores = ncores)$output1
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/R/TotalSpellTimeExceedingThreshold.R b/R/TotalSpellTimeExceedingThreshold.R
index 14507036558203064f334eeef869b8cc4dd51214..2f420b83446bebc97f4f3036fbc100e99443f1c3 100644
--- a/R/TotalSpellTimeExceedingThreshold.R
+++ b/R/TotalSpellTimeExceedingThreshold.R
@@ -115,7 +115,7 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
                                             ncores = ncores)
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -130,10 +130,10 @@ CST_TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -435,6 +435,7 @@ TotalSpellTimeExceedingThreshold <- function(data, threshold, spell, op = '>',
                      ncores = ncores)$output1
     }
   }
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total) 
 }  
 
diff --git a/R/TotalTimeExceedingThreshold.R b/R/TotalTimeExceedingThreshold.R
index 15b822c6ee37a8aea19c9f0c12d4cd91710cde44..70fc188db5657fac78f11c4a41dfbb750268b79a 100644
--- a/R/TotalTimeExceedingThreshold.R
+++ b/R/TotalTimeExceedingThreshold.R
@@ -119,7 +119,7 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
                                        ncores = ncores)
   data$data <- total
   data$dims <- dim(total)
-  data$coords[[time_dim]] <- NULL
+  data$coords[[time_dim]] <- 1 : length(data$dims[[time_dim]])
 
   if (!is.null(Dates)) {
     if (!is.null(start) && !is.null(end)) {
@@ -134,10 +134,10 @@ CST_TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
       # Create time_bounds
       time_bounds <- NULL
       time_bounds$start <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
-                                                 indices = 1, drop = 'selected')
+                                                 indices = 1, drop = FALSE)
       time_bounds$end <- ClimProjDiags::Subset(x = Dates, along = time_dim, 
                                                indices = dim(Dates)[time_dim], 
-                                               drop = 'selected')
+                                               drop = FALSE)
 
       # Add Dates in attrs
       data$attrs$Dates <- time_bounds$start
@@ -431,6 +431,7 @@ TotalTimeExceedingThreshold <- function(data, threshold, op = '>',
                      ncores = ncores)$output1
     }
   }
+  dim(total) <- c(dim(total), setNames(1, time_dim))
   return(total)
 }
 
diff --git a/tests/testthat/test-AccumulationExceedingThreshold.R b/tests/testthat/test-AccumulationExceedingThreshold.R
index 2b557767f2bcf082de89fa72b19198c2c8b1e7fb..2e579b2ad551efcf1a755b9497f6bff807c0362f 100644
--- a/tests/testthat/test-AccumulationExceedingThreshold.R
+++ b/tests/testthat/test-AccumulationExceedingThreshold.R
@@ -141,36 +141,36 @@ test_that("1. Input checks", {
 test_that("2. Output checks", {
   expect_equal(
     AccumulationExceedingThreshold(dat1, 10),
-    155
+    array(155, c(time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat2_1, 10, time_dim = 'time'), 
-    array(c(375, 390), c(x = 2))
+    array(c(375, 390), c(x = 2, time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat2_1, thres2_1, time_dim = 'time'), 
-    array(c(375, 390), c(x = 2))
+    array(c(375, 390), c(x = 2, time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat2_1, thres2_1, time_dim = 'x'),
-    array(c(rep(0,5), seq(23, 79, 4)), c(time = 20))
+    array(c(rep(0, 5), seq(23, 79, 4)), c(time = 20, x = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat2_1, 10, time_dim = 'time'),
-    array(c(375, 390), c(x = 2))
+    array(c(375, 390), c(x = 2, time = 1))
   )
   # dimensions
   expect_equal(
     dim(AccumulationExceedingThreshold(dat2_3, thres2_3)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(AccumulationExceedingThreshold(dat2_3, thres2_4)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(AccumulationExceedingThreshold(dat2_4, thres2_4, time_dim = 'time')),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
 
 })
@@ -180,27 +180,27 @@ test_that("3. Output checks", {
 
   expect_equal(
     dim(AccumulationExceedingThreshold(dat3_1, c(55,58), c('<', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     AccumulationExceedingThreshold(dat3_1, c(55,58), c(">", "<")),
-    array(c(rep(0,11),113), dim = c(fyear = 3, sdate = 2, lat = 2))
+    array(c(rep(0, 11), 113), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat3_1, c(55,58), c(">=", "<=")),
-    array(c(rep(0,10),55,171), dim = c(fyear = 3, sdate = 2, lat = 2))
+    array(c(rep(0, 10), 55, 171), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat3_2, c(46, 35), op = c("<", ">"), time_dim = 'time'), 
-    array(c(76, 114), c(x = 2))
+    array(c(76, 114), c(x = 2, time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(dat3_2, c(7,11), op = c('>=', '<='), time_dim = 'time'), 
-    array(c(27, 18), c(x = 2))
+    array(c(27, 18), c(x = 2, time = 1))
   )
   expect_equal(
     dim(AccumulationExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), op = c('>=', '<'))),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
 
 })
@@ -211,7 +211,7 @@ test_that("4. Output checks", {
 
   expect_equal(
     dim(AccumulationExceedingThreshold(dat4, list(thres4_2, thres4_1), c('<=', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     as.vector(AccumulationExceedingThreshold(dat4, list(thres4_1, thres4_2), c(">", "<="))[1:3]),
@@ -244,7 +244,7 @@ test_that("5. Seasonal forecasts", {
   res <- CST_AccumulationExceedingThreshold(exp, threshold = 280, time_dim = 'ftime')
 
   expect_equal(
-    round(res$data[, 2, 2, 2]),
+    round(res$data[, 2, 2, 2, 1]),
     c(0, 280, 281, 281)
   )
 
@@ -267,12 +267,12 @@ test_that("5. Seasonal forecasts", {
                                         start = list(1, 4), end = list(31, 10), na.rm = TRUE)
   
   expect_equal(
-    round(GDD[,1,1,1]),
+    round(GDD[, 1, 1, 1, 1]),
     c(549, 387, 125, 554, 245, 282)
   )
   expect_equal(
     dim(GDD),
-    c(member = 6, sdate = 3, lat = 4, lon = 4)
+    c(member = 6, sdate = 3, lat = 4, lon = 4, time = 1)
   )
   expect_error(
     AccumulationExceedingThreshold(exp - 17, threshold = 0, dates = Dates, start = list(1, 4), end = list(31, 10), time_dim = 'ftime'),
@@ -291,19 +291,19 @@ test_that("5. Seasonal forecasts", {
 
   expect_equal(
     AccumulationExceedingThreshold(input_1, threshold_1, diff = TRUE),
-    153
+    array(153, c(time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(input_1, threshold_1),
-    204
+    array(204, c(time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(input_2, threshold_2, op = '<'),
-    -105
+    array(-105, c(time = 1))
   )
   expect_equal(
     AccumulationExceedingThreshold(input_2, threshold_2, op = '<', diff = TRUE),
-    -55
+    array(-55, c(time = 1))
   )
 })
 
@@ -321,7 +321,7 @@ test_that("6. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-PeriodAccumulation.R b/tests/testthat/test-PeriodAccumulation.R
index 9dcbcf9c44cb96b10cbe0c520e875e7ad70dfd8e..3aa47ea1a7da387fbaf5d35c5ea3f0b9f71a8485 100644
--- a/tests/testthat/test-PeriodAccumulation.R
+++ b/tests/testthat/test-PeriodAccumulation.R
@@ -72,12 +72,12 @@ test_that("1. Initial checks", {
   )
   expect_equal(
     PeriodAccumulation(1:10), 
-    55
+    array(55, c(time = 1))
   )
   data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
   expect_equal(
     PeriodAccumulation(data),
-    array(c(9, 12, 27, 30, 45, 48, 63, 66), c(sdate = 2, lon = 4))
+    array(c(9, 12, 27, 30, 45, 48, 63, 66), c(sdate = 2, lon = 4, time = 1))
   )
   # Test dates warning
   expect_warning(
@@ -104,7 +104,7 @@ test_that("2. Seasonal", {
   output$data <- array(c(sum(exp$data[1,1,21:82,1]), sum(exp$data[1,2,21:82,1]),
                          sum(exp$data[1,3,21:82,1]), sum(exp$data[1,1,21:82,2]),
                          sum(exp$data[1,2,21:82,2]), sum(exp$data[1,3,21:82,2])),
-                       c(memb = 1, sdate = 3, lon = 2))
+                       c(memb = 1, sdate = 3, lon = 2, ftime = 1))
 
   expect_equal(
     CST_PeriodAccumulation(exp, start = list(21, 4), end = list(21, 6), 
@@ -127,16 +127,16 @@ test_that("3. Subset Dates and check time_bounds", {
   )
   expect_equal(
     dim(res2$data), 
-    dim(exp$data)[-which(names(dim(exp$data)) == 'ftime')]
+    c(dim(exp$data)[-which(names(dim(exp$data)) == 'ftime')], ftime = 1)
   )
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   expect_equal(
     dim(res2$data)['sdate'], 
-    dim(res2$attrs$Dates)
+    dim(res2$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
@@ -196,29 +196,29 @@ test_that("4. Rolling", {
   # Output checks
   expect_equal(
     PeriodAccumulation(data = dat2, rollwidth = -2, dates = dates2, frequency = 'daily'),
-    array(c(4,6,8, 10, NA, NA), dim = c(sdate = 2, time = 3, member = 1))
+    array(c(4, 6, 8, 10, NA, NA), dim = c(sdate = 2, time = 3, member = 1, time = 1))
   )
   expect_equal(
     PeriodAccumulation(data = dat2, rollwidth = 3, dates = dates2, frequency = 'daily'),
-    array(c(rep(NA, 4), 9, 12), dim = c(sdate = 2, time = 3, member = 1))
+    array(c(rep(NA, 4), 9, 12), dim = c(sdate = 2, time = 3, member = 1, time = 1))
   )
   dat2_1 <- dat2
   dat2_1[1,1,1] <- NA
   expect_equal(
     PeriodAccumulation(data = dat2_1, rollwidth = 2, dates = dates2, na.rm = FALSE, 
                        frequency = 'daily'), 
-    array(c(rep(NA, 3), 6,8,10), dim = c(sdate = 2, time = 3, member = 1))
+    array(c(rep(NA, 3), 6, 8, 10), dim = c(sdate = 2, time = 3, member = 1, time = 1))
   )
   
   # Test rolling with start and end
   expect_equal(
     PeriodAccumulation(data = dat2, rollwidth = 1, dates = dates2, 
                        start = list(1, 4), end = list(2, 4), frequency = 'daily'),
-    array(c(1, 2, 3, 4), dim = c(sdate = 2, time = 2, member = 1))
+    array(c(1, 2, 3, 4), dim = c(sdate = 2, time = 2, member = 1, time = 1))
   )
   expect_equal(
     PeriodAccumulation(data = dat2, rollwidth = 2, dates = dates2, 
                        start = list(1, 4), end = list(2, 4), frequency = 'daily'),
-    array(c(NA, NA, 4, 6), dim = c(sdate = 2, time = 2, member = 1))
+    array(c(NA, NA, 4, 6), dim = c(sdate = 2, time = 2, member = 1, time = 1))
   )
 })
diff --git a/tests/testthat/test-PeriodMax.R b/tests/testthat/test-PeriodMax.R
index 1d7437a40f90f272e7740ef0c4a838bb77cb72ea..ae29bb58e6c6b934f38a0168d8c2db5598050b3a 100644
--- a/tests/testthat/test-PeriodMax.R
+++ b/tests/testthat/test-PeriodMax.R
@@ -23,7 +23,7 @@ test_that("1. Sanity Checks", {
   suppressWarnings(
     expect_equal(
       PeriodMax(array(1, c(x = 1)), time_dim = 'x'), 
-      1
+      array(1, c(x = 1))
     )
   )
   expect_error(
@@ -33,7 +33,7 @@ test_that("1. Sanity Checks", {
   suppressWarnings(
     expect_equal(
       PeriodMax(array(1:10, c(time = 10))), 
-      10
+      array(10, c(time = 1))
     )
   )
   data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
@@ -41,7 +41,7 @@ test_that("1. Sanity Checks", {
     expect_equal(
       PeriodMax(data),
       array(c(5, 6, 11, 12, 17, 18, 23, 24), 
-      c(sdate = 2, lon = 4))
+            c(sdate = 2, lon = 4, time = 1))
     )
   )
   # Test dates warning
@@ -81,7 +81,7 @@ test_that("2. Seasonal", {
   output$data <- array(c(max(exp$data[1,1,21:82,1]), max(exp$data[1,2,21:82,1]),
                          max(exp$data[1,3,21:82,1]), max(exp$data[1,1,21:82,2]),
                          max(exp$data[1,2,21:82,2]), max(exp$data[1,3,21:82,2])),
-                       c(memb = 1, sdate = 3, lon = 2))
+                       c(memb = 1, sdate = 3, lon = 2, time = 1))
   expect_equal(
     CST_PeriodMax(exp, start = list(21, 4), end = list(21, 6))$data,
     output$data
@@ -101,7 +101,7 @@ test_that("3. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-PeriodMean.R b/tests/testthat/test-PeriodMean.R
index 9f8c4cf15faa16229c2ccf223ee3df350c1e72ee..6eab8eed610afc529a1d8d62caf6ac07c00afa87 100644
--- a/tests/testthat/test-PeriodMean.R
+++ b/tests/testthat/test-PeriodMean.R
@@ -10,7 +10,7 @@ test_that("1. Sanity Checks", {
   suppressWarnings(
     expect_equal(
       PeriodMean(array(1, c(x = 1)), time_dim = 'x'), 
-      1
+      array(1, c(x = 1))
     )
   )
   expect_error(
@@ -33,7 +33,7 @@ test_that("1. Sanity Checks", {
   suppressWarnings(
     expect_equal(
       PeriodMean(array(1:10, c(time = 10))), 
-      5.5
+      array(5.5, c(time = 1))
     )
   )
   data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
@@ -41,7 +41,7 @@ test_that("1. Sanity Checks", {
     expect_equal(
       PeriodMean(data),
       array(c(3, 4, 9, 10, 15, 16, 21, 22), 
-      c(sdate = 2, lon = 4))
+            c(sdate = 2, lon = 4, time = 1))
     )
   )
   # Test dates warning
@@ -81,7 +81,7 @@ test_that("2. Seasonal", {
   output$data <- array(c(mean(exp$data[1,1,21:82,1]), mean(exp$data[1,2,21:82,1]),
                          mean(exp$data[1,3,21:82,1]), mean(exp$data[1,1,21:82,2]),
                          mean(exp$data[1,2,21:82,2]), mean(exp$data[1,3,21:82,2])),
-                       c(memb = 1, sdate = 3, lon = 2))
+                       c(memb = 1, sdate = 3, lon = 2, time = 1))
   expect_equal(
     CST_PeriodMean(exp, start = list(21, 4), end = list(21, 6))$data,
     output$data
@@ -101,7 +101,7 @@ test_that("3. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-PeriodMin.R b/tests/testthat/test-PeriodMin.R
index fb97fc206681fcbcb773e42dab4aa5f0af3b9449..51e671b7b87c131034e0e391a42db13712045df7 100644
--- a/tests/testthat/test-PeriodMin.R
+++ b/tests/testthat/test-PeriodMin.R
@@ -10,7 +10,7 @@ test_that("1. Sanity Checks", {
   suppressWarnings(
     expect_equal(
       PeriodMin(array(1, c(x = 1)), time_dim = 'x'), 
-      1
+      array(1, c(x = 1))
     )
   )
   expect_error(
@@ -32,13 +32,13 @@ test_that("1. Sanity Checks", {
   )
   expect_equal(
     PeriodMin(array(1:10, c(time = 10))), 
-    1
+    array(1, c(time = 1))
   )
   data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
   expect_equal(
     PeriodMin(data),
     array(c(1, 2, 7, 8, 13, 14, 19, 20), 
-    c(sdate = 2, lon = 4))
+          c(sdate = 2, lon = 4, time = 1))
   )
   # Test dates warning
   expect_warning(
@@ -77,7 +77,7 @@ test_that("2. Seasonal", {
   output$data <- array(c(min(exp$data[1,1,21:82,1]), min(exp$data[1,2,21:82,1]),
                          min(exp$data[1,3,21:82,1]), min(exp$data[1,1,21:82,2]),
                          min(exp$data[1,2,21:82,2]), min(exp$data[1,3,21:82,2])),
-                       c(memb = 1, sdate = 3, lon = 2))
+                       c(memb = 1, sdate = 3, lon = 2, time = 1))
   expect_equal(
     CST_PeriodMin(exp, start = list(21, 4), end = list(21, 6))$data,
     output$data
@@ -97,7 +97,7 @@ test_that("3. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-PeriodVariance.R b/tests/testthat/test-PeriodVariance.R
index e1de0327cd029d4283e443a5e31780939c50b074..e5c72754a8ad832a5c08939814ead2cb52710a6f 100644
--- a/tests/testthat/test-PeriodVariance.R
+++ b/tests/testthat/test-PeriodVariance.R
@@ -9,7 +9,7 @@ test_that("1. Sanity Checks", {
   )
   expect_equal(
     PeriodVariance(array(1:2, c(x = 2)), time_dim = 'x'), 
-    0.5
+    array(0.5, c(x = 1))
   )
   expect_error(
     PeriodVariance(data = NULL), 
@@ -30,14 +30,14 @@ test_that("1. Sanity Checks", {
   )
   expect_equal(
     PeriodVariance(array(1:10, c(time = 10))), 
-    9.166667,
+    array(9.166667, c(time = 1)),
     tolerance = 0.001
   )
   data <- array(1:24, c(sdate = 2, time = 3, lon = 4))
   expect_equal(
     PeriodVariance(data),
     array(rep(4, 8), 
-    c(sdate = 2, lon = 4))
+          c(sdate = 2, lon = 4, time = 1))
   )
   # Test dates warning
   expect_warning(
@@ -76,7 +76,7 @@ test_that("2. Seasonal", {
   output$data <- array(c(var(exp$data[1,1,21:82,1]), var(exp$data[1,2,21:82,1]),
                          var(exp$data[1,3,21:82,1]), var(exp$data[1,1,21:82,2]),
                          var(exp$data[1,2,21:82,2]), var(exp$data[1,3,21:82,2])),
-                       c(memb = 1, sdate = 3, lon = 2))
+                       c(memb = 1, sdate = 3, lon = 2, time = 1))
   expect_equal(
     CST_PeriodVariance(exp, start = list(21, 4), end = list(21, 6))$data,
     output$data
@@ -96,7 +96,7 @@ test_that("3. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
index e65ec1e8523a3be3317378d895dc2643465d8719..d9a226e47e9710c9f4ef57f04023b78c2fe83d1f 100644
--- a/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalSpellTimeExceedingThreshold.R
@@ -156,36 +156,36 @@ test_that("1. Sanity checks", {
 test_that("2. Output checks", {
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat1, thres1, spell = 2, time_dim = 'time'),  
-    array(c(0,rep(2,9)), c(lat = 10))
+    array(c(0, rep(2, 9)), c(lat = 10, time = 1))
   )    
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat1_2, 10, spell = 2), 
-    array(c(15, 15), c(x = 2))
+    array(c(15, 15), c(x = 2, time = 1))
   )
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat1_2,threshold1_2, spell = 2), 
-    array(c(15, 15), c(x = 2, member = 1))
+    array(c(15, 15), c(x = 2, member = 1, time = 1))
   )
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat2_1, thres2_1, spell = 10), 
-    array(c(15, 15), c(x = 2))
+    array(c(15, 15), c(x = 2, time = 1))
   )
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat2_1, thres2_1, spell = 2, time_dim = 'x'),
-    array(c(rep(0,5), rep(2,15)), c(time = 20))
+    array(c(rep(0, 5), rep(2, 15)), c(time = 20, x = 1))
   )
   # dimensions
   expect_equal(
     dim(TotalSpellTimeExceedingThreshold(dat2_3, thres2_3, spell = 3)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(TotalSpellTimeExceedingThreshold(dat2_3, thres2_4, spell = 3)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(TotalSpellTimeExceedingThreshold(dat2_4, thres2_4, spell = 3, time_dim = 'time')),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
 })
 
@@ -194,28 +194,28 @@ test_that("2. Output checks", {
 test_that("3. Output checks", {
 
   expect_equal(
-    dim(TotalSpellTimeExceedingThreshold(dat3_1, c(55,58), spell = 3, c('<', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    dim(TotalSpellTimeExceedingThreshold(dat3_1, c(55, 58), spell = 3, c('<', '>'))),
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
-    TotalSpellTimeExceedingThreshold(dat3_1, c(30,60), spell = 3, c(">", "<")),
-    array(c(rep(0,6),rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2))
+    TotalSpellTimeExceedingThreshold(dat3_1, c(30, 60), spell = 3, c(">", "<")),
+    array(c(rep(0, 6),rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
-    TotalSpellTimeExceedingThreshold(dat3_1, c(55,58), spell = 3, c(">=", "<=")),
-    array(c(rep(0,11),3), dim = c(fyear = 3, sdate = 2, lat = 2))
+    TotalSpellTimeExceedingThreshold(dat3_1, c(55, 58), spell = 3, c(">=", "<=")),
+    array(c(rep(0, 11),3), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat3_2, c(46, 35), spell = 3, op = c("<", ">"), time_dim = 'time'), 
-    array(c(0, 3), c(x = 2))
+    array(c(0, 3), c(x = 2, time = 1))
   )
   expect_equal(
     TotalSpellTimeExceedingThreshold(dat3_2, c(7,11), spell = 3, op = c('>=', '<='), time_dim = 'time'), 
-    array(c(3, 0), c(x = 2))
+    array(c(3, 0), c(x = 2, time = 1))
   )
   expect_equal(
     dim(TotalSpellTimeExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), spell = 3, op = c('>=', '<'))),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
 
 })
@@ -226,7 +226,7 @@ test_that("4. Output checks", {
 
   expect_equal(
     dim(TotalSpellTimeExceedingThreshold(dat4, list(thres4_2, thres4_1), spell = 3, c('<=', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     as.vector(TotalSpellTimeExceedingThreshold(dat4, list(thres4_1, thres4_2), spell = 3, c(">", "<="))[1:3]),
@@ -254,20 +254,20 @@ test_that("5. Seasonal Forecasts", {
   exp$data <- exp$data[1,1:3,1:3,,,]
   res <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = 260, spell = 2, time_dim = 'ftime')
   expect_equal(
-    res$data[,,1,1],
+    res$data[, , 1, 1, 1],
     array(c(3,3,3,3,3,3,3,3,3), c(member = 3, sdate = 3))
   )
   # compare with percentile
   thresholdP <- Threshold(exp$data, threshold = 0.9)
   WSDI <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = thresholdP, spell = 2, time_dim = 'ftime')
   expect_equal(
-    WSDI$data[3,3,3,],
+    WSDI$data[3, 3, 3, ,1],
     c(rep(0, 6), rep(2, 2), 0, 2, 0, 2, rep(0, 29), 2, rep(0, 11))
   )
   thresholdP1 <- thresholdP[1,,]
   WSDI1 <- CST_TotalSpellTimeExceedingThreshold(exp, threshold = thresholdP1, spell = 2, time_dim = 'ftime')
   expect_equal(
-    WSDI1$data[3,3,3,],
+    WSDI1$data[3, 3, 3, ,1],
     c(rep(0, 53)))
 })
 
@@ -285,7 +285,7 @@ test_that("6. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(
diff --git a/tests/testthat/test-TotalTimeExceedingThreshold.R b/tests/testthat/test-TotalTimeExceedingThreshold.R
index cba27ae03b2ddb5be192a701453f7a8e2cf6b336..53c8a84954f34c8bbb4547a8c00e4e0461c4a33a 100644
--- a/tests/testthat/test-TotalTimeExceedingThreshold.R
+++ b/tests/testthat/test-TotalTimeExceedingThreshold.R
@@ -149,36 +149,36 @@ test_that("1. Sanity checks", {
 test_that("2. Output checks", {
   expect_equal(
     TotalTimeExceedingThreshold(dat1, thres1, time_dim = 'time'),  
-    array(c(0,rep(2,9)), c(lat = 10))
+    array(c(0, rep(2, 9)), c(lat = 10, time = 1))
   )    
   expect_equal(
     TotalTimeExceedingThreshold(dat1_2, 10), 
-    array(c(15, 15), c(x = 2))
+    array(c(15, 15), c(x = 2, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat1_2,threshold1_2), 
-    array(c(15, 15), c(x = 2, member = 1))
+    array(c(15, 15), c(x = 2, member = 1, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat2_1, thres2_1), 
-    array(c(15, 15), c(x = 2))
+    array(c(15, 15), c(x = 2, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat2_1, thres2_1, time_dim = 'x'),
-    array(c(rep(0,5), rep(2,15)), c(time = 20))
+    array(c(rep(0, 5), rep(2, 15)), c(time = 20, x = 1))
   )
   # dimensions
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat2_3, thres2_3)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat2_3, thres2_4)),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat2_4, thres2_4, time_dim = 'time')),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
 })
 
@@ -187,27 +187,27 @@ test_that("2. Output checks", {
 test_that("3. Output checks", {
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat3_1, c(55,58), c('<', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat3_1, c(30,60), c(">", "<")),
-    array(c(rep(0, 6), rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2))
+    array(c(rep(0, 6), rep(5, 5), 4), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat3_1, c(55, 58), c(">=", "<=")),
-    array(c(rep(0, 10), 1, 3), dim = c(fyear = 3, sdate = 2, lat = 2))
+    array(c(rep(0, 10), 1, 3), dim = c(fyear = 3, sdate = 2, lat = 2, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat3_2, c(46, 35), op = c("<", ">"), time_dim = 'time'), 
-    array(c(2, 3), c(x = 2))
+    array(c(2, 3), c(x = 2, time = 1))
   )
   expect_equal(
     TotalTimeExceedingThreshold(dat3_2, c(7, 11), op = c('>=', '<='), time_dim = 'time'), 
-    array(c(3, 2), c(x = 2))
+    array(c(3, 2), c(x = 2, time = 1))
   )
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat2_3, list(thres2_4, thres2_4+1), op = c('>=', '<'))),
-    c(sdate = 2, lat = 2)
+    c(sdate = 2, lat = 2, time = 1)
   )
 })
 
@@ -216,7 +216,7 @@ test_that("3. Output checks", {
 test_that("4. Output checks", {
   expect_equal(
     dim(TotalTimeExceedingThreshold(dat4, list(thres4_2, thres4_1), c('<=', '>'))),
-    c(fyear = 3, sdate = 2, lat = 2)
+    c(fyear = 3, sdate = 2, lat = 2, time = 1)
   )
   expect_equal(
     as.vector(TotalTimeExceedingThreshold(dat4, list(thres4_1, thres4_2), c(">", "<="))[1:3]),
@@ -245,7 +245,7 @@ test_that("5. Seasonal forecasts", {
   exp$data <- exp$data[1, 1:3, , , , ] - 247
   SU35_NoP <- CST_TotalTimeExceedingThreshold(exp, threshold = 35, time_dim = 'ftime')$data  
   expect_equal(
-    SU35_NoP[1, , 15, 3], 
+    SU35_NoP[1, , 15, 3, 1], 
     c(0, 1, 1, 1, 0, 0)
   )
   # convert to percentile
@@ -255,7 +255,7 @@ test_that("5. Seasonal forecasts", {
   data$data <- exp_percentile
   SU35_P <- CST_TotalTimeExceedingThreshold(data, threshold = obs_percentile, time_dim = 'ftime')$data
   expect_equal(
-    SU35_P[2, , 5, 5], 
+    SU35_P[2, , 5, 5, 1], 
     c(3, 3, 3, 3, 3, 3)
   )
 })
@@ -274,7 +274,7 @@ test_that("6. Subset Dates and check time_bounds", {
   # Check Dates
   expect_equal(
     dim(res$data)['sdate'], 
-    dim(res$attrs$Dates)
+    dim(res$attrs$Dates)['sdate']
   )
   # Check time_bounds
   expect_equal(