;   Example script to plot a field from metgrid files
;   Loop over all available metgrid files and plot sea level pressure from all files
;   Change input field from Pa to hPa and ensure all plots span same values
;   November 2008

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

begin

; Make a list of all files we are interested in
  DATADir = "./"
  FILES = systemfunc (" ls -1 " + DATADir + "met_em.d01* ")
  numFILES = dimsizes(FILES)


; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_metgrid_4")


  res = True                           ; Set up some basic plot resources
  res@MainTitle = "METGRID FILES"
  res@Footer = False

  pltres = True
  mpres = True


  opts = res
  opts@cnFillOn = True
  opts@ContourParameters = (/ 985., 1025., 5. /)  ; ensure plot levels are always the same

  do ifil = 0,numFILES-1

    a = addfile(FILES(ifil)+".nc","r")

    slp = wrf_user_getvar(a,"PMSL",0)    ; Get seal level pressure from file

    slp = slp * 0.01                     ; We want slp in Pa
    slp@units = "hPa"                    ; Make sure the units match

    contour = wrf_contour(a,wks,slp,opts)
    plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

  end do
   
end