;----------------------------------------------------------------------
; wrf_wind_reproj.ncl
;----------------------------------------------------------------------
; Concepts illustrated:
;   - using uvmet to rotate WRF grid-relative winds to earth-relative
;     and confirming this by reprojecting WRF output.
;   - Using gsn_csm_contour_map to plot WRF-ARW data
;   - Drawing a WRF lat/lon grid using gsn_coordinates
;----------------------------------------------------------------------
; This example is similar to wrf_gsn_3.ncl, except it zooms in further
; on the map and draws the WRF grid as points.
; Contributed March 2016 by David Ovens at
; http://www.atmos.washington.edu/~ovens/wrfwinds.html
;----------------------------------------------------------------------; 
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.
  if (.not. isvar("f")) then      ; is f on command line?
      f = "./wrfout_d1.uvrewrite.nc";
  end if
  a = addfile(f,"r")   ; Open a file

;---Read terrain height and lat/lon off file.
  it        = 0     ; first time step
  lats      = wrf_user_getvar(a,"XLAT",it)   ; latitude/longitude
  lons      = wrf_user_getvar(a,"XLONG",it)  ; latitude/longitude

  if(isatt(a,"STAND_LON")) then
     cenlon = a@STAND_LON
  else
     cenlon = a@CEN_LON
  end if
  print("cenlon = "+cenlon)

  u  = wrf_user_getvar(a,"U10",it)        ; u10 already on mass points
  v  = wrf_user_getvar(a,"V10",it)        ; v10 already on mass points
  t2m  = wrf_user_getvar(a,"T2",it)       ; temp at 2 meters in kelvin
  dims = dimsizes(t2m)                     ;
  print("dims = "+dims)
  t2m  = t2m - 273.15                     ; now in celsius
  nd = dimsizes(dims);
  dimU = dims(nd-1) ;
  dimV = dims(nd-2) ;
  print("dimU = "+dimU+" dimV = "+dimV)

  lllat = lats(0,0) ;
  lllon = lons(0,0) ;
  urlat = lats(dimV-1,dimU-1) ;
  urlon = lons(dimV-1,dimU-1) ;
  print("ll latlon ="+lllat+lllon+" ur latlon ="+urlat+urlon)

; NOTE: a lot of people assume that cen_lon is correct, but truly
; you want to use STAND_LON.  This is because you can set your
; longitude that goes exactly vertical in your WRF domain to be off
; the center.  STAND_LON is what you want.
;
; no need to do this:
;  uvm = wrf_uvmet(u,v,hgt@lat2d,hgt@lon2d,-121.,0.716)
; NCL can do the transformation internally:
  uvm = wrf_user_getvar(a,"uvmet10",it) ; get it 
; uc and vc calculated the long way:
;  cosalpha  = wrf_user_getvar(a,"COSALPHA",it)        ; 
;  sinalpha  = wrf_user_getvar(a,"SINALPHA",it)        ; 
;  uc = 1.94386 * (u*cosalpha - v*sinalpha)
;  vc = 1.94386 * (v*cosalpha + u*sinalpha)
; are identical to u,v calculated the easy way:
  u = uvm(0,:,:)*1.94386     ; kts
  v = uvm(1,:,:)*1.94386     ; kts
  u@units = "kts"
  v@units = "kts"
  t2m@lat2d = lats           ;  required for plotting
  t2m@lon2d = lons;          ;  required for plotting
  t2m@description = "2m Temperature"
  t2m@units = "C"
  u@lat2d = lats           ;  required for plotting
  u@lon2d = lons;          ;  required for plotting
  v@lat2d = lats           ;  required for plotting
  v@lon2d = lons;          ;  required for plotting
;  uc@lat2d = lats           ;  required for plotting
;  uc@lon2d = lons;          ;  required for plotting
;  vc@lat2d = lats           ;  required for plotting
;  vc@lon2d = lons;          ;  required for plotting

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

;---Set some basic plot options
  res               = True

  res@gsnMaximize   = True   ; maximize plot in frame
  res@gsnDraw       = False  ; don't draw
  res@gsnFrame      = False  ; don't advance frame
  res@gsnMaximize   = True   ; use full page
  res@cnInfoLabelOn = False  ; turn off cn info label

  res@mpProjection  = "CylindricalEquidistant"    ; The default

; Usually, when data is placed onto a map, it is TRANSFORMED to the
; specified projection. Since this model is already on a native lambert
; conformal grid, we want to turn OFF the tranformation.
;
  res@tfDoNDCOverlay           = False

;---Zoom in on plot
  res@mpMinLatF     = urlat-5
  res@mpMaxLatF     = urlat+3.
  res@mpMinLonF     = urlon-10
  res@mpMaxLonF     = urlon

  res@gsnAddCyclic             = False            ; regional data 
  res@vcRefMagnitudeF          = 10.0             ; define vector ref mag
  res@vcRefLengthF             = 0.04            ; define length of vec ref
  res@vcGlyphStyle             = "WindBarb"       ; WindBarb or CurlyVector
  res@vcPositionMode           = "ArrowHead"      ; stick tail at value
  res@vcMinDistanceF           = 0.001
  res@gsnScalarContour        = True
  res@cnFillOn                = True
  res@cnFillDrawOrder         = "Draw"
  res@cnFillPalette           = "gui_default"     ; set color map
  res@cnLinesOn               = False             ; turn off contour lines
  res@cnLineLabelsOn          = False             ; turn off contour lines

;---Additional resources desired
  res@pmTickMarkDisplayMode = "Always"   ; nicer tickmarks

  res@mpDataBaseVersion     = "MediumRes"       ; better and more map outlines
  res@mpDataSetName         = "Earth..4"
  res@mpOutlineBoundarySets = "AllBoundaries"
  res@mpOutlineOn           = True
  res@lbLabelBarOn          = False
  res@lbTitleOn             = False

  res@gsnCenterString          = "NCL:  Rotation handled with: ~C~ ~Z75~ uvm = wrf_user_getvar(a,uvmet10,it)"
  res@gsnLeftString          = " ";
  res@gsnRightString          = " ";
  plot = gsn_csm_vector_scalar_map(wks,u,v,t2m,res)
;  plot = gsn_csm_vector_scalar_map(wks,uc,vc,t2m,res)

  mkres                 = True
  mkres@gsMarkerSizeF  = 0.008
  mkres@gsMarkerColor  = "blue"
  gsn_coordinates(wks,plot,t2m,mkres)

end