diff --git a/R/SelectorChecker.R b/R/SelectorChecker.R
index 1bf2681c917abfbe6d2f3283e5af457f3df8ac9d..141cae2178fb7c25f18d1171f470f4bf3cdc15ad 100644
--- a/R/SelectorChecker.R
+++ b/R/SelectorChecker.R
@@ -212,7 +212,7 @@ SelectorChecker <- function(selectors, var = NULL, return_indices = TRUE,
   } else {
     if (!is.null(var)) {
       if (is.list(selectors)) {
-        if (length(selectors != 2)) {
+        if (length(selectors) != 2) {
           stop("'selectors' provided in a wrong format.")
         } else {
           var[selectors[[1]]:selectors[[2]]]
diff --git a/R/Start.R b/R/Start.R
index ef9db3e3d02150faa5beff19afd9b412c74b3cee..eec20acd0b4e9cd583964d27fa6a3be99633931f 100644
--- a/R/Start.R
+++ b/R/Start.R
@@ -1983,10 +1983,19 @@ print("-> THE INNER DIMENSION DOES NOT GO ACROSS ANY FILE DIMENSION OR IT DOES B
                                                                            tolerance_params[[inner_dim]]
                                                                          })
                   }
-                  sub_array_of_indices <- sub_array_of_indices[chunk_indices(length(sub_array_of_indices),
-                                                                             chunks[[inner_dim]]['chunk'],
-                                                                             chunks[[inner_dim]]['n_chunks'],
-                                                                             inner_dim)]
+                  if (!is.list(sub_array_of_indices)) {
+                    sub_array_of_indices <- sub_array_of_indices[chunk_indices(length(sub_array_of_indices), 
+                                                                               chunks[[inner_dim]]["chunk"], 
+                                                                               chunks[[inner_dim]]["n_chunks"], 
+                                                                               inner_dim)]
+                  } else {
+                    tmp <- chunk_indices(length(sub_array_of_indices[[1]]:sub_array_of_indices[[2]]),
+                                                chunks[[inner_dim]]["chunk"], chunks[[inner_dim]]["n_chunks"],
+                                                inner_dim)
+                    start_pt <- sub_array_of_indices[[1]]
+                    sub_array_of_indices[[1]] <- start_pt + tmp[1] - 1
+                    sub_array_of_indices[[2]] <- start_pt + tmp[length(tmp)] - 1
+                  }
                   # The sub_array_of_indices now contains numeric indices of the values to be taken.
 
 #Check if all the files have the selectors assigned (e.g., region = 'Grnland') _20191015
diff --git a/inst/doc/usecase.md b/inst/doc/usecase.md
index a584f1e10a0a45d986ee7b9c2756bd77a897e512..2bfb989349453b8882a9aff242c2cc4450ac7a67 100644
--- a/inst/doc/usecase.md
+++ b/inst/doc/usecase.md
@@ -17,3 +17,4 @@ In this document, you can link to the example scripts for various demands. For t
    5. [Using experimental and (date-corresponding) observational data](inst/doc/usecase/ex2_5_exp_and_obs.R)
    6. [Use external parameters in atomic function](inst/doc/usecase/ex2_6_ext_param_func.R)
    7. [Seasonal forecast verification on cca](inst/doc/usecase/ex2_7_seasonal_forecast_verification.R)
+   8. [Use CSTools Calibration function](inst/doc/usecase/ex2_8_calibration.R)
diff --git a/inst/doc/usecase/ex2_8_calibration.R b/inst/doc/usecase/ex2_8_calibration.R
new file mode 100644
index 0000000000000000000000000000000000000000..13b0739bb49c5e7c0d456a249c6723d02a3be3ba
--- /dev/null
+++ b/inst/doc/usecase/ex2_8_calibration.R
@@ -0,0 +1,85 @@
+# Load startR
+library(startR)
+
+# Define a region
+lons.min <- 0
+lons.max <- 10
+lats.min <- 0
+lats.max <- 10
+
+# Declaration of data sources
+exp <- Start(dat = '/esarchive/exp/ecmwf/system4_m1/monthly_mean/$var$_f6h/$var$_$sdate$.nc',
+             var = 'tas',
+             sdate = paste0(2000:2010, '0101'),
+             ensemble = 'all',
+             time = indices(1),
+             latitude = values(list(lats.min, lats.max)),
+             latitude_reorder = Sort(),
+             longitude = values(list(lons.min, lons.max)),
+             return_vars = list(latitude = 'dat',
+                                longitude = 'dat',
+                                time = c('sdate')),
+             retrieve = FALSE)
+
+obs <- Start(dat = '/esarchive/recon/ecmwf/erainterim/monthly_mean/$var$_f6h/$var$_$sdate$.nc',
+             var = 'tas',
+             sdate = paste0(2000:2010, '02'),
+             time = indices(1),
+             latitude = values(list(lats.min, lats.max)),
+             latitude_reorder = Sort(),
+             longitude = values(list(lons.min, lons.max)),
+             return_vars = list(latitude = 'dat',
+                                longitude = 'dat',
+                                time = c('sdate')),
+             split_multiselected_dims = TRUE,
+             merge_across_dims = TRUE,
+             retrieve = FALSE)
+
+# Define of the workflow
+
+# Function
+wrap_cal <- function(obs, exp) {
+    obs <- s2dverification::InsertDim(obs, 1, 1)
+    names(dim(obs))<- c('member', 'sdate')
+    exp <- t(exp)
+    names(dim(exp))<- c('member', 'sdate')
+    calibrated <- CSTools:::.cal(var_obs = obs, var_exp = exp) # CSTools version 1.0.1 or earlier
+    # calibrated <- CSTools:::.cal(mod = exp, obs = obs,
+    #                             cal.method = "mse_min",
+    #                             eval.method = "leave-one-out",
+    #                             multi.model = FALSE)
+    return(calibrated)
+}
+
+step <- Step(wrap_cal,
+             target_dims = list(obs = c('sdate'), exp = c('sdate', 'ensemble')),
+             output_dims = c('ensemble', 'sdate'))
+
+# workflow of operations
+wf <- AddStep(list(obs = obs, exp = exp), step)
+
+
+# Execution
+res <- Compute(wf, chunks = list(latitude = 2, longitude = 2),
+               threads_load = 2, threads_compute = 4)
+## Check output:
+#summary(res$output1)
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+#  294.1   298.8   300.5   300.6   302.2   306.7 
+
+
+# Declaration of HPC and execution
+## ECFlow is required
+res_fat2 <- Compute(wf,chunks = list(latitude = 2, longitude = 2),
+                    threads_load = 2,
+                    threads_compute = 4,
+                    cluster = list(queue_host = "bsceslogin01.bsc.es",
+                                   cores_per_job = 2,
+                                   max_jobs = 4, job_wallclock = '00:10:00'),
+                    ecflow_suite_dir = "/esarchive/scratch/nperez/ecflow") # your path!
+## Check output:
+#summary(res_fat2$output1)
+#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
+#  294.1   298.8   300.5   300.6   302.2   306.7 
+
+