;***********************************************
; 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
;
;***********************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
; 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
;************************************************
  a    = addfile("precip.mon.total.v2018.nc","r")
  prc  = a->precip                                 ; (time,lat,lon)
  time = a->time                                ; (time) ==> YYMM
  ntim = dimsizes(time)

  months = (/"January", "February", "March", "April" \
            ,"May", "June", "July", "August"         \
            ,"September", "October", "November"      \
            ,"December" /)
;************************************************
; Compute the climatology using a function in contributed.ncl
;************************************************
  prcClm = clmMonTLL (prc)            ; monthly climatology
  prcStd = stdMonTLL (prc)            ; monthly interannual variability
;************************************************
; create colors
;*************************************************
  wks = gsn_open_wks("png","climo")             ; send graphics to PNG file

  colors = (/ "azure1","beige","lavender" \
            ,"PaleGreen","SeaGreen3","LightYellow" \
            ,"Yellow","HotPink","Red"/)         ; choose colors for color map
;************************************************
; create panel plots
;*************************************************
  plot   = new (4 , graphic)                    ; create graphical array

  res                      = True               ; plot options desired

  res@cnFillOn             = True               ; turn on color fill
  res@cnFillPalette        = colors             ; set color map
  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.5, 1.0, 2.0 \ ; set unequal contour levels
                              , 3.0, 4.0, 5.0 \
                              , 7.5,10.0 /)

  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@gsnPanelMainString  = "CPC Merged Prc: Climatology"   ; common title
  resP@gsnMaximize         = True               ; maximize image
 ;resP@gsnPanelLabelBar    = True               ; Add common label bar 

  i = -1                                        ; Climatologies
  do nmo=0,11,3                                 ; loop over the months
     i = i+1
     res@gsnCenterString   = months(nmo)+":"+time(0)/100 +"-"+ time(ntim-1)/100
     plot(i) = gsn_csm_contour_map(wks,prcClm(nmo,:,:), res)  ; create plot
  end do

  gsn_panel(wks,plot,(/2,2/),resP)
;************************************************
; create second panel plot
;*************************************************
  i = -1                                        ; Std Deviations
  do nmo=0,11,3                                 ; loop over the months
     i = i+1
     res@gsnCenterString   = months(nmo)+":"+time(0)/100 +"-"+ time(ntim-1)/100
     plot(i) = gsn_csm_contour_map(wks,prcStd(nmo,:,:), res)  ; create plot
  end do

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

end