;   Example script to plot some 3D fields from a single metgrid file
;   Note the change in map and wind vector colors and the addition of PlotLevel info
;   November 2008

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

begin
;
  f = systemfunc("ls -1 /home/bchaimite/lustre/WRFChem/WRFV3/run/wrfout_oldfiles/wrfout_d01*")  ; Open a file
  a = addfiles(f,"r")
  ListSetType(a,"cat")
  
  f1 = systemfunc("ls -1 /mnt/lustre/users/bchaimite/WRFCHEMBB0/WRFChem/WRFV3/run/wrfout_d01*")
  a1 = addfiles(f1,"r")
  ListSetType(a1,"cat")

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


  res = True                           ; Set up some basic plot resources
  res@MainTitle = ""
  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"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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


    t = wrf_user_getvar(a,"tc",-1)  ; Calculate tc from TT 
    tc = dim_avg_n_Wrap(t,0)
    ;u = wrf_user_getvar(a,"ua",-1)  ; Get U on mass points
    ;ua = dim_avg_n_Wrap(u,0)
    ;v = wrf_user_getvar(a,"va",-1)  ; Get V on mass points
    ;va = dim_avg_n_Wrap(v,0)

    t1 = wrf_user_getvar(a1,"tc",-1)
    tc1= dim_avg_n_Wrap(t1,0)
   
    tc_diff = tc -tc1
    tc_diff@description = "Temperature Effect"
    printVarSummary(tc_diff)
    printMinMax(tc_diff,False)

     p  = wrf_user_getvar(a,"z",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

    ;dims = dim_avg_n_Wrap(tc,0)
    ;nd = dimsizes(dims)
    ;nl = dims
    
    do level = 1,nl-1,5      ; LOOP OVER LEVELS, plot every 5th (surface is 0, lets start one level up)
 
       res@PlotLevelID = p(level,0,0)*.001 + "km"    ; Add level info to plot
 
     ; Temperature
       opts = res
       opts@cnFillOn = True
       opts@ContourParameters = (/ 0.5 /)
       opts@lbBoxLinesOn = True           ; Move lines around label bar
       contour = wrf_contour(a[0],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[0],wks,ua(level,:,:),va(level,:,:),opts)
       ;delete(opts)

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

    end do      ; END OF LEVEL LOOP

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

end