diff --git a/.github/ArgumentUsage.xlsx b/.github/ArgumentUsage.xlsx
index 6a64457..9fb1cf4 100644
Binary files a/.github/ArgumentUsage.xlsx and b/.github/ArgumentUsage.xlsx differ
diff --git a/.gitignore b/.gitignore
index c34e6e9..6a85405 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,4 +4,7 @@
 .httr-oauth
 .DS_Store
 .quarto
-.vscode
\ No newline at end of file
+.vscode
+
+# ignore KiN nc files for now as most are too big to ship the package with
+ExampleData/*
\ No newline at end of file
diff --git a/NAMESPACE b/NAMESPACE
index 743e589..adfd710 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -10,6 +10,7 @@ export(Discovery_Library)
 export(Discovery_QuickFacts)
 export(Discovery_Read)
 export(Discovery_Variables)
+export(Metrics_BootstrapQuantiles)
 export(Metrics_ETCCDI)
 export(NC_Read)
 export(NC_Write)
@@ -51,6 +52,8 @@ importFrom(ncdf4,ncatt_put)
 importFrom(ncdf4,ncdim_def)
 importFrom(ncdf4,ncvar_def)
 importFrom(ncdf4,ncvar_put)
+importFrom(ncdfCF,as_CF)
+importFrom(ncdfCF,create_ncdf)
 importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(progress,progress_bar)
diff --git a/R/Discovery_Functionality.r b/R/Discovery_Functionality.r
index 82a5b47..01b0d72 100644
--- a/R/Discovery_Functionality.r
+++ b/R/Discovery_Functionality.r
@@ -48,18 +48,24 @@ Discovery_DataSet.Check <- function(dataSet) {
 #'
 #' @export
 Discovery_Read <- function(dataSet = "NULL") {
-    owner <- "Clim-Hub"
-    repo <- "ClimHub"
-    branch <- "master"
-    path <- "product-metadata"
+    metadata_dir <- system.file("product-metadata", package = "ClimHub")
+    if (metadata_dir == "") {
+        stop("Could not locate 'product-metadata' directory in package 'ClimHub'.")
+    }
 
     Discovery_DataSet.Check(dataSet)
 
-    # URL of the raw JSON file
-    URL <- paste0("https://raw.githubusercontent.com/", owner, "/", repo, "/", branch, "/", path, "/", paste0(dataSet, ".json"))
+    json_file <- file.path(metadata_dir, paste0(dataSet, ".json"))
+
+    # URL-based fallback kept for reference:
+    # owner <- "Clim-Hub"
+    # repo <- "ClimHub"
+    # branch <- "master"
+    # path <- "product-metadata"
+    # URL <- paste0("https://raw.githubusercontent.com/", owner, "/", repo, "/", branch, "/", path, "/", paste0(dataSet, ".json"))
 
-    # Read JSON content
-    json_data <- jsonlite::fromJSON(URL)
+    # Read JSON content from package files.
+    json_data <- jsonlite::fromJSON(json_file)
 
     # View the data
     json_data
@@ -129,28 +135,21 @@ Discovery_Citation <- function(dataSet = "NULL") {
 #'
 #' @export
 Discovery_Library <- function() {
-    owner <- "Clim-Hub"
-    repo <- "ClimHub"
-    path <- "product-metadata"
-
-    # Construct the URL
-    URL <- paste0("https://github.com/", owner, "/", repo, "/tree/master/", path)
-
-    # Read the page
-    page <- rvest::read_html(URL)
-
-    # Extract the embedded JSON metadata
-    raw_metadata <- page %>%
-        rvest::html_nodes("script") %>% # Extract all <script> tags
-        rvest::html_text() %>% # Get the text content
-        grep("payload", ., value = TRUE) %>% # Locate the script containing "payload"
-        jsonlite::fromJSON(simplifyVector = FALSE) # Parse JSON metadata
+    metadata_dir <- system.file("product-metadata", package = "ClimHub")
+    if (metadata_dir == "") {
+        stop("Could not locate 'product-metadata' directory in package 'ClimHub'.")
+    }
 
-    # Navigate to the relevant "tree" section
-    tree_items <- raw_metadata$payload$codeViewFileTreeLayoutRoute$fileTree$`product-metadata`$items
+    tree_items <- list.files(path = metadata_dir, pattern = "\\.json$", full.names = FALSE)
 
-    # make into vector of names of files
-    tree_items <- unlist(lapply(tree_items, "[[", "name"))
+    # GitHub Contents API approach kept for reference:
+    # owner <- "Clim-Hub"
+    # repo <- "ClimHub"
+    # branch <- "master"
+    # path <- "product-metadata"
+    # URL <- paste0("https://api.github.com/repos/", owner, "/", repo, "/contents/", path, "?ref=", branch)
+    # tree_items <- jsonlite::fromJSON(URL, simplifyVector = TRUE)
+    # tree_items <- tree_items$name
 
     # subset for .json files
     tree_items <- tree_items[grep(pattern = "\\.json$", tree_items)]
diff --git a/R/Helper_ETCCDIDailyBootstraps.r b/R/Helper_ETCCDIDailyBootstraps.r
new file mode 100644
index 0000000..ee474f2
--- /dev/null
+++ b/R/Helper_ETCCDIDailyBootstraps.r
@@ -0,0 +1,49 @@
+#' @title Prepare Bootstrap Intervals of Daily Time Series across Years
+#'
+#' @description Create bootstrap intervals of daily time series across years for a given time sequence. This function is needed for baseperiod quantile calculations for `Metrics_ETCCDI`.
+#'
+#' @param dates A character vector of dates in the format "YYYY-MM-DDTHH:MM:SS" (e.g., "1961-01-01T00:00:00") representing the time series for which bootstrap intervals should be created. Same as is obtained with the `dimnames(ds[["time"]])` of a `CFDataset` object called `ds`.
+#' @param bootstrapWindow Numeric, uneven integer. Length of bootstrapping interval around each target date.
+#'
+#' @return A list of character vectors, each containing the dates that fall within the bootstrap window around each target date. The names of the list correspond to the day of year (in "MM-DD" format) for which the bootstrap intervals are created. Interval creation assumes a leap year (i.e., 366 days) to ensure that the same day of year is always included in the bootstrap intervals, even for February 29th.
+#'
+#' @author Erik Kusch
+#'
+#' @examples
+#' time_series <- seq(
+#'     as.POSIXct("1961-01-01 00:00:00", tz = "UTC"),
+#'     as.POSIXct("1999-12-31 23:00:00", tz = "UTC"),
+#'     by = "1 day"
+#' )
+#' Time_seq <- format(time_series, "%Y-%m-%dT%H:%M:%S")
+#' Helper_ETCCDIDailyBootstraps(dates = Time_seq)
+#'
+Helper_ETCCDIDailyBootstraps <- function(dates, bootstrapWindow = 5) {
+    ## input checking
+    if (bootstrapWindow %% 2 != 1) {
+        stop("Bootstrap window must be an odd number.")
+    }
+
+    ## transforming time sequence into dates
+    Date_seq <- format(as.POSIXct(dates, tryFormats = c("%Y-%m-%dT%H:%M:%OS"), tz = "UTC"), "%Y-%m-%d")
+    dates <- as.Date(Date_seq) # Convert Date_seq to Date objects
+    years <- sort(unique(as.integer(format(dates, "%Y")))) # Get unique years
+
+    ## actual grouping creation
+    half_window <- (bootstrapWindow - 1) / 2
+    offsets <- c(-half_window:-1, 0, 1:half_window)
+    leap_year_dates <- as.Date("2000-01-01") + (0:365) # Compute MM-DD names for each day of year using a leap year
+
+    # For each day of year
+    grouping_ls <- lapply(1:366, function(doy) {
+        centre_dates <- as.Date(paste0(years, "-01-01")) + (doy - 1)
+        target_dates <- do.call(c, lapply(centre_dates, function(centre_date) {
+            neigh_dates <- centre_date + offsets
+            neigh_str <- format(neigh_dates, "%Y-%m-%d")
+            neigh_str
+        }))
+        target_dates[target_dates %in% Date_seq] # Keep only dates that exist in the data
+    })
+    names(grouping_ls) <- format(leap_year_dates, "%m-%d")
+    grouping_ls
+}
diff --git a/R/Helper_ETCCDIGSL.r b/R/Helper_ETCCDIGSL.r
new file mode 100644
index 0000000..4737b78
--- /dev/null
+++ b/R/Helper_ETCCDIGSL.r
@@ -0,0 +1,187 @@
+#' @title Calculate Growing Season Length (GSL)
+#'
+#' @description Calculates Growing Season Length as defined by ETCCDI with different treatment of northern and southern hemisphere. For the northern hemisphere, GSL is defined as the number of days between the first occurrence of 6 consecutive days with mean temperature above 5°C after Jan 1st and the first occurrence of 6 consecutive days with mean temperature below 5°C after Jul 1st. For the southern hemisphere, GSL is defined as the number of days between the first occurrence of 6 consecutive days with mean temperature above 5°C after Jul 1st and the first occurrence of 6 consecutive days with mean temperature below 5°C after Jan 1st of the following year. If no valid start or end date is found, GSL is set to NA for that year and grid cell.
+#'
+#' @param TM A CFVariable containing daily mean temperatures (or similar) to be used for GSL calculation.
+#'
+#' @importFrom ncdfCF as_CF
+#'
+#' @return A CFVariable.
+#'
+#' @author Erik Kusch
+#'
+#' @examples
+#' \dontrun{
+#' TX <- NC_Read("inst/extdata/KiN_tasmax_2050.nc")[["tasmax"]]
+#' TN <- NC_Read("inst/extdata/KiN_tasmin_2050.nc")[["tasmin"]]
+#' TM <- (TN + TX) / 2
+#' GSL <- Helper_ETCCDIGSL(TM)
+#' }
+Helper_ETCCDIGSL <- function(TM) {
+    threshold <- 273.15 + 5
+    TM_array <- TM$raw()
+    dn <- dimnames(TM_array)
+    if (is.null(dn) || is.null(dn$lat)) {
+        stop("Array must have dimnames with a 'lat' entry.")
+    }
+    lat <- as.numeric(dn$lat)
+    if (anyNA(lat)) {
+        stop("Latitude dimnames could not be converted to numeric.")
+    }
+
+    north_idx <- which(lat >= 0)
+    south_idx <- which(lat < 0)
+    GSL_ls <- list(
+        north = TM_array[, north_idx, , drop = FALSE],
+        south = TM_array[, south_idx, , drop = FALSE]
+    )
+
+    first_run_start <- function(x, idx, conditionFun, runLength = 6) {
+        if (length(idx) < runLength) {
+            return(NA_integer_)
+        }
+        v <- conditionFun(x[idx])
+        v[is.na(v)] <- FALSE
+        r <- rle(v)
+
+        hit <- which(r$values & r$lengths >= runLength)
+        if (length(hit) == 0) {
+            return(NA_integer_)
+        }
+
+        relStart <- if (hit[1] == 1) 1L else (sum(r$lengths[seq_len(hit[1] - 1)]) + 1L)
+        idx[relStart]
+    }
+
+
+    GSL_calcs <- lapply(names(GSL_ls), function(hemi) {
+        arr <- GSL_ls[[hemi]]
+        # Skip if no latitude cells
+        if (is.null(arr) || length(dim(arr)) != 3 || dim(arr)[2] == 0) {
+            return(NULL)
+        }
+
+        ## figure out dates and years
+        dn <- dimnames(arr)
+        if (is.null(dn) || is.null(dn$time)) stop("Time dimnames are required.")
+        dates <- as.Date(dn$time)
+        yrs <- sort(unique(as.integer(format(dates, "%Y"))))
+
+        ## makke arrays for filling with information about start/end dates and GSL days
+        nlon <- dim(arr)[which(names(dn) == "lon")]
+        nlat <- dim(arr)[which(names(dn) == "lat")]
+        ny <- length(yrs)
+        startDate <- array(as.Date(NA),
+            dim = c(nlon, nlat, ny),
+            dimnames = list(lon = dn$lon, lat = dn$lat, time = as.character(yrs))
+        )
+        endDate <- array(as.Date(NA),
+            dim = c(nlon, nlat, ny),
+            dimnames = list(lon = dn$lon, lat = dn$lat, time = as.character(yrs))
+        )
+        gslDays <- array(NA_real_,
+            dim = c(nlon, nlat, ny),
+            dimnames = list(lon = dn$lon, lat = dn$lat, time = as.character(yrs))
+        )
+
+        ## iterate over years and grid cells to find start/end dates and GSL days
+        for (yIdx in seq_along(yrs)) {
+            y <- yrs[yIdx]
+
+            ## define windows in which to look for runs
+            if (hemi == "north") {
+                # Start window (north): Jan 1 -> Dec 31 of year y
+                startWin <- which(dates >= as.Date(sprintf("%d-01-01", y)) &
+                    dates <= as.Date(sprintf("%d-12-31", y)))
+
+                # End window (north): Jul 1 -> Dec 31 of year y
+                endWin <- which(dates >= as.Date(sprintf("%d-07-01", y)) &
+                    dates <= as.Date(sprintf("%d-12-31", y)))
+
+                # Fallback if no 6-day <= threshold run
+                endFallback <- as.Date(sprintf("%d-12-31", y))
+            } else {
+                # Start window (south): Jul 1 -> Dec 31 of year y
+                startWin <- which(dates >= as.Date(sprintf("%d-07-01", y)) &
+                    dates <= as.Date(sprintf("%d-12-31", y)))
+
+                # End window (south): Jan 1 -> Jun 30 of year y+1
+                # (June 31 is not a valid date, so Jun 30 is used)
+                yNext <- y + 1L
+                hasNextYear <- any(as.integer(format(dates, "%Y")) == yNext)
+
+                if (hasNextYear) {
+                    endWin <- which(dates >= as.Date(sprintf("%d-01-01", yNext)) &
+                        dates <= as.Date(sprintf("%d-06-30", yNext)))
+                    endFallback <- as.Date(sprintf("%d-06-30", yNext))
+                } else {
+                    endWin <- integer(0)
+                    endFallback <- as.Date(NA)
+                }
+            }
+
+            ## iterate over grid cells to find start/end dates and GSL days
+            for (i in seq_len(nlon)) {
+                for (j in seq_len(nlat)) {
+                    ts <- arr[i, j, ]
+
+                    sIdx <- first_run_start(ts, startWin, function(v) v > threshold)
+
+                    if (is.na(sIdx)) {
+                        # no valid season start -> leave NA
+                        next
+                    }
+
+                    # Find end run; if not found, use fallback date (except missing next year case in south)
+                    eIdx <- first_run_start(ts, endWin, function(v) v < threshold)
+
+                    sDate <- dates[sIdx]
+                    if (is.na(eIdx)) {
+                        eDate <- endFallback
+                    } else {
+                        eDate <- dates[eIdx]
+                    }
+
+                    startDate[i, j, yIdx] <- sDate
+                    endDate[i, j, yIdx] <- eDate
+
+                    if (!is.na(eDate)) {
+                        gslDays[i, j, yIdx] <- as.numeric(eDate - sDate)
+                        # If you want inclusive count, use: as.numeric(eDate - sDate) + 1
+                    }
+                }
+            }
+        }
+
+        # list(
+        # startDate = startDate,
+        # endDate = endDate,
+        # gslDays =
+        gslDays
+        # ,
+        # years = yrs
+        # )
+    })
+    names(GSL_calcs) <- names(GSL_ls)
+
+    ## fusing north and south-hemisphere arrays
+    gsl_array <- abind::abind(GSL_calcs$south, GSL_calcs$north, along = 2)
+
+    ## transforming back to CFVariable
+    names(dimnames(gsl_array)) <- names(dimnames(TM_array))
+    dimnames(gsl_array)[1:2] <- dimnames(TM_array)[1:2]
+    dimnames(gsl_array)$time <- paste0(dimnames(gsl_array)$time, "-07-02T12:00:00")
+    gsl_array_tll <- aperm(gsl_array, c(3, 1, 2))
+
+    # optional: fix dimname order explicitly
+    dn <- dimnames(gsl_array)
+    dimnames(gsl_array_tll) <- list(
+        time = dn$time,
+        lon  = dn$lon,
+        lat  = dn$lat
+    )
+    GSL <- as_CF("GSL", gsl_array_tll)
+
+    ## report back
+    return(GSL)
+}
diff --git a/R/Helper_Threshold.r b/R/Helper_Threshold.r
index b1123e7..59c1456 100644
--- a/R/Helper_Threshold.r
+++ b/R/Helper_Threshold.r
@@ -1,56 +1,100 @@
-#' @title Apply tresholding to raster layers
+#' @title Apply tresholding to CFVariable contents
 #'
-#' @description Loops over layers in a SpatRaster object and evaluates threshold criteria defined by user.
+#' @description Evaluates threshold criteria defined by user against a CFVariable and reports summaries through the `$summarise` method of `ncdfCF`.
 #'
-#' @param spatRaster A SpatRaster within which thresholding should be applied
-#' @param operator Character. One of ">", "<", ">=", "<=", "==", "!=" to evaluate spatRaster cells against threshold value.
-#' @param threshold Numeric or Character. threshold to evaluate spatRaster cells against using operator value.
-#' @param returnValues Logical. Whether to return spatRaster of values matching thresholding (TRUE) or only logical raster indicating where thresholding is matched (FALSE).
-#' @param verbose Logical. If progress should be displayed in the console.
+#' @param CFVariable A CFVariable within which thresholding should be applied
+#' @param operator Character. One of ">", "<", ">=", "<=", "==", "!=" to evaluate CFVariable cells against threshold value.
+#' @param threshold Numeric or Character. threshold to evaluate CFVariable cells against using operator value.
+#' @param threshMode Character, optional. If supplying a CFVariable containing baseline quantiles for ETCCDI calculation, setting `threshMode = "ETCCDIQuantiles"` will internally apply location and day-of-year specific thresholds to `CFVariable` input. Default is `NULL` avoiding this behaviour.
+#' @param returnValues Logical. Whether to return CFVariable of values matching thresholding (TRUE) or only logical raster indicating where thresholding is matched (FALSE).
+#' @param returnSummary Function. Summary function for tresholded CFVariable.
+#' @param returnTResolution Character. Temporal resolution of summary of tresholded CFVariable.
 #'
-#' @importFrom terra nlyr
+#' @importFrom ncdfCF as_CF
+#'
+#' @return A CFVariable.
 #'
-#' @return A SpatRaster.
-#' 
 #' @author Erik Kusch
 #'
 #' @examples
-#' Data_rast <- terra::rast(system.file("extdata", "KiN_rast.nc", package = "ClimHub"))[[1:31]]
+#' Data_CF <- NC_Read("inst/extdata/KiN_tas.nc")[["tas"]]
 #' Helper_Threshold(
-#'     spatRaster = Data_rast,
+#'     CFVariable = Data_CF,
 #'     operator = "<",
-#'     threshold = 2650
+#'     threshold = 273.15
 #' )
 #' Helper_Threshold(
-#'     spatRaster = Data_rast,
+#'     CFVariable = Data_CF,
 #'     operator = "<",
-#'     threshold = 2650,
-#'     returnValues = TRUE
+#'     threshold = 273.15,
+#'     returnValues = TRUE,
+#'     returnSummary = mean
 #' )
-Helper_Threshold <- function(spatRaster, operator, threshold, returnValues = FALSE, verbose = TRUE) {
+Helper_Threshold <- function(CFVariable, operator, threshold, threshMode = NULL, returnValues = FALSE, returnSummary, returnTResolution = "year") {
     ## input check, needs rwriting to Helper.InputChecker
     if (!operator %in% c(">", "<", ">=", "<=", "==", "!=")) {
         stop("operator must be one of: '>', '<', '>=', '<=', '==', '!='")
     }
 
-    ## progress bar
-    pb <- Helper_Progress(iterLength = nlyr(spatRaster), text = "Applying Thresholding")
-
-    ## going layer by layer to avoid error reported in https://github.com/rspatial/terra/issues/1556
-    ThresholdedLayers_ls <- lapply(1:nlyr(spatRaster), FUN = function(Iter) {
-        expr <- paste0("spatRaster[[Iter]] ", operator, " ", threshold)
-        ret_rast <- eval(parse(text = expr))
-        if (verbose) {
-            pb$tick(tokens = list(layer = Iter))
-        }
-        if (returnValues) {
-            ret_rast2 <- spatRaster[[Iter]]
-            ret_rast2[!ret_rast] <- NA
-            ret_rast <- ret_rast2
-        }
-        ret_rast
-    })
-
-    ## return fused layers as SpatRaster
-    do.call(c, ThresholdedLayers_ls)
+    ## Get the raw data array
+    raw_array <- CFVariable$raw()
+    array_dims <- dimnames(raw_array)
+
+    ## reformatting treshold to array if necessary
+    if (class(threshold)[1] == "CFVariable") {
+        threshold <- threshold$raw()
+    }
+
+    ## reformatting of threshold if ETCCDI Baseline is specified
+    if (identical(threshMode, "ETCCDIQuantiles")) {
+        # Extract the time labels from dimnames
+        raw_dates <- as.Date(sub("T.*$", "", dimnames(raw_array)$time)) # drops the "T12:00:00"
+        thresh_dates <- as.Date(dimnames(threshold)$time)
+        # Create a month-day key for matching
+        raw_key <- format(raw_dates, "%m-%d")
+        thresh_key <- format(thresh_dates, "%m-%d")
+        # match dates
+        idx <- match(raw_key, thresh_key)
+        # Now expand `threshold` to a 3D array matching raw_array time length:
+        threshold <- threshold[, , idx, drop = FALSE]
+        dimnames(threshold) <- array_dims
+    }
+
+    ## apply thresholding
+    if (class(threshold) == "array") { # this is comparing to arrays and setting every exceeding value in the former to NA
+        expr <- paste("raw_array", operator, "threshold")
+    } else { # this is setting every value that does NOT match the operator threshold to NA
+        expr <- paste0("raw_array", operator, threshold)
+    }
+    mask <- eval(parse(text = expr))
+    raw_array[!mask] <- NaN
+
+    ## binarise if return of values is not desired
+    if (!returnValues) {
+        raw_array[!is.na(raw_array)] <- 1 # these are those not masked due to threshold
+        raw_array[is.nan(raw_array)] <- 0 # these are those masked due to threshold
+    }
+    dimnames(raw_array) <- array_dims
+
+    # Create a new CFVariable from the filtered data
+    Thresh_CF <- ncdfCF::as_CF("Thresholded", raw_array)
+
+    # The original attributes - drop "actual_range" because it is no longer accurate
+    atts <- CFVariable$attributes[CFVariable$attributes$name != "actual_range", ]
+
+    # Loop over the attributes and set them in the new CFVariable
+    apply(atts, 1, function(a) Thresh_CF$set_attribute(a$name, a$type, a$value))
+
+    ## return summary only when explicitly requested by caller
+    ## (run-based ETCCDI logic needs unsummarised daily output)
+    if (!missing(returnSummary)) {
+        Thresh_CF <- Thresh_CF$summarise(
+            "Tresholded",
+            returnSummary,
+            returnTResolution
+        )
+    }
+
+    ## return object
+    Thresh_CF
 }
diff --git a/R/Metrics_BootstrapQuantiles.r b/R/Metrics_BootstrapQuantiles.r
new file mode 100644
index 0000000..5be402d
--- /dev/null
+++ b/R/Metrics_BootstrapQuantiles.r
@@ -0,0 +1,129 @@
+#' @title Bootstrap Quantile Lookup for Daily Time Series
+#'
+#' @description Computes quantiles across years for each day‑of‑year using a rolling bootstrap window to be used in `Metrics_ETTCDI` calculation of ETCCDI. Internally the function uses \code{Helper_ETCCDIDailyBootstraps()} to determine which dates belong in the window centred on each calendar day (assuming a leap year for consistency).
+#'
+#' This is primarily used by the ETCCDI metric code to obtain base period thresholds but can be called directly by users. Dates in the input `CFVariable` must present individual days.
+#'
+#' @param CFVariable `CFVariable` object for which quantiles for each bootstrap day interval ought to be calculated.
+#' @param probs Numeric vector of probabilities (between 0 and 1) for which quantiles will be returned.  Default is \code{c(0.1,0.9)}.
+#' @param bootstrapWindow Odd integer giving the width of the window in days used for the bootstrap sample around each target date. The window is centred on the day of year of interest.  Defaults to 5.
+#'
+#' @importFrom ncdfCF as_CF
+#' @importFrom ncdfCF create_ncdf
+#'
+#' @return A numeric array with dimensions \code{lon × lat × dayOfYear × probs}.  The third dimension has 366 levels ("01‑01" through "12‑31", including Feb 29) and the fourth dimension is named by the supplied \code{probs}.  Values are the requested quantiles computed over the set of years that fall inside the bootstrap window.  If no valid observations exist for a given cell/day the result is NA.
+#'
+#' @author Erik Kusch
+#'
+#' @examples
+#' \dontrun{
+#' TX <- NC_Read("inst/extdata/KiN_tasmax.nc")[["tasmax"]]
+#' TX_Quant <- Metrics_BootstrapQuantiles(CFVariable = TX, probs = c(0.1, 0.9), bootstrapWindow = 5)
+#' TX_Quant
+#' }
+#' @export
+Metrics_BootstrapQuantiles <- function(
+    CFVariable,
+    probs = c(0.1, 0.9),
+    bootstrapWindow = 5) {
+    ## transform input for internal use
+    data_array <- CFVariable$raw()
+    dates <- as.Date(dimnames(data_array)[[3]])
+
+    ## input checks
+    if (!is.numeric(data_array) || length(dim(data_array)) != 3) {
+        stop("data_array must be a numeric array with three dimensions")
+    }
+    if (any(probs < 0 | probs > 1)) {
+        stop("probs must be between 0 and 1")
+    }
+
+    ## extract years from the dates for later use (leap-year detection)
+    years <- sort(unique(as.integer(format(dates, "%Y"))))
+
+    ## create grouping information
+    group_ls <- Helper_ETCCDIDailyBootstraps(dates, bootstrapWindow = bootstrapWindow)
+    nday <- length(group_ls) # 366
+    nl <- dim(data_array)[1]
+    nt <- dim(data_array)[2]
+    nprob <- length(probs)
+
+    ## construct a time dimension that reflects a leap year; use the first leap year found in the input dates if present, otherwise fall back to 2000
+    leap_years <- years[years %% 4 == 0 & (years %% 100 != 0 | years %% 400 == 0)]
+    if (length(leap_years) > 0) {
+        base_year <- leap_years[1]
+    } else {
+        base_year <- 2000
+    }
+    time_dim <- paste0(base_year, "-", names(group_ls), "T00:00:00")
+
+    ## create empty array which is to be filled (third dim now actual timestamps)
+    out <- array(NA_real_,
+        dim = c(nl, nt, nday, nprob),
+        dimnames = list(
+            dimnames(data_array)[[1]],
+            dimnames(data_array)[[2]],
+            time_dim,
+            as.character(probs)
+        )
+    )
+
+    ## progress bar
+    pb <- Helper_Progress(iterLength = nday, text = "Calculating Quantiles")
+
+    ## iterate over the 366 days
+    for (doy in seq_len(nday)) {
+        ## select data for the relevant days in the base period
+        idx <- match(group_ls[[doy]], dates)
+        idx <- idx[!is.na(idx)]
+        if (length(idx) == 0) next
+        slice <- data_array[, , idx, drop = FALSE] # lon × lat × nobs
+
+        ## iterate over locations and calculate quantiles
+        qmat <- apply(slice, 1:2, function(v) {
+            if (all(is.na(v))) {
+                rep(NA_real_, nprob)
+            } else {
+                as.numeric(stats::quantile(v,
+                    probs = probs,
+                    na.rm = TRUE, type = 7
+                ))
+            }
+        })
+        dim(qmat) <- c(nprob, nl, nt)
+
+        ## assign quantiles to pre-defined array
+        for (ip in seq_len(nprob)) {
+            out[, , doy, ip] <- qmat[ip, , ]
+        }
+
+        ## update progress bar
+        pb$tick(tokens = list(layer = doy))
+    }
+
+    ## convert each probability slice into its own CFVariable and collect in a CFDataset
+    # create CFDataset for reporting back
+    ds <- ncdfCF::create_ncdf()
+
+    # attributes to copy (exclude actual_range since values have changed)
+    atts <- CFVariable$attributes[CFVariable$attributes$name != "actual_range", ]
+
+    for (ip in seq_len(nprob)) {
+        prob <- probs[ip]
+        slice <- out[, , , ip, drop = TRUE] # lon × lat × time
+        names(dimnames(slice)) <- names(dimnames(data_array))
+
+        # variable name uses probability so user can see it later
+        varname <- paste0("quantile_", gsub("\\.", "p", as.character(prob)))
+        new_var <- ncdfCF::as_CF(varname, slice)
+
+        # copy attributes
+        apply(atts, 1, function(a) new_var$set_attribute(a$name, a$type, a$value))
+
+        ## add data
+        ds$add_variable(new_var)
+    }
+
+    ## return CFDataset
+    ds
+}
diff --git a/R/Metrics_ETCCDI.r b/R/Metrics_ETCCDI.r
index 0cff643..e8be7e1 100644
--- a/R/Metrics_ETCCDI.r
+++ b/R/Metrics_ETCCDI.r
@@ -1,232 +1,764 @@
 #' @title Calculate ETCCDI Indices from List Input
-#' 
-#' @description This function calculates \href{https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/jgrd.50203}{ETCCDIs} from a named list of SpatRaster objects. Currently, only the first four ETCCDI are supported.
 #'
-#' @param ls Optional, list. List of SpatRasters. Names of elements must be "TX" and "TN", holding maximum and minimum daily air temperature respectively.
+#' @description This function calculates \href{https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/jgrd.50203}{ETCCDIs} from a named list of CFVariable objects. Currently, only some ETCCDI are supported.
 #'
-#' @importFrom terra units
-#' @importFrom terra varnames
-#' @importFrom terra longnames
-#' @importFrom terra time
+#' @details Input requirements are conditional on the selected ETCCDI indices in the `indices` argument (see table below). To reduce memory usage and speed up calculations, supply only the inputs required for the selected indices.
 #'
-#' @return A named list of SpatRasters with each element corresponding to an ETCCDI.
+#' | **ETCCDI** | **Name** | **Required Data** | **Description** | **Unit** | **Default `TResolution`** |
+#' |---|---|---|---|---|---|
+#' | FD | Number of frost days | `TN` | Count of days when TN (daily minimum temperature) < 0°C. | days | year |
+#' | SU | Number of summer days | `TX` | Count of days when TX (daily maximum temperature) > 25°C. | days | year |
+#' | ID | Number of icing days | `TX` | Count of days when TX (daily maximum temperature) < 0°C. | days | year |
+#' | TR | Number of tropical nights | `TN` | Count of days when TN (daily minimum temperature) > 20°C. | days | year |
+#' | GSL | Growing Season Length | `TX`, `TN` | Annual count of days between the first occurrence of at least 6 consecutive days with (TN+TX)/2 > 5°C and the first occurrence after 1st July (Northern Hemisphere) or 1st January (Southern Hemisphere) of at least 6 consecutive days with (TN+TX)/2 < 5°C | days | year |
+#' | TXx | Max of Daily Max Temp | `TX` | Maximum daily maximum temperature. | K | month |
+#' | TNx | Max of Daily Min Temp | `TN` | Maximum daily minimum temperature. | K | month |
+#' | TXn | Min of Daily Max Temp | `TX` | Minimum daily maximum temperature. | K | month |
+#' | TNn | Min of Daily Min Temp | `TN` | Minimum daily minimum temperature. | K | month |
+#' | TN10p | Percent Days TN < 10th Percentile | `TN`, `TN_Base` | Percent of days where TN < 10th percentile of base period. | % | year |
+#' | TX10p | Percent Days TX < 10th Percentile | `TX`, `TX_Base` | Percent of days where TX < 10th percentile of base period. | % | year |
+#' | TN90p | Percent Days TN > 90th Percentile | `TN`, `TN_Base` | Percent of days where TN > 90th percentile of base period. | % | year |
+#' | TX90p | Percent Days TX > 90th Percentile | `TX`, `TX_Base` | Percent of days where TX > 90th percentile of base period. | % | year |
+#' | WSDI | Warm Spell Duration Index | `TX`, `TX_Base` | Count of days contained within runs of 6+ consecutive days when TX > 90th percentile of base period. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. | days | year |
+#' | CSDI | Cold Spell Duration Index | `TN`, `TN_Base` | Count of days contained within runs of 6+ consecutive days when TN < 10th percentile of base period. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. | days | year |
+#' | DTR | Daily Temperature Range | `TX`, `TN` | Mean difference between daily max (TX) and min (TN) temperatures. | K | month |
+#' | Rx1day | Max 1-day Precipitation | `RR` | Maximum precipitation in a single day. | mm | month |
+#' | Rx5day | Max 5-day Precipitation | `RR` | Maximum precipitation over any 5 consecutive days. | mm | month |
+#' | SDII | Simple Precipitation Intensity Index | `RR` | Mean precipitation amount on wet days (RR ≥ 1mm). | mm/day | year |
+#' | R10mm | Days with Precip ≥ 10mm | `RR` | Count of days with precipitation ≥ 10mm. | days | year |
+#' | R20mm | Days with Precip ≥ 20mm | `RR` | Count of days with precipitation ≥ 20mm. | days | year |
+#' | Rnnmm | Days with Precip ≥ user-defined threshold | `RR` | Count of days with precipitation ≥ nnmm. | days | year |
+#' | CDD | Consecutive Dry Days | `RR` | Maximum number of consecutive days with RR < 1mm. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. | days | year |
+#' | CWD | Consecutive Wet Days | `RR` | Maximum number of consecutive days with RR ≥ 1mm. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. | days | year |
+#' | R95pTOT | Total Precipitation from RR > 95th Percentile | `RR`, `RR_Base` | Total precipitation from wet days (RR > 95th percentile of base period). | mm | year |
+#' | R99pTOT | Total Precipitation from RR > 99th Percentile | `RR`, `RR_Base` | Total precipitation from wet days (RR > 99th percentile of base period). | mm | year |
+#' | RCPTOT | Total Precipitation on Wet Days | `RR` | Sum of precipitation on wet days (RR ≥ 1mm). | mm | year |
+#'
+#' @param projectionList List. List of `CFVariable` objects required by the selected ETCCDI indices via the `indices` argument. Include only named elements that are needed for `indices`: "TX" (daily maximum air temperature in Kelvin), "TN" (daily minimum air temperature in Kelvin), and/or "RR" (daily total precipitation in mm). See details for required data input per index.
+#' @param baseLineList Optional, list. List of `CFDataset` objects containing baseline quantiles required only for quantile-based ETCCDI indices. Include only named elements needed for `indices`: "TX_Base", "TN_Base", and/or "RR_Base" (note that these must be in Kelvin, Kelvin, and mm, respectively). If no quantile-based ETCCDI is selected, this argument can be omitted. See details for required quantile baselines per index.
+#' @param indices Optional, character. Character vector of ETCCDI abbreviations to calculate (see first column of the details table). If missing, all supported indices in this function are calculated. See details for supported indices.
+#' @param TResolution Optional, character. Temporal resolution for ETCCDI calculation. Supports "year", "month" and "season". If omitted, each selected index is calculated using its default temporal resolution from the details table. If provided, the same resolution is used for all selected indices and output indices whose default resolution differs are prefixed with "ALT_". Note that this is not applied to GSL which is always calculated at annual resolution due to its specific logic.
+#' @param RRThreshold Numeric. Custom threshold for daily precipiation in mm for calculation of Rnnmm. Defaults to 42.
+#' @param fileName Character, optional. Character. A file name for the produced file, including path and ".nc" file ending. If no value is supplied, the dataset is not written to disk but returned as a `CFDataset` object in memory. If a file name is supplied and a file with that name already exists, the function will attempt to load and return that file instead of recalculating the indices.
+#'
+#' @importFrom ncdfCF as_CF
+#' @importFrom ncdfCF create_ncdf
+#'
+#' @return A `CFDataset`. If all calculated ETCCDI indices have the same temporal resolution (e.g., when `TResolution` is specified or a subset of `indices` is selected that share a temporal resolution default), the returned `CFDataset` contains all variables in the root group. If indices span multiple temporal resolutions (e.g., when `TResolution` is omitted and some selected `indices` default to "year" while others default to "month"), the returned `CFDataset` contains one subgroup per temporal resolution, each containing only variables of that temporal resolution.
+#'
+#' Each variable is named by its ETCCDI acronym and has a `long_name` attribute describing the index. Acronyms of indices calculated with a non-default temporal resolution are prefixed with "ALT_" to indicate that these are alternative calculations.
 #'
 #' @author Erik Kusch
 #'
 #' @examples
 #' \dontrun{
-#' TX_rast <- terra::rast("inst/extdata/Jotunheimen_TX.nc") # terra::rast(system.file("extdata", "Jotunheimen_TX.nc", package = "ClimHub"))
-#' TN_rast <- terra::rast("inst/extdata/Jotunheimen_TN.nc")
-#' RR_rast <- terra::rast("inst/extdata/Jotunheimen_RR.nc")
-#' terra::units(RR_rast) <- "mm"
-#' # BASEPeriod_rast <- terra::rast("inst/extdata/Jotunheimen_BASEPeriod.nc")
-#'
-#' Metrics_ETCCDI(ls = list(TX = TX_rast, TN = TN_rast, RR = RR_rast
-#' # , BASEPeriod = BASEPeriod_rast)
+#' ## Directory for Data and Output -----------
+#' Dir.Data <- file.path(getwd(), "ExampleData")
+#' dir.create(Dir.Data, showWarnings = FALSE)
+#'
+#' ## Downloading and preparing data -----------
+#' ### Variables for which we need data
+#' vars <- c("minimum_air_temperature", "maximum_air_temperature", "precipitation_flux")
+#' names(vars) <- c("tasmin", "tasmax", "pr")
+#'
+#' ### Data (Down-)Loading and Preparation
+#' Data_ls <- lapply(seq_along(vars), FUN = function(i) {
+#'     print(as.character(vars[i]))
+#'     ## these need to be made into quantiles to define base-period for ETCCDI calculation
+#'     KiN_Base <- Access_KlimaiNorge2100(
+#'         variable = as.character(vars[i]),
+#'         dateStart = "1971-01-01", # note that the standard base period starts 1961, but KiN does start in 1971 so we use that as the start date
+#'         dateStop = "2000-12-31",
+#'         extent = c(9, 11, 59, 61),
+#'         model = "noresm-r1i1p1f1-hclim",
+#'         method = "eqm",
+#'         scenario = "ssp370",
+#'         fileName = file.path(Dir.Data, paste0("KiN_", names(vars)[i], ".nc"))
+#'     )
+#'
+#'     ## these are the files for calculation of ETCCDI metrics in future
+#'     KiN_2090 <- Access_KlimaiNorge2100(
+#'         variable = as.character(vars[i]),
+#'         dateStart = "2090-01-01",
+#'         dateStop = "2099-12-31",
+#'         extent = c(9, 11, 59, 61),
+#'         model = "noresm-r1i1p1f1-hclim",
+#'         method = "eqm",
+#'         scenario = "ssp370",
+#'         fileName = file.path(Dir.Data, paste0("KiN_", names(vars)[i], "_2090.nc"))
+#'     )
+#'
+#'     ## these are the files for baseline quantiles
+#'     QuantF <- file.path(Dir.Data, paste0("KiN_", names(vars)[i], "_BaseLineQuantiles.nc"))
+#'     if (!file.exists(QuantF)) {
+#'         if (names(vars)[i] == "pr") {
+#'             probs_vec <- c(0.95, 0.99)
+#'             input <- KiN_Base[[names(vars)[i]]] * 86400 # to get from mm/day to kg m-2 s-1
+#'             input$set_attribute("units", "NC_CHAR", "mm")
+#'         } else {
+#'             probs_vec <- c(0.1, 0.9)
+#'             input <- KiN_Base[[names(vars)[i]]]
+#'         }
+#'         Quant <- Metrics_BootstrapQuantiles(
+#'             CFVariable = input,
+#'             probs = probs_vec,
+#'             bootstrapWindow = 5
+#'         )
+#'         Quant$save(QuantF)
+#'     } else {
+#'         Quant <- NC_Read(QuantF)
+#'     }
+#'
+#'     return(list(projection = KiN_2090, baseLine = Quant))
+#' })
+#' names(Data_ls) <- names(vars)
+#'
+#' ### Data Object listing
+#' TX <- Data_ls[["tasmax"]][["projection"]][["tasmax"]]
+#' TN <- Data_ls[["tasmin"]][["projection"]][["tasmin"]]
+#' RR <- Data_ls[["pr"]][["projection"]][["pr"]] * 86400 # to get from mm/day to kg m-2 s-1
+#' RR$set_attribute("units", "NC_CHAR", "mm")
+#'
+#' TX_Base <- Data_ls[["tasmax"]][["baseLine"]]
+#' TN_Base <- Data_ls[["tasmin"]][["baseLine"]]
+#' RR_Base <- Data_ls[["pr"]][["baseLine"]]
+#'
+#' rm(Data_ls) # remove list to free up memory
+#'
+#' ## ETCCDI calculation -----------
+#' ### All ETCCDI
+#' Metrics_ETCCDI(
+#'     projectionList = list(TX = TX, TN = TN, RR = RR),
+#'     baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+#'     fileName = file.path(Dir.Data, "ETCCDI_Metrics.nc")
 #' )
-#' }
 #'
+#' ### Select ETCCDI
+#' Metrics_ETCCDI(
+#'     projectionList = list(TX = TX, TN = TN, RR = RR),
+#'     baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+#'     indices = c("FD", "TXx", "R95pTOT")
+#' )
+#'
+#' ### Select ETCCDI with custom temporal resolutions
+#' Metrics_ETCCDI(
+#'     projectionList = list(RR = RR),
+#'     # baseLineList = list(RR_Base = RR_Base), # we do not need to supply the baseline quantiles if we do not select any indices that depend on them
+#'     indices = c("Rnnmm"),
+#'     RRThreshold = 21, # this is the threshold for the Rnnmm metric, in mm
+#'     TResolution = "season"
+#' )
+#'
+#' Metrics_ETCCDI(
+#'     projectionList = list(TX = TX, TN = TN, RR = RR),
+#'     baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+#'     indices = c("FD", "SU", "ID", "TR"),
+#'     TResolution = "month" # note that this is not the default for these indices, so the returned object will bear the prefix "ALT_" for each index to indicate that this is an alternative calculation with a different temporal resolution
+#' )
+#' }
 #' @export
-Metrics_ETCCDI <- function(ls) {
-    ## validate specification
-    InCheck_ls <- list(
-        Unit_TX = list(
-            Input = unique(units(ls$TX)),
-            Allowed = c("K"),
-            Operator = "in"
-        ),
-        Unit_TN = list(
-            Input = unique(units(ls$TN)),
-            Allowed = c("K"),
-            Operator = "in"
-        ),
-        Unit_RR = list(
-            Input = unique(units(ls$RR)),
-            Allowed = c("mm"),
-            Operator = "in"
-        )
+Metrics_ETCCDI <- function(projectionList, baseLineList, indices, TResolution, RRThreshold = 42, fileName) {
+    ## fileName handling
+    if (!missing(fileName)) {
+        fileName <- normalizePath(fileName, mustWork = FALSE)
+    }
+
+    ## index selection handling (moved before FCheck)
+    supportedIndices <- c(
+        "FD", "SU", "ID", "TR", "GSL", "TXx", "TNx", "TXn", "TNn", "TN10p",
+        "TX10p", "TN90p", "TX90p", "WSDI", "CSDI", "DTR", "Rx1day", "Rx5day",
+        "SDII", "R10mm", "R20mm", "Rnnmm", "CDD", "CWD", "R95pTOT", "R99pTOT", "RCPTOT"
+    )
+    if (missing(indices)) {
+        selectedIndices <- supportedIndices
+    } else {
+        if (!is.character(indices)) {
+            stop("'indices' must be a character vector of ETCCDI abbreviations.")
+        }
+        if (length(indices) == 0) {
+            stop("'indices' cannot be an empty vector.")
+        }
+        unknownIndices <- setdiff(indices, supportedIndices)
+        if (length(unknownIndices) > 0) {
+            stop(
+                paste0(
+                    "Unsupported ETCCDI abbreviation(s) in 'indices': ",
+                    paste(unknownIndices, collapse = ", "),
+                    ". Supported abbreviations are: ",
+                    paste(supportedIndices, collapse = ", "),
+                    "."
+                )
+            )
+        }
+        selectedIndices <- unique(indices)
+    }
+
+    ## ETCCDI default temporal resolution handling (moved before FCheck)
+    defaultTResolution <- c(
+        FD = "year",
+        SU = "year",
+        ID = "year",
+        TR = "year",
+        GSL = "year",
+        TXx = "month",
+        TNx = "month",
+        TXn = "month",
+        TNn = "month",
+        TN10p = "year",
+        TX10p = "year",
+        TN90p = "year",
+        TX90p = "year",
+        WSDI = "year",
+        CSDI = "year",
+        DTR = "month",
+        Rx1day = "month",
+        Rx5day = "month",
+        SDII = "year",
+        R10mm = "year",
+        R20mm = "year",
+        Rnnmm = "year",
+        CDD = "year",
+        CWD = "year",
+        R95pTOT = "year",
+        R99pTOT = "year",
+        RCPTOT = "year"
     )
-    Helper_InputChecker(inputCheck = InCheck_ls)
+    userProvidedTResolution <- !missing(TResolution)
+    if (userProvidedTResolution) {
+        if (!is.character(TResolution) || length(TResolution) != 1) {
+            stop("'TResolution' must be a single character value: 'year', 'month', or 'season'.")
+        }
+        providedTResolution <- tolower(TResolution)
+        if (is.na(providedTResolution) || !providedTResolution %in% c("year", "month", "season")) {
+            stop("'TResolution' must be one of 'year', 'month', or 'season'.")
+        }
+    } else {
+        providedTResolution <- NULL
+    }
+    get_index_t_resolution <- function(indexName) {
+        if (userProvidedTResolution) {
+            providedTResolution
+        } else {
+            defaultTResolution[[indexName]]
+        }
+    }
+
+    ## Create tResolution_ls and identify unique resolutions (for FCheck)
+    tResolution_ls <- list()
+    for (idx in selectedIndices) {
+        tResolution_ls[[idx]] <- get_index_t_resolution(idx)
+    }
+    unique_resolutions <- unique(unlist(tResolution_ls))
+
+    ## File Check with awareness of temporal resolution structure =========
+    if (!missing(fileName)) {
+        FCheck <- Helper_FileCheck(fileName = fileName, loadFun = NC_Read, load = TRUE, verbose = TRUE)
+        if (!is.null(FCheck)) {
+            return(FCheck)
+        }
+    }
+
+    ## Setting up Progress Bar
+    pb <- Helper_Progress(iterLength = length(selectedIndices), text = "ETCCDI Calculation")
+
+    ## summary functions
+    ### assign full run lengths to run end-points on the daily axis
+    ### this lets temporal aggregation assign a boundary-crossing run to the period where it ends
+    run_end_weighted_series <- function(binary_cf, minRunLength = 1) {
+        raw_array <- binary_cf$raw()
+        dn <- dimnames(raw_array)
+
+        nlon <- dim(raw_array)[which(names(dn) == "lon")]
+        nlat <- dim(raw_array)[which(names(dn) == "lat")]
+
+        # Allocate output array with the same shape as the input, filled with zeros. Only run end-point positions will be set to non-zero values.
+        out_array <- array(0, dim = dim(raw_array), dimnames = dn)
+
+        # Loop over every spatial grid cell, identical in style to Helper_ETCCDIGSL.
+        for (i in seq_len(nlon)) {
+            for (j in seq_len(nlat)) {
+                # Extract the full daily time series for this lon/lat cell.
+                x <- raw_array[i, j, ]
+
+                # If the entire time series is NA (e.g. ocean mask), propagate NAs and move on.
+                if (all(is.na(x))) {
+                    out_array[i, j, ] <- NA_real_
+                    next
+                }
+
+                # Treat NAs as non-events (0) so they break consecutive runs rather than silently bridging two separate runs into one.
+                x[is.na(x)] <- 0
+
+                # Run-length encode the full daily binary series (1 = event, 0 = non-event). rle() returns $lengths (number of consecutive days) and $values (TRUE/FALSE).
+                runs <- rle(x == 1)
+
+                # cumsum of lengths gives the time-axis index of the LAST day of each run. Example: lengths c(3, 4, 2) → end-points at positions 3, 7, 9. We use this as an index for where a run ends so we can use it for assignment of the run length value to the correct day in the output array.
+                run_ends <- cumsum(runs$lengths)
+
+                # Keep only truly event runs (values == TRUE) that are long enough to qualify.
+                # minRunLength = 1 keeps all runs (used by CDD and CWD).
+                # minRunLength = 6 keeps only spells of 6+ consecutive days (used by WSDI and CSDI).
+                qualifying_runs <- runs$values & runs$lengths >= minRunLength
+
+                if (any(qualifying_runs)) {
+                    # Write the full run length at the run's end-point index. All other days stay 0. This encodes the boundary-attribution rule: a run that crosses a period boundary (e.g. 4 days in year 1 + 3 days in year 2, total 7) will place its value (7) at the end-point day in year 2. When the temporal summary then sums (WSDI/CSDI) or takes the max (CDD/CWD) per period, year 1 contributes nothing and year 2 receives the full run length.
+                    out_array[i, j, run_ends[qualifying_runs]] <- runs$lengths[qualifying_runs]
+                }
+            }
+        }
+
+        # Wrap the result back into a CFVariable so it can be passed to $summarise().
+        run_cf <- ncdfCF::as_CF("Thresholded", out_array)
+        # Copy metadata attributes from the source, excluding actual_range which is now stale.
+        atts <- binary_cf$attributes[binary_cf$attributes$name != "actual_range", ]
+        if (nrow(atts) > 0) {
+            apply(atts, 1, function(a) run_cf$set_attribute(a$name, a$type, a$value))
+        }
+        run_cf
+    }
+
+    ### compute maximum of 5 day intervals
+    max_sum_over_5 <- function(x) {
+        if (length(na.omit(x)) == 0) {
+            return(NA_real_)
+        }
+        x[is.na(x)] <- 0 # keep the NA‑to‑0 convention if you like
+        n <- length(x)
+        if (n < 5) {
+            return(NA_real_)
+        }
+        best <- -Inf
+        for (i in seq_len(n - 5 + 1)) {
+            best <- max(best, sum(x[i:(i + 4)], na.rm = TRUE))
+        }
+        best
+    }
+
+    ### percentage of exceeding values
+    percentage_non_na <- function(x) {
+        x <- na.omit(x)
+        if (length(x) == 0) {
+            return(NA_real_)
+        } else {
+            return((sum(x) / length(x)) * 100)
+        }
+    }
+
+    ### regular functions ignoring NAs
+    sum_non_na <- function(x) {
+        x <- x[!is.na(x)]
+        if (length(x) == 0) {
+            NA_real_
+        } else {
+            sum(x, na.rm = TRUE)
+        }
+    }
+    mean_non_na <- function(x) {
+        x <- x[!is.na(x)]
+        if (length(x) == 0) {
+            NA_real_
+        } else {
+            mean(x, na.rm = TRUE)
+        }
+    }
+    max_non_na <- function(x) {
+        x <- x[!is.na(x)]
+        if (length(x) == 0) {
+            NA_real_
+        } else {
+            max(x, na.rm = TRUE)
+        }
+    }
+    min_non_na <- function(x) {
+        x <- x[!is.na(x)]
+        if (length(x) == 0) {
+            NA_real_
+        } else {
+            min(x, na.rm = TRUE)
+        }
+    }
+
+    ## get list contents
+    list2env(projectionList, env = environment())
+    if (!missing(baseLineList)) {
+        list2env(baseLineList, env = environment())
+    }
+
+    ## validate specification
+    # skipping this for the time being
+    # InCheck_ls <- list(
+    #     Unit_Projection_TX = list(
+    #         Input = unlist(TX$attributes$value[TX$attributes$name == "units"]),
+    #         Allowed = c("K"),
+    #         Operator = "in"
+    #     ),
+    #     Unit_Projection_TN = list(
+    #         Input = unlist(TN$attributes$value[TN$attributes$name == "units"]),
+    #         Allowed = c("K"),
+    #         Operator = "in"
+    #     ),
+    #     Unit_Projection_RR = list(
+    #         Input = unlist(RR$attributes$value[RR$attributes$name == "units"]),
+    #         Allowed = c("mm"),
+    #         Operator = "in"
+    #     )
+    # )
+    # Helper_InputChecker(inputCheck = InCheck_ls)
     # should also check here for:
     #  1. temporal resolution being days
     #  2. time ranges being neat years
     #  3. time across all inputs in Rasters being the same
     #  4. Check that CRS is EPS:4326
-
-    #stop("Add extra checks")
-    #stop("do quantile calculations for percentile based ETCCDI")
-
-    ## split into years
-    AnnualRasters_ls <- lapply(ls, Helper_TimeSplit, tResolution = "Year")
-    Dates <- as.POSIXct(paste0(unique(names(AnnualRasters_ls[[1]])), "-01-01"), tz = "UTC") ## merging lists into SpatRaster does not like numeric names of list elements
-
-    ## ETCCDI
-    ### Frost Days; Number of frost days: Annual count of days when TN (daily minimum temperature) < 0°C.
-    message("===== FD - Number of frost days =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$TN, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 273.15, operator = "<"))
-    })
-    names(Thresholded_ls) <- NULL
-    FD_rast <- do.call(c, Thresholded_ls)
-
-    ### Summer Days; Number of summer days: Annual count of days when TX (daily maximum temperature) > 25°C.
-    message("===== SU - Number of summer days =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$TX, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 273.15 + 25, operator = ">"))
-    })
-    names(Thresholded_ls) <- NULL
-    SU_rast <- do.call(c, Thresholded_ls)
-
-    ### Icing Days; Number of icing days: Annual count of days when TX (daily maximum temperature) < 0°C.
-    message("===== ID - Number of icing days =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$TX, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 273.15, operator = "<"))
-    })
-    names(Thresholded_ls) <- NULL
-    ID_rast <- do.call(c, Thresholded_ls)
-
-    ### Tropical Nights; Number of tropical nights: Annual count of days when TN (daily minimum temperature) > 20°C.
-    message("===== TR - Number of tropical nights =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$TN, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 273.15 + 20, operator = ">="))
-    })
-    names(Thresholded_ls) <- NULL
-    TR_rast <- do.call(c, Thresholded_ls)
-
-    ### GSL - Growing Season Length: Count the number of days between the first occurrence of at least 6 consecutive days with (TN+TX)/2 > 5°C and the first occurrence after 1st July (Northern Hemisphere) or 1st January (Southern Hemisphere) of at least 6 consecutive days with (TN+TX)/2 < 5°C
-    message("===== GSL - Growing Season Length =====")
-    message("Not implemented yet")
-
-    ### TXx - Monthly Max of Daily Max Temp: Maximum daily maximum temperature in each month.
-    message("===== TXx - Monthly Max of Daily Max Temp =====")
-    TXx_rast <- Temporal_Aggregation(spatRaster = ls$TX, tResolution = "month", tStep = 1, fun = max)
-
-    ### TNx - Monthly Max of Daily Min Temp: Maximum daily minimum temperature in each month.
-    message("===== TNx - Monthly Max of Daily Min Temp =====")
-    TNx_rast <- Temporal_Aggregation(spatRaster = ls$TN, tResolution = "month", tStep = 1, fun = max)
-
-    ### TXn - Monthly Min of Daily Max Temp: Minimum daily maximum temperature in each month.
-    message("===== TXn - Monthly Min of Daily Max Temp =====")
-    TXn_rast <- Temporal_Aggregation(spatRaster = ls$TX, tResolution = "month", tStep = 1, fun = min)
-
-    ### TNn - Monthly Min of Daily Min Temp: Minimum daily minimum temperature in each month.
-    message("===== TNn - Monthly Min of Daily Min Temp =====")
-    TNn_rast <- Temporal_Aggregation(spatRaster = ls$TN, tResolution = "month", tStep = 1, fun = min)
-
-    ### TN10p - Percent Days TN < 10th Percentile: Percent of days, per year, where TN < 10th percentile of base period.
-    message("===== TN10p - Percent Days TN < 10th Percentile =====")
-    message("Not implemented yet")
-
-    ### TX10p - Percent Days TX < 10th Percentile: Percent of days, per year, where TX < 10th percentile of base period.
-    message("===== TX10p - Percent Days TX < 10th Percentile =====")
-    message("Not implemented yet")
-
-    ### TN90p - Percent Days TN > 90th Percentile: Percent of days, per year, where TN > 90th percentile of base period.
-    message("===== TN90p - Percent Days TN > 90th Percentile =====")
-    message("Not implemented yet")
-
-    ### TX90p - Percent Days TX > 90th Percentile: Percent of days, per year, where TX > 90th percentile of base period.
-    message("===== TX90p - Percent Days TX > 90th Percentile =====")
-    message("Not implemented yet")
-
-    ### WSDI - Warm Spell Duration Index: Annual count of days with 6+ consecutive days when TX > 90th percentile of base period.
-    message("===== WSDI - Warm Spell Duration Index =====")
-    message("Not implemented yet")
-
-    ### CSDI - Cold Spell Duration Index: Annual count of days with 6+ consecutive days when TN < 10th percentile of base period.
-    message("===== CSDI - Cold Spell Duration Index =====")
-    message("Not implemented yet")
-
-    ### DTR - Daily Temperature Range: Monthly mean difference between daily max (TX) and min (TN) temperatures.
-    message("===== DTR - Daily Temperature Range =====")
-    DTR_rast <- Temporal_Aggregation(spatRaster = ls$TX - ls$TN, tResolution = "month", tStep = 1, fun = mean)
-
-    ### Rx1day - Max 1-day Precipitation per Month: Maximum precipitation in a single day each month.
-    message("===== Rx1day - Max 1-day Precipitation per Month =====")
-    Rx1day_rast <- Temporal_Aggregation(spatRaster = ls$RR, tResolution = "month", tStep = 1, fun = max)
-
-    ### Rx5day - Max 5-day Precipitation per Month: Maximum precipitation over any 5 consecutive days in each month.
-    message("===== Rx5day - Max 5-day Precipitation per Month =====")
-    message("Not implemented yet")
-
-    ### SDII - Simple Precipitation Intensity Index: Mean precipitation amount on wet days (RR ≥ 1mm).
-    message("===== SDII - Simple Precipitation Intensity Index =====")
-    message("Not yet tested")
-    SDII_rast <- ls$RR
-    SDII_rast[SDII_rast < 1] <- NA
-    SDII_rast <- Temporal_Aggregation(spatRaster = SDII_rast, tResolution = "month", tStep = 1, fun = mean)
-
-    ### R10mm - Days with Precip ≥ 10mm: Annual count of days with precipitation ≥ 10mm.
-    message("===== R10mm - Days with Precip ≥ 10mm =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$RR, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 10, operator = ">="))
-    })
-    names(Thresholded_ls) <- NULL
-    R10mm_rast <- do.call(c, Thresholded_ls)
-
-    ### R20mm - Days with Precip ≥ 20mm: Annual count of days with precipitation ≥ 20mm.
-    message("===== R20mm - Days with Precip ≥ 20mm =====")
-    Thresholded_ls <- lapply(AnnualRasters_ls$RR, FUN = function(RastIter) {
-        sum(Helper_Threshold(spatRaster = RastIter, threshold = 20, operator = ">="))
-    })
-    names(Thresholded_ls) <- NULL
-    R20mm_rast <- do.call(c, Thresholded_ls)
-
-    ### Rnnmm - Days with Precip ≥ user-defined threshold: Annual count of days with precipitation ≥ nnmm.
-    message("===== Rnnmm - Days with Precip ≥ user-defined threshold =====")
-    message("Not implemented yet")
-
-    ### CDD - Consecutive Dry Days: Maximum number of consecutive days with RR < 1mm.
-    message("===== CDD - Consecutive Dry Days =====")
-    message("Not implemented yet")
-
-    ### CWD - Consecutive Wet Days: Maximum number of consecutive days with RR ≥ 1mm.
-    message("===== CWD - Consecutive Wet Days =====")
-    message("Not implemented yet")
-
-    ### R95pTOT - Annual Precip from RR > 95th Percentile: Total precipitation from wet days (RR > 95th percentile of base period).
-    message("===== R95pTOT - Annual Precip from RR > 95th Percentile =====")
-    message("Not implemented yet")
-
-    ### R99pTOT - Annual Precip from RR > 99th Percentile: Total precipitation from wet days (RR > 99th percentile of base period).
-    message("===== R99pTOT - Annual Precip from RR > 99th Percentile =====")
-    message("Not implemented yet")
-
-    ### PRCPTOT - Annual Total Precipitation on Wet Days: Sum of precipitation on wet days (RR ≥ 1mm) over a year.
-    message("===== PRCPTOT - Annual Total Precipitation on Wet Days =====")
-    message("Not yet tested")
-    PRCPTOT_rast <- ls$RR
-    PRCPTOT_rast[PRCPTOT_rast < 1] <- 0
-    PRCPTOT_rast <- Temporal_Aggregation(spatRaster = PRCPTOT_rast, tResolution = "year", tStep = 1, fun = sum)
-
-    ## return
-    ## Fusing objects
-    Return_ls <- list(
-        FD = FD_rast,
-        SU = SU_rast,
-        ID = ID_rast,
-        TR = TR_rast,
-        R10 = R10mm_rast,
-        R20 = R20mm_rast
+    #  5. same as above for elements of baseLineList
+
+    # stop("Add extra checks")
+    # stop("do quantile calculations for percentile based ETCCDI")
+
+    ## ETCCDI calculators
+    indexCalculators <- list(
+        FD = function() {
+            Helper_Threshold(
+                CFVariable = TN,
+                operator = "<",
+                threshold = 273.15,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("FD")
+            )
+        },
+        SU = function() {
+            Helper_Threshold(
+                TX,
+                operator = ">",
+                threshold = 273.15 + 25,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("SU")
+            )
+        },
+        ID = function() {
+            Helper_Threshold(
+                TX,
+                operator = "<",
+                threshold = 273.15,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("ID")
+            )
+        },
+        TR = function() {
+            Helper_Threshold(
+                TN,
+                operator = ">",
+                threshold = 273.15 + 20,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("TR")
+            )
+        },
+        GSL = function() {
+            TM <- (TN + TX) / 2
+            Helper_ETCCDIGSL(TM)
+        },
+        TXx = function() {
+            TX$summarise("Tresholded", max_non_na, get_index_t_resolution("TXx"))
+        },
+        TNx = function() {
+            TN$summarise("Tresholded", max_non_na, get_index_t_resolution("TNx"))
+        },
+        TXn = function() {
+            TX$summarise("Tresholded", min_non_na, get_index_t_resolution("TXn"))
+        },
+        TNn = function() {
+            TN$summarise("Tresholded", min_non_na, get_index_t_resolution("TNn"))
+        },
+        TN10p = function() {
+            Helper_Threshold(
+                TN,
+                operator = "<",
+                threshold = TN_Base[[names(TN_Base)[1]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnSummary = percentage_non_na,
+                returnTResolution = get_index_t_resolution("TN10p")
+            )
+        },
+        TX10p = function() {
+            Helper_Threshold(
+                TX,
+                operator = "<",
+                threshold = TX_Base[[names(TX_Base)[1]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnSummary = percentage_non_na,
+                returnTResolution = get_index_t_resolution("TX10p")
+            )
+        },
+        TN90p = function() {
+            Helper_Threshold(
+                TN,
+                operator = ">",
+                threshold = TN_Base[[names(TN_Base)[2]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnSummary = percentage_non_na,
+                returnTResolution = get_index_t_resolution("TN90p")
+            )
+        },
+        TX90p = function() {
+            Helper_Threshold(
+                TX,
+                operator = ">",
+                threshold = TX_Base[[names(TX_Base)[2]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnSummary = percentage_non_na,
+                returnTResolution = get_index_t_resolution("TX90p")
+            )
+        },
+        WSDI = function() {
+            # Calculate spells on the full daily axis, then aggregate by run end period.
+            binary_cf <- Helper_Threshold(
+                TX,
+                operator = ">",
+                threshold = TX_Base[[names(TX_Base)[2]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnTResolution = get_index_t_resolution("WSDI")
+            )
+            run_end_weighted_series(binary_cf, minRunLength = 6)$summarise(
+                "Tresholded",
+                sum_non_na,
+                get_index_t_resolution("WSDI")
+            )
+        },
+        CSDI = function() {
+            # CSDI uses TN against TN baseline 10th percentile.
+            binary_cf <- Helper_Threshold(
+                TN,
+                operator = "<",
+                threshold = TN_Base[[names(TN_Base)[1]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = FALSE,
+                returnTResolution = get_index_t_resolution("CSDI")
+            )
+            run_end_weighted_series(binary_cf, minRunLength = 6)$summarise(
+                "Tresholded",
+                sum_non_na,
+                get_index_t_resolution("CSDI")
+            )
+        },
+        DTR = function() {
+            dtr <- TX - TN
+            dtr$summarise("Tresholded", mean_non_na, get_index_t_resolution("DTR"))
+        },
+        Rx1day = function() {
+            RR$summarise("Tresholded", max_non_na, get_index_t_resolution("Rx1day"))
+        },
+        Rx5day = function() {
+            RR$summarise("Tresholded", max_sum_over_5, get_index_t_resolution("Rx5day"))
+        },
+        SDII = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = 1,
+                returnValues = TRUE,
+                returnSummary = mean_non_na,
+                returnTResolution = get_index_t_resolution("SDII")
+            )
+        },
+        R10mm = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = 10,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("R10mm")
+            )
+        },
+        R20mm = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = 20,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("R20mm")
+            )
+        },
+        Rnnmm = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = RRThreshold,
+                returnValues = FALSE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("Rnnmm")
+            )
+        },
+        CDD = function() {
+            binary_cf <- Helper_Threshold(
+                RR,
+                operator = "<",
+                threshold = 1,
+                returnValues = FALSE,
+                returnTResolution = get_index_t_resolution("CDD")
+            )
+            # Endpoint-weighted runs ensure boundary-crossing runs are assigned at their end.
+            run_end_weighted_series(binary_cf, minRunLength = 1)$summarise(
+                "Tresholded",
+                max_non_na,
+                get_index_t_resolution("CDD")
+            )
+        },
+        CWD = function() {
+            binary_cf <- Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = 1,
+                returnValues = FALSE,
+                returnTResolution = get_index_t_resolution("CWD")
+            )
+            run_end_weighted_series(binary_cf, minRunLength = 1)$summarise(
+                "Tresholded",
+                max_non_na,
+                get_index_t_resolution("CWD")
+            )
+        },
+        R95pTOT = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">",
+                threshold = RR_Base[[names(RR_Base)[1]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = TRUE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("R95pTOT")
+            )
+        },
+        R99pTOT = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">",
+                threshold = RR_Base[[names(RR_Base)[2]]],
+                threshMode = "ETCCDIQuantiles",
+                returnValues = TRUE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("R99pTOT")
+            )
+        },
+        RCPTOT = function() {
+            Helper_Threshold(
+                RR,
+                operator = ">=",
+                threshold = 1,
+                returnValues = TRUE,
+                returnSummary = sum_non_na,
+                returnTResolution = get_index_t_resolution("RCPTOT")
+            )
+        }
     )
 
-    ### Ascribing proper Names and Dates
-    VarNames_ls <- list(
-        FD = list("FD", "Number of frost days"),
-        SU = list("SU", "Number of summer days"),
-        ID = list("ID", "Number of icing days"),
-        TR = list("TR", "Number of tropical nights"),
-        R10 = list("R10", "Days with Precip ≥ 10mm"),
-        R20 = list("R20", "Days with Precip ≥ 20mm")
+    ## calculate only selected ETCCDI objects (tResolution_ls already created before FCheck)
+    Return_ls <- list()
+    for (i in seq_along(selectedIndices)) {
+        idx <- selectedIndices[[i]]
+        Return_ls[[idx]] <- indexCalculators[[idx]]()
+        pb$tick(tokens = list(layer = i))
+    }
+
+    ## Combine into CFDataset(s) with appropriate variable names + long_name attributes
+    long_names <- c(
+        FD = "Number of frost days",
+        SU = "Number of summer days",
+        ID = "Number of icing days",
+        TR = "Number of tropical nights",
+        GSL = "Growing Season Length",
+        TXx = "Max of Daily Max Temp",
+        TNx = "Max of Daily Min Temp",
+        TXn = "Min of Daily Max Temp",
+        TNn = "Min of Daily Min Temp",
+        TN10p = "Percent Days TN < 10th Percentile",
+        TX10p = "Percent Days TX < 10th Percentile",
+        TN90p = "Percent Days TN > 90th Percentile",
+        TX90p = "Percent Days TX > 90th Percentile",
+        WSDI = "Warm Spell Duration Index",
+        CSDI = "Cold Spell Duration Index",
+        DTR = "Daily Temperature Range",
+        Rx1day = "Max 1-day Precipitation",
+        Rx5day = "Max 5-day Precipitation",
+        SDII = "Simple Precipitation Intensity Index",
+        R10mm = "Days with Precip ≥ 10mm",
+        R20mm = "Days with Precip ≥ 20mm",
+        Rnnmm = "Days with Precip ≥ user-defined threshold",
+        CDD = "Consecutive Dry Days",
+        CWD = "Consecutive Wet Days",
+        R95pTOT = "Total Precipitation from RR > 95th Percentile",
+        R99pTOT = "Total Precipitation from RR > 99th Percentile",
+        RCPTOT = "Total Precipitation on Wet Days"
     )
 
-    Return_ls <- lapply(1:length(Return_ls), FUN = function(Iter) {
-        terra::varnames(Return_ls[[Iter]]) <- VarNames_ls[[Iter]][[1]]
-        terra::longnames(Return_ls[[Iter]]) <- VarNames_ls[[Iter]][[2]]
-        terra::time(Return_ls[[Iter]]) <- Dates
-        Return_ls[[Iter]]
-    })
-    names(Return_ls) <- names(VarNames_ls)
-
-    return(Return_ls)
-}
\ No newline at end of file
+    # Determine what to return and save (unique_resolutions already created before FCheck)
+    if (length(unique_resolutions) == 1) {
+        # Single temporal resolution - base dataset is sufficient
+        return_ds <- ncdfCF::create_ncdf()
+        for (nm in names(Return_ls)) {
+            cfvar <- Return_ls[[nm]]
+            raw <- cfvar$raw()
+            outputName <- nm
+            ## adding variable with "ALT_" prefix if user provided a temporal resolution that differs from the default for this index
+            if (userProvidedTResolution && defaultTResolution[[nm]] != providedTResolution) {
+                outputName <- paste0("ALT_", nm)
+            }
+            new_var <- ncdfCF::as_CF(outputName, raw)
+            if (nm %in% names(long_names)) {
+                new_var$set_attribute("long_name", "NC_CHAR", long_names[[nm]])
+            }
+            return_ds$add_variable(new_var)
+        }
+    } else {
+        # Multiple temporal resolutions - create one dataset with one subgroup per resolution
+        return_ds <- ncdfCF::create_ncdf()
+        # List of shared objects that go into the root group (.. = parent group of the sub-group, i.e. the root group)
+        latlonnames <- c(
+            if (exists("TX")) names(TX$axes) else character(0),
+            if (exists("TN")) names(TN$axes) else character(0),
+            if (exists("RR")) names(RR$axes) else character(0)
+        )
+        latlonnames <- unique(latlonnames)
+        latlonnames <- latlonnames[!(latlonnames %in% "time")]
+        latlon <- as.list(rep("..", length(latlonnames) * 2 + 1))
+        names(latlon) <- c(latlonnames, paste0(latlonnames, "_bnds"), "height")
+
+        for (tres in unique_resolutions) {
+            subgroup_name <- as.character(tres)
+            resolution_group <- return_ds$root$create_subgroup(subgroup_name)
+            var_names_for_this_res <- names(tResolution_ls)[unlist(tResolution_ls) == tres]
+
+            for (nm in var_names_for_this_res) {
+                cfvar <- Return_ls[[nm]]
+                raw <- cfvar$raw()
+                outputName <- nm
+                if (userProvidedTResolution && defaultTResolution[[nm]] != providedTResolution) {
+                    outputName <- paste0("ALT_", nm)
+                }
+                new_var <- ncdfCF::as_CF(outputName, raw)
+                if (nm %in% names(long_names)) {
+                    new_var$set_attribute("long_name", "NC_CHAR", long_names[[nm]])
+                }
+                resolution_group$add_variable(new_var, locations = latlon)
+            }
+        }
+    }
+
+    ## saving if fileName provided, otherwise just return dataset object
+    if (!missing(fileName)) {
+        return_ds$save(fileName)
+        return_ds <- NC_Read(fileName)
+    }
+
+    ## return dataset(s) to user
+    return_ds
+}
diff --git a/data-raw/DATASET.R b/data-raw/DATASET.R
index e9aec46..a99d8c0 100644
--- a/data-raw/DATASET.R
+++ b/data-raw/DATASET.R
@@ -1,36 +1,14 @@
 # Load Package ------------------
-setwd("C:/Users/erikkus/OneDrive - CICERO senter for klimaforskning/Documents/31569 - ClimHub/CLIMHUB_GitHub")
-
-document()
-# build()
+# document()
+# # build()
 load_all()
 
-# Subsetted NORA3 file -----------------
-NORA3 <- Access_NORA3(
-  variable = "T2M", dateStart = "2001-08-01 00", dateStop = "2001-08-01 21",
-  extent = c(0, 30, 50, 75), leadTimeHour = 3,
-  fileName = file.path(getwd(), "inst/extdata", "NORA3.nc")
-)
-
 # Make NORA3 raw file ------------------
 NORA3_rast <- Access_NORA3(
   variable = "T2M", dateStart = "2001-08-01 00", dateStop = "2001-08-01 00",
   leadTimeHour = 3,
-  fileName = file.path(getwd(), "inst/extdata", "NORA3_FULL.nc")
+  fileName = file.path(getwd(), "inst/extdata", "NORA3.nc")
 )
-# usethis::use_data(NORA3_rast)
-
-# # Make Klima i Norge raw file ------------------
-# KiN_rast <- Access_KlimaiNorge2100(
-#   variable = "mean_air_temperature",
-#   dateStart = "1995-01-01",
-#   dateStop = "1995-01-10",
-#   model = "noresm-r1i1p1f1-hclim",
-#   method = "eqm",
-#   scenario = "ssp370",
-#   fileName = file.path(getwd(), "inst/extdata", "KiN_rast.nc")
-# )
-
 
 # # Jotunheimen boundary as spatialfeatureobject ------------------
 # Jotunheimen_sf <- sf::st_read("data-raw/Shape/Shape-polygon.shp")
@@ -41,128 +19,3 @@ NORA3_rast <- Access_NORA3(
 # Peaks_df$Longitude <- as.numeric(stringr::str_trim(Peaks_df$Longitude))
 # Nor2K_sf <- sf::st_as_sf(Peaks_df, coords = c("Longitude", "Latitude"), crs = sf::st_crs(Jotunheimen_sf))
 # usethis::use_data(Nor2K_sf, overwrite = TRUE)
-
-# # # Metrics.ETCCDI objects ------------------
-# # ## Raw Data -----
-# # TX_rast <- Access_KlimaiNorge2100(
-# #   variable = "Maximum Air Temperature",
-# #   dateStart = "2050-01-01",
-# #   dateStop = "2051-12-31",
-# #   scenario = "rcp85",
-# #   model = "CNRM_CCLM",
-# #   fileName = file.path(getwd(), "inst/extdata", "KiN_TX.nc"),
-# #   writeFile = FALSE
-# # )
-
-# # TN_rast <- Access_KlimaiNorge2100(
-# #   variable = "Minimum Air Temperature",
-# #   dateStart = "2050-01-01",
-# #   dateStop = "2051-12-31",
-# #   scenario = "rcp85",
-# #   model = "CNRM_CCLM",
-# #   fileName = file.path(getwd(), "inst/extdata", "KiN_TN.nc"),
-# #   writeFile = FALSE
-# # )
-
-# # BP_TM_rast <- Access_KlimaiNorge2100(
-# #   variable = "Mean Air Temperature",
-# #   dateStart = "1971-01-01",
-# #   dateStop = "2000-12-31",
-# #   model = "CNRM_CCLM",
-# #   fileName = file.path(getwd(), "inst/extdata", "KiN_BP_TM.nc"),
-# #   writeFile = FALSE
-# # )
-
-# # RR_rast <- Access_KlimaiNorge2100(
-# #   variable = "Precipitation",
-# #   dateStart = "2050-01-01",
-# #   dateStop = "2051-12-31",
-# #   scenario = "rcp85",
-# #   model = "CNRM_CCLM",
-# #   fileName = file.path(getwd(), "inst/extdata", "KiN_RR.nc"),
-# #   writeFile = FALSE
-# # )
-# # # RR_rast <- RR_rast / 10
-
-# # BP_RR_rast <- Access_KlimaiNorge2100(
-# #   variable = "Precipitation",
-# #   dateStart = "1971-01-01",
-# #   dateStop = "2000-12-31",
-# #   model = "CNRM_CCLM",
-# #   fileName = file.path(getwd(), "inst/extdata", "KiN_BP_RR.nc"),
-# #   writeFile = FALSE
-# # )
-# # # BP_RR_rast <- BP_RR_rast / 10
-
-# # ## Cropping -----
-# # datals <- list(
-# #   TX = TX_rast, TN = TN_rast, RR = RR_rast,
-# #   BP_TM = BP_TM_rast, BP_RR = BP_RR_rast
-# # )
-
-# # x <- lapply(names(datals), function(name) {
-# #   # name <- "X"
-# #   print(name)
-# #   FNAME <- file.path(getwd(), "inst/extdata", paste0("Sognefjord_", name, ".nc"))
-# #   if (file.exists(FNAME)) {
-# #     ret <- terra::rast(FNAME)
-# #   } else {
-# #     x <- datals[[name]]
-# #     # y <- Spatial_Reproject(x, Jotunheimen_sf)
-# #     z <- Spatial_Limit(spatRaster = x, shape = ext(c(1e4, 1.5e5, 6.76e6, 6.86e+06)))
-# #     if (unique(terra::units(x)) == "K") {
-# #       # z <- z / 100
-# #     } else {
-# #       z <- z / 10
-# #       terra::units(z) <- "mm"
-# #     }
-# #     # terra::plot(z[[1]])
-# #     ret <- NC_Write(
-# #       spatRaster = z,
-# #       fileName = FNAME,
-# #       varName = unique(terra::varnames(x)),
-# #       longName = unique(terra::longnames(x)),
-# #       unit = unique(terra::units(x)),
-# #       meta = setNames(terra::metags(x)[, 2], terra::metags(x)[, 1]),
-# #       compression = 9
-# #     )
-# #   }
-# #   ret
-# # })
-# # unlink(list.files(pattern = "TEMP_"))
-
-# # ## Percentiles -----
-# # stop("Make percentiles for BASEPeriods")
-# # ## make Reference into percentiles
-# # library(terra)
-# # BASEPeriod <- terra::rast("inst/extdata/Jotunheimen_BASEPeriod.nc")
-# # percentiles <- app(BASEPeriod, fun = function(x) {
-# #   quantile(x, probs = c(0.1, 0.9), na.rm = TRUE)
-# # })
-# # # Optional: give meaningful layer names
-# # names(percentiles) <- c("p10", "p90")
-# # terra::time(percentiles) <- rep(time(BASEPeriod[[1]]), 2)
-
-# # BASEPeriod_rast <- NC_Write(
-# #   spatRaster = percentiles / 100,
-# #   fileName = file.path(getwd(), "inst/extdata", "Jotunheimen_BP_TM.nc"),
-# #   varName = "Air Temperature Percentiles",
-# #   longName = "Air Temperature Percentiles for 1971-2000",
-# #   unit = "K",
-# #   meta = setNames(terra::metags(BASEPeriod_rast)[, 2], terra::metags(BASEPeriod_rast)[, 1]),
-# #   compression = 9
-# # )
-
-# # BASEPeriod_rast <- terra::rast("inst/extdata/Jotunheimen_BASEPeriod.nc")
-
-
-# # unlink(list.files(pattern = "TEMP_", full.names = TRUE))
-
-# # ### Sufficiently large file for testing of NC_Write ------------------
-# # NORA3 <- Access_NORA3(
-# #   variable = "TS (Surface temperature)", # which variable
-# #   dateStart = "1971-01-01 00", dateStop = "1971-05-31 18", # time-window
-# #   leadTimeHour = 3, cores = 1,
-# #   fileName = "NORA3.nc", compression = NA, # file storing
-# #   removeTemporary = TRUE
-# )
diff --git a/inst/extdata/NORA3.nc b/inst/extdata/NORA3.nc
index 9b6f7b4..b34759f 100644
Binary files a/inst/extdata/NORA3.nc and b/inst/extdata/NORA3.nc differ
diff --git a/inst/extdata/NORA3_FULL.nc b/inst/extdata/NORA3_FULL.nc
deleted file mode 100644
index a61a00f..0000000
Binary files a/inst/extdata/NORA3_FULL.nc and /dev/null differ
diff --git a/inst/figures/CodeBadge.png b/inst/figures/CodeBadge.png
index fa6e38a..68df1df 100644
Binary files a/inst/figures/CodeBadge.png and b/inst/figures/CodeBadge.png differ
diff --git a/inst/figures/PackageWorkflow.png b/inst/figures/PackageWorkflow.png
index 68e51c1..5e204ce 100644
Binary files a/inst/figures/PackageWorkflow.png and b/inst/figures/PackageWorkflow.png differ
diff --git a/man/Helper_ETCCDIDailyBootstraps.Rd b/man/Helper_ETCCDIDailyBootstraps.Rd
new file mode 100644
index 0000000..c7ca7c1
--- /dev/null
+++ b/man/Helper_ETCCDIDailyBootstraps.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Helper_ETCCDIDailyBootstraps.r
+\name{Helper_ETCCDIDailyBootstraps}
+\alias{Helper_ETCCDIDailyBootstraps}
+\title{Prepare Bootstrap Intervals of Daily Time Series across Years}
+\usage{
+Helper_ETCCDIDailyBootstraps(dates, bootstrapWindow = 5)
+}
+\arguments{
+\item{dates}{A character vector of dates in the format "YYYY-MM-DDTHH:MM:SS" (e.g., "1961-01-01T00:00:00") representing the time series for which bootstrap intervals should be created. Same as is obtained with the \code{dimnames(ds[["time"]])} of a \code{CFDataset} object called \code{ds}.}
+
+\item{bootstrapWindow}{Numeric, uneven integer. Length of bootstrapping interval around each target date.}
+}
+\value{
+A list of character vectors, each containing the dates that fall within the bootstrap window around each target date. The names of the list correspond to the day of year (in "MM-DD" format) for which the bootstrap intervals are created. Interval creation assumes a leap year (i.e., 366 days) to ensure that the same day of year is always included in the bootstrap intervals, even for February 29th.
+}
+\description{
+Create bootstrap intervals of daily time series across years for a given time sequence. This function is needed for baseperiod quantile calculations for \code{Metrics_ETCCDI}.
+}
+\examples{
+time_series <- seq(
+    as.POSIXct("1961-01-01 00:00:00", tz = "UTC"),
+    as.POSIXct("1999-12-31 23:00:00", tz = "UTC"),
+    by = "1 day"
+)
+Time_seq <- format(time_series, "\%Y-\%m-\%dT\%H:\%M:\%S")
+Helper_ETCCDIDailyBootstraps(dates = Time_seq)
+
+}
+\author{
+Erik Kusch
+}
diff --git a/man/Helper_ETCCDIGSL.Rd b/man/Helper_ETCCDIGSL.Rd
new file mode 100644
index 0000000..619f0b7
--- /dev/null
+++ b/man/Helper_ETCCDIGSL.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Helper_ETCCDIGSL.r
+\name{Helper_ETCCDIGSL}
+\alias{Helper_ETCCDIGSL}
+\title{Calculate Growing Season Length (GSL)}
+\usage{
+Helper_ETCCDIGSL(TM)
+}
+\arguments{
+\item{TM}{A CFVariable containing daily mean temperatures (or similar) to be used for GSL calculation.}
+}
+\value{
+A CFVariable.
+}
+\description{
+Calculates Growing Season Length as defined by ETCCDI with different treatment of northern and southern hemisphere. For the northern hemisphere, GSL is defined as the number of days between the first occurrence of 6 consecutive days with mean temperature above 5°C after Jan 1st and the first occurrence of 6 consecutive days with mean temperature below 5°C after Jul 1st. For the southern hemisphere, GSL is defined as the number of days between the first occurrence of 6 consecutive days with mean temperature above 5°C after Jul 1st and the first occurrence of 6 consecutive days with mean temperature below 5°C after Jan 1st of the following year. If no valid start or end date is found, GSL is set to NA for that year and grid cell.
+}
+\examples{
+\dontrun{
+TX <- NC_Read("inst/extdata/KiN_tasmax_2050.nc")[["tasmax"]]
+TN <- NC_Read("inst/extdata/KiN_tasmin_2050.nc")[["tasmin"]]
+TM <- (TN + TX) / 2
+GSL <- Helper_ETCCDIGSL(TM)
+}
+}
+\author{
+Erik Kusch
+}
diff --git a/man/Helper_Threshold.Rd b/man/Helper_Threshold.Rd
index 9ca7ed7..1861e73 100644
--- a/man/Helper_Threshold.Rd
+++ b/man/Helper_Threshold.Rd
@@ -2,45 +2,52 @@
 % Please edit documentation in R/Helper_Threshold.r
 \name{Helper_Threshold}
 \alias{Helper_Threshold}
-\title{Apply tresholding to raster layers}
+\title{Apply tresholding to CFVariable contents}
 \usage{
 Helper_Threshold(
-  spatRaster,
+  CFVariable,
   operator,
   threshold,
+  threshMode = NULL,
   returnValues = FALSE,
-  verbose = TRUE
+  returnSummary,
+  returnTResolution = "year"
 )
 }
 \arguments{
-\item{spatRaster}{A SpatRaster within which thresholding should be applied}
+\item{CFVariable}{A CFVariable within which thresholding should be applied}
 
-\item{operator}{Character. One of ">", "<", ">=", "<=", "==", "!=" to evaluate spatRaster cells against threshold value.}
+\item{operator}{Character. One of ">", "<", ">=", "<=", "==", "!=" to evaluate CFVariable cells against threshold value.}
 
-\item{threshold}{Numeric or Character. threshold to evaluate spatRaster cells against using operator value.}
+\item{threshold}{Numeric or Character. threshold to evaluate CFVariable cells against using operator value.}
 
-\item{returnValues}{Logical. Whether to return spatRaster of values matching thresholding (TRUE) or only logical raster indicating where thresholding is matched (FALSE).}
+\item{threshMode}{Character, optional. If supplying a CFVariable containing baseline quantiles for ETCCDI calculation, setting \code{threshMode = "ETCCDIQuantiles"} will internally apply location and day-of-year specific thresholds to \code{CFVariable} input. Default is \code{NULL} avoiding this behaviour.}
 
-\item{verbose}{Logical. If progress should be displayed in the console.}
+\item{returnValues}{Logical. Whether to return CFVariable of values matching thresholding (TRUE) or only logical raster indicating where thresholding is matched (FALSE).}
+
+\item{returnSummary}{Function. Summary function for tresholded CFVariable.}
+
+\item{returnTResolution}{Character. Temporal resolution of summary of tresholded CFVariable.}
 }
 \value{
-A SpatRaster.
+A CFVariable.
 }
 \description{
-Loops over layers in a SpatRaster object and evaluates threshold criteria defined by user.
+Evaluates threshold criteria defined by user against a CFVariable and reports summaries through the \verb{$summarise} method of \code{ncdfCF}.
 }
 \examples{
-Data_rast <- terra::rast(system.file("extdata", "KiN_rast.nc", package = "ClimHub"))[[1:31]]
+Data_CF <- NC_Read("inst/extdata/KiN_tas.nc")[["tas"]]
 Helper_Threshold(
-    spatRaster = Data_rast,
+    CFVariable = Data_CF,
     operator = "<",
-    threshold = 2650
+    threshold = 273.15
 )
 Helper_Threshold(
-    spatRaster = Data_rast,
+    CFVariable = Data_CF,
     operator = "<",
-    threshold = 2650,
-    returnValues = TRUE
+    threshold = 273.15,
+    returnValues = TRUE,
+    returnSummary = mean
 )
 }
 \author{
diff --git a/man/Metrics_BootstrapQuantiles.Rd b/man/Metrics_BootstrapQuantiles.Rd
new file mode 100644
index 0000000..62b7736
--- /dev/null
+++ b/man/Metrics_BootstrapQuantiles.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Metrics_BootstrapQuantiles.r
+\name{Metrics_BootstrapQuantiles}
+\alias{Metrics_BootstrapQuantiles}
+\title{Bootstrap Quantile Lookup for Daily Time Series}
+\usage{
+Metrics_BootstrapQuantiles(
+  CFVariable,
+  probs = c(0.1, 0.9),
+  bootstrapWindow = 5
+)
+}
+\arguments{
+\item{CFVariable}{\code{CFVariable} object for which quantiles for each bootstrap day interval ought to be calculated.}
+
+\item{probs}{Numeric vector of probabilities (between 0 and 1) for which quantiles will be returned.  Default is \code{c(0.1,0.9)}.}
+
+\item{bootstrapWindow}{Odd integer giving the width of the window in days used for the bootstrap sample around each target date. The window is centred on the day of year of interest.  Defaults to 5.}
+}
+\value{
+A numeric array with dimensions \code{lon × lat × dayOfYear × probs}.  The third dimension has 366 levels ("01‑01" through "12‑31", including Feb 29) and the fourth dimension is named by the supplied \code{probs}.  Values are the requested quantiles computed over the set of years that fall inside the bootstrap window.  If no valid observations exist for a given cell/day the result is NA.
+}
+\description{
+Computes quantiles across years for each day‑of‑year using a rolling bootstrap window to be used in \code{Metrics_ETTCDI} calculation of ETCCDI. Internally the function uses \code{Helper_ETCCDIDailyBootstraps()} to determine which dates belong in the window centred on each calendar day (assuming a leap year for consistency).
+
+This is primarily used by the ETCCDI metric code to obtain base period thresholds but can be called directly by users. Dates in the input \code{CFVariable} must present individual days.
+}
+\examples{
+\dontrun{
+TX <- NC_Read("inst/extdata/KiN_tasmax.nc")[["tasmax"]]
+TX_Quant <- Metrics_BootstrapQuantiles(CFVariable = TX, probs = c(0.1, 0.9), bootstrapWindow = 5)
+TX_Quant
+}
+}
+\author{
+Erik Kusch
+}
diff --git a/man/Metrics_ETCCDI.Rd b/man/Metrics_ETCCDI.Rd
index 96317ac..4420b8a 100644
--- a/man/Metrics_ETCCDI.Rd
+++ b/man/Metrics_ETCCDI.Rd
@@ -4,30 +4,174 @@
 \alias{Metrics_ETCCDI}
 \title{Calculate ETCCDI Indices from List Input}
 \usage{
-Metrics_ETCCDI(ls)
+Metrics_ETCCDI(
+  projectionList,
+  baseLineList,
+  indices,
+  TResolution,
+  RRThreshold = 42,
+  fileName
+)
 }
 \arguments{
-\item{ls}{Optional, list. List of SpatRasters. Names of elements must be "TX" and "TN", holding maximum and minimum daily air temperature respectively.}
+\item{projectionList}{List. List of \code{CFVariable} objects required by the selected ETCCDI indices via the \code{indices} argument. Include only named elements that are needed for \code{indices}: "TX" (daily maximum air temperature in Kelvin), "TN" (daily minimum air temperature in Kelvin), and/or "RR" (daily total precipitation in mm). See details for required data input per index.}
+
+\item{baseLineList}{Optional, list. List of \code{CFDataset} objects containing baseline quantiles required only for quantile-based ETCCDI indices. Include only named elements needed for \code{indices}: "TX_Base", "TN_Base", and/or "RR_Base" (note that these must be in Kelvin, Kelvin, and mm, respectively). If no quantile-based ETCCDI is selected, this argument can be omitted. See details for required quantile baselines per index.}
+
+\item{indices}{Optional, character. Character vector of ETCCDI abbreviations to calculate (see first column of the details table). If missing, all supported indices in this function are calculated. See details for supported indices.}
+
+\item{TResolution}{Optional, character. Temporal resolution for ETCCDI calculation. Supports "year", "month" and "season". If omitted, each selected index is calculated using its default temporal resolution from the details table. If provided, the same resolution is used for all selected indices and output indices whose default resolution differs are prefixed with "ALT_". Note that this is not applied to GSL which is always calculated at annual resolution due to its specific logic.}
+
+\item{RRThreshold}{Numeric. Custom threshold for daily precipiation in mm for calculation of Rnnmm. Defaults to 42.}
+
+\item{fileName}{Character, optional. Character. A file name for the produced file, including path and ".nc" file ending. If no value is supplied, the dataset is not written to disk but returned as a \code{CFDataset} object in memory. If a file name is supplied and a file with that name already exists, the function will attempt to load and return that file instead of recalculating the indices.}
 }
 \value{
-A named list of SpatRasters with each element corresponding to an ETCCDI.
+A \code{CFDataset}. If all calculated ETCCDI indices have the same temporal resolution (e.g., when \code{TResolution} is specified or a subset of \code{indices} is selected that share a temporal resolution default), the returned \code{CFDataset} contains all variables in the root group. If indices span multiple temporal resolutions (e.g., when \code{TResolution} is omitted and some selected \code{indices} default to "year" while others default to "month"), the returned \code{CFDataset} contains one subgroup per temporal resolution, each containing only variables of that temporal resolution.
+
+Each variable is named by its ETCCDI acronym and has a \code{long_name} attribute describing the index. Acronyms of indices calculated with a non-default temporal resolution are prefixed with "ALT_" to indicate that these are alternative calculations.
 }
 \description{
-This function calculates \href{https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/jgrd.50203}{ETCCDIs} from a named list of SpatRaster objects. Currently, only the first four ETCCDI are supported.
+This function calculates \href{https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/jgrd.50203}{ETCCDIs} from a named list of CFVariable objects. Currently, only some ETCCDI are supported.
+}
+\details{
+Input requirements are conditional on the selected ETCCDI indices in the \code{indices} argument (see table below). To reduce memory usage and speed up calculations, supply only the inputs required for the selected indices.\tabular{llllll}{
+   \strong{ETCCDI} \tab \strong{Name} \tab \strong{Required Data} \tab \strong{Description} \tab \strong{Unit} \tab \strong{Default \code{TResolution}} \cr
+   FD \tab Number of frost days \tab \code{TN} \tab Count of days when TN (daily minimum temperature) < 0°C. \tab days \tab year \cr
+   SU \tab Number of summer days \tab \code{TX} \tab Count of days when TX (daily maximum temperature) > 25°C. \tab days \tab year \cr
+   ID \tab Number of icing days \tab \code{TX} \tab Count of days when TX (daily maximum temperature) < 0°C. \tab days \tab year \cr
+   TR \tab Number of tropical nights \tab \code{TN} \tab Count of days when TN (daily minimum temperature) > 20°C. \tab days \tab year \cr
+   GSL \tab Growing Season Length \tab \code{TX}, \code{TN} \tab Annual count of days between the first occurrence of at least 6 consecutive days with (TN+TX)/2 > 5°C and the first occurrence after 1st July (Northern Hemisphere) or 1st January (Southern Hemisphere) of at least 6 consecutive days with (TN+TX)/2 < 5°C \tab days \tab year \cr
+   TXx \tab Max of Daily Max Temp \tab \code{TX} \tab Maximum daily maximum temperature. \tab K \tab month \cr
+   TNx \tab Max of Daily Min Temp \tab \code{TN} \tab Maximum daily minimum temperature. \tab K \tab month \cr
+   TXn \tab Min of Daily Max Temp \tab \code{TX} \tab Minimum daily maximum temperature. \tab K \tab month \cr
+   TNn \tab Min of Daily Min Temp \tab \code{TN} \tab Minimum daily minimum temperature. \tab K \tab month \cr
+   TN10p \tab Percent Days TN < 10th Percentile \tab \code{TN}, \code{TN_Base} \tab Percent of days where TN < 10th percentile of base period. \tab \% \tab year \cr
+   TX10p \tab Percent Days TX < 10th Percentile \tab \code{TX}, \code{TX_Base} \tab Percent of days where TX < 10th percentile of base period. \tab \% \tab year \cr
+   TN90p \tab Percent Days TN > 90th Percentile \tab \code{TN}, \code{TN_Base} \tab Percent of days where TN > 90th percentile of base period. \tab \% \tab year \cr
+   TX90p \tab Percent Days TX > 90th Percentile \tab \code{TX}, \code{TX_Base} \tab Percent of days where TX > 90th percentile of base period. \tab \% \tab year \cr
+   WSDI \tab Warm Spell Duration Index \tab \code{TX}, \code{TX_Base} \tab Count of days contained within runs of 6+ consecutive days when TX > 90th percentile of base period. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. \tab days \tab year \cr
+   CSDI \tab Cold Spell Duration Index \tab \code{TN}, \code{TN_Base} \tab Count of days contained within runs of 6+ consecutive days when TN < 10th percentile of base period. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. \tab days \tab year \cr
+   DTR \tab Daily Temperature Range \tab \code{TX}, \code{TN} \tab Mean difference between daily max (TX) and min (TN) temperatures. \tab K \tab month \cr
+   Rx1day \tab Max 1-day Precipitation \tab \code{RR} \tab Maximum precipitation in a single day. \tab mm \tab month \cr
+   Rx5day \tab Max 5-day Precipitation \tab \code{RR} \tab Maximum precipitation over any 5 consecutive days. \tab mm \tab month \cr
+   SDII \tab Simple Precipitation Intensity Index \tab \code{RR} \tab Mean precipitation amount on wet days (RR ≥ 1mm). \tab mm/day \tab year \cr
+   R10mm \tab Days with Precip ≥ 10mm \tab \code{RR} \tab Count of days with precipitation ≥ 10mm. \tab days \tab year \cr
+   R20mm \tab Days with Precip ≥ 20mm \tab \code{RR} \tab Count of days with precipitation ≥ 20mm. \tab days \tab year \cr
+   Rnnmm \tab Days with Precip ≥ user-defined threshold \tab \code{RR} \tab Count of days with precipitation ≥ nnmm. \tab days \tab year \cr
+   CDD \tab Consecutive Dry Days \tab \code{RR} \tab Maximum number of consecutive days with RR < 1mm. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. \tab days \tab year \cr
+   CWD \tab Consecutive Wet Days \tab \code{RR} \tab Maximum number of consecutive days with RR ≥ 1mm. Note that periods of consecutive days are attributed to the time window in which the sequence terminates. \tab days \tab year \cr
+   R95pTOT \tab Total Precipitation from RR > 95th Percentile \tab \code{RR}, \code{RR_Base} \tab Total precipitation from wet days (RR > 95th percentile of base period). \tab mm \tab year \cr
+   R99pTOT \tab Total Precipitation from RR > 99th Percentile \tab \code{RR}, \code{RR_Base} \tab Total precipitation from wet days (RR > 99th percentile of base period). \tab mm \tab year \cr
+   RCPTOT \tab Total Precipitation on Wet Days \tab \code{RR} \tab Sum of precipitation on wet days (RR ≥ 1mm). \tab mm \tab year \cr
+}
 }
 \examples{
 \dontrun{
-TX_rast <- terra::rast("inst/extdata/Jotunheimen_TX.nc") # terra::rast(system.file("extdata", "Jotunheimen_TX.nc", package = "ClimHub"))
-TN_rast <- terra::rast("inst/extdata/Jotunheimen_TN.nc")
-RR_rast <- terra::rast("inst/extdata/Jotunheimen_RR.nc")
-terra::units(RR_rast) <- "mm"
-# BASEPeriod_rast <- terra::rast("inst/extdata/Jotunheimen_BASEPeriod.nc")
-
-Metrics_ETCCDI(ls = list(TX = TX_rast, TN = TN_rast, RR = RR_rast
-# , BASEPeriod = BASEPeriod_rast)
+## Directory for Data and Output -----------
+Dir.Data <- file.path(getwd(), "ExampleData")
+dir.create(Dir.Data, showWarnings = FALSE)
+
+## Downloading and preparing data -----------
+### Variables for which we need data
+vars <- c("minimum_air_temperature", "maximum_air_temperature", "precipitation_flux")
+names(vars) <- c("tasmin", "tasmax", "pr")
+
+### Data (Down-)Loading and Preparation
+Data_ls <- lapply(seq_along(vars), FUN = function(i) {
+    print(as.character(vars[i]))
+    ## these need to be made into quantiles to define base-period for ETCCDI calculation
+    KiN_Base <- Access_KlimaiNorge2100(
+        variable = as.character(vars[i]),
+        dateStart = "1971-01-01", # note that the standard base period starts 1961, but KiN does start in 1971 so we use that as the start date
+        dateStop = "2000-12-31",
+        extent = c(9, 11, 59, 61),
+        model = "noresm-r1i1p1f1-hclim",
+        method = "eqm",
+        scenario = "ssp370",
+        fileName = file.path(Dir.Data, paste0("KiN_", names(vars)[i], ".nc"))
+    )
+
+    ## these are the files for calculation of ETCCDI metrics in future
+    KiN_2090 <- Access_KlimaiNorge2100(
+        variable = as.character(vars[i]),
+        dateStart = "2090-01-01",
+        dateStop = "2099-12-31",
+        extent = c(9, 11, 59, 61),
+        model = "noresm-r1i1p1f1-hclim",
+        method = "eqm",
+        scenario = "ssp370",
+        fileName = file.path(Dir.Data, paste0("KiN_", names(vars)[i], "_2090.nc"))
+    )
+
+    ## these are the files for baseline quantiles
+    QuantF <- file.path(Dir.Data, paste0("KiN_", names(vars)[i], "_BaseLineQuantiles.nc"))
+    if (!file.exists(QuantF)) {
+        if (names(vars)[i] == "pr") {
+            probs_vec <- c(0.95, 0.99)
+            input <- KiN_Base[[names(vars)[i]]] * 86400 # to get from mm/day to kg m-2 s-1
+            input$set_attribute("units", "NC_CHAR", "mm")
+        } else {
+            probs_vec <- c(0.1, 0.9)
+            input <- KiN_Base[[names(vars)[i]]]
+        }
+        Quant <- Metrics_BootstrapQuantiles(
+            CFVariable = input,
+            probs = probs_vec,
+            bootstrapWindow = 5
+        )
+        Quant$save(QuantF)
+    } else {
+        Quant <- NC_Read(QuantF)
+    }
+
+    return(list(projection = KiN_2090, baseLine = Quant))
+})
+names(Data_ls) <- names(vars)
+
+### Data Object listing
+TX <- Data_ls[["tasmax"]][["projection"]][["tasmax"]]
+TN <- Data_ls[["tasmin"]][["projection"]][["tasmin"]]
+RR <- Data_ls[["pr"]][["projection"]][["pr"]] * 86400 # to get from mm/day to kg m-2 s-1
+RR$set_attribute("units", "NC_CHAR", "mm")
+
+TX_Base <- Data_ls[["tasmax"]][["baseLine"]]
+TN_Base <- Data_ls[["tasmin"]][["baseLine"]]
+RR_Base <- Data_ls[["pr"]][["baseLine"]]
+
+rm(Data_ls) # remove list to free up memory
+
+## ETCCDI calculation -----------
+### All ETCCDI
+Metrics_ETCCDI(
+    projectionList = list(TX = TX, TN = TN, RR = RR),
+    baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+    fileName = file.path(Dir.Data, "ETCCDI_Metrics.nc")
 )
-}
 
+### Select ETCCDI
+Metrics_ETCCDI(
+    projectionList = list(TX = TX, TN = TN, RR = RR),
+    baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+    indices = c("FD", "TXx", "R95pTOT")
+)
+
+### Select ETCCDI with custom temporal resolutions
+Metrics_ETCCDI(
+    projectionList = list(RR = RR),
+    # baseLineList = list(RR_Base = RR_Base), # we do not need to supply the baseline quantiles if we do not select any indices that depend on them
+    indices = c("Rnnmm"),
+    RRThreshold = 21, # this is the threshold for the Rnnmm metric, in mm
+    TResolution = "season"
+)
+
+Metrics_ETCCDI(
+    projectionList = list(TX = TX, TN = TN, RR = RR),
+    baseLineList = list(TX_Base = TX_Base, TN_Base = TN_Base, RR_Base = RR_Base),
+    indices = c("FD", "SU", "ID", "TR"),
+    TResolution = "month" # note that this is not the default for these indices, so the returned object will bear the prefix "ALT_" for each index to indicate that this is an alternative calculation with a different temporal resolution
+)
+}
 }
 \author{
 Erik Kusch