;
; Plots long reconstructions (Age-banded, Jones, Mann(v2) and
; Yamal+Torne+Taimyr), and correlates them too!
;
;ploton,1
;
doseas=1
if n_elements(doseas) eq 0 then doseas=0  ; 0=annual, 1=apr-sep, 2=oct-mar
seasname=['annual','aprsep','octmar']
mtit=['ANNUAL (Jan-Dec)','WARM SEASON (Apr-Sep)','COLD SEASON (Oct-Mar)']
dobore=1  ; 0=no boreholes, 1=overlay borehole-based NH estimate
doerr=0   ; 0=no errors, 1=Briffa age errors
dost=850 ; start year
docorr=0        ; 1=compute correlations  0=don't
thalf=30.       ; filter for high or low-pass
dobigy=0        ; 0=compact yrange  1=big yrange
if doseas ne 1 then dobigy=1.4
doinstrwarm=0   ; 1=use warm-season instrumental series as predictor too!
;
; Have to restore the NEW age-banded results now to avoid overwriting key
; variables later
;
; For regional calibrated stuff
;restore,filename='bandtemp_calibrated.idlsave'
; Gets: calregts,nreg,nyr,regname,tempnyr,temptimey,tempregts,timey
;if regname(nreg-1) ne 'NH' then message,'Re-calibrate the age-banded!'
;newagetime=timey
;newagets=calregts(*,nreg-1)
;
; For allsites calibrated stuff
;restore,filename='bandalltemp_calibrated.idlsave'
; Gets: nyr,nhtit,timey,nhtemp
;if nhtit ne 'NH' then message,'Re-calibrate the age-banded!'
;newagetime=timey
;newagets=nhtemp
;
; For PCR calibrated stuff. The sm50 gets the smoothed series and appropriate
; erros - no need to smooth them!!!!  Except we want correlations between
; unsmoothed series, and between high-passed series, so we need to unsmoothed
; time series but the smoothed error bands!!
;
restore,filename='bandtempNH_calmultipcr.idlsave'
newagetime=yrmxd
newagets=prednh
restore,filename='bandtempNHsm50_calmultipcr.idlsave'
newagelow=prednh
newagese=predse
;restore,filename='bandtempNHsm50_calmultipcr_NSIBhug.idlsave'
; Gets: nyr,nhtit,yrmxd,prednh,fullnh,predse
if nhtit ne 'NH' then message,'Re-calibrate the age-banded!'
;
lcol=[10,11,12,13,14,16,15,17,50]
lthi=[4,4,12,4,4,4,4,4]
;
loadct,39
multi_plot,nrow=3,layout='large'
if !d.name eq 'X' then begin
  wintim,ysize=900
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,bold=0,font_size=16
endelse
def_1color,10,color='red'
def_1color,11,color='lpurple'
def_1color,12,color='mlblue'
def_1color,13,color='green'
def_1color,14,color='mdyellow'
def_1color,15,color='black'
def_1color,16,color='orange'
def_1color,17,color='magenta'
def_1color,18,color='white'
def_1color,22,color='vlblue'
def_1color,50,color='mgrey'
def_smearcolor,fromto=[18,22]
;
ndo=8
nyr=2002
alltime=findgen(nyr)
alldo=fltarr(ndo,nyr)*!values.f_nan
lowdo=fltarr(ndo,nyr)*!values.f_nan
namedo=['Jones ','Mann  ','Briffa','NChron','Overpeck','Crowley & Lowery',$
  'Instr','Esper']
iord=[3,0,4,5,1,7,2,6]
;
for jdo = 0 , ndo-1 do begin
  ido=iord(jdo)
  fac=1.
  ;
  case ido of
    0: begin        ; Phil's reconstruction
      alltit="Jones' Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000
      peren=1992
      openr,1,'../tree5/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)
      ; Convert from normalised values to deg C relative to 1961-1990
