diff --git a/R/Corr.R b/R/Corr.R
index 38d3901529d185a92e53b9ec7e944728df6427d9..a74725f14be1990f9bb5352201f5eb3aa8d936f3 100644
--- a/R/Corr.R
+++ b/R/Corr.R
@@ -10,15 +10,20 @@
 #'available at all leadtimes.\cr 
 #'The confidence interval is computed by the Fisher transformation and the 
 #'significance level relies on an one-sided student-T distribution.\cr 
+#'If the dataset has more than one member, ensemble mean is necessary necessary
+#'before using this function since it only allows one dimension 'dat_dim' to
+#'have inconsistent length between 'exp' and 'obs'. If all the dimensions of
+#''exp' and 'obs' are identical, you can simply use apply() and cor() to 
+#'compute the correlation.
 #'
 #'@param exp A named numeric array of experimental data, with at least two 
-#'  dimensions 'time_dim' and 'memb_dim'.
+#'  dimensions 'time_dim' and 'dat_dim'.
 #'@param obs A named numeric array of observational data, same dimensions as  
-#'  parameter 'exp' except along memb_dim.
+#'  parameter 'exp' except along dat_dim.
 #'@param time_dim A character string indicating the name of dimension along  
 #'  which the correlations are computed. The default value is 'sdate'.
-#'@param memb_dim A character string indicating the name of member (nobs/nexp) 
-#'  dimension. The default value is 'member'.
+#'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
+#'  dimension. The default value is 'dataset'.
 #'@param comp_dim A character string indicating the name of dimension along which
 #'  obs is taken into account only if it is complete. The default value
 #'  is NULL.
@@ -38,8 +43,8 @@
 #'@return 
 #'A list containing the numeric arrays with dimension:\cr 
 #'  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-#'nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-#'number of observation (i.e., memb_dim in obs).\cr
+#'nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+#'number of observation (i.e., dat_dim in obs).\cr
 #'\item{$corr}{
 #'  The correlation coefficient. 
 #'}
@@ -57,7 +62,7 @@
 #'# Load sample data as in Load() example: 
 #'example(Load) 
 #'clim <- Clim(sampleData$mod, sampleData$obs) 
-#'corr <- Corr(clim$clim_exp, clim$clim_obs, time_dim = 'ftime')
+#'corr <- Corr(clim$clim_exp, clim$clim_obs, time_dim = 'ftime', dat_dim = 'member')
 #'# Renew the example when Ano and Smoothing is ready
 #'
 #'@rdname Corr
@@ -65,7 +70,7 @@
 #'@importFrom ClimProjDiags Subset
 #'@importFrom stats cor pt qnorm 
 #'@export
-Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', 
+Corr <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset', 
                  comp_dim = NULL, limits = NULL,
                  method = 'pearson', pval = TRUE, conf = TRUE,
                  conf.lev = 0.95, ncores = NULL) {
@@ -80,7 +85,7 @@ Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   }
   if (is.null(dim(exp)) | is.null(dim(obs))) {
     stop(paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
-                "containing time_dim and memb_dim."))
+                "containing time_dim and dat_dim."))
   }
   if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
      any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
@@ -97,12 +102,12 @@ Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
     stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
   }
