; Basic plot of WRF temperature overlaid with wind vectors

begin
;---Open WRF output file
  filename = "/home/rbehnke/WRF/WRFV3/test/em_real/wrfout_d02_2017-06-01_03:00:00"
  a        = addfile(filename,"r")

;---Read several WRF variables at first time step
  it  = -1
  slp_1 = wrf_user_getvar(a,"PSFC",it)   ; surface pressure
  slp_1 = slp_1/100.0
  wrf_smooth_2d(slp_1, 3)             ; smooth slp
  t2_1  = wrf_user_getvar(a,"T2",it)    ; surface temperature
  t2_1  = t2_1 - 273.15	
  u10_1   = wrf_user_getvar(a,"U10",it)    ; surface U
  v10_1   = wrf_user_getvar(a,"V10",it)    ; surface V
  pblh_1  = wrf_user_getvar(a,"PBLH",it)   ; planetary bdry layer height	

  slp = slp_1(0,:,:)
  t2 = t2_1(0,:,:) 
  u10 = u10_1(0,:,:)
  v10 = v10_1(0,:,:)
  pblh = pblh_1(0,:,:)

;---Change the metadata
  t2@description = "Surface Temperature"
  t2@units       = "deg C"
  u10@units      = "m/s"
  v10@units      = "m/s"
  pblh@units     = "m"
  slp@units      = "mb"

  wks = gsn_open_wks("png","wrfout_d02_2017-06-01_03Z")

;---Set common resources for all plots
  res                = True
  res@gsnFrame       = False
  res@gsnDraw        = False
  res@gsnLeftString   = ""
  res@gsnRightString  = ""
  res@gsnCenterString = ""

;---Necessary for contours to be overlaid correctly on WRF projection
  res@tfDoNDCOverlay = True
   
;---Temperature filled contour plot
  t2_res                      = res
  t2_res@cnFillOn             = True  
  t2_res@cnLevelSelectionMode = "ExplicitLevels"
  t2_res@cnLevels             = ispan(8,15,1)
  t2_res@lbLabelFontHeightF   = 0.015
  t2_res@lbOrientation        = "Vertical"
  t2_res@pmLabelBarOrthogonalPosF = -0.005

  contour_tf = gsn_csm_contour(wks,t2,t2_res)

;---SLP line contour plot
  levels      = ispan(700,1000,50)
  info_string = "Sea level pressure contours from 700 to 1000 by 50"

  slp_res                                = res
  slp_res@cnLineColor                    = "NavyBlue"
  slp_res@cnLevelSelectionMode           = "ExplicitLevels"
  slp_res@cnLevels                       = levels
  slp_res@cnLineLabelBackgroundColor     = -1    ; transparent
  slp_res@cnLineThicknessF               = 2.5
  slp_res@cnHighLabelsOn                 = True
  slp_res@cnLowLabelsOn                  = True
  slp_res@cnHighLabelBackgroundColor     = -1
  slp_res@cnLowLabelBackgroundColor      = -1
  slp_res@cnInfoLabelString              = info_string
  slp_res@cnInfoLabelFontColor           = "NavyBlue"
  slp_res@cnInfoLabelPerimOn             = False

; contour_psl = gsn_csm_contour(wks,slp,slp_res)

;---Wind vector plot
  vec_res                  = res
  vec_res@vcMinDistanceF   = 0.02
  vec_res@vcRefLengthF     = 0.06 
  vec_res@vcMinFracLengthF = 0.25
  vec_res@vcGlyphStyle     = "CurlyVector"
  vec_res@vcRefAnnoOn      = False
  vec_res@vcLineArrowColor        = "black"           ; change vector color
  vec_res@vcLineArrowThicknessF   = 1.5               ; change vector thickness

  vector = gsn_csm_vector(wks,u10,v10,vec_res)
  
;---Map plot
  map_res               = True
  map_res@gsnFrame      = False
  map_res@gsnDraw       = False
  map_res@tiMainString  = "Surface Temperature and Wind"
  map_res@gsnLeftString = t2@description + " (" + t2@units + ")~C~" + \
                          ; slp@description + " (" + slp@units + ")~C~" + \
                          "Wind (" + u10@units + ")"
  map_res@gsnLeftStringFontHeightF = 0.01

;---Set map resources based on projection on WRF output file
  map_res = wrf_map_resources(a,map_res)
  map = gsn_csm_map(wks,map_res)

;---Overlay plots on map and draw.
  overlay(map,contour_tf)
  ; overlay(map,contour_psl)
  overlay(map,vector)

;---Attach the shapefile polylines using files read off gadm.org/country.
  states_shp = "/home/rbehnke/Desktop/cb_2015_us_state_500k/cb_2015_us_state_500k.shp"
  county_shp = "/home/rbehnke/Desktop/cb_2015_us_county_500k/cb_2015_us_county_500k.shp"
  ;roads = "/home/rbehnke/Desktop/ne_10m_roads/ne_10m_roads.shp"
  lakes = "/home/rbehnke/Desktop/Lakes/hydrography_p_lakes_v2.shp"
  ;rivers = "/home/rbehnke/Desktop/Rivers/hydrography_l_rivers_v2.shp"

  lnres                  = True
  lnres@gsLineColor      = "black"
  lnres@gsLineThicknessF = 2.0   

  states_id = gsn_add_shapefile_polylines(wks,map,states_shp,lnres)
  county_id = gsn_add_shapefile_polylines(wks,map,county_shp,lnres)
  ;roads_id = gsn_add_shapefile_polylines(wks,map,roads,lnres)
  lakes_id = gsn_add_shapefile_polylines(wks,map,lakes,lnres)
  ;rivers_id = gsn_add_shapefile_polylines(wks,map,rivers,lnres)

;---Add point
  xyres             = True
  xyres@gsMarkerColor = "black"
  xyres@gsMarkerIndex = 1
  xyres@gsMarkerSizeF = 40.
  y = 47.725            ; Latitude of point
  x = -114.22           ; Longitude of point

  xy_id = gsn_add_polymarker(wks, map, x, y, xyres)

  draw(map)   ; This will draw all overlaid plots and the map
  frame(wks)
end