xpareto <- function(array3d,p=0.9,upper=TRUE) {

# Compute shape and scale parameters of a Generalized Pareto Distribution for
# a given three-dimensional array with first two space dimensions and third
# time dimension.
#
# Description:
#
#      Returns an array with the shape and scale parameters of a Generalized 
#      Pareto Distribution for the exceedances above a given quantile that
#      defines the threshold. 
#
# Usage:
#
#      xpareto(array3d,p=0.9, upper=TRUE)
#
# Input:
#
#      array3d: three-dimensional array with p longitude points and q latitude
#               points as the first 2 dimensions and n as the third time 
#               dimension 
#
#      p:     quantile to be computed (Default is 0.9). 
#
#      upper: Logic argument. Default is TRUE (examines upper tail of the 
#             distribution). If FALSE lower tail is examined.
#
# Output:
#
#      $out: an array of the computed statistics. First dimension contains
#            the estimated parameters. (e.g. $out[1,,] is the scale parameter
#            and $out[2,,] is the shape parameter). Second and third dimensions
#            are the same space dimensions as of array3d. 
#
#
# Authors:
#
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 16 June 2005
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 
#      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)
#      xpareto(data)
#      xpareto(data,p=0.95)
#      xpareto(data,p=0.1,upper=FALSE)
#


estim.par <- function(y,p) {

if(upper) fpot(y, quantile(y, p,na.rm=T), "gpd",std.err = FALSE)$est
else fpot(-y, quantile(-y, 1-p,na.rm=T), "gpd",std.err = FALSE)$est

}

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

}