diff --git a/R/RMSSS.R b/R/RMSSS.R
index 3d1de43df26095c4195694c443000281be7fc036..cbd8aa8d5babadd5e323525faeb6733e6f7b6a6d 100644
--- a/R/RMSSS.R
+++ b/R/RMSSS.R
@@ -132,27 +132,6 @@ RMSSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL,
       any(is.null(names(dim(obs)))) | any(nchar(names(dim(obs))) == 0)) {
     stop("Parameter 'exp' and 'obs' must have dimension names.")
   }
-  if (!is.null(ref)) {
-    if (!is.numeric(ref)) {
-      stop("Parameter 'ref' must be numeric.")
-    }
-    if (is.array(ref)) {
-      if (any(is.null(names(dim(ref)))) | any(nchar(names(dim(ref))) == 0)) {
-        stop("Parameter 'ref' must have dimension names.")
-      }
-    } else if (length(ref) != 1 | !all(ref %in% c(0, 1))) {
-      stop("Parameter 'ref' must be a numeric array or number 0 or 1.")
-    }
-  } else {
-    ref <- 0
-    .warning("If a reference dataset is not provided (ref = NULL), the default ", 
-            "value for the climatology is 0 and RMSSS results will only be ",
-            "correct if 'exp' and 'obs' are anomalies. Provide a non-null ",
-            "'ref' for full-field data.")
-  }
-  if (!is.array(ref)) { # 0 or 1
-    ref <- array(data = ref, dim = dim(exp))
-  }
 
   ## time_dim
   if (!is.character(time_dim) | length(time_dim) > 1) {
@@ -264,22 +243,45 @@ RMSSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL,
     stop("Parameter 'exp' and 'obs' must have same length of ",
          "all dimensions except 'dat_dim' and 'memb_dim'.")
   }
-  
+  ## ref
+  if (!is.null(ref)) {
+    if (!is.numeric(ref)) {
+      stop("Parameter 'ref' must be numeric.")
+    }
+    if (is.array(ref)) {
+      if (any(is.null(names(dim(ref)))) | any(nchar(names(dim(ref))) == 0)) {
+        stop("Parameter 'ref' must have dimension names.")
+      }
+    } else {
+      if (length(ref) != 1) {
+        stop("Parameter 'ref' must be a numeric array or a single numeric value.")
+      }
+      ref <- array(data = ref, dim = dim(exp)) # 'ref' is a single numeric value
+    }
+  } else {
+    ref <- MeanDims(data = obs, dims = time_dim, na.rm = TRUE, drop = TRUE)
+    if (!is.array(ref)) {
+      ref <- array(data = ref, dim = dim(exp))
+    } else {
+      ref <- InsertDim(data = ref, posdim = which(names(dim(exp)) == time_dim),
+                       lendim = as.numeric(dim(exp)[time_dim]), name = time_dim)
+    }
+  }
   name_ref <- sort(names(dim(ref)))
   if (!is.null(memb_dim) && memb_dim %in% name_ref) {
     name_ref <- name_ref[-which(name_ref == memb_dim)]
   }
   if (!is.null(dat_dim)) {
     if (dat_dim %in% name_ref) {
-      if (!identical(dim(exp)[dat_dim], dim(ref)[dat_dim])) {
+      if (!identical(dim(obs)[dat_dim], dim(ref)[dat_dim])) {
         stop("If parameter 'ref' has dataset dimension, it must be ", 
-             "equal to dataset dimension of 'exp'.")
+             "equal to dataset dimension of 'obs'.")
       }
       name_ref <- name_ref[-which(name_ref == dat_dim)]
     }
   }
