;************************************************** ; 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" ;*********************************************** ; 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/05082009/RUC/" ;****************CHANGE FOR NEW CASE******************* dir = "/atmomounts/home/grad/mserino/Downloads/" fil = "rap_130_20140407_1900_000.grb2" ;****************CHANGE FOR NEW CASE******************* ; --- Specify one or more locations lat = (/35.24750752118764, 35.337641535679936, 35.24745764464769, 35.15737216143987, 35.24745764464769/) ;****************CHANGE FOR NEW CASE******************* lon = (/-76.53402951501053, -76.53402951501053, -76.42415449744674, -76.53402951501053, -76.64390453257433/) ;****************CHANGE FOR NEW CASE******************* inflow_elevation = 0.3048 ;meters ; 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 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@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) ; "ground" pressure => pressure cape = f->CAPE_P0_L1_GLC0(nt,nn,mm) ; change units and calculate needed variables psfc = psfc/100 psfc@units = "hPa" p = p/100 p@units = "hPa" 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" lv := ind(p.le.psfc) ; indices all standard levels < psfc (extrapolated); ':=' allows reassignment of variables without having to delete them first klv = dimsizes(lv) wspd = sqrt(u^2 + v^2) wdir = wind_direction(u,v,0) ;print(wspd+" "+u+" "+v+" "+wdir) if (n.eq.0) then metersAGL := z(lv) - inflow_elevation end if itime = toint(time) CAPE = conform(p,cape,-1) ; convenience for printing ; write data to file if (n.eq.0) then fields0 = (/"Pressure(hPa)", "Height(gpm)", "Temperature(C)", "Dew Point(C)", "Relative Humidity(%)", "Wind Speed(m/s)", "Wind Direction(deg)", "CAPE(J/kg)", "Height(mAGL)"/) fields0 = ""+fields0+"" header0 = [/str_join(fields0," ")/] data0 = [/p(lv),z(lv),tc(lv),tdc(lv),rh(lv),wspd(lv),wdir(lv),CAPE,metersAGL/] write_table(n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv", "w", header0, "%s") ;write header write_table(n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv", "a", data0, "%f %f %f %f %f %f %f %f %f") ;fill in data, tab delimited ;print("====== Contents of '" +"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv" + "' ======") system("cat " + n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv") else fields = (/"Pressure(hPa)", "Height(gpm)", "Temperature(C)", "Dew Point(C)", "Relative Humidity(%)", "Wind Speed(m/s)", "Wind Direction(deg)", "CAPE(J/kg)"/) fields = ""+fields+"" header = [/str_join(fields," ")/] data = [/p(lv),z(lv),tc(lv),tdc(lv),rh(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, "%f %f %f %f %f %f %f %f") ;fill in data, tab delimited ;print("====== Contents of '" +"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv" + "' ======") system("cat " + n+"Analysis13_"+lat(n)+"_"+lon(n)+"_"+itime+".csv") end if 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