;----------------------------------------------------------------------
; This script shows how to create a basic plot of temperature (degC)
; calculated from a WRF file.
;
; This script is meant to be compared to one that uses wrf_xxxx
; functions to do the plotting. 
;----------------------------------------------------------------------

 filename = "wrfout_d02_2013-02-08_15:30:00"
 f = addfile(filename,"r")
  lats = (/  39.00,  42.50 /)
  lons = (/ -76.0, -70.0/)
; bounding box for February 08 2013


; radar = wrf_user_getvar(f,"radar",0)    ; Temperature (C)
 xlat2d = wrf_user_getvar(f,"XLAT",0)    ; Temperature (C)
 xlon2d = wrf_user_getvar(f,"XLONG",0)    ; Temperature (C)
  it = 5
  radar = wrf_user_getvar(f,"fft_zh",it)
  zdr = wrf_user_getvar(f,"fft_zdr",it)
  crs = wrf_user_getvar(f,"fft_crs",it)
  kdp = wrf_user_getvar(f,"fft_kdp",it)
;   refl_max = dim_max_n(radar,0)
; z    = wrf_user_getvar(f, "z",it)      ; grid point height
; refl_zkm=wrf_user_intrp3d(radar,z,"h",height,0.,False)


;printVarSummary(radar)

 wks = gsn_open_wks("pdf","wrf_demo_plot_radar_gsn_mod")

  res                = True
  res@pmTickMarkDisplayMode = "Always"
  res@mpFillOn              = False           ; turn off map fill
  res@mpOutlineDrawOrder    = "PostDraw"      ; draw continental outline last
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
  res@mpLimbLineColor        = "Black"
  res@mpPerimLineColor       = "Black"
  res@mpPerimLineThicknessF  = 10
  res@mpUSStateLineThicknessF = 10.0
  res@mpGeophysicalLineThicknessF = 10.0
  res@mpGridAndLimbOn = False

  res@cnFillOn       = True  
  res@cnLevelSelectionMode = "ExplicitLevels" ; use explicit levels
  res@cnLevels             = ispan(-10,55,5) ; set the contour levels
  cmap = read_colormap_file("NCV_jet")
  cmap(0,:) = 1.   ; Make first color white
  res@cnFillPalette  = cmap

; If you want the first color to be transparent rather than white, then use:

; cmap(0,3) = 0,0   ; Make first color transparent

 res@cnLinesOn      = False

;--- These are required for gsn_csm_contour_map plotting
;res = wrf_map_resources(f,res)
 res@tfDoNDCOverlay = True
 res@gsnAddCyclic   = False

 nl = 5     ; level index
;plot = gsn_csm_contour_map(wks,radar(n1:,:),res)
 dims3d = dimsizes(radar) 
 print(dims3d)
; loc(0,;) is west-east (x) ; loc(1,:) is south-north (y)
; subtract one since we want to use it as an index in NCL
  loc = wrf_user_ll_to_ij(f, lons, lats, True)
  x_start = loc(0,0) - 1
  x_end   = loc(0,1) - 1
  y_start = loc(1,0) - 1
  y_end   = loc(1,1) - 1
 print("y_start = " + y_start)
 print("x_start = " + x_start)
 print("y_end = " + y_end)
 print("x_end = " + x_end)
 res@mpMinLatF       = xlat2d(y_start,x_start)
 res@mpMaxLatF       = xlat2d(y_end,x_end)
 res@mpMinLonF       = xlon2d(y_start,x_start)
 res@mpMaxLonF       = xlon2d(y_end,x_end)
 radar_plane = radar(nl,y_start:y_end,x_start:x_end)
 xlat2d_sub = xlat2d(y_start:y_end,x_start:x_end)
 xlon2d_sub = xlon2d(y_start:y_end,x_start:x_end)
 radar_plane@lat2d = xlat2d_sub    ; Needed for plotting
 radar_plane@lon2d = xlon2d_sub


 plot = gsn_csm_contour_map(wks,radar_plane,res)