applyfield <- function(array3d,fun = mean, ...) {

# Compute basic statistics (e.g. mean, variance, min, max, skewness, quantile)
# for a given three-dimensional array with first two space dimensions (e.g. 
# longitude and latitude) and third time dimension 
#
# Description:
#
#      Returns a matrix with the statistics of interest, specified by the user
#      (see below the list of statistics that can be computed by this funtion),
#      which is computed over the time dimension of a three-dimensional array
#      for all space points of the array.
#
# Usage:
#
#      applyfield(array3d,fun)
#
# Input:
#
#      array3d: a three-dimensional array with p longitude points and q latitude
#               points as the first two dimensions and n as the third time 
#               dimension 
#
#      fun: Name of the statistics to be computed
#           that must be one of the following options: mean, var,
#           min, max, momentskew, ykskew or quantil, 
#           where  momentskew is a moment measure of skewness and
#           ykskew is the Yule-Kendall skewness statistics
#      
#      ...: Additional argument passed to `fun' (e.g. the quantile p to be 
#           computed). The default quantile p to be computed is 0.5 (i.e. 
#           the median).
#
# Output:
#
#           A matrix of the computed statistics
#
# Authors:
#
#      Dag Johan Steinskog <dag.johan.steinskog@nersc.no> 9 June 2005
#      Caio Coelho <c.a.d.s.coelho@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)
#      applyfield(data, mean)
#      applyfield(data, var)
#      applyfield(data, min)
#      applyfield(data, max)
#      applyfield(data, ykskew)
#      applyfield(data, quantil) # compute the median field
#      applyfield(data, quantil,p=0.9) # compute the 90th quantile field

apply(array3d,c(1,2),fun, ...)

}

momentskew <- function(z) {mean(((z-mean(z,na.rm = T))/sd(z,na.rm = T))^3,na.rm = T)}

ykskew <- function(z) {q<-quantile(z,na.rm = T) 
(q[2]-2*q[3]+q[4])/(q[4]-q[2]) }

quantil <- function(w,p=0.5) {quantile(w,p,na.rm = T)}