;   Example 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_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
 simname = "NARR2d_UCM_Snudge_WangRec53wave_CU11"
 simname_short = "NARR2d_UCM_Snudge_WangRec53wave_CU11"
; Great Plains
;  ixs = 111
  ixs = 121
  ixe = 145
  jys = 61
  jye = 112


;; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
 type = "eps"
; type = "ncgm"
  folder = "wrf"+simname+".2005060100"
  idomain =1
  files = systemfunc("ls ../"+folder+"/wrfout_d0"+idomain+"*ncea.nc")
do ihours = 0, dimsizes(files)-1
  ifile_in_name = ihours
  figurename =  "wrfout_d0"+idomain+"_skewt_"+ifile_in_name; + 25 ; follow balloon
  filename = files(ihours)+".nc"
  a=addfile(filename,"r")

  wks = gsn_open_wks(type,figurename)

  gsn_define_colormap(wks,"WhViBlGrYeOrReWh")

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

    tc  = wrf_user_getvar(a,"tc",0)         ; T in C
    td  = wrf_user_getvar(a,"td",0)         ; dew point temperature
    p   = wrf_user_getvar(a, "pressure",0)  ; grid point pressure
    z   = wrf_user_getvar(a, "z",0)         ; grid point height
    uvm = wrf_user_getvar(a,"uvmet",0)      ; umet and vmet averaged to mass points
					     ; This is a 4D array where
					     ; uvm(0,:,:,:) is umet, and
					     ; uvm(1,:,:,:) is vmet, and
					     ; This function rotate winds to earth coord.

    ; 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 



; 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       = False   ; draw lines in black
  skewtOpts@DrawColAreaFill   = True    ; 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


; Get the skew-T background
  skewtOpts@tiMainString = "Sounding at GPs " +  chartostring(a->Times(0, :))  + " UTC JJA2005 "+ simname
  skewt_bkgd             = skewT_BackGround (wks, skewtOpts)
  draw (skewt_bkgd)


; Draw the skew-T plot
  dataOpts           = True
  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]
  skewT_data = skewT_PlotData(wks, skewt_bkgd, dim_avg_Wrap(dim_avg_Wrap(p(:,jys:jye,ixs:ixe))), \
                                              dim_avg_Wrap(dim_avg_Wrap(tc(:,jys:jye,ixs:ixe))), \
                                              dim_avg_Wrap(dim_avg_Wrap(td(:,jys:jye,ixs:ixe))), \
                                               dim_avg_Wrap(dim_avg_Wrap(z(:,jys:jye,ixs:ixe))), \
                                               dim_avg_Wrap(dim_avg_Wrap(u(:,jys:jye,ixs:ixe))), \
                                               dim_avg_Wrap(dim_avg_Wrap(v(:,jys:jye,ixs:ixe))), \
                                              dataOpts)

; Close the frame
  frame(wks)
  print("finish "+figurename+".eps")

  system("convert -trim -density 300 -resize 1000x1000 "+figurename+".eps /nsftor/xhu/public_html/WRF-UCM/FocusON_OKC/WRFV3.7.1/with_rr_fixed/"+folder+"/"+figurename+".png")
 system("rm "+figurename+".eps")

  end do ; 
end