applyfieldnew <- function(lon,lat,array3d,fun = mean, ...,nonmissing=0.5) {

# 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. 
#
# Note: This function allows the user to specify through the parameter
#       nonmissing (defaul is 0.5) the acceptable percentage of nonmissing
#       values in the time series of each grid point of the three-dimensional
#       array. This is the main difference between this function and the
#       function applyfield.r
#
# 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:
#
#      applyfieldnew(lon,lat,array3d,fun,nonmissing)
#
# Input:
#
#       lon: vector with p longitude values
#
#       lat: vector with q latitude values 
# 
#   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 
#
#nonmissing: Only grid points with fraction given by 'nonmissing' 
#            (between 0 and 1) of non-missing values are used to compute 
#            the statistics below. Default is 0.5.
#
#      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)
#      applyfieldnew(x,y,data, mean)
#      applyfieldnew(x,y,data, var)
#      applyfieldnew(x,y,data, min)
#      applyfieldnew(x,y,data, max)
#      applyfieldnew(x,y,data, ykskew)
#      applyfieldnew(x,y,data, quantil) # compute the median field
#      applyfieldnew(x,y,data, quantil,p=0.9) # compute the 90th quantile field



# reshape three-dimensional array into a data matrix with
# time as first dimension and space as sencond dimension
data <- reshapefield(lon,lat,array3d)$out

# compute percentage of non-missing values at each grid point
aux <- apply(data,2,function(x)sum(!is.na(x))/(length(x)))

# identify grid points with less than 50% missing values
index <- (1:length(aux))[aux >= nonmissing]

out<-apply(data[,index],2,fun, ...)

out1 <- rep(NA,dim(data)[2])
out1[index]<-out

reshapefield(lon,lat,t(as.matrix(out1)))$out[,,1,drop=T]

}

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)}