; ; Creates regional timeseries based on a limited sample of MXD series, ; variance adjusted with a time-dependent rbar, and outputs them with ; regional instrumental records which are also computed here and variance ; adjusted. The instrumental records are for each month and from ; temperature and precipitation. ; ; 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). ; ; Output timeseries are for Ed Cook to do some Kalman filtering with. ; ;----------------------------------------------------------------- ; ; Specify options ; corrper=[1881,1990] ; 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 ; ; 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 monname=['J','F','M','A','M','J','J','A','S','O','N','D'] nmon=12 ; ; First read in the MXD data (and truncate as necessary) ; print,'Reading MXD data' restore,filename='allmxd.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 MXD local correlations' restore,filename='mxd_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) ; ; Define arrays ; mxdts=fltarr(nreg,nyr)*!values.f_nan temts=fltarr(nreg,nmon,nyr)*!values.f_nan prets=fltarr(nreg,nmon,nyr)*!values.f_nan ; ; 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 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 temts(ireg,*,*)=latmon(*,*) ; ; 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 prets(ireg,*,*)=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) regfra=fraction(*,tl) regcor=reform(locr(tl,18,1)) ; 18=Apr-Sep 1=temperature ; ; Next 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 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 endif else begin meanmxd=fltarr(nyr)*!values.f_nan endelse mxdts(ireg,*)=meanmxd(*) ; ; Repeat for next region ; endfor print ; ; Output all the data ; missval=-99.9 ml=where(finite(mxdts) eq 0,nmiss) if nmiss gt 0 then mxdts(ml)=missval ml=where(finite(temts) eq 0,nmiss) if nmiss gt 0 then temts(ml)=missval ml=where(finite(prets) eq 0,nmiss) if nmiss gt 0 then prets(ml)=missval for ireg = 0 , nreg-1 do begin ; fn='reg'+regname(ireg)+'_mxd.dat' print,fn openw,1,fn printf,1,'Year MXD' for iyr = corrper(0) , corrper(1) do begin printf,1,iyr,mxdts(ireg,iyr-analper(0)),format='(I4,F8.2)' endfor close,1 ; fn='reg'+regname(ireg)+'_tem.dat' print,fn openw,1,fn printf,1,'Year Temp J F M A M J J A S O N D' for iyr = corrper(0) , corrper(1) do begin printf,1,iyr,temts(ireg,*,iyr-analper(0)),format='(I4,12F8.2)' endfor close,1 ; fn='reg'+regname(ireg)+'_pre.dat' print,fn openw,1,fn printf,1,'Year Prec J F M A M J J A S O N D' for iyr = corrper(0) , corrper(1) do begin printf,1,iyr,prets(ireg,*,iyr-analper(0)),format='(I4,12F8.2)' endfor close,1 ; endfor ; end