; NCL script to read in ETa monthly and create cumulative ETa

;*******************************************************
; 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
kma = 12 ;   number of months to accumulate
print("number of months to accumulate = " +kma)

datadir = "/Projects/CMIP5C/OBS_DATA/SSEBOP/SSEBOP_NETCDF/MonthlyET/"   ; This directory is the ROOT for input and output files
fname = "modisSSEBopETactual.2000-2016.all.nc"


; Open ETa file 
eta_file_name = datadir + fname
print("file = " + eta_file_name)
etain = addfile(eta_file_name,"r") ;  Open ETa monthly (all data) file
eta = etain->ETa    ; Read in ETa Monthly file
time = etain->time       ;read in time (referenced to 2000-01-01)

d=dimsizes(eta)
nmo = d(0);
eta_accum = eta   ; Initialize array for output (this could be big! -- May need to run on machine with > 32 GB memory

do m = 0,nmo-1   ; loop over years
  m1 = dim_max((/m-kma+1,0/));
  print ("m1 = " + m1 + " m = " +m )
  eta_accum(m,:,:)=kma*dim_avg_n_Wrap(eta(m1:m,:,:),0) ;  sum over nmo timesteps. USe average * kma so that an occasional missing value will return a reasonable magnitude value. 
end do

eta_accum(0:kma-2,:,:) = eta_accum@_FillValue   ; set the first kma-1 months to _FillValue because they are incomplee accumulations

; Set netcdf file format and compression level

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

;  Open the Output File
fnameout = "modisSSEBopETactual.2000-2016." + kma + "month.nc"
etaout = addfile(datadir + fnameout,"c")

; Set global attributes
globalAtt = True
globalAtt@description = "Actual Evapotranspiration from G. Senay(USGS) SSEBOP method:  " + kma + " month accumulation"
fileattdef(etaout,globalAtt)

; Output the data
etaout->ETa = eta_accum

end