#setwd('C:/Users/quant/Desktop/temp/BCSD/precipitation/imodel')
#setwd('/mnt/c/Users/quant/Desktop/temp/BCSD/precipitation/imodel')

setwd('/work/users/hmo/remoclic/Quan/BC/precipitation/imodel')

rm(list = ls())
library(ncdf4)
library(qmap)

for (m in 1:12) {
  if (m<10) {
    obs <- nc_open(print(paste("APH_0",m,".nc",sep="")))
    his <- nc_open(print(paste("imodel.his_0",m,".nc",sep="")))
    fu45 <- nc_open(print(paste("imodel.fu85_0",m,".nc",sep="")))
    fu85 <- nc_open(print(paste("imodel.fu45_0",m,".nc",sep="")))
  } else {
    obs <- nc_open(print(paste("APH_",m,".nc",sep="")))
    his <- nc_open(print(paste("imodel.his_",m,".nc",sep="")))
    fu45 <- nc_open(print(paste("imodel.fu85_",m,".nc",sep="")))
    fu85 <- nc_open(print(paste("imodel.fu45_",m,".nc",sep="")))
  }
  obs.pr <- ncvar_get(obs, "pr")
  his.pr <- ncvar_get(his, "pr")
  fu45.pr <- ncvar_get(fu45, "pr")
  fu85.pr <- ncvar_get(fu85, "pr")
  
  qm <- function(obscell,hiscell,fucell){
    QM=fitQmap(obscell,hiscell,method=c("QUANT"))
    bias.corrected <- doQmap(fucell,QM)
    return(bias.corrected)
  }
  
  host45 <- array(data=NA,dim=dim(fu45.pr))
  host85 <- array(data=NA,dim=dim(fu85.pr))
  
  for (i in 1:dim(obs.pr)[1]) {
    for (j in 1:dim(obs.pr)[2]) {
      if (is.na(obs.pr[i,j,1])) {
        host45[i,j,] <- NA
        host85[i,j,] <- NA
      } else {
        host45[i,j,] <- qm(obs.pr[i,j,],his.pr[i,j,],fu45.pr[i,j,])
        host85[i,j,] <- qm(obs.pr[i,j,],his.pr[i,j,],fu85.pr[i,j,])
        
      }
    }
  }
  
  #### WRITE NETCDF FILE OUTPUT ####
  #Define dimension of output
  print(paste("Writing the precipitation bias correction using QM"))
  lon.vals <- ncvar_get(fu45, "lon")
  lat.vals <- ncvar_get(fu45, "lat")
  time.vals <- ncvar_get(fu45, "time")
  nx <- length(lon.vals)
  ny <- length(lat.vals)
  nt <- length(time.vals)
  lon1 <- ncdim_def("longitude", "degrees_east", lon.vals, unlim=FALSE, create_dimvar=TRUE)
  lat2 <- ncdim_def("latitude", "degrees_north", lat.vals, unlim=FALSE, create_dimvar=TRUE)
  if (m<10) {
    time <- ncdim_def("Time",paste("years since 2006-0",m,"-16 12:00:00",sep=""), 0:93, unlim=TRUE,create_dimvar=TRUE)
  } else {
    time <- ncdim_def("Time",paste("years since 2006-",m,"-16 12:00:00",sep=""), 0:93, unlim=TRUE,create_dimvar=TRUE)
  }
  #Define variable to write to the output
  var_rain.bc <- ncvar_def("pr", "mm/day", list(lon1, lat2, time), missval=NA,longname="Biases Corrected Precipitation using Quantile Mapping", prec="float", verbose=FALSE)
  
  #Create the blank netcdf file 
  fname <- print(paste("pr.RCP45.imodel.",m,".nc",sep=""))
  new.nc <- nc_create(fname, list(var_rain.bc),force_v4=FALSE, verbose=FALSE)
  ncvar_put(new.nc,var_rain.bc,host45,start=c(1,1,1), count=c(nx,ny,nt), verbose=FALSE)
  nc_close(new.nc)
  
  fname <- print(paste("pr.RCP85.imodel.",m,".nc",sep=""))
  new.nc <- nc_create(fname, list(var_rain.bc),force_v4=FALSE, verbose=FALSE)
  ncvar_put(new.nc,var_rain.bc,host85,start=c(1,1,1), count=c(nx,ny,nt), verbose=FALSE)
  nc_close(new.nc)
  
}