;   Functioning script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;   Create a single skewT plot


;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "/grid-scratch/kshourd/skewtBackground.ncl"      ; Updated bckgrnd for better visualization and readability
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl" <-- original skewT bckgrnd

begin

;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.

; YOUR FILE HERE
  a = addfile("/ScenarioB/Part1/WyomingWeaMod/WRF361/WRFV3/test/em_real/2005/wrfout_d03_2004-12-15_00:00:00.nc","r")
  
  ;FirstTime = True
  ;times  = wrf_user_list_times(a)  ; get times in the file
  ;ntimes = dimsizes(times)         ; number of times in the file
  ;print("There are " +ntimes+ " times in this file")


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

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

 ; print("Working on time: " + times(it) )

 ; res           = True
 ; res@TimeLabel = times(it)   ; Set Valid time to use on plots

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

; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  type = "png"
  wks_type = type
  ;wks_type@wkWidth = 2500
  ;wks_type@wkHeight = 2500




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need at time it
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    it   = 0                                    
    tk   = wrf_user_getvar(a, "tk", it)         ; T in K
    tc   = tk - 273.15
    tc@units = "degC"
    td   = wrf_user_getvar(a, "td", it)                       ; dew point temperature in C
    p    = wrf_user_getvar(a, "pressure", it)        ; grid point pressure
    z    = wrf_user_getvar(a, "z", it)          ; grid point height
    u    = wrf_user_getvar(a, "U", it)
    v    = wrf_user_getvar(a, "V", it)
    ;uvm  = wrf_user_getvar(a, "uvm", it)         ; umet and vmet averaged to mass points, lambert conformal grid


  wks = gsn_open_wks(wks_type,"laramieSkew")
  gsn_define_colormap(wks,"WhViBlGrYeOrReWh")

    ; extract u and v from uvm array, and turn wind into kts
    ; u   = uvm(0,:,:,:)*1.94386 
    ; v   = uvm(1,:,:,:)*1.94386

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Get ij point in model domain for our location of interest
; loc(1) is south-north (y) and loc(0) is west-east (x)
; Ensure this point is in the domain, as we are not checking
; We are interested in NCL array pointers, so subtract 1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 

  loc = wrf_user_ll_to_ij(a, -105.87, 42.87, True)
  locX = loc(0) - 1
  locY = loc(0) - 1


; Define a few skew-T plotting options
  skewtOpts                   = True 
  skewtOpts@DrawHeightScale   = True    ; plot height scale on side
  skewtOpts@DrawHeightScaleFt = False   ; plot height scale in km
  skewtOpts@DrawStandardAtm   = True    ; draw standard atm on plot
  skewtOpts@vpXF              = 0.12    ; controls off-set from left
  skewtOpts@vpYF              = 0.87    ; controls off-set from top
  skewtOpts@vpWidthF          = 0.75    ; controls size of plot
  skewtOpts@vpHeightF         = 0.75    ; controls size of plot
  skewtOpts@DrawFahrenheit    = False   ; use deg C scale
  skewtOpts@tiMainFontHeightF = 0.015   ; change height of main title
  skewtOpts@DrawColLine       = True    ; draw lines in black
  skewtOpts@DrawColAreaFill   = False    ; color on background plot
 ;skewtOpts@DrawColAreaColor  = "Green" ; final color may depend on the color table used
 ;skewtOpts@DrawColAreaColor  = 53      ; Light Green for WhViBlGrYeOrReWh color table
  skewtOpts@PrintOpts         = False   ; do not print options out
  skewtOpts@WindBarbLineThicknessF = 6.0
  skewtOpts@WindBarbMarkerSize = 5.0

; Get the skew-T background
  skewtOpts@tiMainString = "Wyoming Laramie Mountains (42.87, -105.87), 12/15/2004 0:00 Z" 
  skewt_bkgd             = skewT_BackGround (wks, skewtOpts)
  draw (skewt_bkgd)


; Draw the skew-T plot
  dataOpts           = True
  dataOpts@colTemperature = "red"
  dataOpts@lineThicknessTemperature = 6
  dataOpts@colDewPt  = 42 
  dataOpts@lineThicknessDewPt = 6
  dataOpts@Parcel    = 1
  dataOpts@WspdWdir  = False  ; wind speed and dir [else: u,v]
  dataOpts@HspdHdir  = True   ; wind speed and dir [else: u,v]
  dataOpts@PlotWindH = False  ; plot wind barbs at h lvls [pibal; special]

  dataOpts@vcWindBarbLineThicknessF = 4.0               ; set the wind barb thickness
  dataOpts@WindBarbLineThicknessF = 4.0           ; set the wind barb thickness



  skewT_data = skewT_PlotData(wks, skewt_bkgd, p(:,locY,locX), \
                                              tc(:,locY,locX), \
                                              td(:,locY,locX), \
                                               z(:,locY,locX), \
                                               u(:,locY,locX), \
                                               v(:,locY,locX), \
                                              dataOpts)


; Close the frame
  frame(wks)

  ;end do 

end