;***********************************************
; climo_3.ncl
;
; Concepts illustrated:
;   - Calculating monthly climatologies
;   - Calculating interannual variability
;   - Calculating standard deviations of monthly means
;   - 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"

begin
;*******************************************************
; Read the file and Compute monthly Climatolgy(36years)
;*******************************************************
  a    = addfile("GPCP.nc","r")
  prec = a->precip(0:431,:,:)                                 ; (time,lat,lon)
  printVarSummary(prec)
  time = a->time                                ; (time) ==> YYMM
  ntim = dimsizes(time)
  precClm = clmMonTLL(prec)
  printVarSummary(precClm)
;****************************************************************************
;Below Selects the year to plot but for this case it is ommited
;************************************************************************
  
  ;ymStrt  = 198001                           ; start yyyymm
   ;ymLast  = 199912                           ; last  
   ;yrStrt  = ymStrt/100
   ;yrLast  = ymLast/100                           ; start yyyymm
;********************************************************************************  
   
   wksType = "x11"
   wksName = "test"                          ; ": "+yrStrt+"_"+yrLast
  
   
;************************************************
; plots ; Compute 3-month climatologies(season)
;************************************************
                                              ; (2,6)
  season   = (/ (/ 2, 3, 4/) \       ; March-May[MAM]
              , (/ 5, 6, 7/)   \      ; June-August [JJA]
	      , (/ 9, 10, 11/) /)         ; Oct-Dec [OND]
  i_season = season - 1                       ; NCL indices
  season_label = (/ "March-May", "June-August","October-December"/)
 plot = new ( 3, "graphic")
  wks = gsn_open_wks(wksType, wksName)        
  gsn_define_colormap(wks,"BlueYellowRed")
  
  print(dimsizes(season))
       
     
;************************************************
; Compute the MAM,JJA,OND climatology using a function in contributed.ncl
;************************************************
 ; precClm = month_to_season(prec,"MAM")            ; monthly climatology
  prec_ts =prec(time|:,lat|:,lon|:)
  prec_ts = runave(prec,3,0)
  ;************************************************
; create colors
;*************************************************
  wks = gsn_open_wks("x11","test")             ; open a workstation
    
      
 colors = (/ (/255,255,255/),(/0,0,0/),(/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)


             
  gsn_define_colormap(wks, colors)              ; generate new color map
;************************************************
; create panel plots
;*************************************************
  plot   = new ( 3, graphic)                    ; create graphical array

  res                      = True               ; plot options desired

  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@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels             = (/0.0,0.5,1.0,5.0,10,15,20,25\ 	; set unequal contour levels
                              ,30,35/)
                             
  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"
  mpOutlineOn		  =True
  mpOutlineSpecifiers     =True
  res@mpMinLonF            = 21.               ; select a subregion
  res@mpMaxLonF            = 54.
  res@mpMinLatF            = -12. 
  res@mpMaxLatF            = 22.5
  
  
  res@mpLandFillColor      = "background"       ; color of land
  res@cnFillDrawOrder      = "Predraw"    
  ;*********************Masking*********************************
  res@mpAreaMaskingOn = True  
   res@mpMaskAreaSpecifiers = "Land" 
   res@mpOceanFillColor = 0 
   res@mpFillDrawOrder = "Postdraw"  
    
  
      resMap= True 
      resMap@mpPerimOn = False 
     resMap@mpGridAndLimbOn = False 
    resMap@gsnTickMarksOn = False 
    resMap@mpDataSetName = "Earth..3" 
   resMap@mpDataBaseVersion = "MediumRes" 
   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
  

  do ns=0,2                                 ; loop over the seasons
         
         res@gsnCenterString   = season(ns,:)+":"+time(0)/100 +"-"+ time(ntim-1)/100
     plot = gsn_csm_contour_map_ce(wks,prec_ts(ns,:,:), res)  ; create plot
  end do

  gsn_panel(wks,plot,(/2,2/),resP)
 
 end