; ; Analyses MXD data in terms of their correlation with local monthly ; temperature and precipitation. All data up until the end of 1984 ; is used (more recent data is not used, because of differing spatial ; temporal coverage. Also computes partial correlations, and correlations ; between a moisture deficit value and MXD. ; ALSO computes correlations after smoothing (to identify the decline). ; OPTIONALLY can replace precipitation with cloud cover! ; ;---------------------------------------------------- ; 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 ; docld=0 ; 0=T & precip, 1=T & cloud trv=1 ; selects tree-ring-variable: 0=MXD 1=TRW 2=MXD-TRW case trv of 0: fnadd='mxd' 1: fnadd='trw' 2: fnadd='mxd-trw' endcase titadd=strupcase(fnadd) ; ; Get chronology locations ; print,'Reading '+titadd+' data' if trv eq 2 then begin restore,filename='alltrw.idlsave' trw=mxd restore,filename='allmxd.idlsave' mxd=mxd-trw endif else begin restore,filename='all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr endelse ; ; Now read in the 1961-90 monthly means of precip and temperature ; print,'Reading precipitation baseline' if docld then cldfn='cld_new_19311990' else cldfn='precip_new_19611990' restore,'/cru/u2/f055/data/obs/grid/surface/'+cldfn+'.idlsave' ; g,nmon,ltm5 pre6190=ltm5 print,'Reading temperature baseline' restore,'/cru/u2/f055/data/obs/grid/surface/lat_new_19611990.idlsave' ; g,nmon,ltm5 lat6190=ltm5 ; ; Now read in the land precipitation dataset. ; Although there is some missing data in Mark New's precip anomalies, all ; MXD boxes have sufficient data, so do not have to use surrounding boxes. ; print,'Reading precipitation data' if docld then cldfn='cldanom' else cldfn='precip' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+cldfn+'_new_19011995.mon.nc') pmon=crunc_rddata(ncid,tst=[1901,0],tend=[1984,11],grid=g) ncdf_close,ncid ; ; Now read in the land air temperature dataset. ; print,'Reading temperature data' print,fntem ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/'+fntem) tmon=crunc_rddata(ncid,tst=[1881,0],tend=[1984,11],grid=gtemp) ncdf_close,ncid ; ; The precipitation series will be padded with missing values at its ; beginning, to make it the same length as the temperature. Find pad length, ; and also find section length of MXD data to extract. ; nmon=12 npad=(gtemp.nt-g.nt)/nmon ppad=replicate(!values.f_nan,nmon*npad) klmxd=where((timey ge 1881) and (timey le 1984),nyrm) timey=timey(klmxd) ; ; ALTHOUGH THE FULL PERIOD IS ALWAYS READ (TO ENSURE CONSISTENT EVALUATION ; OF WHICH TEMPERATURE BOXES DO NOT HAVE 300 OR MORE MONTHS OF NON-MISSING ; 1881-1964 DATA) THE FOLLOWING LINE ALLOWS SELECTION OF A SMALLER SUB-PERIOD ; FOR ANALYSIS. ; klperiod=where((timey ge 1881) and (timey le 1984),nyrper) print,'LENGTH OF ANALYSIS PERIOD (YEARS) =',nyrper,' *************' ; ; 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 A-S M-A JJA) nvar=2 ; 0=MXD,precip, 1=MXD,temperature nyrp=g.nt/nmon nyrt=gtemp.nt/nmon if (nyrp+npad) ne nyrm then message,'Ooops!!!' if nyrt ne nyrm then message,'Ooops!!!' nyr=nyrm z=!values.f_nan allr=fltarr(nchron,ncorr,nvar)*z ; correlations with P or T lowr=fltarr(nchron,ncorr,nvar)*z ; low freq correlation with P or T moir=fltarr(nchron)*z ; corr with moisture deficit moir2=fltarr(nchron)*z ; corr Apr-Sep T vs moist deficit allp=fltarr(nchron,ncorr,nvar)*z ; partial correlations with P or T lowp=fltarr(nchron,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'] ; ; Analyse each chronology in turn ; runlen=0. runnum=0. for i = 0 , nchron-1 do begin print,i,format='($,I4)' ; ; Extract MXD timeseries ; mxd1=reform(mxd(*,i)) mxd1=mxd1(klmxd) ; ; Locate grid box containing site and extract monthly timeseries ; x=statlon(i) y=statlat(i) statval=[1.] fdout=gridit(g.nx,g.ny,g.x,g.y,x,y,statval) loclist=where(finite(fdout),nloc) if nloc ne 1 then message,'Ooops!!!' locx=loclist(0) mod g.nx locy=fix(loclist(0)/g.nx) ; ; Simple single box mean for precipitation, but pad with missing ; pre12=[ppad,reform(pmon(locx,locy,*))] pre12=reform(pre12,nmon,nyr) ; ; For temperature, if less than 300 months of 1881-1964 data use an ; average of the nine surrounding boxes (variance adjusted etc.) ; tem12=reform(tmon(locx,locy,*),nmon,nyr) kl=where(gtemp.year le 1964) dummy=where(finite(tem12(*,kl)),ngot) if ngot lt 300 then begin print ; extract timeseries from 9 surrounding boxes (NB. zonally cyclic) locy=[locy-1,locy,locy+1] locx=[locx-1,locx,locx+1] mod g.nx neglist=where(locx lt 0,nneg) if nneg gt 0 then locx(neglist)=locx(neglist)+g.nx temboxes=tmon(locx,*,*) temboxes=temboxes(*,locy,*) temboxes=reform(temboxes,9,nmon,nyr) ; average them and adjust variance on a month by month basis ; I've altered mkeffective to adjust variance to that of a single ; input timeseries for imon = 0 , nmon-1 do begin print,imon,format='($,I4)' tsin=reform(temboxes(*,imon,*)) tsout=var_adjust(tsin,/no_anomaly,/td_rbar,/mkrbar,/mkcorr) tem12(imon,*)=tsout(*) endfor print endif ; ; Compute the monthly PET year by year ; ; Define PET array monpet=fltarr(nmon,nyr)*!values.f_nan ; Get absolute temperature & precip for 1961-90 tabs=reform(lat6190(locx,locy,*)) pabs=reform(pre6190(locx,locy,*)) ; Repeat for each year for j = 0 , nyr-1 do begin onet=reform(tem12(*,j))+tabs(*) if finite(total(onet)) then begin onepet=thornthwaite(onet,g.y(locy)) monpet(*,j)=onepet(*) endif endfor ; ; Compute April to July sum of PET ; j2apet=monpet(3:6,*) j2apet=total(j2apet,1) ; ; Compute November to July precipitation (TOTAL not MEAN) ; pre1=mkseason(pre12,10,6,datathresh=8) ; need 8 out of 9 for mean pre1=pre1*9. ; convert to 9 month total ; ; Moisture deficit at end of July ; mdef=j2apet-pre1 ; save,filename='datastore/'+savtem+'example_chron'+string(i,format='(I3.3)')+'.idlsave',$ mxd1,pre12,tem12,nmon,nyr,timey,i,location ; ; Extract period for analysis ; mxd1=mxd1(klperiod) pre12=pre12(*,klperiod) tem12=tem12(*,klperiod) mdef=mdef(klperiod) ; ; Correlate moisture with MXD ; kl=where(finite(mxd1+mdef),nkeep) if nkeep gt 5 then begin filter_cru,10.,/nan,tsin=mxd1,tshigh=th1 filter_cru,10.,/nan,tsin=mdef,tshigh=th2 ; moir(i)=correlate(mxd1(kl),mdef(kl)) ; do raw series moir(i)=correlate(th1(kl),th2(kl)) ; do high-pass series endif ; ; Repeat for precipitation then for temperature ; for ivar = 0 , nvar-1 do begin if ivar eq 0 then begin var12=pre12 ; variable to analyse con12=tem12 ; variable to hold constant endif else begin var12=tem12 con12=pre12 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 pre1=reform(var12(imon,*)) pre1=shift(pre1,1) pre1(0)=!values.f_nan kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=reform(con12(imon,*)) con1=shift(con1,1) con1(0)=!values.f_nan kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 endfor ; ; (2) With current year's months (January to September) ; for imon = 0 , 8 do begin ; Full correlation pre1=reform(var12(imon,*)) kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=reform(con12(imon,*)) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 endfor ; ; (3) With Oct to Sep annual mean and Oct to Mar winter mean ; ; Full correlation pre1=mkseason(var12,9,8,datathresh=10) kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=mkseason(con12,9,8,datathresh=10) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 ; ; Full correlation pre1=mkseason(var12,9,2,datathresh=5) kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=mkseason(con12,9,2,datathresh=5) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 ; ; (4) April-September, May-August, and June-August means ; ; Full correlation pre1=mkseason(var12,3,8,datathresh=5) kl=where(finite(mxd1+pre1),nkeep) runlen=runlen+nkeep runnum=runnum+1. if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=mkseason(con12,3,8,datathresh=5) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 ; ; Moisture deficit at end of July and correlate with Apr-Sep temp ; if ivar eq 1 then begin kl=where(finite(pre1+mdef),nkeep) if nkeep gt 5 then moir2(i)=correlate(pre1(kl),mdef(kl)) endif ; ; Full correlation pre1=mkseason(var12,4,7,datathresh=4) kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=mkseason(con12,4,7,datathresh=4) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 ; ; Full correlation pre1=mkseason(var12,5,7,datathresh=2) ; NEWEST!!! USE JUNE-AUGUST ; pre1=mkseason(var12,6,7,datathresh=2) ; NEW!!! USE JULY-AUGUST ; pre1=mkseason(var12,5,6,datathresh=2) ; OLD!!! USE JUNE-JULY kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin allr(i,icorr,ivar)=correlate(mxd1(kl),pre1(kl)) filter_cru,thalf,/nan,tsin=mxd1,tslow=tl1 filter_cru,thalf,/nan,tsin=pre1,tslow=tl2 lowr(i,icorr,ivar)=correlate(tl1(kl),tl2(kl)) endif ; Partial correlation con1=mkseason(con12,5,7,datathresh=2) kl=where(finite(mxd1+pre1+con1),nkeep) if nkeep gt 5 then begin allp(i,icorr,ivar)=p_correlate(pre1(kl),mxd1(kl),reform(con1(kl),1,nkeep)) filter_cru,thalf,/nan,tsin=pre1,tslow=tl1 filter_cru,thalf,/nan,tsin=mxd1,tslow=tl2 filter_cru,thalf,/nan,tsin=con1,tslow=tl3 lowp(i,icorr,ivar)=p_correlate(tl1(kl),tl2(kl),reform(tl3(kl),1,nkeep)) endif ; icorr=icorr+1 ; ; Repeat for next variable ; endfor ; ; Repeat for next site ; endfor print ; runlen=runlen/runnum print,'Average Apr-Sep vs '+titadd+' correlation sample length=',runlen ; ; Save all the data ; if docld eq 0 then fnsave='datastore/'+savtem+fnadd+'_moncorr.idlsave' $ else fnsave='datastore/'+savtem+fnadd+'_moncorr_cldtem.idlsave' print,fnsave save,filename=fnsave,$ allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,$ allp,moir,moir2,lowr,lowp,thalf ; end