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

begin
  dir      = "./"
  filename = dir + "wrfout_d01_2007-03-03_00:00:00"
  filename = filename + ".nc"
  a = addfile(filename,"r")

  hgt   = wrf_user_getvar(a, "HGT", 0)
  lat2d = wrf_user_getvar(a, "XLAT", 0)
  lon2d = wrf_user_getvar(a, "XLONG", 0)
  hgt@lat2d = lat2d
  hgt@lon2d = lon2d

  printVarSummary(hgt)
  printMinMax(hgt@lat2d,0)
  printMinMax(hgt@lon2d,0)

  wks = gsn_open_wks("x11","test")

  res = True
  res = wrf_map_resources(a,res)

  res@cnFillOn    = True
  res@cnLinesOn   = False
  res@gsnMaximize = True

;---Change some of the defaults that wrf_map_resources used.
  res@mpGeophysicalLineColor      = "Black"
  res@mpNationalLineColor         = "Black"
  res@mpUSStateLineColor          = "Black"
  res@mpGeophysicalLineThicknessF = 1.5
  res@mpNationalLineThicknessF    = 1.5
  res@mpUSStateLineThicknessF     = 1.5

  plot = gsn_csm_contour_map(wks,hgt,res)

;---Pick arbitrary area to zoom in on.
  y_start = 100
  x_start = 100
  y_end   = 175
  x_end   = 175

  hgt_dom = hgt(y_start:y_end,x_start:x_end)
  lat_dom = lat2d(y_start:y_end,x_start:x_end)
  lon_dom = lon2d(y_start:y_end,x_start:x_end)
  nlat = dimsizes(lat_dom(:,0))
  nlon = dimsizes(lon_dom(0,:))
  printMinMax(lat_dom,0)
  printMinMax(lon_dom,0)
;
; Draw the plot again with the corners of the new
; domain plotted as filled dots.
;
  lat_pts = (/lat_dom(0,0),lat_dom(nlat-1,nlon-1)/)
  lon_pts = (/lon_dom(0,0),lon_dom(nlat-1,nlon-1)/)

  mkres               = True
  mkres@gsMarkerIndex = 16
  mkres@gsMarkerColor = "yellow"
  mkres@gsMarkerSizeF = 20
  id = gsn_add_polymarker(wks,plot,lon_pts,lat_pts,mkres)
  draw(plot)
  frame(wks)

;---Create a new plot of the new zoomed in domain area.
  hgt_dom@lat2d := lat_dom
  hgt_dom@lon2d := lon_dom
  printVarSummary(hgt_dom)

  res@mpLeftCornerLonF  = lon_dom(0,0)
  res@mpLeftCornerLatF  = lat_dom(0,0)
  res@mpRightCornerLonF = lon_dom(nlat-1,nlon-1)
  res@mpRightCornerLatF = lat_dom(nlat-1,nlon-1)

  plot = gsn_csm_contour_map(wks,hgt_dom,res)
end