;   Example script to plot some 3D fields from a number of metgrid input files
;   Plot some 3D fields - over a number of input files
;   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_7")


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

  pltres = True
  mpres = True
  mpres@mpGeophysicalLineColor      = "Black"   ; Overwrite map settings
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpPerimLineColor            = "Black"
  mpres@mpUSStateLineColor          = "Black"


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Loop over all input files

  do ifil = 0,numFILES-1

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        
; Note: we do not have "tc", "ua", "va" in the input field - but we
;       know how to calculate them

    tc = wrf_user_getvar(a,"tc",0)  ; Calculate tc from TT 
    ua = wrf_user_getvar(a,"ua",0)  ; Get U on mass points
    va = wrf_user_getvar(a,"va",0)  ; Get V on mass points

    p  = wrf_user_getvar(a,"PRES",0) 

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    dims = dimsizes(tc)          ; Get dims for tc
    nd = dimsizes(dims)          ; Find out how many dimensions there are
    nl = dims(nd-3)              ; We know 3rd dim from right is nz

    i = 0
    do level = 1,nl-1,5      ; LOOP OVER LEVELS, plot every 5th (surface is 0, lets start one level up)
       i = i + 1
 
       res@PlotLevelID = p(level,0,0)*0.01 + "hPa"    ; Add level info to plot

 
     ; Temperature
       opts = res
     ; Always ensure levels span same range
       if ( i .eq. 1 ) then
         opts@ContourParameters = (/ -22., 24., .5 /)
       end if
       if ( i .eq. 2 ) then
         opts@ContourParameters = (/ -22., 24., .5 /)
       end if
       if ( i .eq. 3 ) then
         opts@ContourParameters = (/ -24., 4., .5 /)
       end if
       if ( i .eq. 4 ) then
         opts@ContourParameters = (/ -42.,-20., .5 /)
       end if
       if ( i .eq. 5 ) then
         opts@ContourParameters = (/ -66.,-44., .5 /)
       end if
       if ( i .eq. 6 ) then
         opts@ContourParameters = (/ -68.,-60., .5 /)
       end if
       opts@cnFillOn = True
       opts@lbBoxLinesOn = False
       contour = wrf_contour(a,wks,tc(level,:,:),opts)
       delete(opts)
    
     ; Wind 
       opts = res
       opts@FieldTitle = "Wind"            ; Overwrite Field Title
       opts@NumVectors = 35                ; wind barb density - higher is more vectors
       opts@vcWindBarbColor = "White"      ; Draw wins barbs in white
       vector = wrf_vector(a,wks,ua(level,:,:),va(level,:,:),opts)
       delete(opts)

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

    end do      ; END OF LEVEL LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do      ; END OF INPUT FILES

end