;----------------------------------------------------------------------
; gpcp_1.ncl
;
; Concepts illustrated:
;   - Reading a CNRM-CM5 and CESM1-BGC historical and RCP8.5file
;   - Calculate daily annual cycle and panel them averaged data over a region
;----------------------------------------------------------------------
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"
;**************************************************************
; User Input
;**************************************************************
   pltName  = "Annual Cycle=Historical+rcp8.5"
   pltType  = "x11"

   pltTitle = "Annual Cycle"
   regName  = "South Asia"

   latS     =  18.0    ; Region
   latN     =  25.0
   lonL     =  75.0
   lonR     = 90.0
;***************************************************************
; End User Input
;***************************************************************
   a            =addfile("/Projects/CMIP5A/Daily/historical/CNRM-CM5/r1i1p1/pr_day_CNRM-CM5_historical_r1i1p1_18500101-200512.nc","r")
   time1        = a->time
   pr1          = a->pr(:,{latS:latN},{lonL:lonR})   ; (time,lat,lon)  (0,1,2)
   clat1        = latRegWgt(pr1&lat, typeof(pr1), 0)   ; create lat weights
   pr1Area      = wgt_areaave_Wrap(pr1, clat1, 1, 0)    ; weighted area avg at each time

   ymdStrt      = 18500101 ; start yyyymmdd
   ymdLast      = 20051231       ;last
   yrStrt       = ymdStrt/10000
   yrLast       = ymdLast/10000
   ymdPlot2     = 20051231
   nhar         = 4  ; no of fourier comp
   var          = "pr"

   b             =addfile("/Projects/CMIP5A/Daily/rcp85/CNRM-CM5/r1i1p1/pr_day_CNRM-CM5_rcp85_r1i1p1_20060101-21001231.nc", "r")
   time2         = b->time
   pr2           = b->pr(:,{latS:latN},{lonL:lonR})   ; (time,lat,lon)  (0,1,2)
   clat2         = latRegWgt(pr2&lat, typeof(pr2), 0)   ; create lat weights
   pr2Area       = wgt_areaave_Wrap(pr2, clat2, 1, 0)    ; weighted area avg at each time

ymdStrt       = 20060101 ; start yyyymmdd
   ymdLast       = 21001231       ;last
   yrStrt        = ymdStrt/10000
   yrLast        = ymdLast/10000
   ymdPlot2      = 21001231

  c              =addfile ( "/Projects/CMIP5A/Daily/historical/CESM1-BGC/r1i1p1/pr_day_CESM1-BGC_historical_r1i1p1_18500101-20051231.nc", \
"r")
   time3         = c->time
   pr3           = c->pr(:,{latS:latN},{lonL:lonR})   ; (time,lat,lon)  (0,1,2)
   clat3         = latRegWgt(pr3&lat, typeof(pr3), 0)   ; create lat weights
   pr3Area       = wgt_areaave_Wrap(pr3, clat3, 1, 0)    ; weighted area avg at each time
   ymdStrt       = 18500101 ; start yyyymmdd
   ymdLast       = 20051231       ;last
   yrStrt        = ymdStrt/10000
   yrLast        = ymdLast/10000
   ymdPlot2      = 20051231

  d             =addfile ( "/Projects/CMIP5A/Daily/rcp85/CESM1-BGC/r1i1p1/pr_day_CESM1-BGC_rcp85_r1i1p1_20060101-21001231.nc", "r")
  time4         = d->time
  pr4           = d->pr(:,{latS:latN},{lonL:lonR})   ; (time,lat,lon)  (0,1,2)
  clat4         = latRegWgt(pr4&lat, typeof(pr4), 0)   ; create lat weights
  pr4Area       = wgt_areaave_Wrap(pr4, clat4, 1, 0)    ; weighted area avg at each time

  ymdStrt       = 20060101 ; start yyyymmdd
  ymdLast       = 21001231       ;last
  yrStrt        = ymdStrt/10000
  yrLast        = ymdLast/10000
  ymdPlot2      = 21001231
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
   nday         = 366
   pday1        = new ( nday, typeof(pr1), pr1@_FillValue)
    do nd=0,nday-2
      iday  = ind(ddd.eq.(nd+1))
      pday1(nd) = avg(pr1Area(iday))
      delete(iday)                                ; may change next ireration
     end do
    pday1(nday-1) = 0.5*(pday1(0)+pday1(nday-2))
    pr1AvgRegion = avg(pr1Area)                    ; regional avg

   nday          = 366
   pday2         = new ( nday, typeof(pr2), pr2@_FillValue)
    do nd=0,nday-2
      iday  = ind(ddd.eq.(nd+1))
      pday2(nd) = avg(pr2Area(iday))
      delete(iday)                                ; may change next ireration
     end do
    pday2(nday-1) = 0.5*(pday2(0)+pday2(nday-2))
    pr2AvgRegion = avg(pr2Area)                    ; regional avg

   nday          = 366
   pday          = new ( nday, typeof(pr3), pr3@_FillValue)
   do nd=0,nday-2
     iday  = ind(ddd.eq.(nd+1))
     pday3(nd) = avg(pr3Area(iday))
     delete(iday)                                ; may change next ireration
   end do
  pday3(nday-1) = 0.5*(pday3(0)+pday3(nday-2))
  pr3AvgRegion  = avg(pr3Area)                    ; regional avg
