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

require(reshape2)
require(ggplot2)
require(grid)
require(gridExtra)
require(plyr)
require(dplyr)
require(data.table)

##############################################################################################
# initialisation variables
#############################################################################################
bo=c(0,0,0,17,18,6,19,8,9,10,20,21,13,14,15,16)
names=c('','','','','P.Sylvestris','P.Pinaster','PiceaSp','','','QuercusIlex','','BetulaSp','FagusSyl.','QuercusRob.','PopulusSp','')
leaftype=c('','','','','','','','','','','','Deciduous','Deciduous','Deciduous','Deciduous','')
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-Obe','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','CH-Dav','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')
site.sp=list(
    P.Sylvestris=c('NL-Loo','FI-Hyy','FI-Sod'),
    P.Pinaster=c('IT-SRo','FR-LBr'),
    PiceaSp=c('DE-Bay','DE-Tha ','DE-Wet','IT-Ren','RU-Fyo','SE-Nor','SE-Fla'),
    QuercusIlex=c('FR-Pue','IT-Cpz','PT-Mi1'),
    BetulaSp=c(''),
    FagusSyl.=c('FR-Hes','DE-Hai','IT-Col','DK-Sor'),
    QuercusRob.=c('IT-Non','FR-Fon','IT-Ro1','IT-Ro2'),
    PopulusSp=c('IT-PT1')
)

en.mo <- c(
    janv.="January",
    févr.="Febuary",
    mars ="March",
    avril="April",
    mai="May",
    juin="June",
    juil.="July",
    août="August",
    sept.="September",
    oct.="October",
    nov.="November",
    déc.="December")

vol.thres <- list(
    PiceaSp=list(
        YC1=list(DE=c(644,886,26,34,645,944,30,36)),
        YC2=list(DE=c(613,842,22,30,785,958,29,31)),
        YC3=list(DE=c(579,795,25,32,689,1021,28,34)),
        YC4=list(DE=c(543,747,24,31,715,1065,27,32)),
        YC5=list(DE=c(505,698,23,30,746,1117,26,31)),
        YC6=list(DE=c(467,648,22,29,780,1173,25,30)),
        YC7=list(DE=c(427,598,21,28,821,1238,24,29)),
        YC8=list(DE=c(387,547,20,26,872,1320,22,27)),
        YC9=list(DE=c(346,496,19,25,934,1418,21,26)),       
        YC10=list(DE=c(304,445,18,24,1011,1541,20,24)),
        YC11=list(DE=c(262,394,16,22,1115,1701,18,23)),
        YC12=list(DE=c(220,343,14,20,1266,1939,20,21)),
        YC13=list(DE=c(177,291,12,18,1496,2309,15,19)),
        YC14=list(DE=c(131,238,10,15,2361,3911,15,19))       
    )
)

############################################################################################
climate.final <- fread('climate.final.txt')
climate.final2 <- fread('climate_final.txt', fill=TRUE, na.string=c("Inf","-Inf","NA"))
print(summary(climate.final2 ))

Info.Site <- fread('site_info.txt')

clim.sum <- NULL; clim.sum2 <-NULL; LAI.sum <- NULL; PS.sum <- NULL; WSM.sum <- NULL
LAI.doy <- NULL; RDI.sum=NULL; GPP.doy=NULL; tot.sum=NULL

