#setwd("/Users/macos/Documents/work/2023/ThayTan")
library(ncdf4); library(ggplot2)
library(reshape2); library(ggpubr)
clead <- c("01","02","03","04","05","06")
model <- "RegCM"
#idir <- "/home/hapt/nhap/thay/"
idir <- "/home/tanpv/Data/562_VNU/0.25_VN/Dat_Reg/"
odir <- "/home/tanpv/Progs/562_VNU/PDF/" 
for (lead in c(clead)){
    ifile <- paste("pr_",model,"_lead_",lead,sep="")
    iname <- paste(idir,ifile,".nc",sep="")
    nc_file <- nc_open(iname)
    cvar.name <- c("obs","mod","cli","ann")    # Order lai, add more or remove
    for (ivar in c(cvar.name)){
         tmp.cmd <- paste(ivar,".dat<-ncvar_get(nc_file,'",ivar,"')",sep="")
         eval(parse(text=tmp.cmd))
    }
    
    nreg <- 7
    nmon <- 12
    itmp <- 0
    for (ivar in c(cvar.name)){
          itmp <- itmp + 1
          tmp.cmd <- paste("tmp.dat<-",ivar,".dat",sep="")
          eval(parse(text=tmp.cmd))
          nsample <- dim(tmp.dat)[4]
          my.dat <- tmp.dat
          dimnames(my.dat) <- list(c(seq(7)),
                                   c(paste("m",seq(7),sep="")),
                                   c(seq(12)),
                                   c(seq(nsample)))
          tmp.mel <- melt(my.dat)
          colnames(tmp.mel) <- c("reg","mem","mon","sample","value")
          tmp.mel$dat <- ivar
          if (itmp == 1) { 
              odat <- tmp.mel
          } else {
              odat <- rbind(odat,tmp.mel)
          }
          remove(tmp.mel,sub.dat,my.dat)
    }
    
    for (ireg in seq(nreg)){
         draw.dat <- subset(odat,odat$reg==ireg)
         a <- data.frame(draw.dat$value,draw.dat$mon,draw.dat$dat)
         colnames(a) <- c("value","mon","dat")
         a$mon <- factor(a$mon, levels = seq(12))
         a$dat <- as.character(a$dat)
         remove(p)
         a$dat <- factor(a$dat, levels=c(cvar.name))
         p <- ggplot(a,aes(x=mon,y=value))+
           geom_boxplot(aes(fill=dat),
                   outlier.size=0.1,outlier.shape = 1,outlier.color ="gray50" )+
           scale_y_continuous(breaks = c(seq(0,1000,200)),limits = c(0,1000))+
           scale_fill_manual(name="",values = c("obs"="#FF0000","mod"="#ffff00","cli"="#0000FF","ann"="#009900"))+ 
           ggtitle(paste("R",ireg,"-",model,"-LT=",lead,sep="")) + xlab("")+ylab("(mm)")+labs(x=NULL)
         tmp.cmd <- paste("p",ireg,"<-p",sep="")
         eval(parse(text=tmp.cmd))+theme(legend.position = "none")
    }
    
    save.image <- ggarrange(p1,p2,p3,p4,p5,p6,p7,ncol = 1, nrow = 7,
                            common.legend = TRUE, legend="bottom")
    save.image

    png_file <- paste(odir,ifile,".png",sep="")
    png(png_file, width = 8.0, height = 14.0, units = 'in', res = 300)

#    pdf(paste(odir,ifile,".pdf",sep=""), width = 4.5, height = 14.0)
#    plot(save.image) 
    print(save.image)
    dev.off()
    remove(save.image)
}