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"

begin

; --- Read RUC GRIB file------------;
	dir = “./“
	;dir = "/strm1/serino/DATA/12312010/RUC/"
    fil = "ruc2anl_130_20101231_1600_000.grb2"
    
    ; force a 'time' dimension
    setfileoption("grb2","SingleElementDimensions","Initial_time")
    f = addfile(dir+fil,"r")

    time = f->initial_time0_encoded

    ; RUC grid point locations
    lat2d = f->gridlat_0
    lon2d = f->gridlon_0
   
    ; --- Specify one or more locations
    lat = (/37.380176905368884, 37.47027868555142, 37.38012300713148, 37.29007374881771, 37.38012300713148/)
    lon = (/-90.6104383196949, -90.6104383196949, -90.49752896972296, -90.6104383196949, -90.72334766966681/)
    
	npts= dimsizes(lat)

;-------Create IndividualplotsPanels----------

    wks = gsn_open_wks ("png", "skewt") ;save as png

    skewtOpts = True
    skewtOpts@DrawColAreaFill = False
    skewtOpts@DrawFahrenheit  = False
    skewtOpts@DrawHeightScale = False
    skewtOpts@vpWidthF        = 0.27 ;0.38   ; skewT default is 0.85
    skewtOpts@vpHeightF       = 0.27 ;0.38   ;                  0.85
    skewtOpts@DrawWind        = True    ; not possible with panel
    skewtOpts@DrawFahrenheit = False  ; default is True 
	skewtOpts@DrawColAreaFill = False   ; default is False
    skewtOpts@tiMainOffsetXF  = 0.0     ; skewT default is -0.1
    skewtOpts@tmXBLabelFontHeightF = 0.0125   ; 0.14 is default
    skewtOpts@tmYLLabelFontHeightF = 0.0125   ; 0.14 is default
    skewtOpts@tiMainFontHeightF    = 0.015 ;0.0150   ; 0.025 is default
	;skewtOpts@tiXAxisFontHeightF = 0.0001
	;skewtOpts@tiYAxisFontHeightF = 0.0001
    skewtOpts@DrawFahrenheit = False  ; default is True 
    skewtOpts@tiYAxisString  = ""
    skewtOpts@tiXAxisString  = ""

    dataOpts = True
    dataOpts@PrintZ = True             ; place geopotential hgt on figure
    dataOpts@colTemperature = "red"
    dataOpts@colDewPt = "green"
    dataOpts@colCape = "black"
    dataOpts@linePatternCape = 8

    nplt = npts
    plot = new (nplt, "graphic")

    nt = 0
    do n=0,4   ; 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"+n+"=("+lat(n)+","+lon(n)+")  grid=("+lat2d(nn,mm)+","+lon2d(nn,mm)+")")

        tk = f->TMP_P0_L100_GLC0(nt,:,nn,mm)
        z  = f->HGT_P0_L100_GLC0(nt,:,nn,mm)
        rh = f->RH_P0_L100_GLC0(nt,:,nn,mm)
        u  = f->UGRD_P0_L100_GLC0(nt,:,nn,mm)
        v  = f->VGRD_P0_L100_GLC0(nt,:,nn,mm)
		p  = f->lv_ISBL0
                             
        psfc = f->PRES_P0_L1_GLC0(nt,nn,mm)      ; scalar; "ground" pressure  => pressure
		cape = f->CAPE_P0_L1_GLC0(nt,nn,mm)      ; scalar

        ; change units and calculate needed variables
		psfc = psfc/100
		psfc@units = "hPa"

		p  = p/100
		p@units = "hPa"
       
		tc = tk - 273.15
       ;copy_VarMeta(tk,tc)
        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)
		
		;print(ind(p.le.psfc))
		lv := ind(p.le.psfc)     ; indices all standard levels < psfc (extrapolated); ':=' allows reassignment of variables without having to delete them first
		;print(lv)
        klv = dimsizes(lv)
	 
        wspd = sqrt(u^2 + v^2)
        wdir = wind_direction(u,v,0)
		;print(wspd+" "+u+" "+v+" "+wdir)

        itime = toint(time)

		; write data to file
		fields = (/"Pressure(hPa)", "Height(gpm)", "Temperature(C)", "Dew Point(C)", "Wind Speed(m/s)", "Wind Direction(deg)", "CAPE(J/kg)"/)
		fields = ""+fields+"" 
		header = [/str_join(fields," ")/]
		CAPE = conform(p,cape,-1)  ; convenience for printing
		data = [/p(lv),z(lv),tc(lv),tdc(lv),wspd(lv),wdir(lv),CAPE/]

		write_table(n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv", "w", header, "%s")		;write header
		write_table(n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv", "a", data, "%g%g%g%g%g%g%g")	;fill in data
		;print("====== Contents of '" +"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv" + "' ======")
	  	system("cat " + n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv")

		if (n.eq.0) then
			skewtOpts@vpXF      = 0.40 ;0.08        ; upper left
			skewtOpts@vpYF      = 0.95
		end if
		
		if (n.eq.1) then
			 skewtOpts@vpXF      = 0.07 ;0.56        ; upper right
			 skewtOpts@vpYF      = 0.75 ;0.95
		end if
		
		if (n.eq.2) then
			skewtOpts@vpXF      = 0.72 ;0.08        ; lower left
			skewtOpts@vpYF      = 0.75 ;0.46
		end if
		
		if (n.eq.3) then
			skewtOpts@vpXF      = 0.07 ;0.57        ; lower Right
			skewtOpts@vpYF      = 0.36 ;0.46
		end if

		if (n.eq.4) then
			skewtOpts@vpXF      = 0.72 ;0.57        ; lower Right
		    skewtOpts@vpYF      = 0.36 ;0.46
		end if

		skewtOpts@tiMainString = itime+": "+sprintf("%5.2f",lat(n))+", "+sprintf("%5.2f",lon(n));+" CAPE="+toint(cape)
		skewt_bkgd = skewT_BackGround (wks, skewtOpts)
		plot(n)   = skewT_PlotData (wks,skewt_bkgd,p(lv),tc(lv),tdc(lv),z(lv),wspd(lv),wdir(lv),dataOpts)

		draw (skewt_bkgd)
		draw (plot(n))
                
	end do
	frame(wks)

end