###
### R function to compute the temporal correlation between the time series of two fields
###

tcorr.f <- function(varname,file1,file2,fileoutcorr,fileoutp) {

library(ncdf4)       # netcdf
library(fields)       # image.plot
library(maptools)     # data_wrld
library(RColorBrewer) # creates nice color schemes

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 corr and p-value

corr <- 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) {
    ok=complete.cases(var2[i, j,],var1[i, j,])
    x=var2[i, j,ok]
    y=var1[i, j,ok]
    n = length(x)
    if (n >= mindat) {
      corr[i, j] = cor(x,y)
      aux = cor.test(x, y)
      sig[i, j] = aux$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))

    # corr
    fillvalue <- 1e32
    dlname <- "corr"
    value <-
      ncvar_def("corr",
                "",
                list(londim, latdim),
                fillvalue,
                dlname,
                prec = "single")
    
    ncfname <- paste0(fileoutcorr,".nc")
    ncout <- nc_create(ncfname, list(value))
    ncvar_put(ncout, value, corr)
    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)


}