; ; Reads in the age-banded regional series produced by Harry's program ; saband10.f. For each region, averages across bands 50yr-100yr, then ; across bands 50yr-300yr. Both these averaging processes are variance- ; adjusted, and from the second one the rbar is used to compute a time- ; dependent EPS. The regional mean series are then averaged to form a ; northern hemisphere series, weighted by their EPS, and variance adjusted. ; ; *** Now averages bands 20-100yr first (not just 50-100yr), although they ; *** still only count as a single band. And also uses bands 300-550yr, ; *** after first averaging them into a single band (weighting & adjusted etc.) ; ; *** Can optionally exclude any band with one 1 tree in it. ; ; *** Does it for fixed grids *** ; fixedyr=1950 fnfix=string(fixedyr,format='(I4)') if fixedyr eq 1950 then fnfix='' ; ; Define region names. Names are different to those I use, but order ; should be the same to make later comparison easier. ; regname=['northwestern.europe',$ 'southwestern.europe',$ 'northern.siberia',$ 'eastern.siberia',$ 'centralsiberia',$ 'tibetianplateau',$ 'western.northamerica',$ 'northwestern.canada',$ 'centraleastern.canada'] nreg=n_elements(regname) ; ; Specify the number of bands used in each region ; ; **These values are now overwritten, because they're lower for the early ; **fixed grids! nband=[19,20,21,24,18,20,22,17,17] ncol=nband*3 + 1 ; allnyr=2000 allyr=findgen(allnyr)+1. kper=where((allyr ge 1400) and (allyr le 1995),nyr) timey=allyr(kper) ; allts=fltarr(nreg,nyr)*!values.f_nan alleps=fltarr(nreg,nyr)*!values.f_nan rbarper=[1700,1994] ; loadct,39 multi_plot,nrow=5 if !d.name eq 'X' then begin window,ysize=850 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=12 endelse ; ; Repeat for each region ; for ireg = 0 , nreg-1 do begin ; ; Read raw data ; fn='mxd.'+regname(ireg)+'.all.BS5T1S1.17002000.sa.dat'+fnfix print,fn ; ; Check the number of columns ; spawn,'head -1 '+fn+' | wc -w',osoutput ncol(ireg)=fix(osoutput) nband(ireg)=(ncol(ireg)-1)/3 ; if nband(ireg) lt 5 then begin print,'Skipping region!' endif else begin ; openr,1,fn headlin=fltarr(ncol(ireg)) readf,1,headlin rawdat=fltarr(ncol(ireg),allnyr) readf,1,rawdat close,1 ; ; Just keep the 1400 to 1994 portions ; rawdat=rawdat(*,kper) ; ; Locate missing data ; ml=where(rawdat eq 9990.,nmiss) rawdat(ml)=!values.f_nan ; colts=indgen(nband(ireg))*3+1 ; ; First stage of averaging: bands 50 to 100 yrs ; print,'BANDS 50 TO 100 YEARS' bandlist=indgen(8)+2 ; was (5)+5 collist=colts(bandlist) tsdat=rawdat(collist,*) ; Optionally remove those series made up from <2 trees ntree=rawdat(collist+2,*) misslist=where(ntree lt 2,nmiss) if nmiss gt 0 then tsdat(misslist)=!values.f_nan ; End of option totdat=total(tsdat,1,/nan) numdat=total(finite(tsdat),1) tsold=totdat/float(numdat) tsnew=var_adjust(tsdat,timey,refperiod=rbarper,rbar=fullrbar,/mkrbar,$ /normalise) ; If 1700-1994 had insufficient data to compute rbar, use full period if finite(fullrbar) eq 0 then begin print,'Re-doing variance adjustment using full period to compute rbar' tsnew=var_adjust(tsdat,timey,rbar=fullrbar,/mkrbar,/normalise) endif ; pause plot,timey,numdat,/xstyle,/ystyle,yrange=[0,8],$ title='rbar='+string(fullrbar,format='(F6.2)') neff=float(numdat)/(1.+(float(numdat)-1.)*fullrbar) oplot,timey,neff,thick=4 ; ; Store in the 90 to 100 band to make following averaging simpler ; Also replace the ntree column of 90-100 to make sure that no data is ; subsequently removed! ; rawdat(colts(9),*)=tsnew rawdat(colts(9)+2,*)=100. ; ; Second stage of averaging: bands 300 to 550 yrs ; print,'BANDS 300 TO 550 YEARS' bandlist=indgen(5)+14 collist=colts(bandlist) tsdat=rawdat(collist,*) ; Optionally remove those series made up from <2 trees ntree=rawdat(collist+2,*) misslist=where(ntree lt 2,nmiss) if nmiss gt 0 then tsdat(misslist)=!values.f_nan ; End of option totdat=total(tsdat,1,/nan) numdat=total(finite(tsdat),1) tsold=totdat/float(numdat) tsnew=var_adjust(tsdat,timey,refperiod=rbarper,rbar=fullrbar,/mkrbar,$ /normalise) ; If 1700-1994 had insufficient data to compute rbar, use full period if finite(fullrbar) eq 0 then begin print,'Re-doing variance adjustment using full period to compute rbar' tsnew=var_adjust(tsdat,timey,rbar=fullrbar,/mkrbar,/normalise) endif ; ; Store in the 300 to 350 band to make following averaging simpler ; rawdat(colts(14),*)=tsnew rawdat(colts(14)+2,*)=100. ; ; Next stage of averaging: bands 50 to 300 yrs ; bandlist=indgen(6)+9 ; was (5)+9 collist=colts(bandlist) nuse=n_elements(collist) ; ; Compute and plot running variance of 5 input series ; print,'COMPUTING INPUT VARIANCE' allvar=fltarr(nuse,nyr)*!values.f_nan wlen=50 for iyr = 0 , nyr-wlen do begin ; print,timey(iyr),format='($,I4)' for iband = 0 , nuse-1 do begin oneseg=reform(rawdat(collist(iband),iyr:iyr+wlen-1)) kl=where(finite(oneseg),nkeep) if nkeep ge 25 then begin ; print,'.',format='($,A1)' onemon=moment(oneseg(kl)) allvar(iband,iyr)=onemon(1) endif endfor endfor ; print totvar=total(allvar,1,/nan) numvar=total(float(allvar),1) meanvar=totvar/float(numvar) plot,timey,totvar,/xstyle,ytitle='Mean input variance' ; print,'BANDS 50 TO 300 YEARS' tsdat=rawdat(collist,*) ; Optionally remove those series made up from <2 trees ntree=rawdat(collist+2,*) misslist=where(ntree lt 2,nmiss) if nmiss gt 0 then tsdat(misslist)=!values.f_nan ; End of option totdat=total(tsdat,1,/nan) numdat=total(finite(tsdat),1) tsold=totdat/float(numdat) tsnew=var_adjust(tsdat,timey,refperiod=rbarper,rbar=fullrbar,/mkrbar,$ /normalise) ; plot,timey,numdat,/xstyle,/ystyle,yrange=[0,7],$ title='rbar='+string(fullrbar,format='(F6.2)') neff=float(numdat)/(1.+(float(numdat)-1.)*fullrbar) oneeps=fullrbar / (fullrbar + (1.-fullrbar)/float(numdat) ) oplot,timey,neff,thick=4 ; plot,timey,oneeps,/xstyle,/ystyle,yrange=[0.,1.] title='EPS' ; ; Also compute a time-dependent correlation matrix, and then rbar(t) ; and then a new eps(t). If rbar becomes very low, or if it becomes ; zero simply because there is only one band with data, then we use ; the oldest positive rbar value. ; tdrbar=fltarr(nyr) print,'Century rbar:' for icent = 0 , 11 do begin stcent=float(icent)*50.+1400. print,stcent,format='($,I5)' kcent=where((timey ge stcent) and (timey le stcent+49.)) centts=tsdat(*,kcent) tdrbar(kcent)=mkrbar(centts) > 0. endfor for iyr = nyr-100 , 0 , -1 do begin if numdat(iyr) le 0. then begin tdrbar(iyr)=!values.f_nan endif else begin if finite(tdrbar(iyr)) eq 0 then tdrbar(iyr)=tdrbar(iyr+1) $ else if tdrbar(iyr) le 0. then tdrbar(iyr)=tdrbar(iyr+1) endelse endfor tdeps=tdrbar / (tdrbar + (1.-tdrbar)/float(numdat) ) ; ; Filter (optional) the EPS series so that it doesn't suddenly change ; just because a tree moves from one age band to another ; print,'Smoothing EPS' filter_cru,20.,/nan,tsin=tdeps,tslow=tdepslow tdeps=tdepslow ; OPTION TO GIVE UNIFORM WEIGHTING OF ALL REGIONS!!!!! ; tdeps(*)=1. ; oplot,timey,tdeps,thick=3 oplot,timey,tdrbar,thick=5 ; allts(ireg,*)=tsnew(*) alleps(ireg,*)=tdeps(*) ; y=reform(allts(ireg,*)) cpl_barts,timey,y,$ /xstyle,xtitle='Year',$ ytitle='Normalised MXD anomaly',$ title=regname(ireg)+fnfix filter_cru,/nan,20.,tsin=y,tslow=ylow oplot,timey,ylow,thick=4 oplot,!x.crange,[0.,0.] ; endelse ; endfor ; ; Average across regions to give NH ; nhtot=total(allts,1,/nan) nhnum=total(finite(allts),1) nhold=nhtot/float(nhnum) ; ; Compute rbar between the regions for each 50-yr period, to measure any ; deteriation. ; rbarvary=fltarr(nyr) print,'Changing rbar:' for icent = 0 , 11 do begin stcent=float(icent)*50.+1400. print,stcent,format='($,I5)' kcent=where((timey ge stcent) and (timey le stcent+49.)) centts=allts(*,kcent) rbarvary(kcent)=mkrbar(centts) endfor for iyr = nyr-100 , 0 , -1 do begin if nhnum(iyr) le 0. then begin rbarvary(iyr)=!values.f_nan endif else begin if finite(rbarvary(iyr)) eq 0 then rbarvary(iyr)=rbarvary(iyr+1) endelse endfor ; nhnew=var_adjust(allts,timey,weight=alleps,$ refperiod=rbarper,rbar=fullrbar,/mkrbar,/normalise,/td_rbar) ; nhts=nhnew ; ; Plot it ; pause plot,timey,rbarvary,thick=3,title='Varying rbar between regions',$ yrange=[-0.2,0.5],/xstyle,/ystyle oplot,!x.crange,[0.,0.],linestyle=1 ; ;plot,timey,timey,/nodata,/xstyle,/ystyle,yrange=[0.,1.],$ ; title='Time-dependent EPS' ;for ireg = 0 , nreg-1 do oplot,timey,alleps(ireg,*),thick=ireg+1 ; plot,timey,fullrbar,/xstyle,/ystyle,yrange=[0.,0.5],$ title='Time-dependent rbar (due to regions dropping out)' ; plot,timey,nhnum,/xstyle,/ystyle,yrange=[0,10] neff=float(nhnum)/(1.+(float(nhnum)-1.)*fullrbar) oplot,timey,neff,thick=4 ; y=nhts cpl_barts,timey,y,$ /xstyle,xtitle='Year',$ ytitle='Normalised MXD anomaly',$ title='NH mean using age-banded MXD, fixed grid: '+fnfix filter_cru,/nan,20.,tsin=y,tslow=ylow oplot,timey,ylow,thick=4 oplot,!x.crange,[0.,0.] ; print,'COMPUTING OUTPUT VARIANCE' nhvar=fltarr(nyr)*!values.f_nan wlen=50 for iyr = 0 , nyr-wlen do begin print,timey(iyr),format='($,I4)' oneseg=nhts(iyr:iyr+wlen-1) kl=where(finite(oneseg),nkeep) if nkeep ge 10 then begin onemon=moment(oneseg(kl)) nhvar(iyr)=onemon(1) endif endfor print ; plot,timey,nhvar,/xstyle,ytitle='NH variance',$ xtitle='Start year of 50-yr sliding window' ; ; Now save the data for subsequent analysis ; save,filename='bandregnh_fixed'+fnfix+'.idlsave',$ timey,nyr,nhts,allts,alleps,nreg,regname ; end