;***********************************************
; 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
 
 wks  = gsn_open_wks("png","shapefiles")  ; send graphics to PNG file


  mres                     = True
  mres@mpLimitMode         = "Corners"     ; corner method of zoom
  mres@mpLeftCornerLatF    = 29            ; left corner
  mres@mpLeftCornerLonF    = 47          ; left corner
  mres@mpRightCornerLatF   = 33            ; right corner
  mres@mpRightCornerLonF   = 52           ; right corner
  ;mres@mpProjection        = "LambertConformal"    ; choose projection
  ;mres@mpLambertParallel1F = 33            ; first parallel
  ;mres@mpLambertParallel2F = 45            ; second parallel
  ;mres@mpLambertMeridianF  = -98           ; meridian
  mres@tfDoNDCOverlay      = True          ; native grid, no transform
; mres@tfDoNDCOverlay      = "NDCViewport" ; NCL 6.5.0 or later
  mres@gsnDraw             = False         ; don't draw yet
  mres@gsnFrame            = False         ; don't advance frame yet
  mres@gsnMaximize = False
  mres@pmTickMarkDisplayMode = "Always"    ; turn on tickmarks
  

  plot = gsn_csm_map(wks,mres)

  plres = True
  plres@gsLineColor = "black"
  lines_id = gsn_add_shapefile_polylines(wks,plot,"IRN_adm2.shp",plres)

;************************************************
; 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({34675:79227},:,:))            ; monthly climatology
 prcStd = stdMonTLL(prc({34674:79227},:,:))   ; 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 (2, 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,1.5, 2.0, 2.5 \ ; set unequal contour levels
                              , 3.0, 3.5, 4.0,4.5, 5.0 \
                              , 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.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 
  res@mpDataBaseVersion    = "MediumRes"
  res@mpDataSetName        = "Earth..4"
  ;res@mpOutlineBoundarySets = "GeophysicalAndIran"
   res@mpOutlineBoundarySets = "National"
  
  res@mpLimitMode           = "LatLon"
  res@mpMaxLatF                   = 33          ; choose subregion           
  res@mpMinLatF                   = 29
  res@mpMaxLonF                   = 52
  res@mpMinLonF                   = 47
  lnres                   = True         ; Set some line options
  lnres@gsLineThicknessF  = 2.0          ; default is 1.0
  ;id = gsn_add_shapefile_polylines(wks,plot,dir+shapefile_name,lnres)

  res@pmTickMarkDisplayMode = "Always"     ; nicer map tickmarks
; For each plot, customize resources for the labelbar and/or tickmarks.
;
  res@lbLabelBarOn   = False
  res@tmYROn         = False
  res@tmXBOn         = False
  res@tmXTOn         = False
; Turn off tickmarks.
   res@gsnTickMarksOn        = False


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

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

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

end