;
; Plots our NH series that is a combination of high-freq Hugershoff and
; low-freq ABD, plus labels extreme cool summers and plots the VEI volcanoes
;
thalf=30
doland=0         ; 0=portrait, 1=landscape
doremove=1       ; 0=do not, 1=do remove a pre-defined list of years from the
                 ; MXD time series prior to computing the low-pass series.
;
loadct,39
def_1color,50,color='red'
def_1color,51,color='deepblue'
def_1color,52,color='green'
def_1color,53,color='dgreen'
def_1color,54,color='orange'
if doland eq 0 then multi_plot,nrow=3,layout='large' $
               else multi_plot,nrow=1,/landscape
if !d.name eq 'X' then begin
  if doland eq 0 then wintim,ysize=700 else wintim,xsize=800
  !p.font=-1
  initx
endif else begin
  !p.font=0
  if doland eq 0 then device,/helvetica,/bold,font_size=16,xsize=17.9 $
                 else device,/helvetica,/bold,font_size=16,ysize=12.,$
                             xoffset=6.
endelse
;
; Get reconstruction timeseries
;
restore,filename='combtempNH_calibrated.idlsave'
; Gets: nhtit,nyr,yrmxd,prednh
meanmxd=prednh
x=yrmxd
;
; Now plot them
;
if doremove eq 0 then begin
  ; Simply low-pass filter all the NH series
  filter_cru,thalf,tsin=meanmxd,tslow=tslow,/nan 
endif else begin
  ; Remove "volcano-influenced" years prior to computing low-pass filtered
  ; series, then interpolate to get values for those years
  mxd2=meanmxd
  korig=where(finite(mxd2))
  remlist=[1453,1587,1601,1641,1643,1816,1817,1884,1912,1992]
  nrem=n_elements(remlist)
  print,'REMOVING SOME YEARS FROM LOW-PASS FILTER:',nrem
  for irem = 0 , nrem-1 do begin
    krem=where(x eq remlist[irem],ngot)
    if irem eq 0 then listi=krem[0] else listi=[listi,krem[0]]
  endfor
  mxd2[listi]=!values.f_nan
  filter_cru,thalf,tsin=mxd2,tslow=tslow2,/nan 
  kl=where(finite(tslow2))
  xkeep=x[kl]
  ykeep=tslow2[kl]
  tslow3=interpol(ykeep,xkeep,x[korig])
  tslow=meanmxd*!values.f_nan
  tslow[korig]=tslow3[*]
endelse
;
kkk=where(x le 1960)
;
cpl_barts,x[kkk],meanmxd[kkk],$
  xrange=[1400,2000],xtitle='Year',/xstyle,$ 
  ytitle='Temperature averaged over land north of 20!Uo!NN!C(differences in !Uo!NC from the 1961-90 average)',$
  zeroline=tslow[kkk],yrange=[-2,0.3],ystyle=1,bar_color=[50,51]
oplot,x[kkk],tslow[kkk],thick=5
oplot,!x.crange,[0.,0.],linestyle=1
;
; Now let's list the warmest & coolest dozen
;
kl=where(finite(meanmxd),nk)
kmxd=meanmxd[kl]
ktim=x[kl]
isort=sort(kmxd)
jsort=reverse(isort)
print,'Coldest and warmest values in the combined ABD/HUG quasi-NH'
print,'Year Temp(K)      Year Temp(K)    (anomalies wrt 1961-90 mean)'
for i = 0 , 11 do begin
  print,ktim[isort[i]],kmxd[isort[i]],ktim[jsort[i]],kmxd[jsort[i]],$
    format='(I4,F8.2,6X,I4,F8.2)'
endfor
;
axis,xaxis=1,/xstyle           ;,xtitle='Year'
;
; Label some individual year spikes
;
yrlab=['1420','1453','1587','1601','1641/2/3','1666',$
  '1675','1695','1698/9','1742','1783','1816/7','1837','1884','1912']
yry  =[-3.00 ,-3.00 ,-3.03 ,-2.63 ,-3.32 ,-6.0  ,-5.       ,-3.42 ,$
  -3.81 ,-4.15 ,-3.52   ,-3.35 ,-3.05 ,-5.03     ,-3.5    ,-3.6  ,-4.1]
yry=yry/3.
yrx  =[1420  ,1453  ,1587  ,1589  ,1642      ,1658  ,$
  1675  ,1695  ,1711    ,1742  ,1783  ,1816    ,1837  ,1884  ,1912]
