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

begin
;---Open WRF output file and read data
  dir      = "./"
  filename = "wrfout_d01_2015-09-16_00:00:00"
  a = addfile(dir + filename,"r")

  times  = wrf_user_list_times(a)  ; get times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;---Just look at first time step.
  it = 0
  print("Working on time: " + times(it) )

;---Read temperature and terrain height off file
  tc    = wrf_user_getvar(a,"tc",it)     ; T in C
  lat2d = wrf_user_getvar(a,"XLAT",it)   ; latitude
  lon2d = wrf_user_getvar(a,"XLONG",it)  ; longitude
  dims = dimsizes(tc)

;---Start the graphics
  wks = gsn_open_wks("x11","wrf_latlonoverlay")

;---Set some basic plot options
  opts = True
  opts@MainTitle                   = "REAL-TIME WRF"

  pltres = True
  mpres  = True

  opts@TimeLabel = times(it)   ; Set valid time to use on plots

;----------------------------------------------------------------------
; Plot full domain first.
;----------------------------------------------------------------------
  opts@cnFillOn = True  

  contour = wrf_contour(a,wks,tc(it,:,:),opts)
  plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

;----------------------------------------------------------------------
; Plot partial domain.
;----------------------------------------------------------------------
  ;opts@sfXArray = lon2d
  ;opts@sfYArray = lat2d

  ;contour = wrf_contour(a,wks,tc(it,:,:),opts)

;---Set special resource to indicate we are using XLAT/XLONG coordinates.
  pltres@LatLonOverlay = True

;---Zoom in on map, which we can do because we're using lat/lon coordinates.
  mpres@mpLeftCornerLatF  =  -12
  mpres@mpRightCornerLatF =  22.5
  mpres@mpLeftCornerLonF  = 21
  mpres@mpRightCornerLonF = 54
  

  plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
end