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

begin
;---Read data
  filename = "chamba_tmean.nc"
  a = addfile(filename,"r")

  temp  = a->t
  dims  = dimsizes(temp)
  ntim  = dims(0)
  nlat  = dims(1)
  nlon  = dims(2)
  ntim5 = ntim/5
; 
; Reshape temperature array so we can easily average groupings of
; five years at a time.
;
  temp_5year_avg = dim_avg_n(reshape(temp,(/ntim5,5,nlat,nlon/)),1)
  temp_5year_avg@long_name = "Average temperature for 5 years"
  temp_5year_avg@units     = "degC"
  copy_VarCoords(temp(0,:,:),temp_5year_avg(0,:,:))
  printVarSummary(temp_5year_avg)

  temp_avg = dim_avg_n(temp,0)
printVarSummary(temp_avg)
exit

;---Get just the years out of the time array for plotting purposes.
  newtime = cd_calendar(temp&time,0)    ; Convert to N x 6 array (year,month,...,second)
  years = newtime(:,0)
  year5 = reshape(years,(/ntim5,5/))    ; This is for constructing a plot title later

  wks = gsn_open_wks("x11","chamba_contour")

  res                   = True

  res@gsnMaximize       = True     ; maximize plot in frame

;---Comment these two lines if you want to see the individual plots before the panel plot.
  res@gsnDraw           = False
  res@gsnFrame          = False

  res@gsnAddCyclic      = False
  res@cnFillOn          = True
  res@cnFillMode        = "RasterFill"
  res@cnLinesOn         = False
  res@cnLineLabelsOn    = False
  res@lbLabelBarOn      = False     ; will add in panel
  res@mpFillOn          = False
  res@mpOutlineOn       = False
;  res@lbOrientation     = "Vertical"

  res@mpMinLatF         = min(temp_5year_avg&lat)
  res@mpMaxLatF         = max(temp_5year_avg&lat)
  res@mpMinLonF         = min(temp_5year_avg&lon)
  res@mpMaxLonF         = max(temp_5year_avg&lon)
  res@pmTickMarkDisplayMode = "Always"

  res@gsnLeftString         = ""
  res@gsnRightString        = ""
  res@tiMainFontHeightF     = 0.03  
  res@pmTitleZone           = 4       ; Moves main title down

  mnmxint = nice_mnmxintvl( min(temp_5year_avg), max(temp_5year_avg), 18, False)
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF       = mnmxint(0)
  res@cnMaxLevelValF       = mnmxint(1)
  res@cnLevelSpacingF      = mnmxint(2)
;  res@cnLevelSpacingF      = mnmxint(2)/2.    ; Decrease spacing for twice as many levels

  lnres = True
  plots   = new(ntim5,graphic)
  ind_adm = new(ntim5,graphic)
  do n=0,ntim5-1
    res@tiMainString = year5(n,0) + " to " + year5(n,4)
    printMinMax(temp_5year_avg(n,:,:),0)
    plots(n)   = gsn_csm_contour_map(wks,temp_5year_avg(n,:,:),res)
;---Add shapefile outlines for India (downloaded from gadm.org)
    ind_adm(n) = gsn_add_shapefile_polylines(wks,plots(n),"IND_adm/IND_adm3.shp",lnres)
  end do

  pres = True
  pres@gsnMaximize        = True
  pres@gsnPanelLabelBar   = True
  pres@gsnPanelMainString = "Average temperature (degC)"
  pres@lbLabelFontHeightF = 0.012
  gsn_panel(wks,plots,(/2,4/),pres)
end