diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R
index 79f51320f60fab62195235ee2b53969896cf7483..e973c4d84dae5973e4f02653ea92747d13c98c1d 100644
--- a/R/CST_Calibration.R
+++ b/R/CST_Calibration.R
@@ -147,33 +147,25 @@
 #'@export
 CST_Calibration <- function(exp, obs, exp_cor = NULL, cal.method = "mse_min", 
                             eval.method = "leave-one-out", multi.model = FALSE, 
-                            na.fill = TRUE, na.rm = TRUE, apply_to = NULL, alpha = NULL,
-                            memb_dim = 'member', sdate_dim = 'sdate', dat_dim = NULL, ncores = NULL) {
+                            na.fill = TRUE, na.rm = TRUE, apply_to = NULL, 
+                            alpha = NULL, memb_dim = 'member', sdate_dim = 'sdate', 
+                            dat_dim = NULL, ncores = NULL) {
   # Check 's2dv_cube'
   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.")
+        stop("Parameter 'exp' and 'obs' must be of the class 's2dv_cube'.")
   }
   if (!is.null(exp_cor)) {
     if (!inherits(exp_cor, "s2dv_cube")) {
-        stop("Parameter 'exp', 'obs' and 'exp_cor' must be of the class 's2dv_cube', ",
-             "as output by CSTools::CST_Load.")
+        stop("Parameter 'exp_cor' must be of the class 's2dv_cube'.")
     }
-    # if exp_cor is provided, it will be calibrated: "calibrate forecast instead of hindcast"
-    # if exp_cor is provided, eval.method is overruled (because if exp_cor is provided, the 
-    # train data will be all data of "exp" and the evalutaion data will be all data of "exp_cor"; 
-    # no need for "leave-one-out" or "in-sample")
-    eval.method = "hindcast-vs-forecast" 
-  }
-
-  if (!missing(multi.model) & !(cal.method == "mse_min")) {
-	  warning(paste0("The multi.model parameter is ignored when using the calibration method ", cal.method))
   }
 
   Calibration <- Calibration(exp = exp$data, obs = obs$data, exp_cor = exp_cor$data, 
-                             cal.method = cal.method, eval.method = eval.method, multi.model =  multi.model, 
-                             na.fill = na.fill, na.rm = na.rm, apply_to = apply_to, alpha = alpha, 
-                             memb_dim = memb_dim, sdate_dim = sdate_dim, dat_dim = dat_dim, ncores = ncores)
+                             cal.method = cal.method, eval.method = eval.method, 
+                             multi.model =  multi.model, na.fill = na.fill, 
+                             na.rm = na.rm, apply_to = apply_to, alpha = alpha, 
+                             memb_dim = memb_dim, sdate_dim = sdate_dim, 
+                             dat_dim = dat_dim, ncores = ncores)
 
   if (is.null(exp_cor)) {
     exp$data <- Calibration
@@ -312,7 +304,8 @@ Calibration <- function(exp, obs, exp_cor = NULL,
                         cal.method = "mse_min", eval.method = "leave-one-out",  
                         multi.model = FALSE, na.fill = TRUE, 
                         na.rm = TRUE, apply_to = NULL, alpha = NULL,
-                        memb_dim = 'member', sdate_dim = 'sdate', dat_dim = NULL, ncores = NULL) {
+                        memb_dim = 'member', sdate_dim = 'sdate', dat_dim = NULL, 
+                        ncores = NULL) {
 
   # Check inputs
   ## exp, obs
@@ -338,7 +331,11 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   }
   ## exp_cor
   if (!is.null(exp_cor)) {
-    eval.method = "hindcast-vs-forecast"
+    # if exp_cor is provided, it will be calibrated: "calibrate forecast instead of hindcast"
+    # if exp_cor is provided, eval.method is overruled (because if exp_cor is provided, the 
+    # train data will be all data of "exp" and the evalutaion data will be all data of "exp_cor"; 
+    # no need for "leave-one-out" or "in-sample")
+    eval.method <- "hindcast-vs-forecast"
     expcordims <- names(dim(exp_cor))
     if (is.null(expcordims)) {
       stop("Parameter 'exp_cor' must have dimension names.")
@@ -449,25 +446,9 @@ Calibration <- function(exp, obs, exp_cor = NULL,
       stop("Parameter 'ncores' must be either NULL or a positive integer.")
     }
   }
-  ## data.set.sufficiently.large
-  data.set.sufficiently.large.out <- Apply(data = list(exp = exp, obs = obs),
-                                           target_dims = list(exp = target_dims_exp, 
-                                           obs = target_dims_obs),
-                                           ncores = ncores,
-                                           fun = .data.set.sufficiently.large)$output1
-  
-  if (!all(data.set.sufficiently.large.out)) {		
-  	if (na.fill) {
-        warning("Some forecast data could not be corrected due to data lack",
-                " and is replaced with NA values")
-  	} else {
-        warning("Some forecast data could not be corrected due to data lack",
-                " and is replaced with uncorrected values")
-  	 }
-  }
   ## na.rm
-  if (!na.rm %in% c(TRUE,FALSE)) {
-    stop("Parameter 'na.rm' must be TRUE or FALSE.")
+  if (!inherits(na.rm, "logical")) {
+    stop("Parameter 'na.rm' must be a logical value.")
   }
   if (length(na.rm) > 1) {
     na.rm <- na.rm[1]
@@ -480,28 +461,39 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   if (cal.method == 'rpc-based') {
     if (is.null(apply_to)) {
       apply_to <- 'sign'
-      warning("'apply_to' cannot be NULL for 'rpc-based' method so it has been set to 'sign', as in Eade et al. (2014).")
+      warning("Parameter 'apply_to' cannot be NULL for 'rpc-based' method so it ", 
+              "has been set to 'sign', as in Eade et al. (2014).")
     } else if (!apply_to %in% c('all','sign')) {
-      stop("'apply_to' must be either 'all' or 'sign' when 'rpc-based' method is used.")
+      stop("Parameter 'apply_to' must be either 'all' or 'sign' when 'rpc-based' ", 
+           "method is used.")
     }
     if (apply_to == 'sign') {
       if (is.null(alpha)) {
         alpha <- 0.1
-        warning("'alpha' cannot be NULL for 'rpc-based' method so it has been set to 0.1, as in Eade et al. (2014).")
+        warning("Parameter 'alpha' cannot be NULL for 'rpc-based' method so it ", 
+                "has been set to 0.1, as in Eade et al. (2014).")
       } else if (!is.numeric(alpha) | alpha <= 0 | alpha >= 1) {
-        stop("'alpha' must be a number between 0 and 1.")
+        stop("Parameter 'alpha' must be a number between 0 and 1.")
       }
     }
   }
   ## eval.method
   if (!any(eval.method %in% c('in-sample', 'leave-one-out', 'hindcast-vs-forecast'))) {
-    stop("Parameter 'eval.method' must be a character string indicating the sampling method used ('in-sample', 'leave-one-out' or 'hindcast-vs-forecast').")
+    stop(paste0("Parameter 'eval.method' must be a character string indicating ", 
+                "the sampling method used ('in-sample', 'leave-one-out' or ", 
+                "'hindcast-vs-forecast')."))
   }
   ## multi.model
-  if (!multi.model %in% c(TRUE,FALSE)) {
-    stop("Parameter 'multi.model' must be TRUE or FALSE.")
+  if (!inherits(multi.model, "logical")) {
+    stop("Parameter 'multi.model' must be a logical value.")
+  }
+  if (multi.model & !(cal.method == "mse_min")) {
+	  warning(paste0("The 'multi.model' parameter is ignored when using the ", 
+                   "calibration method '", cal.method, "'."))
   }
 
+  warning_shown <- FALSE
+
   if (is.null(exp_cor)) {
     calibrated <- Apply(data = list(exp = exp, obs = obs), dat_dim = dat_dim, 
                         cal.method = cal.method, eval.method = eval.method, multi.model = multi.model,
@@ -519,7 +511,7 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   }
   if (!is.null(dat_dim)) {
     pos <- match(c(names(dim(exp))[-which(names(dim(exp)) == dat_dim)], 'nexp', 'nobs'), 
-                  names(dim(calibrated)))
+                 names(dim(calibrated)))
     calibrated <- aperm(calibrated, pos)
   } else {
     pos <- match(c(names(dim(exp))), names(dim(calibrated)))
@@ -529,8 +521,6 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   if (exp_cor_remove_memb) {
     dim(calibrated) <- dim(calibrated)[-which(names(dim(calibrated)) == memb_dim)]
   }
-
-
   return(calibrated)
 }
 
@@ -576,13 +566,13 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   }
 
   if (is.null(exp_cor)) {
-    # generate a copy of exp
+    # generate a copy of exp so that the same function can run for both cases
     exp_cor <- exp
     cor_dat_dim <- TRUE
   } else {
-    if (length(dim(exp_cor)) == 2) { # ref: [memb, sdate]
+    if (length(dim(exp_cor)) == 2) { # exp_cor: [memb, sdate]
       cor_dat_dim <- FALSE
-    } else {
+    } else {                         # exp_cor: [memb, sdate, dat]
       cor_dat_dim <- TRUE
     }  
   }
@@ -593,119 +583,39 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   sdate <- expdims[2] # sdate
   sdate_cor <- expdims_cor[2]
 
-  return_data_na <- FALSE
   var.cor.fc <- array(dim = c(dim(exp_cor)[1:2], nexp = nexp, nobs = nobs))
+  
   for (i in 1:nexp) {
     for (j in 1:nobs) {
-      if (!.data.set.sufficiently.large(exp = exp[, , i, drop = FALSE], obs = obs[, j, drop = FALSE])) {
-        return_data_na <- TRUE
-        if (cor_dat_dim) {
-          var.cor.fc[, , i, j] <- NA * exp_cor[, , i]
-        } else {
-          var.cor.fc[, , i, j] <- NA * exp_cor[, i]
-        }
+      if (!.data.set.sufficiently.large(exp = exp[, , i, drop = FALSE], 
+                                        obs = obs[, j, drop = FALSE])) {
         if (!na.fill) {
-          var.cor.fc[, , i] <- exp[, , i]
-        }
-      }
-      obs_data <- as.vector(obs[, j])
-      exp_data <- exp[, , i]
-      dim(exp_data) <- dim(exp)[1:2]
-      if (cor_dat_dim) {
-        expcor_data <- exp_cor[, , i]
-        dim(expcor_data) = dim(exp_cor)[1:2]
-      } else {
-        expcor_data <- exp_cor
-      }
-      eval.train.dexeses <- .make.eval.train.dexes(eval.method = eval.method, amt.points = sdate, 
-                                                   amt.points_cor = sdate_cor)
-      amt.resamples <- length(eval.train.dexeses)
-      for (i.sample in seq(1, amt.resamples)) {
-        # defining training (tr) and evaluation (ev) subsets
-        # fc.ev is used to evaluate (not train; train should be done with exp (hindcast))
-        eval.dexes <- eval.train.dexeses[[i.sample]]$eval.dexes
-        train.dexes <- eval.train.dexeses[[i.sample]]$train.dexes
-        fc.ev <- expcor_data[ , eval.dexes, drop = FALSE]
-        fc.tr <- exp_data[ , train.dexes]
-        obs.tr <- obs_data[train.dexes , drop = FALSE] 
-        
-        if (cal.method == "bias") {
-          var.cor.fc[ , eval.dexes, i, j] <- fc.ev + mean(obs.tr, na.rm = na.rm) - mean(fc.tr, na.rm = na.rm)
-          # forecast correction implemented
-        } else if (cal.method == "evmos") {
-          # forecast correction implemented
-          # ensemble and observational characteristics
-          quant.obs.fc.tr <- .calc.obs.fc.quant(obs = obs.tr, fc = fc.tr, na.rm = na.rm)
-          # calculate value for regression parameters
-          init.par <- c(.calc.evmos.par(quant.obs.fc.tr, na.rm = na.rm))
-          # correct evaluation subset
-          var.cor.fc[ , eval.dexes, i, j] <- .correct.evmos.fc(fc.ev , init.par, na.rm = na.rm)
-        } else if (cal.method == "mse_min") {
-          quant.obs.fc.tr <- .calc.obs.fc.quant(obs = obs.tr, fc = fc.tr, na.rm = na.rm)
-          init.par <- .calc.mse.min.par(quant.obs.fc.tr, multi.model, na.rm = na.rm)
-          var.cor.fc[ , eval.dexes, i, j] <- .correct.mse.min.fc(fc.ev , init.par, na.rm = na.rm)      
-        } else if (cal.method == "crps_min") {
-          quant.obs.fc.tr <- .calc.obs.fc.quant.ext(obs = obs.tr, fc = fc.tr, na.rm = na.rm)
-          init.par <- c(.calc.mse.min.par(quant.obs.fc.tr, na.rm = na.rm), 0.001)
-          init.par[3] <- sqrt(init.par[3])
-          # calculate regression parameters on training dataset
-          optim.tmp <- optim(par = init.par, fn = .calc.crps.opt, gr = .calc.crps.grad.opt, 
-                              quant.obs.fc = quant.obs.fc.tr, na.rm = na.rm, method = "BFGS")
-          mbm.par <- optim.tmp$par
-          var.cor.fc[ , eval.dexes, i, j] <- .correct.crps.min.fc(fc.ev , mbm.par, na.rm = na.rm)
-        } else if (cal.method == 'rpc-based') {
-          # Ensemble mean
-          ens_mean.ev <- Apply(data = fc.ev, target_dims = names(memb), fun = mean, na.rm = na.rm)$output1
-          ens_mean.tr <- Apply(data = fc.tr, target_dims = names(memb), fun = mean, na.rm = na.rm)$output1
-          # Ensemble spread
-          ens_spread.tr <- Apply(data = list(fc.tr, ens_mean.tr), target_dims = names(sdate), fun = "-")$output1
-          # Mean (climatology)
-          exp_mean.tr <- mean(fc.tr, na.rm = na.rm) 
-          # Ensemble mean variance
-          var_signal.tr <- var(ens_mean.tr, na.rm = na.rm) 
-          # Variance of ensemble members about ensemble mean (= spread)
-          var_noise.tr <- var(as.vector(ens_spread.tr), na.rm = na.rm)
-          # Variance in the observations 
-          var_obs.tr <- var(obs.tr, na.rm = na.rm)
-          # Correlation between observations and the ensemble mean 
-          r.tr <- cor(x = ens_mean.tr, y = obs.tr, method = 'pearson', 
-                      use = ifelse(test = isTRUE(na.rm), yes = "pairwise.complete.obs", no = "everything"))
-          if ((apply_to == 'all') || (apply_to == 'sign' && 
-              cor.test(ens_mean.tr, obs.tr, method = 'pearson', alternative = 'greater')$p.value < alpha)) {
-            ens_mean_cal <- (ens_mean.ev - exp_mean.tr) * r.tr * sqrt(var_obs.tr) / sqrt(var_signal.tr) + exp_mean.tr
-            var.cor.fc[ , eval.dexes, i, j] <- Reorder(data = Apply(data = list(exp = fc.ev, ens_mean = ens_mean.ev, 
-                                                                    ens_mean_cal = ens_mean_cal), 
-                                                                    target_dims = names(sdate), fun = .CalibrationMembersRPC, 
-                                                                    var_obs = var_obs.tr, var_noise = var_noise.tr, r = r.tr)$output1, 
-                                                        order = names(expdims)[1:2])
-          } else {
-            # no significant -> replacing with observed climatology
-            var.cor.fc[ , eval.dexes, i, j] <- array(data = mean(obs.tr, na.rm = na.rm), dim = dim(fc.ev))
+          exp_subset <- exp[, , i]
+          var.cor.fc[, , i, j] <- exp_subset
+          if (!warning_shown) {
+            warning("Some forecast data could not be corrected due to data lack",
+                    " and is replaced with uncorrected values.")
+            warning_shown <<- TRUE
           }
-        } else {
-          stop("unknown calibration method: ", cal.method)
+        } else if (!warning_shown) {
+          warning("Some forecast data could not be corrected due to data lack",
+                  " and is replaced with NA values.")
+          warning_shown <<- TRUE
         }
-      }
-    }
-  }
-
-
-
-
-  if (!return_data_na) {
-    var.cor.fc <- array(dim = c(dim(exp_cor)[1:2], nexp = nexp, nobs = nobs))
-    for (i in 1:nexp) {
-      for (j in 1:nobs) {
+      } else {
+        # Subset data for dataset dimension
         obs_data <- as.vector(obs[, j])
         exp_data <- exp[, , i]
         dim(exp_data) <- dim(exp)[1:2]
         if (cor_dat_dim) {
           expcor_data <- exp_cor[, , i]
-          dim(expcor_data) = dim(exp_cor)[1:2]
+          dim(expcor_data) <- dim(exp_cor)[1:2]
         } else {
           expcor_data <- exp_cor
         }
-        eval.train.dexeses <- .make.eval.train.dexes(eval.method = eval.method, amt.points = sdate, 
+
+        eval.train.dexeses <- .make.eval.train.dexes(eval.method = eval.method, 
+                                                     amt.points = sdate, 
                                                      amt.points_cor = sdate_cor)
         amt.resamples <- length(eval.train.dexeses)
         for (i.sample in seq(1, amt.resamples)) {
@@ -713,12 +623,12 @@ Calibration <- function(exp, obs, exp_cor = NULL,
           # fc.ev is used to evaluate (not train; train should be done with exp (hindcast))
           eval.dexes <- eval.train.dexeses[[i.sample]]$eval.dexes
           train.dexes <- eval.train.dexeses[[i.sample]]$train.dexes
-          fc.ev <- expcor_data[ , eval.dexes, drop = FALSE]
-          fc.tr <- exp_data[ , train.dexes]
-          obs.tr <- obs_data[train.dexes , drop = FALSE] 
+          fc.ev <- expcor_data[, eval.dexes, drop = FALSE]
+          fc.tr <- exp_data[, train.dexes]
+          obs.tr <- obs_data[train.dexes, drop = FALSE] 
           
           if (cal.method == "bias") {
-            var.cor.fc[ , eval.dexes, i, j] <- fc.ev + mean(obs.tr, na.rm = na.rm) - mean(fc.tr, na.rm = na.rm)
+            var.cor.fc[, eval.dexes, i, j] <- fc.ev + mean(obs.tr, na.rm = na.rm) - mean(fc.tr, na.rm = na.rm)
             # forecast correction implemented
           } else if (cal.method == "evmos") {
             # forecast correction implemented
@@ -727,11 +637,11 @@ Calibration <- function(exp, obs, exp_cor = NULL,
             # calculate value for regression parameters
             init.par <- c(.calc.evmos.par(quant.obs.fc.tr, na.rm = na.rm))
             # correct evaluation subset
-            var.cor.fc[ , eval.dexes, i, j] <- .correct.evmos.fc(fc.ev , init.par, na.rm = na.rm)
+            var.cor.fc[, eval.dexes, i, j] <- .correct.evmos.fc(fc.ev , init.par, na.rm = na.rm)
           } else if (cal.method == "mse_min") {
             quant.obs.fc.tr <- .calc.obs.fc.quant(obs = obs.tr, fc = fc.tr, na.rm = na.rm)
             init.par <- .calc.mse.min.par(quant.obs.fc.tr, multi.model, na.rm = na.rm)
-            var.cor.fc[ , eval.dexes, i, j] <- .correct.mse.min.fc(fc.ev , init.par, na.rm = na.rm)      
+            var.cor.fc[, eval.dexes, i, j] <- .correct.mse.min.fc(fc.ev , init.par, na.rm = na.rm)      
           } else if (cal.method == "crps_min") {
             quant.obs.fc.tr <- .calc.obs.fc.quant.ext(obs = obs.tr, fc = fc.tr, na.rm = na.rm)
             init.par <- c(.calc.mse.min.par(quant.obs.fc.tr, na.rm = na.rm), 0.001)
@@ -740,7 +650,7 @@ Calibration <- function(exp, obs, exp_cor = NULL,
             optim.tmp <- optim(par = init.par, fn = .calc.crps.opt, gr = .calc.crps.grad.opt, 
                                quant.obs.fc = quant.obs.fc.tr, na.rm = na.rm, method = "BFGS")
             mbm.par <- optim.tmp$par
-            var.cor.fc[ , eval.dexes, i, j] <- .correct.crps.min.fc(fc.ev , mbm.par, na.rm = na.rm)
+            var.cor.fc[, eval.dexes, i, j] <- .correct.crps.min.fc(fc.ev , mbm.par, na.rm = na.rm)
           } else if (cal.method == 'rpc-based') {
             # Ensemble mean
             ens_mean.ev <- Apply(data = fc.ev, target_dims = names(memb), fun = mean, na.rm = na.rm)$output1
@@ -761,14 +671,14 @@ Calibration <- function(exp, obs, exp_cor = NULL,
             if ((apply_to == 'all') || (apply_to == 'sign' && 
                 cor.test(ens_mean.tr, obs.tr, method = 'pearson', alternative = 'greater')$p.value < alpha)) {
               ens_mean_cal <- (ens_mean.ev - exp_mean.tr) * r.tr * sqrt(var_obs.tr) / sqrt(var_signal.tr) + exp_mean.tr
-              var.cor.fc[ , eval.dexes, i, j] <- Reorder(data = Apply(data = list(exp = fc.ev, ens_mean = ens_mean.ev, 
-                                                                      ens_mean_cal = ens_mean_cal), 
-                                                                      target_dims = names(sdate), fun = .CalibrationMembersRPC, 
-                                                                      var_obs = var_obs.tr, var_noise = var_noise.tr, r = r.tr)$output1, 
-                                                         order = names(expdims)[1:2])
+              var.cor.fc[, eval.dexes, i, j] <- Reorder(data = Apply(data = list(exp = fc.ev, ens_mean = ens_mean.ev, 
+                                                                     ens_mean_cal = ens_mean_cal), 
+                                                                     target_dims = names(sdate), fun = .CalibrationMembersRPC, 
+                                                                     var_obs = var_obs.tr, var_noise = var_noise.tr, r = r.tr)$output1, 
+                                                        order = names(expdims)[1:2])
             } else {
               # no significant -> replacing with observed climatology
-              var.cor.fc[ , eval.dexes, i, j] <- array(data = mean(obs.tr, na.rm = na.rm), dim = dim(fc.ev))
+              var.cor.fc[, eval.dexes, i, j] <- array(data = mean(obs.tr, na.rm = na.rm), dim = dim(fc.ev))
             }
           } else {
             stop("unknown calibration method: ", cal.method)
@@ -782,7 +692,7 @@ Calibration <- function(exp, obs, exp_cor = NULL,
     dim(var.cor.fc) <- dim(exp_cor)[1:2]
   }
 
-  return(var.cor.fc) 
+  return(var.cor.fc)
 }
 
 # Function to calculate different quantities of a series of ensemble forecasts and corresponding observations
@@ -861,7 +771,7 @@ Calibration <- function(exp, obs, exp_cor = NULL,
 }
 
 # Extended function to calculate different quantities of a series of ensemble forecasts
-.calc.fc.quant.ext <- function(fc, na.rm){ 
+.calc.fc.quant.ext <- function(fc, na.rm) { 
   amt.mbr <- dim(fc)[1]
   repmat1.tmp <- InsertDim(fc, posdim = 1, lendim = amt.mbr, name = 'amt.mbr')
   repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3))
@@ -888,7 +798,7 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   return(par.out)
 }
 
-.calc.evmos.par <- function(quant.obs.fc, na.rm){
+.calc.evmos.par <- function(quant.obs.fc, na.rm) {
   par.out <- rep(NA, 2)
   par.out[2] <- with(quant.obs.fc, obs.sd / fc.sd)
   par.out[1] <- with(quant.obs.fc, obs.av - par.out[2] * fc.ens.av.av, na.rm = na.rm)
@@ -927,17 +837,17 @@ Calibration <- function(exp, obs, exp_cor = NULL,
   quant.fc.mp <- .calc.fc.quant(fc = fc, na.rm = na.rm)
   return(with(quant.fc.mp, par[1] + par[2] * fc))
 }
-.correct.mse.min.fc <- function(fc, par, na.rm){
+.correct.mse.min.fc <- function(fc, par, na.rm) {
   quant.fc.mp <- .calc.fc.quant(fc = fc, na.rm = na.rm)
   return(with(quant.fc.mp, par[1] + par[2] * fc.ens.av.per.ens + fc.dev * par[3]))
 }
-.correct.crps.min.fc <- function(fc, par, na.rm){
+.correct.crps.min.fc <- function(fc, par, na.rm) {
   quant.fc.mp <- .calc.fc.quant.ext(fc = fc, na.rm = na.rm)
   return(with(quant.fc.mp, par[1] + par[2] * fc.ens.av.per.ens + fc.dev * abs((par[3])^2 + par[4] / spr.abs)))
 }
 
 # Function to calibrate the individual members with the RPC-based method
-.CalibrationMembersRPC <- function(exp, ens_mean, ens_mean_cal, var_obs, var_noise, r){
+.CalibrationMembersRPC <- function(exp, ens_mean, ens_mean_cal, var_obs, var_noise, r) {
   member_cal <- (exp - ens_mean) * sqrt(var_obs) * sqrt(1 - r^2) / sqrt(var_noise) + ens_mean_cal
   return(member_cal)
 }
diff --git a/tests/testthat/test-CST_Calibration.R b/tests/testthat/test-CST_Calibration.R
index 898acbed1054b42b5a2914e94cb22ec2a7de1e09..e6aef6f641559e7a35a379bb0bd24d294a6110fc 100644
--- a/tests/testthat/test-CST_Calibration.R
+++ b/tests/testthat/test-CST_Calibration.R
@@ -64,13 +64,11 @@ test_that("1. Input checks", {
   # s2dv_cube
   expect_error(
     CST_Calibration(exp = 1),
-    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
-         "as output by CSTools::CST_Load.")
+    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube'.")
   )
   expect_error(
     CST_Calibration(exp = exp, obs = 1),
-    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube', ",
-             "as output by CSTools::CST_Load.")
+    paste0("Parameter 'exp' and 'obs' must be of the class 's2dv_cube'.")
   )
   expect_error(
     CST_Calibration(exp = exp1),
@@ -78,8 +76,7 @@ test_that("1. Input checks", {
   )
   expect_error(
     CST_Calibration(exp = exp1, obs = obs1, exp_cor = 1),
-    paste0("Parameter 'exp', 'obs' and 'exp_cor' must be of the class 's2dv_cube', ",
-             "as output by CSTools::CST_Load.")
+    paste0("Parameter 'exp_cor' must be of the class 's2dv_cube'.")
   )
   # exp and obs
   expect_error(
@@ -172,9 +169,9 @@ test_that("1. Input checks", {
   )
   expect_error(
     Calibration(exp = array(1:6, c(sdate = 3, member = 2, dataset = 2, lon = 1)), 
-                   obs = array(1:6, c(sdate = 3, member = 1, dataset = 3, lon = 1)),
-                   exp_cor = array(1:6, c(sdate = 3, member = 1, dataset = 3, lon = 1)), 
-                   dat_dim = 'dataset'),
+                obs = array(1:6, c(sdate = 3, member = 1, dataset = 3, lon = 1)),
+                exp_cor = array(1:6, c(sdate = 3, member = 1, dataset = 3, lon = 1)), 
+                dat_dim = 'dataset'),
     paste0("If parameter 'exp_cor' has dataset dimension, it must be", 
            " equal to dataset dimension of 'exp'.")
   )
@@ -191,27 +188,59 @@ test_that("1. Input checks", {
     "Parameter 'ncores' must be either NULL or a positive integer."
   )
   # na.rm
+  expect_error(
+    CST_Calibration(exp = exp, obs = obs, na.rm = 1),
+    "Parameter 'na.rm' must be a logical value."
+  )
   expect_warning(
     CST_Calibration(exp = exp, obs = obs, na.rm = c(T,F)),
     "Paramter 'na.rm' has length greater than 1, and only the fist element is used."
   )
   # cal.method
   expect_error(
-    as.vector(Calibration(exp4, obs4, cal.method = 'biass')),
+    Calibration(exp4, obs4, cal.method = 'biass'),
     "Parameter 'cal.method' must be a character string indicating the calibration method used."
   )
-})
-##############################################
-
-test_that("2. Output checks: dat1", {
+  expect_warning(
+    Calibration(exp4, obs4, cal.method = 'rpc-based'),
+    paste0("Parameter 'apply_to' cannot be NULL for 'rpc-based' method so it ", 
+           "has been set to 'sign', as in Eade et al. (2014)."), 
+    fixed = TRUE
+  )
+  expect_warning(
+    Calibration(exp4, obs4, cal.method = 'rpc-based', apply_to = 'sign'),
+    paste0("Parameter 'alpha' cannot be NULL for 'rpc-based' method so it ", 
+           "has been set to 0.1, as in Eade et al. (2014)."), 
+    fixed = TRUE
+  )
   expect_error(
-    CST_Calibration(exp = 1),
-    "Parameter 'exp' and 'obs' must be of the class 's2dv_cube', "
+    Calibration(exp4, obs4, cal.method = 'rpc-based', apply_to = 'sign', alpha = 'a'),
+    "Parameter 'alpha' must be a number between 0 and 1."
   )
+  # eval.method
   expect_error(
-    CST_Calibration(obs = 1),
-    c("argument \"exp\" is missing, with no default")
+    Calibration(exp4, obs4, eval.method = 'biass'),
+    paste0("Parameter 'eval.method' must be a character string indicating ", 
+           "the sampling method used ('in-sample', 'leave-one-out' or ", 
+           "'hindcast-vs-forecast')."), 
+    fixed = TRUE
   )
+  # multi.model
+  expect_error(
+    Calibration(exp4, obs4, multi.model = 'biass'),
+    "Parameter 'multi.model' must be a logical value."
+  )
+  expect_warning(
+    Calibration(exp4, obs4, multi.model = TRUE, cal.method = 'bias'),
+    paste0("The 'multi.model' parameter is ignored when using the ", 
+           "calibration method 'bias'.")
+  )
+})
+
+##############################################
+
+test_that("2. Output checks: dat1", {
+  # dat_dim = NULL
   cal <- CST_Calibration(exp = exp, obs = obs)
   expect_equal(
     length(cal), 
@@ -251,6 +280,8 @@ test_that("2. Output checks: dat1", {
 ##############################################
 
 test_that("3. Output checks: dat3", {
+  # dat_dim = 'dataset'
+  # exp_cor = NULL
   expect_equal(
     dim(Calibration(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')),
     c(member = 10, sdate = 10, lat = 2, nexp = 2, nobs = 3)
@@ -269,11 +300,21 @@ test_that("3. Output checks: dat3", {
     c(0.4908992, 0.3054365, 0.3673404, 1.5218074, -1.0928551, 1.6905885),
     tolerance = 0.0001
   )
+  expect_equal(
+    as.vector(Calibration(exp3[, , , 1], obs3[, , 1], memb_dim = 'member'))[1:5],
+    c(-0.4799875, 0.3896246, -0.7045297, 1.9049701, 0.5462052), 
+    tolerance = 0.0001
+  )
   expect_equal(
     as.vector(Calibration(exp3[, , , 1], obs3[, , 1], memb_dim = 'member'))[1:5],
     as.vector(Calibration(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')[, , , 1, 1][1:5]),
     tolerance = 0.0001
   )
+  expect_equal(
+    as.vector(Calibration(exp3[, , , 1], obs3[, , 1], exp_cor3, eval.method = "hindcast-vs-forecast"))[1:5],
+    c(-0.80521700, -0.05578454, 0.56143657, -1.01815286, 0.49089916), 
+    tolerance = 0.0001
+  )
   expect_equal(
     as.vector(Calibration(exp3[, , , 1], obs3[, , 1], exp_cor3, memb_dim = 'member'))[1:5],
     as.vector(Calibration(exp3, obs3, exp_cor3, memb_dim = 'member', dat_dim = 'dataset')[, , , 1, 1][1:5])
@@ -288,6 +329,7 @@ test_that("3. Output checks: dat3", {
 ##############################################
 
 test_that("4. Output checks: dat3", {
+  # Check ca.method
   expect_equal(
     as.vector(Calibration(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', cal.method = 'bias'))[1:7],
     c(-0.3984805, 0.4116167, -0.6076553, 1.8232541, 0.5574811, -0.5924950, 0.7154024),
@@ -313,10 +355,10 @@ test_that("4. Output checks: dat3", {
     c(0.4178384, 0.2323757, 0.2942796, 1.4487467, -1.1659159, 1.6175277),
     tolerance = 0.0001
   )
-
 })
 
 ##############################################
+
 test_that("5. Output checks: dat4", {
   expect_equal(
     dim(Calibration(exp4, obs4)),
@@ -332,15 +374,13 @@ test_that("5. Output checks: dat4", {
     c(-0.88657225, -0.08128503, 0.58193721, -1.11537810, 0.50614269),
     tolerance = 0.0001
   )
+
   expect_equal(
-    dim(Calibration(exp4_1, obs4_1, dat_dim = 'dataset')),
-    c(member = 10, sdate = 20, nexp = 1, nobs = 1)
-  )
-  expect_equal(
-    as.vector(Calibration(exp4_1, obs4_1, dat_dim = 'dataset')),
-    as.vector(Calibration(exp4, obs4)),
+    as.vector(Calibration(exp4_1, obs4_1))[1:6],
+    c(-0.6372505, 0.3261267, -0.8860036, 2.0048627, 0.4995904, -0.8679749),
     tolerance = 0.0001
   )
+  # exp_cor4 doesn't have dat_dim
   expect_equal(
     as.vector(Calibration(exp4_1, obs4_1, exp_cor4, dat_dim = 'dataset')),
     as.vector(Calibration(exp4, obs4, exp_cor4)),
@@ -366,9 +406,21 @@ test_that("5. Output checks: dat4", {
     c(0.4583975, 0.2645130, 0.3292279, 1.5361189, -1.1972745, 1.7125642),
     tolerance = 0.0001
   )
+  expect_equal(
+    dim(Calibration(exp4_1, obs4_1, dat_dim = 'dataset')),
+    c(member = 10, sdate = 20, nexp = 1, nobs = 1)
+  )
+  # exp_cor5 with dat_dim
+  expect_equal(
+    dim(Calibration(exp3, obs3, exp_cor5, dat_dim = 'dataset')),
+    c(member = 10, sdate = 10, lat = 2, nexp = 2, nobs = 3)
+  ) 
 })
+
 ##############################################
+
 test_that("6. Output checks: dat4", {
+  # dat_dim = NULL
   expect_equal(
     as.vector(Calibration(exp4, obs4, cal.method = 'bias'))[1:7],
     c(-0.4039514, 0.4061457, -0.6131262, 1.8177832, 0.5520101, -0.5979660, 0.7099314),
@@ -407,9 +459,67 @@ test_that("6. Output checks: dat4", {
 })
 
 ##############################################
-test_that("6. Output checks: dat4", {
-  expect_equal(
-    dim(Calibration(exp3, obs3, exp_cor5, dat_dim = 'dataset')),
-    c(member = 10, sdate = 10, lat = 2, nexp = 2, nobs = 3)
+
+test_that("7. Output checks: dat4", {
+  expect_warning(
+    Calibration(exp = array(1, dim = c(member = 1, sdate = 1)), 
+                obs = array(1, dim = c(sdate = 1))), 
+    "Some forecast data could not be corrected due to data lack and is replaced with NA values."
+  )
+  expect_warning(
+    Calibration(exp = array(1, dim = c(member = 1, sdate = 1)), 
+                obs = array(1, dim = c(sdate = 1)), na.fill = FALSE), 
+    "Some forecast data could not be corrected due to data lack and is replaced with uncorrected values."
+  )
+  suppressWarnings(
+    expect_equal(
+      dim(Calibration(exp = array(1, dim = c(dataset = 1, member = 1, sdate = 1)), 
+                      obs = array(1, dim = c(dataset = 1, sdate = 1)), 
+                      dat_dim = 'dataset', na.fill = FALSE)), 
+      c(member = 1, sdate = 1, nexp = 1, nobs = 1)
+    )
+  )
+  suppressWarnings(
+    expect_equal(
+      dim(Calibration(exp = array(1, dim = c(dataset = 3, member = 1, sdate = 1)), 
+                      obs = array(1, dim = c(dataset = 2, sdate = 1)), 
+                      exp_cor = array(1, dim = c(dataset = 3, member = 1, sdate = 1)),
+                      dat_dim = 'dataset', na.fill = FALSE)), 
+      c(member = 1, sdate = 1, nexp = 3, nobs = 2)
+    )
+  )
+  suppressWarnings(
+    expect_equal(
+      dim(Calibration(exp = array(1, dim = c(dataset = 1, member = 1, sdate = 1)), 
+                      obs = array(1, dim = c(dataset = 2, sdate = 1)), 
+                      exp_cor = array(1, dim = c(member = 1, sdate = 1)),
+                      dat_dim = 'dataset', na.fill = FALSE)), 
+      c(member = 1, sdate = 1, nexp = 1, nobs = 2)
+    )
+  )
+  # Check values
+  suppressWarnings(
+    expect_equal(
+      Calibration(exp = array(1:6, dim = c(member = 2, sdate = 3)), 
+                  obs = array(1:3, dim = c(sdate = 3)), 
+                  na.fill = FALSE), 
+      array(1:6, dim = c(member = 2, sdate = 3))
+    )
+  )
+  suppressWarnings(
+    expect_equal(
+      Calibration(exp = array(1:6, dim = c(member = 2, sdate = 3)), 
+                  obs = array(1:3, dim = c(sdate = 3))), 
+      array(dim = c(member = 2, sdate = 3))
+    )
+  )
+  suppressWarnings(
+    expect_equal(
+      Calibration(exp = array(1:6, dim = c(member = 2, sdate = 3)), 
+                  obs = array(1:3, dim = c(sdate = 3)), 
+                  exp_cor = array(6:12, dim = c(member = 2, sdate = 3)),
+                  na.fill = FALSE), 
+      array(1:6, dim = c(member = 2, sdate = 3))
+    )
   )
 })