# exp <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, member = 10, year = 50))
# obs <- array(rnorm(1000), dim = c(dat = 1, lat = 3, lon = 5, year = 50))
# bias <- Bias(exp = exp, obs = obs, time_dim = 'year', memb_dim = 'member', na.rm = FALSE, ncores = 1)

Bias <- function(exp, obs, time_dim = 'sdate', memb_dim = NULL, na.rm = FALSE, ncores = NULL) {
  
  ## Checks
  
  
  ## Ensemble mean
  if (!is.null(memb_dim)) {
    exp <- MeanDims(exp, memb_dim, na.rm = na.rm)
  }
  
  ## Mean bias
  bias <- multiApply::Apply(data = list(exp, obs),
                            target_dims = time_dim,
                            fun = .Bias,
                            ncores = ncores)$output1
  ## Return the mean bias
  bias <- MeanDims(bias, time_dim, na.rm = na.rm)
  
  return(bias)
}

.Bias <- function(exp, obs) {
  
  bias <- exp - obs
  
  return(bias)
}