diff --git a/modules/Downscaling/tmp/Intbc.R b/modules/Downscaling/tmp/Intbc.R
index 6237840472b57c4d9e5ba552302674844b6e32a7..d11b5e1ac9e6ccfc4f83cfc97e3befdcffd75d67 100644
--- a/modules/Downscaling/tmp/Intbc.R
+++ b/modules/Downscaling/tmp/Intbc.R
@@ -73,7 +73,8 @@
 
 CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_method = NULL, 
                       points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", 
-                      sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = NULL, ...)
+                      sdate_dim = "sdate", member_dim = "member", time_dim = "time",
+                      region = NULL, ncores = NULL, ...)
 {
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -91,7 +92,8 @@ CST_Intbc <- function(exp, obs, exp_cor = NULL, target_grid, bc_method, int_meth
                source_file_exp = exp$attrs$source_files[1], 
                source_file_obs = obs$attrs$source_files[1], 
                method_point_interp = method_point_interp, lat_dim = lat_dim, lon_dim = lon_dim, 
-               sdate_dim = sdate_dim, member_dim = member_dim, region = region, ncores = ncores, ...)
+               sdate_dim = sdate_dim, member_dim = member_dim, 
+               time_dim = time_dim, region = region, ncores = ncores, ...)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   obs$data <- res$obs
@@ -363,46 +365,106 @@ Intbc <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
   #    which(names(dim(obs_ref)) == sdate_dim), 'sdate')
   #}
 