nday          = 366
  pday4         = new ( nday, typeof(pr4), pr4@_FillValue)
  do nd=0,nday-2
     iday  = ind(ddd.eq.(nd+1))
     pday4(nd) = avg(pr4Area(iday))
     delete(iday)                                ; may change next ireration
  end do
  pday4(nday-1) = 0.5*(pday4(0)+pday4(nday-2))

  pr4AvgRegion  = avg(pr4Area)                    ; regional avg
;************************************************
; Create smooth annual cycle
;************************************************

   cf1                = ezfftf1(pday1)
   cf1(:,3 )          = 0.75*cf1(:,3)                        ; arbitrary wgts
   cf1(:,4 )          = 0.50*cf1(:,4)
   cf1(:,5 )          = 0.25*cf1(:,5)
   cf1(:,6:)          = 0.0                                 ; no contribution
   pday1_smth         = ezfftb1 (cf1, cf1@xbar)                ; reconstruct

   cf2                = ezfftf2(pday2)
   cf2(:,3 )          = 0.75*cf2(:,3)                        ; arbitrary wgts
   cf2(:,4 )          = 0.50*cf2(:,4)
   cf2(:,5 )          = 0.25*cf2(:,5)
   cf2(:,6:)          = 0.0                                 ; no contribution
   pday2_smth         = ezfftb2 (cf2, cf2@xbar)                ; reconstruct

   cf3                = ezfftf3(pday3)
   cf3(:,3 )          = 0.75*cf3(:,3)                        ; arbitrary wgts
   cf3(:,4 )          = 0.50*cf3(:,4)
   cf3(:,5 )          = 0.25*cf3(:,5)
   cf3(:,6:)          = 0.0                                 ; no contribution
   pday3_smth         = ezfftb3 (cf3, cf3@xbar)                ; reconstruct

   cf4                = ezfftf4(pday4)
   cf4(:,3 )          = 0.75*cf4(:,3)                        ; arbitrary wgts
   cf4(:,4 )          = 0.50*cf4(:,4)
   cf4(:,5 )          = 0.25*cf4(:,5)
   cf4(:,6:)          = 0.0                                 ; no contribution
   pday4_smth         = ezfftb4(cf4, cf4@xbar)                ; reconstruct
