xgev <- function(array3d,upper=TRUE,n=365) {

# Compute location, shape and scale parameters of a Generalized Extreme Value
# Distribution for block annual maxima or minima of a given three-dimensional 
# array with first two space dimensions and third time dimension.
#
# Description:
#
#      Returns an array with the location, shape and scale parameters of 
#      a Generalized Extreme Value Distribution for block annual maxima or 
#      minima of a given three-dimensional array with first two space 
#      dimensions and third time dimension.
#
# Usage:
#
#      xgev(array3d,upper=TRUE,n=12)
#
# Input:
#
#      array3d: a 3-dimensional array with p longitude points and q latitude
#               points as the first 2 dimensions and n as the third time 
#               dimension 
#
#      upper: Logic argument. Default is TRUE (examines block annual maxima).
#             If FALSE block annual minima is examined.
#
#      n: Length of the time block. Default is 365 (e.g for daily data 
#         annual maxima/minima)
#

# Output:
#
#      $out: an array of the computed statistics. First dimention contains
#            the estimated parameters. (e.g. $out[1,,] is the location parameter
#            $out[2,,] is the scale parameterand and $out[3,,] is the shape 
#            parameter)
#
#
# Authors:
#
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 14 Dec 2005
#      Christopher 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), 1000)
#      data <- array(rnorm(prod(dim)), dim)
#      xgev(data)
#      xgev(data,n=50)
#      xgev(data,upper=FALSE)
#


d<-dim(array3d)

block.max <- function(x){
apply(matrix(x, n, d[3]/n), 2, max)
}

bmax<-aperm(apply(array3d,c(1,2),block.max),c(2,3,1))

estim.par <- function(y) {

if(upper) fgev(y, std.err = FALSE)$est
else fgev(-y, std.err = FALSE)$est

}

out<-apply(bmax,c(1,2),estim.par)
invisible(list(out=out))

}