;***********************************************************
; climo_6.ncl
;
; Concepts illustrated:
;   - Calculating climatologies spanning user-specified years
;   - Getting the indices where data falls in a particular range
;   - Using "landsea_mask" to create a land/sea mask for your dataset
;   - Using "mask" to set land or ocean values in your data to missing
;   - Converting "string" time values using cd_calendar
;   - Converting "short" data to "float"
;   - Selecting a different color map
;   - Masking contour lines so they don't go through contour labels
;   - Paneling three plots vertically on a page
;   - Overlaying line contours on filled contours
;
;***********************************************************
; Generate Sample Climatologies as in Fig 1
;***********************************************************
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
   plev    = 850                              ; 850hPa

   ymStrt  = 198001                           ; start yyyymm
   ymLast  = 199912                           ; last  
   yrStrt  = ymStrt/100
   yrLast  = ymLast/100

   wksType = "ps"
   wksName = "x11"                          ; ": "+yrStrt+"_"+yrLast
  
   diri    = "/home/anyuola/"         ; input directory   


;***********************************************************
; Read year-month PRC; Create climatology for desired period
;***********************************************************

   fili    = "GPCP.nc"
   f       = addfile (diri+fili , "r")                          

   time    = f->time                          ; days since ...
   yyyymm  = cd_calendar(time, -1)            ; convert to yyyymm

   ntStrt  = ind(yyyymm.eq.ymStrt)            ; index start
   ntLast  = ind(yyyymm.eq.ymLast)            ; index last 
   prc     = f->precip(ntStrt:ntLast,:,:)     ; year-month
                                              
   prcClm  = clmMonTLL( prc)                  ; Climatology
                                              ; (12,lat,lon)
   delete(yyyymm)
   delete(time)
   delete(prc)

   printVarSummary( prcClm )
   printMinMax( prcClm, True )


;************************************************
; plots ; Compute 6-month climatologies
;************************************************
                                              ; (2,6)
  season   = (/ (/ 5, 6, 7, 8, 9,10/) \       ; May-Oct [summer]
              , (/ 1, 2, 3, 4,11,12/) /)      ; Nov-Apr [winter]
  i_season = season - 1                       ; NCL indices
  season_label = (/ "May-Oct Climatology", "Nov-Apr Climatology"/)

  plot = new ( 1, "graphic")
  wks = gsn_open_wks(wksType, wksName)        
  gsn_define_colormap(wks,"amwg")

;************************************************
; resource list for first (color) data array
;************************************************
  res                      = True
  res@gsnDraw              = False              ; don't draw
  res@gsnFrame             = False              ; don't advance frame
  res@gsnStringFontHeightF = 0.0125             ; make larger than default

  res@lbOrientation        = "Vertical"         ; vertical label bar
  res@lbLabelFontHeightF   = 0.01               ; make labels larger
  res@pmLabelBarOrthogonalPosF = -0.025         ; move closer to plot
  res@lbLabelAutoStride    = True               ; optimal label stride

;---This resource not needed in NCL V6.1.0
  res@gsnSpreadColors      = True               ; use full range of colors

  res@mpCenterLonF         = 180.               ; center plot at 180		res@mpLimitMode         = "LatLon" 
  res@mpMinLonF            = 21.               ; select a subregion
  res@mpMaxLonF            = 54.
  res@mpMinLatF            = -12. 
  res@mpMaxLatF            =  23

  res@mpLandFillColor      = "background"       ; color of land
  res@cnFillDrawOrder      = "Predraw"         

;***********************************************
; resource list for panel
;************************************************
  resP                     = True                ; modify the panel plot
  resP@gsnMaximize         = True                ; make large
 ;resP@gsnPaperOrientation = "portrait"          ; force "portrait"

 
 

  do ns=0,1                                     ; 2 seasons

     res@cnLevelSpacingF      = 20.0               ; set contour spacing
     

     res@gsnLeftString        = "GPCP PRC: ("+prcClm@units+")"
     res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
    ;res1@cnLevels             = (/1,2,4,6,8,10,12,14/)
     res@cnLevels             = ispan(1,14,1)             
     prcSeaClm = dim_avg_Wrap( prcClm(lat|:,lon|:,month|i_season(ns,:)) )
    ; plot   = gsn_csm_contour_map_overlay(wks, prcSeaClm,res)
     delete(res@cnLevelSelectionMode)
     delete(res@cnLevels)

     resP@txString  = season_label(ns) +": "+ yrStrt+"-"+yrLast
     gsn_panel(wks,plot,(/1,1/),resP)               ; now draw as one plot
   end do
end