;:::: Calculate the regional mean from daily data
;::::Load Libraries::::::::::::::::::::::::::::::::::::::::::::::::
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"
;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
;Begin
   regName  = "South Asia"

   lat1     =  0.0    ; Region
   lat2     =  40.0
   lon1     =  60.0    ;
   lon2     =  100.0    ;
  
;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

 ;::::Open file and read data
  dir       = "/home/ahaque/"
  filename  = "GPCP_1DD_v1.2_199610-201312.nc"
  a         = addfile(dir+filename,"r")
  prec      = a->PREC;
  ymdStrt = 19961001 ; start yyyymmdd
  ymdLast = 20131231       ;last
  ymdPlot2= 20131231
  yyyymmdd= (a->time)/100
iday    = ind(yyyymmdd.eq.ymdPlot2)   ; (time, lat, lon)(0,1,2)
prec     = a->PREC;(iday:iday,:,:)
time    = a->time   ; time:units = days since
TIME    = cd_calendar(time,0)   ; type of float
year    = floattointeger( TIME(:,0) )
month   = floattointeger( TIME(:,1) )
day     = floattointeger( TIME(:,2) )
ddd     = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd

 b=addfile("/Projects/CMIP5A/Daily/historical/CNRM-CM5/r1i1p1/pr_day_CNRM-CM5_historical_r1i1p1_18500101-20051231.nc","r")
  pr0     = b->pr ;
time1    = b->time   ; time:units = days since
TIME    = cd_calendar(time,0)   ; type of float
year    = floattointeger( TIME(:,0) )
month   = floattointeger( TIME(:,1) )
day     = floattointeger( TIME(:,2) )
ddd     = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd

 c=addfile("/Projects/CMIP5A/Daily/historical/CESM1-BGC/r1i1p1/pr_day_CESM1-BGC_historical_r1i1p1_18500101-20051231.nc", "r")
  pr1     = c->pr;
time2   = c->time   ; time:units = days since
TIME    = cd_calendar(time,0)   ; type of float
year    = floattointeger( TIME(:,0) )
month   = floattointeger( TIME(:,1) )
day     = floattointeger( TIME(:,2) )
ddd     = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd

 d=addfile ( "/Projects/CMIP5A/Daily/historical/NorESM1-M/r1i1p1/pr_day_NorESM1-M_historical_r1i1p1_18500101-20051231.nc", "r")
   pr2     = d->pr;  
time3    = d->time   ; time:units = days since
TIME    = cd_calendar(time,0)   ; type of float
year    = floattointeger( TIME(:,0) )
month   = floattointeger( TIME(:,1) )
day     = floattointeger( TIME(:,2) )
ddd     = day_of_year(year, month, day)
yyyyddd = year*1000 + ddd

;::::Convert units kg/(m2-s) ;to mm/day;
  prec=prec*1.
  prec@units= "mm/day" ;to mm/day;
  printVarSummary(prec)
  printMinMax(prec,1)
  pr0 = pr0*86400.
  pr0@units = "mm/day"
  printVarSummary(pr0)
  printMinMax(pr0,1)
  pr1 = pr1*86400.
  pr1@units = "mm/day"
  printVarSummary(pr1)
  printMinMax(pr1,1)
  pr2 = pr2*86400.
  pr2@units = "mm/day"
  printVarSummary(pr2)
  printMinMax(pr2,1)

;Compute regional mean precipitation
precAvg = wgt_areaave_Wrap(prec(:,{lat1:lat2},{lon1:lon2}),1.0,1.0,0)
pr0Avg = wgt_areaave_Wrap(pr0(:,{lat1:lat2},{lon1:lon2}),1.0,1.0,0)
pr1Avg = wgt_areaave_Wrap(pr1(:,{lat1:lat2},{lon1:lon2}),1.0,1.0,0)
pr2Avg = wgt_areaave_Wrap(pr2(:,{lat1:lat2},{lon1:lon2}),1.0,1.0,0)

printVarSummary(precAvg) 
printVarSummary(pr0Avg) 
printVarSummary(pr1Avg)
printVarSummary(pr2Avg) 

;***************************************************************
;:::Setting resources for Plotting 
wks = gsn_open_wks("x11","Regional Mean Pr")
gsn_define_colormap(wks,"rainbow")
res1                 = True            ; Variable to hold plot options
res1@gsnMaximize     = True            ; Maximize plot in frame.
res1@gsnSpreadColors = True            ; Span full color map 
res1@gsnLeftString   = ""
res1@gsnRightString  = ""
res1@gsnAddCyclic    = False
res1@gsnDraw         =  False                   ; do not draw the plot
res1@gsnFrame        =  False                   ; do not advance the frame 

res1@cnLinesOn       = False           ; turn off contour line color
res1@cnLineLabelsOn  = False
res1@cnFillOn        = True            ; Contour fill on
res1@lbLabelBarOn    = False ; added procedure to do a single label bar for the panel

res1@mpLimitMode     = "LatLon" 
res1@mpMinLatF       = 0 ;min(lat2d)
res1@mpMaxLatF       = 40 ;max(lat2d)
res1@mpMinLonF       = 60 ;min(lon2d)
res1@mpMaxLonF       = 100 ;max(lon2d)

res1@pmTickMarkDisplayMode = "Always" ; turn on tick marks
;:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: 
;plotting procedure
 plot = new(4,graphic)
 do i=0,3
      res1@gsnLeftString = "GPCP: PR"
      res1@gsnRightString = "Regional Mean:("+precAvg@units+")"      
      res1@gsnLeftStringFontHeightF = 0.015     
      res1@gsnRightStringFontHeightF = 0.015
      plot(0) = gsn_csm_contour_map_ce(wks,precAvg,res1)
   
      res1@gsnLeftString = "CNRM-CM5: PR"
      res1@gsnRightString = "Regional Mean:("+pr0Avg@units+")"
      plot(1) = gsn_csm_contour_map_ce(wks,pr0Avg,res1) 

      res1@gsnLeftString = "CESM1-BGC: PR"
      res1@gsnRightString = "Regional Mean:("+pr1Avg@units+")"
      plot(2) = gsn_csm_contour_map_ce(wks,pr1Avg,res1) 
 
      res1@gsnLeftString = "NorESM1-M: PR"
      res1@gsnRightString = "Regional Mean:("+pr2Avg@units+")"
      plot(3) = gsn_csm_contour_map_ce(wks,pr2Avg,res1)
 end do

;::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
; panelling 
 
  resP = True 
  resP@txString = " Regional Mean Precipitation by Historical in South Asia from CMIP 5 Models"
  resP@gsnFrame= False
  resP@txFontHeightF = 0.0165
  resP@gsnPanelLabelBar = True   ; add common colorbar
  resP@lbOrientation   = "Vertical"      ; Move labelbar to side
  resP@gsnMaximize    = True                ; maximize plots
  resP@lbTitleString = "mm/d" 
  resP@lbLabelAutoStride   = True
  resP@lbLabelFontHeightF = 0.015
  resP@lbBoxLinesOn = False
  gsn_panel(wks,plot,(/2,2/),resP) ; not draw as one plot
  txres =True
  txres@gsnFrame= False  ; do not advance the frame
  txres@txFontHeightF = 0.0321
  gsn_text_ndc(wks,"Source: Azizul Haque, NOAA ESRL",087,024, txres)
 frame(wks)
end