;----------------------------------------------------------------------
; Procedure that adds lat/lon lines to given NCL plot.
;----------------------------------------------------------------------
procedure add_latlon_lines(wks,plot,lat[*][*],lon[*][*],color)
local lnres
begin
  lnres                 = True
  lnres@gsLineColor      = color
  lnres@gsLineThicknessF = 2.0
  lnres@gsnCoordsLat     = lat
  lnres@gsnCoordsLon     = lon
  lnres@gsnCoordsAsLines = True    ; draw grid as lines, not markers
  lnres@gsnCoordsAttach  = True    ; lines will be attached, and not drawn (yet)

  x = 1   ; dummy variable
  gsn_coordinates(wks,plot,x,lnres)
end


;----------------------------------------------------------------------
; Procedure that adds a marker at a lat, lon location on the given
; NCL plot.
;----------------------------------------------------------------------
procedure add_marker(wks,plot,lat,lon,color,size)
local mkres
begin
  mkres = True
  mkres@gsMarkerIndex = 16    ; filled dot
  mkres@gsMarkerColor = color
  mkres@gsMarkerSizeF = size
  tmpstr = unique_string("marker")
  plot@$tmpstr$ = gsn_add_polymarker(wks,plot,lon,lat,mkres)
end


;----------------------------------------------------------------------
; Main driver code
;----------------------------------------------------------------------
begin
;---Open WRF output file
  filename = "wrfout_d01_2005-12-14_13:00:00"
  a        = addfile(filename,"r")
  hgt = wrf_user_getvar(a,"HGT",0)
  lat = wrf_user_getvar(a,"XLAT",0)
  lon = wrf_user_getvar(a,"XLONG",0)

  wks = gsn_open_wks("png","wrf_gsn")

;---Resources for filled contour plot
  res                      = True

  res@gsnMaximize          = True
  res@gsnDraw              = False
  res@gsnFrame             = False

  res@cnFillOn             = True  
  res@cnLinesOn            = False
  res@cnFillPalette        = "OceanLakeLandSnow"
  res@cnFillOpacityF       = 0.5
 
;---Necessary for contours to be overlaid correctly on WRF projection when using native projection
  res = wrf_map_resources(a,res)
  res@tfDoNDCOverlay       = True
  res@gsnAddCyclic         = False    ; don't add longitude cyclic point

;---Set map resources based on projection on WRF output file
  res = wrf_map_resources(a,res)

;---Create contours over map (won't get drawn yet)
  plot = gsn_csm_contour_map(wks,hgt,res)

;---Get lon, lat at particular i,j location
  i = 27                                  ; indexes go from 1 to n
  j = 27
  loc  = wrf_user_ij_to_ll (a,i,j,res)    ; lon,lat returned
  lat1 = lat(j-1,i-1)                     ; indexes go from 0 to n-1
  lon1 = lon(j-1,i-1)
  print("lon,lat   = " + loc(0) + "," + loc(1))
  print("lon1,lat1 = " + lon1 + "," + lat1)

  add_latlon_lines(wks,plot,lat,lon,"gray35")
  add_marker(wks,plot,lat1,lon1,"navyblue",20)
  add_marker(wks,plot,loc(1),loc(0),"red",10)

  draw(plot)
  frame(wks)
  delete(res)

;----------------------------------------------------------------------
; Draw a new plot that is zoomed in on the area of interest. Won't
; use the native WRF map projection in this plot.
;----------------------------------------------------------------------

  hgt@lat2d = lat    ; allows us to zoom in easily
  hgt@lon2d = lon

  res                   = True

  res@gsnMaximize       = True
  res@gsnDraw           = False
  res@gsnFrame          = False

  res@cnFillOn          = True  
  res@cnLinesOn         = False
  res@cnFillPalette     = "OceanLakeLandSnow"
  res@cnFillOpacityF    = 0.5
  res@gsnAddCyclic      = False    ; don't add longitude cyclic point

;---Zoom in on plot
  res@mpMinLatF     = 32
  res@mpMaxLatF     = 37
  res@mpMinLonF     = -95
  res@mpMaxLonF     = -87
  res@mpDataBaseVersion = "MediumRes"

  res@pmTickMarkDisplayMode = "Always"  ; better map tickmarks

  plot_zoom = gsn_csm_contour_map(wks,hgt,res)

  add_latlon_lines(wks,plot_zoom,lat,lon,"gray35")
  add_marker(wks,plot_zoom,lat1,lon1,"navyblue",20)
  add_marker(wks,plot_zoom,loc(1),loc(0),"red",10)

  draw(plot_zoom)
  frame(wks)

end