; load libraries, twice for good measure load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl" ;load "/home/torn/ncl/wrf_map_info.ncl" ;load "/home/torn/ncl/array_stat.ncl" ;load "/home/torn/ncl/image_conv.ncl" load "heat_stress.ncl" ;load "ETo_FAO56.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/wind_rose.ncl" ;load "/home/torn/ncl/wrf_map_info.ncl" ;load "/home/torn/ncl/array_stat.ncl" ;load "/home/torn/ncl/image_conv.ncl" ;load "wrf_dart_util.ncl" load "heat_stress.ncl" ;load "ETo_FAO56.ncl" begin ;;;;; set all the stuff that doesn't change from station to station ; time stuff times = fspan(int2flt(289)*-5/60, 0, 289) currentday = systemfunc("date -u +%d") currenthour = systemfunc("date -u +%H") past24hour = systemfunc("date -u -d -24hours +%H") currentminute = systemfunc("date -u +%M") currenttime = currenthour+currentminute value = -999.0 ; fill value delim = "," quote = inttochar(34) ; now loop over all stations names=(/"SCHU","WHIT","WATE","CSQR","SPRA","COLD","JORD","WEST","BERK","BELD","ESSX","CHES","STEP","OTIS","LAUR","GABR","COBL", \ "TICO","WFMB","COPE","BATA","SBRI","SHER","DELE","WOLC","JOHN","LOUI","TULL","WARW","WALT","CINC","FRED","BRAN","WELL", \ "CAPE","BELL","WARS","BELM","COHO","BURD","TYRO","PENN","CLYM","RAND","OLEA","BURT","BROC","SCHA","REDF","HARR","OSCE", \ "WALL","MORR","CAMD","ELMI","PHIL","POTS","OSWE","CHAZ","SARA","RUSH","TUPP","CLIF","FAYE","EDIN","MEDI","GROV","OWEG", \ "NBRA","DUAN","ILAK","OPPE","REDH","COPA","ANDE","HARP","BROO","EAUR","YORK","HART","HERK","ADDI","SPRI","MALO","VOOR", \ "HAMM","EDWA","SUFF","BSPA","BEAC","DEPO","MEDU","ROXB","ELDR","GFLD","WGAT","ELLE","KIND","WBOU","OLDF","NHUD","SCIP", \ "BING","NEWC","CLAR","ONTA","CROG","RAQU","WANT","SOUT","HFAL","DOVE","GROT","SCHO","TANN","SOME","BUFF","GFAL","PISE",\ "BREW","QUEE","STAT","BKLN","MANH","BRON","STON"/) numnames = dimsizes(names)-1 day1string = new((/dimsizes(names),25/), string) day1string = "-999" ; this loops over the entire remainder of the script and makes a meteogram for each of the above stations do q=0,numnames,1 ; read in station metadata meta = asciiread(names(q)+".ascii", -1, "string") print(" ") print("Attempting "+meta(1)) print(" ") ; relative = ""+meta(1)+","+meta(0)+","+tofloat(meta(2))+","+tofloat(meta(3))+","+tofloat(meta(4)) ; system("echo "+quote+relative(0)+quote+" >> metadata.csv") ; define files, previously scp'ed and cat'ed datafile = names(q)+"data" snowdatafile = names(q)+"snowdata" ;---Read in files as array of strings so we can parse each line tnew = asciiread(datafile,-1,"string") snownew = asciiread(snowdatafile,-1,"string") ;size = 288 ; this is the size for one full day ;sitename = str_get_field(tnew,1,delim) datetime = str_get_field(tnew,2,delim) temp = tofloat(str_get_field(tnew,9,delim)) relh = tofloat(str_get_field(tnew,11,delim)) precip = tofloat(str_get_field(tnew,12,delim)) snow := tofloat(str_get_field(snownew,3,delim)) windspeed = tofloat(str_get_field(tnew,7,delim)) windmax = tofloat(str_get_field(tnew,8,delim)) windspeedp = tofloat(str_get_field(tnew,4,delim)) windmaxp = tofloat(str_get_field(tnew,5,delim)) pres = tofloat(str_get_field(tnew,14,delim)) windspeed = where(windspeed.eq.-999.0,windspeed,windspeed*2.23694) windmax = where(windmax.eq.-999.0,windmax,windmax*2.23694) windspeedp = where(windspeedp.eq.-999.0,windspeedp,windspeedp*2.23694) windmaxp = where(windmaxp.eq.-999.0,windmaxp,windmaxp*2.23694) temp = where(temp.eq.-999.0,temp,temp*9/5 + 32) ;solar = where(solar.eq.-999.0,solar,solar) precip = where(precip.eq.-999.0,precip,precip*.0393701) snow = where(snow.eq.-999.0,snow,snow*39.3701) pres = where(pres.eq.-999.0,pres,pres/100) relh@_FillValue = value precip@_FillValue = value snow@_FillValue = value pres@_FillValue = value ; this part makes sure there isn't fake snow assumulating with temps well above freezing do i = 0,dimsizes(temp)-1 snow(i+12) = where(temp(i).lt.36,snow(i+12),snow(i+12-1)) end do temp@_FillValue = value dewp = (dewtemp_trh((((temp-32)*5/9)+273.15),relh)-273.15)*9/5 + 32 ;temps = new((/3,dimsizes(temp)/), float) ;temps(0,:) = value temps = temp ;temps(2,:) = dewp dewp@_FillValue = value temps@_FillValue = value winds = new((/2,dimsizes(temp)/), float) winds(0,:) = windmax winds(1,:) = windspeed winds@_FillValue = value windsprop = new((/2,dimsizes(temp)/), float) windsprop(0,:) = windmaxp windsprop(1,:) = windspeedp windsprop@_FillValue = value if (any(ismissing(winds(0,:)))) winds(0,:)=where(ismissing(winds(0,:)),windsprop(0,:),winds(0,:)) end if if (any(ismissing(winds(1,:)))) winds(1,:)=where(ismissing(winds(1,:)),windsprop(1,:),winds(1,:)) end if winds(0,:)=where(temps.lt.(-10.0),windsprop(0,:),winds(0,:)) ;;; tries to get rid of bad sonic data in super cold temps winds(1,:)=where(temps.lt.(-10.0),windsprop(1,:),winds(1,:)) ;;; tries to get rid of bad sonic data in super cold temps gusts = winds(0,:) gusts = where(windmaxp.gt.winds(0,:),windmaxp,winds(0,:)) gusts@_FillValue = value if (any(ismissing(temps))) replace_ieeenan (temps, value, 0) temps@_FillValue = value end if temps=where(temps.gt.(-80),temps,temps@_FillValue) temps=where(temps.lt.(115),temps,temps@_FillValue) pressure = pres pressure = value if .not.all(ismissing(pres)) if .not.all(ismissing(temps)) temp = (temps-32)*5/9 elev = 0.0065*(1.5+stringtofloat(meta(4))) pressure=pres*((1-(elev/(temp+elev+273.15)))^-5.257) end if end if pressure@_FillValue = value if (any(ismissing(precip))) replace_ieeenan(precip, 0.0, 0) precip@_FillValue = -999.0 end if t0 = (/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0/) daycount = 1 do n=1,dimsizes(precip)-1,1 time = stringtochar(datetime(n)) if time(11:15).eq."00:05" t0(daycount(0))=precip(n-1) if ismissing(precip(n-1)) t0(daycount(0))=precip(n-2) end if daycount = daycount + 1 end if precip(n)=precip(n)+t0(daycount(0)-1); sum(t0) end do precip=precip-precip(0) ; This part figures out what the "feels like" temperature(s) are hi = temps if .not.all(ismissing(relh)) if any(temps.gt.70) hi = heat_index_nws(temps,relh,(/2,2/),False) end if end if wc = temps if any(temps.lt.60) if .not.(all(ismissing(winds(1,:)))) newwindday = winds(1,:)^.16 wc = 35.74+(.6215*temps)-(35.75*newwindday)+(.4275*temps*newwindday) end if end if ; make snowfall "totals" snowfall = new((/24*7/), float) snowfall = -999.0 last = dimsizes(snow)-1 increment = 12 ; this is one hour do i=0,((24*7)-1) ; hours in week if .not.all(ismissing(snow(last-(increment*(i+2)):last-(increment*(i+0))))) currentsnow = avg(snow(last-(increment*(i+1)):last-(increment*(i+0)))) pastsnow = avg(snow(last-(increment*(i+2)):last-(increment*(i+1)))) snowfall(i) = currentsnow-pastsnow end if end do snowfall@_FillValue = value accumulation=snowfall if .not.all(ismissing(snowfall)) accumulation = where(snowfall.lt.0,0,snowfall) accumulation = cumsum(accumulation,2) ; note that this counts backward in time, so that it's easier to output past 1 day end if accumulation@_FillValue = value ;realday = tempday(ind(.not.ismissing(tempday))) last = dimsizes(temp)-1 ; should be 2016 for 1 week ;print(datetime(last-288)) ;print(datetime(last-(2*288))) ; temp section ;system("echo "+quote+meta(1)+quote+" >> day1.txt") ;system("echo "+quote+meta(1)+quote+" >> day2.txt") ;system("echo "+quote+meta(1)+quote+" >> day3.txt") ;system("echo "+quote+meta(1)+quote+" >> day5.txt") ;system("echo "+quote+meta(1)+quote+" >> day7.txt") day1string(q,0) = meta(1) if num(ismissing(temps((last-288):last))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(temps((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxtemp1day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(temps((last-288):last)) system("echo "+quote+relative(0)+quote+" >> mintemp1day.dat") day1string(q,1) = tostring(max(temps((last-288):last))) day1string(q,2) = datetime(maxind(temps((last-288):last))) day1string(q,3) = tostring(min(temps((last-288):last))) day1string(q,4) = datetime(minind(temps((last-288):last))) end if if num(ismissing(temps((last-(2*288)):last))).lt.12 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(temps((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxtemp2day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(temps((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> mintemp2day.dat") end if if num(ismissing(temps((last-(3*288)):last))).lt.18 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(temps((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxtemp3day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(temps((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> mintemp3day.dat") end if if num(ismissing(temps((last-(5*288)):last))).lt.30 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(temps((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxtemp5day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(temps((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> mintemp5day.dat") end if if num(ismissing(temps)).lt.42 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(temps) system("echo "+quote+relative(0)+quote+" >> maxtemp7day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(temps) system("echo "+quote+relative(0)+quote+" >> mintemp7day.dat") end if ; pressure section if num(ismissing(pressure((last-288):last))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(pressure((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxpres1day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(pressure((last-288):last)) system("echo "+quote+relative(0)+quote+" >> minpres1day.dat") day1string(q,5) = tostring(max(pressure((last-288):last))) day1string(q,6) = datetime(maxind(pressure((last-288):last))) day1string(q,7) = tostring(min(pressure((last-288):last))) day1string(q,8) = datetime(minind(pressure((last-288):last))) end if if num(ismissing(pressure((last-(2*288)):last))).lt.12 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(pressure((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxpres2day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(pressure((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> minpres2day.dat") end if if num(ismissing(pressure((last-(3*288)):last))).lt.18 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(pressure((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxpres3day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(pressure((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> minpres3day.dat") end if if num(ismissing(pressure((last-(5*288)):last))).lt.30 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(pressure((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxpres5day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(pressure((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> minpres5day.dat") end if if num(ismissing(pressure)).lt.42 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(pressure) system("echo "+quote+relative(0)+quote+" >> maxpres7day.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(pressure) system("echo "+quote+relative(0)+quote+" >> minpres7day.dat") end if ; gust section if .not.all(ismissing(gusts((last-288):last))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(gusts((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxwind1day.dat") day1string(q,9) = tostring(max(gusts((last-288):last))) day1string(q,10) = datetime(maxind(gusts((last-288):last))) end if if .not.all(ismissing(gusts((last-(2*288)):last))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(gusts((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwind2day.dat") end if if .not.all(ismissing(gusts((last-(3*288)):last))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(gusts((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwind3day.dat") end if if .not.all(ismissing(gusts((last-(5*288)):last))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(gusts((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwind5day.dat") end if if .not.all(ismissing(gusts)) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(gusts) system("echo "+quote+relative(0)+quote+" >> maxwind7day.dat") end if ; heat index section if num(ismissing(hi((last-288):last))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(hi((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxhi1day.dat") day1string(q,11) = tostring(max(hi((last-288):last))) day1string(q,12) = datetime(maxind(hi((last-288):last))) end if if num(ismissing(hi((last-(2*288)):last))).lt.12 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(hi((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxhi2day.dat") end if if num(ismissing(hi((last-(3*288)):last))).lt.18 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(hi((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxhi3day.dat") end if if num(ismissing(hi((last-(5*288)):last))).lt.30 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(hi((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxhi5day.dat") end if if num(ismissing(hi)).lt.42 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(hi) system("echo "+quote+relative(0)+quote+" >> maxhi7day.dat") end if ; wind chill section if num(ismissing(wc((last-288):last))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(wc((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxwc1day.dat") day1string(q,13) = tostring(max(wc((last-288):last))) day1string(q,14) = datetime(maxind(wc((last-288):last))) end if if num(ismissing(wc((last-(2*288)):last))).lt.12 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(wc((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwc2day.dat") end if if num(ismissing(wc((last-(3*288)):last))).lt.18 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(wc((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwc3day.dat") end if if num(ismissing(wc((last-(5*288)):last))).lt.30 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(wc((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxwc5day.dat") end if if num(ismissing(wc)).lt.42 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+min(wc) system("echo "+quote+relative(0)+quote+" >> maxwc7day.dat") end if ; dew point section if num(ismissing(dewp((last-288):last))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(dewp((last-288):last)) system("echo "+quote+relative(0)+quote+" >> maxdewp1day.dat") day1string(q,15) = tostring(max(dewp((last-288):last))) day1string(q,16) = datetime(maxind(dewp((last-288):last))) end if if num(ismissing(dewp((last-(2*288)):last))).lt.12 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(dewp((last-(2*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxdewp2day.dat") end if if num(ismissing(dewp((last-(3*288)):last))).lt.18 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(dewp((last-(3*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxdewp3day.dat") end if if num(ismissing(dewp((last-(5*288)):last))).lt.30 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(dewp((last-(5*288)):last)) system("echo "+quote+relative(0)+quote+" >> maxdewp5day.dat") end if if num(ismissing(dewp)).lt.42 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+max(dewp) system("echo "+quote+relative(0)+quote+" >> maxdewp7day.dat") end if ; precip part if .not.any(ismissing(precip(last))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+precip(last) system("echo "+quote+relative(0)+quote+" >> precip7day.dat") if .not.any(ismissing(precip(last-288))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(precip(last)-precip(last-288)) system("echo "+quote+relative(0)+quote+" >> precip1day.dat") day1string(q,17) = tostring((precip(last)-precip(last-288))) system("sed -i.bak 's/"+meta(0)+"/"+meta(0)+"\&\#013;"+(.01*round(100*(precip(last)-precip(last-288)),3))+"\&\#34/' MRMS1dlocations.html") end if if .not.any(ismissing(precip(last-(2*288)))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(precip(last)-precip(last-(2*288))) system("echo "+quote+relative(0)+quote+" >> precip2day.dat") system("sed -i.bak 's/"+meta(0)+"/"+meta(0)+"\&\#013;"+(.01*round(100*(precip(last)-precip(last-(2*288))),3))+"\&\#34/' MRMS2dlocations.html") end if if .not.any(ismissing(precip(last-(3*288)))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(precip(last)-precip(last-(3*288))) system("echo "+quote+relative(0)+quote+" >> precip3day.dat") system("sed -i.bak 's/"+meta(0)+"/"+meta(0)+"\&\#013;"+(.01*round(100*(precip(last)-precip(last-(3*288))),3))+"\&\#34/' MRMS3dlocations.html") end if if .not.any(ismissing(precip(last-(5*288)))) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(precip(last)-precip(last-(5*288))) system("echo "+quote+relative(0)+quote+" >> precip5day.dat") end if end if ; snowfall part if num(ismissing(snowfall)).lt.8 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+accumulation(168-1) system("echo "+quote+relative(0)+quote+" >> snow7day.dat") if num(ismissing(snowfall(0:23))).lt.2 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(accumulation(23)) system("echo "+quote+relative(0)+quote+" >> snow1day.dat") day1string(q,18) = tostring(accumulation(23)) end if if num(ismissing(snowfall(0:47))).lt.3 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(accumulation(47)) system("echo "+quote+relative(0)+quote+" >> snow2day.dat") end if if num(ismissing(snowfall(0:71))).lt.4 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(accumulation(71)) system("echo "+quote+relative(0)+quote+" >> snow3day.dat") end if if num(ismissing(snowfall(0:119))).lt.6 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(accumulation(119)) system("echo "+quote+relative(0)+quote+" >> snow5day.dat") end if end if end do asciiwrite ("dayone.txt",day1string) end