-  ## memb_dim
-  if (!is.character(memb_dim) | length(memb_dim) > 1) {
-    stop("Parameter 'memb_dim' must be a character string.")
+  ## dat_dim
+  if (!is.character(dat_dim) | length(dat_dim) > 1) {
+    stop("Parameter 'dat_dim' must be a character string.")
   }
-  if (!memb_dim %in% names(dim(exp)) | !memb_dim %in% names(dim(obs))) {
-    stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
+  if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
   }
   ## comp_dim
   if (!is.null(comp_dim)) {
@@ -150,11 +155,11 @@ Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   ## exp and obs (2)
   name_exp <- sort(names(dim(exp)))
   name_obs <- sort(names(dim(obs)))
-  name_exp <- name_exp[-which(name_exp == memb_dim)]
-  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  name_exp <- name_exp[-which(name_exp == dat_dim)]
+  name_obs <- name_obs[-which(name_obs == dat_dim)]
   if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
     stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
-                "all dimension expect 'memb_dim'."))
+                "all dimension expect 'dat_dim'."))
   }
   if (dim(exp)[time_dim] < 3) {
     stop("The length of time_dim must be at least 3 to compute correlation.")
@@ -185,8 +190,8 @@ Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   }
 
  res <- Apply(list(exp, obs), 
-              target_dims = list(c(time_dim, memb_dim), 
-                                 c(time_dim, memb_dim)),
+              target_dims = list(c(time_dim, dat_dim), 
+                                 c(time_dim, dat_dim)),
               fun = .Corr, 
               time_dim = time_dim, method = method,
               pval = pval, conf = conf, conf.lev = conf.lev, 
@@ -197,19 +202,19 @@ Corr <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
 .Corr <- function(exp, obs, time_dim = 'sdate', method = 'pearson',
                   conf = TRUE, pval = TRUE, conf.lev = 0.95) {
 
-  # exp: [sdate, member_exp]
-  # obs: [sdate, member_obs]
-  n_exp <- as.numeric(dim(exp)[2])
-  n_obs <- as.numeric(dim(obs)[2])
+  # exp: [sdate, dat_exp]
+  # obs: [sdate, dat_obs]
+  nexp <- as.numeric(dim(exp)[2])
+  nobs <- as.numeric(dim(obs)[2])
   
-  CORR <- array(dim = c(n_exp = n_exp, n_obs = n_obs))
-  eno_expand <- array(dim = c(n_exp = n_exp, n_obs = n_obs))
-  p.val <- array(dim = c(n_exp = n_exp, n_obs = n_obs))
+  CORR <- array(dim = c(nexp = nexp, nobs = nobs))
+  eno_expand <- array(dim = c(nexp = nexp, nobs = nobs))
+  p.val <- array(dim = c(nexp = nexp, nobs = nobs))
 
   # ens_mean
-  for (i in 1:n_obs) {
+  for (i in 1:nobs) {
 
-    CORR[, i] <- sapply(1:n_exp, 
+    CORR[, i] <- sapply(1:nexp, 
                           function(x) {
     if (any(!is.na(exp[, x])) && sum(!is.na(obs[, i])) > 2) {
 cor(exp[, x], obs[, i],
@@ -222,8 +227,8 @@ cor(exp[, x], obs[, i],
   }
 
 #  if (pval) {
-#    for (i in 1:n_obs) {
-#      p.val[, i] <- try(sapply(1:n_exp,
+#    for (i in 1:nobs) {
+#      p.val[, i] <- try(sapply(1:nexp,
 #                           function(x) {(cor.test(exp[, x], obs[, i],
 #                                         use = "pairwise.complete.obs",
 #                                         method = method)$p.value)/2}), silent = TRUE)
@@ -241,7 +246,7 @@ cor(exp[, x], obs[, i],
     } else if (method == "pearson") {
       eno <- Eno(obs, time_dim)  
     }
-    for (i in 1:n_exp) {
+    for (i in 1:nexp) {
       eno_expand[i, ] <- eno
     }
   }
diff --git a/R/GMST.R b/R/GMST.R
index 5e6f4af1b83daca1d4f14dfb566e0af5f9b927bb..8f6f6e4acaa1149977a3595cf35ea39b3dd946be 100644
--- a/R/GMST.R
+++ b/R/GMST.R
@@ -79,7 +79,9 @@
 #' sea_value <- 1
 #' lat <- seq(-90, 90, 10)
 #' lon <- seq(0, 360, 10)
-#' index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value)
+#' index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat,
+#'                   data_lons = lon, type = 'obs', 
+#'                   mask_sea_land = mask_sea_land, sea_value = sea_value)
 #' 
 #' ## Historical simulations
 #' hist_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
@@ -88,7 +90,9 @@
 #' sea_value <- 1
 #' lat <- seq(-90, 90, 10)
 #' lon <- seq(0, 360, 10)
-#' index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value)
+#' index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, 
+#'                    data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, 
+#'                    sea_value = sea_value)
 #' 
 #' ## Decadal predictions
 #' dcpp_tas <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
@@ -97,7 +101,9 @@
 #' sea_value <- 1
 #' lat <- seq(-90, 90, 10)
 #' lon <- seq(0, 360, 10)
-#' index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value)
+#' index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, 
+#'                    data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land,
+#'                    sea_value = sea_value)
 #'
 #'@author Carlos Delgado-Torres, \email{carlos.delgado@bsc.es}
 #'@author Roberto Bilbao, \email{roberto.bilbao@bsc.es}
diff --git a/R/RMS.R b/R/RMS.R
index e9038edc1dc78293d955f359864daf69f87c4c98..c059c8f7ba36bb83d1915c171f6f3855b4c5d4b2 100644
--- a/R/RMS.R
+++ b/R/RMS.R
@@ -11,13 +11,13 @@
 #'The confidence interval is computed by the chi2 distribution.\cr
 #'
 #'@param exp A named numeric array of experimental data, with at least two 
-#'  dimensions 'time_dim' and 'memb_dim'.
+#'  dimensions 'time_dim' and 'dat_dim'.
 #'@param obs A named numeric array of observational data, same dimensions as  
-#'  parameter 'exp' except along memb_dim.
+#'  parameter 'exp' except along dat_dim.
 #'@param time_dim A character string indicating the name of dimension along  
 #'  which the correlations are computed. The default value is 'sdate'.
-#'@param memb_dim A character string indicating the name of member (nobs/nexp) 
-#'  dimension. The default value is 'member'.
+#'@param dat_dim A character string indicating the name of member (nobs/nexp) 
+#'  dimension. The default value is 'dataset'.
 #'@param comp_dim A character string indicating the name of dimension along which
 #'  obs is taken into account only if it is complete. The default value
 #'  is NULL.
@@ -33,8 +33,8 @@
 #'@return 
 #'A list containing the numeric arrays with dimension:\cr 
 #'  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-#'nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-#'number of observation (i.e., memb_dim in obs).\cr
+#'nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+#'number of observation (i.e., dat_dim in obs).\cr
 #'\item{$rms}{
 #'  The root mean square error. 
 #'}
@@ -48,9 +48,9 @@
 #'@examples
 #'# Load sample data as in Load() example:
 #'  set.seed(1)
-#'  exp1 <- array(rnorm(120), dim = c(member = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
+#'  exp1 <- array(rnorm(120), dim = c(dataset = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
 #'  set.seed(2)
-#'  obs1 <- array(rnorm(80),  dim = c(member = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
+#'  obs1 <- array(rnorm(80),  dim = c(dataset = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
 #'  set.seed(2)
 #'  na <- floor(runif(10, min = 1, max = 80))
 #'  obs1[na] <- NA
@@ -62,7 +62,7 @@
 #'@importFrom ClimProjDiags Subset
 #'@importFrom stats qchisq
 #'@export
-RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
+RMS <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
                 comp_dim = NULL, limits = NULL, 
                 conf = TRUE, conf.lev = 0.95, ncores = NULL) {       
   # Check inputs 
@@ -75,7 +75,7 @@ RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   }
   if (is.null(dim(exp)) | is.null(dim(obs))) {
     stop(paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
-                "containing time_dim and memb_dim."))
+                "containing time_dim and dat_dim."))
   }
   if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
      any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
@@ -92,12 +92,12 @@ RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
     stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
   }
-  ## memb_dim
-  if (!is.character(memb_dim) | length(memb_dim) > 1) {
-    stop("Parameter 'memb_dim' must be a character string.")
+  ## dat_dim
+  if (!is.character(dat_dim) | length(dat_dim) > 1) {
+    stop("Parameter 'dat_dim' must be a character string.")
   }
-  if (!memb_dim %in% names(dim(exp)) | !memb_dim %in% names(dim(obs))) {
-    stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
+  if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
   }
   ## comp_dim
   if (!is.null(comp_dim)) {
@@ -137,11 +137,11 @@ RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   ## exp and obs (2)
   name_exp <- sort(names(dim(exp)))
   name_obs <- sort(names(dim(obs)))
-  name_exp <- name_exp[-which(name_exp == memb_dim)]
-  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  name_exp <- name_exp[-which(name_exp == dat_dim)]
+  name_obs <- name_obs[-which(name_obs == dat_dim)]
   if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
     stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
-                "all dimension expect 'memb_dim'."))
+                "all dimension expect 'dat_dim'."))
   }
   if (dim(exp)[time_dim] < 2) {
     stop("The length of time_dim must be at least 2 to compute RMS.")
@@ -172,50 +172,50 @@ RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   }
   
   res <- Apply(list(exp, obs), 
-               target_dims = list(c(time_dim, memb_dim), 
-                                  c(time_dim, memb_dim)),
+               target_dims = list(c(time_dim, dat_dim), 
+                                  c(time_dim, dat_dim)),
                fun = .RMS, 
-               time_dim = time_dim, memb_dim = memb_dim,
+               time_dim = time_dim, dat_dim = dat_dim,
                conf = conf, conf.lev = conf.lev, ncores = ncores)
   return(res)
 }
 
