;*************************************************
; WRF: RPRECIPITATION: Total, Cumulus and non-cumulus prc
;************************************************
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/wrf/WRF_contributed.ncl"
begin
;************************************************
; open file 
; Read Cumulus (rinc) and Non-cumulus (rainnc) prc
;************************************************
  f       = addfile("wrfout_d01_000000.nc", "r")
  rainc   = f->RAINC                 ; (Time, south_north, west_east)
  rainnc  = f->RAINNC     
                                    
  times   = chartostring(f->Times)  ; convert to type string [plot]
  ntim    = dimsizes(times)         ; # time steps
;************************************************
; Use NCL operator > to make sure all values >=0.0
; Sum components and assign attributes
;************************************************
  rainc   = rainc  > 0.0         
  rainnc  = rainnc > 0.0
  rainTot = rainc + rainnc          
  rainTot@description = "Total Precipitation"
  rainTot@units       =  rainc@units

;************************************************
; create plots: create colormap using named colors
;               unequal contour levels
;************************************************
  wks = gsn_open_wks("ps" ,"WRF_lc")           ; ps,pdf,x11,ncgm,eps
  colors = (/"white","black"          \        ; {for/back}ground
            ,"white","azure"          \ 
            ,"green","palegreen","yellowgreen", "greenyellow" \
            ,"yellow","goldenrod","orange","orangered"        \
            ,"red","deeppinK", "violet","darkviolet"          \
            ,"blueviolet","blue"                              /)
  gsn_define_colormap(wks, colors)  
  
  res                       = True             ; plot mods desired
;;res@gsnMaximize           = True             ; uncomment to maximize size
  res@gsnSpreadColors       = True             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels
  
  res@cnLevelSelectionMode  = "ExplicitLevels" ; explicit [unequal] cn levels
  res@cnLevels              = (/0,0.1,1,2.5,5,7.5,10,15,20,25,37.5,50,75,100,125,150/)
  res@cnFillMode            = "RasterFill"
  res@lbOrientation         = "Vertical"       ; default is horizontal
;************************************************
; Use WRF_contributed procedure to set map resources
;************************************************
  WRF_map_c(f, res, 0)                   ; reads info from file

;************************************************
; if appropriate, set True for native (direct) mapping (faster)
;************************************************
   res@tfDoNDCOverlay       = True

;************************************************
; associate the 2-dimensional coordinates to the variables for plotting
; only necessary if it is not a native projection [res@tfDoNDCOverlay=False]
;************************************************
  if (.not.res@tfDoNDCOverlay) then
      lat2d  = f->XLAT(0,:,:)                 ; need for map limits        
      lon2d  = f->XLONG(0,:,:)      
    
      rainc@lat2d   = lat2d
      rainc@lon2d   = lon2d
      rainnc@lat2d  = lat2d
      rainnc@lon2d  = lon2d
      rainTot@lat2d = lat2d
      rainTot@lon2d = lon2d
  end if

;****************************************************************************
; create panel of different components
;****************************************************************************
  plts = new (3 , "graphic")                   ; 1d array to hold plots

  res@gsnDraw               = False            ; (a) do not draw 
  res@gsnFrame              = False            ; (b) do not advance 'frame'
  res@lbLabelBarOn          = False            ; (c) turn off individual lb's 
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP                      = True             ; modify the panel plot
  resP@gsnMaximize          = True             ; maximize panel area
  resP@gsnPanelRowSpec      = True             ; specify 1 top, 2 lower level  
  resP@gsnPanelLabelBar     = True             ; add common colorbar
  resP@pmLabelBarWidthF     = 0.85             ; make label wider 
  resP@lbLabelFontHeightF   = 0.015            ; default 0.02 [demo make smaller]

  nt = 12                                      ; demo only one time 
;;do nt=0,ntim-1                               ; uncomment to loop over all times
     plts(0) = gsn_csm_contour_map(wks,rainTot(nt,:,:),res)
     plts(1) = gsn_csm_contour_map(wks,rainnc(nt,:,:),res)
     plts(2) = gsn_csm_contour_map(wks,rainc(nt,:,:),res)

     resP@txString = f@TITLE+": "+times(nt)      
     gsn_panel(wks,plts,(/1,2/),resP)               ; now draw as one plot
;;end do
end