diff --git a/R/CST_BiasCorrection.R b/R/CST_BiasCorrection.R
index 20b51082ccc7a164213e26d4844b36ea05dafe48..6b9f53989bd207cdcd65d8c7a6d13d5160605ef0 100644
--- a/R/CST_BiasCorrection.R
+++ b/R/CST_BiasCorrection.R
@@ -19,6 +19,8 @@
 #'  dimension. By default, it is set to 'member'.
 #'@param sdate_dim A character string indicating the name of the start date 
 #'  dimension. By default, it is set to 'sdate'.
+#'@param window_dim A character string indicating the name of the dimension to 
+#' enlarge the sample data. By default NULL.
 #'@param ncores An integer that indicates the number of cores for parallel 
 #'  computations using multiApply function. The default value is NULL.
 #'@return An object of class \code{s2dv_cube} containing the bias corrected 
@@ -46,7 +48,7 @@
 #'@export
 CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE, 
                                memb_dim = 'member', sdate_dim = 'sdate',
-                               ncores = NULL) {
+                               window_dim = NULL, ncores = NULL) {
   if (!inherits(exp, 's2dv_cube') || !inherits(obs, 's2dv_cube')) {
     stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
          "as output by CSTools::CST_Load.")
@@ -63,6 +65,7 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
 
   BiasCorrected <- BiasCorrection(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, 
                                   memb_dim = memb_dim, sdate_dim = sdate_dim,
+                                  window_dim = window_dim,
                                   na.rm = na.rm, ncores = ncores)
   
   pos <- match(dimnames, names(dim(BiasCorrected)))
@@ -105,6 +108,8 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
 #'  dimension. By default, it is set to 'member'.
 #'@param sdate_dim A character string indicating the name of the start date 
 #'  dimension. By default, it is set to 'sdate'.
+#'@param window_dim A character string indicating the name of the dimension to 
+#' enlarge the sample data. By default NULL.
 #'@param ncores An integer that indicates the number of cores for parallel 
 #'  computations using multiApply function. The default value is NULL.
 #'
@@ -127,7 +132,7 @@ CST_BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
 #'@export
 BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
                            memb_dim = 'member', sdate_dim = 'sdate',
-                           ncores = NULL) {
+                           window_dim = NULL, ncores = NULL) {
   # Check inputs
   ## exp, obs
   if (!is.array(exp) || !is.numeric(exp)) {
@@ -196,42 +201,67 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
   if (memb_dim %in% names(dim(obs))) {
     target_dims_obs <- c(memb_dim, target_dims_obs)
   }
-  
-  if (is.null(exp_cor)) {
-    BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp),
-                           target_dims = list(target_dims_obs,
+  if (is.null(window_dim)) {
+    if (is.null(exp_cor)) {
+      BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp),
+                             target_dims = list(target_dims_obs,
                                               c(memb_dim, sdate_dim)),
-                           fun = .sbc, 
-                           na.rm = na.rm, ncores = ncores)$output1
-  } else {
-    BiasCorrected <- Apply(data = list(var_obs = obs,
-                                       var_exp = exp,
-                                       var_cor = exp_cor),
+                             fun = .sbc, 
+                             na.rm = na.rm, ncores = ncores)$output1
+    } else {
+      BiasCorrected <- Apply(data = list(var_obs = obs,
+                                        var_exp = exp,
+                                        var_cor = exp_cor),
                            target_dims = list(target_dims_obs,
                                               c(memb_dim, sdate_dim),
                                               c(memb_dim, sdate_dim)),
                            fun = .sbc, 
                            output_dims = c(memb_dim, sdate_dim),
                            na.rm = na.rm, ncores = ncores)$output1
-  }
+    }
+  } else {
+     if (is.null(exp_cor)) {
+      BiasCorrected <- Apply(data = list(var_obs = obs, var_exp = exp),
+                             target_dims = list(c(target_dims_obs, window_dim),
+                                              c(memb_dim, sdate_dim, window_dim)),
+                             fun = .sbc,
+                             na.rm = na.rm, ncores = ncores)$output1
+    } else {
+      BiasCorrected <- Apply(data = list(var_obs = obs,
+                                        var_exp = exp,
+                                        var_cor = exp_cor),
+                           target_dims = list(c(target_dims_obs, window_dim),
+                                              c(memb_dim, sdate_dim, window_dim),
+                                              c(memb_dim, sdate_dim, window_dim)),
+                           fun = .sbc,
+                           output_dims = c(memb_dim, sdate_dim),
+                           na.rm = na.rm, ncores = ncores)$output1
+    }
+  }  
   return(BiasCorrected)
 }
 
 .sbc <- function(var_obs, var_exp , var_cor = NULL, na.rm = FALSE) {
-
+  #var_obs: [member, sdate, (swindow)] ---> subseasonal: member, syear, sday
+  #var_exp: [member, sdate, (swindow)] --> subseasonal: member, syear, sday
+  #var_cor: [member, sdate, (swindow)] --> subsesonal: member, syear, sday
   ntime <- dim(var_exp)[2]
-  corrected <- NA * var_exp
 
   if (is.null(var_cor)) {
+     corrected <- NA * var_exp
     for (t in 1:ntime) {
       # parameters
-      sd_obs <- sd(var_obs[-t], na.rm = na.rm)
-      sd_exp <- sd(var_exp[, -t], na.rm = na.rm)
-      clim_exp <- mean(var_exp[, -t], na.rm = na.rm)
-      clim_obs <- mean(var_obs[-t], na.rm = na.rm)
-      
-      # bias corrected forecast
-      corrected[, t] <- ((var_exp[, t] - clim_exp) * (sd_obs / sd_exp)) + clim_obs
+      sd_obs <- sd(Subset(var_obs, along = 2, indices = -t), na.rm = na.rm)
+      sd_exp <- sd(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm)
+      clim_exp <- mean(Subset(var_exp, along = 2, indices = -t), na.rm = na.rm)
+      clim_obs <- mean(Subset(var_obs, along = 2, indices = -t), na.rm = na.rm)
+      if (length(dim(var_exp)) == 2) {
+        corrected[,t] <- ((Subset(var_exp, along = 2, indices = t) - clim_exp) *
+                      (sd_obs/sd_exp)) + clim_obs
+      } else {
+        corrected[,t,] <- ((Subset(var_exp, along = 2, indices = t) - clim_exp) *
+                      (sd_obs/sd_exp)) + clim_obs
+      }
     }
   } else {
     # parameters
@@ -239,9 +269,14 @@ BiasCorrection <- function(exp, obs, exp_cor = NULL, na.rm = FALSE,
     sd_exp <- sd(var_exp, na.rm = na.rm)
     clim_exp <- mean(var_exp, na.rm = na.rm)
     clim_obs <- mean(var_obs, na.rm = na.rm)
-    
-    # bias corrected forecast
-    corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs
+    if (length(dim(var_cor)) == 2) {
+      # bias corrected forecast
+      corrected <- ((var_cor - clim_exp) * (sd_obs / sd_exp)) + clim_obs
+    } else {
+      # Todo: flexible position of the window subset?
+      corrected <- ((Subset(var_cor, along = 3, indices = 2, drop = 'selected') - clim_exp) *
+                    (sd_obs / sd_exp)) + clim_obs
+    }
   }
 
   return(corrected)
diff --git a/man/BiasCorrection.Rd b/man/BiasCorrection.Rd
index d944a32fa7f61e738e907bdc3c95e1971912f707..20ff7260981d361341ea7025548df16f8fc166ed 100644
--- a/man/BiasCorrection.Rd
+++ b/man/BiasCorrection.Rd
@@ -11,6 +11,7 @@ BiasCorrection(
   na.rm = FALSE,
   memb_dim = "member",
   sdate_dim = "sdate",
+  window_dim = NULL,
   ncores = NULL
 )
 }
@@ -35,6 +36,9 @@ dimension. By default, it is set to 'member'.}
 \item{sdate_dim}{A character string indicating the name of the start date 
 dimension. By default, it is set to 'sdate'.}
 
