; ; Creates composite regional timeseries based on a limited sample of MXD ; series, with regional series variance adjusted with a time-dependent rbar, ; and then averaged with EPS-weighting and themselves variance adjusted, and ; compares them with composite regional and near-hemispheric instrumental ; temperature series which are also computed here (BUT NOT variance adjusted ; since for the instrumental series I do not want to go through a 2-stage ; averaging process, and the number of boxes involved in a direct average ; is too many to variance adjust - and the adjustment would be extremely small ; anyway!). ; ; 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. ; In addition, the near-hemispheric mean of all LAND north of 20N is used. ; ; 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, all sites regardless ; of which sub-region they are in, that have correlations (with local Apr-Sep ; temperatures) >= CUTOFF are used, where CUTOFF is 0.6, 0.5, 0.4, ; 0.3, 0.22, 0.1, 0.0, -1. ; ; Comparison is made by computing correlations over the period 1881-1960, ; both filtered and unfiltered. ; ;----------------------------------------------------------------- ; dotem=1 ; 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 end 1: begin fnadd='trw' iseas=20 ; use local r with Jun-Aug temp to select chronologies seasst=5 seasen=7 dthr=2 end endcase titadd=strupcase(fnadd) ; ; Specify options ; fixedyr=1988 ; year of fixed grid 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.6,0.5,0.4,0.3,0.22,0.1,0.0,-1.] ; cutoffs for local correlations ncut=n_elements(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 ; ; Now define composite regions ; ncomp=3 compname=['HILAT','LOLAT','ALL','NH'] avgnum=[5,4,9] avgreg=intarr(ncomp,max(avgnum)) avgreg(0,0:avgnum(0)-1)=[0,2,3,7,8] avgreg(1,0:avgnum(1)-1)=[1,4,5,6] avgreg(2,0:avgnum(2)-1)=indgen(nreg) ; ; Next read in the correlations between MXD and local climate ; print,'Reading '+titadd+' local correlations' fnlocal='datastore/'+savtem+fnadd+'_moncorr.idlsave' print,fnlocal restore,filename=fnlocal ; 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 nsites=fltarr(nreg,ncut) ; no. sites with data back to fixed yr & local r compr=fltarr(ncomp+1,ncut,3)*z ; composite raw high low correls for cutoffs compmxd=fltarr(anallen,ncomp,2)*z ; MXD series: locr >= 0.22; all comptemp=fltarr(anallen,ncomp+1)*z ; composite temperature timeseries & NH ; ; Now analyse each region in turn ; multi_plot,nrow=4,ncol=2,layout='large',/landscape loadct,39 def_1color,20,color='red' def_1color,21,color='blue' if !d.name eq 'X' then window,xsize=1000,ysize=850 for icomp = 0 , ncomp-1 do begin print,compname(icomp) ; ; Extract monthly instrumental data for current COMPOSITE region ; ; Choose mask (co-located data only, or all temperatures in region) ireg=avgreg(icomp,0:avgnum(icomp)-1) if doreg eq 0 then fdmask=reform(boxlists(*,*,ireg)) $ else fdmask=reform(boxregs(*,*,ireg)) fdmask=total(finite(fdmask),3) ; if any regional mask has data, comp has too ; Count boxes to use blist=where(fdmask gt 0,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) wt=fltarr(nbox,nyrtemp) for iyr = 0 , nyrtemp-1 do wt(*,iyr)=regwt(*) ; Compute a regional average (latitude-weighted, but no variance adjustment) totval=total(regseas*wt,1,/nan) totnum=total(float(finite(regseas))*wt,1,/nan) meanlat=totval/totnum comptemp(*,icomp)=meanlat(*) mknormal,meanlat,yrtemp,refperiod=rbarper ; ; For the 'ALL' region, also compute the northern hemisphere >20N land series ; if icomp eq 2 then begin ; First extract >20N rows kl=where(g.y ge 20.) ylat=g.y(kl) tsnorth=tsmon(*,kl,*) ; Compute latitude-weighted mean nhmon=globalmean(tsnorth,ylat) ; Compute Apr-Sep mean nhmon=reform(nhmon,nmon,nyrtemp) nhseas=mkseason(nhmon,seasst,seasen,datathresh=dthr) comptemp(*,icomp+1)=nhseas(*) mknormal,nhseas,yrtemp,refperiod=rbarper endif ; ; Repeat regional mean for different subsets according to local skills ; !p.multi(0)=0 pause for icut = 0 , ncut-1 do begin print,icut,cutr(icut) ; ; For the MXD, we need to do each sub-region individually, before composite ; averaging ; ndoreg=avgnum(icomp) indts=fltarr(ndoreg,anallen)*!values.f_nan indeps=fltarr(ndoreg,anallen)*!values.f_nan for idoreg = 0 , ndoreg-1 do begin ireg=avgreg(icomp,idoreg) print,regname(ireg) ; ; 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)) ; 18/20=Apr-Sep/Jun-Aug 1=temperature xxx=statlon(tl) yyy=statlat(tl) ; ; Select those sites with local skill >= current cutoff ; kl=where(regcor ge cutr(icut),nkeep) print,'Skill-based subset',nkeep if nkeep gt 0 then begin regmxd=regmxd(*,kl) regfra=regfra(*,kl) regcor=regcor(kl) xxx=xxx(kl) yyy=yyy(kl) ; ; 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,icut)=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) print,nyr,' : THIS MUST BE NON-ZERO!' regmxd=regmxd(kl2,*) regfra=regfra(kl2,*) yrmxd=timey(kl2) ; ; Transpose to give space then time ; regmxd=transpose(regmxd) regfra=transpose(regfra) ; ; Compute correlation matrix and full rbar ; fullrbar=mkrbar(regmxd,grandr=corrmat) if nx eq 1 then fullrbar=1. ; ; Compute weighted regional mean with appropriate variance adjustment ; meanmxd=var_adjust(regmxd,yrmxd,weight=regfra,corrmat=corrmat,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,yrmxd,refperiod=rbarper ; ; Store series for later compositing ; indts(idoreg,*)=meanmxd(*) ; ; Compute time-dependent EPS for later weighting (need time-dependent ; sample size and full rbar) ; nsamp=float(total(finite(regmxd),1)) indeps(idoreg,*)=fullrbar / (fullrbar + (1.-fullrbar)/nsamp(*) ) ; endif endif ; ; Repeat for next region ; endfor ; ; Check in case we have no series at all (all failed the local r cutoff) ; if total(finite(indts)) eq 0 then begin domiss=1 meanmxd=fltarr(anallen)*!values.f_nan endif else begin domiss=0 ; ; Now compute an EPS-weighted, variance-adjusted composite regional mean ; meanmxd=var_adjust(indts,yrmxd,weight=indeps,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,yrmxd,refperiod=rbarper ; endelse ; ; Store composite regional mean series if required ; if cutr(icut) eq 0.22 then begin compmxd(*,icomp,0)=meanmxd(*) endif if cutr(icut) eq -1. then begin compmxd(*,icomp,1)=meanmxd(*) endif ; ; Correlate MXD with temperature ; if domiss eq 0 then begin r3=mkcorrelation(meanmxd,meanlat,yrmxd,filter=10,refperiod=corrper,$ datathresh=10) compr(icomp,icut,*)=r3(*) plot,yrmxd,meanmxd,/xstyle,title=string(cutr(icut))+' '+string(r3(0)) oplot,yrtemp,meanlat,color=20 endif else begin compr(icomp,icut,*)=!values.f_nan endelse ; ; If region is 'ALL' then show and correlate with NH temperature too ; if icomp eq 2 then begin oplot,yrtemp,nhseas,color=21 r3=mkcorrelation(meanmxd,nhseas,yrmxd,filter=10,refperiod=corrper,$ datathresh=10) compr(icomp+1,icut,*)=r3(*) endif ; ; Repeat for next skill cutoff ; endfor ; ; Repeat for next composite region ; endfor ; ; Now save results ; timey=yrmxd if doreg eq 0 then fnadd2='' else fnadd2='full' save,filename='datastore/'+savtem+fnadd2+'compbest_'+fnadd+'_fixed'+$ string(fixedyr,format='(I4)')+'.idlsave',$ ncomp,compname,nreg,regname,nsites,compr,compmxd,comptemp,timey,$ fixedyr,corrper,doreg,nchron,cutr,ncut,seasst,seasen ; end