-.RMS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
+.RMS <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset',
                  conf = TRUE, conf.lev = 0.95) { 
 
-  # exp: [sdate, member_exp]
-  # obs: [sdate, member_obs]
-  n_exp <- as.numeric(dim(exp)[2])
-  n_obs <- as.numeric(dim(obs)[2])
-  n_sdate <- as.numeric(dim(exp)[1])
+  # exp: [sdate, dat_exp]
+  # obs: [sdate, dat_obs]
+  nexp <- as.numeric(dim(exp)[2])
+  nobs <- as.numeric(dim(obs)[2])
+  nsdate <- as.numeric(dim(exp)[1])
 
-  dif <- array(dim = c(sdate = n_sdate, n_exp = n_exp, n_obs = n_obs))
-  chi <- array(dim = c(nexp = n_exp, nobs = n_obs))
+  dif <- array(dim = c(sdate = nsdate, nexp = nexp, nobs = nobs))
+  chi <- array(dim = c(nexp = nexp, nobs = nobs))
   if (conf) {
     conflow <- (1 - conf.lev) / 2
     confhigh <- 1 - conflow
-    conf.lower <- array(dim = c(nexp = n_exp, nobs = n_obs))
-    conf.upper <- array(dim = c(nexp = n_exp, nobs = n_obs))
+    conf.lower <- array(dim = c(nexp = nexp, nobs = nobs))
+    conf.upper <- array(dim = c(nexp = nexp, nobs = nobs))
   }
 
   # dif
-  for (i in 1:n_obs) {
-    dif[, , i] <- sapply(1:n_exp, function(x) {exp[, x] - obs[, i]})
+  for (i in 1:nobs) {
+    dif[, , i] <- sapply(1:nexp, function(x) {exp[, x] - obs[, i]})
   }
-  rms <- apply(dif^2, c(2, 3), mean, na.rm = TRUE)^0.5  #array(dim = c(n_exp, n_obs))
+  rms <- apply(dif^2, c(2, 3), mean, na.rm = TRUE)^0.5  #array(dim = c(_exp, nobs))
 
   if (conf) {
-    #eno <- Eno(dif, 1) #count effective sample along sdate. dim = c(n_exp, n_obs)
+    #eno <- Eno(dif, 1) #count effective sample along sdate. dim = c(nexp, nobs)
     eno <- Eno(dif, time_dim) #change to this line when Eno() is done
 
     # conf.lower
-    chi <- sapply(1:n_obs, function(i) {
+    chi <- sapply(1:nobs, function(i) {
                              qchisq(confhigh, eno[, i] - 1)
                            })
     conf.lower <- (eno * rms ** 2 / chi) ** 0.5
 
     # conf.upper
-    chi <- sapply(1:n_obs, function(i) {
+    chi <- sapply(1:nobs, function(i) {
                              qchisq(conflow, eno[, i] - 1)
                            })
     conf.upper <- (eno * rms ** 2 / chi) ** 0.5
diff --git a/R/RMSSS.R b/R/RMSSS.R
index 5967bd00ebcead4a20a80382b445187484058a28..a00606666df797dbb441e87191f7d7d4a0ac12e1 100644
--- a/R/RMSSS.R
+++ b/R/RMSSS.R
@@ -2,7 +2,7 @@
 #'
 #'Compute the root mean square error skill score (RMSSS) between an array of 
 #'forecast 'exp' and an array of observation 'obs'. The two arrays should 
-#'have the same dimensions except along memb_dim, where the length can be 
+#'have the same dimensions except along dat_dim, where the length can be 
 #'different, with the number of experiments/models (nexp) and the number of 
 #'observational datasets (nobs).\cr
 #'RMSSS computes the root mean square error skill score of each jexp in 1:nexp 
@@ -13,17 +13,17 @@
 #'The p-value is optionally provided by an one-sided Fisher test.\cr
 #'
 #'@param exp A named numeric array of experimental data which contains at least
-#'  two dimensions for memb_dim and time_dim. It can also be a vector with the 
+#'  two dimensions for dat_dim and time_dim. It can also be a vector with the 
 #'  same length as 'obs', then the vector will automatically be 'time_dim' and 
-#'  'memb_dim' will be 1.
+#'  'dat_dim' will be 1.
 #'@param obs A named numeric array of observational data which contains at least
-#'  two dimensions for memb_dim and time_dim. The dimensions should be the same 
-#'  as paramter 'exp' except the length of 'memb_dim' dimension. The order of 
+#'  two dimensions for dat_dim and time_dim. The dimensions should be the same 
+#'  as paramter 'exp' except the length of 'dat_dim' dimension. The order of 
 #'  dimension can be different. It can also be a vector with the same length as
-#'  'exp', then the vector will automatically be 'time_dim' and 'memb_dim' will
+#'  'exp', then the vector will automatically be 'time_dim' and 'dat_dim' will
 #'  be 1.
-#'@param memb_dim A character string indicating the name of member (nobs/nexp) 
-#'  dimension. The default value is 'member'.
+#'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
+#'  dimension. The default value is 'dataset'.
 #'@param time_dim A character string indicating the name of dimension along  
 #'  which the RMSSS are computed. The default value is 'sdate'.
 #'@param pval A logical value indicating whether to compute or not the p-value 
@@ -35,8 +35,8 @@
 #'@return
 #'A list containing the numeric arrays with dimension:\cr 
 #'  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-#'nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-#'number of observation (i.e., memb_dim in obs).\cr
+#'nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+#'number of observation (i.e., dat_dim in obs).\cr
 #'\item{$rmsss}{
 #'  The root mean square error skill score. 
 #'}
@@ -46,16 +46,16 @@
 #'
 #'@examples
 #' set.seed(1)
-#' exp <- array(rnorm(15), dim = c(dat = 1, time = 3, member = 5))
+#' exp <- array(rnorm(30), dim = c(dataset = 2, time = 3, memb = 5))
 #' set.seed(2)
-#' obs <- array(rnorm(6), dim = c(time = 3, member = 2, dat = 1))
+#' obs <- array(rnorm(15), dim = c(time = 3, memb = 5, dataset = 1))
 #' res <- RMSSS(exp, obs, time_dim = 'time')
 #'
 #'@rdname RMSSS
 #'@import multiApply
 #'@importFrom stats pf
 #'@export
-RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', 
+RMSSS <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset', 
                   pval = TRUE, ncores = NULL) {
   
   # Check inputs 
@@ -69,16 +69,16 @@ RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   if (is.null(dim(exp)) & is.null(dim(obs))) {  #is vector
     if (length(exp) == length(obs)) {
       exp <- array(exp, dim = c(length(exp), 1))
-      names(dim(exp)) <- c(time_dim, memb_dim)
+      names(dim(exp)) <- c(time_dim, dat_dim)
       obs <- array(obs, dim = c(length(obs), 1))
-      names(dim(obs)) <- c(time_dim, memb_dim)
+      names(dim(obs)) <- c(time_dim, dat_dim)
     } else {
     stop(paste0("Parameter 'exp' and 'obs' must be array with as least two ",
-                "dimensions time_dim and memb_dim, or vector of same length."))
+                "dimensions time_dim and dat_dim, or vector of same length."))
     }
   } else if (is.null(dim(exp)) | is.null(dim(obs))) {
     stop(paste0("Parameter 'exp' and 'obs' must be array with as least two ",
-                "dimensions time_dim and memb_dim, or vector of same length."))
+                "dimensions time_dim and dat_dim, or vector of same length."))
   }
   if(any(is.null(names(dim(exp))))| any(nchar(names(dim(exp))) == 0) |
      any(is.null(names(dim(obs))))| any(nchar(names(dim(obs))) == 0)) {
@@ -95,12 +95,12 @@ RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   if (!time_dim %in% names(dim(exp)) | !time_dim %in% names(dim(obs))) {
     stop("Parameter 'time_dim' is not found in 'exp' or 'obs' dimension.")
   }
-  ## memb_dim
-  if (!is.character(memb_dim) | length(memb_dim) > 1) {
-    stop("Parameter 'memb_dim' must be a character string.")
+  ## dat_dim
+  if (!is.character(dat_dim) | length(dat_dim) > 1) {
+    stop("Parameter 'dat_dim' must be a character string.")
   }
-  if (!memb_dim %in% names(dim(exp)) | !memb_dim %in% names(dim(obs))) {
-    stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
+  if (!dat_dim %in% names(dim(exp)) | !dat_dim %in% names(dim(obs))) {
+    stop("Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension.")
   }
   ## pval
   if (!is.logical(pval) | length(pval) > 1) {
@@ -116,11 +116,11 @@ RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   ## exp and obs (2)
   name_exp <- sort(names(dim(exp)))
   name_obs <- sort(names(dim(obs)))
-  name_exp <- name_exp[-which(name_exp == memb_dim)]
-  name_obs <- name_obs[-which(name_obs == memb_dim)]
+  name_exp <- name_exp[-which(name_exp == dat_dim)]
+  name_obs <- name_obs[-which(name_obs == dat_dim)]
   if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
     stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
-                "all dimension expect 'memb_dim'."))
+                "all dimension expect 'dat_dim'."))
   }
   if (dim(exp)[time_dim] <= 2) {
     stop("The length of time_dim must be more than 2 to compute RMSSS.")
@@ -139,47 +139,47 @@ RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   # Calculate RMSSS
   
   res <- Apply(list(exp, obs), 
-               target_dims = list(c(time_dim, memb_dim), 
-                                  c(time_dim, memb_dim)),
+               target_dims = list(c(time_dim, dat_dim), 
+                                  c(time_dim, dat_dim)),
                fun = .RMSSS, 
-               time_dim = time_dim, memb_dim = memb_dim,
+               time_dim = time_dim, dat_dim = dat_dim,
                pval = pval, #conf = conf, conf.lev = conf.lev,
                ncores = ncores)
   
   return(res)
 }
 
-.RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member', pval = TRUE) {
-  # exp: [sdate, member_exp]
-  # obs: [sdate, member_obs]
-  n_exp <- as.numeric(dim(exp)[2])
-  n_obs <- as.numeric(dim(obs)[2])
-  n_sdate <- as.numeric(dim(exp)[1])
+.RMSSS <- function(exp, obs, time_dim = 'sdate', dat_dim = 'dataset', pval = TRUE) {
+  # exp: [sdate, dat_exp]
+  # obs: [sdate, dat_obs]
+  nexp <- as.numeric(dim(exp)[2])
+  nobs <- as.numeric(dim(obs)[2])
+  nsdate <- as.numeric(dim(exp)[1])
   
-  p_val <- array(dim = c(nexp = n_exp, nobs = n_obs))
-  dif1 <- array(dim = c(n_sdate, n_exp, n_obs))
+  p_val <- array(dim = c(nexp = nexp, nobs = nobs))
+  dif1 <- array(dim = c(nsdate, nexp, nobs))
   names(dim(dif1)) <- c(time_dim, 'nexp', 'nobs')
 
 #  if (conf) {
 #    conflow <- (1 - conf.lev) / 2
 #    confhigh <- 1 - conflow
-#    conf_low <- array(dim = c(nexp = n_exp, nobs = n_obs))
-#    conf_high <- array(dim = c(nexp = n_exp, nobs = n_obs))
+#    conf_low <- array(dim = c(nexp = nexp, nobs = nobs))
+#    conf_high <- array(dim = c(nexp = nexp, nobs = nobs))
 #  }
 
   # dif1
-  for (i in 1:n_obs) {
-    dif1[, , i] <- sapply(1:n_exp, function(x) {exp[, x] - obs[, i]})
+  for (i in 1:nobs) {
+    dif1[, , i] <- sapply(1:nexp, function(x) {exp[, x] - obs[, i]})
   }
 
   # rms1 and eno1
-  rms1 <- apply(dif1^2, c(2, 3), mean, na.rm = TRUE)^0.5  #array(dim = c(n_exp, n_obs))
+  rms1 <- apply(dif1^2, c(2, 3), mean, na.rm = TRUE)^0.5  #array(dim = c(nexp, nobs))
   # rms2 and eno2
-  rms2 <- array(colMeans(obs^2, na.rm = TRUE)^0.5, dim = c(n_obs = n_obs))
+  rms2 <- array(colMeans(obs^2, na.rm = TRUE)^0.5, dim = c(nobs = nobs))
   rms2[which(abs(rms2) <= (max(abs(rms2), na.rm = TRUE) / 1000))] <- max(abs(
     rms2), na.rm = TRUE) / 1000
   #rms2 above: [nobs]
-  rms2 <- array(rms2, dim = c(nobs = n_obs, nexp = n_exp))
+  rms2 <- array(rms2, dim = c(nobs = nobs, nexp = nexp))
   #rms2 above: [nobs, nexp]
   rms2 <- Reorder(rms2, c(2, 1))
   #rms2 above: [nexp, nobs]
@@ -191,7 +191,7 @@ RMSSS <- function(exp, obs, time_dim = 'sdate', memb_dim = 'member',
   if (pval) {
     eno1 <- Eno(dif1, time_dim)  
     eno2 <- Eno(obs, time_dim)  
-    eno2 <- array(eno2, dim = c(nobs = n_obs, nexp = n_exp))
+    eno2 <- array(eno2, dim = c(nobs = nobs, nexp = nexp))
     eno2 <- Reorder(eno2, c(2, 1))
   }
 
diff --git a/man/Corr.Rd b/man/Corr.Rd
index 45eb166197dd3fef4ece125678d266915607cfb4..bf5575e443e13ff5a9a7d69eb166ec5045042ab3 100644
--- a/man/Corr.Rd
+++ b/man/Corr.Rd
@@ -4,22 +4,22 @@
 \alias{Corr}
 \title{Compute the correlation coefficient between an array of forecast and their corresponding observation}
 \usage{
-Corr(exp, obs, time_dim = "sdate", memb_dim = "member", comp_dim = NULL,
+Corr(exp, obs, time_dim = "sdate", dat_dim = "dataset", comp_dim = NULL,
   limits = NULL, method = "pearson", pval = TRUE, conf = TRUE,
   conf.lev = 0.95, ncores = NULL)
 }
 \arguments{
 \item{exp}{A named numeric array of experimental data, with at least two 
-dimensions 'time_dim' and 'memb_dim'.}
+dimensions 'time_dim' and 'dat_dim'.}
 
 \item{obs}{A named numeric array of observational data, same dimensions as
-parameter 'exp' except along memb_dim.}
+parameter 'exp' except along dat_dim.}
 
 \item{time_dim}{A character string indicating the name of dimension along
 which the correlations are computed. The default value is 'sdate'.}
 
-\item{memb_dim}{A character string indicating the name of member (nobs/nexp) 
-dimension. The default value is 'member'.}
+\item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
+dimension. The default value is 'dataset'.}
 
 \item{comp_dim}{A character string indicating the name of dimension along which
 obs is taken into account only if it is complete. The default value
@@ -46,8 +46,8 @@ computation. The default value is NULL.}
 \value{
 A list containing the numeric arrays with dimension:\cr 
  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-number of observation (i.e., memb_dim in obs).\cr
+nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+number of observation (i.e., dat_dim in obs).\cr
 \item{$corr}{
  The correlation coefficient. 
 }
@@ -71,13 +71,18 @@ limits[1] and limits[2]. This option can be activated if the user wants to
 account only for the forecasts which the corresponding observations are 
 available at all leadtimes.\cr 
 The confidence interval is computed by the Fisher transformation and the 
-significance level relies on an one-sided student-T distribution.\cr
+significance level relies on an one-sided student-T distribution.\cr 
+If the dataset has more than one member, ensemble mean is necessary necessary
+before using this function since it only allows one dimension 'dat_dim' to
+have inconsistent length between 'exp' and 'obs'. If all the dimensions of
+'exp' and 'obs' are identical, you can simply use apply() and cor() to 
+compute the correlation.
 }
 \examples{
 # Load sample data as in Load() example: 
 example(Load) 
 clim <- Clim(sampleData$mod, sampleData$obs) 
-corr <- Corr(clim$clim_exp, clim$clim_obs, time_dim = 'ftime')
+corr <- Corr(clim$clim_exp, clim$clim_obs, time_dim = 'ftime', dat_dim = 'member')
 # Renew the example when Ano and Smoothing is ready
 
 }
diff --git a/man/GMST.Rd b/man/GMST.Rd
index 3c9ace76f88fc47780cbf65b12c88ae5524134d2..be2d479eb788603207305eb38808e3b8ed6d2f37 100644
--- a/man/GMST.Rd
+++ b/man/GMST.Rd
@@ -107,7 +107,9 @@ mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
 sea_value <- 1
 lat <- seq(-90, 90, 10)
 lon <- seq(0, 360, 10)
-index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat, data_lons = lon, type = 'obs', mask_sea_land = mask_sea_land, sea_value = sea_value)
+index_obs <- GMST(data_tas = obs_tas, data_tos = obs_tos, data_lats = lat,
+                  data_lons = lon, type = 'obs', 
+                  mask_sea_land = mask_sea_land, sea_value = sea_value)
 
 ## Historical simulations
 hist_tas <- array(1:100, dim = c(year = 5, lat = 19, lon = 37, month = 12, member = 5))
@@ -116,7 +118,9 @@ mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
 sea_value <- 1
 lat <- seq(-90, 90, 10)
 lon <- seq(0, 360, 10)
-index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, sea_value = sea_value)
+index_hist <- GMST(data_tas = hist_tas, data_tos = hist_tos, data_lats = lat, 
+                   data_lons = lon, type = 'hist', mask_sea_land = mask_sea_land, 
+                   sea_value = sea_value)
 
 ## Decadal predictions
 dcpp_tas <- array(1:100, dim = c(sdate = 5, lat = 19, lon = 37, fmonth = 24, member = 5))
@@ -125,7 +129,9 @@ mask_sea_land <- array(c(1,0,1), dim = c(lat = 19, lon = 37))
 sea_value <- 1
 lat <- seq(-90, 90, 10)
 lon <- seq(0, 360, 10)
-index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land, sea_value = sea_value)
+index_dcpp <- GMST(data_tas = dcpp_tas, data_tos = dcpp_tos, data_lats = lat, 
+                   data_lons = lon, type = 'dcpp', monini = 1, mask_sea_land = mask_sea_land,
+                   sea_value = sea_value)
 
 }
 \author{
diff --git a/man/RMS.Rd b/man/RMS.Rd
index ac546867c422cdfac3e703279e437934267c1e9b..3ac1359ea069fd63b1765be06c91d31bebeb5f5e 100644
--- a/man/RMS.Rd
+++ b/man/RMS.Rd
@@ -4,21 +4,21 @@
 \alias{RMS}
 \title{Compute root mean square error}
 \usage{
-RMS(exp, obs, time_dim = "sdate", memb_dim = "member", comp_dim = NULL,
+RMS(exp, obs, time_dim = "sdate", dat_dim = "dataset", comp_dim = NULL,
   limits = NULL, conf = TRUE, conf.lev = 0.95, ncores = NULL)
 }
 \arguments{
 \item{exp}{A named numeric array of experimental data, with at least two 
-dimensions 'time_dim' and 'memb_dim'.}
+dimensions 'time_dim' and 'dat_dim'.}
 
 \item{obs}{A named numeric array of observational data, same dimensions as
-parameter 'exp' except along memb_dim.}
+parameter 'exp' except along dat_dim.}
 
 \item{time_dim}{A character string indicating the name of dimension along
 which the correlations are computed. The default value is 'sdate'.}
 
-\item{memb_dim}{A character string indicating the name of member (nobs/nexp) 
-dimension. The default value is 'member'.}
+\item{dat_dim}{A character string indicating the name of member (nobs/nexp) 
+dimension. The default value is 'dataset'.}
 
 \item{comp_dim}{A character string indicating the name of dimension along which
 obs is taken into account only if it is complete. The default value
@@ -39,8 +39,8 @@ computation. The default value is NULL.}
 \value{
 A list containing the numeric arrays with dimension:\cr 
  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-number of observation (i.e., memb_dim in obs).\cr
+nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+number of observation (i.e., dat_dim in obs).\cr
 \item{$rms}{
  The root mean square error. 
 }
@@ -65,9 +65,9 @@ The confidence interval is computed by the chi2 distribution.\cr
 \examples{
 # Load sample data as in Load() example:
  set.seed(1)
- exp1 <- array(rnorm(120), dim = c(member = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
+ exp1 <- array(rnorm(120), dim = c(dataset = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
  set.seed(2)
- obs1 <- array(rnorm(80),  dim = c(member = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
+ obs1 <- array(rnorm(80),  dim = c(dataset = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
  set.seed(2)
  na <- floor(runif(10, min = 1, max = 80))
  obs1[na] <- NA
diff --git a/man/RMSSS.Rd b/man/RMSSS.Rd
index 0939a94730c08ec0e297064dcdccdd9f61d85297..1b8274f41c159812f639b8e580190f2ad3e0dfac 100644
--- a/man/RMSSS.Rd
+++ b/man/RMSSS.Rd
@@ -4,27 +4,27 @@
 \alias{RMSSS}
 \title{Compute root mean square error skill score}
 \usage{
-RMSSS(exp, obs, time_dim = "sdate", memb_dim = "member", pval = TRUE,
+RMSSS(exp, obs, time_dim = "sdate", dat_dim = "dataset", pval = TRUE,
   ncores = NULL)
 }
 \arguments{
 \item{exp}{A named numeric array of experimental data which contains at least
-two dimensions for memb_dim and time_dim. It can also be a vector with the 
+two dimensions for dat_dim and time_dim. It can also be a vector with the 
 same length as 'obs', then the vector will automatically be 'time_dim' and 
-'memb_dim' will be 1.}
+'dat_dim' will be 1.}
 
 \item{obs}{A named numeric array of observational data which contains at least
-two dimensions for memb_dim and time_dim. The dimensions should be the same 
-as paramter 'exp' except the length of 'memb_dim' dimension. The order of 
+two dimensions for dat_dim and time_dim. The dimensions should be the same 
+as paramter 'exp' except the length of 'dat_dim' dimension. The order of 
 dimension can be different. It can also be a vector with the same length as
-'exp', then the vector will automatically be 'time_dim' and 'memb_dim' will
+'exp', then the vector will automatically be 'time_dim' and 'dat_dim' will
 be 1.}
 
 \item{time_dim}{A character string indicating the name of dimension along
 which the RMSSS are computed. The default value is 'sdate'.}
 
-\item{memb_dim}{A character string indicating the name of member (nobs/nexp) 
-dimension. The default value is 'member'.}
+\item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
+dimension. The default value is 'dataset'.}
 
 \item{pval}{A logical value indicating whether to compute or not the p-value 
 of the test Ho: RMSSS = 0. If pval = TRUE, the insignificant RMSSS will 
@@ -36,8 +36,8 @@ computation. The default value is NULL.}
 \value{
 A list containing the numeric arrays with dimension:\cr 
  c(nexp, nobs, all other dimensions of exp except time_dim).\cr
-nexp is the number of experiment (i.e., memb_dim in exp), and nobs is the 
-number of observation (i.e., memb_dim in obs).\cr
+nexp is the number of experiment (i.e., dat_dim in exp), and nobs is the 
+number of observation (i.e., dat_dim in obs).\cr
 \item{$rmsss}{
  The root mean square error skill score. 
 }
@@ -48,7 +48,7 @@ number of observation (i.e., memb_dim in obs).\cr
 \description{
 Compute the root mean square error skill score (RMSSS) between an array of 
 forecast 'exp' and an array of observation 'obs'. The two arrays should 
-have the same dimensions except along memb_dim, where the length can be 
+have the same dimensions except along dat_dim, where the length can be 
 different, with the number of experiments/models (nexp) and the number of 
 observational datasets (nobs).\cr
 RMSSS computes the root mean square error skill score of each jexp in 1:nexp 
@@ -60,9 +60,9 @@ The p-value is optionally provided by an one-sided Fisher test.\cr
 }
 \examples{
 set.seed(1)
-exp <- array(rnorm(15), dim = c(dat = 1, time = 3, member = 5))
+exp <- array(rnorm(30), dim = c(dataset = 2, time = 3, memb = 5))
 set.seed(2)
-obs <- array(rnorm(6), dim = c(time = 3, member = 2, dat = 1))
+obs <- array(rnorm(15), dim = c(time = 3, memb = 5, dataset = 1))
 res <- RMSSS(exp, obs, time_dim = 'time')
 
 }
diff --git a/tests/testthat/test-Corr.R b/tests/testthat/test-Corr.R
index ca694e04fed15e88c5745e97718c2e30b54c550c..4a06f8247fc378cfba5b28abc51f536be8ccfdf0 100644
--- a/tests/testthat/test-Corr.R
+++ b/tests/testthat/test-Corr.R
@@ -3,11 +3,11 @@ context("s2dv::Corr tests")
 ##############################################
   # dat1
   set.seed(1)
-  exp1 <- array(rnorm(240), dim = c(dataset = 1, member = 2, sdate = 5,
+  exp1 <- array(rnorm(240), dim = c(member = 1, dataset = 2, sdate = 5,
                                    ftime = 3, lat = 2, lon = 4))
 
   set.seed(2)
-  obs1 <- array(rnorm(120), dim = c(dataset = 1, member = 1, sdate = 5,
+  obs1 <- array(rnorm(120), dim = c(member = 1, dataset = 1, sdate = 5,
                                    ftime = 3, lat = 2, lon = 4))
   set.seed(2)
   na <- floor(runif(10, min = 1, max = 120))
@@ -27,7 +27,7 @@ test_that("1. Input checks", {
   expect_error(
   Corr(c(1:10), c(2:4)),
   paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
-         "containing time_dim and memb_dim.")
+         "containing time_dim and dat_dim.")
   )
   expect_error(
   Corr(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
@@ -38,12 +38,12 @@ test_that("1. Input checks", {
   "Parameter 'exp' and 'obs' must have same dimension name"
   )
   expect_error(
-  Corr(exp1, obs1, memb_dim = 1),
-  "Parameter 'memb_dim' must be a character string."
+  Corr(exp1, obs1, dat_dim = 1),
+  "Parameter 'dat_dim' must be a character string."
   )
   expect_error(
-  Corr(exp1, obs1, memb_dim = 'a'),
-  "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension."
+  Corr(exp1, obs1, dat_dim = 'a'),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
   )
   expect_error(
   Corr(exp1, obs1, time_dim = c('sdate', 'a')),
@@ -91,13 +91,13 @@ test_that("1. Input checks", {
   "Parameter 'ncores' must be a positive integer."
   )
   expect_error(
-  Corr(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)),
-       obs = array(1:4, dim = c(a = 1, sdate = 2, member = 2))),
-  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'memb_dim'."
+  Corr(exp = array(1:10, dim = c(sdate = 1, dataset = 5, a = 1)),
+       obs = array(1:4, dim = c(a = 1, sdate = 2, dataset = 2))),
+  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'dat_dim'."
   )
   expect_error(
-  Corr(exp = array(1:10, dim = c(sdate = 2, member = 5, a = 1)),
-       obs = array(1:4, dim = c(a = 1, sdate = 2, member = 2))),
+  Corr(exp = array(1:10, dim = c(sdate = 2, dataset = 5, a = 1)),
+       obs = array(1:4, dim = c(a = 1, sdate = 2, dataset = 2))),
   "The length of time_dim must be at least 3 to compute correlation."
   )
 
@@ -108,7 +108,7 @@ test_that("2. Output checks: dat1", {
 
   expect_equal(
     dim(Corr(exp1, obs1)$corr),
-    c(n_exp = 2, n_obs = 1, dataset = 1, ftime = 3, lat = 2, lon = 4)
+    c(nexp = 2, nobs = 1, member = 1, ftime = 3, lat = 2, lon = 4)
   )
   expect_equal(
     Corr(exp1, obs1)$corr[1:6],
diff --git a/tests/testthat/test-RMS.R b/tests/testthat/test-RMS.R
index b5a1b994bdfe55f9e4f81d366ec5295d9ef66a90..7b3e7b973cf09cf8ad5840945a434af75ffb71c7 100644
--- a/tests/testthat/test-RMS.R
+++ b/tests/testthat/test-RMS.R
@@ -3,10 +3,10 @@ context("s2dv::RMS tests")
 ##############################################
   # dat1
   set.seed(1)
-  exp1 <- array(rnorm(120), dim = c(member = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
+  exp1 <- array(rnorm(120), dim = c(dataset = 3, sdate = 5, ftime = 2, lon = 1, lat = 4))
 
   set.seed(2)
-  obs1 <- array(rnorm(80),  dim = c(member = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
+  obs1 <- array(rnorm(80),  dim = c(dataset = 2, sdate = 5, ftime = 2, lon = 1, lat = 4))
   set.seed(2)
   na <- floor(runif(10, min = 1, max = 80))
   obs1[na] <- NA
@@ -25,7 +25,7 @@ test_that("1. Input checks", {
   expect_error(
   RMS(c(1:10), c(2:4)),
   paste0("Parameter 'exp' and 'obs' must be at least two dimensions ",
-         "containing time_dim and memb_dim.")
+         "containing time_dim and dat_dim.")
   )
   expect_error(
   RMS(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
@@ -36,12 +36,12 @@ test_that("1. Input checks", {
   "Parameter 'exp' and 'obs' must have same dimension name"
   )
   expect_error(
-  RMS(exp1, obs1, memb_dim = 1),
-  "Parameter 'memb_dim' must be a character string."
+  RMS(exp1, obs1, dat_dim = 1),
+  "Parameter 'dat_dim' must be a character string."
   )
   expect_error(
-  RMS(exp1, obs1, memb_dim = 'a'),
-  "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension."
+  RMS(exp1, obs1, dat_dim = 'a'),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
   )
   expect_error(
   RMS(exp1, obs1, time_dim = c('sdate', 'a')),
@@ -81,13 +81,13 @@ test_that("1. Input checks", {
   "Parameter 'ncores' must be a positive integer."
   )
   expect_error(
-  RMS(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)),
-       obs = array(1:4, dim = c(a = 1, sdate = 2, member = 2))),
-  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'memb_dim'."
+  RMS(exp = array(1:10, dim = c(sdate = 1, dataset = 5, a = 1)),
+       obs = array(1:4, dim = c(a = 1, sdate = 2, dataset = 2))),
+  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'dat_dim'."
   )
   expect_error(
-  RMS(exp = array(1:5, dim = c(sdate = 1, member = 5, a = 1)),
-       obs = array(1:2, dim = c(a = 1, sdate = 1, member = 2))),
+  RMS(exp = array(1:5, dim = c(sdate = 1, dataset = 5, a = 1)),
+       obs = array(1:2, dim = c(a = 1, sdate = 1, dataset = 2))),
   "The length of time_dim must be at least 2 to compute RMS."
   )
 
@@ -100,7 +100,7 @@ test_that("2. Output checks: dat1", {
 
   expect_equal(
     dim(RMS(exp1, obs1)$rms),
-    c(n_exp = 3, n_obs = 2, ftime = 2, lon = 1, lat = 4)
+    c(nexp = 3, nobs = 2, ftime = 2, lon = 1, lat = 4)
   )
   expect_equal(
     RMS(exp1, obs1)$rms[1:6],
diff --git a/tests/testthat/test-RMSSS.R b/tests/testthat/test-RMSSS.R
index 7d5cc8d4223b2b0f1007e16d62a35ce561e91602..7d2a16da9da52c492dad263985f71c042783927b 100644
--- a/tests/testthat/test-RMSSS.R
+++ b/tests/testthat/test-RMSSS.R
@@ -3,9 +3,9 @@ context("s2dv::RMSSS tests")
 ##############################################
  # case 0
  set.seed(1)
- exp0 <- array(rnorm(15), dim = c(sdate = 3, member = 5))
+ exp0 <- array(rnorm(15), dim = c(sdate = 3, dataset = 5))
  set.seed(2)
- obs0 <- array(rnorm(6), dim = c(sdate = 3, member = 2))
+ obs0 <- array(rnorm(6), dim = c(sdate = 3, dataset = 2))
 
  # case 1
  set.seed(1)
@@ -15,9 +15,9 @@ context("s2dv::RMSSS tests")
  
  # case 2
  set.seed(3)
- exp2 <- array(rnorm(120), dim = c(sdate = 10, dat = 1, lon = 3, lat = 2, member = 2))
+ exp2 <- array(rnorm(120), dim = c(sdate = 10, dat = 1, lon = 3, lat = 2, dataset = 2))
  set.seed(4)
- obs2 <- array(rnorm(60), dim = c(dat = 1, sdate = 10, member = 1, lat = 2, lon = 3))
+ obs2 <- array(rnorm(60), dim = c(dat = 1, sdate = 10, dataset = 1, lat = 2, lon = 3))
 
  # case 3: vector
  set.seed(5)
@@ -40,7 +40,7 @@ test_that("1. Input checks", {
   expect_error(
   RMSSS(c(1:10), c(2:4)),
   paste0("Parameter 'exp' and 'obs' must be array with as least two dimensions ",
-         "time_dim and memb_dim, or vector of same length.")
+         "time_dim and dat_dim, or vector of same length.")
   )
   expect_error(
   RMSSS(array(1:10, dim = c(2, 5)), array(1:10, dim = c(2, 5))),
@@ -59,12 +59,12 @@ test_that("1. Input checks", {
   "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
   )
   expect_error(
-  RMSSS(exp0, obs0, memb_dim = NA),
-  "Parameter 'memb_dim' must be a character string."
+  RMSSS(exp0, obs0, dat_dim = NA),
+  "Parameter 'dat_dim' must be a character string."
   )
   expect_error(
-  RMSSS(exp0, obs0, memb_dim = 'memb'),
-  "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension."
+  RMSSS(exp0, obs0, dat_dim = 'memb'),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
   )
   expect_error(
   RMSSS(exp0, obs0, pval = c(T, T)),
@@ -75,13 +75,13 @@ test_that("1. Input checks", {
   "Parameter 'ncores' must be a positive integer."
   )
   expect_error(
-  RMSSS(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)),
-        obs = array(1:4, dim = c(a = 1, sdate = 2, member = 2))),
-  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'memb_dim'."
+  RMSSS(exp = array(1:10, dim = c(sdate = 1, dataset = 5, a = 1)),
+        obs = array(1:4, dim = c(a = 1, sdate = 2, dataset = 2))),
+  "Parameter 'exp' and 'obs' must have same length of all dimension expect 'dat_dim'."
   )
   expect_error(
-  RMSSS(exp = array(1:10, dim = c(sdate = 1, member = 5, a = 1)),
-        obs = array(1:4, dim = c(a = 1, sdate = 1, member = 2))),
+  RMSSS(exp = array(1:10, dim = c(sdate = 1, dataset = 5, a = 1)),
+        obs = array(1:4, dim = c(a = 1, sdate = 1, dataset = 2))),
   "The length of time_dim must be more than 2 to compute RMSSS."
   )
 })
@@ -89,7 +89,7 @@ test_that("1. Input checks", {
 ##############################################
 test_that("2. Output checks: case 1", {
 
-  res1_1 <- RMSSS(exp1, obs1, time_dim = 'time', memb_dim = 'memb')
+  res1_1 <- RMSSS(exp1, obs1, time_dim = 'time', dat_dim = 'memb')
   expect_equal(
     dim(res1_1$rmsss),
     c(nexp = 5, nobs = 2)
@@ -108,7 +108,7 @@ test_that("2. Output checks: case 1", {
   exp1_2[2:4] <- NA
   obs1_2 <- obs1
   obs1_2[1:2] <- NA
-  res1_2 <- RMSSS(exp1_2, obs1_2, time_dim = 'time', memb_dim = 'memb', pval = TRUE)
+  res1_2 <- RMSSS(exp1_2, obs1_2, time_dim = 'time', dat_dim = 'memb', pval = TRUE)
 
   expect_equal(
     length(res1_2$rmsss[which(is.na(res1_2$rmsss))]),