#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)

require(ncdf4)
require(reshape2)
require(plyr)
library(dplyr)
library(data.table)
library(lubridate)

print(args)

##############################################################################################
# initialisation variables
#############################################################################################
var <- c("NPP","GPP","IND","DIAMETER_DOM","HEIGHT_DOM","LAI","MOISTRESS","PLANT_STATUS","RDI","CCTRW_001","CCTRW_002","CCTRW_003")
age=c(160,37,37,96,80,80,150,105,47,44,30,135,150,12,30,63,40,20,105,50,50)
monthOrder <- c("janv.","févr.","mars","avril","mai","juin","juil.","août","sept.","oct.","nov.","déc.")
LatOrderSite <- c('RU-Cok','FI-Kaa','FI-Sod','SE-Deg','SE-Fla','IS-Gun','FI-Hyy','RU-Zot','SE-Sk2','SE-SK1','SE-Nor','UK-Gri','DK-Fou','RU-Fyo','SE-Faj','UK-ESa','UK-EBu','UK-AMo','DK-Lva','DK-Ris','DK-Sor','RU-Ha2','RU-Ha1','RU-Ha3','NL-Lut','IE-Ca1','PL-wet','NL-Loo','NL-Hor','NL-Haa','IE-Dri','NL-Ca1','NL-Lan','UK-Her','NL-Mol','UK-PL3','BE-Bra','DE-Meh','UK-Tad','UK-Ham','DE-Geb','DE-Hai','DE-Tha','DE-Gri','DE-Kli','BE-Jal','BE-Lon','DE-Wet','BE-Vie','DE-Bay','CZ-BK1','SK-Tat','CZ-wet','FR-Gri','FR-Hes','FR-Fon','HU-Mat','CH-Lae','CH-Oe2','CH-Oe1','AT-Neu','HU-Bug','IT-Ren','IT-Mal','IT-MBo','IT-Lav','FR-Lq1','FR-Lq2','IT-LMa','IT-PT1','IT-Cas','FR-LBr','IT-Non','FR-Pue','IT-SRo','IT-Lec','IT-Pia','IT-Ro1','IT-Ro2','ES-VDA','IT-Amp','IT-Col','IT-Cpz','IT-Cp2','IT-BCi','ES-LMa','ES-ES1','ES-ES2','PT-Esp','PT-Mi1','PT-Mi2')

#########################################################################################
#   extract result from netcdf data
##########################################################################################
Table.Summary <- NULL
h=0
for (l in var) {
    h <- h+1
    j=0
    Table.final <- NULL
    for (i in LatOrderSite ) {
        j <- j+1
        file_path <- paste0(args[1],i,'_ALL.nc')
        if (file.exists(file_path)==FALSE) next
        else print(paste('OK pour :',i))
        grid.nc <- nc_open(file_path)    
        lon <- ncvar_get(grid.nc,"lon")
        lat <- ncvar_get(grid.nc,"lat")
        time <- ncvar_get(grid.nc,"time_counter")
        if (args[5]=='D')tunits <- seq(as.Date("2000/1/1"), by = "day", length.out = length(time))
        else tunits <- seq(as.Date("2000/1/1"), by = "month", length.out = length(time))
        PFT  <- ncvar_get(grid.nc,"PFT")
        tmp_array <- t(ncvar_get(grid.nc,l))    
        tmp_array <- melt(tmp_array, value.name = l)
        tmp_array[,"Date"] <- tunits[tmp_array[,"Var1"]]
        tmp_array[,"Var1"] <- NULL
        tmp_array[,"Lat"] <- round(lat,4)
        tmp_array[,"Lon"] <- lon
        tmp_array[,"Site"] <- substr(i,1,6)
        if(j==1) { Table.final <- as.data.frame(tmp_array) }
        else { Table.final <- bind_rows(Table.final,as.data.frame(tmp_array)) }    
        nc_close(grid.nc)   
    }
    if(h==1) { Table.Summary <- Table.final }
    else { Table.Summary <- left_join(Table.Summary,Table.final) }
}
if(args[5]=='D') {
   Table.Summary$Day <- factor(format(Table.Summary$Date, "%d"))
   Table.Summary$DOY <- factor(strftime(Table.Summary$Date, format = "%j"))  
}
Table.Summary$Month <- factor(format(Table.Summary$Date, "%b"),levels=monthOrder)
Table.Summary$Year <- as.numeric(format(Table.Summary$Date , "%Y"))
print(summary(Table.Summary))
fwrite(Table.Summary,paste0(args[2],'Table.Summary_',args[7],'.txt'), sep=' ',quote=FALSE)
############################################################################################

