load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin
setvalues NhlGetWorkspaceObjectId() 
        "wsMaximumSize" : 100000000000 
end setvalues
;Input and output directory
 diri="/Users/ahaque/Desktop/project/"
 imgdir="/Users/ahaque/Desktop/project/"

;Reading the file for 21 Century Precipitation

infile=addfile("/Users/ahaque/Desktop/project/pr_day_CCSM4_rcp85_r1i1p1_20510101-21001231_sa.nc","r")

;Reading file for 20 century precipitation
infile2=addfile("/Users/ahaque/Desktop/project/pr_day_CCSM4_historical_r1i1p1_19510101-20001231_sa.nc","r")

;;365day calendar!
j1=152; jun1
j2=274; Sept30
;;time indices associated with julian day params
t1=j1-1 ;; 8xdaily will be reduced to 1 daily avg
t2=j2-1 ;; 8xdaily will be reduced to 1 daily avg
nyrs=50

pr21c=infile->pr
pr20c=infile2->pr
delete(pr21c@lat)
delete(pr21c@lon)
delete(pr20c@lat)
delete(pr20c@lon)
printVarSummary(pr20c)
printVarSummary(pr21c)

tsiz=nyrs*(t2-t1+1);;number of times after seasonal subset (JJA)
latsiz=dimsizes(infile->lat)
lonsiz=dimsizes(infile->lon)

lat=infile->lat
lon=infile->lon
lon=infile->lon
lat=infile->lat

percentile=fspan(0.,100.,tsiz)
percentile!0="percentile"
percentile&percentile=percentile
percentile@units="%"
percentile@long_name="percentile of JJAS daily precipitation distribution"
skip=0 ;count leap days

inpr=new((/tsiz,latsiz,lonsiz/),"float")
inpr!0="time"
inpr!1="lat"
inpr!2="lon"
inpr@lat=lat
inpr@lon=lon

inpr20c=inpr


;compute daily total precip -- using 06z to 06z to define the daily sum
tmp21c=runave_Wrap(pr21c(lat|:,lon|:,time|:),8,0)  ;8 times step running mean 
prdy21c=tmp21c(time|2::8,lat|:,lon|:) ;; 3rd time step is 12z --> this is the 06z-06z time avg
prdy21c=prdy21c*8 ;;mm/24hrs

tmp20c=runave_Wrap(pr20c(lat|:,lon|:,time|:),8,0)  ;8 times step running mean 
prdy20c=tmp20c(time|2::8,lat|:,lon|:) ;; 3rd time step is 12z --> this is the 06z-06z time avg
prdy20c=prdy20c*8 ;;mm/24hrs

skip=0 ;count leap days

;pull out summer only
do iyr=0,nyrs-1,1

; no leap 365 day calendar
 istrt=floattoint(iyr*365+t1)  ;input time indices
 iend=istrt+t2-t1
 ostrt=iyr*(1+t2-t1) ;;output time indices
 oend=ostrt+t2-t1
 print(istrt+" "+iend+" "+ostrt+" "+oend)
 inpr(ostrt:oend,:,:)=prdy21c(istrt:iend,:,:); mm/3hrs
 inpr20c(ostrt:oend,:,:)=prdy20c(istrt:iend,:,:); mm/3hrs
end do


printVarSummary(inpr)

;sort precip data
inprsorted21c=inpr21c
sort21c=dim_pqsort_n(inprsorted21c,2,0)
inprsorted21c!0="percentile"
inprsorted21c&percentile=percentile


inprsorted20c=inpr20c
sort20c=dim_pqsort_n(inprsorted20c,2,0)
inprsorted20c!0="percentile"
inprsorted20c&percentile=percentile

printVarSummary(inprsorted20c)
printVarSummary(inprsorted)

outname=diri+"pr.CMIP5.pct.CCSM4.all.nc"
system("rm "+outname)
outfile=addfile(outname,"c")
outfile->lat=lat
outfile->lon=lon
outfile->pr20c=inprsorted20c
outfile->pr21c=inprsorted21c
delete(outfile)

exit
end