;   Example script to produce dbz plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   November 2008

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("../wrfout_d01_2000-01-24_12:00:00.nc","r")

; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_dbz2")
  gsn_define_colormap(wks,"WhViBlGrYeOrReWh")

; Set some basic resources
  res = True
  res@MainTitle                   = "REAL-TIME WRF"

  pltres = True
  mpres = True

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Which times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 1,ntimes-1               ; TIME LOOP
 
     print("Working on time: " + times(it) )
     res@TimeLabel = times(it)   ; Set Valid time to use on plots

; First get the variables we will need        
; Both dbz and mdbz will be calculated using intercept parameters
; for rain, snow, and graupel, which are consistent with
; Thompson, Rasmussen, and Manning (2004, Monthly Weather Review,
; Vol. 132, No. 2, pp. 519-542.)
;	First "1" in wrf_user_getvar
; Frozen particles that are at a temperature above freezing will be
; assumed to scatter as a liquid particle.
;	Second "1" in wrf_user_getvar

     mdbz = wrf_user_getvar(a,(/"mdbz","1","1"/),it) 
     dbz = wrf_user_getvar(a,(/"dbz","1","1"/),it) 

     opts = res                         
     opts@cnFillOn = True  
     opts@ContourParameters = (/ 5., 75., 5./)
     contour = wrf_contour(a,wks,dbz(1,:,:),opts)    ; plot only lowest level
     plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
 
     contour = wrf_contour(a,wks,mdbz,opts)
     plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
 

   end do        ; END OF TIME LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
end