;*************************************************
; WRF: plot data with "missing_values"
;************************************************
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 
; (1) Read SST and SMOIS at specified time/level
; (2) Read character variable Times; Convert to string for plots
;************************************************
  f       = addfile ("wrfout_d01_000000.nc", "r")
  sm      = f->SMOIS(12,0,:,:)      ; (Time, soil_layers_stag, south_north, west_east )
  sm@_FillValue = 1.0               ; manually set _FillValue 
  sst     = f->SST(12,:,:)          ; (Time, south_north, west_east )
  sst@_FillValue = 0.0              ; manually set _FillValue

  times   = chartostring(f->Times)  ; convert to type string
  ntim    = dimsizes(times)         ; # time steps

;************************************************
; create plots  
;************************************************
  wks = gsn_open_wks("ps"  ,"WRF_lc")          ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

  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
  res@lbOrientation         = "Vertical"       ; default is horizontal 
;************************************************
; 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)
;************************************************
  res@tfDoNDCOverlay       = True      

;************************************************
; associate the 2-dimensional coordinates to the variable for plotting
;************************************************
  if (.not.res@tfDoNDCOverlay) then
      lat2d     = f->XLAT(0,:,:)              
      lon2d     = f->XLONG(0,:,:)      
    
      sm@lat2d  = lat2d
      sm@lon2d  = lon2d
      sst@lat2d = lat2d
      sst@lon2d = lon2d
  end if

;************************************************
; For individual plots 
;************************************************
  res@gsnDraw               = False            ; do not draw 
  res@gsnFrame              = False            ; do not advance 'frame'
;************************************************
; Allocate array to store plots: specify time step
;************************************************
  plts = new (2 , "graphic")                   ; 1d array to hold plots

;************************************************
; NCL contouring does not like _FillValue=0.0
;************************************************
  sst@_FillValue = -999.                       ; change to -999

;************************************************
; create panel: panel plots have their own set of resources
; for demo: use Raster fo soil moisture 
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@gsnMaximize      = True                   ; maximize panel area

  res@cnFillMode     = "RasterFill"              ; turn raster on         
  plts(0) = gsn_csm_contour_map(wks,sm,res)
  delete(res@cnFillMode)                     ; delete raster mode
  plts(1) = gsn_csm_contour_map(wks,sst,res)

  resP@txString      = times(12)      
  gsn_panel(wks,plts,(/2,1/),resP)               ; now draw as one plot
end