; ; Creates regional timeseries based on a limited sample of MXD series, ; variance adjusted with a time-dependent rbar, and compares them with ; regional instrumental records which are also computed here and variance ; adjusted. The instrumental records are for each month (including previous ; year) and from temperature and precipitation. The comparison is made via ; correlations and partial correlations. ; ; 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. Large-scale composite means are not done here. ; ; The MXD sites can be limited by specifying a minimum correlation that ; they should have with local Apr-Sep temperature (already computed, using ; the period 1881-1984 where data is available). ; ; Comparison is made by computing correlations over a specified period, ; both filtered and unfiltered. ; ; CAN SPECIFY IN ADVANCE corrper=[xxxx,xxxx] ; AND addsave='_early' or '_late' ; ;----------------------------------------------------------------- ; trv=1 ; selects tree-ring-variable: 0=MXD 1=TRW case trv of 0: begin fnadd='mxd' iseas=18 ; use local r with Apr-Sep temperature for chronology select end 1: begin fnadd='trw' iseas=20 ; use local r with Jun-Aug temperature for chronology select end endcase titadd=strupcase(fnadd) ; ; Specify options ; if n_elements(corrper) ne 2 then $ corrper=[1881,1984] ; period over which comparison with temperature is made rbarper=[1881,1960] ; period for computing correlation matrix 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 if n_elements(addsave) eq 0 then addsave='' ; ; 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. ; All data (MXD, temp, and prec) must be this length (either truncated or ; padded with missing data). ; analper=[1860,1994] anallen=analper(1)-analper(0)+1 nmon=12 ; ; First read in the MXD data (and truncate as necessary) ; print,'Reading '+titadd+' data' restore,filename='all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr kl=where((timey ge analper(0)) and (timey le analper(1)),nyr) mxd=mxd(kl,*) fraction=fraction(kl,*) timey=timey(kl) ; ; 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=fnadd+'_moncorr.idlsave' ; allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2 locr=allr ; ; Read in the instrumental land temperatures for the required period ; print,'Reading temperatures' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_jones_18511998.mon.nc') tmon=crunc_rddata(ncid,tst=[analper(0),0],tend=[analper(1),11],grid=gtemp) ncdf_close,ncid ntemp=gtemp.nt nyrtemp=ntemp/nmon if nyrtemp ne nyr then message,'Temperature and MXD data are different lengths' yrtemp=reform(gtemp.year,nmon,nyrtemp) yrtemp=reform(yrtemp(0,*)) if total(abs(timey-yrtemp)) ne 0 then message,'Year do not match' tmon=reform(tmon,gtemp.nx,gtemp.ny,nmon,nyr) ; ; Read in the instrumental land precipitation for the required period, and ; pad with missing data ; print,'Reading precipitation' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/precip_new_19011995.mon.nc') pmon1=crunc_rddata(ncid,tst=[1901,0],tend=[analper(1),11],grid=gprec) ncdf_close,ncid nprec=gprec.nt nyrprec=nprec/nmon npad=(nyr-nyrprec)*nmon pmon=fltarr(gprec.nx,gprec.ny,nyr*nmon) pmon(*,*,0:npad-1)=!values.f_nan pmon(*,*,npad:(nyr*nmon)-1)=pmon1(*,*,*) pmon1=0 ; free up the memory now we've finished with this variable pmon=reform(pmon,gprec.nx,gprec.ny,nmon,nyr) ; ; Now define composite regions ; ncomp=4 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) ; ; Define arrays ; thalf=10. ; filter for smoothed series, lower the better cos of sample size ncorr=21 ; 16 months (JJASOND-JFMAMJJAS) & 5 seasons (O-S O-M AS MA JJA) nvar=2 ; 0=MXD,precip, 1=MXD,temperature z=!values.f_nan allr=fltarr(ncomp,ncorr,nvar)*z ; correlations with P or T lowr=fltarr(ncomp,ncorr,nvar)*z ; low freq correlation with P or T allp=fltarr(ncomp,ncorr,nvar)*z ; partial correlations with P or T lowp=fltarr(ncomp,ncorr,nvar)*z ; low freq partial correlations with P or T compt=fltarr(nmon,nyr,ncomp)*z compp=fltarr(nmon,nyr,ncomp)*z compm=fltarr(nyr,ncomp)*z varname=['Prec','Temp'] corrname=['J','J','A','S','O','N','D','J','F','M','A','M','J','J','A','S',$ 'Oct-Sep','Oct-Mar','Apr-Sep','May-Aug','Jun-Aug'] ; ; Read in the regional breakdown of grid boxes ; print,'Reading instrumental regions' restore,filename='reg_mxdboxes.idlsave' ; nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp ; ; Make and analyse each composite regional mean in turn ; for icomp = 0 , ncomp-1 do begin print,compname(icomp) ; ; Extract monthly temperature data for current COMPOSITE region ; if icomp ne 3 then begin ; 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 an array to hold latitude weighting locy=fix(blist/gtemp.nx) regwt=cos(gtemp.y(locy)*!dtor) wt=fltarr(nbox,nyr) for iyr = 0 , nyr-1 do wt(*,iyr)=regwt(*) ; Repeat for each month for imon = 0 , nmon-1 do begin print,'Regional temperature for month',imon ; Define array to hold monthly grid box timeseries from region regts=fltarr(nbox,nyr) ; Fill the array for iyr = 0 , nyr-1 do begin fd1mon=reform(tmon(*,*,imon,iyr)) regts(*,iyr)=fd1mon(blist) endfor ; Compute regional average (latitude-weighted, but no variance adjustment) totval=total(regts*wt,1,/nan) totnum=total(float(finite(regts))*wt,1,/nan) meanlat=totval/totnum mknormal,meanlat,timey,refperiod=rbarper compt(imon,*,icomp)=meanlat(*) endfor endif else begin ; ; 'NH' region is the northern hemisphere >20N land series ; ; First extract >20N rows kl=where(gtemp.y ge 20.,nynh) ylat=gtemp.y(kl) tsnorth=reform(tmon(*,kl,*,*),gtemp.nx,nynh,nmon*nyr) ; Compute latitude-weighted mean nhmon=globalmean(tsnorth,ylat) ; Separate into monthly series nhmon=reform(nhmon,nmon,nyr) for imon = 0 , nmon-1 do begin onets=reform(nhmon(imon,*)) mknormal,onets,timey,refperiod=rbarper compt(imon,*,icomp)=onets(*) endfor endelse ; ; Extract monthly precipitation data for current COMPOSITE region ; if icomp ne 3 then begin ; Choose mask (co-located data only, or all precipitation 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 an array to hold latitude weighting locy=fix(blist/gprec.nx) regwt=cos(gprec.y(locy)*!dtor) wt=fltarr(nbox,nyr) for iyr = 0 , nyr-1 do wt(*,iyr)=regwt(*) ; Repeat for each month for imon = 0 , nmon-1 do begin print,'Regional precipitation for month',imon ; Define array to hold monthly grid box timeseries from region regts=fltarr(nbox,nyr) ; Fill the array for iyr = 0 , nyr-1 do begin fd1mon=reform(pmon(*,*,imon,iyr)) regts(*,iyr)=fd1mon(blist) endfor ; Compute regional average (latitude-weighted, but no variance adjustment) totval=total(regts*wt,1,/nan) totnum=total(float(finite(regts))*wt,1,/nan) meanlat=totval/totnum mknormal,meanlat,timey,refperiod=rbarper compp(imon,*,icomp)=meanlat(*) endfor endif else begin ; ; 'NH' region is the northern hemisphere >20N land series ; ; First extract >20N rows kl=where(gprec.y ge 20.,nynh) ylat=gprec.y(kl) tsnorth=reform(pmon(*,kl,*,*),gprec.nx,nynh,nmon*nyr) ; Compute latitude-weighted mean nhmon=globalmean(tsnorth,ylat) ; Separate into monthly series nhmon=reform(nhmon,nmon,nyr) for imon = 0 , nmon-1 do begin onets=reform(nhmon(imon,*)) mknormal,onets,timey,refperiod=rbarper compp(imon,*,icomp)=onets(*) endfor endelse ; ; For the MXD, we need to do each sub-region individually, before composite ; averaging ; if icomp ne 3 then begin 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(locr(tl,iseas,1)) ; 18=Apr-Sep 1=temperature ; ; K keep only those whose local r with Apr-Sep temperature is good enough ; kskill=where(regcor ge cutr,nskill) if nskill gt 0 then begin nx=nskill regmxd=regmxd(*,kskill) regfra=regfra(*,kskill) ; ; Transpose to give space then time ; regmxd=transpose(regmxd) regfra=transpose(regfra) ; ; Compute correlation matrix and full rbar ; fullrbar=mkrbar(regmxd,grandr=corrmat) ; ; Compute weighted regional mean with appropriate variance adjustment ; meanmxd=var_adjust(regmxd,timey,weight=regfra,corrmat=corrmat,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,timey,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 ; ; Repeat for next region ; endfor ; ; Now compute an EPS-weighted, variance-adjusted composite regional mean ; meanmxd=var_adjust(indts,timey,weight=indeps,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,timey,refperiod=rbarper compm(*,icomp)=meanmxd(*) ; endif else begin compm(*,icomp)=compm(*,icomp-1) ; for MXD, NH=ALL endelse ; ; Repeat for precipitation then for temperature ; for ivar = 0 , nvar-1 do begin if ivar eq 0 then begin var12=reform(compp(*,*,icomp)) ; variable to analyse con12=reform(compt(*,*,icomp)) ; variable to hold constant endif else begin var12=reform(compt(*,*,icomp)) con12=reform(compp(*,*,icomp)) endelse ; ; Now correlate the MXD with monthly data ; (1) With previous year's months (June to December) ; icorr=0 for imon = 5 , nmon-1 do begin ; Full correlation var1=reform(var12(imon,*)) var1=shift(var1,1) var1(0)=!values.f_nan r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=reform(con12(imon,*)) con1=shift(con1,1) con1(0)=!values.f_nan r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 endfor ; ; (2) With current year's months (January to September) ; for imon = 0 , 8 do begin ; Full correlation var1=reform(var12(imon,*)) r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=reform(con12(imon,*)) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 endfor ; ; (3) With Oct to Sep annual mean and Oct to Mar winter mean ; ; Full correlation var1=mkseason(var12,9,8,datathresh=10) r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=mkseason(con12,9,8,datathresh=10) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 ; ; Full correlation var1=mkseason(var12,9,2,datathresh=5) r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=mkseason(con12,9,2,datathresh=5) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 ; ; (4) April-September, May-August, and July-August means ; ; Full correlation var1=mkseason(var12,3,8,datathresh=5) r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=mkseason(con12,3,8,datathresh=5) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 ; ; Full correlation var1=mkseason(var12,4,7,datathresh=4) r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=mkseason(con12,4,7,datathresh=4) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 ; ; Full correlation var1=mkseason(var12,5,7,datathresh=2) ; NEWEST!!! USE JUNE-AUGUST ; var1=mkseason(var12,6,7,datathresh=2) ; NEW!!! USE JULY-AUGUST ; var1=mkseason(var12,5,6,datathresh=2) ; OLD!!! USE JUNE-JULY r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allr(icomp,icorr,ivar)=r3(0) lowr(icomp,icorr,ivar)=r3(2) ; Partial correlation con1=mkseason(con12,5,7,datathresh=2) r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) allp(icomp,icorr,ivar)=r3(0) lowp(icomp,icorr,ivar)=r3(2) icorr=icorr+1 ; ; Repeat for next variable ; endfor ; ; Repeat for next region ; endfor print ; ; Save all the data ; save,filename='comp'+fnadd+'_moncorr'+addsave+'.idlsave',$ allr,ncorr,nvar,ncomp,compname,varname,corrname,$ allp,lowr,lowp,thalf,compt,compp,compm ; end