;*************************************************
; WRF: Panel an individual variable at 6 hour time steps 
;************************************************
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/WRF_contributed.ncl"
begin
;************************************************
; get all desired files
;************************************************
  diri = "./"
  fili = systemfunc("cd "+diri +"; ls *2003*")
  fili = fili +".nc"                           ; explicitly add file extension

  F    = addfiles (diri+fili, "r")             ; refers to *all* files      
  f    = addfile  (diri+fili(0), "r")          ; refers to 1st file only
;************************************************
; open file and read in data
; (1) Read Q2 at all times: convert to g/kg
; (2) Read character variable Times; Convert to string for plots
;************************************************
 
  x       = addfiles_GetVar (F, diri+fili,"Q2"); read Q@ across all files
  x       = x*1000.                            ; convert units
  x@units = "g/kg"                             ; upgrade attribute

  times   = chartostring(F[:]->Times)          ; built-in function
  ntim    = dimsizes(times)                    ; # of time steps
;************************************************
; create plots  
;************************************************
  wks = gsn_open_wks("ps" ,"WRF_me")           ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

  res                       = True             ; plot mods desired
  res@gsnSpreadColors       = True             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels

;************************************************
; Use WRF_contributed procedure to set map resources
;************************************************
  WRF_map_c(f, res, 0)                   ; reads info from file

;************************************************
; if appropriate, set True for native (direct) mapping (faster)
; set False otherwise
;************************************************
  res@tfDoNDCOverlay       = True 

;************************************************
; associate the 2-dimensional coordinates to the variable for plotting
; only if non-native plot
;************************************************
  if (.not.res@tfDoNDCOverlay) then
      x@lat2d = f->XLAT(0,:,:)                      ; direct assignment
      x@lon2d = f->XLONG(0,:,:)
  end if

;************************************************
; For individual plots 
;************************************************
  res@gsnDraw               = False            ; do not draw 
  res@gsnFrame              = False            ; do not advance 'frame'
  res@lbLabelBarOn          = False            ; turn off individual lb's 
;************************************************
; common label bar for panel requires that
; all plots should be set to the same interval
; Use built in function to determine "nice" limits
;************************************************
  mnmxint                   = nice_mnmxintvl( min(x), max(x), 14, False)
  res@cnLevelSelectionMode  = "ManualLevels"   
  res@cnMinLevelValF        = mnmxint(0)
  res@cnMaxLevelValF        = mnmxint(1)
  res@cnLevelSpacingF       = mnmxint(2)
;************************************************
; Allocate array to store plots: specify time step
;************************************************
  tstep                     = 3                ; time step to plot          
  plts = new (ntim/tstep , "graphic")          ; 1d array to hold plots

;************************************************
; loop over each forecast 
;************************************************
  n    = -1                                    ; counter
  do nt=1,ntim-1,tstep                         ; plot every 6 hours 
     res@gsnLeftString         = times(nt)
     n  = n+1
     plts(n)                   = gsn_csm_contour_map(wks,x(nt,:,:),res)
  end do
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@txString         = x@description 
  resP@gsnPanelLabelBar = True                   ; add common colorbar
  resP@gsnMaximize      = True                   ; maximize panel area
  resP@lbLabelAutoStride= True                   ; let NCL figure lb stride
  gsn_panel(wks,plts,(/2,2/),resP)               ; now draw as one plot
end