;----------------------------------------------------------------------
; This is a basic NCL template for creating contours over maps using
; WRF data. This template plots just one timestep from a WRF file.
;
; This script does NOT use the wrf_xxxx scripts to do the plotting.
; It shows how to use gsn_csm_xxxx scripts to do the plotting.
;----------------------------------------------------------------------

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

begin
;---Open file. You may need to include ".nc" at the end.
  a = addfile("wrfout_d01_2005-12-14_13:00:00.GWATC_FCST.nc","r")

;---Read variables directly or use wrf_user_getvar
  hgt = a->HGT(0,:,:)       ; 0 is the first time step
  lat = a->XLAT(0,:,:)
  lon = a->XLONG(0,:,:)

;---Open workstation and change color map
  wks = gsn_open_wks("x11","wrf_contour_map")   ; "ps", "pdf", "png"
  gsn_define_colormap(wks,"OceanLakeLandSnow")

;---Set some resources
  res                     = True

  res@gsnMaximize         = True     ; maximize plot in frame

  res@cnFillOn            = True     ; turn on contour fill

  res@cnLinesOn           = False    ; turn off contour lines
  res@cnLineLabelsOn      = False    ; turn off line labels

;---These values will likely need to be changed
  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF       = 1
  res@cnMaxLevelValF       = 2000
  res@cnLevelSpacingF      = 40

;---These two resources are not needed in V6.1.0
  res@gsnSpreadColors      = True     ; span full color map
  res@lbLabelAutoStride    = True     ; nice spacing for labelbar labels

  res@tiMainString         = "Variable 'HGT'"

  res@gsnAddCyclic      = False  ; Set to False if plotting regional data

;
; This function looks at the global attributes on the WRF output 
; file and sets some map resources based on their values.
; 
  res = wrf_map_resources(a,res)
; print(res)     ; If you want to see what resources were set.

;---Necessary to put data on map correctly.
  res@sfXArray          = lon
  res@sfYArray          = lat

  plot = gsn_csm_contour_map(wks,hgt,res)
end