-  if (bc_method == 'qm' | bc_method == 'quantile_mapping') {
-   
-    res <- QuantileMapping(exp = exp_interpolated$data, obs = obs_ref, 
-                           exp_cor = exp_cor_interpolated$data,
-                           na.rm = TRUE,  memb_dim = member_dim, sdate_dim = sdate_dim, 
-                           ncores = ncores, ...)
+
+  k_out <- loop <- 1 # the data is NOT hourly/daily/weekly, i.e., time = 1. Apply leave one data out.
+  daily <- FALSE # for the case at which time > 1 or for the case of forecast dwn (exp_cor) 
+  eval.method <- "leave-one-out"
+  # if time > 1, apply leave-one-year-out
+  if (  time_dim %in% names(dim(obs_ref)) & dim(obs_ref)[time_dim] > 1 & is.null(exp_cor)) {
+    daily <- TRUE
+    eval.method <- "in-sample"
+    sdate_num <- as.numeric (dim(obs_ref)[sdate_dim])
+    k_out <- as.numeric (dim(obs_ref)[time_dim]) ## in case loocv = TRUE and the data is daily, leave one-year data out
+    obs_ref <-  MergeDims(obs_ref, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
+    exp_interpolated$data <- MergeDims(exp_interpolated$data, 
+                                       merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
+    loop <- dim(obs_ref)[sdate_dim]/k_out
+  } else if (!is.null(exp_cor)) {
+    eval.method <- "in-sample"
   }
-  else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') {
-    # Dynamical bias correction is not yet prepared to handle hindcast-forecast data
-    # It will return only the hindcast downscaled
-    if (!is.null(exp_cor)) {
-      warning("For the dynamical bias correction only the hindcast downscaled will be returned.")
-    }
-    # the temporal dimension must be only one dimension called "time"
-    if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) {
-      exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), 
-        fun = as.vector, output_dims = "time", ncores = ncores)$output1
-    } 
-    if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) {
-      obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, 
-        output_dims = "time", ncores = ncores)$output1
+
+  sdate_num <- as.numeric (dim(obs_ref)[sdate_dim])
+  res <- list()
+
+  for (yy in 1:loop) {  # with cv purposes. leave-one-year-out
+
+    if (loop == 1) {
+      window <- -(1:sdate_num)
+      EXP_COR <- NULL
+      if (!is.null(exp_cor)) {
+        EXP_COR <- exp_cor_interpolated$data
+      }
+    } else {
+      window <- ((yy-1) * k_out+1):((yy-1) * k_out + k_out)
+      EXP_COR <- ClimProjDiags::Subset(exp_interpolated$data, along = sdate_dim,
+                                       indices = (1:sdate_num)[window], drop = "selected")
     }
-    # REMEMBER to add na.rm = T in colMeans in .proxiesattractor 
-    res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ncores = ncores, ...)
-  } else {
-    if (dim(exp_interpolated$data)[member_dim] == 1) {
-      stop('Calibration must not be used with only one ensemble member.')
+
+    OBS_REF <- ClimProjDiags::Subset(obs_ref, along = sdate_dim, 
+                                     indices = (1:sdate_num)[-window], drop = "selected")
+    EXP <- ClimProjDiags::Subset(exp_interpolated$data, along = sdate_dim, 
+                                 indices = (1:sdate_num)[-window], drop = "selected")
+
+    if (bc_method == 'qm' | bc_method == 'quantile_mapping') {
+      dum <- QuantileMapping(exp = EXP, obs = OBS_REF, 
+                             exp_cor = EXP_COR,
+                             na.rm = TRUE,  memb_dim = member_dim, sdate_dim = sdate_dim, 
+                             ncores = ncores, eval.method = eval.method, ...)
+
+      res[[yy]] <- dum
+      rm(dum)
     }
-    if (dim(obs_ref)[sdate_dim] == 1) {
-      warning('Simple Bias Correction should not be used with only one observation. Returning NA.')
+    else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') {
+      # Dynamical bias correction is not yet prepared to handle hindcast-forecast data
+      # It will return only the hindcast downscaled
+      # Dynamical bias correction is not yet prepared to handle leave-one-year-out for daily case
+      # It will properly work only with the monthly or coarser temporal scale data. 
+      if (!is.null(exp_cor)) {
+        warning("For the dynamical bias correction only the hindcast downscaled will be returned.")
+      }
+      # the temporal dimension must be only one dimension called "time"
+      if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) {
+        exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), 
+        fun = as.vector, output_dims = "time", ncores = ncores)$output1
+      } 
+      if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) {
+        obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, 
+          output_dims = "time", ncores = ncores)$output1
+      }
+      # REMEMBER to add na.rm = T in colMeans in .proxiesattractor 
+      dum <- DynBiasCorrection(exp = EXP, obs = OBS_REF, ncores = ncores, ...)
+      res[[yy]] <- dum
+      rm(dum)
+    } else {
+      if (dim(exp_interpolated$data)[member_dim] == 1) {
+        stop('Calibration must not be used with only one ensemble member.')
+      }
+      if (dim(obs_ref)[sdate_dim] == 1) {
+        warning('Simple Bias Correction should not be used with only one observation. Returning NA.')
+      }
+      dum <- Calibration(exp = EXP, obs = OBS_REF, 
+                         exp_cor = EXP_COR, 
+                         memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, 
+                         cal.method = bc_method,
+	  	         eval.method = eval.method) 
+      res[[yy]] <- dum
+      rm(dum)
     }
-    res <- CSTools::Calibration(exp = exp_interpolated$data, obs = obs_ref,
-                                exp_cor = exp_cor_interpolated$data, 
-                                memb_dim = member_dim, sdate_dim = sdate_dim,
-                                ncores = ncores, 
-                                cal.method = bc_method) 
   } 
- 
+
+  require (abind)  
+  dimnames <- names(dim(res[[1]]))
+  res <- do.call(abind, c(res, along = which(names(dim(res[[1]])) == "syear")))
+  names(dim(res)) <- dimnames
+
+  if (daily) {
+    res <- SplitDim(res, split_dim = sdate_dim, new_dim_name = time_dim,
+                    indices = rep(1:k_out, sdate_num/k_out))
+    obs_ref <- SplitDim(obs_ref, split_dim = sdate_dim, new_dim_name = time_dim,
+                        indices = rep(1:k_out, sdate_num/k_out))
+  }
+    
   # Return a list of three elements
   res <- list(data = res, obs = obs_ref, lon = exp_interpolated$lon, lat = exp_interpolated$lat)
