;**************************************************
; Concepts illustrated:
;   - Reading RUC (Rapid Update Cycle) GRIB2 data
;   - Using getind_latlon2d to determine grid locations
;   - Drawing Skew-T plots at nearest grid locations to user specified locations
;**************************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/skewt_func.ncl"
;***********************************************
; RUC Data downloaded from:
;     http://nomads.ncdc.noaa.gov/data.php?name=access#hires_weather_datasets
; RUC-ANL: http
;     Specifically: http://nomads.ncdc.noaa.gov/data/rucanl/201205/20120501/
;***********************************************
; The GRIB file's contents can be examined via:
;     ncl_filedump -itime ruc2anl_130_20120501_2200_000.grb2  | less
;***********************************************
begin

; --- Read RUC GRIB file------------;
    dir = "/strm1/serino/DATA/06052009/RUC/" ; ****************CHANGE FOR NEW CASE*******************
    fil = "ruc2anl_130_20090605_1200_000.grb2" ; ****************CHANGE FOR NEW CASE*******************
    
    ; force a 'time' dimension
    setfileoption("grb2","SingleElementDimensions","Initial_time")
    f = addfile(dir+fil,"r")

    p = (f->lv_ISBL0)
    time = f->initial_time0_encoded

    ; RUC grid point locations
    lat2d = f->gridlat_0
    lon2d = f->gridlon_0
   
    ; --- Specify one or more locations
    lat = (/39.781/) ; lat = (/41.147, 43.25/) <-- FOR MULTIPLE LOCATIONS ****************CHANGE FOR NEW CASE*******************
    lon = (/-104.54/) ; lon = (/-104.801, -101.75/) <-- FOR MULTIPLE LOCATIONS ****************CHANGE FOR NEW CASE*******************
    npts = dimsizes(lat)

;*************************
;   create plot(s)
;*************************

    skewtOpts = True
    skewtOpts@DrawColAreaFill = False
    skewtOpts@DrawFahrenheit = False
    skewtOpts@DrawHeightScale = True

    dataOpts = True
    dataOpts@PrintZ = True
    dataOpts@colTemperature = "red"
    dataOpts@colDewPt = "green"
    dataOpts@colCape = "black"
    dataOpts@linePatternCape = 8

    do n=0,npts-1           ; loop over each grid pt

        ; find grid point nearest the user specified location
        nm  = getind_latlon2d (lat2d,lon2d, lat(n), lon(n))
        nn  = nm(0,0)
        mm  = nm(0,1)
        print("location=("+lat(n)+","+lon(n)+")  grid=("+lat2d(nn,mm)+","+lon2d(nn,mm)+")")

        tk = f->TMP_P0_L100_GLC0(0,:,nn,mm)
        z = f->HGT_P0_L100_GLC0(0,:,nn,mm)
        rh = f->RH_P0_L100_GLC0(0,:,nn,mm)
        u = f->UGRD_P0_L100_GLC0(0,:,nn,mm)
        v = f->VGRD_P0_L100_GLC0(0,:,nn,mm)

        ; change units and calculate needed variables
        tc = tk - 273.15
        tc@units = "degC"

	p = p / 100.
	p = floattointeger(p)
	p@units = "hPa"

        q = mixhum_ptrh (p, tk, rh, 2)
        q@units = "kg/kg"

        tdc = dewtemp_trh(tk,rh) - 273.15
;	print(tc+","+tdc)
        tdc@units = "degC"

        wspd = sqrt(u^2 + v^2)
        wdir = wind_direction(u,v,0)
;       print(wspd+" "+u+" "+v+" "+wdir)
	 
        print(sprintf("%6.1f", p)+sprintf("%6.1f",  tc))

        itime= toint(time)
        skewtOpts@tiMainString = "R2A13: "+itime+" - "+lat+", "+lon+""
       
        ; each location will have a different file name
        wks = gsn_open_wks ("png", "ruc2anl13_"+lat+"_"+lon+"_"+itime+"") ;"_"+sprinti("%0.3i",n))
        skewt_bkgd = skewT_BackGround (wks, skewtOpts)
        skewt_data = skewT_PlotData (wks, skewt_bkgd, p,tc,tdc,z, wspd,wdir, dataOpts)
        draw (skewt_bkgd)
        draw (skewt_data)
        frame(wks)
    end do
end