; ; 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=3 ; selects tree-ring-variable: 0=MXD 1=TRW 3=TRW|MXD case trv of 0: begin fnadd='mxd' lseas=18 ; select chronologies on the basis of local Apr-Sep corr end 1: begin fnadd='trw' lseas=20 ; select chronologies on the basis of local Jun-Aug corr end 3: begin fnadd='trwmxd' lseas=-1 ; do not subset the chronologies at this stage end endcase titadd=strupcase(fnadd) ; ; Specify options ; if n_elements(corrper) ne 2 then $ corrper=[1881,1960] ; 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' if fnadd eq 'trwmxd' then begin restore,filename='alltrw.idlsave' trw=mxd restore,filename='allmxd.idlsave' endif else begin restore,filename='all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr endelse kl=where((timey ge analper(0)) and (timey le analper(1)),nyr) mxd=mxd(kl,*) if fnadd eq 'trwmxd' then trw=trw[kl,*] fraction=fraction(kl,*) timey=timey(kl) ; ; Next read in the regional breakdown ; print,'Reading '+titadd+' 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 dummy=where(locr lt 0.22,nlnl) & print,'sites failed=',nlnl ; ; Read in the instrumental land temperatures for the required period ; print,'Reading temperatures' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/crutem2v_18512001.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) ; ; 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(nreg,ncorr,nvar)*z ; correlations with P or T lowr=fltarr(nreg,ncorr,nvar)*z ; low freq correlation with P or T allp=fltarr(nreg,ncorr,nvar)*z ; partial correlations with P or T lowp=fltarr(nreg,ncorr,nvar)*z ; low freq partial correlations with P or T 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 regional mean in turn ; for ireg = 0 , nreg-1 do begin print,regname(ireg) ; ; Extract monthly temperature data for current region ; ; Define array to hold final monthly timeseries latmon=fltarr(nmon,nyr) ; 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 an array to hold latitude weighting locy=fix(blist/gtemp.nx) regwt=cos(gtemp.y(locy)*!dtor) ; 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 a regional average (adjust: sample size & optionally input var) latone=var_adjust(regts,timey,weight=regwt,/td_rbar,$ refperiod=rbarper,/no_anomaly,/td_sdev) mknormal,latone,timey,refperiod=rbarper ; For WNA and NWNA, early values are contaminated by warm anomalies, prob. ; due to the use of north facing walls instead of Stephenson screens. Must ; remove the contaminated data. if ireg eq 6 then latone(where(timey le 1871))=!values.f_nan if ireg eq 7 then latone(where(timey le 1886))=!values.f_nan ; Store timeseries and repeat for next month latmon(imon,*)=latone(*) endfor ; ; Extract monthly precipitation data for current region ; ; Define array to hold final monthly timeseries premon=fltarr(nmon,nyr) ; Choose mask (co-located data only, or all precipitation 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 an array to hold latitude weighting locy=fix(blist/gtemp.nx) regwt=cos(gprec.y(locy)*!dtor) ; 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 a regional average (adjust: sample size & optionally input var) preone=var_adjust(regts,timey,weight=regwt,/td_rbar,$ refperiod=rbarper,/no_anomaly,/td_sdev) mknormal,preone,timey,refperiod=rbarper ; Store timeseries and repeat for next month premon(imon,*)=preone(*) endfor ; ; Store the regional monthly instrumental series for later plotting ; save,filename='datastore/regobs_'+regname(ireg)+'.idlsave',$ nmon,nyr,timey,latmon,premon ; ; Now compute MXD regional timeseries ; ; First select sites in region ; nx=ntree(ireg) print,'Maximum number of sites',nx tl=treelist(0:nx-1,ireg) regmxd=mxd(*,tl) if fnadd eq 'trwmxd' then regtrw=trw[*,tl] regfra=fraction(*,tl) if lseas ge 0 then begin regcor=reform(locr(tl,lseas,1)) ; 18=Apr-Sep 1=temperature endif else begin regcor=replicate(1.,nx) ; select everything endelse dummy=where(locr lt 0.22,nlnl) & print,'sites failed=',nlnl ; ; Next keep only those whose local r with Apr-Sep temperature is good enough ; kskill=where(regcor ge cutr,nskill) help,regmxd,regcor,kskill if nskill gt 0 then begin nx=nskill regmxd=regmxd(*,kskill) if fnadd eq 'trwmxd' then regtrw=regtrw[*,kskill] regfra=regfra(*,kskill) ; ; Transpose to give space then time ; regmxd=transpose(regmxd) if fnadd eq 'trwmxd' then regtrw=transpose(regtrw) regfra=transpose(regfra) ; ; Compute weighted regional mean with appropriate variance adjustment ; meanmxd=var_adjust(regmxd,timey,weight=regfra,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meanmxd,timey,refperiod=rbarper ; if fnadd eq 'trwmxd' then begin meantrw=var_adjust(regtrw,timey,weight=regfra,$ refperiod=rbarper,/no_anomaly,/td_rbar) mknormal,meantrw,timey,refperiod=rbarper endif ; endif else begin meanmxd=fltarr(nyr)*!values.f_nan if fnadd eq 'trwmxd' then meantrw=fltarr(nyr)*!values.f_nan endelse ; ; Repeat for precipitation then for temperature ; for ivar = 0 , nvar-1 do begin if ivar eq 0 then begin var12=premon ; variable to analyse con12=latmon ; variable to hold constant endif else begin var12=latmon con12=premon endelse ; ; Now correlate the MXD with monthly data and seasonal data ; for icorr = 0 , ncorr-1 do begin ; Select month (possibly lagged) and season of instrumental data if icorr le 6 then begin imon=icorr+5 ; Full correlation var1=reform(var12(imon,*)) var1=shift(var1,1) var1(0)=!values.f_nan ; Partial correlation con1=reform(con12(imon,*)) con1=shift(con1,1) con1(0)=!values.f_nan endif if (icorr ge 7) and (icorr le 15) then begin imon=icorr-7 ; Full correlation var1=reform(var12(imon,*)) ; Partial correlation con1=reform(con12(imon,*)) endif if icorr ge 16 then begin stseas=[9,9,3,4,5] enseas=[8,2,8,7,7] dtseas=[10,5,5,3,2] iseas=icorr-16 var1=mkseason(var12,stseas[iseas],enseas[iseas],datathresh=dtseas[iseas]) con1=mkseason(con12,stseas[iseas],enseas[iseas],datathresh=dtseas[iseas]) endif ; Now do the correlations if fnadd eq 'trwmxd' then begin r3=mkpcorrelation(meantrw,var1,meanmxd,timey,filter=thalf,refperiod=corrper,$ datathresh=10) endif else begin r3=mkcorrelation(meanmxd,var1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) endelse allr(ireg,icorr,ivar)=r3(0) lowr(ireg,icorr,ivar)=r3(2) ; Now do the partial correlations if fnadd eq 'trwmxd' then begin if n_elements(con1) ne n_elements(meanmxd) then message,'Mismatch!!' reminf=transpose(reform([con1,meanmxd],nyr,2)) r3=mkp2correlation(meantrw,var1,reminf,timey,filter=thalf,refperiod=corrper,$ datathresh=10) endif else begin r3=mkpcorrelation(var1,meanmxd,con1,timey,filter=thalf,refperiod=corrper,$ datathresh=10) endelse allp(ireg,icorr,ivar)=r3(0) lowp(ireg,icorr,ivar)=r3(2) ; endfor ; ; Repeat for next variable ; endfor ; ; Repeat for next region ; endfor print ; ; Save all the data ; save,filename='reg'+fnadd+'_moncorr'+addsave+'.idlsave',$ allr,ncorr,nvar,nreg,regname,varname,corrname,$ allp,lowr,lowp,thalf ; end