;      ts=ts*0.521 -0.1134       ; via Phil's variance matching (vs. NH all)
      case doseas of
        0: ts=ts*0.3690-0.2051       ; via regression vs. NH land>20 ann
        1: ts=ts*0.3369-0.1209       ; via regression vs. NH land>20 warm
        2: ts=ts*0.3949-0.2989       ; via regression vs. NH land>20 cold
      endcase
    end
    1: begin        ; Mike Mann's reconstruction (now use his 1000 yr one)
                    ; ALSO, NOW HAVE A LAND-NORTH-OF-20N ANNUAL RECON
      alltit="Mann's Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=1000
      peren=1980
      doln20=1
      if doln20 eq 0 then begin
        openr,1,'../tree5/mann_nhrecon1000.dat'
        nyr=981
        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]
;        ts=ts-0.12            ; convert to oC wrt 1961-90 (vs. NH all)
;        ts=ts*1.0641-0.0764   ; convert to oC wrt 1961-90 (vs. NH land>20)
        case doseas of
          0: ts=ts*1.1106-0.1645       ; via regression vs. NH land>20 ann
          1: ts=ts*0.9532-0.0886       ; via regression vs. NH land>20 warm
          2: ts=ts*1.2834-0.2480       ; via regression vs. NH land>20 cold
        endcase
      endif else begin
        openr,1,'../tree5/mannarea_all.dat'
        nyr=981
        headdat=' '
        rawdat=fltarr(11,nyr)
        readf,1,headdat
        readf,1,rawdat           ;,format='(I6,F11.7)'
        close,1
        timey=reform(rawdat(0,*))
        ts=reform(rawdat(10,*))
        kl=where((timey ge perst) and (timey le peren),nyr)
        timey=timey(kl)
        ts=ts[kl]
;OLD VERSION ts=ts*0.7246-0.0945   ; convert to oC wrt 1961-90 (vs. NH land>20)
        case doseas of
          0: ts=ts*0.8270-0.1771       ; via regression vs. NH land>20 ann
          1: ts=ts*0.6354-0.1061       ; via regression vs. NH land>20 warm
          2: ts=ts*0.9820-0.2602       ; via regression vs. NH land>20 cold
        endcase
      endelse
    end
    2: begin        ; Age-banded MXD
      alltit="Age-banded density NH growing-season reconstruction"
      ; Period to consider
      perst=1402
      peren=1960
;      fac=0.0       ; do not smooth it any further!!!
;      restore,filename='../treeharry/densadj_all(330).idlsave'
;      timey=x
;          ; CONVERSION FACTORS FOR AGE-BANDED MXD, BY REGRESSION ON INSTR.
;      ts=densadj*0.156525   ; converts it from density to temperature anom
      timey=newagetime
      ts=newagets
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
;      ts=ts(kl)-0.140369   ; to convert it oC wrt 1961-90
;
; Alternative overwrite with classically calibrated series
;
openr,10,'bandtemp_classiccalibrated_nh.dat'
headst=''
readf,10,headst
readf,10,headst
readf,10,nyr
readf,10,headst
readf,10,headst
readf,10,headst
rawdat=fltarr(2,nyr)
readf,10,rawdat
close,10
timey=reform(rawdat[0,*])
ts=reform(rawdat[1,*])
;
      case doseas of
        0: ts=ts*0.8996-0.1079       ; via regression vs. NH land>20 ann
        1: ts=ts*0.9281-0.0151       ; via regression vs. NH land>20 warm
        ; this isn't 1.0-0.0 because new instr data modified it
        2: ts=ts*0.7619-0.2271       ; via regression vs. NH land>20 cold
      endcase
    end
    3: begin         ; Torn+Yama+Taim
      perst=1
      peren=1993
      openr,2,'tornyamataim.ave'
      readf,2,nnn
      rawdat=fltarr(2,nnn)
      readf,2,rawdat
      close,2
      timey=reform(rawdat(0,*))
      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts[kl]