;************************************************
; Create plot
;************************************************
   ntim              = dimsizes(time)
   PDAY1             = new( (/2,nday/), typeof(pday1), "No_FillValue")
   PDAY1(0,:)        = pday1
   PDAY1(1,:)        = pday1_smth
   PDAY2             = new( (/2,nday/), typeof(pday2), "No_FillValue")
   PDAY2(0,:)        = pday2
   PDAY2(1,:)        = pday2_smth
   PDAY3             = new( (/2,nday/), typeof(pday3), "No_FillValue")
   PDAY3(0,:)        = pday3
   PDAY3(1,:)        = pday3_smth
   PDAY4             = new( (/2,nday/), typeof(pday4), "No_FillValue")
   PDAY4(0,:)        = pday4
   PDAY4(1,:)        = pday4_smth

   pltPath= pltDir+pltName+"."+pltType

   wks               = gsn_open_wks("x11", "Annual Cycle")

   res               = True        ; plot mods desired
   res@gsnMaximize   = True
  ;res@vpHeightF     = 0.4         ; change aspect ratio of plot
   res@vpWidthF      = 0.8
res@trXMinF       = 0           ; max value on x-axis
   res@trXMaxF       = nday+1      ; max value on x-axis
   res@xyLineThicknesses = (/2.0, 2.0/)
   res@xyLineColors      = (/"blue","red"/)  ; change line color
   res@xyMonoDashPattern = True

   res@vpXF              = 0.1         ; start plot at x ndc coord
   res@tmYLFormat        = "f"

   res@tiXAxisString     = "day of year"     ; x

   ;res@gsnYRefLine       = prcAvgRegion           ; create a reference line
   ;res@gsnRightString    = "Areal Mean="+sprintf("%4.2f", prcAvgRegion)+" mm/day"
   ;res@gsnLeftString     = regName

;plotting
   res@gsnLeftString            = "South Asia"
   res@gsnLeftString                    = "CNRM-CM5-Hist"
   res@gsnRightString                   = "Annual Cycle:("+pr1Avg@units+")"
   res@gsnLeftStringFontHeightF         = 0.015
   res@gsnRightStringFontHeightF        = 0.015
   res@gsnCenterStringFontHeightF       = 0.015
   PDAY1(0)                             = gsn_csm_xy(wks, PDAY1, res)
	res@gsnLeftString     	= "South Asia"
   res@gsnRightString    	= "CNRM-CM5-RCP85"
   res@tiYAxisString     	= "Regional: pr (mm/day)"    ; y-axis label
   res@gsnYRefLine       	= avg(pr2Avg)   ; create a reference line
   res@gsnCenterString   	= "Areal Mean="+sprintf("%4.2f", avg(pr2Avg))+" mm/day"
   PDAY2(1)                	= gsn_csm_xy(wks, PDAY2,res)

   res@gsnLeftString     	= "South Asia"
   res@gsnRightString    	= "CESM1-BGC-Hist"
   res@tiYAxisString     	= "Regional: pr (mm/day)"    ; y-axis label
   res@gsnYRefLine       	= avg(pr3Avg)   ; create a reference line
   res@gsnCenterString   	= "Areal Mean="+sprintf("%4.2f", avg(pr3Avg))+" mm/day"
   PDAY3(2)              	= gsn_csm_xy(wks,PDAY3,res)

   res@gsnLeftString     	=  "South Asia"
   res@gsnRightString    	=  "CESM1-BGC-RCP8.5"
   res@tiYAxisString     	=  "Regional: pr (mm/day)"    ; y-axis label
   res@gsnYRefLine        	=  avg(pr4Avg)   ; create a reference line
   res@gsnCenterString   	=  "Areal Mean="+sprintf("%4.2f", avg(pr4Avg))+" mm/day"
   PDAY4(3)                     	=  gsn_csm_xy(wks,PDAY4,res)
 end do
; panelling

   resP = True
   resP@txString               = " Differences of Annual Cycle of Mean Precipitation by Historical and RCP8.5 from CMIP 5 Models"
   resP@gsnFrame           = False
   ;resP@txFontHeightF   = 0.0165
   resP@gsnMaximize      = True                ; maximize pl
   gsn_panel(wks, PDAY,(/2,2/),resP) ; not draw as one plot
end