; NCL script to read in Gridmet PET climatology and ET fraction and create ETactual files according to the method of G. Senay (USGS).  
;
;*******************************************************
; Code written by Joe Barsugli joseph.barsugli@noaa.gov
; from the NOAA/ESRL PSD
;********************************************************
;
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"         	
; +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
begin

; Define beginning and end years for loop over years
begyr = 2000
endyr = 2017

datadir = "/Projects/CMIP5C/OBS_DATA/SSEBOP/SSEBOP_NETCDF/"   ; This directory is the ROOT for input and output files

; Read in Gridmet PET 8-day climatology netCDF file
gridmet_file_name = datadir + "Gridmet_PET/" + "gridmet_pet_8day_climatology.nc"
gmc = addfile(gridmet_file_name,"r") ;  Open Gridmet Climatology PET file (8 day -- 46 entries)
pet_climo = gmc->pet_pm_8d    ; Read in PET climatology 


do yr = begyr,endyr   ; loop over years
print("Year:" + yr)

etf_file_name = datadir + "ETfractions_8day/" + "etfca.8d." + yr + ".nc"  ; open ET fraction file for year yr
etfin = addfile(etf_file_name,"r")
etf = etfin->ETf
time = etfin->time       ;read in time (referenced to 2000-01-01)

  
 vdims=getfilevardimsizes(etfin,"ETf")
 nlats = vdims(1)
 nlons = vdims(2)
print(vdims)
 
daybeg=etfin->daybeg     ; read in auxiliary time coordinates
dayend = etfin->dayend
YMDbeg = etfin->YMDbeg
YMDend = etfin->YMDend

eta = etf * pet_climo   ; Calculation of ETa:  Note that the grids of etf and pet_climo are aligned.
 print("multiplication done")
copy_VarCoords(etf,eta);    Copy over the coordinates and dimensions 
eta@units = "mm"
eta@long_name = "Actual Evapotranspiration (mm)"
daybeg@long_name = "day at beginning of 8-day period"
daybeg@units =  "days since 2000-1-1 0:0:0" 
dayend@long_name = "day at end of 8-day period"
dayend@units =  "days since 2000-1-1 0:0:0" 
YMDbeg@long_name = "Year, Month, Day  at beginning of 8-day period"
YMDend@long_name = "Year, Month, Day  at end of 8-day period"

; Set netcdf file format and compression level

setfileoption("nc","Format","NetCDF4")
setfileoption("nc","CompressionLevel",3)

;  Open the Output File
etaout = addfile(datadir + "ETa_8day/" + "test/" + "ETactual.8day." + yr + ".nc","c")
print("addfile done")
; Set global attributes
globalAtt = True
globalAtt@description = "Actual Evapotranspiration computed from 8-day ET fraction multiplied by 8-day climatology of PET from Gridmet data"
fileattdef(etaout,globalAtt)

; define variable in output file and set chunk sizes (NOte the default NCL behavior was ridiculous and broke downstream applications)
filedimdef(etaout,(/"time","lat","lon"/),vdims,(/True,False,False/))
print("filedimdef done")
filevardef(etaout, "ETa", "float", (/  "time", "lat", "lon" /) )
print("filevardef done")
chunkSizes = (/ 1, 569, 1324 /)
filevarchunkdef(etaout, "ETa", chunkSizes)

print("chunking done")
; Output the data
etaout->ETa = eta
etaout->daybeg = daybeg
etaout->dayend = dayend
etaout->YMDbeg = YMDbeg
etaout->YMDend = YMDend


end do  ; end of year loop

end