PRO wetmodel,variables ; IDL program to plot figures for Lost Creek publication ; ;variables: ; no entry = reload large variable datasets ; any number/character(s); E.g., '1' or 'y' = reload small variable ; dataset print,'Loading data...' IF n_elements(variables) EQ 0 THEN begin print,'Restoring original variables...' ;need to restore variables in the following files: restore,'/eddy/s2/bcook/variables/par.sav' restore,'/eddy/s2/bcook/variables/filled_nee.sav' restore,'/eddy/s2/bcook/variables/soilt_h2o.sav' ; restore,'/eddy/s2/bcook/variables/afm.sav' restore,'/eddy/s2/bcook/variables/daily.sav' restore,'/eddy/s2/bcook/variables/met_obs.sav' restore,'/eddy/s2/bcook/variables/h2o_vpd_wind.sav' ;open phenology files (for triming data to estimated leaf on/off dates) loadsmall,'/eddy/s1/cheas/phenology/wcreek_leafon_dates',wc_leafon loadsmall,'/eddy/s1/cheas/phenology/aspen_leafon_dates',asp_leafon ;save variables save,filename='/eddy/s2/bcook/variables/wet_mss.sav',par_in_fill,wcreek_filled_nee,soilt_5cm_wc_fill,wc_leafon,asp_leafon,daily_met,mcanopy_wc_fill,max_par,solarzenith,wcreek_met,wc_vpd_top_fill ENDIF ELSE BEGIN print,'Restoring wet_mss variables...' restore,'/eddy/s2/bcook/variables/wet_mss.sav' ENDELSE stop ;hourly k computed from solar zenith (assuming spherical orientation, x=1) zenithk=0.5/cos(solarzenith*!pi/180) ;solar zenith with nighttime removed sz_nonight=solarzenith sz_nonight(where(solarzenith ge 90))=0. ;VARIABLES ht=[0,2,7.6,12.2,18.3,24]; belowcanopy sensor heights ;name daily variables year=daily_met(0,*); year month=daily_met(1,*); month jday=daily_met(2,*); jday obs24=daily_met(28,*); observed incoming par obs18=daily_met(31,*); observed 18.3m obs12=daily_met(32,*);observed 12.2m obs7=daily_met(33,*);observed 7.6m obs2=daily_met(34,*);observed 2m obs_asp=daily_met(29,*); observed aspen, 2.6m ;yearly indexes i2000=where(year EQ 2000); entire year i2001=where(year EQ 2001) i2002=where(year EQ 2002) i2003=where(year EQ 2003) i2004=where(year EQ 2004) i00_02=where(year GE 2000 AND year LE 2002) ;CALCULATIONS ;daily GPP, ER (use modeled ER and gap-filled NEE; g/m2,d), deltaG++, incoming PAR (mol/m2,d), and filled soil T gpp=fltarr(n_elements(wcreek_filled_nee(0,*))/48.); er=gpp; modeled er nee=gpp; gap-filled nee par=gpp; filled PAR mpar=gpp; maximum PAR szk=gpp; solar zenith k sz=gpp; solar zenith angle sz10=gpp; solar zenith angle at 10 am LST mvpd=gpp; maximum VPD avpd=gpp; average daylight VPD soil_t=gpp; 5cm soil T can_t=gpp; mid-canopy T min_t=gpp; minimum air T neevals=gpp FOR i=0l,n_elements(gpp)-1 DO BEGIN IF year(i) GE 2000 THEN begin rows=where(wc_leafon(0,*) EQ year(i)) leafout=wc_leafon(1,rows) leafoff=wc_leafon(6,rows) ; IF (jday(i) GE leafout AND jday(i) LE leafoff) THEN $; (modeled er)-(observed nee) ; gpp(i)=total([wcreek_filled_nee(8,i*48.:i*48.+47.)-wcreek_filled_nee(16,i*48.:i*48.+47.)])*1e-6*12.011*(60.*30.) gpp(i)=total([wcreek_filled_nee(13,i*48.:i*48.+47.)])*1e-6*12.011*(60.*30.); modeled gpp er(i)=total([wcreek_filled_nee(8,i*48.:i*48.+47.)])*1e-6*12.011*(60.*30.) nee(i)=total([wcreek_filled_nee(16,i*48:i*48+47.)])*1e-6*12.011*(60.*30.) neevals(i)=n_elements(where(wcreek_filled_nee(14,i*48:i*48+47) GT -900)); number of nee observations ENDIF par(i)=total([par_in_fill(i*48.:i*48.+47.)])*(60.*30.)/1000000. mpar(i)=total([max_par(i*48.:i*48.+47.)])*(60.*30.)/1000000.; maximum incoming par mvpd(i)=max([wc_vpd_top_fill(i*48.:i*48.+47.)]) avpd(i)=mean([wc_vpd_top_fill(i*48+where(par_in_fill(i*48.:i*48.+47.) GT 0))]) soil_t(i)=mean([soilt_5cm_wc_fill(i*48:i*48+47)])+273.15; K can_t(i)=mean([mcanopy_wc_fill(i*48:i*48+47)])+273.15 min_t(i)=min([mcanopy_wc_fill(i*48:i*48+47)])+273.15 ; szk(i)=total(zenithk(i*48:i*48+47))/n_elements(where(zenithk(i*48:i*48+47) GT 0)) ; sz(i)= total(sz_nonight(i*48:i*48+47))/n_elements(where(sz_nonight(i*48:i*48+47) GT 0)) ; weighted averages (using incident PAR) ; szk(i)=total(zenithk(i*48:i*48+47)*par_in_fill(i*48:i*48+47))/total(par_in_fill(i*48:i*48+47)) sz(i)= total(solarzenith(i*48:i*48+47)*par_in_fill(i*48:i*48+47))/total(par_in_fill(i*48:i*48+47)) ; sz10(i)=sz_nonight(i*48+20) ENDFOR ;fill day 366 with 365 index=where(jday EQ 366) gpp(index)=gpp(index-1) er(index)=er(index-1) nee(index)=nee(index-1) neevals(index)=neevals(index-1) ;yearly indexes (>50% real observations) o2000=where(year EQ 2000 AND neevals GE 24) o2001=where(year EQ 2001 AND neevals GE 24) o2002=where(year EQ 2002 AND neevals GE 24) o200x=where(year GE 2000 AND year LE 2002 AND neevals GE 24) gpp(where(gpp Le 0))=0. gpp(where(year LT 2000))=-999. ;cloudiness coefficient ccoef=1-(par/mpar) ;Transmitted PAR (Iz/Io) t18=obs18/obs24 index=where((obs18 GT obs24) OR obs24 EQ -999. OR obs18 EQ -999.) IF index(0) NE -1 THEN t18(index)=-999. t12=obs12/obs24 index=where((obs12 GT obs24) OR obs24 EQ -999. OR obs12 EQ -999.) IF index(0) NE -1 THEN t12(index)=-999. t7=obs7/obs24 index=where((obs7 GT obs24) OR obs24 EQ -999. OR obs7 EQ -999.) IF index(0) NE -1 THEN t7(index)=-999. t2=obs2/obs24 index=where((obs2 GT obs24) OR obs24 EQ -999. OR obs2 EQ -999.) IF index(0) NE -1 THEN t2(index)=-999. tasp=obs_asp/daily_met(28,*) index=where((obs_asp GT daily_met(28,*)) OR daily_met(28,*) EQ -999. OR obs_asp EQ -999.) IF index(0) NE -1 THEN tasp(index)=-999. ;Background/Leafoff (Iz/Io and Wood Area Index, WAI) - 2 weeks before leaf on, 2 weeks after leaf off ;2000 data row=where(wc_leafon(0,*) EQ 2000.); WCreek leafon00=wc_leafon(1,row) leafoff00=wc_leafon(6,row) index=where(year EQ 2000 AND ((jday GE leafon00(0)-14 and jday lt leafon00(0)) OR (jday gt leafoff00(0) AND jday LE leafoff00(0)+14)) AND t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) b18_00=mean([t18(index)]) b12_00=mean([t12(index)]) b7_00=mean([t7(index)]) b2_00=mean([t2(index)]) row=where(asp_leafon(0,*) EQ 2000.); Aspen leafon_asp00=asp_leafon(1,row) leafoff_asp00=asp_leafon(6,row) index=where(year EQ 2000 AND ((jday GE leafon_asp00(0)-14 and jday lt leafon_asp00(0)) OR (jday gt leafoff_asp00(0) AND jday LE leafoff_asp00(0)+14)) AND tasp NE -999.) basp_00=mean([tasp(index)]) ;2001 data row=where(wc_leafon(0,*) EQ 2001.); WCreek leafon01=wc_leafon(1,row) leafoff01=wc_leafon(6,row) index=where(year EQ 2001 AND ((jday GE leafon01(0)-14 and jday lt leafon01(0)) OR (jday GT leafoff01(0) AND jday LE leafoff01(0)+14)) AND t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) b18_01=mean([t18(index)]) b12_01=mean([t12(index)]) b7_01=mean([t7(index)]) b2_01=mean([t2(index)]) row=where(asp_leafon(0,*) EQ 2001.); Aspen leafon_asp01=asp_leafon(1,row) leafoff_asp01=asp_leafon(6,row) index=where(year EQ 2001 AND ((jday GE leafon_asp01(0)-14 and jday lt leafon_asp01(0)) OR (jday gt leafoff_asp01(0) AND jday LE leafoff_asp01(0)+14)) AND tasp NE -999.) basp_01=mean([tasp(index)]) ;2002 data row=where(wc_leafon(0,*) EQ 2002.); WCreek leafon02=wc_leafon(1,row) leafoff02=wc_leafon(6,row) index=where(year EQ 2002 AND ((jday GE leafon02(0)-14 and jday lt leafon02(0)) OR (jday gt leafoff02(0) AND jday LE leafoff02(0)+14)) AND t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) b18_02=mean([t18(index)]) b12_02=mean([t12(index)]) b7_02=mean([t7(index)]) b2_02=mean([t2(index)]) row=where(asp_leafon(0,*) EQ 2002.); Aspen leafon_asp02=asp_leafon(1,row) leafoff_asp02=asp_leafon(6,row) index=where(year EQ 2002 AND ((jday GE leafon_asp02(0)-14 and jday lt leafon_asp02(0)) OR (jday gt leafoff_asp02(0) AND jday LE leafoff_asp02(0)+14)) AND tasp NE -999.) basp_02=mean([tasp(index)]) ;2003 data row=where(wc_leafon(0,*) EQ 2003.); WCreek leafon03=wc_leafon(1,row) leafoff03=wc_leafon(6,row) index=where(year EQ 2003 AND ((jday GE leafon03(0)-14 and jday lt leafon03(0)) OR (jday gt leafoff03(0) AND jday LE leafoff03(0)+14)) AND t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) b18_03=mean([t18(index)]) b12_03=mean([t12(index)]) b7_03=mean([t7(index)]) b2_03=mean([t2(index)]) row=where(asp_leafon(0,*) EQ 2003.); Aspen leafon_asp03=asp_leafon(1,row) leafoff_asp03=asp_leafon(6,row) index=where(year EQ 2003 AND ((jday GE leafon_asp03(0)-14 and jday lt leafon_asp03(0)) OR (jday gt leafoff_asp03(0) AND jday LE leafoff_asp03(0)+14)) AND tasp NE -999.) basp_03=mean([tasp(index)]) ;2004 data row=where(wc_leafon(0,*) EQ 2004.); WCreek leafon04=wc_leafon(1,row) leafoff04=wc_leafon(6,row) index=where(year EQ 2004 AND ((jday GE leafon04(0)-14 and jday lt leafon04(0)) OR (jday gt leafoff04(0) AND jday LE leafoff04(0)+14)) AND t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) b18_04=mean([t18(index)]) b12_04=mean([t12(index)]) b7_04=mean([t7(index)]) b2_04=mean([t2(index)]) row=where(asp_leafon(0,*) EQ 2004.); Aspen leafon_asp04=asp_leafon(1,row) leafoff_asp04=asp_leafon(6,row) index=where(year EQ 2004 AND ((jday GE leafon_asp04(0)-14 and jday lt leafon_asp04(0)) OR (jday gt leafoff_asp04(0) AND jday LE leafoff_asp04(0)+14)) AND tasp NE -999.) basp_04=mean([tasp(index)]) ;PAR extinction coefficent (k) ; - use July/Aug mean as peak leaf-on period ; - good LAI and (Iz/Io) only available for 2000 and 2002 ;2000 data index=where(year EQ 2000 AND (month GE 7 AND month LE 8) AND obs24 NE -999. AND obs2 NE -999.); WCreek p2_00=mean([t2(index)]) index=where(year EQ 2000 AND (month GE 7 AND month LE 8) AND daily_met(28,*) NE -999. AND tasp NE -999.); Aspen pasp_00=mean([tasp(index)]) ;2002 data index=where(year EQ 2002 AND (month GE 7 AND month LE 8) AND obs24 NE -999. AND obs2 NE -999.); WCreek p2_02=mean([t2(index)]) index=where(year EQ 2002 AND (month GE 7 AND month LE 8) AND daily_met(28,*) NE -999. AND tasp NE -999.); Aspen pasp_02=mean([tasp(index)]) ;k kfit=linfit([0,0,5.26,6.00],[alog([b2_00,b2_02,p2_00,p2_02])]) print,'WILLOW CREEK k = '+string(kfit(1)) kfit_00=linfit([0,5.26],[alog([b2_00,p2_00])]) print,' 2000 = '+string(kfit_00(1)) kfit_02=linfit([0,6.00],[alog([b2_02,p2_02])]) print,' 2002 = '+string(kfit_02(1)) print,' ' kfit_asp=linfit([0,0,4.30,4.73],[alog([basp_00,basp_02,pasp_00,pasp_02])]) print,'ASPEN k = '+string(kfit_asp(1)) kfit_asp00=linfit([0,4.30],[alog([basp_00,pasp_00])]) print,' 2000 = '+string(kfit_asp00(1)) kfit_asp02=linfit([0,4.73],[alog([basp_02,pasp_02])]) print,' 2002 = '+string(kfit_asp02(1)) ;Defoliation period LAI (jday 150-180) mean (for WC profile plot) ;2000 data index=where(year EQ 2000 AND (jday GE 150 and jday le 180) AND $ t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) d18_00=mean([t18(index)]) d12_00=mean([t12(index)]) d7_00=mean([t7(index)]) d2_00=mean([t2(index)]) ;2001 data index=where(year EQ 2001 AND (jday GE 150 and jday le 180) AND $ t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) d18_01=mean([t18(index)]) d12_01=mean([t12(index)]) d7_01=mean([t7(index)]) d2_01=mean([t2(index)]) ;2002 data index=where(year EQ 2002 AND (jday GE 150 and jday le 180) AND $ t18 NE -999. AND t12 NE -999. AND t7 NE -999. AND t2 NE -999.) d18_02=mean([t18(index)]) d12_02=mean([t12(index)]) d7_02=mean([t7(index)]) d2_02=mean([t2(index)]) ;LAI calculation (method #1 - contrained by litter fall) ; - used avg k for 2001, but annual leaf-off intercept) lai=alog(t2); WCreek lai(where(year LT 2000))=-999. lai(i2000)=(lai(i2000)-kfit_00(0))/kfit_00(1) lai(i2001)=(lai(i2001)-kfit(0))/kfit(1) lai(i2002)=(lai(i2002)-kfit_02(0))/kfit_02(1) lai(i2003)=(lai(i2003)-kfit(0))/kfit(1) lai(i2004)=(lai(i2004)-kfit(0))/kfit(1) lai(where(t2 EQ -999. OR year LT 2000))=-999. lai(where(lai LT 0.1 and lai ne -999.))=-999.; not sensitive below 0.1 LAI (this screen cuts out noise) leafage=jday; array containing leaf age leafage(where(year LT 2000))=0. lai_asp=alog(tasp); Aspen lai_asp(where(year LT 2000))=-999. lai_asp(i2000)=(lai_asp(i2000)-kfit_asp00(0))/kfit_asp00(1) lai_asp(i2001)=(lai_asp(i2001)-kfit_asp(0))/kfit_asp(1) lai_asp(i2002)=(lai_asp(i2002)-kfit_asp02(0))/kfit_asp02(1) lai_asp(i2003)=(lai_asp(i2003)-kfit_asp(0))/kfit_asp(1) lai_asp(i2004)=(lai_asp(i2004)-kfit_asp(0))/kfit_asp(1) lai_asp(where(tasp EQ -999. OR year LT 2000))=-999. lai_asp(where(lai_asp LT 0.1 and lai_asp ne -999.))=-999.; not sensitive below 0.1 LAI (this screen cuts out noise) leafage_asp=jday; array containing leaf age leafage_asp(where(year LT 2000))=0. ;LAI (WC and aspen), and IPAR and fPAR (WC only) time series ; - gpp used to signal leaf out, par used to signal leaf off ; - interpolation follow by smoothing ;2000 lai00=lai(i2000); WCreek fipar00=1.0-t2(i2000) gppstart00=min(jday(where(gpp(i2000) GT 0.2))) noleaves=where(jday(i2000) LT gppstart00(0) OR jday(i2000) GT leafoff00(0)) leafage(i2000)=leafage(i2000)-gppstart00+1 leafage(i2000(noleaves))=0. lai00(noleaves)=0.; set leafoff periods to LAI=0 fipar00(noleaves)=1.0-b2_00 lai00(gppstart00-1)=0.1; mark leaf on day lai00(leafoff00(0)-1)=0.1; mark leaf off day goodvals=where(lai00 ne -999. AND fipar00 LE 1.0) lai00_int=interpol(lai00(goodvals),jday(i2000(goodvals)),jday(i2000)) lai00_smooth=lai00_int lai00_smooth=smooth(lai00_smooth,7) fapar00_smooth=1.0-(exp(lai00_smooth*kfit_00(1))) fipar00_int=interpol(fipar00(goodvals),jday(i2000(goodvals)),jday(i2000)) fipar00_smooth=fipar00_int fipar00_smooth=smooth(fipar00_smooth,7) lai_asp00=lai_asp(i2000); Aspen noleaves=where(jday(i2000) LT gppstart00(0) OR jday(i2000) GT leafoff_asp00(0)) leafage_asp(i2000)=leafage(i2000)-gppstart00+1 leafage_asp(i2000(noleaves))=0. lai_asp00(noleaves)=0.; set leafoff periods to LAI=0 lai_asp00(gppstart00-1)=0.1; mark leaf on day lai_asp00(leafoff_asp00(0)-1)=0.1; mark leaf off day goodvals=where(lai_asp00 NE -999.) lai_asp00_int=interpol(lai_asp00(goodvals),jday(i2000(goodvals)),jday(i2000)) lai_asp00_smooth=lai_asp00_int lai_asp00_smooth=smooth(lai_asp00_smooth,7) ;2001 lai01=lai(i2001); WCreek fipar01=1.0-t2(i2001) gppstart01=min(jday(where(gpp(i2001) GT 0.2))) noleaves=where(jday(i2001) LT gppstart01(0) OR jday(i2001) GT leafoff01(0)) leafage(i2001)=leafage(i2001)-gppstart01+1 leafage(i2001(noleaves))=0. lai01(noleaves)=0.; set leafoff periods to LAI=0 fipar01(noleaves)=1.0-b2_01 lai01(gppstart01-1)=0.1; mark leaf on day lai01(leafoff01(0)-1)=0.1; mark leaf off day goodvals=where(lai01 NE -999. AND fipar01 LE 1.0) lai01_int=interpol(lai01(goodvals),jday(i2001(goodvals)),jday(i2001)) lai01_smooth=lai01_int lai01_smooth=smooth(lai01_smooth,7) fapar01_smooth=1.0-(exp(lai01_smooth*kfit(1))) fipar01_int=interpol(fipar01(goodvals),jday(i2001(goodvals)),jday(i2001)) fipar01_smooth=fipar01_int fipar01_smooth=smooth(fipar01_smooth,7) lai_asp01=lai_asp(i2001); Aspen noleaves=where(jday(i2001) LT gppstart01(0) OR jday(i2001) GT leafoff_asp01(0)) leafage_asp(i2001)=leafage_asp(i2001)-gppstart01+1 leafage_asp(i2001(noleaves))=0. lai_asp01(noleaves)=0.; set leafoff periods to LAI=0 lai_asp01(gppstart01-1)=0.1; mark leaf on day lai_asp01(leafoff_asp01(0)-1)=0.1; mark leaf off day goodvals=where(lai_asp01 NE -999.) lai_asp01_int=interpol(lai_asp01(goodvals),jday(i2001(goodvals)),jday(i2001)) lai_asp01_smooth=lai_asp01_int lai_asp01_smooth=smooth(lai_asp01_smooth,7) ;2002 lai02=lai(i2002); WCreek fipar02=1.0-t2(i2002) gppstart02=min(jday(where(gpp(i2002) GT 0.2))) noleaves=where(jday(i2002) LT gppstart02(0) OR jday(i2002) GT leafoff02(0)) leafage(i2002)=leafage(i2002)-gppstart02+1 leafage(i2002(noleaves))=0. lai02(noleaves)=0.; set leafoff periods to LAI=0 fipar02(noleaves)=1.0-b2_02 lai02(gppstart02-1)=0.1; mark leaf on day lai02(leafoff02(0)-1)=0.1; mark leaf off day goodvals=where(lai02 NE -999. AND fipar02 LE 1.0) lai02_int=interpol(lai02(goodvals),jday(i2002(goodvals)),jday(i2002)) lai02_smooth=lai02_int lai02_smooth=smooth(lai02_smooth,7) fapar02_smooth=1.0-(exp(lai02_smooth*kfit_02(1))) fipar02_int=interpol(fipar02(goodvals),jday(i2002(goodvals)),jday(i2002)) fipar02_smooth=fipar02_int fipar02_smooth=smooth(fipar02_smooth,7) lai_smooth=lai; concatonate smoothed arrays lai_smooth(i2000)=lai00_smooth lai_smooth(i2001)=lai01_smooth lai_smooth(i2002)=lai02_smooth lai_asp02=lai_asp(i2002); Aspen noleaves=where(jday(i2002) LT gppstart02(0) OR jday(i2002) GT leafoff_asp02(0)) leafage_asp(i2002)=leafage(i2002)-gppstart02+1 leafage_asp(i2002(noleaves))=0. lai_asp02(noleaves)=0.; set leafoff periods to LAI=0 lai_asp02(gppstart02-1)=0.1; mark leaf on day lai_asp02(leafoff_asp02(0)-1)=0.1; mark leaf off day goodvals=where(lai_asp02 NE -999.) lai_asp02_int=interpol(lai_asp02(goodvals),jday(i2002(goodvals)),jday(i2002)) lai_asp02_smooth=lai_asp02_int lai_asp02_smooth=smooth(lai_asp02_smooth,7) ;2003 lai03=lai(i2003); WCreek fipar03=1.0-t2(i2003) gppstart03=min(jday(where(gpp(i2003) GT 0.2))) noleaves=where(jday(i2003) LT gppstart03(0) OR jday(i2003) GT leafoff03(0)) leafage(i2003)=leafage(i2003)-gppstart03+1 leafage(i2003(noleaves))=0. lai03(noleaves)=0.; set leafoff periods to LAI=0 fipar03(noleaves)=1.0-b2_03 lai03(gppstart03-1)=0.1; mark leaf on day lai03(leafoff03(0)-1)=0.1; mark leaf off day goodvals=where(lai03 NE -999. AND fipar03 LE 1.0) lai03_int=interpol(lai03(goodvals),jday(i2003(goodvals)),jday(i2003)) lai03_smooth=lai03_int lai03_smooth=smooth(lai03_smooth,7) fapar03_smooth=1.0-(exp(lai03_smooth*kfit(1))) fipar03_int=interpol(fipar03(goodvals),jday(i2003(goodvals)),jday(i2003)) fipar03_smooth=fipar03_int fipar03_smooth=smooth(fipar03_smooth,7) lai_asp03=lai_asp(i2003); Aspen noleaves=where(jday(i2003) LT gppstart03(0) OR jday(i2003) GT leafoff_asp03(0)) leafage_asp(i2003)=leafage_asp(i2003)-gppstart03+1 leafage_asp(i2003(noleaves))=0. lai_asp03(noleaves)=0.; set leafoff periods to LAI=0 lai_asp03(gppstart03-1)=0.1; mark leaf on day lai_asp03(leafoff_asp03(0)-1)=0.1; mark leaf off day goodvals=where(lai_asp03 NE -999.) lai_asp03_int=interpol(lai_asp03(goodvals),jday(i2003(goodvals)),jday(i2003)) lai_asp03_smooth=lai_asp03_int lai_asp03_smooth=smooth(lai_asp03_smooth,7) ;2004 lai04=lai(i2004); WCreek fipar04=1.0-t2(i2004) gppstart04=min(jday(where(gpp(i2004) GT 0.2))) noleaves=where(jday(i2004) LT gppstart04(0) OR jday(i2004) GT leafoff04(0)) leafage(i2004)=leafage(i2004)-gppstart04+1 leafage(i2004(noleaves))=0. lai04(noleaves)=0.; set leafoff periods to LAI=0 fipar04(noleaves)=1.0-b2_04 lai04(gppstart04-1)=0.1; mark leaf on day lai04(leafoff04(0)-1)=0.1; mark leaf off day goodvals=where(lai04 NE -999. AND fipar04 LE 1.0) lai04_int=interpol(lai04(goodvals),jday(i2004(goodvals)),jday(i2004)) lai04_smooth=lai04_int lai04_smooth=smooth(lai04_smooth,7) fapar04_smooth=1.0-(exp(lai04_smooth*kfit(1))) fipar04_int=interpol(fipar04(goodvals),jday(i2004(goodvals)),jday(i2004)) fipar04_smooth=fipar04_int fipar04_smooth=smooth(fipar04_smooth,7) lai_asp04=lai_asp(i2004); Aspen noleaves=where(jday(i2004) LT gppstart04(0) OR jday(i2004) GT leafoff_asp04(0)) leafage_asp(i2004)=leafage_asp(i2004)-gppstart04+1 leafage_asp(i2004(noleaves))=0. lai_asp04(noleaves)=0.; set leafoff periods to LAI=0 lai_asp04(gppstart04-1)=0.1; mark leaf on day lai_asp04(leafoff_asp04(0)-1)=0.1; mark leaf off day goodvals=where(lai_asp04 NE -999.) lai_asp04_int=interpol(lai_asp04(goodvals),jday(i2004(goodvals)),jday(i2004)) lai_asp04_smooth=lai_asp04_int lai_asp04_smooth=smooth(lai_asp04_smooth,7) ;estimated peak reduction in LAI print,'LAI during defoliation peak (jday 175):' print,' 2000 = '+string(lai00_smooth(174)) print,' 2001 = '+string(lai01_smooth(174)) print,' 2002 = '+string(lai02_smooth(174)) print,' Fraction lost ='+ $ string(1.-(lai_smooth(i2001(174))/mean([lai_smooth(i2000(174)),lai_smooth(i2002(174))]))) ;NEED TO CORRECT TERMINOLOGY--THIS IS "LEAF ONLY" APAR (not fIPAR) SINCE CALC FROM LAI W/O INTERCEPT ;(Iz/Io) time series ; - calculated from gap-filled LAI: (Iz/Io)=exp(LAI*k) izio=lai_smooth; concatonated years index=where(izio NE -999.) izio(index)=exp(izio(index)*kfit(1)) ;fIPAR = 1.-(Iz/Io) fipar=izio fipar(index)=1.-fipar(index) ;IPAR (MJ/m2,d) = (incoming par,mol/m2,d)*fIPAR*(0.218 MJ/mol photons) ipar=fipar ipar(index)=ipar(index)*par(index)*0.218 ;eg = GPP/IPAR eg=gpp eg(index)=gpp(index)/ipar(index) ;PLOTS ;Define colors: ; 0=black ; 1=red ; 2=green ; 3=blue ; 4=cyan ; 5=magenta ; 6=yellow ; 7=white red=[0,1,0,0,0,1,1,1]*255 green=[0,0,1,0,1,0,1,1]*255 blue=[0,0,0,1,1,1,0,1]*255 tvlct,red,green,blue ;create open circle for ploting a=findgen(17)*(!pi*2/16.) usersym,cos(a),sin(a) ;Plot LAI vs ht (2000/02 vs 2001) during defoliation b_00=[b2_00,b7_00,b12_02,b18_02]; 2000 d_00=[d2_00,d7_00,d12_02,d18_02] dpoints_00=(alog(d_00)-alog(b_00))/kfit(1) dpoints_00=[dpoints_00(0),dpoints_00,0.] dprofile_00=spline(ht,dpoints_00,findgen(25)) b_01=[b2_01,b7_01,b12_01,b18_01]; 2001 d_01=[d2_01,d7_01,d12_01,d18_01] dpoints_01=(alog(d_01)-alog(b_01))/kfit(1); defoliated dpoints_01=[dpoints_01(0),dpoints_01,0.] dprofile_01=spline(ht,dpoints_01,findgen(25)) b_02=[b2_02,b7_02,b12_02,b18_02]; 2002 d_02=[d2_02,d7_02,d12_02,d18_02] dpoints_02=(alog(d_02)-alog(b_02))/kfit(1); defoliated dpoints_02=[dpoints_02(0),dpoints_02,0.] dprofile_02=spline(ht,dpoints_02,findgen(25)) window,1,title='LAI vs Canopy ht',xsize=400,ysize=400 plot,dpoints_00(1:4),ht(1:4),xtitle='LAI (m!u2!n m!u-2!n)',ytitle='Canopy height (m)',xrange=[0,5.5],yrange=[0,24],psym=1 oplot,dprofile_00,findgen(25) oplot,dpoints_01(1:4),ht(1:4),psym=5 oplot,dprofile_01,findgen(25),linestyle=1 oplot,dpoints_02(1:4),ht(1:4),psym=8 oplot,dprofile_02,findgen(25),color=150 legend,['2000','2001','2002'],psym=[-1,-5,-8],color=[0,0,0],linestyle=[0,1,0],/right,box=0 ;Plot LAI and GPP vs DOY window,2,title='LAI, GPP vs DOY',xsize=400,ysize=400 ;plot graphs (note:0 empty fields,1 column, 5 rows, 0 z-dim plots, starts in top-left then down) !p.multi=[0,1,2,0,1] nolabels=[' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ',' '] ;LAI plot,jday(i2000),lai_smooth(i2000),min_val=-998,xrange=[0,366],ytitle='LAI (m!u2!n m!u-2!n)',xtickname=nolabels,ymargin=[0,0],xmargin=[4,4],yrange=[0,7] oplot,jday(i2001),lai_smooth(i2001),linestyle=1 oplot,jday(i2002),lai_smooth(i2002),color=150 legend,['2000','2001','2002'],color=[0,0,150],linestyle=[0,1,0],box=0,/left ;GPP plot,jday(i2000),smooth(gpp(i2000),7,/edge_truncate),min_val=-998,xrange=[0,366],ytitle='GPP (g C m!u-2!n d!u-1!n)',xtitle='Day of year',ymargin=[0,0],xmargin=[4,4] oplot,jday(i2001),smooth(gpp(i2001),7,/edge_truncate),min_val=-998.,linestyle=1 oplot,jday(i2002),smooth(gpp(i2002),7,/edge_truncate),min_val=-998.,color=150 !p.multi=[0,1,1,0,1] !p.multi=[0,1,2,0,1] ;plot potential GPP light use efficiency (Eg) vs DOY ;Note: the following screens were used to eliminate low light and high VPD conditions: ; - max VPD lt 1.6 ; - clear skies (par/max_par > 75%) bestvals=where(gpp GT 0 AND ipar GT 0 AND ((year EQ 2002 AND lai NE -999.) OR year EQ 2000) AND mvpd LT 1.5 AND par/mpar GT 0.75) vals00=where(gpp GT 0 AND ipar GT 0 AND year eq 2000 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals01=where(gpp GT 0 AND ipar GT 0 AND year eq 2001 AND mvpd LT 1.5 AND par/mpar GT 0.75 AND lai NE -999.) vals02=where(gpp GT 0 AND ipar GT 0 AND year eq 2002 AND mvpd LT 1.5 AND par/mpar GT 0.75 AND lai NE -999.) window,3,title='Eg vs DOY',xsize=400,ysize=400 plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00),psym=1,xrange=[0,200],yrange=[0,1.4],ytitle='!7e!X!dg!n (g C m!u-2!n ground MJ!u-1!n PAR)',xtickname=nolabels,ymargin=[0,0],xmargin=[4,4] oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01),psym=-5,color=150,linestyle=1 oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02),psym=6 legend,['2000','2001','2002'],psym=[1,5,6],color=[0,150,0],/right egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals),2); quadratic fit oplot,egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) LUEmax=egfit(0)+(egfit(1)*leafage)+(egfit(2)*leafage^2) LUEmax(where(lai_smooth le 0))=0. print,'' print,'Emax from leaf age quadratic: '+string(max([egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2)])) ;plot potential GPP light use efficiency (Eg) vs DOY ***WITH LAI ADJUSTMENT*** bestvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND gpp/ipar/lai_smooth LT 0.25 AND year GE 2000 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals00=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2000 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals01=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2001 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals02=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2002 AND mvpd LT 1.5 AND par/mpar GT 0.75) ;plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='!7e!X!dg!n (g C m!u-2!n leaf MJ!u-1!n PAR))',yrange=[0,.27],ymargin=[0,0],xmargin=[4,4] ;oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals00)/lai_smooth(vals01),psym=5,color=150 ;oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02)/lai_smooth(vals02),psym=6 ;egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals)/lai_smooth(bestvals),2); quadratic fit ;oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) plot,jday(vals00)-gppstart00+1,gpp(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='!7e!X!dg!n (g C m!u-2!n leaf))',ymargin=[0,0],yrange=[0,2.8],xmargin=[4,4] oplot,jday(vals01)-gppstart01+1,gpp(vals01)/lai_smooth(vals01),psym=5,color=150 oplot,jday(vals02)-gppstart02+1,gpp(vals02)/lai_smooth(vals02),psym=6 egfit=poly_fit(leafage(bestvals),gpp(bestvals)/lai_smooth(bestvals),2); quadratic fit oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) !p.multi=[0,1,1,0,1] !p.multi=[0,1,2,0,1] ;plot potential GPP light use efficiency (Eg) vs TMIN window,12,title='Epsilon vs. TMIN',xsize=400,ysize=400 tmin_vals00=where(gpp GT 0 AND ipar GT 0 AND year eq 2000 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75) tmin_vals01=where(gpp GT 0 AND ipar GT 0 AND year eq 2001 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75) tmin_vals02=where(gpp GT 0 AND ipar GT 0 AND year eq 2002 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75) regvals=where(gpp GT 0 AND ipar GT 0 AND (year eq 2000 OR year EQ 2002) AND par/mpar GT 0.75 AND (min_t GE (0)+273.15 AND min_t LE 7.94+273.15)) reg_fit=linear_fit(min_t(regvals)-273.15,gpp(regvals)/ipar(regvals)); linear fit data between TMINmin and TMINmax (BPLUT) tmin_min=-(reg_fit(0))/reg_fit(1) maxvals=where(gpp GT 0 AND ipar GT 0 AND year GE 2000 AND avpd LT 0.65 AND par/mpar GT 0.75 AND min_t GT 7.94+273.15) emax=mean(gpp(maxvals)/ipar(maxvals)); mean LUE above TMINmax from BPLUT tmin_max=(emax-reg_fit(0))/reg_fit(1) plot,min_t(tmin_vals00)-273.15,gpp(tmin_vals00)/ipar(tmin_vals00),psym=1,ytitle='!7e!X!dg!n (g C MJ!u-1!n PAR))',xrange=[-15,25],yrange=[-0.1,1.5],xtickname=nolabels,ymargin=[0,0],xmargin=[4,4] oplot,min_t(tmin_vals01)-273.15,gpp(tmin_vals01)/ipar(tmin_vals01),psym=5,color=150 oplot,min_t(tmin_vals02)-273.15,gpp(tmin_vals02)/ipar(tmin_vals02),psym=6 oplot,[-15,tmin_min],[0,0] oplot,[tmin_min,tmin_max],[0,emax] oplot,[tmin_max,25],[emax,emax] oplot,[-15,-8,7.94,25],[0,0,1.044,1.044],linestyle=2 ;plot potential GPP light use efficiency (Eg) vs TMIN ***WITH LAI ADJUSTMENT*** tmin_vals00=where(gpp Gt 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2000 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) tmin_vals01=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2001 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) tmin_vals02=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2002 AND (avpd LT 0.65 OR min_t LT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) ;regvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND (year eq 2000 OR year EQ 2002) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25 AND (min_t GE (0)+273.15 AND min_t LE 7.94+273.15)) regvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND (year ge 2000 and year le 2002) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25 AND (min_t GE (0)+273.15 AND min_t LE 7.94+273.15)) reg_fit=linear_fit(min_t(regvals)-273.15,gpp(regvals)/ipar(regvals)/lai_smooth(regvals)); linear fit data between TMINmin and TMINmax (BPLUT) tmin_min=-(reg_fit(0))/reg_fit(1) maxvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year GE 2000 AND avpd LT 0.65 AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth Le 0.25 AND min_t GT 7.94+273.15) emax=mean(gpp(maxvals)/ipar(maxvals)/lai_smooth(maxvals)); mean LUE above TMINmax from BPLUT laimax=mean(lai_smooth(maxvals)) tmin_max=(emax-reg_fit(0))/reg_fit(1) print,' ' print,'Earea,max = '+string(emax) print,'TMINmin = '+string(tmin_min) print,'TMINmax = '+string(tmin_max) print,'T response = '+string((emax*laimax)/(tmin_max-tmin_min)) tmin_scalar=((min_t-273.15)*reg_fit(1)+reg_fit(0))/emax tmin_scalar(where(min_t LT tmin_min+273.15))=0. tmin_scalar(where(min_t GT tmin_max+273.15))=1 plot,min_t(tmin_vals00)-273.15,gpp(tmin_vals00)/ipar(tmin_vals00)/lai_smooth(tmin_vals00),psym=1,xtitle='T!dmin!u',ytitle='!7e!X!dg!n (g C m!u-2!n leaf MJ!u-1!n PAR))',xrange=[-15,25],yrange=[-0.05,0.3],ymargin=[0,0],xmargin=[4,4] oplot,min_t(tmin_vals01)-273.15,gpp(tmin_vals01)/ipar(tmin_vals01)/lai_smooth(tmin_vals01),psym=5,color=150 oplot,min_t(tmin_vals02)-273.15,gpp(tmin_vals02)/ipar(tmin_vals02)/lai_smooth(tmin_vals02),psym=6 oplot,[-15,tmin_min],[0,0] oplot,[tmin_min,tmin_max],[0,emax] oplot,[tmin_max,25],[emax,emax] oplot,[-15,-8,7.94,25],[0,0,1.044,1.044]/laimax,linestyle=2 !p.multi=[0,1,1,0,1] !p.multi=[0,1,2,0,1] ;plot potential GPP light use efficiency (Eg) vs VPD window,13,title='Epsilon vs. VPD',xsize=400,ysize=400 vpd_vals00=where(gpp GT 0 AND ipar GT 0 AND year eq 2000 AND (avpd LT 0.65 OR min_t gt 7.94+273.15) AND par/mpar GT 0.75) vpd_vals01=where(gpp GT 0 AND ipar GT 0 AND year eq 2001 AND (avpd gT 0.65 OR min_t gT 7.92+273.15) AND par/mpar GT 0.75) vpd_vals02=where(gpp GT 0 AND ipar GT 0 AND year eq 2002 AND (avpd gT 0.65 OR min_t gT 7.94+273.15) AND par/mpar GT 0.75) regvals=where(gpp GT 0 AND ipar GT 0 AND (year eq 2000 OR year EQ 2002) AND par/mpar GT 0.75 AND (avpd GE (0.65) AND avpd LE 2.5)) reg_fit=linear_fit(avpd(regvals),gpp(regvals)/ipar(regvals)); linear fit data between TMINmin and TMINmax (BPLUT) vpd_min=-(reg_fit(0))/reg_fit(1) maxvals=where(gpp GT 0 AND ipar GT 0 AND year GE 2000 AND avpd LT 0.65 AND par/mpar GT 0.75 AND min_t GT 7.94+273.15) emax=mean(gpp(maxvals)/ipar(maxvals)); mean LUE above VPDmax from BPLUT vpd_max=(emax-reg_fit(0))/reg_fit(1) plot,avpd(vpd_vals00),gpp(vpd_vals00)/ipar(vpd_vals00),psym=1,ytitle='!7e!X!dg!n (g C MJ!u-1!n PAR))',xrange=[0,3],yrange=[-0.1,1.5],xtickname=nolabels,ymargin=[0,0],xmargin=[4,4] oplot,avpd(vpd_vals01),gpp(vpd_vals01)/ipar(vpd_vals01),psym=5,color=150 oplot,avpd(vpd_vals02),gpp(vpd_vals02)/ipar(vpd_vals02),psym=6 oplot,[vpd_min,3],[0,0] oplot,[vpd_min,vpd_max],[0,emax] oplot,[0,vpd_max],[emax,emax] oplot,[0,0.65,2.5,3],[1.044,1.044,0,0],linestyle=2 ;plot potential GPP light use efficiency (Eg) vs VPD ***WITH LAI ADJUSTMENT*** vpd_vals00=where(gpp Gt 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2000 AND (avpd gT 0.65 OR min_t gT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) vpd_vals01=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2001 AND (avpd gT 0.65 OR min_t gT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) vpd_vals02=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2002 AND (avpd gT 0.65 OR min_t gT 7.94+273.15) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25) regvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND (year eq 2000 OR year EQ 2002) AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth LE 0.25 AND (avpd GE 0.65 AND avpd LE 2.5)) reg_fit=linear_fit(avpd(regvals),gpp(regvals)/ipar(regvals)/lai_smooth(regvals)); linear fit data between TMINmin and TMINmax (BPLUT) vpd_min=-(reg_fit(0))/reg_fit(1) maxvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year GE 2000 AND avpd LT 0.65 AND par/mpar GT 0.75 AND gpp/ipar/lai_smooth Le 0.25 AND min_t GT 7.94+273.15) emax=mean(gpp(maxvals)/ipar(maxvals)/lai_smooth(maxvals)); mean LUE above VPDmax from BPLUT laimax=mean(lai_smooth(maxvals)) vpd_max=(emax-reg_fit(0))/reg_fit(1) print,'VPDmin = '+string(vpd_min) print,'VPDmax = '+string(vpd_max) print,'VPD response = '+string((emax*laimax)/(vpd_max-vpd_min)) print,' ' vpd_scalar=(avpd*reg_fit(1)+reg_fit(0))/emax ;vpd_scalar(where(avpd gt vpd_min))=0.; note: there were no values gt vpd_min! vpd_scalar(where(avpd GT vpd_max))=1 plot,avpd(vpd_vals00),gpp(vpd_vals00)/ipar(vpd_vals00)/lai_smooth(vpd_vals00),psym=1,xtitle='VPD (kPa)',ytitle='!7e!X!dg!n (g C m!u-2!n leaf MJ!u-1!n PAR))',xrange=[0,3],yrange=[-0.05,0.3],ymargin=[0,0],xmargin=[4,4] oplot,avpd(vpd_vals01),gpp(vpd_vals01)/ipar(vpd_vals01)/lai_smooth(vpd_vals01),psym=5,color=150 oplot,avpd(vpd_vals02),gpp(vpd_vals02)/ipar(vpd_vals02)/lai_smooth(vpd_vals02),psym=6 oplot,[vpd_min,3],[0,0] oplot,[vpd_min,vpd_max],[0,emax] oplot,[0,vpd_max],[emax,emax] oplot,[0,0.65,2.5,3],[1.044,1.044,0,0]/laimax,linestyle=2 !p.multi=[0,1,1,0,1] ;cloudiness coefficient relationship window,5,title='Eg vs Cloudiness',xsize=400,ysize=400 ;gsindex=where(month GE 7 AND month LE 8 AND gpp/ipar/lai_smooth GT 0.1) gsindex=where(month GE 6 AND month LE 9 AND gpp/ipar/lai_smooth GT 0.1 AND year NE 2001) ;gsindex=where(month GE 7 AND month LE 8 AND gpp/ipar/lai_smooth GT 0.1 AND year NE 2001) ;e_norm=(gpp/ipar/lai_smooth)/(egfit(0)+(egfit(1)*leafage)+(egfit(2)*leafage^2)); normalized eg e_norm=(gpp/ipar/lai_smooth)/(emax*tmin_scalar*vpd_scalar); normalized epsilon plot,ccoef(gsindex),e_norm(gsindex),psym=1,xrange=[0,1],xtitle='Cloudiness coefficient',ytitle='Relative LUE',yrange=[0.5,3.5] enorm_fit=linear_fit(alog(1-ccoef(gsindex)),e_norm(gsindex)); ln fit oplot,findgen(100)/100.,enorm_fit(0)+alog(1-findgen(100)/100.)*enorm_fit(1) enorm_fit2=linear_fit(alog(1-ccoef(gsindex)),e_norm(gsindex)-1,fita=[1,0]); ln fit, intercept=1 oplot,findgen(100)/100.,(alog(1-findgen(100)/100.)*enorm_fit2(1))+1,color=1 enorm_fit3=linear_fit(ccoef(gsindex)^2,e_norm(gsindex)-1,fita=[1,0]) oplot,findgen(100)/100.,((findgen(100)/100.)^2*enorm_fit3(1))+1,color=2; power eq, intercept=1 legend,['ln','ln w/intercept=1','power eq'],color=[0,1,2],linestyle=[0,0,0],box=0,/left ;compute gpp ;;print,'Eg parameters (quadratic fit):' ;;print,' '+string(egfit) ;e_coef=egfit(0)+(egfit(1)*leafage)+(egfit(2)*leafage^2); eg coeficient (clear skies, low vpd) ;;e_cloud=e_coef*(enorm_fit(0)+alog(1-ccoef)*enorm_fit(1)); eg corrected for cloudiness (ln fit) ;;e_cloud=e_coef*((alog(1-ccoef)*enorm_fit2(1))+1); eg corrected for cloudiness (ln fit, intercept=1) ;e_cloud=e_coef*((ccoef^2*enorm_fit3(1))+1); eg corrected for cloudiness (power eq, intercept=1) ;;index=where(e_cloud LT e_coef); use cloudiness correction only if LUE would increase ;;IF index(0) NE -1 THEN e_cloud(index)=e_coef(index) ;LUEobs=e_cloud*lai_smooth ;gpp_model=LUEobs*ipar cloud_scalar=(ccoef^2*enorm_fit3(1))+1 gpp_model=emax*tmin_scalar*vpd_scalar*cloud_scalar*lai_smooth*ipar gpp_model(where(lai_smooth EQ -999))=-999. ;plot cumulative GPP window,6,title='Cumulative GPP',xsize=400,ysize=400 plot,jday(i2000),total(gpp(i2000),/cumulative),xtitle='Day of year',ytitle='Cumulative GPP (g C m!u-2!n)',xrange=[0,366],yrange=[0,1500] oplot,jday(i2000),total(gpp_model(i2000),/cumulative),linestyle=1 oplot,jday(i2001),total(gpp(i2001),/cumulative),color=1 oplot,jday(i2001),total(gpp_model(i2001),/cumulative),linestyle=1,color=1 oplot,jday(i2002),total(gpp(i2002),/cumulative),color=150 oplot,jday(i2002),total(gpp_model(i2002),/cumulative),linestyle=1,color=150 legend,['2000','2001','2002','previous model','this model'],color=[0,1,150,0,0],linestyle=[0,0,0,0,1],/left,box=0 ;create running total for ecosystem C ; - start on 1 Jan 2000 with Jon's biometry estimate ; - includes biomass, litter layer, and leaf mass soil C) ; - add cumulative nee totalc=nee totalc(where(year LT 2000))=0. totalc=18228.5-total(totalc,/cumulative); g C m-2 ;eliminated biomass_c (i.e., tree C) b/c no way to estimate soil storage ;biomass_c=totalc-433.7-10262.2; totalc minus litter AND soilc, and w/o canopy c (added below to create stemleaf_c variable) ;estimate canopy mass (total, maple, and susceptible pools) maple_sla=267.439; cm^2/g from 2000/02 litter baskets maple_lai=lai_smooth maple_lai(i2000)=maple_lai(i2000)*.7238; fraction of lai maple (from litter baskets) maple_lai(i2002)=maple_lai(i2002)*.6791 maple_lai(i2001)=(maple_lai(i2000(0:364))+maple_lai(i2002))/2.; use mean 2000/02 for 2001 lai index=where(maple_lai GT lai_smooth); delete where 2000/02 averaged values > observed LAI for 2001 IF index(0) NE -1 THEN maple_lai(index)=lai_smooth(index) maple_mass=maple_lai/maple_sla*10000. sus_sla=mean([190.65,168.59]); weighted means for susceptible species (basswood, green ash, red oak) during 2000 and 2002 sus_lai=lai_smooth-maple_lai sus_mass=sus_lai/sus_sla*10000. leaf_mass=maple_mass+sus_mass; g DW per m^2 leaf_c=leaf_mass*0.5; g C per m^2 froot_c=leaf_c*1.1; MODIS BPLUT froot_c2=(1.92*(leaf_c*0.5))+130; Raiche and Nadelhoffer (1989) ;create and degrade frass pool sus_c=sus_mass(i2001)*0.5; susceptible leaf C (assume 50% C in leaf tissue) ref_pt=sus_c(149) leafc_lost=fltarr(365); gC fras_pool=fltarr(365); gC herb=fltarr(365); gC lost to FTC herbivores exit_rate=0.0214680; g/(g,d); from Luo (2003) for metabolic litter fras_co2=fltarr(365); gC/d herb_co2=fltarr(365);gC/d FOR i=150,364 DO begin IF (i GE 150 AND i LE 176) AND (sus_c(i) LT ref_pt) THEN BEGIN leafc_lost(i)=ref_pt-sus_c(i) ref_pt=sus_c(i) ENDIF herb(i)=leafc_lost(i)*0.17; 17% of leaf consumed --> FTC biomass herb_co2(i)=leafc_lost(i)*.33; 33% of leaf consumed --> CO2 fras_co2(i)=(fras_pool(i-1)*exit_rate) fras_pool(i)=fras_pool(i-1)-fras_co2(i)+(leafc_lost(i)*.5); 50% of leaf consumed --> frass endfor ;compute leaf-off Rh+m (heterotrophic and wood maintenance respiration) temperature response ; - in the absence of plant growth and when soils < 10 deg to exclude spring root growth ; - linear fit using 2000 through 2002 data index=where((gpp le 0 OR lai_smooth EQ 0) AND soil_t LT (10+273.15) AND year Ge 2000 AND year LE 2002) rmlo_fit=linfit(soil_t(index),er(index)/totalc(index)); y=mx+b rmlo=(((soil_t)*rmlo_fit(1))+rmlo_fit(0))*(totalc-leaf_c-froot_c); note: leaf and fine root mr calculated below rmlo(where(year LT 2000 OR year GT 2002))=-999. window,7,title='Rm vs T (leaf-off)',xsize=400,ysize=400 plot,soil_t(index)-273.15,er(index)/totalc(index),xrange=[-1,10],xtitle='Soil T!d5 cm!n (!uo!nC)',ytitle='ER (g CO!d2!n-C g!u-1!n TOC d!u-1!n)',psym=1,symsize=.7 oplot,findgen(300)-273.15-1,(findgen(300)-1)*rmlo_fit(1)+rmlo_fit(0) xyouts,.2,.9,'r!dh+m, leaf-off!n = '+strtrim(string(rmlo_fit(0)),2)+' g g!u-1!n d!u-1!n',/normal xyouts,.2,.85,'k!dT!n = '+strtrim(string(rmlo_fit(1)),2)+' g g!u-1!n d!u-1!n K!u-1!n',/normal ;estimate heterotrophic respiration ;tfrass=total(fras_co2) cwd=79;coarse woody debris estimate (J Martin, unpublished data) ;Note on root detritus estimate: used Raich and Nadelhoffer eq to estimate C allocation to roots; assume 50% goes to production, 25% to Rg and 25% to Rm (Ryan, 1991) rootd00=((1.92*(218.5*0.5))+130)*0.5 rootd01=((1.92*(max(leaf_mass(i2001))*0.5))+130)*0.5 rootd02=((1.92*(252.25*0.5))+130)*0.5 rh2000=cwd+(226.4*0.5)+rootd00; 226.4=leaf mass from previous year's litter baskets (1999) ;rh2001=cwd+(218.5*0.5)+total(fras_co2)+total(herb_co2)+rootd01 rh2001=cwd+(218.5*0.5)+rootd01; add FTC CO2 separately rh2002=cwd+(max(leaf_mass(i2001))*0.5)+rootd02; max leaf mass during 2001 (computed from spp. sla) print,' ' print,'Leaf+CWD+Root Detritus=total' print,' 2000: '+string(226.4*0.5)+' + '+string(cwd)+' + '+string(rootd00)+' = '+string(rh2000) print,' 2001: '+string(218.5*0.5)+' + '+string(cwd)+' + '+string(rootd01)+' = '+string(rh2001) print,' 2002: '+string(max(leaf_mass(i2001))*0.5)+' + '+string(cwd)+' + '+string(rootd02)+' = '+string(rh2002) ;estimate daily heterotropic and woody respiration (including leaf decomposition) ;rh_rmlo=(rh2000+rh2002)/(total(rmlo(i2000))+total(rmlo(i2002))); rh as fraction of leaf-off respiration ;rh=rmlo*rh_rmlo ;print,'' ;print,'Rh as fraction of Rh+Rw during non-defoliated years: '+string(rh_rmlo*100.)+'%' ;rh(where(year LT 2000 or year GT 2002))=-999. rh=rmlo*0.402439; fraction of ER reported by Bolstad et al ;rh=rmlo/(totalc-leaf_c-froot_c) ;rh(where(year lt 2000 or year GT 2002))=-999. ;rh(i2000)=rh2000*rh(i2000)/total(rh(i2000)) ; weight annual rh by observed rates ;;rh(i2001)=((rh2001-herb_co2-fras_co2)*rh(i2001)/total(rh(i2001)))+herb_co2+fras_co2 ;rh(i2001)=rh2001*rh(i2001)/total(rh(i2001)) ;rh(i2002)=rh2002*rh(i2002)/total(rh(i2002)) wood_rm=rmlo-rh wood_rm(where(year lt 2000 or year GT 2002))=-999. ;estimate daily leaf growth ;dleafmass=leaf_mass ;dleafmass(*)=0. ;refpt=0. ;FOR i=1,n_elements(leaf_mass)-2 DO BEGIN ; IF jday(i) EQ 1 THEN refpt=0. ; IF leaf_mass(i-1) GT 0 and leaf_mass(i) GT 0 and leaf_mass(i) GT refpt THEN begin ; dleafmass(i)=leaf_mass(i)-refpt ; refpt=leaf_mass(i) ; endif ;endfor ;respriation calcs based on MODIS logic and Bolstad observations bs_wood_rmf=(3.90*12./1e9*(60.*60.*24.))*2.99^((can_t-273.15-20.)/10.); Bolstad Q10 factors gr_ash_rmf=(3.81*12./1e9*(60.*60.*24.))*2.9^((can_t-273.15-20.)/10.) rd_oak_rmf=(3.21*12./1e9*(60.*60.*24.))*2.30^((can_t-273.15-20.)/10.) sg_map_rmf=(3.89*12./1e9*(60.*60.*24.))*2.98^((can_t-273.15-20.)/10.) maple_rm=maple_mass*sg_map_rmf sus_rm=sus_mass*((0.38*bs_wood_rmf)+(0.39*gr_ash_rmf)+(0.23*rd_oak_rmf)) leaf_rm=maple_rm+sus_rm froot_rm=froot_c*0.00519*2^((soil_t-273.15-20.)/10.); fine root Q10 from MODIS BPLUT rm_total=wood_rm+leaf_rm+froot_rm; sum of maintenance respiration rm_total(where(year lt 2000 or year GT 2002))=-999. rg2000=(total(er(i2000))-total(rh(i2000))-total(rm_total(i2000)))/total(gpp(i2000)); annual growth respiration rate (g ER g-1 GPP d-1) rg2001=(total(er(i2001))-total(rh(i2001))-total(rm_total(i2001)))/total(gpp(i2001)) rg2002=(total(er(i2002))-total(rh(i2002))-total(rm_total(i2002)))/total(gpp(i2002)) rg_mean=mean([rg2000,rg2001,rg2002]) rg=gpp*rg_mean rg(where(year lt 2000 or year GT 2002))=-999. ra=rm_total+rg ra(where(year lt 2000 or year GT 2002))=-999. ;calculate NEE, NPP, CUE, and modeled ER and NEP npp=gpp_model-ra npp(where(year lt 2000 or year GT 2002))=-999. cue=npp/gpp_model cue(where(year lt 2000 or year GT 2002))=-999. er_model=rh+ra er(where(year LT 2000 OR year GT 2002))=-999. nep=npp-rh nep(where(year LT 2000 OR year GT 2002))=-999. plot,total(GPP(i2000),/cumulative),yrange=[-1000,1500] oplot,total(GPP_model(i2000),/cumulative),linestyle=1 oplot,total(npp(i2000),/cumulative),color=1 oplot,total(nep(i2000),/cumulative),color=2 oplot,-(total(nee(i2000),/cumulative)),color=2,linestyle=1 oplot,-(total(ER(i2000),/cumulative)),color=3 oplot,-(total(er_model(i2000),/cumulative)),linestyle=1,color=3 oplot,-(total(rh(i2000),/cumulative)),color=4 oplot,-(total(rm_total(i2000),/cumulative)),color=5 oplot,-(total(rg(i2000),/cumulative)),color=6 oplot,findgen(365)*0,color=150,linestyle=1 legend,['Tower GPP','Model GPP','NPP','Model NEP','Tower NEE','Tower ER','Model ER','Rh','Rm','Rg'],color=[0,0,1,2,2,3,3,4,5,6],linestyle=[0,1,0,0,1,0,1,0,0,0],/left,/top,box=0 stop ;compute leaf-on Rm and Rg, and evaluate T effect ;rmlo(i2001)=rmlo(i2001)-fras_co2-herb_co2; substract out FTC decomposition from rg + rl regression ;;stemleaf_c=biomass_c+leaf_c; add leaf C to total biomass ;window,8,title='Leaf-On Rm and Rg',xsize=400,ysize=400 ;index=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (10+273.15) AND year ge 2000 AND year LE 2002) ;;index01=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (10+273.15) AND year eq 2001) ;;index10_12=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (10+273.15) AND soil_t LE (12+273.15) AND year GE 2000) ;;index15_17=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (15+273.15) AND soil_t LE (17+273.15) AND year GE 2000) ;;index20_22=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (20+273.15) AND soil_t LE (22+273.15) AND year GE 2000) ;rl_fit=linfit((gpp(index)-rmlo(index))/stemleaf_c(index),(er(index)-rmlo(index))/stemleaf_c(index));/totalc(index)) ;;rl10_12fit=linfit((gpp(index10_12)-rmlo(index10_12))/totalc(index10_12),(er(index10_12)-rmlo(index10_12))/totalc(index10_12),sigma=sigma10_12) ;;rl15_17fit=linfit((gpp(index15_17)-rmlo(index15_17))/totalc(index15_17),(er(index15_17)-rmlo(index15_17))/totalc(index15_17),sigma=sigma15_17) ;rml=rl_fit(0)*stemleaf_c;totalc; leaf maintenance respiration ;rml(where(year LT 2000))=-999. ;rml(where(gpp-rmlo lt 0))=0. ;rg=((gpp-rmlo)/stemleaf_c)*rl_fit(1)*stemleaf_c; growth respiration ;;rg=((gpp-rmlo)/totalc)*rl_fit(1)*totalc; growth respiration ;rg(where(gpp-rmlo lt 0))=0. ;rg(where(year LT 2000))=-999. ;;plot,(gpp(index)-rmlo(index))/stemleaf_c(index),(er(index)-rmlo(index))/stemleaf_c(index),psym=1,symsize=.7,xtitle='Relative growth rate (g C g!u-1!n TOC d!u-1!n)',ytitle='ER - R!dleafless!n (g CO!d2!n-C g!u-1!n TOC d!u-1!n)',yrange=[-.00005,.0006] ;plot,gpp(index)/leaf_c(index),(er(index)-rmlo(index))/leaf_c(index),psym=1,symsize=.7,xtitle='Relative growth rate (g C g!u-1!n TOC d!u-1!n)',ytitle='ER - R!dleafless!n (g CO!d2!n-C g!u-1!n TOC d!u-1!n)',yrange=[-.00005,.0006] ;;plot,(gpp(index)-rmlo(index))/totalc(index),(er(index)-rmlo(index))/totalc(index),psym=1,symsize=.7,xtitle='Relative growth rate (g C g!u-1!n TOC d!u-1!n)',ytitle='ER - R!dleafless!n (g CO!d2!n-C g!u-1!n TOC d!u-1!n)',yrange=[-.00005,.0003] ;;oplot,(gpp(index01)-rmlo(index01))/totalc(index01),(er(index01)-rmlo(index01))/totalc(index01),psym=1,symsize=.7,color=150 ;oplot,findgen(3)-1,(findgen(3)-1)*rl_fit(1)+rl_fit(0) ;xyouts,.2,.9,'r!dm,leaf-on!n = '+strtrim(string(rl_fit(0)),2)+' g g!u-1!n d!u-1!n',/normal ;xyouts,.2,.85,'r!dg!n = '+strtrim(string(rl_fit(1)),2),/normal ;;total heterotrophic/maintenance respiration ;rhm=rmlo+rml ;rhm(where(year LT 2000))=-999. ;;compute total ER (modeled) ;er_model=rhm+rg ;er_model(where(year LT 2000))=-999. ;add herbivore and frass co2 to er er_frass=er_model er_frass(i2001)=er_frass(i2001)+fras_co2+herb_co2 ;plot cumulative ER window,9,title='Cumulative ER',xsize=400,ysize=400 plot,jday(i2000),total(er(i2000),/cumulative),xtitle='Day of year',ytitle='Cumulative ER (g C m!u-2!n)',xrange=[0,366],yrange=[0,1500] oplot,jday(i2000),total(er_model(i2000),/cumulative),linestyle=1 oplot,jday(i2001),total(er(i2001),/cumulative),color=1 ;oplot,jday(i2001),total(er_model(i2001),/cumulative),linestyle=1,color=1 oplot,jday(i2001),total(er_frass(i2001),/cumulative),linestyle=1,color=1 oplot,jday(i2002),total(er(i2002),/cumulative),color=150 oplot,jday(i2002),total(er_model(i2002),/cumulative),linestyle=1,color=150 legend,['2000','2001','2002','previous model','this model'],color=[0,1,150,0,0],linestyle=[0,0,0,0,1],/left,box=0 stop ;compute/plot modeled vs previous modeled/gapfilled NEE model_nee=er_model-gpp_model model_nee(where(year LT 2000))=-999. model_frass=er_frass-gpp_model model_frass(where(year LT 2000))=-999. window,10,title='Cumulative NEE',xsize=400,ysize=400 plot,total(model_nee(i2000),/cumulative),xrange=[0,366],yrange=[-600,200],linestyle=1,xtitle='Day of year',ytitle='Cumulative NEE (g C m!u-2!n)'; this model w/o frass oplot,total(er(i2000)-gpp(i2000),/cumulative); gap-filling model oplot,total(nee(i2000),/cumulative),thick=3; gap-filled observations oplot,total(model_nee(i2001),/cumulative),color=1,linestyle=1; oplot,total(model_frass(i2001),/cumulative),color=1,linestyle=5,thick=3; this model w/frass oplot,total(er(i2001)-gpp(i2001),/cumulative),color=1 oplot,total(nee(i2001),/cumulative),thick=3,color=1 oplot,total(model_nee(i2002),/cumulative),color=150,linestyle=1 oplot,total(er(i2002)-gpp(i2002),/cumulative),color=150 oplot,total(nee(i2002),/cumulative),thick=3,color=150 legend,['2000','2001','2001 (w/frass)','2002','previous model','gap-filled','this model'],color=[0,1,1,150,0,0,0],linestyle=[0,0,5,0,0,0,1],thick=[0,0,3,0,0,3,0],/left,/bottom,box=0 window,11,title='C allocated to growth processes',xsize=400,ysize=400 growthc=(rg-nee)/gpp*100. growthc(where(growthc LT 0))=0. growthc(where(year LT 2000))=-999. plot,jday(i2000),smooth(growthc(i2000),7),yrange=[0,100],xtitle='Day of year',ytitle='Carbon allocated to growth (%)' oplot,jday(i2001),smooth(growthc(i2001),7),linestyle=1 oplot,jday(i2002),smooth(growthc(i2002),7),nsum=7,color=150 legend,['2000','2001','2002'],color=[0,0,150],linestyle=[0,1,0],/left,box=0 print,'NEE Model error' print,' daily avg, all years (%) = '+string(mean((model_nee(o200x)-nee(o200x))/nee(o200x)*100.)) print,' 2000 (gC/m2,y) = '+string(total(model_nee(i2000))-total(nee(i2000))) print,' 2001 (gC/m2,y) = '+string(total(model_frass(i2001))+total(herb)-total(nee(i2001))) print,' 2002 (gC/m2,y) = '+string(total(model_nee(i2002))-total(nee(i2002))) ;print annual sums for various fractions ;;estimate Ra ;ra=rmlo-rh+rml+rg ;estimate daily NPP rate (note: NPP will be negative in winter!) npp=gpp_model index=where(gpp_model NE -999.) npp(index)=gpp_model(index)-ra(index) stop print,' ' print,'***** 2000 *****' print,'GPP = '+string(total(gpp_model(i2000))) print,'Rh = '+string(rh2000) print,'Rm = '+string(total(rmlo(i2000))-rh2000+total(rml(i2000))) print,' Rm,non-leaf = '+string(total(rmlo(i2000))-rh2000) print,' Rm,leaf = '+string(total(rml(i2000))) print,'Rg = '+string(total(rg(i2000))) print,'Ra = '+string(total(rmlo(i2000))-rh2000+total(rml(i2000))+total(rg(i2000))) print,'ER = '+string(total(er_model(i2000))) print,'NEE = '+string(total(er_model(i2000))-total(gpp_model(i2000))) print,'NPP = '+string(total(npp(i2000))) print,'CUE = '+string(total(npp(i2000))/total(gpp_model(i2000))); note: includes understory also print,'R:P = '+string(total(ra(i2000))/total(gpp_model(i2000))) print,' ' print,'***** 2001 *****' print,'GPP = '+string(total(gpp_model(i2001))) print,'Rh = '+string(rh2001) print,'FTC herbivory = '+string(total(herb)) print,'Rm = '+string(total(rmlo(i2001))-(rh2001-total(herb_co2+fras_co2))+total(rml(i2001))) print,' Rm,non-leaf = '+string(total(rmlo(i2001))-(rh2001-total(herb_co2+fras_co2))) print,' Rm,leaf = '+string(total(rml(i2001))) print,'Rg = '+string(total(rg(i2001))) print,'Ra = '+string(total(rmlo(i2001))-rh2001+total(rml(i2001))+total(rg(i2001))) print,'ER = '+string(total(er_frass(i2001))) print,'NEE = '+string(total(er_frass(i2001))-total(gpp_model(i2001))-total(herb)) print,'NPP = '+string(total(npp(i2001))) print,'CUE = '+string(total(npp(i2001))/total(gpp_model(i2001))) print,'R:P = '+string(total(ra(i2001))/total(gpp_model(i2001))) print,' ' print,'***** 2002 *****' print,'GPP = '+string(total(gpp_model(i2002))) print,'Rh = '+string(rh2002) print,'Rm = '+string(total(rmlo(i2002))-rh2002+total(rml(i2002))) print,' Rm,non-leaf = '+string(total(rmlo(i2002))-rh2002) print,' Rm,leaf = '+string(total(rml(i2002))) print,'Rg = '+string(total(rg(i2002))) print,'Ra = '+string(total(rmlo(i2002))-rh2002+total(rml(i2002))+total(rg(i2002))) print,'ER = '+string(total(er_model(i2002))) print,'NEE = '+string(total(er_model(i2002))-total(gpp_model(i2002))) print,'NPP = '+string(total(npp(i2002))) print,'CUE = '+string(total(npp(i2002))/total(gpp_model(i2002))) print,'R:P = '+string(total(ra(i2002))/total(gpp_model(i2002))) ;print annual ipar differences anpp=npp/1.2 print,' ' print,'Annual IPAR (MJ/year), LUEgpp, LUEnpp, LUEanpp' print,' 2000: '+string(total(ipar(i2000)))+' , '+string(total(gpp(i2000))/total(ipar(i2000)))+' , '+string(total(npp(i2000))/total(ipar(i2000)))+' , '+string(total(anpp(i2000))/total(ipar(i2000))) print,' 2001: '+string(total(ipar(i2001)))+' , '+string(total(gpp(i2001))/total(ipar(i2001)))+' , '+string(total(npp(i2001))/total(ipar(i2001)))+' , '+string(total(anpp(i2001))/total(ipar(i2001))) print,' 2002: '+string(total(ipar(i2002)))+' , '+string(total(gpp(i2002))/total(ipar(i2002)))+' , '+string(total(npp(i2002))/total(ipar(i2002)))+' , '+string(total(anpp(i2002))/total(ipar(i2002))) ;estimate daily and growing season CUE cue=npp/gpp_model cue(where(gpp_model EQ 0))=0.; note: CUE can be negative, i.e., growing but loosing mass print,' ' print,'Mean mid-growing season CUE (Jul and Aug)' gs2000=where(year EQ 2000 AND month GE 7 AND month LE 8) print,' 2000: '+string(mean([cue(gs2000)])) gs2001=where(year EQ 2001 AND month GE 7 AND month LE 8) print,' 2001: '+string(mean([cue(gs2001)])) gs2002=where(year EQ 2002 AND month GE 7 AND month LE 8) print,' 2002: '+string(mean([cue(gs2002)])) gsall=where((year ge 2000 AND year LE 2002) AND (month GE 7 AND month LE 8)) meancue=mean([cue(gsall)]) print,' 2000-03: '+string(meancue) print,' ' ;print NPP with constant CUE=0.47 (Warring et al, 1998) print,'NPP with constant CUE=0.47 (Warring et al, 1998) - gC/m2,y (% error)' print,' 2000 = '+string(total(gpp_model(i2000)*0.47))+' ('+string(abs((total(gpp_model(i2000)*0.47))-total(npp(i2000)))/total(npp(i2000))*100.)+')' print,' 2001 = '+string(total(gpp_model(i2001)*0.47))+' ('+string(abs((total(gpp_model(i2001)*0.47))-total(npp(i2002)))/total(npp(i2001))*100.)+')' print,' 2002 = '+string(total(gpp_model(i2002)*0.47))+' ('+string(abs((total(gpp_model(i2002)*0.47))-total(npp(i2002)))/total(npp(i2002))*100.)+')' ;Postscript plots for manuscript !p.charthick=3 !p.thick=3 !x.thick=3 !y.thick=3 !p.charsize=2 set_plot,'ps' ;PAR profile during defoliation period device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/lai_profile.ps',/inches,xoffset=1,yoffset=1,xsize=6.5,ysize=9,/helvetica plot,dpoints_00(1:4),ht(1:4),xtitle='LAI (m!u2!n m!u-2!n)',ytitle='Canopy height (m)',xrange=[0,5.5],yrange=[0,24],psym=1,symsize=1.5 oplot,dprofile_00,findgen(25) oplot,dpoints_01(1:4),ht(1:4),psym=5,symsize=1.5 oplot,dprofile_01,findgen(25),linestyle=1 oplot,dpoints_02(1:4),ht(1:4),psym=8,symsize=1.5,color=125 oplot,dprofile_02,findgen(25),color=125 legend,['2000','2001','2002'],psym=[-1,-5,-8],color=[0,0,125],linestyle=[0,1,0],/right,box=0,charsize=1.5 device,/close; close ps file ;LAI/GPP time series device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/lai_gpp.ps',/inches,xoffset=1.5,yoffset=1.5,xsize=6.5,ysize=9,/helvetica !p.multi=[0,1,2,0,1] plot,jday(i2000),lai_smooth(i2000),min_val=-998,xrange=[0,366],ytitle='LAI (m!u2!n m!u-2!n)',xtickname=nolabels,ymargin=[0,0],xmargin=[4,4],yrange=[0,7]; LAI oplot,jday(i2001),lai_smooth(i2001),linestyle=1 oplot,jday(i2002),lai_smooth(i2002),color=125 legend,['2000','2001','2002'],color=[0,0,125],linestyle=[0,1,0],box=0,/left,charsize=1.5 plot,jday(i2000),smooth(gpp(i2000),7,/edge_truncate),min_val=-998,xrange=[0,366],ytitle='GPP (g C m!u-2!n d!u-1!n)',xtitle='Day of year',ymargin=[0,0],xmargin=[4,4]; GPP oplot,jday(i2001),smooth(gpp(i2001),7,/edge_truncate),min_val=-998.,linestyle=1 oplot,jday(i2002),smooth(gpp(i2002),7,/edge_truncate),min_val=-998.,color=125 oplot,[0,366],[0,0]; to make x-axis black (overwrite color when y=0) !p.multi=[0,1,1,0,1] device,/close; close ps file ;LUE vs leaf age (expressed per leaf area only) device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/lue.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='LUE (g C MJ!u-1!n PAR m!u-2!n leaf)',yrange=[0,.27] oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01)/lai_smooth(vals01),psym=5 oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02)/lai_smooth(vals02),psym=8,color=125 egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals)/lai_smooth(bestvals),2); quadratic fit oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) legend,['2000','2001','2002'],psym=[1,5,8],color=[0,0,125],/right,charsize=1.5 device,/close; close ps file ;LUE vs leaf age; a) ground vs b) leaf area device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/lue2.ps',/inches,xoffset=1.5,yoffset=1.5,xsize=6.5,ysize=9,/helvetica !p.multi=[0,1,2,0,1] bestvals=where(gpp GT 0 AND ipar GT 0 AND ((year EQ 2002 AND lai NE -999.) OR year EQ 2000) AND mvpd LT 1.5 AND par/mpar GT 0.75) vals00=where(gpp GT 0 AND ipar GT 0 AND year eq 2000 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals01=where(gpp GT 0 AND ipar GT 0 AND year eq 2001 AND mvpd LT 1.5 AND par/mpar GT 0.75 AND lai NE -999.) vals02=where(gpp GT 0 AND ipar GT 0 AND year eq 2002 AND mvpd LT 1.5 AND par/mpar GT 0.75 AND lai NE -999.) plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00),psym=1,xrange=[0,200],yrange=[0,1.4],ytitle='!7e!X (g C MJ!u-1!n)',xtickname=nolabels,ymargin=[0,0],xmargin=[4,4] oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01),psym=-5,linestyle=1 oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02),psym=8,color=125 legend,['2000','2001','2002'],psym=[1,5,8],color=[0,0,125],/right,charsize=1.5 egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals),2); quadratic fit oplot,egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) ;plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00),psym=1,xrange=[0,200],yrange=[0,1.4],ytitle='!4e',xtickname=nolabels,ymargin=[0,0],xmargin=[4,4],ytickformat='(f5.2)' ;oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01),psym=-5,linestyle=1 ;oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02),psym=8,color=125 ;legend,['2000','2001','2002'],psym=[1,5,8],color=[0,0,125],/right,charsize=1.5 ;egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals),2); quadratic fit ;oplot,egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) bestvals=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND gpp/ipar/lai_smooth LT 0.25 AND year GE 2000 AND year LE 2002 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals00=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2000 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals01=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2001 AND mvpd LT 1.5 AND par/mpar GT 0.75) vals02=where(gpp GT 0 AND ipar GT 0 AND lai_smooth GT 0 AND year eq 2002 AND mvpd LT 1.5 AND par/mpar GT 0.75) ;plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='!4e!X (g C MJ!u-1!n m!u-2!n)',yrange=[0,.27],ymargin=[0,0],xmargin=[4,4] ;oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01)/lai_smooth(vals01),psym=5 ;oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02)/lai_smooth(vals02),psym=6,color=125 ;egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals)/lai_smooth(bestvals),2); quadratic fit ;oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) plot,jday(vals00)-gppstart00+1,gpp(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='GPP!darea!n (g C d!u-1!n m!u-2!n)',ymargin=[0,0],yrange=[0,2.8],xmargin=[4,4] oplot,jday(vals01)-gppstart01+1,gpp(vals01)/lai_smooth(vals01),psym=5 oplot,jday(vals02)-gppstart02+1,gpp(vals02)/lai_smooth(vals02),psym=8,color=125 egfit=poly_fit(leafage(bestvals),gpp(bestvals)/lai_smooth(bestvals),2); quadratic fit oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) ;plot,jday(vals00)-gppstart00+1,gpp(vals00)/ipar(vals00)/lai_smooth(vals00),psym=1,xrange=[0,200],xtitle='Leaf age (d)',ytitle='!4e!darea!n',yrange=[0,.27],ymargin=[0,0],xmargin=[4,4] ;oplot,jday(vals01)-gppstart01+1,gpp(vals01)/ipar(vals01)/lai_smooth(vals01),psym=5 ;oplot,jday(vals02)-gppstart02+1,gpp(vals02)/ipar(vals02)/lai_smooth(vals02),psym=8,color=125 ;egfit=poly_fit(leafage(bestvals),gpp(bestvals)/ipar(bestvals)/lai_smooth(bestvals),2); quadratic fit ;oplot,findgen(201),egfit(0)+(egfit(1)*findgen(201))+(egfit(2)*findgen(201)^2) ;xyouts,-.14,.55,'Light use efficiency (LUE)',orientation=90,alignment=0.5,/normal !p.multi=[0,1,1,0,1] device,/close; close ps file ;TMIN scalar device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/TMINscalar.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,min_t(tmin_vals00)-273.15,gpp(tmin_vals00)/ipar(tmin_vals00)/lai_smooth(tmin_vals00)/emax,psym=1,xtitle='T!dmin!u',ytitle='f!dT!n',xrange=[-15,25],yrange=[-0.1,1.3] oplot,min_t(tmin_vals01)-273.15,gpp(tmin_vals01)/ipar(tmin_vals01)/lai_smooth(tmin_vals01)/emax,psym=5 oplot,min_t(tmin_vals02)-273.15,gpp(tmin_vals02)/ipar(tmin_vals02)/lai_smooth(tmin_vals02)/emax,psym=8,color=125 oplot,[-15,tmin_min],[0,0] oplot,[tmin_min,tmin_max],[0,1] oplot,[tmin_max,25],[1,1] oplot,[-15,-8,7.94,25],[0,0,1,1],linestyle=2 legend,['2000','2001','2002','Observed fit','MODIS'],psym=[1,5,8,0,0],color=[0,0,125,0,0],linestyle=[0,0,0,0,2],box=0,/right,/bottom,charsize=1. ;need legend for lines device,/close; close ps file ;VPD scalar device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/VPDscalar.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,avpd(vpd_vals00),gpp(vpd_vals00)/ipar(vpd_vals00)/lai_smooth(vpd_vals00)/emax,psym=1,xtitle='Daylight VPD',ytitle='f!dVPD!n',xrange=[0,3],yrange=[-0.1,1.3] oplot,avpd(vpd_vals01),gpp(vpd_vals01)/ipar(vpd_vals01)/lai_smooth(vpd_vals01)/emax,psym=5 oplot,avpd(vpd_vals02),gpp(vpd_vals02)/ipar(vpd_vals02)/lai_smooth(vpd_vals02)/emax,psym=8,color=125 oplot,[vpd_min,3],[0,0] oplot,[vpd_min,vpd_max],[0,1] oplot,[0,vpd_max],[1,1] oplot,[0,0.65,2.5,3],[1,1,0,0],linestyle=2 ;legend,['2000','2001','2002'],psym=[1,5,8],color=[0,0,125],/right,charsize=1.5 legend,['2000','2001','2002','Observed fit','MODIS'],psym=[1,5,8,0,0],color=[0,0,125,0,0],linestyle=[0,0,0,0,2],box=0,/right,charsize=1. ;need legend for lines device,/close; close ps file ;Relative Eg vs cloudiness device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/cloudiness.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,ccoef(gsindex),e_norm(gsindex),psym=1,xrange=[0,1],xtitle='Cloudiness index (CI)',yrange=[0.5,2.7],ytitle='f!dCI!n';ytitle='!7e!X!darea!n/!7e!X!dcloudless!n',yrange=[0.5,2.7] ;oplot,findgen(100)/100.,enorm_fit(0)+alog(1-findgen(100)/100.)*enorm_fit(1); ln fit (variable intercept oplot,findgen(100)/100.,((findgen(100)/100.)^2*enorm_fit3(1))+1; power eq, intercept=1 ;xyouts,.3,.85,'y = '+strmid(strtrim(string(enorm_fit(0)),2),0,5)+' + ln(1 - x)',/normal,charsize=1.4 xyouts,.4,.85,'f!dCI!n = '+strmid(strtrim(string(enorm_fit3(1)),2),0,5)+'CI!u2!n',/normal,charsize=1.4 device,/close; close ps file ;Leaf-off rm, kT device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/leafoff_r.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica index=where((gpp le 0 OR lai_smooth EQ 0) AND soil_t LT (10+273.15) AND year Ge 2000) index_snow=where((gpp le 0 OR lai_smooth EQ 0) AND soil_t LT (10+273.15) AND year Ge 2000 AND daily_met(85,*) Ge 0.1) index_nosnow=where((gpp le 0 OR lai_smooth EQ 0) AND soil_t LT (10+273.15) AND year Ge 2000 AND daily_met(85,*) lt 0.1) plot,soil_t(index)-273.15,er(index)/totalc(index),xrange=[-1,10],xtitle='T!ds!n (!uo!nC)',ytitle='r!dh!n + r!dw!n (d!u-1!n)',psym=1,symsize=.7 oplot,soil_t(index_snow)-273.15,er(index_snow)/totalc(index_snow),psym=1,symsize=.7,color=125 oplot,findgen(12)-1,(findgen(12)+273.15-1)*(rmlo_fit(1))+(rmlo_fit(0)); range of fit (solid line) xyouts,.3,.87,'r!dh+w!n = '+strmid(strtrim(string(rmlo_fit(1)),2),0,4)+'x10!u-6!nT!ds!n + '+strmid(strtrim(string(rmlo_fit(0)+(273.15*rmlo_fit(1))),2),0,4)+'x10!u-4!n',/normal,charsize=1.4 ;xyouts,.3,.87,'r!dh+m!n = '+strmid(strtrim(string(1e3*rmlo_fit(0)+(273.15*1e3*rmlo_fit(1))),2),0,6)+' g kg!u-1!n d!u-1!n',/normal,charsize=1.4 ;xyouts,.3,.82,'k!dT!n = '+strmid(strtrim(1e3*string(rmlo_fit(1)),2),0,7)+' !uo!nC!u-1!n',/normal,charsize=1.4 device,/close; close ps file ;Leaf-on rm, rg device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/leafon_r.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica index=where((gpp-rmlo GT 0 and lai_smooth GT 0) AND soil_t GT (10+273.15) AND year GE 2000) plot,(gpp(index)-rmlo(index))/stemleaf_c(index),(er(index)-rmlo(index))/stemleaf_c(index),psym=1,symsize=.7,xtitle='RGR (d!u-1!n)',ytitle='r!dg!n + r!dl!n (d!u-1!n)',yrange=[-.00002,.0006],xrange=[0.0001,0.0012],ystyle=1, $ xtickname=['2x10!u-4!n',' ','6x10!u-4!n',' ','1x10!u-3!n',' ','1.4x10!u-3!n'],ytickname=['0','1x10!u-4!n','2x10!u-4!n','3x10!u-4!n','4x10!u-4!n','5x10!u-5!n','6x10!u-4!n'] oplot,(findgen(11)-1),(findgen(11)-1)*rl_fit(1)+rl_fit(0) ;xyouts,0,.55,'ER - R!dleafless!n',alignment=0.5,orientation=90,/normal ;xyouts,.6,-0.02,'(g C kg!u-1!n d!u-1!n)',alignment=0.5,/normal xyouts,.3,.87,'r!dl!n = '+strmid(strtrim(string(rl_fit(0)),2),0,4)+'x10!u-5!n d!u-1!n',/normal,charsize=1.4 xyouts,.3,.82,'r!dg!n = '+strmid(strtrim(string(rl_fit(1)),2),0,5),/normal,charsize=1.4 device,/close; close ps file ;Cumulative NEE device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/nee.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,total(model_nee(i2000),/cumulative),xrange=[0,366],yrange=[-600,200],xtitle='Day of year',ytitle='Cumulative NEE (g C m!u-2!n)' oplot,total(nee(i2000),/cumulative),thick=10 oplot,total(model_frass(i2001),/cumulative),linestyle=1 oplot,total(nee(i2001),/cumulative),linestyle=1,thick=10 oplot,total(model_nee(i2002),/cumulative),color=125 oplot,total(nee(i2002),/cumulative),color=125,thick=10 legend,['2000','2001','2002'],color=[0,0,125],linestyle=[0,1,0],/left,/bottom,box=0,charsize=1.5 device,/close; close ps file ;Rm (fraction of Ra) ;device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/rm_fraction.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica ;plot,jday(i2000),smooth((rhm(i2000)-rh(i2000))/(rhm(i2000)-rh(i2000)+rg(i2000)),7),xtitle='Day of year',ytitle='R!dM!n (fraction of R!dA!n)',yrange=[0.5,1] ;oplot,jday(i2001),smooth((rhm(i2001)-rh(i2001))/(rhm(i2001)-rh(i2001)+rg(i2001)),7),linestyle=1 ;oplot,jday(i2002),smooth((rhm(i2002)-rh(i2002))/(rhm(i2002)-rh(i2002)+rg(i2002)),7),color=125 ;oplot,[0,366],[1,1]; to make x-axis black (overwrite color when y=1) ;legend,['2000','2001','2002'],color=[0,0,125],linestyle=[0,1,0],/left,/bottom,box=0,charsize=1.5 ;device,/close; close ps file ;Forest respiration plot (to obtain 3-year average rm, rg, and max CUE) ;Notes: rg units = g glucose/g DM ; rg * (0.5 gC/g glucose) = g/gDM index=where((gpp-rh GT 0 and lai_smooth GT 0) AND soil_t GT (10+273.15) AND year GE 2000 AND year LE 2002) forestfit=linfit((gpp(index)-rh(index))/(biomass_c(index)*2.)*1e3,(er(index)-rh(index))*6./(biomass_c(index)*2)*1e3,yfit=yfit) device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/forest_r.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,(gpp(index)-rh(index))/(biomass_c(index)*2.)*1e3,(er(index)-rh(index))*6./(biomass_c(index)*2.)*1e3,psym=1,symsize=.7,xtitle='Relative growth rate',ytitle='(g glucose kg!u-1!n DM d!u-1!n)';,xrange=[0.,0.4],yrange=[0.2,1.3] oplot,(gpp(index)-rh(index))/(biomass_c(index)*2.)*1e3,yfit xyouts,.07,.55,'Specific forest respiraiton',alignment=0.5,orientation=90,/normal xyouts,.6,-0.02,'(g kg!u-1!n d!u-1!n)',alignment=0.5,/normal xyouts,.3,.87,'r!dm!n = '+strmid(strtrim(string(forestfit(0)),2),0,5)+' g kg!u-1!n DM d!u-1!n',/normal,charsize=1.4 xyouts,.3,.82,'r!dg!n = '+strmid(strtrim(string(forestfit(1)),2),0,4),/normal,charsize=1.4 device,/close; close ps file ;CUE vs DOY device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/cue.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica cue00_smooth=smooth(cue(i2000),7) plot,jday(i2000),cue00_smooth,nsum=7,xtitle='Day of year',ytitle='CUE',yrange=[0,0.8] cue01_smooth=smooth(cue(i2001),7) oplot,jday(i2001),cue01_smooth,nsum=7,linestyle=1 cue02_smooth=smooth(cue(i2002),7) oplot,jday(i2002),cue02_smooth,nsum=7,color=125 oplot,[0,366],[0,0]; to make x-axis black (overwrite color when y=0) legend,['2000','2001','2002'],color=[0,0,125],linestyle=[0,1,0],/left,box=0,charsize=1.5 device,/close; close ps file ;NPP relative errors w/constant CUE=0.47 (Warring et al., 1998) device,file='/eddy/s2/bcook/wcreek/mss_figs/caterpillar/npp_err.ps',/inches,xoffset=1,yoffset=4,xsize=6.5,ysize=6.5,/helvetica plot,jday(i2000),abs((gpp_model(i2000)*0.47)-npp(i2000))/npp(i2000)*100.,xtitle='Day of year',ytitle='Residual NPP (g m!u-2!n d!u-1!n)',yrange=[0,100],psym=1 oplot,jday(i2001),abs((gpp_model(i2001)*0.47)-npp(i2001))/npp(i2001)*100.,linestyle=1,psym=5 oplot,jday(i2002),abs((gpp_model(i2002)*0.47)-npp(i2002))/npp(i2002)*100.,color=125,psym=8 legend,['2000','2001','2002'],color=[0,0,125],linestyle=[0,1,0],psym=[1,5,8],/left,box=0,charsize=1.5 device,/close; close ps file ;return plotting to unix set_plot,'x' !p.charthick=1 !p.thick=1 !x.thick=1 !y.thick=1 !p.charsize=1 ;output LAI estimates for Jon (year,jday,WC,aspen) lai_output=fltarr(4,n_elements([jday(i2000),jday(i2001),jday(i2002),jday(i2003),jday(i2004)])) lai_output(0,*)=[year(i2000),year(i2001),year(i2002),year(i2003),year(i2004)] lai_output(1,*)=[jday(i2000),jday(i2001),jday(i2002),jday(i2003),jday(i2004)] lai_output(2,*)=[lai00_smooth,lai01_smooth,lai02_smooth,lai03_smooth,lai04_smooth] lai_output(3,*)=[lai_asp00_smooth,lai_asp01_smooth,lai_asp02_smooth,lai_asp03_smooth,lai_asp04_smooth] fname='lai_output.txt' print,'Saving '+fname+'...' openw,lun1,filepath(fname,root_dir='/eddy/s2/bcook/caterpillar/'),/get_lun printf,lun1,lai_output close,lun1 free_lun,lun1 ;output Willow Creek LAI, fIPAR, fAPAR, LUEmax, LUEobs, CUE for Faith Ann ;wc_output=fltarr(8,n_elements([jday(i2000),jday(i2001),jday(i2002)])) ;wc_output(0,*)=[year(i2000),year(i2001),year(i2002)]; year ;wc_output(1,*)=[jday(i2000),jday(i2001),jday(i2002)]; jday ;wc_output(2,*)=[lai00_smooth,lai01_smooth,lai02_smooth]; LAI ;wc_output(3,*)=[fipar00_smooth,fipar01_smooth,fipar02_smooth]; fIPAR (does not account for stems, bole, etc) ;wc_output(4,*)=[fapar00_smooth,fapar01_smooth,fapar02_smooth]; fAPAR (only intercepted by leaf) ;wc_output(5,*)=[luemax(i2000),luemax(i2001),luemax(i2002)]; LUEmax (low VPD, clear skies, no LAI factor) ;wc_output(6,*)=[LUEobs(i2000),LUEobs(i2001),LUEobs(i2002)]; LUEobs (corrected for LAI, cloudiness) ;wc_output(7,*)=[cue00_smooth,cue01_smooth,cue02_smooth]; CUE ;fname='wc_output.txt' ;print,'Saving '+fname+'...' ;openw,lun1,filepath(fname,root_dir='/eddy/s2/bcook/caterpillar/'),/get_lun,width=200 ;printf,lun1,wc_output ;close,lun1 ;free_lun,lun1 ;create and output aspen input for SAIL model sail=fltarr(7,23*2) sail(*,*)=-999. row=0 FOR j=0,1 DO BEGIN i=0 FOR i=0,22 DO BEGIN index=where(year EQ 2000+j AND jday Gt i*16 AND jday le i*16+16) tindex=where(year EQ 2000+j AND jday GT i*16 AND jday LE i*16+16 AND tasp NE -999.) sail(0,row)=mean([year(index)]) sail(1,row)=mean([jday(index)]) sail(2,row)=mean([sz10(index)]) IF j EQ 0 THEN BEGIN IF i LT 22 THEN BEGIN sail(3,row)=mean([lai_asp00_smooth(i*16:i*16+15)]) sail(4,row)=stdev([lai_asp00_smooth(i*16:i*16+15)]) ENDIF ELSE begin sail(3,row)=mean([lai_asp00_smooth(i*16:n_elements(lai_asp00_smooth)-1)]) sail(4,row)=stdev([lai_asp00_smooth(i*16:n_elements(lai_asp00_smooth)-1)]) endelse endif ELSE BEGIN IF i LT 22 THEN BEGIN sail(3,row)=mean([lai_asp01_smooth(i*16:i*16+15)]) sail(4,row)=stdev([lai_asp01_smooth(i*16:i*16+15)]) ENDIF ELSE begin sail(3,row)=mean([lai_asp01_smooth(i*16:n_elements(lai_asp01_smooth)-1)]) sail(4,row)=stdev([lai_asp01_smooth(i*16:n_elements(lai_asp01_smooth)-1)]) endelse endelse IF i GE 5 AND i LE 19 THEN begin sail(5,row)=mean([1.-tasp(tindex)]) sail(6,row)=stdev([1.-tasp(tindex)])/sqrt(n_elements(tindex)) ENDIF row=row+1 endfor endfor fname='asp_SAIL.txt' print,'Saving '+fname+'...' openw,lun1,filepath(fname,root_dir='/eddy/s2/bcook/caterpillar/'),/get_lun,width=120 ;thisformat='(i4,tr3,i3,tr3,4(d12.4,tr1))' printf,lun1,sail close,lun1 free_lun,lun1 stop return end