for(i in args[-c(1:6)]){

    Table.Summary <- fread(paste0(args[2],'Table.Summary_',i,'.txt')) 
############################################################################################

new=NULL
for (j in LatOrderSite ) {
        if (file.exists(paste0(args[1],j,'_ALL.nc'))==FALSE) next
        else new <- c(new,j)
}
LatOrderSite <- new
print(LatOrderSite)
###########################################################################################
#                    REshape data to montly value only
############################################################################################

if(args[5]=="M"){

attach(Table.Summary)
ob <-melt(tapply(GPP,list(Month,Var2,Site),quantile,probs=0.75), value.name ='GPP75')
ob2 <-melt(tapply(GPP,list(Month,Var2,Site),quantile,probs=0.25), value.name ='GPP25')
obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
names(obf) <- c('Month','PFT','Site','GPP75','GPP25')
obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
clim <- merge(obf,climate.final,c('Month','Site'))
clim$Month <- factor(clim$Month,levels=monthOrder)
clim$Version <- i
detach(Table.Summary)
clim.sum <- rbind(clim.sum,clim)

Table.Sub <- subset(Table.Summary, Month %in% c("juin","juil.","août"))
attach(Table.Sub)
ob <-melt(tapply(NPP/GPP,list(Var2,Site),quantile,probs=0.75,na.rm=TRUE), value.name ='Ratio75')
ob2 <-melt(tapply(NPP/GPP,list(Var2,Site),quantile,probs=0.25,na.rm=TRUE), value.name ='Ratio25')
obf <- merge(ob,ob2,c('Var1','Var2'))
names(obf) <- c('PFT','Site','Ratio75','Ratio25')
obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
clim2 <- merge(obf,climate.final,c('Site'))
clim2 <- subset(clim2, Month=='mai')
clim2$Version <- i
detach(Table.Sub)
clim.sum2 <- rbind(clim.sum2,clim2)

attach(Table.Summary)
ob <-melt(tapply(LAI,list(Month,Var2,Site),max,na.rm=TRUE), value.name ='LAI75')
ob2 <-melt(tapply(LAI,list(Month,Var2,Site),min,na.rm=TRUE), value.name ='LAI25')
obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
names(obf) <- c('Month','PFT','Site','LAI75','LAI25')
obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
LAII <- merge(obf,climate.final,c('Month','Site'))
LAII$Month <- factor(LAII$Month,levels=monthOrder)
LAII$Version <- i
detach(Table.Summary)
LAI.sum <- rbind(LAI.sum,LAII)

}

if(args[5]=="D"){

Table.Sub <- Table.Summary
attach(Table.Sub)
ob <-melt(tapply(LAI,list(Var2,Site,DOY),max,na.rm=TRUE), value.name ='LAIMax')
ob2 <-melt(tapply(LAI,list(Var2,Site,DOY),min,na.rm=TRUE), value.name ='LAIMin')
obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
names(obf) <- c('PFT','Site','DOY','LAIMax','LAIMin')
LAII <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
LAII <- merge(obf,climate.final2,c('Site','DOY'))
LAII <- merge(LAII,Info.Site,'Site')
LAII$Version <- i
detach(Table.Sub)
LAI.doy <- rbind(LAI.doy,LAII)


Table.Sub <- Table.Summary
attach(Table.Sub)
ob <-melt(tapply(GPP,list(Var2,Site,DOY),max,na.rm=TRUE), value.name ='GPPMax')
ob2 <-melt(tapply(GPP,list(Var2,Site,DOY),min,na.rm=TRUE), value.name ='GPPMin')
ob3 <-melt(tapply(GPP,list(Var2,Site,DOY),median,na.rm=TRUE), value.name ='GPPMed')
obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
obf <- merge(obf,ob3,c('Var1','Var2','Var3'))
names(obf) <- c('PFT','Site','DOY','GPPMax','GPPMin','GPPMed')
LAII <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
LAII <- merge(obf,climate.final2,c('Site','DOY'))
LAII <- merge(LAII,Info.Site,c('Site'))
LAII$Version <- i
detach(Table.Sub)
GPP.doy <- rbind(GPP.doy,LAII)


TS <- Table.Summary[,c('Year','DOY','PLANT_STATUS','Var2','Site')]
TS <- subset(TS, (PLANT_STATUS==5 & DOY>150) | PLANT_STATUS==3 )
TS1 <- TS %>% group_by(Var2,Site,Year) %>% slice(which(PLANT_STATUS==3)[1])
TS1[,'Phenology'] <- 'budburst'
TS2 <- TS %>% group_by(Var2,Site,Year) %>% slice(which(PLANT_STATUS==5)[1])
TS2[,'Phenology'] <- 'leaf senescence'
TS <- rbind(TS1,TS2)


attach(TS)
ob <-melt(tapply(DOY,list(Phenology,Var2,Site),max,na.rm=TRUE), value.name ='PSmax')
ob2 <-melt(tapply(DOY,list(Phenology,Var2,Site),min,na.rm=TRUE), value.name ='PSmin')
obf <- merge(ob,ob2,c('Var1','Var2','Var3'))
names(obf) <- c('Phenology','PFT','Site','PSmax','PSmin')
obf <- subset(obf, PFT %in% c(5,6,7,10,12,13,14,15,16,17,18,19,20,21))
PS <- merge(obf,climate.final,c('Site'))
PS <- unique(PS[,c('Phenology','PFT','Site','PSmax','PSmin','Lat')])
PS[,'Site'] <- factor(PS[,'Site'],levels=rev(LatOrderSite))
PS$Version <- i
detach(TS)
PS.sum <- rbind(PS.sum,PS)

tot.tab = merge(Table.Summary,Info.Site,'Site')
tot.tab$Version <- i
tot.sum <- rbind(tot.sum,tot.tab)

}

}
print('Reformating done !')