monthOrder <- c("janv.","févr.","mars","avril","mai","juin","juil.","août","sept.","oct.","nov.","déc.")
reloc <- time_length(interval(start = ymd("1970-01-01"),end = ymd("2000-01-01") ),unit = "seconds")
climfile2 <- c("FR-Hes_1997-2006.nc","FI-Sod_2000-2006.nc","IT-Non_2001-2006.nc","DE-Bay_1996-1999.nc","DE-Wet_2002-2006.nc","SE-Nor_1996-2005.nc","SE-Fla_1996-2002.nc")
var.f=c('sum','mean')
var.o=c("GPPalt")
var.oo=c("GPP")
var.conv=c(1,1,1800)
###########################################################################################
# extract climate and GPP obs from netcdf data
############################################################################################
climate.final <- NULL
for(l in c("FR","IT","DE","DK","FI","SE","NL","BE","CH","RU","CZ","OLD")){
    if(l=='OLD') climfile=climfile2
    else climfile=list.files('/home/surface5/nraoult/FLUXNET2015/done/',pattern = l)
    for(i in climfile){
    k=0;s=1
    Site <- substr(i,1,6)
    clim.t=data.frame(DOY=1:366)
    if(l=='OLD') obs.nc <- nc_open(paste0('/home/orchideeshare/igcmg/IGCM/BC/OOL/OL2/FLUXNET/',i))
    else obs.nc <- nc_open(paste0('/home/surface5/nraoult/FLUXNET2015/done/',i))
    if(l=='OLD') var.final=var.oo
    else var.final=var.o
    for (j in var.final){
        k=k+1
        func=var.f[k]
        conv=var.conv[k]
        var.nc <- ncvar_get(obs.nc, j)
        var.nc <- ifelse(var.nc>10e30,NA,var.nc)
        lat <- ncvar_get(obs.nc,"lat")
        lon <- ncvar_get(obs.nc,"lon")
        time <- ncvar_get(obs.nc, "tstep")+ reloc
        class(time) = c('POSIXt','POSIXct')
        df <- data.frame(date=time,obs=var.nc*conv)
        df$Month <- factor(format(df$date, "%b"),levels = monthOrder)
        df$Year <- as.numeric(format(df$date , "%Y"))
        df$DOY <- factor(strftime(df$date, format = "%j"))
        ob <- tapply(df[,'obs'],list(df[,'DOY'],df[,'Year']),func)
        clim.t <- data.frame(clim.t,apply(ob,1,max,na.rm=TRUE),
            apply(ob,1,min,na.rm=TRUE),apply(ob,1,median,na.rm=TRUE))
        names(clim.t)[c(s+1,s+2,s+3)]=c(paste0('GPP.max'),paste0('GPP.min'),paste0('GPP.med'))
        s=s+3
    }
    clim.t$Site <- Site
    clim.t$Lat <- round(lat,4)
    clim.t$Lon <- round(lon,4)
    climate.final <- rbind(climate.final,clim.t)
    nc_close(obs.nc)     
    }
}
#fwrite(climate.final,'climate_final.txt', sep=' ',quote=FALSE, na.string=NA)
##########################################################################################