begin

a = addfile("wrfout_d03_2007-09-24_05-07:00:00.nc","r")

; type = "x11"
  type = "pdf"
; type = "ps"
; type = "ncgm"

wks = gsn_open_wks(type,"d04_AVO_loop_05-07")

  pltres                             = True
  pltres@NoTitles                    = True 
  mpres                              = True
  mpres@mpOutlineOn                  = True
  mpres@gsnMaximize                  = True
  mpres@mpFillOn                     = False
  mpres@mpGeophysicalLineThicknessF  = 2.0
  mpres@mpGeophysicalLineColor       = "Black"

;***********************************************
; plot controls
;***********************************************

  res          = True
  res@gsnDraw  = False
  res@gsnFrame = False 
   
;  cmap                     = read_colormap_file("vorticity12")
;  cmap(2:8,:)              = 0.0
;  cmap(0:4,3)              = 0.3
;  cmap(7:8,3)              = 0.3

;  res@cnFillColors         = cmap
;  res@cnFillOn             = True               ; turns on the color
;  res@cnLinesOn            = False              ; turn off contour lines

  res@cnFillOn             = False
  res@cnLinesOn            = True

  res@cnLevelSelectionMode = "ManualLevels"
  res@cnMinLevelValF       = 0.005
  res@cnMaxLevelValF       = 0.040
  res@cnLevelSpacingF      = 0.005
  res@cnLineThicknessF     = 0.8
  res@cnInfoLabelOn        = False    

;  res@lbTitleOn            =  True                          ; turn on title
;  res@lbTitleString        = "Absolute Vorticity (1/s)"     ; title string
;  res@pmLabelBarOrthogonalPosF = -0.05                      ; distances label bar

xlat  = wrf_user_getvar(a, "XLAT",0)     
xlon  = wrf_user_getvar(a, "XLONG",0)      
lat   = xlat(:,0)                                
lon   = xlon(0,:)        

;**********************************************
; Begin loop
;**********************************************

times = wrf_user_list_times(a)     ; get times in the file
ntimes = dimsizes(times)           ; number of times in the file
contour = new(ntimes,graphic)

do it = 0,ntimes-1                 ; loop

print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots

avo  = wrf_user_getvar(a,"avo",it)
z     = wrf_user_getvar(a,"z",it)

height = 500.
avo_plane = wrf_user_intrp3d(avo,z,"h",height,0.,False)    ; Interpolates avo to height
avo_plane1 = avo_plane/100000

printVarSummary(avo_plane1)
printMinMax(avo_plane1,0)

contour(it) = wrf_contour(a,wks,avo_plane1,res)

end do

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

end