;      ts=ts*0.1342-0.2761    ; to convert it oC wrt 1961-90
      case doseas of
        0: ts=ts*0.1308-0.3650       ; via regression vs. NH land>20 ann
        1: ts=ts*0.1188-0.2664       ; via regression vs. NH land>20 warm
        2: ts=ts*0.1452-0.4746       ; via regression vs. NH land>20 cold
      endcase
    end
    4: begin         ; Overpeck
      perst=1600
      peren=1990
      fac=0.18       ; multiply filter by this because record is already 5-yrs
      openr,2,'overpeck_arctic.dat'
      headst=strarr(2)
      readf,2,headst
      readf,2,nnn
      rawdat=fltarr(2,nnn)
      readf,2,rawdat
      close,2
      timey=reverse(reform(rawdat(0,*)))
      ts=reverse(reform(rawdat(1,*)))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts[kl]
;      ts=ts(kl)*0.3352-0.0802    ; to convert it oC wrt 1961-90
      case doseas of
        0: ts=ts*0.3900-0.1593       ; via regression vs. NH land>20 ann
        1: ts=ts*0.2737-0.0947       ; via regression vs. NH land>20 warm
        2: ts=ts*0.4812-0.2303       ; via regression vs. NH land>20 cold
      endcase
    end
    5: begin         ; Crowley & Lowery
      fac=0.8        ; multiply filter by this 'cos record is already smooth
      restore,'/cru/u2/f055/data/paleo/crowley/crowley_lowery.'+$
        seasname[doseas]+'.idlsave'
      ; Gets crownyr,crowtimey,crowts
      timey=crowtimey
      ts=crowts
    end
    6: begin         ; Instrumental NH land > 20N, Apr-Sep (or other season!)
      perst=1871
      peren=1993
      print,'Reading temperatures'
      ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2_18512001.mon.nc')
      tsmon=crunc_rddata(ncid,tst=[1851,0],tend=[2001,11],grid=gtemp)
      ncdf_close,ncid
      nmon=12
      ntemp=gtemp.nt
      nyrtemp=ntemp/nmon
      timey=reform(gtemp.year,nmon,nyrtemp)
      timey=reform(timey(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
        0: ts=mkseason(nhmon,0,11,datathresh=6)
        1: ts=mkseason(nhmon,3,8,datathresh=3)
        2: ts=mkseason(nhmon,9,2,datathresh=3)
      endcase
      warmyr=timey
      warmts=mkseason(nhmon,3,8,datathresh=3)
      case doseas of
        0: warmts=warmts*0.8929-0.1058   ; via regression vs. NH land>20 ann
        1: warmts=warmts*1.0000-0.0000   ; via regression vs. NH land>20 warm
        2: warmts=warmts*0.7365-0.2285   ; via regression vs. NH land>20 cold
      endcase
;      openr,2,'../treeharry/nhland20_amjjas.dat'
;      readf,2,nnn
;      rawdat=fltarr(2,nnn)
;      readf,2,rawdat
;      close,2
;      print,rawdat(0,0)
;      timey=reform(rawdat(0,*))
;      ts=reform(rawdat(1,*))
      kl=where((timey ge perst) and (timey le peren),nyr)
      timey=timey(kl)
      ts=ts(kl)
      warmyr=warmyr[kl]
      warmts=warmts[kl]
    end
    7: begin        ; Esper's reconstruction
      alltit="Esper et al.'s Northern Hemisphere temperature reconstruction"
      ; Period to consider
      perst=831
      peren=1992
      ; Read raw reconstruction
      openr,1,'/cru/u2/f055/data/paleo/esper2002/esper.txt'
      readf,1,nyr
      headst=''
      readf,1,headst
      rawdat=fltarr(7,nyr)
      readf,1,rawdat
      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)
      ; Convert from raw values to deg C relative to 1961-1990
      case doseas of
        0: ts=ts*1.6664-2.1539       ; via regression vs. NH land>20 ann
        1: ts=ts*1.4724-1.8442       ; via regression vs. NH land>20 warm
        2: ts=ts*2.0851-2.7289       ; via regression vs. NH land>20 cold
      endcase
    end
  endcase
  ;
  ; Now plot the results
  ;
  if fac gt 0. then filter_cru,thalf*fac,/nan,tsin=ts,tslow=tslow $
               else tslow=ts
  ;
  if jdo eq 0 then begin
    yrange=[-0.75-0.2*dobore-(doseas eq 2)*0.2,0.4+dobigy*0.05+$
        (doseas eq 2)*0.1]
    if thalf eq 100. then yrange=[-0.95,0.2]
;    yrange=[-0.9,0.6]
    yrange=[-1.35,0.5]
    plot,timey,tslow,/nodata,$
      xstyle=1,xrange=[dost,2001],xtitle='Year  (AD)',$
      /ystyle,ytitle='Temperature anomaly  (!Uo!NC wrt 1961-90)',$
      yrange=yrange
    ttt=newagetime
    tts=newagelow
    tse=newagese
    tkk=where((ttt ge 1402) and (ttt le 1960))
    ttt=ttt(tkk) & tts=tts(tkk) & tse=tse(tkk)
    if doerr eq 1 then begin
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse*2.,reverse(tts+tse*2.)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
      polyfill,xfill,yfill,color=19,noclip=0
      xfill=[ttt,reverse(ttt)]
      yfill=[tts-tse,reverse(tts+tse)]
      kl=where(finite(yfill),nkeep)
      yfill=yfill(kl)
      xfill=xfill(kl)
      polyfill,xfill,yfill,color=21,noclip=0
    endif
    oplot,!x.crange,[0.,0.],linestyle=1
;
; Optionally overlay pre-computed borehole temperature estimates
;
if dobore eq 1 then begin
  restore,'/cru/u2/f055/data/paleo/borehole/borehole_NH.idlsave'
  ; Gets: boretime, boresite, boregrid
  ; Compute and apply an offset so that the borehole line crosses zero
  ; in 1975.5 (i.e., mid-way through the 1961-1990 reference period)
  v=boresite[0:1]
  voff=(v[0]-v[1])*(2000.-1975.5)/(2000.-1900.)-v[0]
;  oplot,boretime,boresite+voff,color=50,thick=6,linestyle=2
  v=boregrid[0:1]
  voff=(v[0]-v[1])*(2000.-1975.5)/(2000.-1900.)-v[0]
  oplot,boretime,boregrid+voff,color=50,thick=6
endif
  ;
  endif
  ;
  oplot,timey,tslow,thick=lthi(ido),color=lcol(ido)
print,min(tslow,/nan)
  ;
  ; For Overpeck, let's interpolate to get annual values from decadal ones
  ;
  if ido eq 4 then begin
    inyr=timey(nyr-1)-timey(0)+1
    itimey=findgen(inyr)+timey(0)
    lowits=interpol(tslow,timey,itimey)
    its=interpol(ts,timey,itimey)
    ts=its
    tslow=lowits
    timey=itimey
  endif
  ;
  ; Store timeseries for later cross-correlation
  ;
  ist=where(alltime eq timey(0))
  alldo(ido,ist(0):ist(0)+n_elements(ts)-1)=ts(*)
  lowdo(ido,ist(0):ist(0)+n_elements(ts)-1)=tslow(*)
  ;
endfor
;
if doinstrwarm eq 1 then begin
  filter_cru,thalf,/nan,tsin=warmts,tslow=warmlow
  oplot,warmyr,warmlow,thick=4
endif
;
case doerr of
  0: terr=''
  1: terr=' (& errors)'
endcase
case dost of
  1: iiipos=50
  850: iiipos=880
  1000: iiipos=1025
  1400: iiipos=1415
  1800: iiipos=1940
endcase
if dost le 1000 then jjjpos=-0.39 else jjjpos=0.35
iiipos=!x.crange[0]
jjjpos=!y.crange[1]
legpos=convert_coord([iiipos],[jjjpos],/data,/to_normal)
legord=[4,0,1,5,3,2,7,6]
;legtxt=['Jones et al. (1998)','Mann et al. (1999)','This study'+terr,$
;    'Briffa (2000)','Overpeck et al. (1997)','Crowley & Lowery (2000)',$
;    'Observations']
;legtxt=['Jones et al. (1998)','Mann et al. (1999)','Briffa et al. (2000)',$
;    'Briffa (2000)','Overpeck et al. (1997)','Crowley & Lowery (2000)',$
;    'Observations','Esper et al. (2002)']
;legend,position=legpos,legtxt(legord),$
;  thick=lthi(legord),color=lcol(legord)
legtxt=['Jones98','Mann99','Briffa01',$
    'Briffa00','Overpeck97','Crowley00',$
    'Obs','Esper02']
legend,position=legpos,legtxt(legord),horiz=1,$
  margin=0.5,pspacing=1,thick=lthi(legord),color=lcol(legord)
xyouts,dost+75,!y.crange[1]-0.27,mtit[doseas]
;
if docorr ne 0 then begin
ndo=6
for kkk = 1 , 1 do begin
  if kkk eq 0 then yrst=1000 else yrst=1402
  print
  print,'Analysis period:',yrst,1960
  print
  print,'TS1     TS2     Raw   <50yr   >50yr'
for ido = 0 , ndo-2 do begin
  for jdo = ido+1 , ndo-1 do begin
    x=reform(alldo(ido,*))
    y=reform(alldo(jdo,*))
    fil2=[thalf,thalf]
    if namedo(ido) eq 'Overpeck' then fil2(0)=thalf*0.8
    if namedo(jdo) eq 'Overpeck' then fil2(1)=thalf*0.8
    if namedo(ido) eq 'Crowley & Lowery' then fil2(0)=thalf*0.8
    if namedo(jdo) eq 'Crowley & Lowery' then fil2(1)=thalf*0.8
    print,fil2
;    allr1=mkcorrelation(x,y,alltime,filter=10.,refperiod=[yrst,1960])
    allr=mkcorrelation(x,y,alltime,filter=fil2,refperiod=[yrst,1960])
    print,namedo(ido),namedo(jdo),allr,format='(2A8,6F8.2)'
  endfor
endfor
endfor
endif
;
openw,1,'multiple_reconstructions.dat'
printf,1,'Data file to accompany Plate 3 of Briffa et al. (2001).'
printf,1
printf,1,'The following reconstructions have been taken from the source'
printf,1,'references listed below, and then RECALIBRATED to obtain estimates'
printf,1,'of April-September mean temperatures from all land regions north'
printf,1,'of 20N.  All series are temperature anomalies in degrees C with'
printf,1,'respect to the 1961-1990 mean.  Note that in Plate 3 of Briffa'
printf,1,'et al. (2001) all the series had been smoothed, while the data'
printf,1,'listed below are unsmoothed (though column 5 never had any'
printf,1,'sub-5-year variability even when unsmoothed).'
printf,1
printf,1,'Full references and details of the recalibration are given in'
printf,1,'Briffa et al. (2001) J Geophys Res 106, 2929-2941.'
printf,1
printf,1,'Columns are:'
printf,1,'  1: Jones et al. (1998) Holocene'
printf,1,'  2: Mann et al. (1999) Geophys Res Lett'
printf,1,'  3: Briffa et al. (2001) J Geophys Res'
printf,1,'  4: Briffa (2000) Quat Sci Rev'
printf,1,'  5: Overpeck et al. (1997) Science'
printf,1,'  6: Crowley & Lowery (2000) Ambio'
printf,1,'  7: Observed temperatures from Jones et al. (1999) Rev Geophys'
printf,1
printf,1,'  Year         1         2         3         4         5         6         7'
ist=where(alltime eq dost)
ist=ist(0)
nyr=n_elements(alltime)
ml=where(finite(alldo) eq 0)
alldo(ml)=-999.999
for i = ist , nyr-1 do begin
  printf,1,alltime(i),alldo(*,i),format='(I6,7F10.3)'
endfor
close,1
;
if (dost eq 850) and (thalf eq 50) then begin
  if doseas eq 0 then fnadd='' else fnadd=seasname[doseas]
  for ijk = 0 , 1 do begin
    if ijk eq 0 then openw,1,'blowinghotandcold.dat'+fnadd $
                else openw,1,'blowinghotandcold.dat'+fnadd+'.unsmoothed'
    printf,1,'Data file to accompany:'
    printf,1,'Briffa & Osborn (2002) Blowing hot and cold. Science 295, 2227-2228.'
    printf,1
    printf,1,'The figure in Briffa & Osborn (2002) was used to place a new'
    printf,1,'reconstruction by Esper et al. (2002) in the context of other'
    printf,1,'reconstructions of past temperature variability over large parts'
    printf,1,'of the Northern Hemisphere.'
    printf,1
    printf,1,'IMPORTANT NOTES:'
    printf,1
    printf,1,'(1) The reconstructions have been taken from the source'
    printf,1,'references listed below, and then RECALIBRATED to obtain estimates'
    printf,1,'of '+seasname[doseas]+' mean temperatures from all land regions north'
    printf,1,'of 20N.  This was because, for our purposes, we wanted to ensure'
    printf,1,'that common data and a common period had been used for calibration (accepting'
    printf,1,'that the selected season and spatial coverage may not be optimal'
    printf,1,'for all records).'
    printf,1
    printf,1,'(2) For series 2, from Mann et al., we computed the mean of their'
    printf,1,'spatially-resolved reconstruction from only land grid boxes north'
    printf,1,'of 20N, and used that series (rather than their published full-NH'
    printf,1,'series) for our calibration.'
    printf,1
    printf,1,'(3) All series are temperature anomalies in degrees C with'
    printf,1,'respect to the observed 1961-1990 mean (which may differ from each'
    printf,1,"series' own 1961-1990 mean)."
    printf,1
    if ijk eq 0 then begin
      printf,1,'(4) All series have been smoothed with a 50-year filter.  Beware of end'
      printf,1,'effects introduced because we take the smoothed values right up to the'
      printf,1,'ends of the record (by padding with the local mean).  Not also that we'
      printf,1,'cut off the instrumental record in 1993 (the last year in which any of'
      printf,1,'the proxy reconstructions have data) in order to prevent post-1993'
      printf,1,'values influencing the smoothed instrumental record in a way that would'
      printf,1,'bias the comparison with the proxy reconstructions.'
    endif else begin
      printf,1,'(4) In the paper, all series had been smoothed with a 50-year filter,'
      printf,1,'BUT IN THIS FILE WE PROVIDE UNSMOOTHED TIMESERIES (though column 5'
      printf,1,'never had any sub-5-year variability even when unsmoothed).'
      printf,1,'Not also that we cut off the instrumental record in 1993 (the last'
      printf,1,'year in which any of the proxy reconstructions have data).'
    endelse
    printf,1
    printf,1,'(5) Our figure also included three other curves which are not'
    printf,1,'included here (two showing reconstructions from borehole temperature'
    printf,1,'profiles, obtained using different averaging processes, and the'
    printf,1,'third was what you get if you try to reconstruct annual temperatures'
    printf,1,'from instrumental warm-season temperatures!).'
    printf,1
    printf,1,'Columns are:'
    printf,1,'  1: Jones et al. (1998) Holocene'
    printf,1,'  2: Mann et al. (1999) Geophys Res Lett'
    printf,1,'  3: Briffa et al. (2001) J Geophys Res'
    printf,1,'  4: Briffa (2000) Quat Sci Rev'
    printf,1,'  5: Overpeck et al. (1997) Science'
    printf,1,'  6: Crowley & Lowery (2000) Ambio'
    printf,1,'  7: Observed temperatures from Jones et al. (1999) Rev Geophys'
    printf,1,'  8: Esper et al. (2002) Science'
    printf,1
    printf,1,'  Year         1         2         3         4         5         6         7         8'
    ;
    if ijk eq 0 then alldat=lowdo else alldat=alldo
    ;
    ist=where(alltime eq dost)
    ist=ist(0)
    nyr=n_elements(alltime)
    ml=where(finite(alldat) eq 0,nmiss)
    if nmiss gt 0 then alldat(ml)=-999.999
    for i = ist , nyr-1 do begin
      printf,1,alltime(i),alldat(*,i),format='(I6,8F10.3)'
    endfor
    close,1
  endfor
endif
;
;plotoff
;spawn,'ps2web'
;
end