;*************************************************
; WRF: panel three different variables at the same time step
;************************************************
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/wrf/WRF_contributed.ncl"
begin
;************************************************
; open file and read in data
;************************************************
  f     = addfile("wrfout_d01_2003-07-13_12:00:00.nc", "r")
;************************************************
; Read character variable Times; Convert to string for plots
; Read vertical coordinate for plot labels
;************************************************
  times = chartostring(f->Times)               ; built-in function
  znu   = f->ZNU(0,:)                          ; (Time, bottom_top)
;************************************************
; Read perturbation geopotential, pressure, water vapor 
; at all times and levels
;************************************************
  ph    = f->PH                     ; (Time, bottom_top, south_north, west_east)
  p     = f->P                      ; (Time, bottom_top, south_north, west_east)
  q2    = f->Q2                     ; (Time, south_north, west_east )
;************************************************
; create plots  
;************************************************
  wks = gsn_open_wks("ps" ,"WRF_me")           ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks ,"BlAqGrYeOrReVi200"); choose colormap

  res                       = True             ; plot mods desired
;;res@gsnMaximize           = True             ; uncomment to maximize size
  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
  res@lbLabelAutoStride     = True             ; let NCL figure lb stride

;************************************************
; Use WRF_contributed procedure to set map resources
;************************************************
  WRF_map_c(f, res, 0)                   ; reads info from file
;************************************************
; set True for native direct mapping (faster)
;************************************************
  res@tfDoNDCOverlay       = true

;************************************************
; associate the 2-dimensional coordinates to variables for plotting
; only if res@tfDoNDCOverlay=False
;************************************************
  if (.not.res@tfDoNDCOverlay) then
      lat2d    = f->XLAT(0,:,:)                 ; need for map limits        
      lon2d    = f->XLONG(0,:,:)      
    
      p@lat2d  = lat2d
      p@lon2d  = lon2d
      ph@lat2d = lat2d
      ph@lon2d = lon2d
      q2@lat2d = lat2d
      q2@lon2d = lon2d
  end if

;************************************************
; allocate array for 3 plots
;************************************************
  plts                      = new (3,"graphic")   
;************************************************
; Specify (arbitrarily chosen) subscripts 
; This could also be done in a do loop or explicitly specified
;************************************************
  nt                        = 3                ; last time step
  kl                        = 3      
;************************************************
; Tell NCL not to draw or advance frame for individual plots
;************************************************
  res@gsnDraw               = False            ; (a) do not draw 
  res@gsnFrame              = False            ; (b) do not advance 'frame'

  res@gsnLeftString         = ph@description+": znu="+znu(kl)
  plts(0)                   = gsn_csm_contour_map(wks,ph(nt,kl,:,:),res)

  res@gsnLeftString         = p@description+": znu="+znu(kl)
  plts(1)                   = gsn_csm_contour_map(wks,p(nt,kl,:,:),res)
  delete(res@gsnLeftString)            

  plts(2)                   = gsn_csm_contour_map(wks,q2(nt,:,:),res)
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@txString         = f@TITLE+": "+times(nt)      
  resP@gsnMaximize      = True                   ; maximize panel area
  resP@gsnPanelRowSpec  = True                   ; specify 1 top, 2 lower level  
  gsn_panel(wks,plts,(/1,2/),resP)               ; now draw as one plot

end