library(ncdf4); library(ggplot2)
library(reshape2); library(ggpubr)
#clead <- c("01","02","03","04","05","06")
clead <- c("04")
model <- "RegCM"
idir <- "/home/tanpv/Data/562_VNU/0.25_VN/Dat_Reg/"
odir <- "/home/tanpv/Progs/562_VNU/PDF/" 
#c.mem <- c("m1","m2","m3","m4","m5","m6","m7")
c.mem <- c("m4")
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)
    }
    
    odat.1 <- odat
    for (imem in c(c.mem)){
         odat <- subset(odat.1,odat.1$mem==imem)
            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,"-",imem,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,"-",imem,".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)
    }
}