nlab=n_elements(yrlab)
for i = 0 , nlab-1 do xyouts,yrx(i),yry(i),/data,yrlab(i),size=0.5,align=0.5
;
; Make VEI volcanic events
;
keepcol=!p.color
!p.color=53
yrvei=[1450,1452,1471,1477,1480,1482,1580,1586,1593,1600,1640,1641,1660,$
  1663,1667,1673,1680,1707,1739,1800,1815,1835,1853,1854,1883,1886,1902,$
  1907,1912,1932,1956,1980,1982,1991]
yrsiz=[5,6,5,5,5,5,6,5,5,6,5,6,6,5,5,5,5,5,5,5,7,5,5,5,6,5,6,5,6,5,5,5,5,6]
nvei=n_elements(yrvei)
ym=!y.crange(0)
ycr=!y.crange(1)-!y.crange(0)
f03=ycr*0.3/12.
f04=ycr*0.4/12.
f05=ycr*0.5/12.
f07=ycr*0.7/12.
f11=ycr*1.1/12.
f14=ycr*1.4/12.
f16=ycr*1.6/12.
f22=ycr*2.2/12.
f25=ycr*2.5/12.
f27=ycr*2.7/12.
for i = 0 , nvei-1 do begin
  z=yrvei(i)
  y=yrsiz(i)
  case y of
    5: plots,[z,z,z-1,z,z+1,z],[ym,ym+f04,ym+f04,ym+f07,ym+f04,ym+f04]
    6: begin
      plots,[z,z],[ym,ym+f11],thick=3
      polyfill,[z-2.,z,z+2.,z-2.],[ym+f11,ym+f16,ym+f11,ym+f11]
      end
    7: begin
      plots,[z,z],[ym,ym+f22],thick=5
      polyfill,[z-3.,z,z+3.,z-3.],[ym+f22,ym+f27,ym+f22,ym+f22]
      end
  endcase
endfor
;
xyouts,1449.,ym+f03,align=1.,'?',size=0.5
xyouts,1579.,ym+f03,align=1.,'?',size=0.5
xyouts,1659.,ym+f03,align=1.,'?',size=0.5
;
xyouts,1420.,ym+f03,orient=90.,'VEI',size=0.75
plots,[1432.,1432.],[ym,ym+f27],thick=3
plots,[1432.,1438.],[ym+f27,ym+f27],thick=3
plots,[1432.,!x.crange(1)],[ym+f27,ym+f27],linestyle=1
plots,[1432.,1438.],[ym+f16,ym+f16],thick=3
plots,[1432.,!x.crange(1)],[ym+f16,ym+f16],linestyle=1
plots,[1432.,1438.],[ym+f07,ym+f07],thick=3
plots,[1432.,!x.crange(1)],[ym+f07,ym+f07],linestyle=1
xyouts,1431.,ym+f25,'7',size=0.6,align=1.
xyouts,1431.,ym+f14,'6',size=0.6,align=1.
xyouts,1431.,ym+f05,'5',size=0.6,align=1.
!p.color=keepcol
;
; Output the reconstruction time series
;
openw,1,'gpc2003_fig5.dat'
printf,1,'Data file to accompany Figure 5 of Briffa, Osborn and Schweingruber'
printf,1,'(2003) Large-scale temperature inferences from tree rings: a review.'
printf,1,'Global and Planetary Change 40, 11-26.
printf,1
printf,1,'These data are reconstructions of the average April-September'
printf,1,'temperature from all land data north of 20N, expressed as'
printf,1,'anomalies in degrees C with respect to the 1961-1990 mean.'
printf,1,'The series is a combination of the 25-year low-pass filtered'
printf,1,'reconstruction of Briffa et al. (2001) with the 25-year high-pass'
printf,1,'filtered reconstruction of Briffa et al. (2002).'
printf,1
printf,1,'Year  Anomaly'
for iyr = 0 , nyr-1 do begin
  if finite(meanmxd[iyr]) then printf,1,x[iyr],meanmxd[iyr],format='(I4,F10.2)'
endfor
close,1
;
; Finally, overplot the instrumental temperatures
;
obsts=crutempave(year=timey,latrange=[20,90],/land)
annts=mkseason(obsts,3,8)
filter_cru,thalf,/nan,tsin=annts,tslow=lowts
kkk=where(timey ge 1881)
oplot,timey[kkk],lowts[kkk],thick=8,color=54
;
end