-
   return(res)
 }
diff --git a/modules/Downscaling/tmp/Interpolation.R b/modules/Downscaling/tmp/Interpolation.R
index 51de9ec0748f9e1b2c575e6ecd3664c14cf03d48..f1bebfa0cc11f5784dc5aee656b0b969e45bca21 100644
--- a/modules/Downscaling/tmp/Interpolation.R
+++ b/modules/Downscaling/tmp/Interpolation.R
@@ -68,7 +68,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
     stop("The name of the latitude/longitude dimensions in 'exp$data' must match the parametres 'lat_dim' and 'lon_dim'")
   }
 
-  res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$attrs[[lon_dim]], 
+  res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], 
                        source_file = exp$attrs$source_files[1], points = points,
                        method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim,
                        lon_dim = lon_dim, region = region, method_point_interp = method_point_interp, ncores = ncores)
diff --git a/modules/Downscaling/tmp/Intlr.R b/modules/Downscaling/tmp/Intlr.R
index beb8d7ee9ca6d8d11f96c025bf8408c42c3b8975..f6d972171b14b152b5adfaaea049ef779175a6e0 100644
--- a/modules/Downscaling/tmp/Intlr.R
+++ b/modules/Downscaling/tmp/Intlr.R
@@ -504,13 +504,16 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
     if (!is.null(exp_cor)) {
       aux_dim <- NULL
       forecast <- exp_cor_interpolated$data
-      target_dims_predictor <- c(sdate_dim, member_dim)
-      target_dims_forecast <- c(sdate_dim, member_dim)
+      if (dim(predictor)[[member_dim]] == dim(forecast)[[member_dim]]) {
+        target_dims_forecast <- c(sdate_dim)
+      } else {
+        target_dims_predictor <- c(sdate_dim, member_dim)
+        target_dims_forecast <- c(sdate_dim, member_dim)
+      }
     } else {
       forecast <- NULL
       target_dims_forecast <- NULL
     }
-
   } 
   
   # (Multi) linear regression with large-scale predictors
@@ -534,8 +537,12 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
       } else {
         forecast <- exp_cor$data
       }
-      target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname, member_dim)
-      target_dims_forecast <- c(sdate_dim, large_scale_predictor_dimname, member_dim)
+      if (dim(predictor)[[member_dim]] == dim(forecast)[[member_dim]]) {
+        target_dims_forecast <- c(sdate_dim, large_scale_predictor_dimname)
+      } else {
+        target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname, member_dim)
+        target_dims_forecast <- c(sdate_dim, large_scale_predictor_dimname, member_dim)
+      }     
     }
   } 
   
@@ -597,27 +604,30 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
     }
 
     target_dims_predictand <- sdate_dim 
+    target_dims_predictor <- c(sdate_dim,'nn')
 
     if (!is.null(exp_cor)) {
-      target_dims_predictor <- c(sdate_dim,'nn', member_dim)
-      target_dims_forecast <- c(sdate_dim,'nn', member_dim)   
-    } else {
-      target_dims_predictor <- c(sdate_dim,'nn')
-    }
-  } else {
+      if (dim(predictor)[[member_dim]] == dim(forecast)[[member_dim]]) {
+        target_dims_forecast <- c(sdate_dim,'nn')   
+      } else {
+        target_dims_predictor <- c(sdate_dim,'nn', member_dim)
+        target_dims_forecast <- c(sdate_dim,'nn', member_dim)
+      }
+    } 
 
+  } else {
     stop(paste0(lr_method, " method is not implemented yet"))
   }
 
   k_out <- 1 ## ##  It is used in constructing window for loocv = TRUE case. k_out is also used to order the dimensions of pred.lm output for both loocv = TRUE and loocv = FALSE cases. In case the data is NOT daily, k_out is 1.
   daily <- FALSE
