Still More on the CRU Data Dump...

Here's one more code file, this one calibrate_nhrecon.pro.  This is another IDL code file, and it appears to be calibrating tree ring data against direct temperature measurement data:

;
; Calibrates, usually via regression, various NH and quasi-NH records
; against NH or quasi-NH seasonal or annual temperatures.
;
; Specify period over which to compute the regressions (stop in 1960 to avoid
; the decline that affects tree-ring density records)
;
perst=1881.  ;note 1
peren=1960.
thalf=10.     ; filter to use to give extra hi & lo pass info
;
; Select season of the temperature data against which to calibrate
;
if n_elements(doseas) eq 0 then doseas=0      ; 0=annual, 1=Apr-Sep, 2=Oct-Mar
seasname=['annual','aprsep','octmar']
;
; Select record to calibrate
;
if n_elements(dorec) eq 0 then dorec=0
recname=['esper','jones','mann','3tree','briffa','iwarm','mj2000','crowley03','rutherford04']
print,recname[dorec]
doland=0      ; for Mann et al. only, 0=use NH mean, 1=use land>20N mean
;
if recname[dorec] eq 'mann' then begin
  if doland eq 0 then openw,1,'recon_mannNH.out'+seasname[doseas] $
                 else openw,1,'recon_mannLN20.out'+seasname[doseas]
endif else begin
  openw,1,'recon_'+recname[dorec]+'.out'+seasname[doseas]
endelse
;
multi_plot,nrow=2
if !d.name eq 'X' then window,ysize=800
;
; Compute the >20N land instrumental temperature timseries and filter
;
datst=1860
daten=1980
print,'Reading temperatures'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2_18512001.mon.nc')
tsmon=crunc_rddata(ncid,tst=[datst,0],tend=[daten,11],grid=gtemp)
ncdf_close,ncid
nmon=12
ntemp=gtemp.nt
nyrtemp=ntemp/nmon
yrtemp=reform(gtemp.year,nmon,nyrtemp)
yrtemp=reform(yrtemp(0,*))
;
; Compute the northern hemisphere >20N land series
;
; First extract >20N rows
kl=where(gtemp.y gt 20.)
ylat=gtemp.y(kl)
tsnorth=tsmon(*,kl,*)
; Compute latitude-weighted mean
nhmon=globalmean(tsnorth,ylat)
; Compute seasonal/annual mean
nhmon=reform(nhmon,nmon,nyrtemp)
case doseas of  ;note 2
  0: lvy=mkseason(nhmon,0,11,datathresh=6)   ; could try 9,8 (Oct-Sep annual)!
  1: lvy=mkseason(nhmon,3,8,datathresh=3)
  2: lvy=mkseason(nhmon,9,2,datathresh=3)
endcase
;
; Filter it
;
filter_cru,thalf,tsin=lvy,tslow=lvylow,tshigh=lvyhi,/nan
;
; Read in record and filter
;
case recname[dorec] of
  'esper': begin
    openr,2,'/cru/u2/f055/data/paleo/esper2002/esper.txt'
    readf,2,nyr
    headst=''
    readf,2,headst
    rawdat=fltarr(7,nyr)
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(1,*))
  end
  'jones': begin
    openr,2,'../tree5/phil_nhrecon.dat'
    nyr=992
    rawdat=fltarr(4,nyr)
    readf,2,rawdat,format='(I5,F7.2,I3,F7.2)'
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(3,*))
  end
  'mann': begin
    if doland eq 0 then begin
      openr,2,'../tree5/mann_nhrecon1000.dat'
      nyr=981
      rawdat=fltarr(2,nyr)
      readf,2,rawdat           ;,format='(I6,F11.7)'
      close,2
      x=reform(rawdat(0,*))
      densall=reform(rawdat(1,*))
    endif else begin
      openr,2,'../tree5/mannarea_all.dat'
      nyr=981
      rawdat=fltarr(11,nyr)
      headdat=' '
      readf,2,headdat
      readf,2,rawdat           ;,format='(I6,F11.7)'
      close,2
      x=reform(rawdat(0,*))
      densall=reform(rawdat(10,*))
    endelse
  end
  '3tree': begin
    openr,2,'../tree6/tornyamataim.ave'
    readf,2,nnn
    rawdat=fltarr(2,nnn)
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(1,*))
  end
  'briffa': begin
    restore,filename='/cru/u2/f055/tree6/bandtempNH_calmultipcr.idlsave'
    x=yrmxd
    densall=prednh
  end
  'iwarm': begin     ; use warm-season instrumental series as the predictor!
    x=yrtemp
    densall=mkseason(nhmon,3,8,datathresh=3)
  end
  'mj2000': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/mann03_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
  'crowley03': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/crowley03_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
  'rutherford04': begin
    openr,2,'/cru/u2/f055/data/paleo/ipccar4/data/rutherford04_orig.dat'
    readf,2,ncol
    readf,2,icol
    readf,2,nyr
    readf,2,nhead
    headst=strarr(nhead)
    rawdat=fltarr(ncol,nyr)
    readf,2,headst
    readf,2,rawdat
    close,2
    x=reform(rawdat(0,*))
    densall=reform(rawdat(icol,*))
  end