+\item{window_dim}{A character string indicating the name of the dimension to 
+enlarge the sample data. By default NULL.}
+
 \item{ncores}{An integer that indicates the number of cores for parallel 
 computations using multiApply function. The default value is NULL.}
 }
diff --git a/man/CST_BiasCorrection.Rd b/man/CST_BiasCorrection.Rd
index fb96babd17b91ea27c8b4c4c0d794f33d8c5a734..4e4bcee8323d463c236bfe36f2154e0ed4a66a31 100644
--- a/man/CST_BiasCorrection.Rd
+++ b/man/CST_BiasCorrection.Rd
@@ -11,6 +11,7 @@ CST_BiasCorrection(
   na.rm = FALSE,
   memb_dim = "member",
   sdate_dim = "sdate",
+  window_dim = NULL,
   ncores = NULL
 )
 }
@@ -35,6 +36,9 @@ dimension. By default, it is set to 'member'.}
 \item{sdate_dim}{A character string indicating the name of the start date 
 dimension. By default, it is set to 'sdate'.}
 
+\item{window_dim}{A character string indicating the name of the dimension to 
+enlarge the sample data. By default NULL.}
+
 \item{ncores}{An integer that indicates the number of cores for parallel 
 computations using multiApply function. The default value is NULL.}
 }