;   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"

Dnum="03"
Dsize="1.33"
case="20141125"
year="2014"
month="11"
day="25"
starthour="18"
model="RAP"
PBL="MYNN2"
micro="THOM"
dom="NE"
SKEBS_string=""

begin
;

; Switch color map env var to use colormap saved locally
system (" export NCARG_COLORMAPS=/D2/ryan/ncl/colormaps:$NCARG_ROOT/lib/ncarg/colormaps ")

diri = "../Research/"+case+"_Case/"+model+"_"+PBL+"_"+micro+"_"+dom+"/wrfprd/"
fils = systemfunc ("csh -c 'cd " + diri + " ; ls *wrfout*d03*'")

do j = 121,dimsizes(fils)-1,1

  fili = fils(j)

; The WRF ARW input file.
  a = addfile(diri+fili,"r")

; We generate plots, but what kind do we prefer?
;  type = "x11"
 type = "png"
; type = "ps"
; type = "ncgm"

; Set some basic resources
  res = True
;  res@MainTitle                   = Dsize+" km WRF, 10 cm Sim Reflec at z = 2 km, "+model+" IC/BC"
  res@MainTitle                   = Dsize+" km dx WRF 1.0 km Sim Reflectivity"

  pltres = True
  mpres = True
  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpProvincialLineColor  = "Black"
  mpres@mpNationalLineColor    = "Black"
  mpres@mpUSStateLineColor     = "Black"
  mpres@mpGridLineColor        = 0
  mpres@mpLimbLineColor        = "Black"
  mpres@mpPerimLineColor       = "Black"
  mpres@mpUSStateLineThicknessF = 4.0
  mpres@mpNationalLineThicknessF = 6.0
  mpres@mpGeophysicalLineThicknessF = 4.0

; Plotting options for limited area consistent across domains
;  mpres@mpLimitMode = "Corners"
;  mpres@mpRightCornerLatF = 45.2
;  mpres@mpLeftCornerLonF = -70.
;  mpres@mpLeftCornerLatF = 40.2
;  mpres@mpRightCornerLonF = -76.

  cmap = (/"gray100", "black", "gray100", "lightskyblue1", "lightskyblue", "cornflowerblue", "royalblue3", "blue3", "lawngreen", "green2", "green3", "green4", "darkgreen", "yellow", "gold", "darkgoldenrod1", "darkorange", "darkorange2", "firebrick1", "firebrick3", "firebrick", "firebrick4", "maroon4", "violet"/) 

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

; 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

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

  it = 0			
 
;  wks = gsn_open_wks(type,"../Research/"+case+"_Case/"+model+"_"+PBL+"_"+micro+"_"+dom+"_"+SKEBS_string+"/wrf_dbz_"+SKEBS_string+"_"+dom+"_"+times(it))
  wks = gsn_open_wks(type,"../Research/"+case+"_Case/"+model+"_"+PBL+"_"+micro+"_"+dom+"/rainbow_1_min_wrf_dbz_"+times(it))
;  wks = gsn_open_wks(type,"../Research/"+case+"_Case/"+model+"_"+PBL+"_"+micro+"_9-3-1_"+SKEBS_string+"/"+Dsize+"_km_wrf_dbz_"+dom+"_"+times(it))
;  gsn_define_colormap(wks,cmap)       ; Overwrite the standard color map
  gsn_define_colormap(wks,"NCV_bright_white")       ; Overwrite the standard color map
  

     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 constant intercept parameters
; (as in early Reisner-2), with values of 8x10^6, 2x10^7, and 4x10^6 m^-4,
; for rain, snow, and graupel, respectively.

     mdbz = wrf_user_getvar(a,"mdbz",it) 
     dbz = wrf_user_getvar(a,"REFL_10CM",it) 
     z = wrf_user_getvar(a,"z",it)
     lat2d = wrf_user_getvar(a,"XLAT",it) 
     lon2d = wrf_user_getvar(a,"XLONG",it)  

     radar_Z = 10.^(dbz/10.)

     height = 1000.	; 1 km
     radar_Z_plane  = wrf_user_intrp3d(radar_Z,z,"h",height,0.,False)
     dbz_plane = 10.*log10(radar_Z_plane)

     opts = res
     pltres@gsnFrame = False               ;;--- Holds framing the plot until called // KS
     pltres@gsnDraw  = False               ;;--- Holds framing the plot until called // KS
     pltres@gsnPanelPlot = True
     opts@sfXArray = lon2d 
     opts@sfYArray = lat2d                         
     opts@cnFillOn = True  
     opts@ContourParameters = (/ 0., 60., 0.5/)
     contour = wrf_contour(a,wks,dbz_plane(:,:),opts)     ; plot only lowest level
     pltres@LatLonOverlay = False
     plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

draw(plot)

; Construct a polyline to show cross-section location

  pltres                 = True
  pltres@gsFillColor     = "black"
  pltres@cnLineDrawOrder      = "Draw"
  pltres@gsLineThicknessF     = 10.0

  lat = (/44.4,  42.9/)
  lon = (/-71.3, -69.3/)
  dum = gsn_add_polyline(wks,plot,lon,lat,pltres)

; draw(plot)
  frame(wks)           ;;--- closes wks AFTER adding the polyline to your plot! // KS

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


  end do        ; END OF TIME LOOP

; Reset NCARG_ROOT for future uses
system (" export NCARG_ROOT=/opt/ncl ")

end