; ; Creates regional timeseries based on a limited sample of MXD series, ; variance adjusted with a time-dependent rbar, and compares them with ; a regional instrumental temperature record which is also computed here ; and variance adjusted. ; ; The instrumental record can be based on either all available land ; observations in the region, or only those with an MXD site within their grid ; box (in the latter case, the region is defined by the FULL sample of MXD ; data, not by the limited sample of MXD sites actually being used. ; The instrumental data is based on Apr-Sep means, although this is adjustable. ; Large-scale composite means are not done here. ; ; The MXD sites are limited first by specifying a year for which each site ; used must have had data for (e.g., specifying 1700 would use only sites ; with data back to at least 1700). After this selection, the available ; sites are sorted (on a region-by-region basis) into descending order ; according to their correlation with local Apr-Sep temperature (already ; computed, using the period 1881-1984 where data is available). Calculations ; are repeated using gradually more and more data. ; ; Comparison is made by computing correlations over the period 1881-1984, ; both filtered and unfiltered. ; ;----------------------------------------------------------------- dotem=0 ; 0=use CRUTEM1, 1=use CRUTEM2v if dotem eq 0 then begin fntem='lat_jones_18512000.mon.nc' savtem='' endif else begin fntem='crutem2v_18512001.mon.nc' savtem='ct2v_' endelse ; trv=0 ; selects tree-ring-variable: 0=MXD 1=TRW case trv of 0: begin fnadd='mxd' iseas=18 ; use local r with Apr-Sep temp to select chronologies seasst=3 seasen=8 dthr=3 ; used to be 3! end 1: begin fnadd='trw' iseas=20 ; use local r with Jun-Aug temp to select chronologies seasst=5 seasen=7 dthr=2 ; used to be 2! end endcase titadd=strupcase(fnadd) ; ; Specify options ; fixedyr=1400 ; year of fixed grid (1950, 1988, 1800, 1700,...,1400) rbarper=[1881,1960] ; period for computing correlation matrix corrper=[1881,1960] ; period over which comparison with temperature is made doreg=0 ; 0=use co-located temperature only, 1=use full region cutr=0.22 ; a 'best' series is based on sites with local r >= cutr ; ; Identify analysis period: ; For speed, analyse the minimum possible. Must cover all of rbarper and ; corrper, plus a bit more either side of corrper for filtering. Does not ; need to cover the fixedyr. ; analper=[1860,2000] anallen=analper(1)-analper(0)+1 ; ; First read in the MXD data ; print,'Reading '+titadd+' data' restore,filename='all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr ; ; Next read in the regional breakdown ; print,'Reading MXD regions' restore,filename='reg_mxdlists.idlsave' ; ntree,treelist,nreg,regname ; ; Next read in the correlations between MXD and local climate ; print,'Reading '+titadd+' local correlations' restore,filename='datastore/'+savtem+fnadd+'_moncorr.idlsave' ; allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2 ; ; Read in the instrumental temperatures for the required period ; print,'Reading temperatures' print,fntem ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+fntem) tsmon=crunc_rddata(ncid,tst=[analper(0),0],tend=[analper(1),11],grid=gtemp) ncdf_close,ncid nmon=12 ntemp=gtemp.nt nyrtemp=ntemp/nmon yrtemp=reform(gtemp.year,nmon,nyrtemp) yrtemp=reform(yrtemp(0,*)) ; ; Read in the regional breakdown of grid boxes ; print,'Reading temperature regions' restore,filename='reg_mxdboxes.idlsave' ; nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp ; ; Define arrays for storing all regional results ; z=!values.f_nan maxntree=max(ntree) nsites=fltarr(nreg) ; number of sites with data back to fixed yr regr=fltarr(maxntree,nreg,3)*z ; regional raw high low correls as extra series locr=fltarr(maxntree,nreg)*z ; local correlations sorted into descending order bestx=fltarr(maxntree,nreg)*z ; location of sites making best raw correlation besty=fltarr(maxntree,nreg)*z bestmxd=fltarr(anallen,nreg,3)*z ; MXD series: best raw regr; cutoff locr; all regtemp=fltarr(anallen,nreg)*z ; regional temperature timeseries rawtemp=fltarr(anallen,nreg)*z ; regional temps unadjusted for variance changes cutsites=fltarr(nreg) cutx=fltarr(maxntree,nreg) cuty=fltarr(maxntree,nreg) ; ; Now analyse each region in turn ; multi_plot,nrow=6,ncol=3,layout='large',/landscape loadct,39 def_1color,20,color='red' if !d.name eq 'X' then window,xsize=1000,ysize=850 for ireg = 0 , nreg-1 do begin print,regname(ireg) ; ; Extract monthly instrumental data for current region ; ; Choose mask (co-located data only, or all temperatures in region) if doreg eq 0 then fdmask=reform(boxlists(*,*,ireg)) $ else fdmask=reform(boxregs(*,*,ireg)) ; Count boxes to use blist=where(finite(fdmask),nbox) ; Define array to hold monthly grid box timeseries from region regts=fltarr(nbox,ntemp) ; Fill the array for imon = 0 , ntemp-1 do begin fd1mon=reform(tsmon(*,*,imon)) regts(*,imon)=fd1mon(blist) endfor ; Separate months from years regts=reform(regts,nbox,nmon,nyrtemp) ; Compute seasonal average (April-September) regseas=mkseason(regts,seasst,seasen,datathresh=dthr) ; Define an array to hold latitude weighting locy=fix(blist/g.nx) regwt=cos(g.y(locy)*!dtor) ; Compute a regional average (make second adjustment for input variance!) meanlat=var_adjust(regseas,yrtemp,weight=regwt,/td_rbar,$ refperiod=rbarper,/no_anomaly,td_sdev=1,rbar=rbar,corrmat=tempcm,$ /mkrbar,/mkcorr) ; For WNA and NWNA, early values are contaminated by warm anomalies, probably ; due to the use of north facing walls instead of Stephenson screens. Must ; remove the contaminated data. if ireg eq 6 then meanlat(where(yrtemp le 1871))=!values.f_nan if ireg eq 7 then meanlat(where(yrtemp le 1886))=!values.f_nan ; Need to store a regional series in oC. meanlat is adjusted but still in ; oC (hopefully), but represents a single series. Need to multiply by the ; time-constant rbar to get regional variance. klper=where((yrtemp ge rbarper(0)) and (yrtemp le rbarper(1))) constrbar=mkrbar(regseas(*,klper)) regtemp(*,ireg)=meanlat(*)*sqrt(constrbar) mknormal,meanlat,yrtemp,refperiod=rbarper ; ; Also compute unadjusted regional mean series ; wt=fltarr(nbox,nyrtemp) for iyr = 0 , nyrtemp-1 do wt(*,iyr)=regwt(*) totval=total(regseas*wt,1,/nan) totnum=total(float(finite(regseas))*wt,1,/nan) rawtemp(*,ireg)=totval/totnum ; ; First select sites in region ; nx=ntree(ireg) print,'Maximum number of sites',nx tl=treelist(0:nx-1,ireg) ; regmxd=mxd(*,tl) regfra=fraction(*,tl) regcor=reform(allr(tl,iseas,1)) ; iseas=18/20=Apr-Sep/Jun-Aug 1=temperature xxx=statlon(tl) yyy=statlat(tl) ; ; Select those sites with data back to required year ; iyr=where(timey eq fixedyr) kl=where(finite(regmxd(iyr,*)),nx) print,'Fixed grid number of sites',nx nsites(ireg)=nx if nx gt 0 then begin regmxd=regmxd(*,kl) regfra=regfra(*,kl) regcor=regcor(kl) xxx=xxx(kl) yyy=yyy(kl) ; ; Reduce time dimension to minimum required ; kl2=where((timey ge analper(0)) and (timey le analper(1)),nyr) regmxd=regmxd(kl2,*) regfra=regfra(kl2,*) yrmxd=timey(kl2) ; ; Transpose to give space then time ; regmxd=transpose(regmxd) regfra=transpose(regfra) ; ; Compute correlation matrix of all retained sites ; fullrbar=mkrbar(regmxd,grandr=corrmat) ; ; Identify decreasing order of local skill ; slist=reverse(sort(regcor)) locr(0:nx-1,ireg)=regcor(slist) ; ; Identify those with significant local correlations (approx.) ; dummy=where(regcor ge cutr,n2cut) cutsites(ireg)=n2cut ; !p.multi(0)=0 pause plot,regcor(slist),/xstyle,title=regname(ireg),psym=-1 ; ; Compute regional series for a range of subsets (based on skill sort) ; print,'Number of sites being used:' bestr=-1. for iuse = 0 , nx-1 do begin print,iuse,format='($,I4)' ; ; Extract skillful subset and subset of the correlation matrix ; submxd=regmxd(slist(0:iuse),*) subfra=regfra(slist(0:iuse),*) subcorr=corrmat(slist(0:iuse),*) subcorr=subcorr(*,slist(0:iuse)) subcorr=reform(subcorr,iuse+1,iuse+1) ; make sure it is 2D ; ; Compute weighted regional mean with appropriate variance adjustment ; meanmxd=var_adjust(submxd,yrmxd,weight=subfra,corrmat=subcorr,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,yrmxd,refperiod=rbarper ; ; Correlate MXD with temperature ; r3=mkcorrelation(meanmxd,meanlat,yrmxd,filter=10,refperiod=corrper,$ datathresh=10) regr(iuse,ireg,*)=r3(*) if (iuse eq 0) or (r3(0) gt bestr) then begin bestr=r3(0) bestmxd(*,ireg,0)=meanmxd(*) bestx(0:iuse,ireg)=xxx(slist(0:iuse)) besty(0:iuse,ireg)=yyy(slist(0:iuse)) endif if iuse+1 eq n2cut then begin bestmxd(*,ireg,1)=meanmxd(*) cutx(0:iuse,ireg)=xxx(slist(0:iuse)) cuty(0:iuse,ireg)=yyy(slist(0:iuse)) endif if iuse eq nx-1 then bestmxd(*,ireg,2)=meanmxd(*) ; if (iuse mod 5) eq 0 then begin plot,yrmxd,meanmxd,/xstyle,title=string(iuse)+string(r3(0)) oplot,yrtemp,meanlat,color=20 endif ; endfor print ; endif ; endfor ; ; Now save results ; timey=yrmxd if doreg eq 0 then fnadd2='' else fnaddi2='full' ; save,filename='datastore/'+savtem+fnadd2+'regbest_'+fnadd+'_fixed'+$ string(fixedyr,format='(I4)')+'.idlsave',$ nreg,regname,nsites,regr,locr,bestx,besty,bestmxd,regtemp,timey,$ fixedyr,corrper,doreg,maxntree,cutx,cuty,cutsites,seasst,seasen,rawtemp ; end