-  if (!identical(length(name_exp), length(name_ref)) |
-      !identical(dim(exp)[name_exp], dim(ref)[name_ref])) {
+  if (!identical(length(name_exp), length(name_ref)) | ## obs?
+      !identical(dim(exp)[name_exp], dim(ref)[name_ref])) { ## obs?
     stop("Parameter 'exp' and 'ref' must have the same length of ",
          "all dimensions except 'memb_dim' and 'dat_dim' if there is ",
          "only one reference dataset.")
@@ -360,8 +362,8 @@ RMSSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', dat_dim = NULL,
   # ref: [sdate, (dat)] or NULL
 
   if (is.null(ref)) {
-    ref <- array(data = 0, dim = dim(obs))
-  } else if (identical(ref, 0) | identical(ref, 1)) {
+    ref <- array(data = mean(obs), dim = dim(obs))
+  } else if (is.numeric(ref) & length(ref)==1) {
     ref <- array(ref, dim = dim(exp))
   }
 
diff --git a/tests/testthat/test-RMSSS.R b/tests/testthat/test-RMSSS.R
index a364b40c50cd917002560c6e0b0613e7e3091a4a..e11f5b0b8bcebb8bde2075cc6e0fa6b6c4379aab 100644
--- a/tests/testthat/test-RMSSS.R
+++ b/tests/testthat/test-RMSSS.R
@@ -5,7 +5,7 @@
  set.seed(2)
  obs1 <- array(rnorm(6), dim = c(sdate = 3, dataset = 2))
  set.seed(3)
- ref1_1 <- array(rnorm(15), dim = c(sdate = 3, dataset = 5))
+ ref1_1 <- array(rnorm(6), dim = c(sdate = 3, dataset = 2))
  ref1_2 <- exp1[, 3]
  dim(ref1_2) <- c(sdate = 3)
 
@@ -116,12 +116,12 @@ test_that("2. Output checks: case 1", {
   )
   expect_equal(
   c(res1_1$rmsss)[3:8],
-  c(0.01693900, -0.02730428, -0.34167869, 0.01588531, -1.55611403, -0.63596896),
+  c(-0.02274264, -0.27334497, -0.39583607, 0.46331207, -2.05011175, 0.10782269),
     tolerance = 0.00001
   )
   expect_equal(
   as.vector(res1_1$p.val)[3:7],
-  c(0.4914588, 0.5134658, 0.6428701, 0.4919943, 0.8672634),
+  c(0.5112420, 0.6185255, 0.6608282, 0.2236229, 0.9029426),
   tolerance = 0.001
   )
 
@@ -135,17 +135,17 @@ test_that("2. Output checks: case 1", {
   )
   expect_equal(
   c(res1_2$rmsss),
-  c(rep(NA, 7), -0.6359690, -0.5877153, -1.1108657),
+  c(rep(NA, 8), -0.8945591, NA),
   tolerance = 0.0001
   )
   expect_equal(
   c(res1_2$p.val),
-  c(rep(NA, 7), 0.7279944, 0.7159769, 0.8167073),
+  c(rep(NA, 8), 0.7821044, NA),
   tolerance = 0.0001
   )
   expect_equal(
   c(res1_2$sign),
-  c(rep(NA, 7), rep(FALSE, 3))
+  c(rep(NA, 8), FALSE, NA)
   )
 
   #ref
@@ -160,12 +160,12 @@ test_that("2. Output checks: case 1", {
   )
   expect_equal(
   as.vector(res1_3$rmsss[2, ]),
-  c(-1.197432, -8.879809),
+  c(-10.600974, -8.879809),
   tolerance = 0.0001
   )
   expect_equal(
   as.vector(res1_3$p.val[2, ]),
-  c(0.8284354, 0.9898591),
+  c(0.9926244, 0.9898591),
   tolerance = 0.0001
   )
   res1_4 <- RMSSS(exp1, obs1, ref = ref1_2, dat_dim = 'dataset', sign = T, alpha = 0.3)
@@ -275,7 +275,7 @@ test_that("4. Output checks: case 3", {
   )
   expect_equal(
   as.vector(RMSSS(exp3, obs3)$rmsss),
-  -0.6289915,
+  -0.6408554,
   tolerance = 0.00001
   )