;   Example script - plot terrain (using some basic WRF scripts)
;   The map background is added via the special wrf_map_overlays function
;   We get the HGT_M field using the special function wrf_user_getvar
;   November 2008

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

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

begin

  a = addfile("./geo_em.d04.nc","r")           ; Open a file


  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_geo_5")         ; Create a plot workstation
  

  opts = True                                  ; Set some Basic Plot options
  opts@MainTitle = "GEOGRID FIELDS"

  ter = wrf_user_getvar(a,"HGT_M",0)           ; Read the variable to memory

  res = opts                                   ; Use basic options for this field
  res@cnFillOn = True                          ; Create a color fill plot
  res@ContourParameters = (/ 0., 1000., 50. /) ; Set the levels

  contour = wrf_contour(a,wks,ter,res)

  pltres = True                                ; Set plot options
  mpres = True                                 ; Set map options
  plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres) ; Plot the field over a map background
 
  pmres = True
  pmres@gsMarkerColor = "Black"
  pmres@gsMarkerIndex = 16
  pmres@gsMarkerSizeF = 0.01

  

  gsn_polymarker(wks,contour,2.53115,36.8569,pmres)
   draw(plot)
  frame(wks)

end