; 1 mm/day    = 28.4 W/m2 

  lonL   = 117.5
  lonR   = 122.5
  latS   = 11.25
  latN   = 18.75

  YRSTRT = 1979    ; sa tellite era
  YRLAST = 2013

;---Open file 
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/addfiles.shtml
  diri   = "/Users/shea/Data/CDC/"
  fili   = "olr.day_19740601_20131231.mean.nc"
  pthi   = diri+fili
  f      = addfile(pthi, "r")

;---What indices correspond to the YRSTRT-YRLAST period
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/cd_calendar.shtml
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/ind.shtml
  TIME_file  = cd_calendar(f->time, 0)
  iTIME_file = ind(TIME_file(:,0).ge.YRSTRT .and. \
                   TIME_file(:,0).le.YRLAST)
  niTIME     = dimsizes(iTIME_file)
  itStrt     = iTIME_file(0)                  
  itLast     = iTIME_file(niTIME-1)
  delete([/ TIME_file, iTIME_file, niTIME/])

;---Read data for specified YRSTRT-YRLAST period only
;---https://www.ncl.ucar.edu/Document/Functions/Contributed/short2flt.shtml
  x    = short2flt(f->olr(itStrt:itLast,{latS:latN},{lonL:lonR}))
  printVarSummary(x)
  printMinMax(x,0)
  print("==========")

;---All times spanning desired period
  TIME  = cd_calendar(x&time, 0)     
  NTIM  = dimsizes(x&time)

;---Extract assorted 'date' information in different formats
;---https://www.ncl.ucar.edu/Document/Functions/Built-in/day_of_year.shtml
  yyyy = toint(TIME(:,0))
  mm   = toint(TIME(:,1))
  dd   = toint(TIME(:,2))
  mmdd = mm*100+dd
  ddd  = day_of_year(yyyy, mm, dd)       ; needed for clmDayTLL
  yyyyddd  = yyyy*1000  + ddd 
  yyyymmdd = yyyy*10000 + mmdd

;---Raw daily climatology
;---http://www.ncl.ucar.edu/Document/Functions/Contributed/clmDayTLL.shtml
  xClmDay    = clmDayTLL(x, yyyyddd)     ; daily climatology at each grid point
  printVarSummary(xClmDay)
  printMinMax(xClmDay,0)
  print("==========")

;---Smooth daily climatologies  
;---http://www.ncl.ucar.edu/Document/Functions/Contributed/smthClmDayTLL.shtml
  xClmDay = smthClmDayTLL(xClmDay, 2)  
  printVarSummary(xClmDay)
  printMinMax(xClmDay,0)
  print("==========")

;---Calculate anomalies from smoothed climatology
;---http://www.ncl.ucar.edu/Document/Functions/Contributed/calcDayAnomTLL.shtml
  xAnom  = calcDayAnomTLL (x, yyyyddd, xClmDay)     
  xAnom@long_name = "Anomalies from Smooth Daily Climatology"
  printVarSummary(xAnom)
  printMinMax(xAnom, 0)
  print("==========")

;---Use straight average. No need to areally weight for so small an area
;---https://www.ncl.ucar.edu/Document/Functions/Contributed/dim_avg_n_Wrap.shtml
  xts     = dim_avg_n_Wrap(xAnom,(/1,2/))
  printVarSummary(xts)
  printMinMax(xts, 0)
  print("==========")

;---Specify yearly daily periods
  mmStrt    = 5               ; May 15
  ddStrt    = 15
  mmLast    = 7               ; July 15
  ddLast    = 15

;---Extract temporal indices for desired daily periods
  mmddStrt  = mmStrt*100 + ddStrt
  mmddLast  = mmLast*100 + ddLast
  dddStrt   = day_of_year(1901, mmStrt, ddStrt)  
  dddLast   = day_of_year(1901, mmLast, ddLast)  
  dddSpan   = dddLast - dddStrt +1               ; # no. of days in segment
  print("Start: mmddStrt="+ mmddStrt+"  dddStrt="+dddStrt)
  print("Last:  mmddLast="+ mmddLast+"  dddLast="+dddLast)
  print("Span:  dddSpan ="+ dddSpan )
  print("==========")

;---Create array to be plotted: (year,day)
  yrStrt = yyyy(0)
  yrLast = yyyy(NTIM-1)
  year   = ispan(yrStrt, yrLast, 1)   ; left axis
  year!0 = "year" 
  printVarSummary(year)
  printMinMax(year,0)
  print("==========")

  day    = fspan(dddStrt, dddLast, dddSpan)
  day!0  = "day"
  day@long_name = "day of year"
  printVarSummary(day)
  printMinMax(day,0)
  print("==========")

  nyear = dimsizes(year)
  data  = new( (/nyear,dddSpan/), typeof(x), getFillValue(x))

;---Fill Array with values
  do nyr=0,nyear-1    ; loop over each year: extract and assign day segment
     yyyymmddStrt = year(nyr)*10000 + mmStrt*100 + ddStrt  ; current year start
     yyyymmddLast = year(nyr)*10000 + mmLast*100 + ddLast  ; current year last
     immdd        = ind(yyyymmdd.ge.yyyymmddStrt .and. \   ; current year indices
                        yyyymmdd.le.yyyymmddLast)
    ;print("nyr="+nyr+"  yyyymmddStrt="+yyyymmddStrt+" yyyymmddLast="+yyyymmddLast\
    ;                +"  dimsizes(immdd)="+dimsizes(immdd))
     data(nyr,:)  = xts(immdd) 
  end do
 ;print("==========")

;---Add meta data
  data!0    = "year"
  data&year =  year
  data!1    = "day" 
  data&day  =  day
  data@long_name = "OLR anomalies"
  data@units     = x@units

  printVarSummary(data)
  printMinMax(data,0)
  print("==========")

;==============================================
; create plot
;============================================== 
  pltType = "png"
  pltDir  = "./"
  pltName = "DayYear.olr"
  pltPath = pltDir+pltName
  
  wks  = gsn_open_wks (pltType, pltPath  ) 

  res               = True            ; Plot mods desired
  res@gsnMaximize   = True            ; Maximize plot in frame
  res@cnFillOn      = True            ; Turn on contour fill
  res@cnFillPalette = "BlueWhiteOrangeRed"
  res@lbOrientation = "vertical"
  res@tiMainString  = YRSTRT+"-"+YRLAST+": mmdd="+       \
                      sprinti("0.4i%",mmddStrt)  + "-" + \ 
                      sprinti("0.4i%",mmddLast) 
  symMinMaxPlt (data,18,False,res)
  res@cnLevelSpacingF = 10
  plot = gsn_csm_contour(wks, data, res ) 

 ;print(mmdd(0:200)+"   "+ddd(0:200))

;---Set some month-day axis labels
  res@tmXBMode              = "Explicit"
  res@tmXBValues            = (/ 135  , 152 , 166  , 182 ,  196 /)                  
  res@tmXBLabels            = (/"5/15","6/1","6/15","7/1","7/15"/)
  plot = gsn_csm_contour(wks, data, res )