;**************************************************
; 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_0000_000.grb2  | less
;***********************************************
begin

; --- Read RUC GRIB file------------;
    dir  = "/atmomounts/home/grad/mserino/Desktop/"
    fil  = "ruc2anl_252_20090605_2000_000.grb"

    ; force a 'time' dimension
    setfileoption("grb","SingleElementDimensions","Initial_time")

    f     = addfile(dir+fil,"r")
    p     = int2flt(f->lv_ISBL2)               ;  ( lv_ISBL2 )
    time  = f->initial_time0_encoded  ; yyyymmddhh.hh_frac

    ; RUC grid point locations
    lat2d = f->gridlat_252              ; ( ygrid_0, xgrid_0 )
    lon2d = f->gridlon_252
    print("lat2d: min="+min(lat2d)+"  ;  max="+max(lat2d))
    print("lon2d: min="+min(lon2d)+"  ;  max="+max(lon2d))

    ;p     = p*0.01                    ; change units
    p@units = "hPa"                   ; skewT, mixhum_ptrh use mb (hPa)

    ; --- Specify one or more locations ***** POSSIBLE TO CREATE MULTIPLE SOUNDINGS FROM MULTIPLE LOCATIONS? *****

    lat   = (/ 37.25    , 43.25 /)
    lon   = (/-107.75 ,-101.75   /)
    npts  = dimsizes(lat)

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

    skewtOpts                 = True
    skewtOpts@DrawColAreaFill = True    ; default is False

    dataOpts   = True
    dataOpts@PrintZ = True

    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(0,:,nn,mm)
        z    = f->HGT_P0_252_ISBL(0,:,nn,mm)
        rh   = f->R_H_P0_252_ISBL(0,:,nn,mm)
        u    = f->U_GRD_252_ISBL(0,:,nn,mm)
        v    = f->V_GRD_252_ISBL(0,:,nn,mm)

        ; change units and calculate needed variables
        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"

        wspd = sqrt(u^2 + v^2)
        wdir = wind_direction(u,v,0)

        itime= toint(time)
        skewtOpts@tiMainString    = "RUC: "+itime+": ("+lat(n)+","+lon(n)+")"

        ; each location will have a different file name
        wks  = gsn_open_wks ("png", "ruc2anl_skewt_"+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