; ; Reads in some long series (single tree chronology, regional tree ; timeseries, or palaeoaverage), and compute the distribution of ; trends of various lengths from the series, to find the 95% threshold ; of trend magnitudes. Only trends that start from near-normal values ; are used. ; ; Lengths of trends to consider ; tlen=[6,8,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100] ntry=n_elements(tlen) maxsize=min(tlen) ; ; Prepare for plotting ; loadct,39 def_1color,20,color='mgrey' multi_plot,nrow=3,ncol=2,layout='large' if !d.name eq 'X' then begin window, ysize=800 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=14 endelse cpl_usersym,/circle,/fill ; ; Now get the timeseries to use ; for i = 0 , 2 do begin case i of 0: begin ; Phil's reconstruction mtit="Jones' Northern Hemisphere temperature reconstruction" ; Period to consider perst=1400 peren=1980 openr,1,'phil_nhrecon.dat' nyr=992 rawdat=fltarr(4,nyr) readf,1,rawdat,format='(I5,F7.2,I3,F7.2)' close,1 timey=reform(rawdat(0,*)) ts=reform(rawdat(3,*)) kl=where((timey ge perst) and (timey le peren),nyr) timey=timey(kl) ts=ts(kl) yr=[0.,0.9] end 1: begin ; Mike Mann's reconstruction mtit="Mann's Northern Hemisphere temperature reconstruction" ; Period to consider perst=1400 peren=1980 openr,1,'mann_nhrecon.dat' nyr=596 headdat=' ' rawdat=fltarr(2,nyr) readf,1,headdat readf,1,rawdat,format='(I6,F11.7)' close,1 timey=reform(rawdat(0,*)) ts=reform(rawdat(1,*)) kl=where((timey ge perst) and (timey le peren),nyr) timey=timey(kl) ts=ts(kl) yr=[0.,0.5] end 2: begin ; NHD1 mtit="NHD1 growing-season temperature reconstruction" ; Period to consider perst=1400 peren=1960 restore,filename='densadj_MEAN.idlsave' timey=x ts=densadj*0.116626 ; converts it from density to temperature anom kl=where((timey ge perst) and (timey le peren),nyr) timey=timey(kl) ts=ts(kl) yr=[0.,0.5] end 3: begin ; Jasper mtit="Jasper temperature reconstruction" ; Period to consider ALL ;perst=1400 ;peren=1980 openr,1,'jasper.dat' readf,1,nyr,format='(1X,I4)' rawdat=fltarr(2,nyr) readf,1,rawdat,format='(1X,I5,F8.2)' close,1 timey=reform(rawdat(0,*)) ts=reform(rawdat(1,*)) ; kl=where((timey ge perst) and (timey le peren),nyr) ; timey=timey(kl) ; ts=ts(kl) yr=[0.,1.5] end 4: begin ; Polar Urals mtit="Polar Urals temperature reconstruction" ; Period to consider ALL ;perst=1400 ;peren=1980 openr,1,'polarural.dat' readf,1,nyr,format='(1X,I4)' rawdat=fltarr(2,nyr) readf,1,rawdat,format='(1X,I5,F8.2)' close,1 timey=reform(rawdat(0,*)) ts=reform(rawdat(1,*)) ; kl=where((timey ge perst) and (timey le peren),nyr) ; timey=timey(kl) ; ts=ts(kl) yr=[0.,1.5] end 5: begin ; Tornetrask mtit="Tornetrask temperature reconstruction" ; Period to consider ALL ;perst=1400 ;peren=1980 openr,1,'tornetrask.dat' readf,1,nyr,format='(I4)' rawdat=fltarr(2,nyr) readf,1,rawdat,format='(I5,F8.2)' close,1 timey=reform(rawdat(0,*)) ts=reform(rawdat(1,*)) ; kl=where((timey ge perst) and (timey le peren),nyr) ; timey=timey(kl) ; ts=ts(kl) yr=[0.,1.5] end endcase ; ; Compute the trend distributions ; ; Only compute them if first year in window is within +-0.5 sd from the ; post-1900 mean ; meanl=where(timey ge 1900,nmean) mean20c=total(ts(meanl))/float(nmean) dummy=moment(ts,sdev=tssd) stmin=mean20c-0.5*tssd stmax=mean20c+0.5*tssd ; tskeep=ts timeykeep=timey nyrkeep=nyr for kkk = 0 , 1 do begin ; ; For kkk=0, do all trends but stop in 1900 (and also compute trend ; since 1901) ; For kkk=1, do only trends starting from 'near-normal' conditions ; if kkk eq 0 then begin kl=where(timeykeep le 1900,nyr) ts=tskeep(kl) timey=timeykeep(kl) kl=where(timeykeep ge 1901) yval=tskeep(kl) xval=timeykeep(kl) kl=where(finite(yval),n19) xval=xval(kl) yval=yval(kl) acoeff=linfit(xval,yval) tren19=acoeff(1)*float(n19) endif else begin ts=tskeep timey=timeykeep nyr=nyrkeep endelse ; tren95=fltarr(ntry) ; for j = 0 , ntry-1 do begin print,j,tlen(j) ; ; ntren=nyr-tlen(j)+1 obstren=fltarr(ntren)*!values.f_nan ; for k = 0 , ntren-1 do begin xval=timey(k:k+tlen(j)-1) yval=ts(k:k+tlen(j)-1) kl=where(finite(yval),nkeep) if nkeep gt 0.66*tlen(j) then begin xval=xval(kl) yval=yval(kl) if ((yval(0) ge stmin) and (yval(0) le stmax)) or (kkk eq 0) then begin acoeff=linfit(xval,yval) obstren(k)=acoeff(1)*float(tlen(j)) endif endif endfor kl=where(finite(obstren),ntren) obstren=obstren(kl) print,tlen(j),ntren ; ; Compute the 95% threshold (easy cos there's no missing data) ; sortobs=obstren(sort(obstren)) i95=ntren*0.95 tren95(j)=sortobs(i95) ; endfor ; if kkk eq 0 then begin plot,tlen,tren95,title=mtit,$ /ystyle,yrange=yr,$ xtitle='Trend length (yr)',$ ytitle='95 centile of trends (!Uo!NC per trend length)',$ psym=1 oplot,tlen,tren95,thick=3 plots,n19,tren19,psym=4,symsize=2 endif else begin oplot,tlen,tren95,psym=8,color=20 oplot,tlen,tren95,color=20,thick=3 endelse ; endfor ; endfor ; end