; NCL script to read in ETa for all years and compute ranking in the range of 1:Nyrs/(Nyrs+1)
;    
;
;*******************************************************
; 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 = 2016   ; This is the last complete year of data

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

; Read in all years of data
;eta_file_name = datadir + "MonthlyET/" + "modisSSEBopETactual." +  begyr + "-" + endyr + ".all.nc"  ; open ET  file that contains all complete years from the range of years
eta_file_name = datadir + "MonthlyET/" + "modisSSEBopETactual." +  begyr + "-" + endyr + ".7month.nc"  ; open ET  file that contains all complete years from the range of years
 etain = addfile(eta_file_name,"r")
 time = etain->time      ;read in time (referenced to 2000-01-01)
 lat = etain->lat
 lon = etain->lon

 Nyrs = endyr-begyr+1
 nmo = 12

 year = ispan(begyr,endyr,1)
 month = ispan(1,nmo,1)
  
 vdims=getfilevardimsizes(etain,"ETa")
 nlats = vdims(1)
 nlons = vdims(2)
 dims3d=(/Nyrs,nlats,nlons/)

 eta = new(dims3d,float)  ; Define new float for the slice that will be processed (a single calendar month)
 eta!0 = "year"
 eta&year = year
 eta!1 = "lat"
 eta&lat = lat
 eta!2 = "lon"
 eta&lon = lon
 
 irank = new((/Nyrs,nmo,nlats,nlons/),"byte");  ; Define output variable
 good = new((/nmo,nlats,nlons/),"byte");  ; Define output variablea
 etastd = new((/Nyrs,nmo,nlats,nlons/),"float");  ; Define output variable
 irank@_FillValue = integertobyte(99)
 irank@_ZeroValue = integertobyte(-10)


 do mm = 0,nmo-1
  print("month " + mm)
  do iyr = 0,16
    imm=iyr*12+mm
    eta(iyr,:,:) = etain->ETa(imm,:,:);
  end do
  
  ip = dim_pqsort_n(eta, 1, 0)  ;; return permuation  along time dimension but do not sort eta.
  ii=ispan(1,Nyrs,1)

  iranktemp = integertobyte(dim_pqsort_n(ip,1,0)+1)
  iranktemp = where(ismissing(eta), irank@_FillValue, iranktemp)
  iranktemp = where(eta .lt. 1, irank@_ZeroValue, iranktemp)

  irank(:,mm,:,:) =  iranktemp
  good(mm,:,:)=integertobyte(dim_num_n(.not.ismissing(eta),0))
  etastd(:,mm,:,:)=dim_standardize_n_Wrap(eta,0,0)

  
  
  ;do i = 0,nlats-1
  ;  do j = 0,nlons-1
  ;    irank(ip(:,i,j),mm,i,j)=integertobyte(ii);
  ;  end do
  ;end do
end do

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

year!0 = "year"
month!0 = "month"


irank!0 = "year"
irank&year = year
irank!1 = "month"
irank&month = month
irank!2 = "lat"
irank&lat = lat
irank!3 = "lon"
irank&lon = lon
irank@description = "Ranking from lowest to highest of the Actual ETa from Senay's SSEBop method; Ranking includes missing values.Only use ranks from 1 to ngood"

good@_FillValue = integertobyte(-99)
good = where(good.lt.0, good@_FillValue, good)
good!0 = "month"
good&month = month
good!1 = "lat"
good&lat = lat
good!2 = "lon"
good&lon = lon
good@description = "Number of non-missing values across years for each calendar month at each pixel"

etastd@_FillValue = 9999
etastd!0 = "year"
etastd&year = year
etastd!1 = "month"
etastd&month = month
etastd!2 = "lat"
etastd&lat = lat
etastd!3 = "lon"
etastd&lon = lon
etastd@description = "Standardized actual ET from Senay SSEBop method"


;copy_VarCoords(eta,ir);    Copy over the coordinates and dimensions 

; first write out year by year 
do k = begyr,endyr
  kk=k-begyr
  
foutyr = addfile("test2/MonthlyET_ranked_" + k + ".6month.nc","c")
foutyr->good = good
foutyr->irank = irank(kk,:,:,:)


foutyr2 = addfile("test2/MonthlyET_standardized_" + k + ".6month.nc","c")
foutyr2->good = good
foutyr2->ETa_std = etastd(kk,:,:,:)

end do

; the write out one big file
;fout = addfile("MonthlyET_ranked_2000-2016.nc","c")
;fout->good = good
;fout->irank = irank


;fout2 = addfile("MonthlyET_standardized_2000-2016.nc","c")
;fout2->good = good
;fout2->ETa_std = etastd



end