;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
;load "./WRFUserARW.ncl"

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("./wrfout_d01_2019-11-07_00:00:00.nc","r")

; We generate plots, but what kind do we prefer?
; type = "x11"
  type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_Surface3")
  gsn_define_colormap(wks,"MPL_Pastel1")

; Set some basic resources
  res = True
  res@MainTitle                   = "WRF MODEL OUTPUT"

  pltres = True
  mpres = True


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; What times and how many time steps are in the data set?
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  do it = 0,ntimes-1,1             ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it)   ; Set Valid time to use on plots


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    slp = wrf_user_getvar(a,"slp",it)    ; slp
      wrf_smooth_2d( slp, 3 )            ; smooth slp
    tc = wrf_user_getvar(a,"tc",it)      ; 3D tc 
    td = wrf_user_getvar(a,"td",it)      ; 3D td 
    u  = wrf_user_getvar(a,"ua",it)      ; 3D U at mass points
    v  = wrf_user_getvar(a,"va",it)      ; 3D V at mass points
    td2 =  wrf_user_getvar(a,"td2",it)   ; Td2 in C
    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
       tc2 = tc2-273.16                  ; T2 in C
    u10 = wrf_user_getvar(a,"U10",it)    ; u at 10 m, mass point
    v10 = wrf_user_getvar(a,"V10",it)    ; v at 10 m, mass point
    hgt = wrf_user_getvar(a,"HGT",it)    ; topography height
    wsp = sqrt(u10*u10 + v10*v10)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ; Plotting options for HGT                
      opts = res                         
      opts@cnFillOn = True  
      opts@cnLevelSelectionMode = "ExplicitLevels"
;     opts@ContourParameters = (/ 0., 1300., 100./)
      opts@cnLevels    = (/ 0.,1.,2.,5.,10.,15.,20.,25.,30.,35.,40.,60./)
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_ht = wrf_contour(a,wks,wsp,opts)
      delete(opts)
    
; Plotting options for Wind Vectors                 
      opts = res         
      opts@FieldTitle = "Wind"       ; overwrite Field Title
      opts@NumVectors = 20           ; density of wind barbs
      opts@vcGlyphStyle = "CurlyVector"     ; turn off curly vectors
;     opts@vcLineArrowWidthF = 0.1    ; thickness of arrow 
      opts@vcMinFracLengthF  = 0.0 
      opts@vcRefMagnitudeF           = 0.0             ; define vector ref mag
      opts@vcRefLengthF              = 0.045           ; define length of vec ref
      opts@vcRefAnnoOn               = True
      vector = wrf_vector(a,wks,u10,v10,opts)
      delete(opts)
 

    ; MAKE PLOTS                                       
      plot = wrf_map_overlays(a,wks,(/contour_ht,vector/),pltres,mpres)
;     mpNationalLineThicknessF = 2.0
;     mpres@mpNationalLineThicknessF = 2.0 ; interior boundaries 
;     mpres@mpDataBaseVersion = "MediumRes"
;     mpres@mpGeophysicalLineThicknessF = 1.0 ; lines separating land/ocean 
;      plot = wrf_map_overlays(a,wks,(/contour_td,vector/),pltres,mpres)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  end do        ; END OF TIME LOOP

end