;----------------------------------------------------------------------
; This script shows how to create a basic plot of temperature (degC)
; calculated from a WRF file.
;
; This script is meant to be compared to one that uses wrf_xxxx
; functions to do the plotting. 
;----------------------------------------------------------------------

 filename = "wrfout_d02_2013-02-08_15:30:00"
 f = addfile(filename,"r")

 tc = wrf_user_getvar(f,"tc",0)    ; Temperature (C)
 xlat2d = wrf_user_getvar(f,"XLAT",0)    ; Temperature (C)
 xlon2d = wrf_user_getvar(f,"XLONG",0)    ; Temperature (C)

 printVarSummary(tc)

 wks = gsn_open_wks("pdf","wrf_demo_plot_tc_gsn")

 res                = True
  res@pmTickMarkDisplayMode = "Always"
  res@mpFillOn              = False           ; turn off map fill
  res@mpOutlineDrawOrder    = "PostDraw"      ; draw continental outline last
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates" ; state boundaries
  res@mpLimbLineColor        = "Black"
  res@mpPerimLineColor       = "Black"
  res@mpPerimLineThicknessF  = 10
  res@mpUSStateLineThicknessF = 10.0
  res@mpGeophysicalLineThicknessF = 10.0
  res@mpGridAndLimbOn = False

 res@cnFillOn       = True  
;res@cnFillPalette  = "WhiteBlueGreenYellowRed"
 res@cnFillPalette  = "NCV_jet"
 res@cnLinesOn      = False

;--- These are required for gsn_csm_contour_map plotting
;res = wrf_map_resources(f,res)
 res@tfDoNDCOverlay = True
 res@gsnAddCyclic   = False

 nl = 0     ; level index
;plot = gsn_csm_contour_map(wks,tc(nl,:,:),res)
 dims3d = dimsizes(tc) 
 print(dims3d)
 y_start = 2*dims3d(1)/5
 y_end = 3*dims3d(1)/5
 x_start = 2*dims3d(2)/5
 x_end = 3*dims3d(2)/5
 res@mpMinLatF       = xlat2d(y_start,x_start)
 res@mpMaxLatF       = xlat2d(y_end,x_end)
 res@mpMinLonF       = xlon2d(y_start,x_start)
 res@mpMaxLonF       = xlon2d(y_end,x_end)
 tc_plane = tc(nl,y_start:y_end,x_start:x_end)
 xlat2d_sub = xlat2d(y_start:y_end,x_start:x_end)
 xlon2d_sub = xlon2d(y_start:y_end,x_start:x_end)
 tc_plane@lat2d = xlat2d_sub    ; Needed for plotting
 tc_plane@lon2d = xlon2d_sub
 print("y_start = " + y_start)
 print("x_start = " + x_start)
 print("y_end = " + y_end)
 print("x_end = " + x_end)


 plot = gsn_csm_contour_map(wks,tc_plane,res)