diff --git a/NAMESPACE b/NAMESPACE
index a5229e3d20dd73f6cd5ca216897e937ad3507c79..d440ec097f07df8d3c0b7e37629587de38d0a9d9 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -2,9 +2,11 @@
 
 export(ACC)
 export(AMV)
+export(AbsBiasSS)
 export(AnimateMap)
 export(Ano)
 export(Ano_CrossValid)
+export(Bias)
 export(BrierScore)
 export(CDORemap)
 export(CRPS)
diff --git a/R/AbsBiasSS.R b/R/AbsBiasSS.R
new file mode 100644
index 0000000000000000000000000000000000000000..0838f251a97e62a4e82f65423a6ccd567b2b6b77
--- /dev/null
+++ b/R/AbsBiasSS.R
@@ -0,0 +1,280 @@
+#'Compute the Absolute Mean Bias Skill Score
+#'
+#'The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks,
+#' 2011) between the ensemble mean forecast and the observations. It measures 
+#'the accuracy of the forecast in comparison with a reference forecast to assess
+#'whether the forecast presents an improvement or a worsening with respect to 
+#'that reference. The Mean Bias Skill Score ranges between minus infinite and 1.
+#'Positive values indicate that the forecast has higher skill than the reference
+#'forecast, while negative values indicate that it has a lower skill. Examples
+#'of reference forecasts are the climatological forecast (average of the 
+#'observations), a previous model version, or another model. It is computed as
+#'\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance
+#'is obtained based on a Random Walk test at the 95% confidence level (DelSole 
+#'and Tippett, 2016). If there is more than one dataset, the result will be 
+#'computed for each pair of exp and obs data.
+#'
+#'@param exp A named numerical array of the forecast with at least time 
+#'  dimension.  
+#'@param obs A named numerical array of the observation with at least time 
+#'  dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 
+#'  'dat_dim'.
+#'@param ref A named numerical array of the reference forecast data with at 
+#'  least time dimension. The dimensions must be the same as 'exp' except
+#'  'memb_dim' and 'dat_dim'. If there is only one reference dataset, it should
+#'  not have dataset dimension. If there is corresponding reference for each
+#'  experiement, the dataset dimension must have the same length as in 'exp'. If
+#'  'ref' is NULL, the climatological forecast is used as reference forecast.
+#'  The default value is NULL.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'sdate'.
+#'@param memb_dim A character string indicating the name of the member dimension
+#'  to compute the ensemble mean; it should be set to NULL if the parameter 'exp'
+#'  and 'ref' are already the ensemble mean. The default value is NULL.
+#'@param dat_dim A character string indicating the name of dataset dimension. 
+#'  The length of this dimension can be different between 'exp' and 'obs'. 
+#'  The default value is NULL.
+#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or
+#'  kept (FALSE) for computation. The default value is FALSE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return
+#'\item{$biasSS}{
+#'  A numerical array of BiasSS with dimensions nexp, nobs and the rest 
+#'  dimensions of 'exp' except 'time_dim' and 'memb_dim'.
+#'}
+#'\item{$sign}{
+#'  A logical array of the statistical significance of the BiasSS
+#'  with the same dimensions as $biasSS. 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). If dat_dim is NULL, nexp and nobs are omitted.
+#'}
+#'
+#'@references 
+#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+#'DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1
+#'
+#'@examples
+#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+#'ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
+#'biasSS1 <- AbsBiasSS(exp = exp, obs = obs, ref = ref, memb_dim = 'member')
+#'biasSS2 <- AbsBiasSS(exp = exp, obs = obs, ref = NULL, memb_dim = 'member')
+#'
+#'@import multiApply
+#'@export
+AbsBiasSS <- function(exp, obs, ref = NULL, time_dim = 'sdate', memb_dim = NULL, 
+                      dat_dim = NULL, na.rm = FALSE, ncores = NULL) {
+  
+  # Check inputs
+  ## exp, obs, and ref (1)
+  if (!is.array(exp) | !is.numeric(exp)) {
+    stop("Parameter 'exp' must be a numeric array.")
+  }
+  if (!is.array(obs) | !is.numeric(obs)) {
+    stop("Parameter 'obs' must be a numeric array.")
+  }
+  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)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  if (!is.null(ref)) {
+    if (!is.array(ref) | !is.numeric(ref))
+      stop("Parameter 'ref' must be a numeric array.")
+    if (any(is.null(names(dim(ref))))| any(nchar(names(dim(ref))) == 0)) {
+      stop("Parameter 'ref' must have dimension names.")
+    }
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1) {
+    stop("Parameter 'time_dim' must be a character string.")
+  }
+  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.")
+  }
+  if (!is.null(ref) & !time_dim %in% names(dim(ref))) {
+    stop("Parameter 'time_dim' is not found in 'ref' dimension.")
+  }
+  ## memb_dim
+  if (!is.null(memb_dim)) {
+    if (!is.character(memb_dim) | length(memb_dim) > 1) {
+      stop("Parameter 'memb_dim' must be a character string.")
+    }
+    if (!memb_dim %in% names(dim(exp))) {
+      stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+    }
+    if (memb_dim %in% names(dim(obs))) {
+      if (identical(as.numeric(dim(obs)[memb_dim]), 1)) {
+        obs <- ClimProjDiags::Subset(x = obs, along = memb_dim, indices = 1, drop = 'selected')
+      } else {
+        stop("Not implemented for observations with members ('obs' can have 'memb_dim', ",
+             "but it should be of length = 1).")
+      }
+    }
+  }
+  ## dat_dim
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim) | length(dat_dim) > 1) {
+      stop("Parameter 'dat_dim' must be a character string.")
+    }
+    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.",
+           " Set it as NULL if there is no dataset dimension.")
+    }
+  }
+  ## exp, obs, and ref (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  if (!is.null(memb_dim)) {
+    name_exp <- name_exp[-which(name_exp == memb_dim)]
+  }
+  if (!is.null(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim)]
+    name_obs <- name_obs[-which(name_obs == dat_dim)]
+  }
+  if (!identical(length(name_exp), length(name_obs)) |
+      !identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimensions except 'memb_dim' and 'dat_dim'."))
+  }
+  if (!is.null(ref)) {
+    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])) {
+          stop(paste0("If parameter 'ref' has dataset dimension, it must be", 
+                      " equal to dataset dimension of 'exp'."))
+        }
+        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])) {
+      stop(paste0("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."))
+    }
+  }
+  ## na.rm
+  if (!is.logical(na.rm) | length(na.rm) > 1) {
+    stop("Parameter 'na.rm' must be one logical value.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be either NULL or a positive integer.")
+    }
+  }  
+ 
+ ############################
+ 
+  ## Ensemble mean
+  if (!is.null(memb_dim)) {
+    exp <- MeanDims(exp, memb_dim, na.rm = na.rm)
+    if (!is.null(ref) & memb_dim %in% names(dim(ref))) {
+      ref <- MeanDims(ref, memb_dim, na.rm = na.rm)
+    }
+  }
+  
+  ## Mean bias skill score
+  if (!is.null(ref)) { # use "ref" as reference forecast
+    if (!is.null(dat_dim) && dat_dim %in% names(dim(ref))) {
+      target_dims_ref <- c(time_dim, dat_dim)
+    } else {
+      target_dims_ref <- c(time_dim)
+    }
+    data <- list(exp = exp, obs = obs, ref = ref)
+    target_dims = list(exp = c(time_dim, dat_dim),
+                       obs = c(time_dim, dat_dim),
+                       ref = target_dims_ref)
+  } else {
+    data <- list(exp = exp, obs = obs)
+    target_dims = list(exp = c(time_dim, dat_dim),
+                       obs = c(time_dim, dat_dim))
+  }
+
+  output <- Apply(data,
+                  target_dims = target_dims,
+                  fun = .AbsBiasSS,
+                  dat_dim = dat_dim, 
+                  na.rm = na.rm, 
+                  ncores = ncores)
+
+  return(output)
+}
+
+.AbsBiasSS <- function(exp, obs, ref = NULL, dat_dim = NULL, na.rm = FALSE) {
+  # exp and obs: [sdate, (dat_dim)]
+  # ref: [sdate, (dat_dim)] or NULL
+
+  # Adjust exp, obs, ref to have dat_dim temporarily
+  if (is.null(dat_dim)) {
+    nexp <- 1
+    nobs <- 1
+    exp <- InsertDim(exp, posdim = 2, lendim = 1, name = 'dataset')
+    obs <- InsertDim(obs, posdim = 2, lendim = 1, name = 'dataset')
+    if (!is.null(ref)) {
+      ref <- InsertDim(ref, posdim = 2, lendim = 1, name = 'dataset')
+    }
+    ref_dat_dim <- FALSE
+  } else {
+    nexp <- as.numeric(dim(exp)[dat_dim])
+    nobs <- as.numeric(dim(obs)[dat_dim])
+    if (length(dim(ref)) == 1) { # ref: [sdate]
+      ref_dat_dim <- FALSE
+    } else {
+      ref_dat_dim <- TRUE
+    }
+  }
+
+  biasSS <- array(dim = c(nexp = nexp, nobs = nobs))
+  sign <- array(dim = c(nexp = nexp, nobs = nobs))
+
+  for (i in 1:nexp) {
+    exp_data <- exp[, i]
+    if (isTRUE(ref_dat_dim)) {
+      ref_data <- ref[, i]
+    } else {
+      ref_data <- ref
+    }
+    for (j in 1:nobs) {
+      obs_data <- obs[, j]
+
+      if (isTRUE(na.rm)) {
+        if (is.null(ref)) {
+          good_values <- !is.na(exp_data) & !is.na(obs_data)
+          exp_data <- exp_data[good_values]
+          obs_data <- obs_data[good_values]
+        } else {
+          good_values <- !is.na(exp_data) & !is.na(ref_data) & !is.na(obs_data)
+          exp_data <- exp_data[good_values]
+          ref_data <- ref_data[good_values]
+          obs_data <- obs_data[good_values]
+        }
+      }
+      ## Bias of the exp
+      bias_exp <- .Bias(exp = exp_data, obs = obs_data, na.rm = na.rm, absolute = TRUE, time_mean = TRUE)
+      ## Bias of the ref
+      if (is.null(ref)) { ## Climatological forecast
+        ref_data <- mean(obs_data, na.rm = na.rm)
+      }
+      bias_ref <- .Bias(exp = ref_data, obs = obs_data, na.rm = na.rm, absolute = TRUE, time_mean = TRUE)
+      ## Skill score and significance
+      biasSS[i, j] <- 1 - bias_exp / bias_ref
+      sign[i, j] <- .RandomWalkTest(skill_A = bias_exp, skill_B = bias_ref)$signif
+    }
+  }
+
+  if (is.null(dat_dim)) {
+    dim(biasSS) <- NULL
+    dim(sign) <- NULL
+  }
+  
+  
+  return(list(biasSS = biasSS, sign = sign))
+}
diff --git a/R/Bias.R b/R/Bias.R
new file mode 100644
index 0000000000000000000000000000000000000000..0319a0f08e23e388046ce895e7a06aa82a0d9a41
--- /dev/null
+++ b/R/Bias.R
@@ -0,0 +1,189 @@
+#'Compute the Mean Bias
+#'
+#'The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference 
+#'between the ensemble mean forecast and the observations. It is a deterministic
+#'metric. Positive values indicate that the forecasts are on average too high
+#'and negative values indicate that the forecasts are on average too low.
+#'It also allows to compute the Absolute Mean Bias or bias without temporal 
+#'mean. If there is more than one dataset, the result will be computed for each
+#'pair of exp and obs data.
+#'
+#'@param exp A named numerical array of the forecast with at least time 
+#'  dimension.  
+#'@param obs A named numerical array of the observation with at least time 
+#'  dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 
+#'  'dat_dim'.
+#'@param time_dim A character string indicating the name of the time dimension.
+#'  The default value is 'sdate'.
+#'@param dat_dim A character string indicating the name of dataset dimension. 
+#'  The length of this dimension can be different between 'exp' and 'obs'. 
+#'  The default value is NULL.
+#'@param memb_dim A character string indicating the name of the member dimension
+#'  to compute the ensemble mean; it should be set to NULL if the parameter 
+#'  'exp' is already the ensemble mean. The default value is NULL.
+#'@param na.rm A logical value indicating if NAs should be removed (TRUE) or
+#'  kept (FALSE) for computation. The default value is FALSE.
+#'@param absolute A logical value indicating whether to compute the absolute 
+#'  bias. The default value is FALSE.
+#'@param time_mean A logical value indicating whether to compute the temporal 
+#'  mean of the bias. The default value is TRUE.
+#'@param ncores An integer indicating the number of cores to use for parallel 
+#'  computation. The default value is NULL.
+#'
+#'@return
+#'A numerical array of bias with dimensions c(nexp, nobs, the rest dimensions of
+#''exp' except 'time_dim' (if time_mean = T) and 'memb_dim'). 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). If dat_dim is NULL, nexp and nobs are omitted.
+#'
+#'@references 
+#'Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+#'
+#'@examples
+#'exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+#'obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
+#'bias <- Bias(exp = exp, obs = obs, memb_dim = 'member')
+#'
+#'@import multiApply
+#'@importFrom ClimProjDiags Subset
+#'@export
+Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, dat_dim = NULL, na.rm = FALSE, 
+                 absolute = FALSE, time_mean = TRUE, ncores = NULL) {
+  
+  # Check inputs
+  ## exp and obs (1)
+  if (!is.array(exp) | !is.numeric(exp))
+    stop("Parameter 'exp' must be a numeric array.")
+  if (!is.array(obs) | !is.numeric(obs))
+    stop("Parameter 'obs' must be a numeric array.")
+  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)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  ## time_dim
+  if (!is.character(time_dim) | length(time_dim) != 1)
+    stop("Parameter 'time_dim' must be a character string.")
+  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.null(memb_dim)) {
+    if (!is.character(memb_dim) | length(memb_dim) > 1) {
+      stop("Parameter 'memb_dim' must be a character string.")
+    }
+    if (!memb_dim %in% names(dim(exp))) {
+      stop("Parameter 'memb_dim' is not found in 'exp' dimension.")
+    }
+    if (memb_dim %in% names(dim(obs))) {
+      if (identical(as.numeric(dim(obs)[memb_dim]), 1)) {
+        obs <- ClimProjDiags::Subset(x = obs, along = memb_dim, indices = 1, drop = 'selected')
+      } else {
+        stop("Not implemented for observations with members ('obs' can have 'memb_dim', ",
+             "but it should be of length = 1).")
+      }
+    }
+  }
+  ## dat_dim
+  if (!is.null(dat_dim)) {
+    if (!is.character(dat_dim) | length(dat_dim) > 1) {
+      stop("Parameter 'dat_dim' must be a character string.")
+    }
+    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.",
+           " Set it as NULL if there is no dataset dimension.")
+    }
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  if (!is.null(memb_dim)) {
+    name_exp <- name_exp[-which(name_exp == memb_dim)]
+  }
+  if (!is.null(dat_dim)) {
+    name_exp <- name_exp[-which(name_exp == dat_dim)]
+    name_obs <- name_obs[-which(name_obs == dat_dim)]
+  }
+  if (!identical(length(name_exp), length(name_obs)) |
+      !identical(dim(exp)[name_exp], dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all dimensions except 'memb_dim' and 'dat_dim'."))
+  }
+  ## na.rm
+  if (!is.logical(na.rm) | length(na.rm) > 1) {
+    stop("Parameter 'na.rm' must be one logical value.")
+  }
+  ## absolute
+  if (!is.logical(absolute) | length(absolute) > 1) {
+    stop("Parameter 'absolute' must be one logical value.")
+  }
+  ## time_mean
+  if (!is.logical(time_mean) | length(time_mean) > 1) {
+    stop("Parameter 'time_mean' must be one logical value.")
+  }
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | ncores %% 1 != 0 | ncores <= 0 |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be either NULL or a positive integer.")
+    }
+  }
+
+  ###############################
+
+  ## Ensemble mean
+  if (!is.null(memb_dim)) {
+    exp <- MeanDims(exp, memb_dim, na.rm = na.rm)
+  }
+  
+  ## (Mean) Bias
+  bias <- Apply(data = list(exp, obs),
+                target_dims = c(time_dim, dat_dim),
+                fun = .Bias, 
+                time_dim = time_dim,
+                dat_dim = dat_dim, 
+                na.rm = na.rm,
+                absolute = absolute,
+                time_mean = time_mean,
+                ncores = ncores)$output1
+  
+  return(bias)
+}
+
+
+.Bias <- function(exp, obs, time_dim = 'sdate', dat_dim = NULL, na.rm = FALSE, 
+                  absolute = FALSE, time_mean = TRUE) {
+  # exp and obs: [sdate, (dat)]
+
+  if (is.null(dat_dim)) {
+    bias <- exp - obs
+    
+    if (isTRUE(absolute)) {
+      bias <- abs(bias)
+    } 
+    
+    if (isTRUE(time_mean)) {
+      bias <- mean(bias, na.rm = na.rm)
+    }
+
+  } else {
+    nexp <- as.numeric(dim(exp)[dat_dim])
+    nobs <- as.numeric(dim(obs)[dat_dim])
+    bias <- array(dim = c(dim(exp)[time_dim], nexp = nexp, nobs = nobs))
+
+    for (i in 1:nexp) {
+      for (j in 1:nobs) {
+        bias[, i, j] <-  exp[, i] - obs[, j]
+      }
+    }
+
+    if (isTRUE(absolute)) {
+      bias <- abs(bias)
+    }
+
+    if (isTRUE(time_mean)) {
+      bias <- MeanDims(bias, time_dim, na.rm = na.rm)
+    }
+  }
+
+  return(bias)
+}
diff --git a/man/AbsBiasSS.Rd b/man/AbsBiasSS.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..029101d06d6c3f31b69d8936ef3852b2702e2c8c
--- /dev/null
+++ b/man/AbsBiasSS.Rd
@@ -0,0 +1,89 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/AbsBiasSS.R
+\name{AbsBiasSS}
+\alias{AbsBiasSS}
+\title{Compute the Absolute Mean Bias Skill Score}
+\usage{
+AbsBiasSS(
+  exp,
+  obs,
+  ref = NULL,
+  time_dim = "sdate",
+  memb_dim = NULL,
+  dat_dim = NULL,
+  na.rm = FALSE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numerical array of the forecast with at least time 
+dimension.}
+
+\item{obs}{A named numerical array of the observation with at least time 
+dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 
+'dat_dim'.}
+
+\item{ref}{A named numerical array of the reference forecast data with at 
+least time dimension. The dimensions must be the same as 'exp' except
+'memb_dim' and 'dat_dim'. If there is only one reference dataset, it should
+not have dataset dimension. If there is corresponding reference for each
+experiement, the dataset dimension must have the same length as in 'exp'. If
+'ref' is NULL, the climatological forecast is used as reference forecast.
+The default value is NULL.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'sdate'.}
+
+\item{memb_dim}{A character string indicating the name of the member dimension
+to compute the ensemble mean; it should be set to NULL if the parameter 'exp'
+and 'ref' are already the ensemble mean. The default value is NULL.}
+
+\item{dat_dim}{A character string indicating the name of dataset dimension. 
+The length of this dimension can be different between 'exp' and 'obs'. 
+The default value is NULL.}
+
+\item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or
+kept (FALSE) for computation. The default value is FALSE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+\item{$biasSS}{
+ A numerical array of BiasSS with dimensions nexp, nobs and the rest 
+ dimensions of 'exp' except 'time_dim' and 'memb_dim'.
+}
+\item{$sign}{
+ A logical array of the statistical significance of the BiasSS
+ with the same dimensions as $biasSS. 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). If dat_dim is NULL, nexp and nobs are omitted.
+}
+}
+\description{
+The Absolute Mean Bias Skill Score is based on the Absolute Mean Error (Wilks,
+2011) between the ensemble mean forecast and the observations. It measures 
+the accuracy of the forecast in comparison with a reference forecast to assess
+whether the forecast presents an improvement or a worsening with respect to 
+that reference. The Mean Bias Skill Score ranges between minus infinite and 1.
+Positive values indicate that the forecast has higher skill than the reference
+forecast, while negative values indicate that it has a lower skill. Examples
+of reference forecasts are the climatological forecast (average of the 
+observations), a previous model version, or another model. It is computed as
+\code{AbsBiasSS = 1 - AbsBias_exp / AbsBias_ref}. The statistical significance
+is obtained based on a Random Walk test at the 95% confidence level (DelSole 
+and Tippett, 2016). If there is more than one dataset, the result will be 
+computed for each pair of exp and obs data.
+}
+\examples{
+exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+ref <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
+biasSS1 <- AbsBiasSS(exp = exp, obs = obs, ref = ref, memb_dim = 'member')
+biasSS2 <- AbsBiasSS(exp = exp, obs = obs, ref = NULL, memb_dim = 'member')
+
+}
+\references{
+Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+DelSole and Tippett, 2016; https://doi.org/10.1175/MWR-D-15-0218.1
+}
diff --git a/man/Bias.Rd b/man/Bias.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2a02f2d52b9d39b8df0c7b8fa06ef02d702b1b11
--- /dev/null
+++ b/man/Bias.Rd
@@ -0,0 +1,73 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Bias.R
+\name{Bias}
+\alias{Bias}
+\title{Compute the Mean Bias}
+\usage{
+Bias(
+  exp,
+  obs,
+  time_dim = "sdate",
+  memb_dim = NULL,
+  dat_dim = NULL,
+  na.rm = FALSE,
+  absolute = FALSE,
+  time_mean = TRUE,
+  ncores = NULL
+)
+}
+\arguments{
+\item{exp}{A named numerical array of the forecast with at least time 
+dimension.}
+
+\item{obs}{A named numerical array of the observation with at least time 
+dimension. The dimensions must be the same as 'exp' except 'memb_dim' and 
+'dat_dim'.}
+
+\item{time_dim}{A character string indicating the name of the time dimension.
+The default value is 'sdate'.}
+
+\item{memb_dim}{A character string indicating the name of the member dimension
+to compute the ensemble mean; it should be set to NULL if the parameter 
+'exp' is already the ensemble mean. The default value is NULL.}
+
+\item{dat_dim}{A character string indicating the name of dataset dimension. 
+The length of this dimension can be different between 'exp' and 'obs'. 
+The default value is NULL.}
+
+\item{na.rm}{A logical value indicating if NAs should be removed (TRUE) or
+kept (FALSE) for computation. The default value is FALSE.}
+
+\item{absolute}{A logical value indicating whether to compute the absolute 
+bias. The default value is FALSE.}
+
+\item{time_mean}{A logical value indicating whether to compute the temporal 
+mean of the bias. The default value is TRUE.}
+
+\item{ncores}{An integer indicating the number of cores to use for parallel 
+computation. The default value is NULL.}
+}
+\value{
+A numerical array of bias with dimensions c(nexp, nobs, the rest dimensions of
+'exp' except 'time_dim' (if time_mean = T) and 'memb_dim'). 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). If dat_dim is NULL, nexp and nobs are omitted.
+}
+\description{
+The Mean Bias or Mean Error (Wilks, 2011) is defined as the mean difference 
+between the ensemble mean forecast and the observations. It is a deterministic
+metric. Positive values indicate that the forecasts are on average too high
+and negative values indicate that the forecasts are on average too low.
+It also allows to compute the Absolute Mean Bias or bias without temporal 
+mean. If there is more than one dataset, the result will be computed for each
+pair of exp and obs data.
+}
+\examples{
+exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, sdate = 50))
+obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, sdate = 50))
+bias <- Bias(exp = exp, obs = obs, memb_dim = 'member')
+
+}
+\references{
+Wilks, 2011; https://doi.org/10.1016/B978-0-12-385022-5.00008-7
+}
diff --git a/tests/testthat/test-AbsBiasSS.R b/tests/testthat/test-AbsBiasSS.R
new file mode 100644
index 0000000000000000000000000000000000000000..b2e70f3e73b9aa1cd0b1ca350b8e3c62e6200c04
--- /dev/null
+++ b/tests/testthat/test-AbsBiasSS.R
@@ -0,0 +1,319 @@
+context("s2dv::AbsBiasSS tests")
+
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+set.seed(3)
+ref1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(60), dim = c(member = 2, sdate = 10, lat = 3))
+set.seed(2)
+obs2 <- array(rnorm(30), dim = c(member = 1, sdate = 10, lat = 3))
+set.seed(3)
+ref2 <- array(rnorm(60), dim = c(member = 2, sdate = 10, lat = 3))
+
+# dat3
+set.seed(1)
+exp3 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2))
+set.seed(2)
+obs3 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3))
+set.seed(3)
+ref3 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat4
+set.seed(1)
+exp4 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2))
+set.seed(2)
+obs4 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3))
+set.seed(3)
+ref4 <- array(rnorm(40), dim = c(sdate = 10, lat = 2, dataset = 2))
+
+# dat5
+set.seed(1)
+exp5 <- array(rnorm(10), dim = c(sdate = 10))
+exp5_1 <- array(exp5, dim = c(sdate = 10, dataset = 1))
+exp5_2 <- exp5_1
+exp5_2[1] <- NA
+set.seed(2)
+obs5 <- array(rnorm(10), dim = c(sdate = 10))
+obs5_1 <- array(obs5, dim = c(sdate = 10, dataset = 1))
+obs5_2 <- obs5_1
+obs5_2[c(1, 3)] <- NA
+set.seed(3)
+ref5 <- array(rnorm(10), dim = c(sdate = 10))
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs (1)
+  expect_error(
+    AbsBiasSS(c()),
+    "Parameter 'exp' must be a numeric array."
+  )
+  expect_error(
+    AbsBiasSS(exp1, c()),
+    "Parameter 'obs' must be a numeric array."
+  )
+  expect_error(
+    AbsBiasSS(exp1, array(rnorm(20), dim = c(10, lat = 2))),
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+    AbsBiasSS(exp1, array(rnorm(20), dim = c(10, lat = 2))),
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, ref = 'ref'),
+    "Parameter 'ref' must be a numeric array."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, ref = array(rnorm(20), dim = c(10, lat = 2))),
+    "Parameter 'ref' must have dimension names."
+  )
+  # time_dim
+  expect_error(
+    AbsBiasSS(exp1, obs1, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, time_dim = 'time'),
+    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, ref = array(rnorm(20), dim = c(lon = 10, lat = 2)), time_dim = 'time'),
+    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # memb_dim
+  expect_error(
+    AbsBiasSS(exp1, obs1, memb_dim = TRUE),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, memb_dim = 'memb'),
+    "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  expect_error(
+    AbsBiasSS(exp2, array(rnorm(20), dim = c(member = 3, sdate = 10, lat = 2)), memb_dim = 'member'),
+    "Not implemented for observations with members ('obs' can have 'memb_dim', but it should be of length = 1).", fixed = TRUE
+  )
+  # dat_dim
+  expect_error(
+    AbsBiasSS(exp1, obs1, dat_dim = TRUE, ),
+    "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, dat_dim = 'dat_dim'),
+    "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension. Set it as NULL if there is no dataset dimension."
+  )
+  # exp, ref, and obs (2)
+  expect_error(
+    AbsBiasSS(exp1, array(1:9, dim = c(sdate = 9))),
+    "Parameter 'exp' and 'obs' must have same length of all dimensions except 'memb_dim' and 'dat_dim'."
+  )
+  expect_error(
+    AbsBiasSS(exp3, obs3, ref = obs3, dat_dim = 'dataset', memb_dim = 'member'),
+    "If parameter 'ref' has dataset dimension, it must be equal to dataset dimension of 'exp'."
+  )
+  expect_error(
+    AbsBiasSS(exp1, obs1, ref = ref2),
+    "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."
+  )
+  # na.rm
+  expect_error(
+    AbsBiasSS(exp2, obs2, memb_dim = 'member', na.rm = 1.5),
+    "Parameter 'na.rm' must be one logical value."
+  )
+  # ncores
+  expect_error(
+    AbsBiasSS(exp2, obs2, memb_dim = 'member', ncores = 1.5),
+    "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+  expect_equal(
+    dim(AbsBiasSS(exp1, obs1)$biasSS),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp1, obs1)$sign),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp1, obs1, ref1)$biasSS),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp1, obs1, ref1)$sign),
+    c(lat = 2)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp1, obs1)$biasSS),
+    c(-0.3103594, 0.0772921),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp1, obs1, ref1)$biasSS),
+    c(-0.07871642, 0.28868904),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp1, obs1)$sign),
+    c(FALSE, FALSE),
+    tolerance = 0.0001
+  )
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+  expect_equal(
+    dim(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS),
+    c(lat = 3)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp2, obs2, ref2, memb_dim = 'member')$biasSS),
+    c(lat = 3)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS),
+    c(-0.4743706, -0.2884077, -0.4064404),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp2, obs2, ref2, memb_dim = 'member')$biasSS),
+    c(-0.07319869, 0.06277502, -0.17321998),
+    tolerance = 0.0001
+  )
+
+})
+
+##############################################
+test_that("4. Output checks: dat3", {
+
+  expect_equal(
+    dim(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS),
+    c(nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset')$biasSS),
+    c(nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10],
+    c(-0.09794912, -0.11814710, -0.28840768, 0.02117376, -0.31282805, -0.08754112),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10],
+    c(0.264602089, 0.251073637, 0.006772886, 0.245427687, 0.168998894, 0.311602249),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS)[1:2],
+    as.vector(AbsBiasSS(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')$biasSS[1,1,])
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp3, obs3, ref3, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[1:5],
+    c(-0.21373395, -0.34444526, 0.11039962, 0.05523797, 0.26460209)
+  )
+
+})
+##############################################
+test_that("5. Output checks: dat4", {
+
+  expect_equal(
+    dim(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset')$biasSS),
+    c(nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset')$biasSS),
+    c(nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[5:10],
+    c(-0.09794912, -0.11814710, -0.28840768, 0.02117376, -0.31282805, -0.08754112),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset')$biasSS)[5:10],
+    c(0.264602089, 0.133379718, 0.006772886, 0.242372951, 0.168998894, 0.272767238),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp2, obs2, memb_dim = 'member')$biasSS)[1:2],
+    as.vector(AbsBiasSS(exp4, obs4, memb_dim = 'member', dat_dim = 'dataset')$biasSS[1,1,])
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp4, obs4, ref4, memb_dim = 'member', dat_dim = 'dataset', na.rm = TRUE)$biasSS)[1:5],
+    c(-0.213733950, -0.214240924, 0.110399615, -0.009733463, 0.264602089)
+  )
+
+})
+
+##############################################
+test_that("6. Output checks: dat5", {
+  expect_equal(
+    dim(AbsBiasSS(exp5, obs5)$biasSS),
+    NULL
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp5, obs5)$sign),
+    NULL
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5, obs5)$biasSS),
+    -0.3103594,
+    tolerance = 0.0001
+  )
+  expect_equal( 
+    as.vector(AbsBiasSS(exp5, obs5)$sign),
+    FALSE
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5, obs5, ref5)$biasSS),
+    -0.07871642,
+    tolerance = 0.0001
+  )
+  #5_1
+  expect_equal(
+    dim(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$biasSS),
+    c(nexp = 1, nobs = 1)
+  )
+  expect_equal(
+    dim(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$sign),
+    c(nexp = 1, nobs = 1)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5_1, obs5_1, dat_dim = 'dataset')$biasSS),
+    as.vector(AbsBiasSS(exp5, obs5)$biasSS)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5_1, obs5_1, ref = ref5, dat_dim = 'dataset')$biasSS),
+    as.vector(AbsBiasSS(exp5, obs5, ref = ref5)$biasSS)
+  )
+  #5_2: NA test
+  expect_equal(
+    as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset')$biasSS),
+    as.numeric(NA)
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5_1, obs5_2, dat_dim = 'dataset', na.rm = T)$biasSS),
+    c(-0.4636772),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(AbsBiasSS(exp5_2, obs5_2, ref = ref5, dat_dim = 'dataset', na.rm = T)$biasSS),
+    c(0.08069355),
+    tolerance = 0.0001
+  )
+})
diff --git a/tests/testthat/test-Bias.R b/tests/testthat/test-Bias.R
new file mode 100644
index 0000000000000000000000000000000000000000..842ecc2c2b626649603f3d964ca20adee2fd22df
--- /dev/null
+++ b/tests/testthat/test-Bias.R
@@ -0,0 +1,240 @@
+context("s2dv::Bias tests")
+
+##############################################
+
+# dat1
+set.seed(1)
+exp1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+set.seed(2)
+obs1 <- array(rnorm(20), dim = c(sdate = 10, lat = 2))
+
+# dat2
+set.seed(1)
+exp2 <- array(rnorm(40), dim = c(member = 2, sdate = 10, lat = 2))
+set.seed(2)
+obs2 <- array(rnorm(20), dim = c(member = 1, sdate = 10, lat = 2))
+
+# dat3
+set.seed(1)
+exp3 <- array(rnorm(80), dim = c(member = 2, sdate = 10, lat = 2, dataset = 2))
+set.seed(2)
+obs3 <- array(rnorm(60), dim = c(sdate = 10, lat = 2, dataset = 3))
+
+# dat4
+set.seed(1)
+exp4 <- array(rnorm(10), dim = c(sdate = 10))
+exp4_1 <- array(exp4, dim = c(sdate = 10, dataset = 1))
+set.seed(2)
+obs4 <- array(rnorm(10), dim = c(sdate = 10))
+obs4_1 <- array(obs4, dim = c(sdate = 10, dataset = 1))
+obs4_2 <- obs4_1
+obs4_2[c(1, 3)] <- NA
+
+##############################################
+test_that("1. Input checks", {
+  # exp and obs (1)
+  expect_error(
+    Bias(c()),
+    "Parameter 'exp' must be a numeric array."
+  )
+  expect_error(
+    Bias(exp1, c()),
+    "Parameter 'obs' must be a numeric array."
+  )
+  expect_error(
+    Bias(exp1, array(rnorm(20), dim = c(10, lat = 2))),
+    "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  # time_dim
+  expect_error(
+    Bias(exp1, obs1, time_dim = 1),
+    "Parameter 'time_dim' must be a character string."
+  )
+  expect_error(
+    Bias(exp1, obs1, time_dim = 'time'),
+    "Parameter 'time_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  # memb_dim
+  expect_error(
+    Bias(exp1, obs1, memb_dim = TRUE),
+    "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+    Bias(exp1, obs1, memb_dim = 'memb'),
+    "Parameter 'memb_dim' is not found in 'exp' dimension."
+  )
+  expect_error(
+    Bias(exp2, array(rnorm(20), dim = c(member = 3, sdate = 10, lat = 2)), memb_dim = 'member'),
+    "Not implemented for observations with members ('obs' can have 'memb_dim', but it should be of length = 1).", fixed = TRUE
+  )
+  # dat_dim
+  expect_error(
+    Bias(exp1, obs1, dat_dim = TRUE, ),
+    "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+    Bias(exp1, obs1, dat_dim = 'dat_dim'),
+    "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension. Set it as NULL if there is no dataset dimension."
+  )
+  # exp, ref, and obs (2)
+  expect_error(
+    Bias(exp1, array(1:9, dim = c(sdate = 9))),
+    "Parameter 'exp' and 'obs' must have same length of all dimensions except 'memb_dim' and 'dat_dim'."
+  )
+  # na.rm
+  expect_error(
+    Bias(exp2, obs2, memb_dim = 'member', na.rm = 1.5),
+    "Parameter 'na.rm' must be one logical value."
+  )
+  # absolute
+  expect_error(
+    Bias(exp2, obs2, memb_dim = 'member', ncores = 1.5),
+    "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+  # time_mean
+  expect_error(
+    Bias(exp2, obs2, memb_dim = 'member', time_mean = 1.5),
+    "Parameter 'time_mean' must be one logical value."
+  )
+  # ncores
+  expect_error(
+    Bias(exp2, obs2, memb_dim = 'member', ncores = 1.5),
+    "Parameter 'ncores' must be either NULL or a positive integer."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+  expect_equal(
+    dim(Bias(exp1, obs1)),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(Bias(exp1, obs1, time_mean = FALSE)),
+    c(sdate = 10, lat = 2)
+  )
+  expect_equal(
+    as.vector(Bias(exp1, obs1)),
+    c(-0.07894886, 0.06907455),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(Bias(exp1, obs1, absolute = TRUE)),
+    c(0.9557288, 0.8169118),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(Bias(exp1, obs1, time_mean = FALSE, na.rm = TRUE))[1:5],
+    c(0.27046074, -0.00120586, -2.42347394, 2.72565648, 0.40975953),
+    tolerance = 0.0001
+  )
+
+})
+
+##############################################
+test_that("3. Output checks: dat2", {
+
+  expect_equal(
+    dim(Bias(exp2, obs2, memb_dim = 'member')),
+    c(lat = 2)
+  )
+  expect_equal(
+    dim(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)),
+    c(sdate = 10, lat = 2)
+  )
+  expect_equal(
+    as.vector(Bias(exp2, obs2, memb_dim = 'member')),
+    c(-0.02062777, -0.18624194),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)[1:2,1:2]),
+    c(0.6755093, 0.1949769, 0.4329061, -1.9391461),
+    tolerance = 0.0001
+  )
+
+})
+
+##############################################
+test_that("4. Output checks: dat3", {
+
+  expect_equal(
+    dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')),
+    c(nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    dim(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)),
+    c(sdate = 10, nexp = 2, nobs = 3, lat = 2)
+  )
+  expect_equal(
+    as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset'))[5:10],
+    c(0.23519286, 0.18346575, -0.18624194, -0.07803352, 0.28918537, 0.39739379),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', absolute = TRUE, time_mean = FALSE))[5:10],
+    c(0.2154482, 0.8183919, 2.1259250, 0.7796967, 1.5206510, 0.8463483),
+    tolerance = 0.0001
+  )
+  expect_equal(
+    as.vector(Bias(exp2, obs2, memb_dim = 'member')),
+    as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset')[1,1,])
+  )
+  expect_equal(
+    as.vector(Bias(exp2, obs2, memb_dim = 'member', time_mean = FALSE)),
+    as.vector(Bias(exp3, obs3, memb_dim = 'member', dat_dim = 'dataset', time_mean = FALSE)[ ,1,1,])
+  )
+
+})
+
+##############################################
+test_that("5. Output checks: dat4", {
+  expect_equal(
+    dim(Bias(exp4, obs4)),
+    NULL
+  )
+  expect_equal(
+    dim(Bias(exp4, obs4, time_mean = F)),
+    c(sdate = 10)
+  )
+  expect_equal(
+    as.vector(Bias(exp4, obs4, time_mean = F)),
+    as.vector(exp4 - obs4)
+  )
+  expect_equal(
+    as.vector(Bias(exp4, obs4, time_mean = F, absolute = T)),
+    abs(as.vector(exp4 - obs4))
+  )
+  expect_equal(
+    as.vector(Bias(exp4, obs4, absolute = T)),
+    mean(abs(as.vector(exp4 - obs4)))
+  )
+
+  expect_equal(
+    dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset')),
+    c(nexp = 1, nobs = 1)
+  )
+  expect_equal(
+    dim(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)),
+    c(sdate = 10, nexp = 1, nobs = 1)
+  )
+  expect_equal(
+    as.vector(Bias(exp4_1, obs4_1, dat_dim = 'dataset', time_mean = F)),
+    as.vector(Bias(exp4, obs4, time_mean = F))
+  )
+  # 4_2: NA
+  expect_equal(
+    as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset')),
+    as.numeric(NA)
+  )
+  expect_equal(
+    as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F))[c(1, 3)],
+    as.numeric(c(NA, NA))
+  )
+  expect_equal(
+    as.vector(Bias(exp4_1, obs4_2, dat_dim = 'dataset', time_mean = F))[c(2, 4:10)],
+    as.vector(Bias(exp4, obs4, time_mean = F))[c(2, 4:10)]
+  )
+})