; --------------------------------------------------------------------------- ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  Program Name: combine_qpe_bias.ncl                                         ;
;                                                                             ;
;  National Water Model (NWM) WRF-hydro forcing engine is developed by        ;
;  National Center for Atmospheric Research (NCAR), under the sponsorship     ;
;  of National Water Center (NWC).                                            ;
;                                                                             ;
;   Team Members:                                                             ;
;     NCAR Staff: Linlin Pan, Wei Yu, and David Gochis                        ;
;      NWC/OWP Staff: Brian Cosgrove, Zhengtao Cui, Cham Pham, and James Taft ;
;                                                                             ;
;  This is a ncl program to perform Precipitation bias correction             ;
;                    and combination of HRRR, RAP and MRMS                    ;
;                                                                             ;
;  Input files: 1) hrrrFile: HRRR rainfall data,                              ;
;               2) rapFile : RAP rainfall data                                ;
;               3) mrmsFile : MRMS rainfall data                              ;
;               4) hrrrWgt : HRRR weighting function                          ;
;               5) rapWgt : RAP weighting function                            ;
;               6) mrmsWgt : MRMS weighting function                          ;
;               7) hrrrBias : HRRR bias                                       ;
;               8) rapBias : RAP bias                                         ;
;               9) mrmsBias : MRMS bias                                       ;
;  Output: combined rainrate                                                  ;
;                                                                             ;
;                                                                             ;
; For non-fatal errors output is witten to $DATA/logs                         ;
;                                                                             ;
; Author(s)/Contact(s): Linlin Pan, lpan@ucar.edu                             ;
; Origination                                                  Oct., 2015     ;
;                                                                             ;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;

begin
           LOCKFILE = getenv("LOCKFILE") ; HRRR data
           hrrrFile = getenv("hrrrFile") ; HRRR data
           rapFile = getenv("rapFile") ; RAP data
           mrmsFile = getenv("mrmsFile") ; MRMS data
           hrrrWgt = getenv("hrrrWgt") ; HRRR wgt
           rapWgt = getenv("rapWgt") ; RAP wgt
           mrmsWgt = getenv("mrmsWgt") ; MRMS wgt
           hrrrBias = getenv("hrrrBias") ; HRRR Bias
           rapBias = getenv("rapBias") ; RAP Bias
           mrmsBias = getenv("mrmsBias") ; MRMS bias
           outFile = getenv("outFile") ; outFile data
           forcingDir = getenv("forcDir") ; outFile data
;           
           fhrrrF = addfile(hrrrFile ,"r") 
           frapF = addfile(rapFile,"r") 
           fmrmsF = addfile(mrmsFile,"r") 
;
           fhrrrW = addfile(hrrrWgt ,"r") 
           frapW = addfile(rapWgt,"r") 
           fmrmsW = addfile(mrmsWgt,"r") 
;
           fhrrrB = addfile(hrrrBias ,"r") 
           frapB = addfile(rapBias,"r") 
           fmrmsB = addfile(mrmsBias,"r") 
;           system("rm -f "+outFile)

           fout = addfile(outFile + ".nc" ,"w")

        ndim=dimsizes(dimsizes(fhrrrF->RAINRATE))
        if(ndim.eq.2)then
           rainH = fhrrrF->RAINRATE(:,:)
        else
           rainH = fhrrrF->RAINRATE(0,:,:)
        end if
        ndim=dimsizes(dimsizes(frapF->RAINRATE))
        if(ndim.eq.2)then
           rainR = frapF->RAINRATE(:,:)
        else
           rainR = frapF->RAINRATE(0,:,:)
        end if
        ndim=dimsizes(dimsizes(fmrmsF->precip_rate))
        if(ndim.eq.2)then
           rainM = fmrmsF->precip_rate(:,:)
        else
           rainM = fmrmsF->precip_rate(0,:,:)
        end if

; ### temporary hardwired
       rainH(3308:3338,:)=rainR(3308:3338,:)
       rainH(209:239,:)=rainR(209:239,:)

          ;rainwH = fhrrrW->APCP_P0_L102_GLC0(:,:)
          if(isfilevar(fhrrrW,"APCP_P0_L102_GLC0")) then
            rainwH = fhrrrW->APCP_P0_L102_GLC0(:,:)
          else
            rainwH = fhrrrW->POP_P0_L102_GLC0(:,:)
          end if


           rainbH = fhrrrB->POP_P0_L102_GLC0(:,:)
           rainbR = frapB->POP_P0_L102_GLC0(:,:)
           rainbM = fmrmsB->POP_P0_L102_GLC0(:,:)
           rainwR = frapW->POP_P0_L102_GLC0(:,:)
           rainwM = fmrmsW->POP_P0_L102_GLC0(:,:)
	   printVarSummary(rainH)
	   printVarSummary(rainwH)
	   printVarSummary(rainbH)
	   printVarSummary(rainR)
	   printVarSummary(rainwR)
	   printVarSummary(rainbR)
	   printVarSummary(rainM)
	   printVarSummary(rainwM)
	   printVarSummary(rainbM)
;	   rainbR=where(ismissing(rainbR),rainR,rainbR)
;	   rainbH=where(ismissing(rainbH),rainR,rainbH)
;	   rainbM=where(ismissing(rainbM),rainR,rainbM)
        if(ndim.eq.2)then
	   rainrr=rainR*rainwR*rainbR
	   rainhh=rainH*rainwH*rainbH
	   rainmm=rainM*rainwM*rainbM
        else
	   rainrr=rainR(0,:,:)*rainwR*rainbR
	   rainhh=rainH(0,:,:)*rainwH*rainbH
	   rainmm=rainM(0,:,:)*rainwM*rainbM
        end if
	   rain=rainrr+rainhh+rainmm
	   rain = where(ismissing(rainmm).and.(.not.ismissing(rainhh)),rainrr+rainhh,rain)
	   rain = where(ismissing(rainhh).and.(.not.ismissing(rainmm)),rainrr+rainmm,rain)
        if(ndim.eq.2)then
	   rain = where(ismissing(rain),rainR,rain)
	   rain = where((rain.gt.1).or.(rain.lt.0),rainR,rain)
        else
	   rain = where(ismissing(rain),rainR(0,:,:),rain)
	   rain = where((rain.gt.1).or.(rain.lt.0),rainR(0,:,:),rain)
        end if


        if(ndim.eq.2)then
	   fout->precip_rate = (/rain(:,:)/)
        else
	   fout->precip_rate(0,:,:) = (/rain(:,:)/)
        end if
	   fout->precip_rate@units="mm s-1"
           fout->precip_rate@long_name="Surface Precipitation Rate"

	   system("mv -f " + outFile + " " + forcingDir + "/.")

           system("rm -f "+LOCKFILE)

exit
end
status_exit(1)