+  output_dims <- c(time_dim, sdate_dim) # for hindcast dwn 
   if (  time_dim %in% names(dim(predictor)) ) {
     daily <- TRUE
     k_out <- as.numeric (dim(predictor)[time_dim]) ## in case loocv = TRUE and the data is daily, leave one-year data out
     sdate_num <- as.numeric (dim(predictand)[sdate_dim]) ## sdate_num of hindcast
     predictor <- MergeDims (predictor, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
     predictand <- MergeDims (predictand, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
-    output_dims <- c(time_dim, sdate_dim) # for hindcast dwn 
     if (!is.null(exp_cor)) {
       sdate_num_fr <- as.numeric (dim(forecast)[sdate_dim]) ## sdate_num of forecast
       forecast <- MergeDims (forecast, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
@@ -634,8 +644,8 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
                                sdate_dim = sdate_dim, k_out = k_out)$output1
 
     if (daily) {
-      res <- SplitDim(data = res, split_dim = sdate_dim, new_dim_name = time_dim,
-                      indices = rep (1:k_out, sdate_num_fr))
+      res <- CSTools::SplitDim(data = res, split_dim = sdate_dim, new_dim_name = time_dim,
+                               indices = rep (1:k_out, sdate_num_fr))
     }
   } 
   ## case hindcast only
@@ -645,16 +655,20 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
                  fun = .intlr, loocv = loocv, ncores = ncores, sdate_dim = sdate_dim, 
                                k_out = k_out,
                                output_dims = output_dims)$output1
+    if (!daily) {
+      res <- ClimProjDiags::Subset(x = res, along = time_dim, indices = 1, drop = 'selected')
+    }
   }
   
   if (daily) {
-    predictand <- SplitDim(data = predictand, split_dim = sdate_dim, new_dim_name = time_dim,
-                           indices = rep (1:k_out, sdate_num))
+    predictand <- CSTools::SplitDim(data = predictand, split_dim = sdate_dim, 
+                                    new_dim_name = time_dim,
+                                    indices = rep (1:k_out, sdate_num))
   }
 
   # restore ensemble dimension in observations if it existed originally
   if (restore_ens) {
-    predictand <- MergeDims (predictand, merge_dims = c(time_dim, sdate_dim), rename_dim = sdate_dim)
+    predictand <- s2dv::InsertDim(predictand, posdim = 1, lendim = 1, name = member_dim)
   }
   
   # Return a list of three elements
@@ -682,9 +696,9 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
 
   require (plyr)
 
-  if (!is.null(f)) {
+  if (!is.null(f) & any(names(dim(f)) == member_dim)) {
     if (!is.null(aux_dim)) {
-      tmp_df <- data.frame(x = adply(x,.margins = 3, .id = NULL), y = y)
+      tmp_df <- data.frame(x = adply(x,.margins = 3, .id = NULL, .fun = as.matrix), y = y)
     } else {
       tmp_df <- data.frame(x = as.vector(x), y = y)
     }
@@ -692,6 +706,8 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
     tmp_df <- data.frame(x = x, y = y)
   }
 
+  colnames(tmp_df) <- c(paste0("x.",1:(ncol(tmp_df)-1)), "y")
+
   # if the data is all NA, force return return NA
   if (all(is.na(tmp_df)) | (sum(apply(tmp_df, 2, function(x) !all(is.na(x)))) == 1)) {
     if (is.null(f)) {  
@@ -704,98 +720,133 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
       }
     }
   } else {
-    
-    ## apply principle components if method is 9nn
-    ## TODO 9NN MIGHT GIVE ERROR FOR FORECAST CASE. REVIEW AFTER DEVELOPING FORECAST STRATEGY!! 
-    if (any(names(dim(x)) == "nn")) { 
-      PR <- prcomp(~ x.1 + x.2 + x.3 + x.4 + x.5 + x.6 + x.7 + x.8 + x.9, 
-                   data = tmp_df[, -ncol(tmp_df)], 
-                   na.action = na.omit, scale = FALSE) ## prcomp between predictors
-      ## 95% of variablity should be explained by the principle components
-      N_predictors <- which (summary(PR)$importance[ "Cumulative Proportion",] > .95)[1] 
-      PR_mat <- PR$x[, 1:N_predictors, drop = FALSE]
-      if (!is.null(PR$na.action)) { ## if there are NAs, relocate them to the correct order
-        dum <- array(NA, dim = c(nrow(tmp_df), ncol(PR_mat) ))
-        dum [as.numeric(rownames(PR_mat)), ] <- PR_mat
-        PR_mat <- dum
-      }
-      tmp_df <- data.frame(x = PR_mat, y = y)
-      colnames(tmp_df) <- c(paste0("x.",1:(ncol(tmp_df)-1)), "y")
-    } 
+
+    # for now, PCA is TRUE only for nn. in future, it might be true for multiple var options and ls.   
+    if (any(names(dim(x)) == "nn")) {
+      pca <- TRUE
+    } else {
+      pca <- FALSE
+    }
 
     # training
-    lm1 <- .train_lm(df = tmp_df, loocv = loocv, k_out = k_out)  
+    lm1 <- .train_lm(df = tmp_df, loocv = loocv, k_out = k_out, pca = pca)  
  
     # prediction
-    res <- .pred_lm(lm1 = lm1, df = tmp_df, f = f, loocv = loocv, aux_dim = aux_dim, k_out = k_out)
+    res <- .pred_lm(lm1 = lm1, df = tmp_df, f = f, loocv = loocv, aux_dim = aux_dim, 
+                    k_out = k_out, pca = pca)
   }
-
-    return(res)
+  return(res)
 }
 
 #-----------------------------------
 # Function to generate the linear regressions.
 # Returns a list
 #-----------------------------------
