; 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" ;load "/ct12/abrammer/ncl6.3.0/lib/ncarg/nclscripts/csm/contributed.ncl" begin ;set history count here, for number of days in period: history = 365 ; 365 is the current maximum value. mindays = round((history*.94),3) ;; minimum number of days required to compute climate info ;; check to see if it's worth running whole script - no longer necessary ;recorddata = asciiread("records.ascii", -1, "string") ;sitedays = stringtofloat(recorddata(52)) ;if sitedays.lt.mindays ; exit() ;end if 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,:)=.1*(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 ;;;;;;; figure out what days to plot daystring=new((/(history+12)/), "string") yearstring=new((/(history+12)/), "string") daynumber=new((/(history+12)/), "integer") official=new((/4,(history+12)/), "float") do i=((history*-1)-1),10,1 daystring(i+(history+1)) = tostring(systemfunc("date +%j -d "+i+"days")) yearstring(i+(history+1))= tostring(systemfunc("date +%y -d "+i+"days")) daynumber(i+(history+1)) = tointeger(systemfunc("date +%j -d "+i+"days")) end do ;;;;;; determine how many days are in this year offset=365 ; normal year first = ((history*-1)-1) ; 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=((history*-1)+1),0,1 official(:,i+history) = asciiread(meta(1)+"/"+meta(1)+yearstring(i+history)+daystring(i+history)+".txt", -1, "float") end do ;;;;;; getting right climo data do i=0,(history),1 if daynumber(i).gt.daynumber((history)) 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((history+11))-1)) official@_FillValue = -996 overlap=new((/2,(history+12)/), "float") print(dimsizes(overlap)) print(dimsizes(window)) 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,:)) 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 precipwindow = precip(:,(offset+daynumber(0)-1):(offset+daynumber((history+11))-1)) precipreal = precipwindow if all(ismissing(official(2,:))) precipreal(1,:)=0 else precipreal(1,:)=cumsum(official(2,:),2) ; plot only real data end if h = official(0,ind(.not.ismissing(official(0,:)))) l = official(1,ind(.not.ismissing(official(0,:)))) ; find expected historical 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) ;w =.1*round(10*max(winds(1,:)),3) ;t =.01*round(100*max(cumsum(precip,2)),3) indices=ind(.not.ismissing(official(0,:))) ;write out historical anomaly, if it exists if dimsizes(indices).ge.round((history*.94),3) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(actualmean-expectedmean) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> specified.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+actualmean system("echo "+quote+relative(0)+quote+" >> specifiedmean.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(h) system("echo "+quote+relative(0)+quote+" >> specifiedhigh.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ avg(l) system("echo "+quote+relative(0)+quote+" >> specifiedlow.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ max(h) system("echo "+quote+relative(0)+quote+" >> specifiedmaximum.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ min(l) system("echo "+quote+relative(0)+quote+" >> specifiedminimum.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(h)-expectedhigh) system("echo "+quote+relative(0)+quote+" >> specifiedhighanomaly.dat") relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+ (avg(l)-expectedlow) system("echo "+quote+relative(0)+quote+" >> specifiedlowanomaly.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+" >> specifiedheating.dat") monthd = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+cooling system("echo "+quote+monthd(0)+quote+" >> specifiedcooling.dat") end if ;;;;; precip output ;write out precip historical anomaly and total, if it exists if dimsizes(pindices).ge.round((history*.94),3) 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+" >> precipspecified.dat") system("echo "+quote+mrelative(0)+quote+" >> preciptotalspecified.dat") end if ;;;;; write out solar data ; write out historical average, if it exists sindices=ind(.not.ismissing(official(3,:))) solar = ((sum(official(3,:))*300)/1000000)/dimsizes(sindices) if dimsizes(sindices).ge.round((history*.94),3) relative = ""+tofloat(meta(2))+","+tofloat(meta(3))+","+(solar) quote = inttochar(34) system("echo "+quote+relative(0)+quote+" >> solarspecified.dat") end if end