; ; Plots scatter graphs of MXD versus Apr-Sep temperature to look for ; evidence of a nonlinear relationship. Only does it for the sites ; with the very best correlations with Apr-Sep, so that the scatter ; is small enough to see any nonlinearity. ; ;---------------------------------------------------- ; trv=0 ; selects tree-ring-variable: 0=MXD 1=TRW 2=MXD-TRW case trv of 0: fnadd='mxd' 1: fnadd='trw' endcase titadd=strupcase(fnadd) ; ; Read in the correlations with Apr-Sep temperature ; restore,filename=fnadd+'_moncorr.idlsave' ; Gets: allr (plus others) ; oner=reform(allr(*,18,1)) ; vs. Apr-Sep temperature cutr=0.5 kcut=where(oner gt cutr,ncut) print,'Number of sites to analyse=',ncut ; ; Get chronology locations ; print,'Reading '+titadd+' data' restore,filename='all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr ; ; Now read in the land air temperature dataset. ; print,'Reading temperature data' ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_jones_18512000.mon.nc') tmon=crunc_rddata(ncid,tst=[1881,0],tend=[1984,11],grid=g) ncdf_close,ncid ; nmon=12 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 ; nyrt=g.nt/nmon if nyrt ne nyrm then message,'Ooops!!!' nyr=nyrm ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=4,ncol=3,layout='large' if !d.name eq 'X' then begin window,ysize=850 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=11 endelse def_1color,10,color='mlgrey' ; ; Define bins to use ; binmin=-2.5 binmax=+4.1 binsize=0.2 nbin=(binmax-binmin)/binsize binx=findgen(nbin)*binsize+binmin biny=fltarr(nbin) bintot=fltarr(nbin) bin2tot=fltarr(nbin) binnum=fltarr(nbin) coefftot=fltarr(2) coeffnum=0. ; ; Analyse each chronology in turn ; for iii = 0 , ncut-1 do begin i=kcut(iii) print,i ; ; 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) ; ; 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(g.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 ; ; Extract period for analysis ; mxd1=mxd1(klperiod) tem12=tem12(*,klperiod) var12=tem12 ; ; April-September means ; ; Full correlation pre1=mkseason(var12,3,8,datathresh=5) kl=where(finite(mxd1+pre1),nkeep) if nkeep gt 5 then begin mxd1=mxd1(kl) pre1=pre1(kl) mknormal,mxd1 mknormal,pre1 r1=correlate(mxd1,pre1) acoeff=linfit(pre1,mxd1) coefftot(*)=coefftot(*)+acoeff(*) coeffnum=coeffnum+1. print,r1,min(pre1),max(pre1),min(mxd1),max(mxd1) pause plot,pre1,mxd1,psym=def_sym(15),$ /xstyle,xrange=[-2.5,4.1],xtitle='Normalised Apr-Sep temp',$ /ystyle,yrange=[-3.8,2.8],ytitle='Normalised '+titadd,$ xmargin=[6,0] oplot,!x.crange,[0,0],linestyle=1 oplot,!x.crange,!x.crange,linestyle=1 oplot,[0,0],!y.crange,linestyle=1 oplot,!x.crange,!x.crange*acoeff(1)+acoeff(0),thick=2 xyouts,1,-3,string(r1,format='(F5.2)') for ibin = 0 , nbin-1 do begin binlist=where((pre1 ge binx(ibin)) and (pre1 lt binx(ibin)+binsize),nl) if nl eq 0 then begin biny(ibin)=!values.f_nan endif else begin mval=total(mxd1(binlist)) m2val=total(mxd1(binlist)^2) biny(ibin)=mval/float(nl) bintot(ibin)=bintot(ibin)+mval bin2tot(ibin)=bin2tot(ibin)+m2val binnum(ibin)=binnum(ibin)+nl endelse endfor plots,binx+0.5*binsize,biny,psym=-def_sym(14),thick=4 endif ; ; Repeat for next site ; endfor ; ; Compute and plot mean binnned curve ; multi_plot,nrow=2 pause binmean=bintot/binnum bin2mean=bin2tot/binnum binsd=bin2mean-binmean^2 ml=where(binnum lt 4,nmiss) coeffmean=coefftot/coeffnum plot,binx+0.5*binsize,binmean,/nodata,$ /xstyle,xrange=[-2.5,4.1],xtitle='Normalised Apr-Sep temp',$ /ystyle,yrange=[-3.8,2.8],ytitle='Normalised '+titadd,$ xmargin=[6,0] oplot,binx+0.5*binsize,binmean,psym=-def_sym(14),thick=4,color=10 if nmiss gt 0 then binmean(ml)=!values.f_nan oplot,binx+0.5*binsize,binmean,psym=-def_sym(14),thick=4 oplot,!x.crange,[0,0],linestyle=1 oplot,!x.crange,!x.crange,linestyle=1 oplot,[0,0],!y.crange,linestyle=1 oplot,!x.crange,!x.crange*coeffmean(1)+coeffmean(0),thick=2 ; ; Now compute uncertainty ranges ; for ibin = 0 , nbin-1 do begin if finite(binmean(ibin)) then begin v=binnum(ibin)-1 t05=t_cvf(0.025,v) halfrange=t05*binsd(ibin)/sqrt(v) oplot,replicate(binx(ibin)+0.5*binsize,2),binmean(ibin)+halfrange*[-1,1],$ thick=2 endif endfor ; end