One of the leaked files from the so-called Climategate scandal

briffa_sep98_e.pro

;
; PLOTS 'ALL' REGION MXD timeseries from age banded and from hugershoff
; standardised datasets.
; Reads Harry's regional timeseries and outputs the 1600-1992 portion
; with missing values set appropriately.  Uses mxd, and just the
; "all band" timeseries
;****** APPLIES A VERY ARTIFICIAL CORRECTION FOR DECLINE*********
;
yrloc=[1400,findgen(19)*5.+1904]
valadj=[0.,0.,0.,0.,0.,-0.1,-0.25,-0.3,0.,-0.1,0.3,0.8,1.2,1.7,2.5,2.6,2.6,$
  2.6,2.6,2.6]*0.75         ; fudge factor
if n_elements(yrloc) ne n_elements(valadj) then message,'Oooops!'
;
loadct,39
def_1color,20,color='red'
plot,[0,1]
multi_plot,nrow=4,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=18
endelse
;
; Get regional tree lists and rbar
;
restore,filename='reglists.idlsave'
harryfn=['nwcan','wnam','cecan','nweur','sweur','nsib','csib','tib',$
  'esib','allsites']
;
rawdat=fltarr(4,2000)
for i = nreg-1 , nreg-1 do begin
  fn='mxd.'+harryfn(i)+'.pa.mean.dat'
  print,fn
  openr,1,fn
  readf,1,rawdat
  close,1
  ;
  densadj=reform(rawdat(2:3,*))
  ml=where(densadj eq -99.999,nmiss)
  densadj(ml)=!values.f_nan
  ;
  x=reform(rawdat(0,*))
  kl=where((x ge 1400) and (x le 1992))
  x=x(kl)
  densall=densadj(1,kl)     ; all bands
  densadj=densadj(0,kl)     ; 2-6 bands
  ;
  ; Now normalise w.r.t. 1881-1960
  ;
  mknormal,densadj,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
  mknormal,densall,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd
;
; APPLY ARTIFICIAL CORRECTION
;
yearlyadj=interpol(valadj,yrloc,x)
densall=densall+yearlyadj
  ;
  ; Now plot them
  ;
  filter_cru,20,tsin=densall,tslow=tslow,/nan 
  cpl_barts,x,densall,title='Age-banded MXD from all sites',$ 
    xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$ 
    zeroline=tslow,yrange=[-7,3]
  oplot,x,tslow,thick=3 
  oplot,!x.crange,[0.,0.],linestyle=1
  ;
endfor
;
; Restore the Hugershoff NHD1 (see Nature paper 2)
;
xband=x
restore,filename='../tree5/densadj_MEAN.idlsave'
; gets: x,densadj,n,neff
;
; Extract the post 1600 part
;
kl=where(x ge 1400)
x=x(kl)
densadj=densadj(kl)
;
; APPLY ARTIFICIAL CORRECTION
;
yearlyadj=interpol(valadj,yrloc,x)
densadj=densadj+yearlyadj
;
; Now plot it too
;
filter_cru,20,tsin=densadj,tslow=tshug,/nan 
cpl_barts,x,densadj,title='Hugershoff-standardised MXD from all sites',$ 
  xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$ 
  zeroline=tshug,yrange=[-7,3],bar_color=20
oplot,x,tshug,thick=3,color=20
oplot,!x.crange,[0.,0.],linestyle=1
;
; Now overplot their bidecadal components
;
plot,xband,tslow,$
  xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$ 
  yrange=[-6,2],thick=3,title='Low-pass (20-yr) filtered comparison'
oplot,x,tshug,thick=3,color=20
oplot,!x.crange,[0.,0.],linestyle=1
;
; Now overplot their 50-yr components
;
filter_cru,50,tsin=densadj,tslow=tshug,/nan 
filter_cru,50,tsin=densall,tslow=tslow,/nan 
plot,xband,tslow,$
  xrange=[1399.5,1994.5],xtitle='Year',/xstyle,$ 
  yrange=[-6,2],thick=3,title='Low-pass (50-yr) filtered comparison'
oplot,x,tshug,thick=3,color=20
oplot,!x.crange,[0.,0.],linestyle=1
;
; Now compute the full, high and low pass correlations between the two
; series
;
perst=1400.
peren=1992.
;
openw,1,'corr_age2hug.out'
thalf=[10.,30.,50.,100.]
ntry=n_elements(thalf)
printf,1,'Correlations between timeseries'
printf,1,'Age-banded vs. Hugershoff-standardised'
printf,1,'     Region    Full   <10   >10   >30   >50  >100'
;
kla=where((xband ge perst) and (xband le peren))
klh=where((x ge perst) and (x le peren))
ts1=densadj(klh)
ts2=densall(kla)
;
r1=correlate(ts1,ts2)
rall=fltarr(ntry)
for i = 0 , ntry-1 do begin
  filter_cru,thalf(i),tsin=ts1,tslow=tslow1,tshigh=tshi1,/nan
  filter_cru,thalf(i),tsin=ts2,tslow=tslow2,tshigh=tshi2,/nan
  if i eq 0 then r2=correlate(tshi1,tshi2)
  rall(i)=correlate(tslow1,tslow2)
endfor
;
printf,1,'ALL SITES',r1,r2,rall,$
  format='(A11,2X,6F6.2)'
;
printf,1,' '
printf,1,'Correlations carried out over the period ',perst,peren
;
close,1
;
end