; ; Synthesises together multiple quasi-NH recons ; if n_elements(doseas) eq 0 then doseas=1 ; 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=0 ; 0=no boreholes, 1=overlay borehole-based NH estimate doerr=0 ; 0=no errors, 1=Briffa age errors dost=1850 ; start year docorr=0 ; 1=compute correlations 0=don't thalf=0. ; filter for high or low-pass dobigy=1 ; 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] if thalf gt 0 then begin lthi=[6,6,6,6,6,6,4,6] endif else begin lthi=[2,2,2,2,2,2,4,2] endelse ; loadct,39 multi_plot,nrow=2,layout='large' if !d.name eq 'X' then begin lthi=lthi*0.5 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='brightblue' def_1color,16,color='orange' def_1color,17,color='magenta' def_1color,50,color='mgrey' ; def_1color,18,color='white' def_1color,22,color='vlblue' def_smearcolor,fromto=[18,22] def_1color,21,color='mdgrey' def_1color,22,color='lgrey' ; ndo=8 allnyr=2002 alltime=findgen(allnyr) alldo=fltarr(ndo,allnyr)*!values.f_nan lowdo=fltarr(ndo,allnyr)*!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 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=1881 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 ; if fac gt 0. then filter_cru,thalf*fac,/nan,tsin=ts,tslow=tslow $ else tslow=ts ; ; 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 ; ; Having read in all series, we want to synthesise them together. ; uselist=[0,1,2,3,4,5,7] ; don't use 6, the instrumental data nuse=n_elements(uselist) ; niter=1 ; 12? for ij = 0 , niter-1 do begin ; ; Now let's re-calibrate each over the period 1000-1850 (or whatever subset ; each has data for) against the mean series. ; if ij gt 0 then begin allnew=fltarr(ndo,allnyr)*!values.f_nan lownew=fltarr(ndo,allnyr)*!values.f_nan for jdo = 0 , nuse-1 do begin ido=uselist[jdo] ts1=reform(lowdo[ido,*]) ts2=reform(alldo[ido,*]) kl=where((alltime ge 1600) and (alltime le 1850) and (finite(ts1) ne 0),$ nkeep) print,'re-calibrating' acoeff=linfit(ts1[kl],alllow[kl]) print,namedo[ido],nkeep,acoeff lownew[ido,*]=acoeff[0]+acoeff[1]*ts1[*] allnew[ido,*]=acoeff[0]+acoeff[1]*ts2[*] endfor lowdo[uselist,*]=lownew[uselist,*] alldo[uselist,*]=allnew[uselist,*] endif ; ; Compute their average and calibrate this against the instrumental data, ; estimating uncertainty/errors too! ; allts=alldo[uselist,*] alltarg=reform(alldo[6,*]) allave=multits_calib(allts,alltime,alltarg,alltime,calibper=[1881,1960],$ stderr=serr,no_anomaly=0,no_varadj=0,smoothlist=[thalf],plotdiag=0) ; ; Optionally keep the uncalibrated average of the multiple reconstructions ; if ij eq 0 then begin print,'keeping average' keepave=total(allts,1,/nan)/float(total(finite(allts),1)) filter_cru,thalf,/nan,tsin=keepave,tslow=keeplow endif ; ; Make smoothed series for plotting ; filter_cru,thalf,/nan,tsin=allave,tslow=alllow !p.multi[0]=0 ; ;print,'Calibrate average vs. observed:',r,acoeff ; ; Now plot the results ; if (ij eq 0) or (ij eq niter-1) then begin for jdo = 0 , ndo-1 do begin ido=iord(jdo) tslow=reform(lowdo[ido,*]) ; 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.8,0.6] ; yrange=[-0.75,0.25] pause plot,alltime,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) ttt=alltime tts=alllow tse=serr if doerr eq 1 then begin print,'PLOTTING ERRORS!' xfill=[ttt,reverse(ttt)] yfill=[tts-tse*2.,reverse(tts+tse*2.)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) kerr=where(xfill ge 1600) polyfill,xfill[kerr],yfill[kerr],color=21,noclip=0 xfill=[ttt,reverse(ttt)] yfill=[tts-tse,reverse(tts+tse)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) kerr=where(xfill ge 1600) polyfill,xfill[kerr],yfill[kerr],color=22,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,noclip=1 endif ; endif ; oplot,alltime,tslow,thick=lthi(ido),color=lcol(ido) 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 950: iiipos=950 1000: iiipos=1025 1400: iiipos=1415 1800: iiipos=1940 else: iiipos=dost 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.17,mtit[doseas] ; if (ij ne niter-1) or (ij eq 0) then begin ; oplot,alltime,alllow,thick=4+4*(thalf gt 0) ; oplot,alltime,keeplow,thick=2+(thalf gt 0) endif endif ; endfor ; end