#'Compute forecast or observation anomalies
#'
#'This function computes anomalies from a multidimensional data array and a 
#'climatology array. 
#'
#'@param data  A numeric array with named dimensions, representing the model or
#'  observational data to be calculated the anomalies. It should involve all
#'  the dimensions in parameter 'clim', and it can have more other dimensions.
#'@param clim A numeric array with named dimensions, representing the 
#'  climatologies to be deducted from parameter 'data'. It can be generated by 
#'  Clim(). The dimensions should all be involved in parameter 'data' with the
#'  same length.
#'@param ncores An integer indicating the number of cores to use for parallel 
#'  computation. The default value is NULL.
#'
#'@return An array with same dimensions as parameter 'data' but with different 
#'  dimension order. The dimensions in parameter 'clim' are ordered first.
#'
#'@examples 
#'# Load sample data as in Load() example:
#'example(Load)
#'clim <- Clim(sampleData$mod, sampleData$obs)
#'ano_exp <- Ano(sampleData$mod, clim$clim_exp)
#'ano_obs <- Ano(sampleData$obs, clim$clim_obs)
#'ano_exp <- Reorder(ano_exp, c(1, 2, 4, 3))
#'ano_obs <- Reorder(ano_obs, c(1, 2, 4, 3))
#'\donttest{
#'PlotAno(ano_exp, ano_obs, startDates, 
#'        toptitle = 'Anomaly', ytitle = c('K', 'K', 'K'), 
#'        legends = 'ERSST', biglab = FALSE, fileout = 'tos_ano.png')
#'}
#'@import multiApply
#'@export
  Ano <- function(data, clim, ncores = NULL) {

  # Check inputs 
  ## data
  if (is.null(data)) {
    stop("Parameter 'data' cannot be NULL.")
  }
  if (!is.numeric(data)) {
    stop("Parameter 'data' must be a numeric array.")
  }
  if (is.null(dim(data))) {  #is vector
    dim(data) <- c(length(data))
    names(dim(data)) <- 'tmp_name'
  }
  if(any(is.null(names(dim(data))))| any(nchar(names(dim(data))) == 0)) {
    stop("Parameter 'data' must have dimension names.")
  }
  ## clim
  if (is.null(clim)) {
    stop("Parameter 'clim' cannot be NULL.")
  }
  if (!is.numeric(clim)) {
    stop("Parameter 'clim' must be a numeric array.")
  }
  if (is.null(dim(clim))) {  #is vector
    dim(clim) <- c(length(clim))
    names(dim(clim)) <- 'tmp_name'
  }
  if (any(is.null(names(dim(clim))))| any(nchar(names(dim(clim))) == 0)) {
    stop("Parameter 'clim' must have dimension names.")
  }
  for (i in 1:length(dim(clim))) {
    if (!(names(dim(clim))[i] %in% names(dim(data)))) {
      stop("Parameter 'data' must have all the dimensions of parameter 'clim'.")
    } else {
      pos <- names(dim(data))[which(names(dim(clim))[i] == names(dim(data)))]
      if ((dim(clim)[i] != dim(data)[pos])) {
        stop("Some dimensions of parameter 'clim' have different length from parameter 'data'.")
      }
    }
  }
  ## 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.")
    }
  }

  ###############################
  # Calculate Ano

  res <- Apply(list(data),
               target_dims = names(dim(clim)),
               output_dims = names(dim(clim)),
               fun = .Ano,
               clim = clim,
               ncores = ncores)$output1

    return(invisible(res))
  }

  .Ano <- function(data, clim) {
    ano <- data - clim
    return(ano)
  }