;**************************************************
; Concepts illustrated:
;   - Reading RUC (Rapid Update Cycle)  GRIB1 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*******************
    dir = "./"
    fil = "ruc2anl_252_20090605_1200_000.grb" ; ****************CHANGE FOR NEW CASE*******************
   
    ; force a 'time' dimension
    setfileoption("grb","SingleElementDimensions","Initial_time")
    f = addfile(dir+fil,"r")

    time = f->initial_time0_encoded

    ; RUC grid point locations
    lat2d = f->gridlat_252
    lon2d = f->gridlon_252
   
    ; --- 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

    nt = 0

    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_252_ISBL(nt,:,nn,mm)
        z  = f->HGT_252_ISBL(nt,:,nn,mm)
        rh = f->R_H_252_ISBL(nt,:,nn,mm)
        u  = f->U_GRD_252_ISBL(nt,:,nn,mm)
        v  = f->V_GRD_252_ISBL(nt,:,nn,mm)
	p  = f->lv_ISBL3                                 ; hPa
print(p)

	psfc = f->PRES_252_SFC(nt,nn,mm)		; "ground" pressure => pressure (Pa)

        ; change units and calculate needed variables
	psfc = psfc/100
	psfc@units = "hPa"
print("psfc="+psfc)

        tc = tk - 273.15
        tc@units = "degC"

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

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

        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))
       ;print(p)
       ;print(psfc)

	lv = ind(p.le.psfc)	; indices all standard levels < psfc (extrapolated)

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