###
### R function to compute the difference between the means of two fields
###

diff.f <- function(varname,file1,file2,fileoutdiff,fileoutp) {

library(ncdf4)       
library(fields)      
library(maptools)    
library(RColorBrewer) 

file2.nc = nc_open(file2)
var2 = ncvar_get(file2.nc,varname)
var2[var2 == "-8.99999987309029e+33"] <- NA
file2.nc$dim$x$vals -> lon
file2.nc$dim$y$vals -> lat
ni = dim(var2)[1]
nj = dim(var2)[2]
nt = dim(var2)[3]

mindat=0.75*nt ## minimum number of data in the series to be considered
  
file1.nc = nc_open(file1)
var1 = ncvar_get(file1.nc,varname)
var1[var1 == "-8.99999987309029e+33"] <- NA

## compute diff and p-value

diff <- matrix(data = NA,
               nrow = ni,
               ncol = nj)
sig <- matrix(data = NA,
              nrow = ni,
              ncol = nj)

for (i in 1:ni) {
  for (j in 1:nj) {
    n1 = length(which(!is.na(var1[i, j,])))
    n2 = length(which(!is.na(var2[i, j,])))
    if (n1 >= mindat & n2 >= mindat) {
      diff[i, j] = mean(var2[i, j,], na.rm = TRUE) - mean(var1[i, j,], na.rm = TRUE)
      dum = t.test(var2[i, j,], var1[i, j,], na.omit = TRUE)
      sig[i, j] = dum$p.value
    }
  }
}
 
## write netcdf output files

  # define dimensions
  londim <- ncdim_def("lon", "degrees_east", as.double(lon))
  latdim <- ncdim_def("lat", "degrees_north", as.double(lat))

    # diff
    fillvalue <- 1e32
    dlname <- "diff"
    value <-
      ncvar_def("diff",
                "",
                list(londim, latdim),
                fillvalue,
                dlname,
                prec = "single")
    
    ncfname <- paste0(fileoutdiff,".nc")
    ncout <- nc_create(ncfname, list(value))
    ncvar_put(ncout, value, diff)
    ncatt_put(ncout, "lon", "axis", "X") 
    ncatt_put(ncout, "lat", "axis", "Y")
    nc_close(ncout)

    # p-value
    fillvalue <- 1e32
    dlname <- "pval"
    value <-
      ncvar_def("pval",
                "",
                list(londim, latdim),
                fillvalue,
                dlname,
                prec = "single")
    
    ncfname <- paste0(fileoutp,".nc")
    ncout <- nc_create(ncfname, list(value))
    ncvar_put(ncout, value, sig)
    ncatt_put(ncout, "lon", "axis", "X") 
    ncatt_put(ncout, "lat", "axis", "Y")
    nc_close(ncout)


}