anomalyfield <- function(data,firstmonth,firstyear,calc="mac") {

# Calculate anomalies by subtracting out mean annual cycle or long-term mean of 
# a three-dimensional p x q x n array of monthly mean values. The structure of this
# array must be as follows: - first dimension p is longitude points
#                           - second dimension q is latitude points
#                           - third dimension n is time (months of the year in 
#                             ascending order, e.g. Jan 1948, Feb 1948, ..., Sep 2005)
#
# Description:
#
#      Returns a p x q x n array containing anomalies
#
# Usage:
#
#      anomalyfield(data,firstmonth,firstyear,calc="mac")
#
# Input:
#
#      data:       p x q x n array (see data structure description above)
#      firstmonth: number of the first month containing data in the 
#                  data array (e.g. 1 (for January)) 
#      firstyear:  first year of data in the data array (e.g. 1948)
#      calc:       string of characters that can be "mac" (Default) to  calculate
#                  anomalies by subtracting the mean annual cycle, or "ltm" to
#                  calculate anomalies by subtracting the long-term mean.
#
# Output:
#
#      $anomaly:   p x q x n array containing anomalies
#
# Authors:
#
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk>  28 October 2005
#      Chris Ferro <c.a.t.ferro@reading.ac.uk>
#
# Examples:
#      
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 100)
#      data <- array(rnorm(prod(dim)), dim)
#      out<-anomalyfield(data,1,1948)
#      out<-anomalyfield(data,5,1970)
#      out<-anomalyfield(data,5,1970,calc="ltm")

if(dim(data)[3] < 12) {stop("Not possible to compute either mean annual cycle or long-term mean. Less than a year of data available")}

if(calc=="mac"){

mean.month <- function(data) {
  apply(matrix(data, nrow = 12), 1, mean)
}
anomaly <- data

if (firstmonth == 1) {
initialmonth <- firstmonth
finalmonth <- floor(dim(data)[3]/12)*12
}

if (firstmonth != 1) {
initialmonth <- 12-firstmonth+2
finalmonth <- (floor((dim(data)[3]-(12-firstmonth+1))/12)*12)+initialmonth-1
}

m <- aperm(apply(data[,,initialmonth:finalmonth],c(1,2),mean.month),c(2,3,1))
anom <- data[,,initialmonth:finalmonth]
for(i in 1:((finalmonth-initialmonth+1)/12)) {
  if(firstmonth == 1)
  anom[,,(((i-1)*12)+1):(12*i)] <- data[,,(((i-1)*12)+1):(12*i)] - m
  else
  anom[,,(((i-1)*12)+1):(12*i)] <- data[,,(((i-1)*12)+initialmonth):((12*i)+initialmonth-1)] - m
}

if (firstmonth == 1 & finalmonth == dim(data)[3]) {
anomaly <- anom
}

if (firstmonth == 1 & finalmonth < dim(data)[3]) {
anomaly[,,firstmonth:finalmonth] <- anom
j <- 1
for(i in (finalmonth+1):dim(data)[3]){
anomaly[,,i] <- data[,,i] - m[,,j]
j <- j+1
}
}

if (firstmonth != 1 & finalmonth == dim(data)[3]) {
j <- 1
for(i in firstmonth:12){
anomaly[,,j] <- data[,,j] - m[,,i]
j <- j+1
}
anomaly[,,initialmonth:finalmonth] <- anom
}

if (firstmonth != 1 &  finalmonth < dim(data)[3]){
j <- 1
for(i in firstmonth:12){
anomaly[,,j] <- data[,,j] - m[,,i]
j <- j+1
}
anomaly[,,initialmonth:finalmonth] <- anom
j <- 1
for(i in (finalmonth+1):dim(data)[3]){
anomaly[,,i] <- data[,,i] - m[,,j]
j <- j+1
}
}


if(firstmonth == 1) y1 <- firstyear
if(firstmonth != 1) y1 <- firstyear+1

y2<- y1+((dim(anom)[3])/12)-1

cat(paste("Mean annual cycle was computed with",as.character(dim(anom)[3]),"months of data","\n",sep=" "))
cat(paste("from Jan",as.character(y1),"to Dec",as.character(y2),"(i.e.",as.character(((dim(anom)[3])/12)),"years)","\n",sep=" "))

}

if(calc=="ltm"){
anomaly <- data
m<-apply(data,c(1,2),mean)
for(i in 1:dim(data)[3]) {
  anomaly[,,i] <- data[,,i] - m
}
cat(paste("Long-term mean was computed with",as.character(dim(data)[3]),"months of data","\n",sep=" "))
}

invisible(list(anomaly=anomaly))

}