############################################################################################
# forestry graphics 
############################################################################################
if(args[5]=="M"){

pdf(paste0('Forestry_Eval_',args[6],'.pdf'),width=16,height=10)
for(i in 7){
    STR.bot <- NULL
    fm='DE'
    ycorder <- c('YC1','YC2','YC3','YC4','YC5','YC6','YC7','YC8','YC9','YC10','YC11','YC12','YC13','YC14')
    for(j in ycorder){
        vol.bot.min <- vol.thres[[names[[i]]]][[j]][[fm]][1]/10000
        vol.bot.max <- vol.thres[[names[[i]]]][[j]][[fm]][2]/10000   
        STR.b <- subset(Table.Summary, Var2==i & WOOD_VOL > vol.bot.min & WOOD_VOL < vol.bot.max & Site %in% c('DE-Hai','DE-Tha','DE-Wet','DE-Bay'))
        if(nrow(STR.b) == 0) next
        STR.b[,"Yield"] <- j
        STR.b[,"ind.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][5]
        STR.b[,"ind.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][6]
        STR.b[,"dg.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][3]
        STR.b[,"dg.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][4]
        STR.b[,"H.min"] <- vol.thres[[names[[i]]]][[j]][[fm]][7]
        STR.b[,"H.max"] <- vol.thres[[names[[i]]]][[j]][[fm]][8]
        STR.bot <- rbind(STR.bot,STR.b)
    } 
    STR.bot[,"Yield"] <- factor(STR.bot[,"Yield"],levels=ycorder)
    x2= ggplot(STR.bot,aes(as.factor(Site), IND*10000, linetype=Version))+
        geom_boxplot(fill="grey50")+ ylim(0,14000)+
        geom_hline(aes(yintercept = ind.max),linetype = 2, colour='red',size=1.5)+
        geom_hline(aes(yintercept = ind.min),linetype = 2, colour='red',size=1.5)+
        xlab("") + ylab("number of individue (tree/ha)") + coord_flip()+ facet_grid(.~Yield)
    x3= ggplot(STR.bot,aes(as.factor(Site), DIAMETER_DOM*100, linetype=Version))+
        geom_boxplot(fill="grey50")+ ylim(0,40) +
        geom_hline(aes(yintercept = dg.max),linetype = 2, colour='red',size=1.5)+
        geom_hline(aes(yintercept = dg.min),linetype = 2, colour='red',size=1.5)+
        xlab("") + ylab("dominant diameter (cm)") + coord_flip()+ facet_grid(.~Yield)
    x4= ggplot(STR.bot,aes(as.factor(Site), HEIGHT_DOM, linetype=Version))+
        geom_boxplot(fill="grey50")+ ylim(0,25) +
        geom_hline(aes(yintercept = H.max),linetype = 2, colour='red',size=1.5)+
        geom_hline(aes(yintercept = H.min),linetype = 2, colour='red',size=1.5)+
        xlab("") + ylab("dominant height (m)") + coord_flip()+ facet_grid(.~Yield)   
    grid.arrange(x2,x3,x4,nrow=3)
}
dev.off()

##############################################################################################
# LAI WSM graphics V2
##############################################################################################
pdf(paste0('LAI-WSM_eval_',args[7],'.pdf'),width=16,height=10)
for(i in c(5,7,12,13,14)){
    subtable <- subset(LAI.sum, PFT==i  & Lat<57 & Lat>=45)
    subtablebo <- subset(LAI.sum, PFT==bo[i] & Lat>=57 )
    subtableso <- subset(LAI.sum, PFT==i  & Lat<45 )
    subtableBE <- subset(subtablebo, PFT==bo[i] & Site %in% site.sp[[names[i]]])
    subtableTE <- subset(subtable, PFT==i &  Site %in% site.sp[[names[i]]])
    subtableME <- subset(subtableso, PFT==i & Site %in% site.sp[[names[i]]])
    x1= ggplot(LAI.sum,aes(as.factor(Lat),y=LAI25, ymin=LAI25,ymax=LAI75,linetype=Version)) +
        theme_bw()+ ylim(0,10)+
        geom_crossbar(data= subtablebo, fill="deepskyblue3",alpha=0.3,fatten=1)+
        geom_crossbar(data= subtable, fill="chartreuse4",alpha=0.3,fatten=1)+
        geom_crossbar(data= subtableso, fill="darkorange3",alpha=0.3,fatten=1) +
        geom_crossbar(data= subtableTE, fill="chartreuse4",fatten=1) +
        geom_crossbar(data= subtableME, fill="darkorange3",fatten=1) +
        geom_crossbar(data= subtableBE, fill="deepskyblue3",fatten=1)+
        xlab("") + ylab("LAI in m2/m2") + scale_x_discrete(labels = rev(LatOrderSite)) +
        facet_grid(.~Month,labeller = labeller(Month=en.mo)) + coord_flip()+
        theme(panel.spacing = unit(0, "lines"),strip.text.x = element_text(face = "bold",size=18))
    subtable <- subset(WSM.sum, PFT==i  & Lat<57 & Lat>=45)
    subtablebo <- subset(WSM.sum, PFT==bo[i] & Lat>=57 )
    subtableso <- subset(WSM.sum, PFT==i  & Lat<45 )
    subtableBE <- subset(subtablebo, PFT==bo[i] & Site %in% site.sp[[names[i]]])
    subtableTE <- subset(subtable, PFT==i &  Site %in% site.sp[[names[i]]])
    subtableME <- subset(subtableso, PFT==i & Site %in% site.sp[[names[i]]])
    x2= ggplot(WSM.sum,aes(as.factor(Lat),y=WSM25, ymin=WSM25,ymax=WSM75,linetype=Version)) +theme_bw()+
        geom_crossbar(data= subtablebo, fill="deepskyblue3",alpha=0.3,fatten=1)+
        geom_crossbar(data= subtable, fill="chartreuse4",alpha=0.3,fatten=1)+
        geom_crossbar(data= subtableso, fill="darkorange3",alpha=0.3,fatten=1) +
        geom_crossbar(data= subtableTE, fill="chartreuse4",fatten=1) +
        geom_crossbar(data= subtableME, fill="darkorange3",fatten=1) +
        geom_crossbar(data= subtableBE, fill="deepskyblue3",fatten=1)+
        xlab("") + ylab("monthly average water stress") + ylim(0,1) + coord_flip() +
        scale_x_discrete(labels = rev(LatOrderSite)) + 
        facet_grid(.~Month,labeller = labeller(Month=en.mo)) 
    grid.arrange(x1,x2,nrow=2,top = textGrob(paste('PFT:',i,'called',names[i], "10% and 90% quantile of a 250 years simulation"),gp=gpar(fontsize=20,font=3)))
}
dev.off()
###################################################################################################
}

if(args[5]=="D" & args[6]=="SBG"){

######################################################################
# Global analysis                                                    #
######################################################################

pdf(paste0('global_eval_',args[7],'.pdf'),width=16,height=10)
for(i in c(10,13,14,5,6,7,18,19)){
datasp=subset(GPP.doy,PFT==i & PFT==SPE,Site)
droplevels(datasp)
nbsite = unique(datasp)
nofacet = dim(nbsite)[1]>1
x1 = ggplot(subset(GPP.doy, PFT==SPE & PFT==i), aes(x=DOY,ymin=GPPMin,ymax=GPPMax,colour=Version)) +
     geom_ribbon(alpha=0.4,size=0.2)+ xlab("DOY") + ylab("GPP") + 
     {if(nofacet)facet_wrap(.~Site, scales = "free")}+ 
     geom_ribbon(aes(ymin=GPP.min,ymax=GPP.max),alpha=0.1,colour='black',size=0.2)+
     geom_smooth(aes(y=GPP.med),colour='black',linetype=2,size=0.2,span = 0.3)+
     geom_smooth(aes(y=GPPMed,colour=Version),linetype=2,size=0.2,span = 0.3)
grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i],"Min and Max of 50 years simulation"),gp=gpar(fontsize=20,font=3)))

x1 = ggplot(subset(LAI.doy, PFT==SPE & PFT==i), aes(x=DOY,ymin=LAIMin,ymax=LAIMax,colour=Version))+
     geom_ribbon(alpha=0.4,size=0.2)+ xlab("DOY") + ylab("LAI") + 
     geom_line(aes(y=LAImax),linetype=2,colour='black',size=0.2)+
     geom_line(aes(y=LAImin),linetype=2,colour='black',size=0.2)+
     {if(nofacet)facet_wrap(.~Site, scales = "free")}
grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i],"Interannual Min and Max"),gp=gpar(fontsize=20,font=3)))

