;----------------------------------------------------------------------
; This script reads a mask array off a file that was created by 
; running another NCL script called create_wrf_mask.ncl.
;
; By reading the mask array off the file, rather than generating
; it every time via "shapefile_mask_data", you save yourself a
; lot of time!
;----------------------------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
; 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/WRFUserARW.ncl"

load "./shapefile_utils.ncl"

begin
;---Open WRF output file
  filename = "wrfout_d01_2005-12-14_13:00:00"
  a        = addfile(filename,"r")

;---Read several WRF variables at first time step
  it    = 0
  tc    = wrf_user_getvar(a,"tc",it)      ; 3D temperature
  slp   = wrf_user_getvar(a,"slp",it)     ; 2D sea level pressure
  ter   = wrf_user_getvar(a,"ter",it)     ; 2D terrain

  lat2d = wrf_user_getvar(a,"lat",it)     ; 2D latitude
  lon2d = wrf_user_getvar(a,"lon",it)     ; 2D longitude

;---Get the lowest (bottommost) level for temperature
  nl  = 0
  tc2 = tc(nl,:,:)

;---This is necessary for plotting only!
  tc2@lat2d = lat2d    
  tc2@lon2d = lon2d
  slp@lat2d = lat2d    
  slp@lon2d = lon2d
  ter@lat2d = lat2d    
  ter@lon2d = lon2d

;---Open file with WRF mask and apply to variables
  f        = addfile("wrf_mask.nc","r")
  marray   = f->WRF_mask
  tc2_mask = mask(tc2,marray,1)
  slp_mask = mask(slp,marray,1)
  ter_mask = mask(ter,marray,1)

;---Copy metadata from original variables to masked variables
  copy_VarMeta(tc2,tc2_mask)  
  copy_VarMeta(slp,slp_mask)
  copy_VarMeta(ter,ter_mask)

;---Plot original data and masked data.
  wks = gsn_open_wks("png","wrf_mask")

  res                = True
  res@gsnMaximize    = True
  res@gsnDraw        = False  ; comment these if you want to see the
  res@gsnFrame       = False  ; individual plots before the panel
  res@cnFillOn       = True  
  res@lbOrientation  = "Vertical"
  res = wrf_map_resources(a,res)    ; use same map projection as defined on WRF file.

;---Plot original data with different color maps
  res@cnFillPalette = "BlueYellowRed"
  plot_orig_tc2 = gsn_csm_contour_map(wks,tc2,res)

  res@cnFillPalette = "BlAqGrYeOrReVi200"
  plot_orig_slp = gsn_csm_contour_map(wks,slp,res)

  res@cnFillPalette = "OceanLakeLandSnow"
  plot_orig_ter = gsn_csm_contour_map(wks,ter,res)

;---Plot masked data with the same color maps as their original data.
  res@cnFillPalette = "BlueYellowRed"
  plot_mask_tc2 = gsn_csm_contour_map(wks,tc2_mask,res)
  
  res@cnFillPalette = "BlAqGrYeOrReVi200"
  plot_mask_slp = gsn_csm_contour_map(wks,slp_mask,res)

  res@cnFillPalette = "OceanLakeLandSnow"
  plot_mask_ter = gsn_csm_contour_map(wks,ter_mask,res)

  pres = True
  pres@gsnPanelMainString = "WRF data masked by Mississippi River Basin shapefile"
  pres@gsnMaximize = True
  gsn_panel(wks,(/plot_orig_tc2,plot_orig_slp,plot_orig_ter,\
                  plot_mask_tc2,plot_mask_slp,plot_mask_ter/),(/2,3/),pres)
end