; 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/wrf/WRFUserARW.ncl" ;load "/home/torn/ncl/wrf_map_info.ncl" ;load "/home/torn/ncl/array_stat.ncl" ;load "/home/torn/ncl/image_conv.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/wrf/WRFUserARW.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" begin fasc = "../*.txt" ; a unique identifier for file DASC = "./" ; output dir FASC = "BIG.txt" ; output file name system ("/bin/rm -f "+DASC+FASC) ; rm any pre-existing file ncol = stringtoint(systemfunc("ls -1 ../*.txt | wc -l")) ; Use UNIX "cat" to concatenate the files into one file. system ("cat "+fasc+" > "+DASC+FASC) ; You can now read the file via "asciiread". data = asciiread(DASC+FASC,(/ncol,365/),"float") system ("/bin/rm -f "+DASC+FASC) ; rm any pre-existing file ;;;;;;;;;Here's the order of data ;sdtmax(0-3),sdtmin(4-7),tavg(8-11),tmax(12-15),tmin(16-19),ytdp(20-23) ; meta files with info meta = asciiread("station.ascii", -1, "string") metadata = stringtofloat(meta(2:5)) ;;;;;; read the water data, and calculate distance from site if needed delim = "," if metadata(3).eq.0 waterfile = "water.dat" waterdata = asciiread(waterfile,-1,"string") waterlats = tofloat(str_get_field(waterdata,3,delim)) waterlons = tofloat(str_get_field(waterdata,4,delim)) wateryes = tofloat(str_get_field(waterdata,5,delim)) distance = new((/dimsizes(waterlats)/),float) distance = 10000 do i = 0,dimsizes(waterlats)-1 distance(i) = gc_latlon(metadata(0),metadata(1),waterlats(i),waterlons(i),2,4) distance(i) = distance(i)*wateryes(i) end do distance@_FillValue = 0 print("Minimum distance to water for "+meta(1)+":"+min(distance)) else distance = metadata(3) end if ;;;;;;;;;Here's the order of data ;sdtmax(0-4),sdtmin(5-9),tavg(10-14),tmax(15-19),tmin(20-24),ytdp(25-28/9) temp=new((/8,365/), "float") temp(0,:)=data(15,:)+metadata(2)*data(16,:)+metadata(0)*data(17,:)+metadata(1)*data(18,:)+min(distance)*data(19,:) ; tmax temp(1,:)=temp(0,:)+(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:) ; tmax+ temp(2,:)=temp(0,:)-(data(0,:)+metadata(2)*data(1,:)+metadata(0)*data(2,:)+metadata(1)*data(3,:))+min(distance)*data(4,:) ; tmax- temp(3,:)=data(20,:)+metadata(2)*data(21,:)+metadata(0)*data(22,:)+metadata(1)*data(23,:)+min(distance)*data(24,:) ; tmin temp(4,:)=temp(3,:)+(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:) ; tmin+ temp(5,:)=temp(3,:)-(data(5,:)+metadata(2)*data(6,:)+metadata(0)*data(7,:)+metadata(1)*data(8,:))+min(distance)*data(9,:) ; tmin- temp(6,:)=data(10,:)+metadata(2)*data(11,:)+metadata(0)*data(12,:)+metadata(1)*data(13,:)+min(distance)*data(14,:) ; tavg temp(7,:)=data(25,:)+metadata(2)*data(26,:)+metadata(0)*data(27,:)+metadata(1)*data(28,:)+min(distance)*data(29,:) ; precip - use last part if you have fixed problem temp = temp*10 ;;;;;; figure out what days to plot daystring=new((/42/), "string") yearstring=new((/42/), "string") daynumber=new((/42/), "integer") official=new((/4,42/), "float") do i=-31,10,1 daystring(i+31) = tostring(systemfunc("date +%j -d "+i+"days")) yearstring(i+31)= tostring(systemfunc("date +%y -d "+i+"days")) daynumber(i+31) = tointeger(systemfunc("date +%j -d "+i+"days")) end do ;;;;;; determine how many days are in this year offset=365 ; normal year first = -31 ; set as minimum i above if mod(tointeger(systemfunc("date +%y -d "+first+"days")),4.0).eq.0 offset = 366 end if ;;;;;;reading in station data official(:,:) = -996 do i=-29,0,1 official(:,i+30) = asciiread(meta(1)+"/"+meta(1)+yearstring(i+30)+daystring(i+30)+".txt", -1, "float") end do ;;;;;; getting right climo data do i=0,41,1 if daynumber(i).gt.daynumber(41) daynumber(i)=daynumber(i)-offset end if end do ;;;;;; getting right climo data days = daynumber+offset temps = table_attach_columns(temp,temp,0) ; make two years worth of data to handle transition across year temps(7,365:729) = temps(7,365:729)+ temp(7,364) ; make sure precip is continuous window = temps(:,(days(0)-1):(days(41)-1)) official@_FillValue = -996 past24hour = systemfunc("date -u -d -24hours +%H") monplus1week = systemfunc("date -d 7days +%b") mon1week = systemfunc("date -d -7days +%b") mon2week = systemfunc("date -d -14days +%b") mon3week = systemfunc("date -d -21days +%b") mon4week = systemfunc("date -d -28days +%b") dayplus1week = systemfunc("date -d 7days +%d") day1week = systemfunc("date -d -7days +%d") day2week = systemfunc("date -d -14days +%d") day3week = systemfunc("date -d -21days +%d") day4week = systemfunc("date -d -28days +%d") wks_type = "png" wks_type@wkWidth = 2048 wks_type@wkHeight = 2048 wks = gsn_open_wks(wks_type,meta(1)) res = True res@gsnFrame = False res@tiMainFontHeightF = 0.018 res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tiXAxisString = "Day" res@tmXBTickSpacingF = 1. res@tmXBTickStartF = 0. res@tmXMajorGrid = True res@tmYLMinorPerMajor = 1 res@tmYMajorGrid = True res@trXMinF = days(0) res@trXMaxF = days(41) res@tmXBTickSpacingF = 5 res@tmXBMode = "Explicit" ; explicitly set X-axis labels ;res@tmXBValues = (/"1","32","60","91","121","152","182","213","244","274","305","335"/) ;res@tmXBLabels = (/"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan"/) res@tmXBValues = (/days(31)-28,days(31)-21,days(31)-14,days(31)-7,days(31)+7/) res@tmXBLabels = (/mon4week+" "+day4week,mon3week+" "+day3week,mon2week+" "+day2week,mon1week+" "+day1week,monplus1week+" "+dayplus1week,"","","","","","","","",""/) ;res@tmXBLabels = (/"-4 Weeks","-3 Weeks","-2 Weeks","-1 Week","+1 Week","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan"/) res@tmXBLabelFontHeightF = .015 res@tmYROn = False res@tiMainString = "Recent Climate Values For " + meta(0) ; add title res@tiYAxisString = "Average Daily High/Low/Mean (Red/Blue/Green) +/- 1 SD" res@xyLineColors = (/"firebrick3","firebrick3","firebrick3","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) ; change line color res@xyLineThicknesses = (/30.0,8.0,8.0,24.0,4.0,4.0,24.0,24.0/) res@tiYAxisFontHeightF = 0.015 res@trYMinF = min(official(1,:))-8; min(window)-30 res@trYMaxF = where((max(official(0,:))+5).gt.max(window(1,0:36)),(max(official(0,:))+5),max(window(1,0:36))); max(window)+20 res@tmYLPrecision = 2 ;res@xyMonoMarkLineMode= False ;res@xyMarkLineModes = (/"Lines","Lines","Lines","Markers","Markers","Markers","Markers"/) res@xyMarkLineModes = "Lines" ;res@xyMonoDashPattern = False ;res@xyDashPattern = (/"SolidLine","SolidLine","SolidLine",4,4,4,4/) res@xyDashPattern = (/0,0,0,4,4,4,4,4/) ;res@xyMarkers = 16 ; some dot ;res@xyMarkerColors = "firebrick3" ;res@xyMarkerSizes = 0.055 res@gsnXYFillColors = 180 ; "Light Orange" res@xyLineColor = -1 plot = gsn_csm_xy(wks,days,window(0:2,:),res); create plot res@xyLineColors = (/"dodgerblue","dodgerblue","dodgerblue","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) ; change line color res@gsnXYFillColors = 50 ; "pink1" res@xyLineColor = -1 plot = gsn_csm_xy(wks,days,window(3:5,:),res); create plot overlap=new((/2,42/), "float") overlap(0,:)=window(2,:) overlap(1,:)=window(4,:) overlap(0,:)=where(overlap(0,:).gt.overlap(1,:),overlap(1,:),overlap(0,:)) overlap(1,:)=where(overlap(0,:).gt.overlap(1,:),overlap(0,:),overlap(1,:)) res@xyLineThicknesses = (/0.001,0.001,1,24,4,4,24,24/) res@xyLineColors = (/"dodgerblue","firebrick3","dodgerblue","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) ; change line color res@gsnXYFillColors = 255 ; "pink1" res@xyLineColor = -1 plot = gsn_csm_xy(wks,days,overlap,res); create plot precip = new((/2,730/), "float") precip(1,:) = temps(7,0:729) do i=1,729,1 precip(0,i) = 0 precip(1,i) = temps(7,i)-temps(7,i-1) end do res@xyLineThicknesses = (/40.0,8.0,8.0,24.0,4.0,4.0,24.0,24.0/) res@xyLineColors = (/"green3","dodgerblue","dodgerblue","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) res@gsnXYFillColors = 255 ; "pink1" res@xyLineColor = -1 res@tiMainString = "" ; add title res2 = res res2@tmYUseLeft = False res2@xyLineThicknesses = (/0.01,0.01,8.0,24.0,4.0,4.0,24.0,24.0/) res2@xyLineColors = (/"black","black","dodgerblue","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) res2@gsnXYFillColors = 85 ; "pink1" res2@tiYAxisString = "Cumulative Precipitation (in), If Available" ;res2@tiYAxisOffsetYF = 0 ; -.15 ; moves y axis label res2@tmYMajorGrid = False res2@xyMarkLineModes = "Lines" ; choose which have markers res2@trYMinF = 0 res2@trYMaxF = 25 res2@tmYRLabelsOn = True ; turn on the YR tick mark labels res2@tmYRMode ="Explicit" res2@tmYRValues = (/"0","4","8","12"/) res2@tmYRLabels = (/"0","4","8","12"/) res2@tmYROn = True res2@tmYRLabelJust = "CenterLeft" precipwindow = precip(:,(offset+daynumber(0)-1):(offset+daynumber(41)-1)) precipreal = precipwindow if all(ismissing(official(2,:))) precipreal(1,:)=0 else precipreal(1,:)=cumsum(official(2,:),2) ; plot only real data end if plot = gsn_csm_xy2(wks,days,window(6,:),precipreal,res,res2) ; create plot delete(res) res = True res@gsnFrame = False res@tiMainFontHeightF = 0.018 res@tiXAxisFontHeightF = 0.015 res@tiYAxisFontHeightF = 0.015 res@tmBorderThicknessF = 40 res@tmXBTickSpacingF = 1. res@tmXBTickStartF = 0. res@tmXMajorGrid = True res@tmYLMinorPerMajor = 1 res@tmYMajorGrid = False res@trXMinF = days(0) res@trXMaxF = days(41) res@tmXBTickSpacingF = 5 res@tmXBMode = "Explicit" ; explicitly set X-axis labels res@tmXBValues = days(31) res@tmXBLabels = "Today" res@tmXBLabelFontHeightF = .015 res@tmYROn = False res@xyLineColors = (/"black","black","firebrick3","dodgerblue","dodgerblue","dodgerblue","green3","green3"/) ; change line color res@xyLineThicknesses = (/15.0,15.0,8.0,24.0,4.0,4.0,24.0,24.0/) ; this seems to control the bar chart too! res@tiYAxisFontHeightF = 0.015 res@trYMinF = min(official(1,:))-8; min(window)-30 res@trYMaxF = where((max(official(0,:))+5).gt.max(window(1,0:36)),(max(official(0,:))+5),max(window(1,0:36))); max(window)+20 res@tmYLPrecision = 2 res@tmXMajorGridThicknessF = 10 resbar = res resbar@gsnXYBarChart = True resbar@gsnXYBarChartBarWidth = .04 ; Default is 1.0 resbar@gsnYRefLineColor = -1 ; reference line resbar@gsnXYBarChart = True ; create bar chart resbar@gsnAboveYRefLineColor = "black" ; above ref line fill red resbar@gsnBelowYRefLineColor = "black" ; below ref line fill blue plot1 = new((/2,days(41)/),graphic) ; get indices of value data, and only plot those indices=ind(.not.ismissing(official(0,:))) do xx = 0,(dimsizes(indices)-1) gg = indices(xx) resbar@gsnXYBarChartColors = "black" resbar@gsnYRefLine = (official(0,gg)+official(1,gg))/2 plot1(0,gg) = gsn_csm_xy(wks,days(gg),official(0,gg),resbar) ; draw each individual max bar plot1(1,gg) = gsn_csm_xy(wks,days(gg),official(1,gg),resbar) ; draw each individual min bar end do res@tmBorderThicknessF = 0 res@xyMarkers = 16 ; choose type of marker res@xyMarkerColor = "Black" ; Marker color res@xyMarkerSizeF = 0.015 res@xyMarkLineMode = "Markers" ; choose to use markers plot = gsn_csm_xy(wks,days,official(0,:),res) ; create plot plot = gsn_csm_xy(wks,days,official(1,:),res) ; create plot res@xyMarkerColor = "gray30" plot = gsn_csm_xy(wks,(/days(31),days(31)/),(/window(3,31),window(0,31)/),res) ; current day's dot tres = True ; text mods desired ;tres@tfPolyDrawOrder = "PostDraw" tres@txFont = "helvetica-bold" tres@txFontHeightF = 0.014 tres@txFontColor = "gray30" gsn_text(wks,plot,(/"Avg:~C~"+.1*round(10*window(3,31),3)+"~S~o~N~F","Avg:~C~"+.1*round(10*window(0,31),3)+"~S~o~N~F","1 SD:~C~~F34~1~F22~"+.1*round(10*(window(1,31)-window(0,31)),3)+"~S~o~N~F","1 SD:~C~~F34~1~F22~"+.1*round(10*(window(4,31)-window(3,31)),3)+"~S~o~N~F"/),(/days(33),days(33),days(33),days(33)/)+.5,(/window(3,31)-2,window(0,31)-2,window(0,31)+4,window(3,31)+4/),tres) ;stationtext = gsn_add_text(wks,plot,(/"T","L"/),(/days(32),days(32)/),(/window(3,31)-2,window(0,31)-2/),tres) ;minmarker = gsn_add_polymarker(wks,plot,days(31),window(3,31),res) ;plot = gsn_csm_xy(wks,days,(official(0,:)+official(1,:))/2.0,res) ; unflag to put in average dot res@xyMarkerColor = "firebrick3" ; Marker color res@xyMarkerSizeF = 0.0075 ; Marker size (default 0.01) plot = gsn_csm_xy(wks,days,official(0,:),res) ; create plot res@xyMarkerColor = "dodgerblue" ; Marker color plot = gsn_csm_xy(wks,days,official(1,:),res) ; create plot ;res@xyMarkerColor = "green3" ; Marker color ; unflag to put in average dot ;plot = gsn_csm_xy(wks,days,(official(0,:)+official(1,:))/2.0,res) ; unflag to put in average dot ; print extremes, etc. txres = True txres@txFontHeightF = .012 ; Set the font height txres@txJust = "CenterLeft" label = "Location: " + meta(0) gsn_text_ndc(wks,label,.09,.12 ,txres) label = "Station elevation: " + .01*round(100*tofloat(meta(4)),3) + " m" gsn_text_ndc(wks,label,.09,.10 ,txres) label = "Station lat/lon: " + .01*round(100*tofloat(meta(2)),3) + "; " + .01*round(100*tofloat(meta(3)),3) gsn_text_ndc(wks,label,.09,.08 ,txres) label = "Figure created on " + systemfunc("date -u +%a") + ", " + systemfunc("date -u +%F") + ", at "+systemfunc("date -u +%R") + " UTC" gsn_text_ndc(wks,label,.09,.06 ,txres) label = "Climate Normals Estimated From 1990-2020" gsn_text_ndc(wks,label,.09,.04 ,txres) ; get actual available highs/lows h = official(0,ind(.not.ismissing(official(0,:)))) l = official(1,ind(.not.ismissing(official(0,:)))) ; find expected monthly precip, but first check to make sure they exist pindices=ind(.not.ismissing(official(2,:))) if dimsizes(pindices).gt.1 pexpvalid = new((/dimsizes(pindices)/), "float") pvalid = new((/dimsizes(pindices)/), "float") do xx = 0,(dimsizes(pindices)-1) pexpvalid(xx) = precipwindow(1,pindices(xx)) pvalid(xx) = official(2,pindices(xx)) end do pflag = 0 else pexpvalid = 0.00 pvalid = 0.00 pflag = 1 end if p =sum(pvalid)-sum(pexpvalid) p =.01*round((100*p(0)),3) ; get various means expectedhigh = avg(window(0,ind(.not.ismissing(official(0,:))))) expectedlow = avg(window(3,ind(.not.ismissing(official(0,:))))) expectedmean = avg(window(6,ind(.not.ismissing(official(0,:))))) actualmean = (avg(h)+avg(l))/2 d = .1*round(10*(actualmean - expectedmean),3) ;t =.01*round(100*max(cumsum(precip,2)),3) txres = True txres@txFontHeightF = .015 ; Set the font height txres@txJust = "CenterRight" label = "Summary Of Past "+(dimsizes(indices))+" Days" gsn_text_ndc(wks,label,.9,.12 ,txres) label = "Warmest/Coldest:" gsn_text_ndc(wks,label,.8,.095 ,txres) label = "Versus Expected Temperature:" gsn_text_ndc(wks,label,.835,.07 ,txres) label = "Versus Expected Precipitation:" gsn_text_ndc(wks,label,.835,.045 ,txres) txres@txFont = "helvetica-bold" label = .1*round(10*max(official(0,:)),3)+"/"+ .1*round(10*min(official(1,:)),3) gsn_text_ndc(wks,label,.9,.095 ,txres) txres@txFont = "helvetica-bold" txres@txFontColor = "dodgerblue" text = "" if d.ge.0 text = "+" txres@txFontColor = "firebrick3" end if label = text(0)+d gsn_text_ndc(wks,label,.9,.07 ,txres) text = "" if pflag.eq.0 txres@txFontColor = "sienna" if p.ge.0 text = "+" txres@txFontColor = "green3" end if label = text(0)+p gsn_text_ndc(wks,label,.9,.046 ,txres) else text = "N/A" txres@txFontColor = "black" label = text(0) gsn_text_ndc(wks,label,.9,.046 ,txres) end if txres@txFontColor = "black" txres@txJust = "CenterLeft" label = ""+(dimsizes(pindices)) gsn_text_ndc(wks,label,.66,.24 ,txres) ;gsn_text_ndc(wks,"Precip",.66,.24 ,txres) if (dimsizes(pindices)-1).eq.1 gsn_text_ndc(wks,"Day",.66,.215 ,txres) else gsn_text_ndc(wks,"Days",.66,.215 ,txres) end if ;;;;;; temperature anomaly output ; write out daily anomaly, if it exists if .not.any(ismissing(official(0,30))) dayh = official(0,30) dayl = official(1,30) dayex= window(6,30) dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(((dayh+dayl)/2)-dayex) quote = inttochar(34) system("echo "+quote+dayd(0)+quote+" >> daily.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+((dayh+dayl)/2) system("echo "+quote+dayd(0)+quote+" >> dailymean.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+dayh system("echo "+quote+dayd(0)+quote+" >> dailyhigh.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+dayl system("echo "+quote+dayd(0)+quote+" >> dailylow.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(dayh-window(0,30)) system("echo "+quote+dayd(0)+quote+" >> dailyhighanomaly.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(dayl-window(3,30)) system("echo "+quote+dayd(0)+quote+" >> dailylowanomaly.dat") highsigma = window(1,30)-window(0,30) dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+((dayh-window(0,30))/highsigma) system("echo "+quote+dayd(0)+quote+" >> dailyhighsigmaanomaly.dat") lowsigma = window(4,30)-window(3,30) dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+((dayl-window(3,30))/lowsigma) system("echo "+quote+dayd(0)+quote+" >> dailylowsigmaanomaly.dat") heating = 0.0 cooling = 0.0 if ((dayh+dayl)/2).gt.65 cooling = ((dayh+dayl)/2) - 65 else if ((dayh+dayl)/2).lt.65 heating = ((dayh+dayl)/2) - 65 end if end if dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+heating*-1 system("echo "+quote+dayd(0)+quote+" >> dailyheating.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+cooling system("echo "+quote+dayd(0)+quote+" >> dailycooling.dat") end if ; write out weekly anomaly, if it exists if .not.any(ismissing(official(0,24:30))) weekh = official(0,24:30) weekl = official(1,24:30) weekex= window(6,24:30) weekd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(((avg(weekh)+avg(weekl))/2)-avg(weekex)) quote = inttochar(34) system("echo "+quote+weekd(0)+quote+" >> weekly.dat") weekd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+((avg(weekh)+avg(weekl))/2) system("echo "+quote+weekd(0)+quote+" >> weeklymean.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+avg(weekh) system("echo "+quote+dayd(0)+quote+" >> weeklyhigh.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+avg(weekl) system("echo "+quote+dayd(0)+quote+" >> weeklylow.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(avg(weekh)-avg(window(0,24:30))) system("echo "+quote+dayd(0)+quote+" >> weeklyhighanomaly.dat") dayd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(avg(weekl)-avg(window(3,24:30))) system("echo "+quote+dayd(0)+quote+" >> weeklylowanomaly.dat") weekhc = (weekh + weekl)/2 heating = 0.0 cooling = 0.0 do i=0,6 if (weekhc(i).gt.65) cooling=cooling+(weekhc(i)-65) else if (weekhc(i).lt.65) heating=heating+(weekhc(i)-65) end if end if end do weekd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+heating*-1 system("echo "+quote+weekd(0)+quote+" >> weeklyheating.dat") weekd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+cooling system("echo "+quote+weekd(0)+quote+" >> weeklycooling.dat") end if ;write out monthly anomaly, if it exists if dimsizes(indices).gt.27 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(actualmean-expectedmean) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> monthly.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+actualmean system("echo "+quote+relative(0)+quote+" >> monthlymean.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(h) system("echo "+quote+relative(0)+quote+" >> monthlyhigh.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(l) system("echo "+quote+relative(0)+quote+" >> monthlylow.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(h)-expectedhigh) system("echo "+quote+relative(0)+quote+" >> monthlyhighanomaly.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(l)-expectedlow) system("echo "+quote+relative(0)+quote+" >> monthlylowanomaly.dat") monthhc = (h + l)/2 heating = 0.0 cooling = 0.0 do i=0,(dimsizes(indices)-1) if (monthhc(i).gt.65) cooling=cooling+(monthhc(i)-65) else if (monthhc(i).lt.65) heating=heating+(monthhc(i)-65) end if end if end do monthd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+heating*-1 system("echo "+quote+monthd(0)+quote+" >> monthlyheating.dat") monthd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+cooling system("echo "+quote+monthd(0)+quote+" >> monthlycooling.dat") end if ;;;;; precip output ;write out precip monthly anomaly and total, if it exists if dimsizes(pindices).gt.27 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(p) mrelative =""+tofloat(meta(2))+","+tofloat(meta(3))+","+ sum(pvalid) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> precipmonthly.dat") system("echo "+quote+mrelative(0)+quote+" >> preciptotalmonthly.dat") end if ;write out precip weekly anomaly and total, if it exists if .not.any(ismissing(official(2,24:30))) pweeklyvalid = sum(official(2,24:30)) pexweekly = sum(precipwindow(1,24:30)) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(pweeklyvalid-pexweekly) mrelative =""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (pweeklyvalid) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> precipweekly.dat") system("echo "+quote+mrelative(0)+quote+" >> preciptotalweekly.dat") end if if .not.any(ismissing(official(2,30))) pdailyvalid = sum(official(2,30)) pdailyex = sum(precipwindow(1,30)) mrelative =""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (pdailyvalid) quote = inttochar(34) system("echo "+quote+mrelative(0)+quote+" >> preciptotaldaily.dat") pdailyanom = pdailyvalid-pdailyex mrelative =""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (pdailyanom) system("echo "+quote+mrelative(0)+quote+" >> precipdaily.dat") end if ;;;;; write out solar data ; write out monthly average, if it exists sindices=ind(.not.ismissing(official(3,:))) solar = ((sum(official(3,:))*300)/1000000)/dimsizes(sindices) if dimsizes(sindices).gt.27 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(solar) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> solarmonthly.dat") end if ; write out weekly average, if it exists if .not.any(ismissing(official(3,24:30))) solar = ((sum(official(3,24:30))*300)/1000000)/7 relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(solar) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> solarweekly.dat") end if ; write out daily total, if it exists if .not.any(ismissing(official(3,30))) solar = ((sum(official(3,30))*300)/1000000) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(solar) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> solardaily.dat") end if frame(wks) end