; ; Combines the multi-PCR models into one reconstruction, including ; combining and interpolating the time-dependent errors. ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=4,layout='large' if !d.name eq 'X' then begin window,ysize=850 !p.font=-1 endif else begin !p.font=0 device,/helvetica,bold=0,font_size=14 endelse def_1color,20,color='lsand' def_1color,21,color='lyellow' def_1color,22,color='red' def_1color,23,color='mblue' def_1color,24,color='lgreen' ; ; ***DECIDE WHETHER TO REPLACE THE NSIB SERIES BY IT'S HUGERSHOFF VERSION, ; ***BECAUSE THE FORMER HAS HIGH VARIANCE ; dohug=0 ; 0=age-banded NSIB 1=hugershoff NSIB if dohug eq 1 then hadd='_NSIBhug' else hadd='' ; for jjj = 0 , 3 do begin case jjj of 0: begin dosmooth=0 thalf=10 end 1: begin dosmooth=1 thalf=10 end 2: begin dosmooth=1 thalf=30 end 3: begin dosmooth=1 thalf=50 end endcase ; regname='NH' ; 'ALL' or 'NH' ;dosmooth=1 ; 0=errors for interannual data 1=for smoothed data ;thalf=25. ; fnhead='bandtemp'+regname if dosmooth eq 1 then smst='sm'+string(thalf,format='(I2.2)') else smst='' fntail=smst+'_calpcr'+hadd+'.idlsave' ; ; Read in the different PCR models for different fixed grids ; perlist=[ 0, 0, 0, 1, 2, 3, 4, 5, 5, 6, 6, 7, 8, 9] fixlist=[1950,1800,1700,1700,1700,1700,1600,1600,1500,1500,1400,1950,1988,1988] nver=n_elements(perlist) ; for jver = 0 , nver-1 do begin iver=perlist(jver) verst=string(iver,format='(I1)') ifix=fixlist(jver) fixst=string(ifix,format='(I4)') ; fn=fnhead+verst+'_fixed'+fixst+fntail print,fn restore,filename=fn ; Gets: nyr,nhtit,yrmxd,prednh,fullnh,predserr,nhtemp,tempyr ; if jver eq 0 then begin allnyr=nyr alltimey=yrmxd allts=fltarr(nyr,nver) allse=fltarr(nyr,nver) endif allts(*,jver)=prednh(*) allse(*,jver)=predserr(*) ; endfor ; ; Make a combined reconstruction using different models for different periods, ; and different fixed grid error estimates ; prednh=fltarr(nyr) & prednh(*)=!values.f_nan predse=fltarr(nyr) & predse(*)=!values.f_nan ; for iyr = 0 , nyr-1 do begin jyr=alltimey(iyr) ; ; Identify which model to use ; iuse=13 if jyr le 1991 then iuse=12 if jyr le 1988 then iuse=11 if jyr le 1981 then iuse=0 if jyr le 1682 then iuse=3 if jyr le 1666 then iuse=4 if jyr le 1624 then iuse=5 if jyr le 1603 then iuse=6 if jyr le 1587 then iuse=7 if jyr le 1479 then iuse=9 ; ; Use the selected model prediction for this year ; if jyr ge 1402 then begin prednh(iyr)=allts(iyr,iuse) ; ; Also use the selected model to provide the standard error, UNLESS we ; are in one of the following periods ; se=allse(iyr,iuse) if (jyr ge 1683) and (jyr le 1700) then begin se=allse(iyr,2) endif if (jyr ge 1701) and (jyr le 1800) then begin se=interpol([allse(iyr,2),allse(iyr,1)],[1700,1800],[jyr]) endif if (jyr ge 1801) and (jyr le 1899) then begin se=interpol([allse(iyr,1),allse(iyr,0)],[1800,1900],[jyr]) endif if (jyr ge 1480) and (jyr le 1499) then begin se=allse(iyr,8) endif if (jyr ge 1500) and (jyr le 1587) then begin se=interpol([allse(iyr,8),allse(iyr,7)],[1500,1600],[jyr]) endif if (jyr ge 1402) and (jyr le 1479) then begin se=interpol([allse(iyr,10),allse(iyr,9)],[1400,1500],[jyr]) endif predse(iyr)=se ; endif ; endfor ; plot,[0,1],/nodata,$ /xstyle,xrange=[1400,1994],xtitle='Year',$ /ystyle,yrange=[-2.2,1.0],ytitle='Temperature (!Uo!NC wrt 1961-90)',$ title=regname+' '+smst+' '+hadd ; filter_cru,thalf,/nan,tsin=nhtemp,tslow=temlow filter_cru,thalf,/nan,tsin=prednh,tslow=agelow if dosmooth eq 1 then begin agets=agelow temts=temlow endif else begin agets=prednh temts=nhtemp endelse ; timey=alltimey serr=predse xfill=[timey,reverse(timey)] yfill=[agets-serr*2.,reverse(agets+serr*2.)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) polyfill,xfill,yfill,color=20,noclip=0 ; xfill=[timey,reverse(timey)] yfill=[agets-serr,reverse(agets+serr)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) polyfill,xfill,yfill,color=21,noclip=0 ; oplot,!x.crange,[0.,0.],linestyle=1 ; oplot,tempyr,temts,thick=1+3*dosmooth,color=22 oplot,timey,agets,thick=1+3.*dosmooth ; ; Save the composite, multi-model reconstruction with errors for later use ; yrmxd=alltimey prednh=agets nhtit=regname nyr=allnyr ;save,filename='bandtemp'+regname+smst+'_calmultipcr'+hadd+'.idlsave',$ ; yrmxd,prednh,nhtit,nyr,predse ; ; Also output the series for MATLAB spectral analysis ; Stop in 1994 cos the Hugershoff doesn't have 1995 in it ; Start in 1402 cos 1400-1 are missing in age-banded multi-PCR ; ;if regname eq 'ALL' then begin ; openw,1,'ts2.matlab_mxdband_all' ; printf,1,prednh(2:nyr-2),format='(F7.2)' ; close,1 ;endif ; ; ; Output the final time series to an ASCII file too ; ;if jjj eq 0 then begin openw,1,regname+'temp'+smst+'_agebandbriffa'+hadd+'.dat' if regname eq 'ALL' then begin printf,1,'Calibrated against observed Apr-Sep temperature over 1881-1960' printf,1,'averaged over all grid boxes containing tree-ring chronologies' printf,1,'(i.e., co-located temperatures).' endif else begin printf,1,'Calibrated against observed Apr-Sep temperature over 1881-1960' printf,1,'averaged over all land grid boxes north of 20N' printf,1 endelse printf,1 if dosmooth then printf,1,'Time series and errors smoothed by '+txt(thalf)+' years' $ else printf,1 printf,1 if dohug eq 1 then begin printf,1,'NSIB series was replaced by the Hugershoff standardised version!' endif printf,1 printf,1,'Year Reconstructed temperature anomaly (degrees C wrt 1961-90) kkk=where(finite(prednh),nkkk) for i = 0 , nkkk-1 do printf,1,yrmxd(kkk[i]),prednh(kkk[i]),predse[kkk[i]],format='(I4,2F10.3)' close,1 ;endif ; endfor ; end