;***********************************************
; climo_3.ncl
;
; Concepts illustrated:
;   - Calculating monthly climatologies
;   - Calculating seasonal totals for 36yrs
;   - Creating a color map using named colors
;   - Paneling four plots on a page
;   - Adding a common labelbar to paneled plots
;   - Explicitly setting contour levels
;   - Changing the center longitude for a cylindrical equidistant projection
;   - Turning off the individual contour labelbar
;   - Adding a common title to paneled plots
;   - Turning off map fill
;   - Creating a center subtitle
;
;***********************************************
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
;*******************************************************
; Read the file for GPCC
;*******************************************************
  ymStrt = 198101 
  ymLast = 199012
  a    = addfile("GPCC.nc","r")
 
  TIME = a->time
  YYYYMM = cd_calendar(TIME, -1)
  iStrt  = ind(YYYYMM.eq.ymStrt)
  iLast  = ind(YYYYMM.eq.ymLast)

  prec = a->precip(iStrt:iLast,:,:)              ; (time,lat,lon)
  delete(prec@units)
  delete(prec@long_name)
  printMinMax(prec,True)
  printVarSummary(TIME)
  
  printVarSummary(prec)
   time = a->time                                ; (time) ==> YYMM
  ntim = dimsizes(time)
  precClm = clmMonTLL(prec)
   printVarSummary(precClm)

;*******************************************************
; Read the file for GPCP
;*******************************************************
  ymStrt = 198101 
  ymLast = 199012
  b    = addfile("GPCP.nc","r")
  TIME := b->time
  ;time =b->time
  ;delete(TIME)
  printVarSummary(TIME)
  yyyymm = cd_calendar(time, -1)
  iStrt  = ind(yyyymm.eq.ymStrt)
  iLast  = ind(yyyymm.eq.ymLast)

  prec1 = b->precip(0:431,:,:)              ; (time,lat,lon)
  delete(prec1@units)
  delete(prec1@long_name)
  printMinMax(prec1,True)
  
  printVarSummary(prec1)
   ;time = b->time                             ; (time) ==> YYMM
  ;ntim = dimsizes(time)
  prec1Clm = clmMonTLL(prec1)
   printVarSummary(prec1Clm)

;************************************************
; Compute the monthly climatology using a function in contributed.ncl
;************************************************
  precMonClm = clmMonTLL(prec)   ; (12,:,:)
  precSeaClm = runave_n_Wrap(precMonClm, 3, -1, 0)  ; unweighted seasonal average
  printVarSummary(precSeaClm)
  
   prec1MonClm = clmMonTLL(prec1)   ; (12,:,:)
  prec1SeaClm = runave_n_Wrap(prec1MonClm, 3, -1, 0)  ; unweighted seasonal average
  printVarSummary(prec1SeaClm)
  
   season_index = (/3,6,10/)  ; MAM, JJA, OND 
  season_label = (/ "MAM", "JJA","OND"/)
  
;************************************************
; create colors
;*************************************************
  wks = gsn_open_wks("x11", "seasons")          ; open  a work station
                                          
   
 colors = (/ (/255,255,255/),(/255,255,255/), (/244,255,244/), \
  (/217,255,217/), (/163,255,163/), (/106,255,106/), \
  (/43,255,106/), (/0,224,0/), (/0,134,0/),(/255,255,0/),\
  (/255,127,0/) /) * 1.0 ; we multiply by 1 to make colors float
  
    colors = colors/255.               ; normalize (required by NCL)

 ;************************************************
; create panel plots
;*************************************************
 plot   = new ( 6, graphic)                    ; create graphical array
 
    res                      = True               ; plot options desired
    
  ;res@cnFillDrawOrder      = "Predraw"    
  res@cnFillOn             = True               ; turn on color fill
  res@cnInfoLabelOn        = False              ; turn off contour info label
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLineLabelsOn       = False              ; turn off line labels
  res@cnFillPalette        =colors
  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
 res@cnLevels        =(/0.1,0.2,0.4,0.8,1.6,3.2,6.4,12.8,25,30,40,45,60,75,83,100/)

                             
  res@mpFillOn             = False              ; turn off gray continents
  res@mpCenterLonF         = 180                ; Centers the plot at 180 
 
  res@lbLabelBarOn         = False              ; No single label bar

  res@gsnDraw              = False
  res@gsnFrame             = False

  resP                     = True               ; panel options
  resP@txString            = "Seasons"   ; common title
  resP@gsnMaximize         = True               ; maximize image
  resP@gsnPanelLabelBar    = True               ; Add common label bar 
  
  ;****************ZOOM INTO YOUR LOCATION*********************
  
  res@mpLimitMode           = "LatLon"
 ;res@mpOutlineOn		  =True
 ;res@mpOutlineSpecifiers     =True
  res@mpMinLonF            = 21.               ; select a subregion
  res@mpMaxLonF            = 54.
  res@mpMinLatF            = -12. 
  res@mpMaxLatF            = 22.5
  
  
  ;res@mpLandFillColor      = "background"       ; color of land
  
  ;*********************Masking*********************************
 

      res@mpFillDrawOrder = "Postdraw"  
    
  
      resMap= True 
      resMap@mpPerimOn = False 
     resMap@mpGridAndLimbOn = False 
    resMap@gsnTickMarksOn = False 
    resMap@mpDataSetName = "Earth..3" 
     resMap@mpOutlineOn = True 
   
  
  res@mpFillOn                    = False        ; turn off gray fill
  res@mpOutlineBoundarySets       = "National"   ; turn on country boundaries
   
  resMap@mpFillOn = True 

  resMap@mpOceanFillColor = "transparent" 
  resMap@mpLandFillColor = "transparent" 
  resMap@mpInlandWaterFillColor = "transparent" 
   
   res@mpDataBaseVersion = "MediumRes"            ; choose higher resolution
  
    i=0 
   do ns=0,dimsizes(season_index)-1
    
    res@gsnCenterString   = season_label(ns)+":"+(ymStrt/100) +"-"+(ymLast/100)
	
    ;res@gsnCenterString   = ""
    plot(i) = gsn_csm_contour_map_ce(wks,precSeaClm(ns,:,:), res)  ; create plot
    
    i=1
    do ns=1,dimsizes(season_index)-1
           
    res@gsnCenterString   = season_label(ns)+":"+(ymStrt/100) +"-"+(ymLast/100)
	
    plot(i+1) = gsn_csm_contour_map_ce(wks,prec1SeaClm(ns,:,:), res)  ; create plot
    
    
   end do
end do
  gsn_panel(wks,plot,(/3,3/),resP)
 
 end