x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=LAI, colour=Version))+
     ylab("LAI") + geom_smooth(aes(linetype=Version),size=0.2,span = 0.1)+
     {if(nofacet)facet_wrap(.~Site, scales = "free")}+ ylim(0,10)
grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))

x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=GPP,colour=Version))+
     ylab("GPP") + geom_smooth(aes(linetype=Version),size=0.2,span = 0.1)+
     {if(nofacet)facet_wrap(.~Site, scales = "free")}
grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))

x1 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=IND,colour=Version))+
     ylab("IND") + geom_line(aes(linetype=Version),size=0.2)+ ylim(0,1.4)+
     {if(nofacet)facet_wrap(.~Site, scales = "free")}
x2 = ggplot(subset(tot.sum, Var2==SPE & Var2==i), aes(x=as.Date(Date),y=RDI,colour=Version))+
     ylab("RDI") + geom_line(aes(linetype=Version),size=0.2)+ ylim(0.25,1)+
     {if(nofacet)facet_wrap(.~Site, scales = "free")}
grid.arrange(x1,x2,nrow=2,top = textGrob(paste('PFT:',i,'called',names[i]),gp=gpar(fontsize=20,font=3)))

#If(leaftype(i)=='Deciduous'){
#   x1 = ggplot(subset(PS.sum,PFT==i),aes(Site, y=PSmin, ymin=PSmin,ymax=PSmax,colour=Version)) +
#                  geom_crossbar(aes(group=Phenology),fill='black',alpha=0.3,fatten=1)+
#                  geom_hline(yintercept = 90,linetype = 2)+ 
#                  geom_hline(yintercept = 130,linetype = 2)+
#                  geom_hline(yintercept = 310,linetype = 2)+
#                  geom_hline(yintercept = 285,linetype = 2)+
#                  xlab("") + ylab("DOY") + ylim(0,366) + coord_flip()
#   grid.arrange(x1,nrow=1,top = textGrob(paste('PFT:',i,'called',names[i], "Min and Max of 50 years simulation"),gp=gpar(fontsize=20,font=3)))
#}

}
dev.off()

}

if(args[5]=="D" & args[6]=="SFR"){

}