-.train_lm <- function(df, loocv, k_out = 1) {
+.train_lm <- function(df, loocv, k_out = 1, pca = FALSE) {
   
   # Remove predictor columns containing only NA's
-  df <- df[ ,apply(as.matrix(df[,colnames(df) != 'y'],nrow(df),ncol(df)-1), 2, 
-                   function(x) !all(is.na(x)))]
-  
-  if (loocv) {
+  df_all<- df[ ,apply(as.matrix(df[,colnames(df) != 'y'],nrow(df),ncol(df)-1), 2, 
+                        function(x) !all(is.na(x)))]
 
-    lm1 <- lapply(1:(nrow(df)/k_out), function(j) {
-      window <- ((j-1)*k_out+1):((j-1)*k_out + k_out) # omit the data of the year including corresponding day
-      if (all(is.na(df[-window ,]$y))) {
+  if (loocv) {
+    lm1 <- lapply(1:(nrow(df_all)/k_out), function(j) {
+      window <- ((j-1)*k_out+1):((j-1)*k_out + k_out) # omit the data of the year including corresponding day  
+      df <- df_all[-window,]
+      ## apply principle components if method is 9nn
+      ## TODO 9NN MIGHT GIVE ERROR FOR FORECAST CASE. REVIEW AFTER DEVELOPING FORECAST STRATEGY!! 
+      if (pca) {
+        PCA <- .pca(dx = df[, -ncol(df)], exp.var = 0.95, center = TRUE, scale = FALSE)
+        df <- data.frame(x = PCA$PR_mat, y = df[,ncol(df)])
+        colnames(df) <- c(paste0("x.",1:(ncol(df)-1)), "y")
+        PCA$rotation <- as.matrix(PCA$rotation)
+        colnames(PCA$rotation) <- c(paste0("x.",1:(ncol(df)-1)))
+      } else {
+        PCA <- NULL
+      }
+      if (all(is.na(df$y))) {
         return(NA)
       } else {
-        return(lm(df[-window,], formula = y ~ .)) 
+        return(list(lm = lm(df, formula = y ~ .), N_predictors = PCA$N_predictors, 
+                    rotation = PCA$rotation)) 
       }})
   } else {
-    
-    lm1 <- ifelse(all(is.na(df$y)), NA, list(lm(data = df, formula = y ~ .)))
+    if (pca) {
+      PCA <- .pca(dx = df_all[,-ncol(df_all)], exp.var = 0.95, center = TRUE, scale = FALSE)
+      df_all <- data.frame(x = PCA$PR_mat, y = df_all[,ncol(df_all)])
+      colnames(df_all) <- c(paste0("x.",1:(ncol(df_all)-1)), "y")
+      PCA$rotation <- as.matrix(PCA$rotation)
+      colnames(PCA$rotation) <- c(paste0("x.",1:(ncol(df_all)-1)))
+    } else {
+      PCA <- NULL
+    }
+    if(all(is.na(df_all$y))) {
+      lm1 <- NA
+    } else { 
+      lm1<- list(lm = lm(data = df_all, formula = y ~ .), N_predictors = PCA$N_predictors, 
+                 rotation = PCA$rotation)
+    }
   }
-  
   return(lm1)
 }
 
 #-----------------------------------
 # Function to apply the linear regressions.
 #-----------------------------------
-.pred_lm <- function(df, lm1, f, loocv, aux_dim, k_out = 1) {
+.pred_lm <- function(df, lm1, f, loocv, aux_dim, k_out = 1, pca = FALSE) {
 
   require (plyr)
-
   
   if (loocv) {
     pred_vals <- sapply(1:length(lm1), function(j) {
-      if (all(is.na(lm1[[j]]))) {
+      if (all(is.na(lm1[[j]]$lm))) {
         return(array(rep(NA, k_out), c(k_out, nrow(df) / k_out)))
       } else {
         window <- ((j-1)*k_out+1):((j-1)*k_out + k_out) # test the daily data of the corresponding year
-        return(predict(lm1[[j]], df[window,]))
+        if (pca) {
+          ## transfer test data to PCA space.
+          pred_data <- as.matrix(df[window,-ncol(df)])%*%lm1[[j]]$rotation
+          pred_data <- as.data.frame (pred_data)
+        } else {
+          pred_data <- df[window,]
+        }
+        return(predict(lm1[[j]]$lm, pred_data))
       }})
     pred_vals <- array (pred_vals, c(k_out, nrow(df) / k_out))
   } else {
-    if (!is.na(lm1)) {
+    # if lm function exist, in other words not NA
+    if (any(!is.na(lm1))) {
       # case to downscale hindcasts
-      if (is.null(f)) { 
-        pred_vals_ls <- lapply(lm1, predict, newdata = df)    
-        pred_vals <- array (unlist(pred_vals_ls), c(k_out, nrow(df) / k_out))
+      if (is.null(f)) {
+        if (pca) {
+          pred_data <- as.matrix(df[,-ncol(df)])%*%lm1$rotation
+          pred_data <- as.data.frame(pred_data)
+        } else {
+          pred_data <- df
+        }
+        pred_vals_ls <- predict(lm1$lm, newdata = pred_data)    
+        pred_vals <- array (pred_vals_ls, c(k_out, nrow(df) / k_out))
       # case to downscale forecasts
       } else {
         if (!is.null(aux_dim)) {
           # 9nn & large-scale
-          fcst_df <- adply(f,.margins = 3, .id = NULL)
+          if (length(dim(f)) == 3) { 
+          # if ens member number is different in hcst and fcst, concatenate members
+            fcst_df <- as.data.frame(matrix(aperm(f, c(1, 3, 2)), 
+                                            nrow = dim(f)[1] * dim(f)[3], ncol = dim(f)[2]))
+          } else {
+            fcst_df <- as.data.frame(f)
+          }
+          if (pca) {
+            ## transfer test data to PCA space.
+            fcst_df <- as.matrix(fcst_df)%*%lm1$rotation
+            fcst_df <- as.data.frame (fcst_df)
+          } 
           colnames(fcst_df) <- paste0('x.',1:ncol(fcst_df))
-          pred_vals_ls <- lapply(lm1, predict, newdata = fcst_df)
-          pred_vals <- unlist(pred_vals_ls)
+          pred_vals <- predict(lm1$lm, newdata = fcst_df)
           dim(pred_vals) <- dim(f)[names(dim(f)) != aux_dim] 
         } else { 
         # basic
-          pred_vals_ls <- lapply(lm1, predict, newdata = data.frame(x = as.vector(f)))
-          pred_vals <- unlist(pred_vals_ls)
+          pred_vals <- predict(lm1$lm, newdata = data.frame(x.1 = as.vector(f)))
           dim(pred_vals) <- dim(f)
         }
       }
@@ -811,7 +862,6 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
       }
     }
   }
- 
   return(pred_vals)
 }
 
@@ -831,9 +881,9 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
 
   names(dim(dist_id)) <- c("nn", lon_dim)
 
-  idh <- SplitDim(dist_id, split_dim = lon_dim,
-                  indices = rep(1:(length(lats_hres)),length(lons_hres)),
-                  new_dim_name = lat_dim)
+  idh <- CSTools::SplitDim(dist_id, split_dim = lon_dim,
+                           indices = rep(1:(length(lats_hres)),length(lons_hres)),
+                           new_dim_name = lat_dim)
 
   idh <- aperm (idh, c(1,3,2))  # dimension: nn, lat_dim, lon_dim
 
@@ -861,7 +911,7 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
                                       "grd" ),   
                    fun = function(x, y) {
                            a <- x[, y[1], y[2], , drop = FALSE]
-                           a <- Subset (a, along = c(lat_dim, lon_dim), 
+                           a <- ClimProjDiags::Subset (a, along = c(lat_dim, lon_dim), 
                                         indices = list (1,1), 
                                         drop = "selected") # drop lat lon dim coming from coar data
                            return(a)
@@ -870,3 +920,19 @@ Intlr <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats, obs_lo
 
   return(nearest)
 }
+
+.pca <- function (dx, exp.var = 0.95, scale = FALSE, center = TRUE, ...) {
+  PR <- suppressWarnings(prcomp (~., dx, center = center,
+                                 na.action = na.omit,
+                                 scale = scale, ...)) ## prcomp between predictors
+
+  N_predictors <- which (summary(PR)$importance[ "Cumulative Proportion",] > exp.var)[1]
+  PR_mat <- PR$x[, 1:N_predictors, drop = FALSE]
+  if (!is.null(PR$na.action)) { ## if there are NAs, relocate them to the correct order
+    dum <- array(NA, dim = c(nrow(dx), ncol(PR_mat) ))
+    dum [as.numeric(rownames(PR_mat)), ] <- PR_mat
+    PR_mat <- dum
+  }
+  return (list(PR_mat = PR_mat, N_predictors = N_predictors, 
+               rotation = PR$rotation[,1:N_predictors]))
+}
diff --git a/modules/Downscaling/tmp/LogisticReg.R b/modules/Downscaling/tmp/LogisticReg.R
index c0d914abd45843b6b23a41dae76b38371635a7ba..2a6259bea18d023bafc773064b5dcae418de658b 100644
--- a/modules/Downscaling/tmp/LogisticReg.R
+++ b/modules/Downscaling/tmp/LogisticReg.R
@@ -484,8 +484,6 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats,
     obs_ref <- obs
   }
 
-  # convert observations to categorical predictands
-
   k_out <- 1 ##  in case loocv = TRUE and the data is NOT daily, leave one data out.
   daily <- FALSE
   if (  time_dim %in% names(dim(obs_ref)) & dim(obs_ref)[time_dim] > 1) {
@@ -500,6 +498,7 @@ LogisticReg <- function(exp, obs, exp_cor = NULL, exp_lats, exp_lons, obs_lats,
     }
   }
   
+  # convert observations to categorical predictands
   obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) {
                      if (!any(!is.na(x))) {
                        rep(NA,length(x))