endcase
;
oopsx=x
oopsy=densall
;
kl=where((x ge datst) and (x le daten))
x=x(kl)
densall=densall(kl)
if total(abs(yrtemp-x)) ne 0 then message,'Incompatible years'
y1=densall
filter_cru,thalf,tsin=y1,tslow=ylow1,tshigh=yhi1,/nan
;
; Now correlate and regress them
;
printf,1,'Correlations and regression coefficients for '+seasname[doseas]
keeplist=where(finite(y1+lvy) and (x ge perst) and (x le peren),nkeep)
x=x(keeplist)
ts1=y1(keeplist)
ts2=lvy(keeplist)
r1=correlate(ts1,ts2)
dum=linfit(ts1,ts2,sigma=err_sigma)
c1=dum
printf,1,'Full timeseries:',r1,c1
plot,x,ts2
oplot,x,ts1*c1[1]+c1[0],thick=2
;
tsa=ts1(0:nkeep-2)
tsb=ts1(1:nkeep-1)
;mxd_ar1=correlate(tsa,tsb)
mxd_ar1=a_correlate(ts1,[-1])
tsa=ts2(0:nkeep-2)
tsb=ts2(1:nkeep-1)
;nht_ar1=correlate(tsa,tsb)
nht_ar1=a_correlate(ts2,[-1])
printf,1,'AR1 for MXD and NHEMI:',mxd_ar1,nht_ar1
;
tspred=c1(0)+c1(1)*ts1
tswant=ts2
plot,tspred,ts2,psym=1,$
  xtitle='Scaled Esper et al. anomaly  (!Uo!NC)',$
  ytitle='Northern Hemisphere temperature anomaly  (!Uo!NC)',$
  /xstyle,xrange=[-0.6,0.3],$
  /ystyle,yrange=[-0.6,0.3]
oplot,!x.crange,[0.,0.],linestyle=1
oplot,[0.,0.],!y.crange,linestyle=1
dum=linfit(tspred,ts2,sigma=se)
oplot,!x.crange,!x.crange*dum(1)+dum(0)
oplot,!x.crange,!x.crange*dum(1)+dum(0)+se(0)
oplot,!x.crange,!x.crange*dum(1)+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)+se(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)+se(1))+dum(0)-se(0)
oplot,!x.crange,!x.crange*(dum(1)-se(1))+dum(0)+se(0)
;
ts1=yhi1(keeplist)
ts2=lvyhi(keeplist)
r2=correlate(ts1,ts2)
dum=linfit(ts1,ts2)
c2=dum
printf,1,'High-pass      :',r2,c2
;
ts1=ylow1(keeplist)
ts2=lvylow(keeplist)
r3=correlate(ts1,ts2)
dum=linfit(ts1,ts2)
c3=dum
printf,1,'Low-pass       :',r3,c3
;
; Now compute the rms error between the reconstruction and the original
; timeseries
;
tserr=tswant-tspred
rmserr=sqrt( total( tserr^2 ) / float(n_elements(tserr)) )
printf,1,'RMS error between land>20Napr-sep temperature and Esper et al. reconstruction'
printf,1,rmserr
printf,1,'Uncertainties surrounding regression coefficients'
printf,1,err_sigma
;
printf,1,' '
printf,1,'Computations carried out over the period ',perst,peren
printf,1,' '
printf,1,'To separate low and high frequency components, a gaussian weighted'
printf,1,'filter was used with a half-width (years) of ',thalf
;
close,1
;
end

Where to start?  Here, I guess:

1.  Another arbitrary-looking choice of limits.  In this case, they've selected the period 1881 through 1960 to correlate the two data series.  The comment above, however, is very suggestive of why the 1960 cutoff was chosen.

2.  This piece of code looks like classic cherry-picking and fudging.  Here the software writer is choosing arbitrary months and data thresholds when sampling the annual mean of the instrumental data (which will be used to calibrate the tree-ring data.  The comment reinforces this notion.

Something that strikes me as I read through all this code: it looks amateurish, not at all like what I'd imagine code written by scientists expecting peer review would produce.  I'm no scientist, so it may well just reflect my naiveté about things scientific.  But I'm quite disappointed in it...