diff --git a/R/CRPSS.R b/R/CRPSS.R
index a6b4a1405a80156149cf16cad5c955b9135428c2..cfd8fdfdd7b80ae60b527822c7f164f904e66763 100644
--- a/R/CRPSS.R
+++ b/R/CRPSS.R
@@ -222,11 +222,20 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member',
   if (is.null(ref)) {
     ## using climatology as reference forecast
     ## all the time steps are used as if they were members
-    ## then, ref dimensions are [sdate, memb], both with length(sdate)
+    ## then, ref dimensions are [sdate, memb]
+    ## memb dimension has length(sdate) - 1 due to cross-validation
 
     obs_time_len <- dim(obs)[time_dim]
     if (is.null(dat_dim)) {
-      ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len))
+      
+      ## Without cross-validation: 
+      ## ref <- array(data = rep(obs, each = obs_time_len), dim = c(obs_time_len, obs_time_len))
+      ## 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 <- .CRPS(exp = ref, obs = obs, time_dim = time_dim, memb_dim = memb_dim,
@@ -237,7 +246,15 @@ CRPSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = 'member',
       crps_ref <- array(dim = c(obs_time_len, nobs))
       names(dim(crps_ref)) <- c(time_dim, 'nobs')
       for (i_obs in 1:nobs) {
-        ref <- array(data = rep(obs[, i_obs], each = obs_time_len), dim = c(obs_time_len, obs_time_len))
+        
+        ## Without cross-validation: 
+        ## ref <- array(data = rep(obs[, i_obs], each = obs_time_len), dim = c(obs_time_len, obs_time_len))
+        ## 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, i_obs]
+        }
+        
         names(dim(ref)) <- c(time_dim, memb_dim)
         crps_ref[, i_obs] <- .CRPS(exp = ref, obs = ClimProjDiags::Subset(obs, dat_dim, i_obs, drop = 'selected'), 
                                    time_dim = time_dim, memb_dim = memb_dim, dat_dim = NULL, Fair = Fair)
diff --git a/tests/testthat/test-CRPSS.R b/tests/testthat/test-CRPSS.R
index 6937edb2332888a8368663a1660c89db944e6749..db0eecd632c1494a332987d6bc0f31f079ae7ed3 100644
--- a/tests/testthat/test-CRPSS.R
+++ b/tests/testthat/test-CRPSS.R
@@ -141,7 +141,7 @@ test_that("2. Output checks: dat1", {
   # ref = NULL
   expect_equal(
     as.vector(CRPSS(exp1, obs1)$crpss),
-    c(-0.1582765, -0.2390707),
+    c(0.061796021, -0.003647287),
     tolerance = 0.0001
   )
   expect_equal(
@@ -150,12 +150,12 @@ test_that("2. Output checks: dat1", {
   )
   expect_equal(
     as.vector(CRPSS(exp1, obs1, Fair = T)$crpss),
-    c(0.07650872, -0.09326681),
+    c(0.2612070, 0.1253865),
     tolerance = 0.0001
   )
   expect_equal(
     as.vector(CRPSS(exp1, obs1)$crpss),
-    c(-0.1582765, -0.2390707),
+    c(0.061796021, -0.003647287),
     tolerance = 0.0001
   )
   # ref = ref
@@ -175,7 +175,7 @@ test_that("2. Output checks: dat1", {
   )
   expect_equal(
     as.vector(CRPSS(exp1, obs1, ref1)$crpss),
-    c( 0.3491793, 0.3379610),
+    c(0.3491793, 0.3379610),
     tolerance = 0.0001
   )
 
@@ -211,7 +211,7 @@ test_that("3. Output checks: dat2", {
   # ref = NULL
   expect_equal(
     as.vector(CRPSS(exp2, obs2)$crpss),
-    c(-0.8209236),
+    -0.4749481,
     tolerance = 0.0001
   )
   expect_equal(
@@ -222,11 +222,11 @@ test_that("3. Output checks: dat2", {
 
   expect_equal(
     as.vector(CRPSS(exp2, obs2)$sign),
-    TRUE,
+    FALSE,
   )
   expect_equal(
     as.vector(CRPSS(exp2, obs2, Fair = T)$crpss),
-    c(-0.468189),
+    c(-0.1745512),
     tolerance = 0.0001
   )
   expect_equal(
@@ -265,12 +265,12 @@ test_that("4. Output checks: dat3", {
   )
   expect_equal(
     mean(CRPSS(exp3, obs3, dat_dim = 'dataset')$crpss),
-    c(-0.7390546),
+    c(-0.4086342),
     tolerance = 0.0001
   )
   expect_equal(
     as.vector(CRPSS(exp3, obs3, dat_dim = 'dataset')$crpss),
-    c(-0.8209236, -0.6270744, -1.0829403, -0.6574485, -0.5970569, -0.6488837),
+    c(-0.4749481, -0.3179303, -0.6871816, -0.3425333, -0.2936161, -0.3355958),
     tolerance = 0.0001
   )
   expect_equal(
@@ -279,12 +279,12 @@ test_that("4. Output checks: dat3", {
   )
   expect_equal(
     mean(CRPSS(exp3, obs3, dat_dim = 'dataset', Fair = T)$crpss),
-    c(-0.5302703),
+    -0.2242162,
     tolerance = 0.0001
   )
   expect_equal(
     as.vector(CRPSS(exp3, obs3, dat_dim = 'dataset', Fair = T)$crpss),
-    c(-0.4681890, -0.3580685, -1.0116628, -0.3730576, -0.3965618, -0.5740819),
+    c(-0.17455119, -0.08645477, -0.60933021, -0.09844606, -0.11724947, -0.25926556),
     tolerance = 0.0001
   )
   # ref = ref3
@@ -310,7 +310,7 @@ test_that("5. Output checks: dat4", {
   )
   expect_equal(
     as.vector(CRPSS(exp4, obs4, dat_dim = 'dataset', Fair = T)$crpss)[1:6],
-    c(-0.4681890, -0.3580685, -1.0116628, -0.3730576, -0.3965618, -0.5740819),
+    c(-0.17455119, -0.08645477, -0.60933021, -0.09844606, -0.11724947, -0.25926556),
     tolerance = 0.0001
   )
   # ref = ref3