diff --git a/NAMESPACE b/NAMESPACE
index 6da8d0c7ee74b297aab584a94a5431777f714d19..182ac412402ec1912df026d42ee4821531ce9500 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(ACC)
 export(AMV)
 export(AnimateMap)
 export(Ano)
@@ -27,6 +28,7 @@ export(LeapYear)
 export(Load)
 export(MeanDims)
 export(Persistence)
+export(PlotACC)
 export(PlotAno)
 export(PlotClim)
 export(PlotEquiMap)
@@ -91,6 +93,8 @@ importFrom(stats,pf)
 importFrom(stats,pt)
 importFrom(stats,qchisq)
 importFrom(stats,qnorm)
+importFrom(stats,qt)
+importFrom(stats,quantile)
 importFrom(stats,rnorm)
 importFrom(stats,sd)
 importFrom(stats,ts)
diff --git a/R/ACC.R b/R/ACC.R
new file mode 100644
index 0000000000000000000000000000000000000000..bf3daa27a2bc4017bc1c7f444671fda7adc9c700
--- /dev/null
+++ b/R/ACC.R
@@ -0,0 +1,618 @@
+#'Compute the anomaly correlation coefficient between the forecast and corresponding observation
+#'
+#'Calculate the anomaly correlation coefficient for the ensemble mean of each 
+#'model and the corresponding references over a spatial domain. It can return a
+#'forecast time series if the data contain forest time dimension, and also the 
+#'start date mean if the data contain start date dimension. 
+#'The domain of interest can be specified by providing the list 
+#'of longitudes/latitudes (lon/lat) of the data together with the corners 
+#'of the domain: lonlatbox = c(lonmin, lonmax, latmin, latmax).
+#'
+#'@param exp A numeric array of experimental anomalies with named dimensions.
+#'  It must have at least 'dat_dim' and 'space_dim'.
+#'@param obs A numeric array of observational anomalies with named dimensions.
+#'  It must have the same dimensions as 'exp' except the length of 'dat_dim' 
+#'  and 'memb_dim'.
+#'@param dat_dim A character string indicating the name of dataset (nobs/nexp) 
+#'  dimension. The default value is 'dataset'.
+#'@param space_dim A character string vector of 2 indicating the name of the
+#'  latitude and longitude dimensions (in order) along which ACC is computed. 
+#'  The default value is c('lat', 'lon').
+#'@param avg_dim A character string indicating the name of the dimension to be
+#'  averaged. It must be one of 'time_dim'. The mean ACC is calculated along 
+#'  averaged. If no need to calculate mean ACC, set as NULL. The default value 
+#'  is 'sdate'.
+#'@param memb_dim A character string indicating the name of the member 
+#'  dimension. If the data are not ensemble ones, set as NULL. The default 
+#'  value is 'member'.
+#'@param lat A vector of the latitudes of the exp/obs grids. Only required when
+#'  the domain of interested is specified. The default value is NULL.
+#'@param lon A vector of the longitudes of the exp/obs grids. Only required when
+#'  the domain of interested is specified. The default value is NULL.
+#'@param lonlatbox A numeric vector of 4 indicating the corners of the domain of
+#'  interested: c(lonmin, lonmax, latmin, latmax). Only required when the domain
+#'  of interested is specified. The default value is NULL.
+#'@param conf A logical value indicating whether to retrieve the confidence 
+#'  intervals or not. The default value is TRUE.
+#'@param conftype A charater string of "parametric" or "bootstrap". 
+#'  "parametric" provides a confidence interval for the ACC computed by a 
+#'  Fisher transformation and a significance level for the ACC from a one-sided
+#'  student-T distribution. "bootstrap" provides a confidence interval for the
+#'  ACC and MACC computed from bootstrapping on the members with 100 drawings 
+#'  with replacement. To guarantee the statistical robustness of the result, 
+#'  make sure that your experiment and observation always have the same number
+#'  of members. "bootstrap" requires 'memb_dim' has value. The default value is
+#'  'parametric'.
+#'@param conf.lev A numeric indicating the confidence level for the 
+#'  regression computation. The default value is 0.95.
+#'@param pval A logical value indicating whether to compute the p-value or not.
+#'  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 list containing the numeric arrays:\cr 
+#'\item{acc}{
+#'  The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
+#'  space_dim 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).
+#'}
+#'\item{conf.lower (if conftype = "parametric") or acc_conf.lower (if 
+#'      conftype = "bootstrap")}{
+#'  The lower confidence interval of ACC with the same dimensions as ACC. Only
+#'  present if \code{conf = TRUE}.
+#'}
+#'\item{conf.upper (if conftype = "parametric") or acc_conf.upper (if 
+#'      conftype = "bootstrap")}{
+#'  The upper confidence interval of ACC with the same dimensions as ACC. Only 
+#'  present if \code{conf = TRUE}.
+#'}
+#'\item{p.val}{
+#'  The p-value with the same dimensions as ACC. Only present if 
+#'  \code{pval = TRUE} and code{conftype = "parametric"}.  
+#'}
+#'\item{macc}{
+#'  The mean anomaly correlation coefficient with dimensions
+#'  c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
+#'  avg_dim). Only present if 'avg_dim' is not NULL.
+#'}
+#'\item{macc_conf.lower}{
+#'  The lower confidence interval of MACC with the same dimensions as MACC. 
+#'  Only present if \code{conftype = "bootstrap"}.
+#'}
+#'\item{macc_conf.upper}{
+#'  The upper confidence interval of MACC with the same dimensions as MACC. 
+#'  Only present if \code{conftype = "bootstrap"}.
+#'}
+#'
+#'@examples
+#'  \dontshow{
+#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+#'                                                c('observation'), startDates,
+#'                                                leadtimemin = 1,
+#'                                                leadtimemax = 4,
+#'                                                output = 'lonlat',
+#'                                                latmin = 27, latmax = 48,
+#'                                                lonmin = -12, lonmax = 40)
+#'  }
+#'sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+#'sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2) 
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+#'acc <- ACC(ano_exp, ano_obs)
+#'acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+#'# Combine acc results for PlotACC
+#'res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
+#'             dim = c(dim(acc$acc), 4))
+#'res_bootstrap <- array(c(acc$acc_conf.lower, acc$acc, acc$acc_conf.upper, acc$p.val),
+#'                       dim = c(dim(acc$acc), 4))
+#'  \donttest{
+#'PlotACC(res, startDates)
+#'PlotACC(res_bootstrap, startDates)
+#'  }
+#'@references Joliffe and Stephenson (2012). Forecast Verification: A 
+#'  Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.
+#'@import multiApply
+#'@importFrom abind abind
+#'@importFrom stats qt qnorm quantile
+#'@importFrom ClimProjDiags Subset
+#'@export
+ACC <- function(exp, obs, dat_dim = 'dataset', space_dim = c('lat', 'lon'),
+                avg_dim = 'sdate', memb_dim = 'member', 
+                lat = NULL, lon = NULL, lonlatbox = NULL, 
+                conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE,
+                ncores = NULL) {
+
+  # Check inputs 
+  ## exp and obs (1)
+  if (is.null(exp) | is.null(obs)) {
+    stop("Parameter 'exp' and 'obs' cannot be NULL.")
+  }
+  if (!is.numeric(exp) | !is.numeric(obs)) {
+    stop("Parameter 'exp' and 'obs' must be a numeric array.")
+  }
+  if (is.null(dim(exp)) | is.null(dim(obs))) {
+    stop(paste0("Parameter 'exp' and 'obs' must have at least dimensions ",
+                "dat_dim and space_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)) {
+    stop("Parameter 'exp' and 'obs' must have dimension names.")
+  }
+  if(!all(names(dim(exp)) %in% names(dim(obs))) |
+     !all(names(dim(obs)) %in% names(dim(exp)))) {
+    stop("Parameter 'exp' and 'obs' must have same dimension names.")
+  }
+  ## 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.")
+  }
+  ## space_dim
+  if (!is.character(space_dim) | length(space_dim) != 2) {
+    stop("Parameter 'space_dim' must be a character vector of 2.")
+  }
+  if (any(!space_dim %in% names(dim(exp))) | any(!space_dim %in% names(dim(obs)))) {
+    stop("Parameter 'space_dim' is not found in 'exp' or 'obs' dimension.")
+  }
+  ## avg_dim
+  if (!is.null(avg_dim)) {
+    if (!is.character(avg_dim) | length(avg_dim) > 1) {
+      stop("Parameter 'avg_dim' must be a character string.")
+    } 
+    if (!avg_dim %in% names(dim(exp)) | !avg_dim %in% names(dim(obs))) {
+      stop("Parameter 'avg_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)) | !memb_dim %in% names(dim(obs))) {
+      stop("Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension.")
+    }
+  }
+  ## lat
+  if (!is.null(lat)) {
+    if (!is.numeric(lat) | length(lat) != dim(exp)[space_dim[1]]) {
+      stop(paste0("Parameter 'lat' must be a numeric vector with the same ",
+                  "length as the latitude dimension of 'exp' and 'obs'."))
+    }
+  }
+  ## lon
+  if (!is.null(lon)) {
+    if (!is.numeric(lon) | length(lon) != dim(exp)[space_dim[2]]) {
+      stop(paste0("Parameter 'lon' must be a numeric vector with the same ",
+                  "length as the longitude dimension of 'exp' and 'obs'."))
+    }
+  }
+  ## lonlatbox
+  if (!is.null(lonlatbox)) {
+    if (!is.numeric(lonlatbox) | length(lonlatbox) != 4) {
+      stop("Parameter 'lonlatbox' must be a numeric vector of 4.")
+    }
+  }
+  ## lat, lon, and lonlatbox
+  if (!is.null(lon) & !is.null(lat) & !is.null(lonlatbox)) {
+    select_lonlat <- TRUE
+  } else if (is.null(lon) & is.null(lat) & is.null(lonlatbox)) {
+    select_lonlat <- FALSE
+  } else {
+    stop(paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
+                "NULL at the same time."))
+  }
+  ## conf
+  if (!is.logical(conf) | length(conf) > 1) {
+    stop("Parameter 'conf' must be one logical value.")
+  }
+  if (conf) {
+    ## conftype 
+    if (!conftype %in% c('parametric', 'bootstrap')) {
+      stop("Parameter 'conftype' must be either 'parametric' or 'bootstrap'.")
+    }
+    if (conftype == 'bootstrap' & is.null(memb_dim)) {
+      stop("Parameter 'memb_dim' cannot be NULL when parameter 'conftype' is 'bootstrap'.")
+    }
+    ## conf.lev
+    if (!is.numeric(conf.lev) | conf.lev < 0 | conf.lev > 1 | length(conf.lev) > 1) {
+      stop("Parameter 'conf.lev' must be a numeric number between 0 and 1.")
+    }
+  }
+  ## pval
+  if (!is.logical(pval) | length(pval) > 1) {
+    stop("Parameter 'pval' 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 a positive integer.")
+    }
+  }
+  ## exp and obs (2)
+  name_exp <- sort(names(dim(exp)))
+  name_obs <- sort(names(dim(obs)))
+  name_exp <- name_exp[-which(name_exp == dat_dim)]
+  name_obs <- name_obs[-which(name_obs == dat_dim)]
+  if (!is.null(memb_dim)) {
+    name_exp <- name_exp[-which(name_exp == memb_dim)]
+    name_obs <- name_obs[-which(name_obs == memb_dim)]
+  }
+  if(!all(dim(exp)[name_exp] == dim(obs)[name_obs])) {
+    stop(paste0("Parameter 'exp' and 'obs' must have same length of ",
+                "all the dimensions expect 'dat_dim' and 'memb_dim'."))
+  }
+
+#-----------------------------------------------------------------
+
+
+  ###############################
+  # Sort dimension
+  name_exp <- names(dim(exp))
+  name_obs <- names(dim(obs))
+  order_obs <- match(name_exp, name_obs)
+  obs <- Reorder(obs, order_obs)
+  ###############################
+
+  # Select the domain 
+  if (select_lonlat) {
+    for (jind in 1:2) {
+      while (lonlatbox[jind] < 0) {
+        lonlatbox[jind] <- lonlatbox[jind] + 360
+      }
+      while (lonlatbox[jind] > 360) {
+        lonlatbox[jind] <- lonlatbox[jind] - 360
+      }
+    }
+    indlon <- which((lon >= lonlatbox[1] & lon <= lonlatbox[2]) | 
+                    (lonlatbox[1] > lonlatbox[2] & (lon > lonlatbox[1] | lon < lonlatbox[2])))
+    indlat <- which(lat >= lonlatbox[3] & lat <= lonlatbox[4])
+
+    exp <- ClimProjDiags::Subset(exp, space_dim, list(indlat, indlon), drop = FALSE)
+    obs <- ClimProjDiags::Subset(obs, space_dim, list(indlat, indlon), drop = FALSE)
+  }
+
+   # Ensemble mean
+  if (!is.null(memb_dim)) {
+    if (conftype == 'bootstrap') {
+      exp_ori <- exp
+      obs_ori <- obs
+    }
+    exp <- MeanDims(exp, memb_dim, na.rm = TRUE, ncores = ncores)
+    obs <- MeanDims(obs, memb_dim, na.rm = TRUE, ncores = ncores)
+  }
+
+  if (is.null(avg_dim)) {
+    res <- Apply(list(exp, obs),
+                 target_dims = list(c(space_dim, dat_dim),
+                                    c(space_dim, dat_dim)),
+                 fun = .ACC,
+                 dat_dim = dat_dim, avg_dim = avg_dim,
+                 conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
+                 ncores_input = ncores,
+                 ncores = ncores)
+
+    if (conftype == 'bootstrap') {
+      res_conf <- Apply(list(exp_ori, obs_ori),
+                 target_dims = list(c(memb_dim, dat_dim, space_dim),
+                                    c(memb_dim, dat_dim, space_dim)),
+                 fun = .ACC_bootstrap,
+                 dat_dim = dat_dim, memb_dim = memb_dim, avg_dim = avg_dim,
+                 conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
+                 ncores_input = ncores,
+                 ncores = ncores)
+      #NOTE: pval?
+      res <- list(acc = res$acc, 
+                  acc_conf.lower = res_conf$acc_conf.lower, 
+                  acc_conf.upper = res_conf$acc_conf.upper,
+                  macc = res$macc,
+                  macc_conf.lower = res_conf$macc_conf.lower, 
+                  macc_conf.upper = res_conf$macc_conf.upper)
+    }
+
+  } else {
+    res <- Apply(list(exp, obs),
+                 target_dims = list(c(space_dim, avg_dim, dat_dim),
+                                    c(space_dim, avg_dim, dat_dim)),
+                 fun = .ACC,
+                 dat_dim = dat_dim, avg_dim = avg_dim,
+                 conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
+                 ncores_input = ncores,
+                 ncores = ncores)
+
+    if (conftype == 'bootstrap') {
+      res_conf <- Apply(list(exp_ori, obs_ori),
+                        target_dims = list(c(memb_dim, dat_dim, avg_dim, space_dim),
+                                           c(memb_dim, dat_dim, avg_dim, space_dim)),
+                        fun = .ACC_bootstrap,
+                        dat_dim = dat_dim, memb_dim = memb_dim, avg_dim = avg_dim,
+                        conftype = conftype, pval = pval, conf = conf, conf.lev = conf.lev,
+                        ncores_input = ncores,
+                        ncores = ncores)
+      res <- list(acc = res$acc, 
+                  acc_conf.lower = res_conf$acc_conf.lower, 
+                  acc_conf.upper = res_conf$acc_conf.upper,
+                  macc = res$macc,
+                  macc_conf.lower = res_conf$macc_conf.lower,
+                  macc_conf.upper = res_conf$macc_conf.upper)
+
+    }
+
+  }
+
+ return(res)
+}
+
+.ACC <- function(exp, obs,  dat_dim = 'dataset', #space_dim = c('lat', 'lon'),
+                 avg_dim = 'sdate', #memb_dim = NULL,
+                 lon = NULL, lat = NULL, lonlatbox = NULL,
+                 conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE,
+                 ncores_input = NULL) {
+
+# if (is.null(avg_dim)) 
+  # exp: [space_dim, dat_exp]
+  # obs: [space_dim, dat_obs]
+# if (!is.null(avg_dim)) 
+  # exp: [space_dim, avg_dim, dat_exp]
+  # obs: [space_dim, avg_dim, dat_obs]
+
+  # .ACC() should use all the spatial points to calculate ACC. It returns [nexp, nobs].
+
+  nexp <- as.numeric(dim(exp)[length(dim(exp))])
+  nobs <- as.numeric(dim(obs)[length(dim(obs))])
+
+  if (is.null(avg_dim)) {
+    acc <- array(dim = c(nexp = nexp, nobs = nobs))
+    if (pval) p.val <- array(dim = c(nexp = nexp, nobs = nobs))
+    if (conf) {
+      conf.upper <- array(dim = c(nexp = nexp, nobs = nobs))
+      conf.lower <- array(dim = c(nexp = nexp, nobs = nobs))
+      if (conftype == 'bootstrap') {
+         ndraw <- 100
+         acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw)) 
+      }
+    }
+    
+  } else {
+        acc <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
+        names(dim(acc))[3] <- avg_dim
+        macc <- array(dim = c(nexp = nexp, nobs = nobs))
+        if (pval) p.val <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
+    if (conf) {
+        conf.upper <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
+        conf.lower <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1]))
+      if (conftype == 'bootstrap') {
+         ndraw <- 100
+         acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[length(dim(exp)) - 1], ndraw))
+         macc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw))
+      }
+    }
+  }  
+  
+  # Per-paired exp and obs. NAs should be in the same position in both exp and obs
+  for (iobs in 1:nobs) {
+    for (iexp in 1:nexp) {
+      exp_sub <- ClimProjDiags::Subset(exp, dat_dim, iexp, drop = 'selected')
+      obs_sub <- ClimProjDiags::Subset(obs, dat_dim, iobs, drop = 'selected')
+      # dim: [space_dim]
+
+      # Variance(iexp) should not take into account any point 
+      # that is not available in iobs and therefore not accounted for 
+      # in covariance(iexp, iobs) and vice-versa 
+      exp_sub[is.na(obs_sub)] <- NA
+      obs_sub[is.na(exp_sub)] <- NA
+
+      if (is.null(avg_dim)) {
+        # ACC
+        top <- sum(exp_sub*obs_sub, na.rm = TRUE)  #a number
+        bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
+        acc[iexp, iobs] <- top/bottom #a number
+        # handle bottom = 0
+        if (is.infinite(acc[iexp, iobs])) acc[iexp, iobs] <- NA
+        # pval and conf
+        if (pval | conf) {
+          if (conftype == "parametric") {
+            # calculate effective sample size along space_dim
+            # combine space_dim into one dim first
+            obs_tmp <- array(obs_sub, dim = c(space = length(obs_sub)))
+            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a number
+            if (pval) {
+              t <- qt(conf.lev, eno - 2)  # a number
+              p.val[iexp, iobs] <- sqrt(t^2 / (t^2 + eno - 2))
+            }
+            if (conf) {
+              conf.upper[iexp, iobs] <- tanh(atanh(acc[iexp, iobs]) + qnorm(1 - (1 - conf.lev) / 2) / sqrt(eno - 3))
+              conf.lower[iexp, iobs] <- tanh(atanh(acc[iexp, iobs]) + qnorm((1 - conf.lev) / 2) / sqrt(eno - 3))
+            }
+          }
+        }
+
+      } else {  #avg_dim is not NULL
+        # MACC
+        top <- sum(exp_sub*obs_sub, na.rm = TRUE)  #a number
+        bottom <- sqrt(sum(exp_sub^2, na.rm = TRUE) * sum(obs_sub^2, na.rm = TRUE))
+        macc[iexp, iobs] <- top/bottom #a number
+        # handle bottom = 0
+        if (is.infinite(macc[iexp, iobs])) macc[iexp, iobs] <- NA
+        # ACC
+        for (i in 1:dim(acc)[3]) {   #NOTE: use sapply!!!
+          exp_sub_i <- ClimProjDiags::Subset(exp_sub, avg_dim, i, drop = 'selected')
+          obs_sub_i <- ClimProjDiags::Subset(obs_sub, avg_dim, i, drop = 'selected')
+          #dim: [space_dim]
+          top <- sum(exp_sub_i*obs_sub_i, na.rm = TRUE)  #a number
+          bottom <- sqrt(sum(exp_sub_i^2, na.rm = TRUE) * sum(obs_sub_i^2, na.rm = TRUE))
+          acc[iexp, iobs, i] <- top/bottom #a number
+          # handle bottom = 0
+          if (is.infinite(acc[iexp, iobs, i])) acc[iexp, iobs, i] <- NA
+        }
+
+        # pval and conf
+        if (pval | conf) {
+          if (conftype == "parametric") {
+            # calculate effective sample size along space_dim 
+            # combine space_dim into one dim first
+            obs_tmp <- array(obs_sub, dim = c(space = prod(dim(obs_sub)[-length(dim(obs_sub))]), 
+                                              dim(obs_sub)[length(dim(obs_sub))]))
+            eno <- Eno(obs_tmp, 'space', ncores = ncores_input)  # a vector of avg_dim
+            if (pval) {
+              t <- qt(conf.lev, eno - 2)  # a vector of avg_dim
+              p.val[iexp, iobs, ] <- sqrt(t^2 / (t^2 + eno - 2))
+            }
+            if (conf) {
+              conf.upper[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + qnorm(1 - (1 - conf.lev) / 2) / sqrt(eno - 3))
+              conf.lower[iexp, iobs, ] <- tanh(atanh(acc[iexp, iobs, ]) + qnorm((1 - conf.lev) / 2) / sqrt(eno - 3))
+            }
+          }
+        }
+
+      }  # if avg_dim is not NULL
+
+    }
+  }
+
+#------------------------------------------------
+  
+
+
+  # Return output
+  if (is.null(avg_dim)) {
+    if (conf & pval) {
+      return(list(acc = acc, conf.lower = conf.lower, conf.upper = conf.upper,
+                  p.val = p.val))
+    } else if (conf & !pval) {
+      return(list(acc = acc, conf.lower = conf.lower, conf.upper = conf.upper,
+                  macc = macc))
+    } else if (!conf & pval) {
+      return(list(acc = acc, p.val = p.val))
+    } else {
+      return(list(acc = acc))
+    }
+  } else {
+    if (conf & pval) {
+      return(list(acc = acc, conf.lower = conf.lower, conf.upper = conf.upper, 
+                  p.val = p.val, macc = macc))
+    } else if (conf & !pval) {
+      return(list(acc = acc, conf.lower = conf.lower, conf.upper = conf.upper,
+                  macc = macc))
+    } else if (!conf & pval) {
+      return(list(acc = acc, p.val = p.val, macc = macc))
+    } else {
+      return(list(acc = acc, macc = macc))
+    }
+  }
+
+}
+
+
+.ACC_bootstrap <- function(exp, obs,  dat_dim = 'dataset', #space_dim = c('lat', 'lon'),
+                           avg_dim = 'sdate', memb_dim = NULL,
+                           lon = NULL, lat = NULL, lonlatbox = NULL,
+                           conf = TRUE, conftype = "parametric", conf.lev = 0.95, pval = TRUE, 
+                           ncores_input = NULL) {
+# if (is.null(avg_dim)) 
+  # exp: [memb_exp, dat_exp, space_dim]
+  # obs: [memb_obs, dat_obs, space_dim]
+# if (!is.null(avg_dim)) 
+  # exp: [memb_exp, dat_exp, avg_dim, space_dim]
+  # obs: [memb_obs, dat_obs, avg_dim, space_dim]
+
+  nexp <- as.numeric(dim(exp)[2])
+  nobs <- as.numeric(dim(obs)[2])
+  nmembexp <- as.numeric(dim(exp)[1])
+  nmembobs <- as.numeric(dim(obs)[1])
+
+  ndraw <- 100
+  if (is.null(avg_dim)) {
+    acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw))
+  } else {
+    acc_draw <- array(dim = c(nexp = nexp, nobs = nobs, dim(exp)[3], ndraw))
+    macc_draw <- array(dim = c(nexp = nexp, nobs = nobs, ndraw))
+  }
+
+  for (jdraw in 1:ndraw) {
+    #choose a randomly member index for each point of the matrix 
+    indexp <- array(sample(nmembexp, size = prod(dim(exp)[-c(length(dim(exp)) - 1, length(dim(exp)))]),
+                           replace = TRUE), 
+                    dim = dim(exp))
+    indobs <- array(sample(nmembobs, size = prod(dim(obs)[-c(length(dim(obs)) - 1, length(dim(obs)))]),               
+                           replace = TRUE),
+                     dim = dim(obs))
+
+      #combine maxtrix of data and random index
+      varindexp <- abind::abind(exp, indexp, along = length(dim(exp)) + 1)
+      varindobs <- abind::abind(obs, indobs, along = length(dim(obs)) + 1)
+
+    #select randomly the members for each point of the matrix
+#    if (is.null(avg_dim)) {
+
+    drawexp <- array( 
+                    apply(varindexp, c(2:length(dim(exp))), function(x) x[,1][x[,2]] ),
+                          dim = dim(exp)) 
+    drawobs <- array(
+                    apply(varindobs, c(2:length(dim(obs))), function(x) x[,1][x[,2]] ),
+                          dim = dim(obs))
+
+    # ensemble mean before .ACC
+    drawexp <- MeanDims(drawexp, memb_dim, na.rm = TRUE, ncores = ncores_input)
+    drawobs <- MeanDims(drawobs, memb_dim, na.rm = TRUE, ncores = ncores_input)
+    # Reorder
+    if (is.null(avg_dim)) {
+      drawexp <- Reorder(drawexp, c(2, 3, 1))
+      drawobs <- Reorder(drawobs, c(2, 3, 1))
+    } else {
+      drawexp <- Reorder(drawexp, c(3, 4, 2, 1))
+      drawobs <- Reorder(drawobs, c(3, 4, 2, 1))
+    } 
+   
+    #calculate the ACC of the randomized field
+    tmpACC <- .ACC(drawexp, drawobs, conf = FALSE, pval = FALSE, avg_dim = avg_dim,
+                   ncores_input = ncores_input)
+    if (is.null(avg_dim)) {
+      acc_draw[, , jdraw] <- tmpACC$acc
+    } else {
+      acc_draw[, , , jdraw] <- tmpACC$acc
+      macc_draw[, , jdraw] <- tmpACC$macc
+    }
+  }
+
+  #calculate the confidence interval
+  if (is.null(avg_dim)) {
+  acc_conf.upper <- apply(acc_draw, c(1, 2),
+                          function (x) {
+                            quantile(x, 1 - (1 - conf.lev) / 2, na.rm = TRUE)})
+  acc_conf.lower <- apply(acc_draw, c(1, 2),
+                          function (x) {
+                            quantile(x, (1 - conf.lev) / 2, na.rm = TRUE)})
+
+  } else {
+  acc_conf.upper <- apply(acc_draw, c(1, 2, 3), 
+                          function (x) {
+                            quantile(x, 1 - (1 - conf.lev) / 2, na.rm = TRUE)})
+  acc_conf.lower <- apply(acc_draw, c(1, 2, 3),
+                          function (x) {
+                            quantile(x, (1 - conf.lev) / 2, na.rm = TRUE)})
+  macc_conf.upper <- apply(macc_draw, c(1, 2),
+                          function (x) {
+                            quantile(x, 1 - (1 - conf.lev) / 2, na.rm = TRUE)})
+  macc_conf.lower <- apply(macc_draw, c(1, 2),
+                          function (x) {
+                            quantile(x, (1 - conf.lev) / 2, na.rm = TRUE)})
+  }
+
+  # Return output
+  if (is.null(avg_dim)) {
+      return(list(acc_conf.lower = acc_conf.lower,
+                  acc_conf.upper = acc_conf.upper))
+  } else {
+      return(list(acc_conf.lower = acc_conf.lower,
+                  acc_conf.upper = acc_conf.upper,
+                  macc_conf.lower = macc_conf.lower,
+                  macc_conf.upper = macc_conf.upper))
+  }
+                  
+}
diff --git a/R/PlotACC.R b/R/PlotACC.R
new file mode 100644
index 0000000000000000000000000000000000000000..3bffa68f279234a7d6881eb1848830df6f7ad3d9
--- /dev/null
+++ b/R/PlotACC.R
@@ -0,0 +1,238 @@
+#'Plot Plumes/Timeseries Of Anomaly Correlation Coefficients
+#'
+#'Plots plumes/timeseries of ACC from an array with dimensions 
+#'(output from \code{ACC()}): \cr
+#'c(nexp, nobs, nsdates, nltime, 4)\cr
+#'where the fourth dimension is of length 4 and contains the lower limit of 
+#'the 95\% confidence interval, the ACC, the upper limit of the 95\% 
+#'confidence interval and the 95\% significance level given by a one-sided 
+#'T-test.
+#'
+#'@param ACC An ACC array with with dimensions:\cr
+#'  c(nexp, nobs, nsdates, nltime, 4)\cr
+#'  with the fourth dimension of length 4 containing the lower limit of the 
+#'  95\% confidence interval, the ACC, the upper limit of the 95\% confidence 
+#'  interval and the 95\% significance level.
+#'@param sdates A character vector of startdates: c('YYYYMMDD','YYYYMMDD').
+#'@param toptitle A character string of the main title, optional.
+#'@param sizetit A multiplicative factor to scale title size, optional.
+#'@param ytitle A character string of the title of Y-axis for each experiment:
+#'  c('', ''), optional.
+#'@param limits A numeric vector c(lower limit, upper limit): limits of the 
+#'  Y-axis, optional.
+#'@param legends A character vector of flags to be written in the legend, 
+#'  optional.
+#'@param freq A integer: 1 = yearly, 12 = monthly, 4 = seasonal, ... Default: 12.
+#'@param biglab A logical value for presentation/paper plot, Default = FALSE.
+#'@param fill A logical value if filled confidence interval. Default = FALSE.
+#'@param linezero A logical value if a line at y=0 should be added. Default = FALSE.
+#'@param points A logical value if points instead of lines. Default = TRUE.\cr
+#'  Must be TRUE if only 1 leadtime.
+#'@param vlines A vector of x location where to add vertical black lines, optional.
+#'@param fileout A character string of the output file name. Extensions allowed: 
+#'  eps/ps, jpeg, png, pdf, bmp and tiff. Default is NULL.
+#'@param width A numeric of the file width, in the units specified in the 
+#'  parameter size_units (inches by default). Takes 8 by default.
+#'@param height A numeric of the file height, in the units specified in the parameter 
+#'  size_units (inches by default). Takes 5 by default.
+#'@param size_units A character string of the units of the size of the device 
+#'  (file or window) to plot in. Inches ('in') by default. See ?Devices and the 
+#'  creator function of the corresponding device.
+#'@param res Resolution of the device (file or window) to plot in. See 
+#'  ?Devices and the creator function of the corresponding device.
+#'@param \dots Arguments to be passed to the method. Only accepts the following
+#'  graphical parameters:\cr
+#'  adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+#'  csi cxy err family fg fig fin font font.axis font.lab font.main font.sub 
+#'  lend lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page 
+#'  plt smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog\cr
+#'  For more information about the parameters see `par`.
+#'
+#'@examples
+#'  \dontshow{
+#'startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+#'sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+#'                                                c('observation'), startDates,
+#'                                                leadtimemin = 1,
+#'                                                leadtimemax = 4,
+#'                                                output = 'lonlat',
+#'                                                latmin = 27, latmax = 48,
+#'                                                lonmin = -12, lonmax = 40)
+#'  }
+#'sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+#'sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2) 
+#'clim <- Clim(sampleData$mod, sampleData$obs)
+#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+
+#'acc <- ACC(ano_exp, ano_obs)
+#'acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+#'# Combine acc results for PlotACC
+#'res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
+#'             dim = c(dim(acc$acc), 4))
+#'res_bootstrap <- array(c(acc$acc_conf.lower, acc$acc, acc$acc_conf.upper, acc$p.val),
+#'                       dim = c(dim(acc$acc), 4))
+#'  \donttest{
+#'PlotACC(res, startDates)
+#'PlotACC(res_bootstrap, startDates)
+#'  }
+#'@importFrom grDevices dev.cur dev.new dev.off 
+#'@importFrom stats ts
+#'@export
+PlotACC <- function(ACC, sdates, toptitle = "", sizetit = 1, ytitle = "", 
+                    limits = NULL, legends = NULL, freq = 12, biglab = FALSE, 
+                    fill = FALSE, linezero = FALSE, points = TRUE, vlines = NULL, 
+                    fileout = NULL, 
+                    width = 8, height = 5, size_units = 'in', res = 100, ...) {
+  # Process the user graphical parameters that may be passed in the call
+  ## Graphical parameters to exclude
+  excludedArgs <- c("cex", "cex.axis", "cex.lab", "cex.main", "col", "lab", "las", "lty", "lwd", "mai", "mgp", "new", "pch", "pin", "ps", "pty")
+  userArgs <- .FilterUserGraphicArgs(excludedArgs, ...)
+
+  # If there is any filenames to store the graphics, process them
+  # to select the right device 
+  if (!is.null(fileout)) {
+    deviceInfo <- .SelectDevice(fileout = fileout, width = width, height = height, units = size_units, res = res)
+    saveToFile <- deviceInfo$fun
+    fileout <- deviceInfo$files
+  }
+
+  #
+  if (length(dim(ACC)) != 5 | dim(ACC)[5] != 4) {
+    stop("5 dim needed : c(nexp, nobs, nsdates, nltime, 4)")
+  }
+  nexp <- dim(ACC)[1]
+  nobs <- dim(ACC)[2]
+  nleadtime <- dim(ACC)[4]
+  nsdates <- dim(ACC)[3]
+  if (is.null(limits) == TRUE) {
+    ll <- min(ACC, na.rm = TRUE)
+    ul <- max(ACC, na.rm = TRUE)
+    if (biglab) {
+      ul <- ul + 0.3 * (ul - ll)
+    } else {
+      ul <- ul + 0.2 * (ul - ll)
+    }
+  } else {
+    ll <- limits[1]
+    ul <- limits[2]
+  }
+  yearinit <- as.integer(substr(sdates[1], 1, 4))
+  moninit <- as.integer(substr(sdates[1], 5, 6))
+  lastyear <- as.integer(substr(sdates[nsdates], 1, 4)) + (moninit + (
+                         nleadtime - 1) * 12 / freq - 1) %/% 12
+  lastmonth <- (moninit + (nleadtime - 1) * (12 / freq) - 1) %% 12 + 1
+  empty_ts <- ts(start = c(yearinit, (moninit - 1) %/% (12 / freq) + 1), 
+                 end = c(lastyear, (lastmonth - 1) %/% (12 / freq) + 1), 
+                 frequency = freq)
+  color <- c("red4", "dodgerblue4", "lightgoldenrod4", "deeppink4", 
+             "mediumpurple4", "green4", "orange4", "lightblue4", "mediumorchid4", 
+             "olivedrab4")
+  colorblock <- c("red1", "dodgerblue1", "lightgoldenrod1", "deeppink1", 
+                  "mediumpurple1", "green1", "orange1", "lightblue1", 
+                  "mediumorchid1", "olivedrab1")
+
+  # Open connection to graphical device
+  if (!is.null(fileout)) {
+    saveToFile(fileout)
+  } else if (names(dev.cur()) == 'null device') {
+    dev.new(units = size_units, res = res, width = width, height = height)
+  }
+
+  # Load the user parameters
+  par(userArgs)
+
+  if (biglab) {
+    par(mai = c(1, 1.1, 0.5, 0), mgp = c(2.8, 0.9, 0))
+    par(cex = 1.3, cex.lab = 2, cex.axis = 1.8)
+    cexmain <- 2.2
+    legsize <- 1.5
+  } else {
+    par(mai = c(0.8, 0.8, 0.5, 0.1), mgp = c(2, 0.5, 0))
+    par(cex = 1.3, cex.lab = 1.5, cex.axis = 1.1)
+    cexmain <- 1.5
+    legsize <- 1
+  }
+  plot(empty_ts, ylim = c(ll, ul), xlab = "Time (years)", ylab = ytitle, 
+       main = toptitle, cex.main = cexmain * sizetit)
+  for (jexp in 1:nexp) {
+    for (jobs in 1:nobs) {
+      numcol <- jobs + (jexp - 1) * nobs
+      for (jdate in 1:nsdates) {
+        year0 <- as.integer(substr(sdates[jdate], 1, 4))
+        mon0 <- as.integer(substr(sdates[jdate], 5, 6))
+        start <- (year0 - yearinit) * freq + 1
+        end <- start + nleadtime - 1
+        var <- array(dim = c(3, length(empty_ts)))
+        var[, start:end] <- t(ACC[jexp, jobs, jdate, , 1:3])
+        if (fill) {
+          par(new = TRUE)
+          bordup <- ACC[jexp, jobs, jdate, , 3]
+          borddown <- ACC[jexp, jobs, jdate, , 1]
+          tmp <- c(start:end)
+          xout <- is.na(bordup + borddown)
+          tmp <- tmp[which(xout == FALSE)]
+          xx <- c(tmp, rev(tmp))
+          bordup <- bordup[which(xout == FALSE)]
+          borddown <- borddown[which(xout == FALSE)]
+          yy <- c(bordup, rev(borddown))
+          if (jdate == 1) {
+            matplot(t(var), type = "l", lty = 1, lwd = 1, ylim = c(ll, ul), 
+                    col = color[numcol], xlab = "", ylab = "", axes = FALSE)
+          }
+          polygon(xx, yy, col = colorblock[numcol], border = NA)
+        }
+        if (points) {
+          par(new = TRUE)
+          plot(var[2, ], type = "p", lty = 1, lwd = 6, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE,
+               cex = 0.6)
+          par(new = TRUE)
+          plot(var[1, ], type = "p", pch = 6, lwd = 3, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE,
+               cex = 0.6)
+          par(new = TRUE)
+          plot(var[3, ], type = "p", pch = 2, lwd = 3, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE,
+               cex = 0.6)
+          par(new = TRUE)
+          for (jind in start:end) {
+            lines(c(jind, jind), var[c(1, 3), jind], lwd = 1, 
+                  ylim = c(ll, ul), col = color[numcol], xlab = "", 
+                  ylab = "", axes = FALSE)
+          }
+        } else {
+          par(new = TRUE)
+          plot(var[2, ], type = "l", lty = 1, lwd = 4, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE)
+          par(new = TRUE)
+          plot(var[1, ], type = "l", lty = 1, lwd = 1, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE)
+          par(new = TRUE)
+          plot(var[3, ], type = "l", lty = 1, lwd = 1, ylim = c(ll, ul), 
+               col = color[numcol], xlab = "", ylab = "", axes = FALSE)
+        }
+      }
+    }
+  }
+  if (linezero) {
+    abline(h = 0, col = "black")
+  }
+  if (is.null(vlines) == FALSE) {
+    for (x in vlines) {
+      abline(v = x, col = "black")
+    }
+  }
+  if (is.null(legends) == FALSE) {
+    if (points) {
+      legend(0, ul, legends[1:(nobs * nexp)], lty = 3, lwd = 10, 
+             col = color[1:(nobs * nexp)], cex = legsize)
+    } else {
+      legend(0, ul, legends[1:(nobs * nexp)], lty = 1, lwd = 4, 
+             col = color[1:(nobs * nexp)], cex = legsize)
+    }
+  }
+
+  # If the graphic was saved to file, close the connection with the device
+  if(!is.null(fileout)) dev.off()
+}
diff --git a/man/ACC.Rd b/man/ACC.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..d48d6b8fa01aca912ecd6ca5dfac319d0402a2f5
--- /dev/null
+++ b/man/ACC.Rd
@@ -0,0 +1,144 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ACC.R
+\name{ACC}
+\alias{ACC}
+\title{Compute the anomaly correlation coefficient between the forecast and corresponding observation}
+\usage{
+ACC(exp, obs, dat_dim = "dataset", space_dim = c("lat", "lon"),
+  avg_dim = "sdate", memb_dim = "member", lat = NULL, lon = NULL,
+  lonlatbox = NULL, conf = TRUE, conftype = "parametric",
+  conf.lev = 0.95, pval = TRUE, ncores = NULL)
+}
+\arguments{
+\item{exp}{A numeric array of experimental anomalies with named dimensions.
+It must have at least 'dat_dim' and 'space_dim'.}
+
+\item{obs}{A numeric array of observational anomalies with named dimensions.
+It must have the same dimensions as 'exp' except the length of 'dat_dim' 
+and 'memb_dim'.}
+
+\item{dat_dim}{A character string indicating the name of dataset (nobs/nexp) 
+dimension. The default value is 'dataset'.}
+
+\item{space_dim}{A character string vector of 2 indicating the name of the
+latitude and longitude dimensions (in order) along which ACC is computed. 
+The default value is c('lat', 'lon').}
+
+\item{avg_dim}{A character string indicating the name of the dimension to be
+averaged. It must be one of 'time_dim'. The mean ACC is calculated along 
+averaged. If no need to calculate mean ACC, set as NULL. The default value 
+is 'sdate'.}
+
+\item{memb_dim}{A character string indicating the name of the member 
+dimension. If the data are not ensemble ones, set as NULL. The default 
+value is 'member'.}
+
+\item{lat}{A vector of the latitudes of the exp/obs grids. Only required when
+the domain of interested is specified. The default value is NULL.}
+
+\item{lon}{A vector of the longitudes of the exp/obs grids. Only required when
+the domain of interested is specified. The default value is NULL.}
+
+\item{lonlatbox}{A numeric vector of 4 indicating the corners of the domain of
+interested: c(lonmin, lonmax, latmin, latmax). Only required when the domain
+of interested is specified. The default value is NULL.}
+
+\item{conf}{A logical value indicating whether to retrieve the confidence 
+intervals or not. The default value is TRUE.}
+
+\item{conftype}{A charater string of "parametric" or "bootstrap". 
+"parametric" provides a confidence interval for the ACC computed by a 
+Fisher transformation and a significance level for the ACC from a one-sided
+student-T distribution. "bootstrap" provides a confidence interval for the
+ACC and MACC computed from bootstrapping on the members with 100 drawings 
+with replacement. To guarantee the statistical robustness of the result, 
+make sure that your experiment and observation always have the same number
+of members. "bootstrap" requires 'memb_dim' has value. The default value is
+'parametric'.}
+
+\item{conf.lev}{A numeric indicating the confidence level for the 
+regression computation. The default value is 0.95.}
+
+\item{pval}{A logical value indicating whether to compute the p-value or not.
+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 list containing the numeric arrays:\cr 
+\item{acc}{
+ The ACC with the dimensions c(nexp, nobs, the rest of the dimension except 
+ space_dim 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).
+}
+\item{conf.lower (if conftype = "parametric") or acc_conf.lower (if 
+     conftype = "bootstrap")}{
+ The lower confidence interval of ACC with the same dimensions as ACC. Only
+ present if \code{conf = TRUE}.
+}
+\item{conf.upper (if conftype = "parametric") or acc_conf.upper (if 
+     conftype = "bootstrap")}{
+ The upper confidence interval of ACC with the same dimensions as ACC. Only 
+ present if \code{conf = TRUE}.
+}
+\item{p.val}{
+ The p-value with the same dimensions as ACC. Only present if 
+ \code{pval = TRUE} and code{conftype = "parametric"}.  
+}
+\item{macc}{
+ The mean anomaly correlation coefficient with dimensions
+ c(nexp, nobs, the rest of the dimension except space_dim, memb_dim, and 
+ avg_dim). Only present if 'avg_dim' is not NULL.
+}
+\item{macc_conf.lower}{
+ The lower confidence interval of MACC with the same dimensions as MACC. 
+ Only present if \code{conftype = "bootstrap"}.
+}
+\item{macc_conf.upper}{
+ The upper confidence interval of MACC with the same dimensions as MACC. 
+ Only present if \code{conftype = "bootstrap"}.
+}
+}
+\description{
+Calculate the anomaly correlation coefficient for the ensemble mean of each 
+model and the corresponding references over a spatial domain. It can return a
+forecast time series if the data contain forest time dimension, and also the 
+start date mean if the data contain start date dimension. 
+The domain of interest can be specified by providing the list 
+of longitudes/latitudes (lon/lat) of the data together with the corners 
+of the domain: lonlatbox = c(lonmin, lonmax, latmin, latmax).
+}
+\examples{
+ \dontshow{
+startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+                                               c('observation'), startDates,
+                                               leadtimemin = 1,
+                                               leadtimemax = 4,
+                                               output = 'lonlat',
+                                               latmin = 27, latmax = 48,
+                                               lonmin = -12, lonmax = 40)
+ }
+sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2) 
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+acc <- ACC(ano_exp, ano_obs)
+acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+# Combine acc results for PlotACC
+res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
+            dim = c(dim(acc$acc), 4))
+res_bootstrap <- array(c(acc$acc_conf.lower, acc$acc, acc$acc_conf.upper, acc$p.val),
+                      dim = c(dim(acc$acc), 4))
+ \donttest{
+PlotACC(res, startDates)
+PlotACC(res_bootstrap, startDates)
+ }
+}
+\references{
+Joliffe and Stephenson (2012). Forecast Verification: A 
+ Practitioner's Guide in Atmospheric Science. Wiley-Blackwell.
+}
+
diff --git a/man/PlotACC.Rd b/man/PlotACC.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..1dbd7cf6579e33460c0c79f83ad5731914ebd8cb
--- /dev/null
+++ b/man/PlotACC.Rd
@@ -0,0 +1,109 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/PlotACC.R
+\name{PlotACC}
+\alias{PlotACC}
+\title{Plot Plumes/Timeseries Of Anomaly Correlation Coefficients}
+\usage{
+PlotACC(ACC, sdates, toptitle = "", sizetit = 1, ytitle = "",
+  limits = NULL, legends = NULL, freq = 12, biglab = FALSE,
+  fill = FALSE, linezero = FALSE, points = TRUE, vlines = NULL,
+  fileout = NULL, width = 8, height = 5, size_units = "in", res = 100,
+  ...)
+}
+\arguments{
+\item{ACC}{An ACC array with with dimensions:\cr
+c(nexp, nobs, nsdates, nltime, 4)\cr
+with the fourth dimension of length 4 containing the lower limit of the 
+95\% confidence interval, the ACC, the upper limit of the 95\% confidence 
+interval and the 95\% significance level.}
+
+\item{sdates}{A character vector of startdates: c('YYYYMMDD','YYYYMMDD').}
+
+\item{toptitle}{A character string of the main title, optional.}
+
+\item{sizetit}{A multiplicative factor to scale title size, optional.}
+
+\item{ytitle}{A character string of the title of Y-axis for each experiment:
+c('', ''), optional.}
+
+\item{limits}{A numeric vector c(lower limit, upper limit): limits of the 
+Y-axis, optional.}
+
+\item{legends}{A character vector of flags to be written in the legend, 
+optional.}
+
+\item{freq}{A integer: 1 = yearly, 12 = monthly, 4 = seasonal, ... Default: 12.}
+
+\item{biglab}{A logical value for presentation/paper plot, Default = FALSE.}
+
+\item{fill}{A logical value if filled confidence interval. Default = FALSE.}
+
+\item{linezero}{A logical value if a line at y=0 should be added. Default = FALSE.}
+
+\item{points}{A logical value if points instead of lines. Default = TRUE.\cr
+Must be TRUE if only 1 leadtime.}
+
+\item{vlines}{A vector of x location where to add vertical black lines, optional.}
+
+\item{fileout}{A character string of the output file name. Extensions allowed: 
+eps/ps, jpeg, png, pdf, bmp and tiff. Default is NULL.}
+
+\item{width}{A numeric of the file width, in the units specified in the 
+parameter size_units (inches by default). Takes 8 by default.}
+
+\item{height}{A numeric of the file height, in the units specified in the parameter 
+size_units (inches by default). Takes 5 by default.}
+
+\item{size_units}{A character string of the units of the size of the device 
+(file or window) to plot in. Inches ('in') by default. See ?Devices and the 
+creator function of the corresponding device.}
+
+\item{res}{Resolution of the device (file or window) to plot in. See 
+?Devices and the creator function of the corresponding device.}
+
+\item{\dots}{Arguments to be passed to the method. Only accepts the following
+graphical parameters:\cr
+adj ann ask bg bty cex.sub cin col.axis col.lab col.main col.sub cra crt 
+csi cxy err family fg fig fin font font.axis font.lab font.main font.sub 
+lend lheight ljoin lmitre mar mex mfcol mfrow mfg mkh oma omd omi page 
+plt smo srt tck tcl usr xaxp xaxs xaxt xlog xpd yaxp yaxs yaxt ylbias ylog\cr
+For more information about the parameters see `par`.}
+}
+\description{
+Plots plumes/timeseries of ACC from an array with dimensions 
+(output from \code{ACC()}): \cr
+c(nexp, nobs, nsdates, nltime, 4)\cr
+where the fourth dimension is of length 4 and contains the lower limit of 
+the 95\% confidence interval, the ACC, the upper limit of the 95\% 
+confidence interval and the 95\% significance level given by a one-sided 
+T-test.
+}
+\examples{
+ \dontshow{
+startDates <- c('19851101', '19901101', '19951101', '20001101', '20051101')
+sampleData <- s2dv:::.LoadSampleData('tos', c('experiment'),
+                                               c('observation'), startDates,
+                                               leadtimemin = 1,
+                                               leadtimemax = 4,
+                                               output = 'lonlat',
+                                               latmin = 27, latmax = 48,
+                                               lonmin = -12, lonmax = 40)
+ }
+sampleData$mod <- Season(sampleData$mod, monini = 11, moninf = 12, monsup = 2)
+sampleData$obs <- Season(sampleData$obs, monini = 11, moninf = 12, monsup = 2) 
+clim <- Clim(sampleData$mod, sampleData$obs)
+ano_exp <- Ano(sampleData$mod, clim$clim_exp)
+ano_obs <- Ano(sampleData$obs, clim$clim_obs)
+acc <- ACC(ano_exp, ano_obs)
+acc_bootstrap <- ACC(ano_exp, ano_obs, conftype = 'bootstrap')
+# Combine acc results for PlotACC
+res <- array(c(acc$conf.lower, acc$acc, acc$conf.upper, acc$p.val), 
+            dim = c(dim(acc$acc), 4))
+res_bootstrap <- array(c(acc$acc_conf.lower, acc$acc, acc$acc_conf.upper, acc$p.val),
+                      dim = c(dim(acc$acc), 4))
+ \donttest{
+PlotACC(res, startDates)
+PlotACC(res_bootstrap, startDates)
+ }
+}
+
diff --git a/tests/testthat/test-ACC.R b/tests/testthat/test-ACC.R
new file mode 100644
index 0000000000000000000000000000000000000000..7631401d0ae4864e0aa418daa1e9b116542719cc
--- /dev/null
+++ b/tests/testthat/test-ACC.R
@@ -0,0 +1,199 @@
+context("s2dv::ACC tests")
+
+##############################################
+##NOTE: bootstrap is not tested because sample() is used inside.
+
+  # dat1
+  set.seed(1)
+  exp1 <- array(rnorm(60), dim = c(dataset = 1, member = 2, sdate = 5,
+                                   ftime = 1, lat = 2, lon = 3))
+
+  set.seed(2)
+  obs1 <- array(rnorm(30), dim = c(dataset = 1, member = 1, sdate = 5,
+                                   ftime = 1, lat = 2, lon = 3))
+  # dat2
+  set.seed(1)
+  exp2 <- array(rnorm(60), dim = c(dataset = 2, sdate = 5,
+                                   ftime = 1, lat = 2, lon = 3))
+  set.seed(2)
+  obs2 <- array(rnorm(30), dim = c(dataset = 1, sdate = 5,
+                                   ftime = 1, lat = 2, lon = 3))
+  set.seed(2)
+  na <- floor(runif(2, min = 1, max = 30))
+  obs2[na] <- NA
+
+##############################################
+test_that("1. Input checks", {
+
+  expect_error(
+  ACC(c(), c()),
+  "Parameter 'exp' and 'obs' cannot be NULL."
+  )
+  expect_error(
+  ACC(c('b'), c('a')),
+  "Parameter 'exp' and 'obs' must be a numeric array."
+  )
+  expect_error(
+  ACC(c(1:10), c(2:4)),
+  paste0("Parameter 'exp' and 'obs' must have at least dimensions ",
+         "dat_dim and space_dim.")
+  )
+  expect_error(
+  ACC(array(1:10, dim = c(2, 5)), array(1:4, dim = c(2, 2))),
+  "Parameter 'exp' and 'obs' must have dimension names."
+  )
+  expect_error(
+  ACC(array(1:10, dim = c(a = 2, c = 5)), array(1:4, dim = c(a = 2, b = 2))),
+  "Parameter 'exp' and 'obs' must have same dimension name"
+  )
+  expect_error(
+  ACC(exp1, obs1, dat_dim = 1),
+  "Parameter 'dat_dim' must be a character string."
+  )
+  expect_error(
+  ACC(exp1, obs1, dat_dim = 'a'),
+  "Parameter 'dat_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+  ACC(exp1, obs1, space_dim = c('lon')),
+  "Parameter 'space_dim' must be a character vector of 2."
+  )
+  expect_error(
+  ACC(exp1, obs1, space_dim = c('lon', 'lev')),
+  "Parameter 'space_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+  ACC(exp1, obs1, avg_dim = 1),
+  "Parameter 'avg_dim' must be a character string."
+  )
+  expect_error(
+  ACC(exp1, obs1, avg_dim = c('lev')),
+  "Parameter 'avg_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+  ACC(exp1, obs1, memb_dim = TRUE),
+  "Parameter 'memb_dim' must be a character string."
+  )
+  expect_error(
+  ACC(exp1, obs1, memb_dim = 'memb'),
+  "Parameter 'memb_dim' is not found in 'exp' or 'obs' dimension."
+  )
+  expect_error(
+  ACC(exp1, obs1, lat = c(1, 2, 3)),
+  paste0("Parameter \'lat\' must be a numeric vector with the same ",
+         "length as the latitude dimension of \'exp\' and \'obs\'.")
+  )
+  expect_error(
+  ACC(exp1, obs1, lon = c(1, 3)),
+  paste0("Parameter \'lon\' must be a numeric vector with the same ",
+         "length as the longitude dimension of \'exp\' and \'obs\'.")
+  )
+  expect_error(
+  ACC(exp1, obs1, lonlatbox = c(-90, 90)),
+  "Parameter 'lonlatbox' must be a numeric vector of 4."
+  )
+  expect_error(
+  ACC(exp1, obs1, lat = c(1, 2), lon = c(1, 2, 3)),
+  paste0("Parameters 'lon', 'lat', and 'lonlatbox' must be used or be ",
+         "NULL at the same time.")
+  )
+  expect_error(
+  ACC(exp1, obs1, conf = 1),
+  "Parameter 'conf' must be one logical value."
+  )
+  expect_error(
+  ACC(exp1, obs1, conftype = 'a'),
+  "Parameter \'conftype\' must be either \'parametric\' or \'bootstrap\'."
+  )
+  expect_error(
+  ACC(exp1, obs1, memb_dim = NULL, conftype = 'bootstrap'),
+  "Parameter 'memb_dim' cannot be NULL when parameter 'conftype' is 'bootstrap'."
+  )
+  expect_error(
+  ACC(exp1, obs1, conf.lev = -1),
+  "Parameter 'conf.lev' must be a numeric number between 0 and 1."
+  )
+  expect_error(
+  ACC(exp1, obs1, pval = 'TRUE'),
+  "Parameter 'pval' must be one logical value."
+  )
+  expect_error(
+  ACC(exp1, obs1, ncores = 1.5),
+  "Parameter 'ncores' must be a positive integer."
+  )
+  expect_error(
+  ACC(exp = array(1:10, dim = c(dataset = 5, member = 1, lat = 2, lon = 1)),
+      obs = array(1:4, dim = c(dataset = 2, member = 2, lat = 1, lon = 1)),
+      avg_dim = NULL),
+  "Parameter 'exp' and 'obs' must have same length of all the dimensions expect 'dat_dim' and 'memb_dim'."
+  )
+
+})
+
+##############################################
+test_that("2. Output checks: dat1", {
+
+  expect_equal(
+    dim(ACC(exp1, obs1)$acc),
+    c(nexp = 1, nobs = 1, sdate = 5, ftime = 1)
+  )
+  expect_equal(
+    names(ACC(exp1, obs1)),
+    c("acc", "conf.lower", "conf.upper", "p.val", "macc")
+  )
+  expect_equal(
+    mean(ACC(exp1, obs1)$acc),
+    -0.2001352,
+    tolerance = 0.00001
+  )
+  expect_equal(
+    as.vector(ACC(exp1, obs1)$p.val),   
+    c(0.7292993, 0.7292993, 0.7292993, 0.7292993, 0.7292993),
+    tolerance = 0.00001    
+  )
+  expect_equal(
+    as.vector(ACC(exp1, obs1)$conf.lower),
+    c(-0.8595534, -0.9644555, -0.9408508, -0.6887500, -0.7619374),
+    tolerance = 0.00001
+  )
+  expect_equal(
+    as.vector(ACC(exp1, obs1)$conf.upper),
+    c(0.7493799, 0.2515608, 0.4759707, 0.8890967, 0.8517117),
+    tolerance = 0.00001
+  )
+  expect_equal(
+    names(ACC(exp1, obs1, avg_dim = NULL)),
+    c("acc", "conf.lower", "conf.upper", "p.val")
+  )
+  expect_equal(
+    dim(ACC(exp1, obs1, dat_dim = 'member', memb_dim = NULL)$acc),
+    c(nexp = 2, nobs = 1, sdate = 5, dataset = 1, ftime = 1)
+  )
+  expect_equal(
+    names(ACC(exp1, obs1, conf = FALSE)),
+    c("acc", "p.val", "macc")
+  )
+  expect_equal(
+    names(ACC(exp1, obs1, pval = FALSE)),
+    c("acc", "conf.lower", "conf.upper", "macc")
+  )
+  expect_equal(
+    names(ACC(exp1, obs1, conf = FALSE, pval = FALSE)),
+    c("acc", "macc")
+  )
+  expect_equal(
+    as.vector(ACC(exp1, obs1, conf = FALSE, avg_dim = NULL, conf.lev = 0.9)$p.val),
+    c(0.6083998, 0.6083998, 0.6083998, 0.6083998, 0.6083998),
+    tolerance = 0.00001
+  )
+  expect_equal(
+    mean(ACC(exp1, obs1, lat = c(10, 20), lon = c(20, 30, 40), lonlatbox = c(20, 30, 10, 20))$acc),
+    -0.1681097,
+    tolerance = 0.00001
+  )
+})
+
+
+
+
+