diff --git a/modules/Skill/R/CRPS_clim.R b/modules/Skill/R/CRPS_clim.R
index 36db4e945fcfc7256ed66c05b98bf9db8c6892e3..b66cab7872bf91102766a9124c82830deba90d24 100644
--- a/modules/Skill/R/CRPS_clim.R
+++ b/modules/Skill/R/CRPS_clim.R
@@ -1,13 +1,27 @@
-CRPS_clim <- function(obs, memb_dim ='ensemble'){  
+# CRPS version for climatology 
+CRPS_clim <- function(obs, memb_dim ='ensemble', return_mean = TRUE, clim.cross.val= TRUE){  
   time_dim <- names(dim(obs))
   obs_time_len <- dim(obs)[time_dim]
-
+  
+  if (isFALSE(clim.cross.val)) { ## Without cross-validation
     ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len))
-    names(dim(ref)) <- c(time_dim, memb_dim)
-    # ref: [sdate, memb]
-    # obs: [sdate]
+  } else if (isTRUE(clim.cross.val)) { ## With cross-validation (excluding the value of that year to create ref for that year)
+    ref <- array(data = NA, dim = c(obs_time_len, obs_time_len - 1))
+    for (i in 1:obs_time_len) {
+      ref[i, ] <- obs[-i]
+    }
+  }
+  
+  names(dim(ref)) <- c(time_dim, memb_dim)
+  # ref: [sdate, memb]
+  # obs: [sdate]
   crps_ref <- s2dv:::.CRPS(exp = ref, obs = obs, time_dim = time_dim, memb_dim = memb_dim,
-                      dat_dim = NULL, Fair = FALSE)
-    # crps_ref should be [sdate]
-  return(mean(crps_ref))
+                           dat_dim = NULL, Fair = FALSE)
+  
+  # crps_ref should be [sdate]
+  if (return_mean == TRUE) {
+    return(mean(crps_ref)) 
+  } else {
+    return(crps_ref)
+  }
 }
diff --git a/modules/Skill/R/RPS_clim.R b/modules/Skill/R/RPS_clim.R
index e8b6452d14f4bfa2972af651b16f56011cf20b13..c93c67476ffcdcbf450f3f2efb37d9f6a2c69b7e 100644
--- a/modules/Skill/R/RPS_clim.R
+++ b/modules/Skill/R/RPS_clim.R
@@ -1,12 +1,12 @@
 # RPS version for climatology 
-RPS_clim <- function(obs, indices_for_clim = NULL, prob_thresholds = c(1/3, 2/3)) { 
+RPS_clim <- function(obs, indices_for_clim = NULL, prob_thresholds = c(1/3, 2/3), cross.val = TRUE) { 
   
   if (is.null(indices_for_clim)){
     indices_for_clim <- 1:length(obs)
   }
   
   obs_probs <- .GetProbs(data = obs, indices_for_quantiles = indices_for_clim, ## temporarily removed s2dv:::
-                                 prob_thresholds = prob_thresholds, weights = NULL, cross.val = T) ## here!
+                                 prob_thresholds = prob_thresholds, weights = NULL, cross.val = cross.val) 
   # clim_probs: [bin, sdate]
   clim_probs <- c(prob_thresholds[1], diff(prob_thresholds), 1 - prob_thresholds[length(prob_thresholds)])
   clim_probs <- array(clim_probs, dim = dim(obs_probs))
diff --git a/modules/Skill/Skill.R b/modules/Skill/Skill.R
index fd54f17f3471cd0ca880712eee91f64f304666cf..04fbc52efb5246e86c601007be1ced893439fc14 100644
--- a/modules/Skill/Skill.R
+++ b/modules/Skill/Skill.R
@@ -89,6 +89,7 @@ Skill <- function(recipe, data, agg = 'global') {
       skill_metrics[[ metric ]] <- skill
       rps_clim <- Apply(list(data$obs$data),
                              target_dims = c(time_dim, memb_dim),
+                             cross.val = cross.val,
                              RPS_clim)$output1
       rps_clim <- .drop_dims(rps_clim)
       skill_metrics[[paste0(metric, "_clim")]] <- rps_clim
@@ -140,7 +141,8 @@ Skill <- function(recipe, data, agg = 'global') {
       skill <- .drop_dims(skill)
       skill_metrics[[ metric ]] <- skill
       crps_clim <- Apply(list(data$obs$data), target_dims = time_dim,
-                         fun = CRPS_clim, memb_dim = memb_dim)$output1
+                         fun = CRPS_clim, memb_dim = memb_dim,
+                         clim.cross.val = cross.val)$output1
       crps_clim <- .drop_dims(crps_clim)
       skill_metrics[['crps_clim']] <- crps_clim
     # CRPSS and FCRPSS