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-19_12: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("png","wrf_test")

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

  ;mpres@mpGeophysicalColor          = "black"
  mpres@mpNationalLineColor         = "black"
  mpres@mpUSStateLineColor          = "black"
  mpres@mpGeophysicalLineThicknessF = 6.5
  mpres@mpNationalLineThicknessF    = 6.5
  mpres@mpUSStateLineThicknessF     = 6.5
  mpOutlineBoundarySets             = "National" 
  mpDataBaseVersion                 = "MediumRes"
  mpFillOn                          = True 
  mpProjection                      = "Mercator"
  resMap               = True  
  resMap@mpDataSetName = "Earth..3"
  

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