;----------------------------------------------------------------------
; NCL script to plot fields from a WRF simulation
;----------------------------------------------------------------------
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"
;----------------------------------------------------------------------
begin
;---Open WRF output file
;  diri = "/home/sfermepin/WRF/SNCF/FR2/RUNS/mpWSM3-raRRTM-blYSU-cuKF/20141023_00/"
; fname = "wrfout_d01_2014-10-23_00:00:00.nc"
  fname = "t2_wrfout_d01_2014-10-23_00:00:00.nc"
  diri = "./"
  frun = diri + fname
  a        = addfile(frun,"r")
  date="20141023"
  dom="03"

;---Read WRF variable
  it   = 0      ; first record
  tc2  = wrf_user_getvar(a,"T2",0);
  tc2 = tc2-273.15

;---Change the metadata
  tc2@description  = "2m Temperature"
  tc2@units        = "degC"

  wks = gsn_open_wks("png","wrf_test")

;---Set common resources for all plots
  res                          = True
  res@gsnMaximize              = True
  res@gsnDraw                  = False
  res@gsnFrame                 = False

  res@cnFillOn                 = True
  res@cnLinesOn                = False
  res@cnLineLabelsOn           = False
  res@cnInfoLabelOn            = False
  res@cnLevelSelectionMode     = "ExplicitLevels"
  res@cnLevels                 =  (/ 0.,2.,4,5,6,7,8,9,10,11,12,13,14,15,20/)

  res@lbLabelBarOn             = False     ; turn labelbar off and draw in panel

  res@gsnLeftString            = "Temperature (" + tc2@units + ")"
  res@gsnRightString           = ""

  res                          = wrf_map_resources(a,res)
  res@gsnAddCyclic             = False
  res@tfDoNDCOverlay           = True

  res@mpGeophysicalLineColor        = "black"    ; So we can see the map
  res@mpGeophysicalLineThicknessF   = 2.0        ; outlines better.

;---Plot using "MediumRes" map database
  res@tiMainString = "MediumRes NCL map database"
  plot_med = gsn_csm_contour_map(wks,tc2,res)

;---Plot using "HighRes" map database
  res@tiMainString      = "HighRes RANGS map database"
  res@mpDataBaseVersion = "HighRes"
  plot_hgh = gsn_csm_contour_map(wks,tc2,res)

;---Plot using shapefile outlines, downloaded from gadm.org/country.
  res@tiMainString = "Shapefile outlines"
  res@mpOutlineOn  = False                             ; turn off NCL outlines so we can add shapefile outlines
  plot_shp = gsn_csm_contour_map(wks,tc2,res)
  id_fra   = gsn_add_shapefile_polylines(wks,plot_shp,"FRA_adm/FRA_adm0.shp",True)    ; France
  id_esp   = gsn_add_shapefile_polylines(wks,plot_shp,"ESP_adm/ESP_adm0.shp",True)    ; Spain

;---Draw all three plots full-sized.
  draw(plot_med)
  frame(wks)
  draw(plot_hgh)
  frame(wks)
  draw(plot_shp)
  frame(wks)

;----------------------------------------------------------------------
; Plot the WRF lat/lon grid for debugging purposes.  This helps
; us see whether the map projection on the file is correct.
;----------------------------------------------------------------------
  tc2@lat2d = a->XLAT(0,:,:)
  tc2@lon2d = a->XLONG(0,:,:)

  gres                  = True
  gres@gsnCoordsAsLines = True
  gres@gsLineColor      = "black"
  gsn_coordinates(wks,plot_med,tc2,gres)

;----------------------------------------------------------------------
; Panel all three plots so we can compare them.
;----------------------------------------------------------------------
  pres                  = True
  pres@gsnPanelLabelBar = True
  pres@gsnPanelRowSpec  = True
  gsn_panel(wks,(/plot_med,plot_hgh,plot_shp/),(/2,1/),pres)
end