; this code is used to decompose the NEE measured at WLEF using ; footprint model ;copied from decompose_nee.pro ; NEE~ temperature at night yyyy='2001' print,'input yyyy' read,yyyy constraint='0' ; 0--no constraint, otherwise, constraint yr2=strmid(yyyy,2,2) hh=18 ;14+4-4 smeanday='28' type='10' meanday=fix(smeanday) adjust=1.0d IF fix(yyyy) EQ 2000 THEN adjust=0.5d; 1.0d ;0.4d IF fix(yyyy) EQ 2001 THEN adjust=0.5d; 0.5d ;0.4d IF fix(yyyy) EQ 2003 THEN adjust=2.d; 2.5d ; 1.0d ;0000000000000000000000000000000000000000000000000 !p.multi=[0,2,2] badval=-999. ;1111111111111111111111111111111111111111111111111111111111111111111111111111111111111 ; Load footprint data - weights get_data,'/eddy/s3/wang/thesis/ft/decompose_flux/results/vegratio_wlef'+yyyy+'_30.txt'+type,data30,column=24,rowmax=10000L ; load WLEF met data get_data,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yyyy+'_met.txt',wlefmet,column=44,rowmax=10000L ;LST ;222222222222222222222222222222222222222222222222222222222222222222222222222222222222222 ;load preferred PAR get_data,'/eddy/s3/wang/thesis/ablbudget/data/pref_par/prefpar'+yr2+'.txt',prepar,column=5,rowmax=10000L ;;; load Willow Creek and Lost creek data for some years ;3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333 IF fix(yyyy) GE 1999 THEN BEGIN ; after 1999, we can use WC data as contraints get_data,'/eddy/s1/cheas/flux/wcreek/wcreek'+yyyy+'_flux.txt',wcreekflux,column=13,rowmax=20000L get_data,'/eddy/s1/cheas/flux/wcreek/wcreek'+yyyy+'_met.txt',wcreekmet,column=69,rowmax=20000L stor_wc=wcreekflux(6,*) ; storage flux eddy_wc=wcreekflux(7,*) ; eddy flux nee_wc=wcreekflux(8,*) ; nee from Willow creek tower rr=where(stor_wc LE -990) & nee_wc(rr)=eddy_wc(rr) ; use eddy flux as NEE due to bad storage flux nee_flag=wcreekflux(12,*) ; quality of flux,flag=0---OK par_wc=wcreekmet(55,*) ; par from Willow creek .... 0.5 hour data jday_wc=wcreekflux(4,*) hour_wc=wcreekflux(3,*) airt_wc=wcreekmet(30,*) ; Temperature at 97ft , C soilt_wc=wcreekmet(40,*) ; temperature 5cm below the surface ttt_wc=airt_wc*0.0+badval rr=where(airt_wc NE badval AND soilt_wc NE badval) ttt_wc(rr)=0.5*(airt_wc(rr)+soilt_wc(rr)) ttt_wc=airt_wc ;9/28/2004 ;;;;ttt_wc=soilt_wc openw,16,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/wc_para'+smeanday+yyyy+'.txt' ENDIF IF fix(yyyy) GE 2001 THEN BEGIN ; after 2001, we can use LC data as contraints get_data,'/eddy/s1/cheas/flux/lcreek/lcreek'+yyyy+'_flux.txt',lcreekflux,column=12,rowmax=20000L get_data,'/eddy/s1/cheas/flux/lcreek/lcreek'+yyyy+'_met.txt',lcreekmet,column=35,rowmax=20000L stor_lc=lcreekflux(6,*) ; storage eddy_lc=lcreekflux(7,*) ; eddy flux nee_lc=lcreekflux(8,*) ; nee from Lost creek tower rr=where(stor_lc LE -990) & nee_lc(rr)=eddy_lc(rr) ; in some cases, too mnay badvals in NEE because of bad storage fluxes flux_lc=lcreekflux(7,*) ; turbulent flux rr=where(nee_lc EQ badval) & nee_lc(rr)=flux_lc(rr) ; storage flux is small. Most winter data, NEE is badval due to no storage fluxes which are so small. par_lc=lcreekmet(32,*) ; par from Lost creek .... 0.5 hour data jday_lc=lcreekflux(4,*) hour_lc=lcreekflux(3,*) airt_lc=lcreekmet(13,*) ; Temperature at 10.2m, C soilt_lc=lcreekmet(22,*) ; temperature 5cm below the surface ttt_lc=airt_lc*0.0+badval rr=where(airt_lc NE badval AND soilt_lc NE badval) ttt_lc(rr)=0.5*(airt_lc(rr)+soilt_lc(rr)) ttt_lc=airt_lc ;9/28/04 ;;;;ttt_lc=soilt_lc openw,18,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/lc_para'+smeanday+yyyy+'.txt' ENDIF ;44444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444 ; load wc and lc, wlef monthly fitted parameters load,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/nee_t_wc_month.txt',wc_para load,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/nee_t_lc_month.txt',lc_para load,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/nee_t_wlef_month.txt',wlef_para data=data30 mon=data(1,*) day=data(2,*) hour=data(3,*) jday=data(4,*) ;;;par=data(10,*) par=prepar(3,*) ;preferred par; here used to judge if it is at night, i.e., Par LE 0 nee=data(11,*) airt=wlefmet(21,*) ; 21-30m-airT C ; 20- 122m soilt=wlefmet(23,*) ; soil temperature 5cn below the surface ttt=soilt*0+badval ;rr=where(soilt NE badval AND airt NE badval) ;ttt(rr)=0.5*(soilt(rr)+airt(rr)) ;;ttt=soilt ttt=airt ; 9/28/2004 changed n=len(hour) nx=col(data) ratio=fltarr(12,n) FOR j=0L,n-1 DO BEGIN ratio(*,j)=data(12:nx-1,j) ENDFOR ;need to regroup, combine some small contributions nk=7-1 ; type of veg,1108 combine mixed con/dec with pine+mixed coniferous ratio1=fltarr(nk,len(ratio))+badval FOR j=0,len(ratio)-1 DO BEGIN IF ratio(0,j) NE badval THEN begin ratio1(0,j)=ratio(3,j)+ratio(4,j)+ratio(7,j) ; pine+mixed coniferous ratio1(1,j)=ratio(5,j) ; aspen ratio1(2,j)=ratio(6,j) ; mixed deciduous ; ratio1(3,j)=ratio(7,j) ; mixed con/dec ratio1(3,j)=ratio(8,j) ; lowland wetland ratio1(4,j)=ratio(9,j) ; forested wetland ratio1(5,j)=ratio(0,j)+ratio(2,j)+ratio(10,j)+ratio(11,j) ; others excluding water ENDIF ENDFOR ;==================================start fitting ======================================= ; using some-day data to fit openw,14,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/decomposedpara'+type+constraint+smeanday+yyyy+'.txt' openw,20,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/wlef_para'+smeanday+yyyy+'.txt' FOR day1=0,365-1-meanday*0 DO BEGIN ; decompose=fltarr(2*nk)+badval ; two parameters flux=a*exp(b*T) for one type result_wc=fltarr(2)+badval result_lc=fltarr(2)+badval result_wlef=fltarr(2)+badval ; added on 7/15/2004 result_mix=fltarr(2)+badval r=where(jday GE day1 AND jday Lt day1+meanday AND ratio(0,*) NE badval AND par le 0 AND abs(nee) LE 30+10 AND nee GE -999 AND ttt NE badval,count) IF count LE 19 THEN BEGIN status=-1 GOTO, contd ENDIF ;find inital values using mixed measurements start=[1.0D,0.05] ;[a,b] x_mix=transpose(ttt(r)) y_mix=nee(r) ratio_mix=ratio1(*,r) help,x_mix help,y_mix err=0.2d*abs(median(y_mix)) pi=replicate({fixed:0.0,limited:[0,0],limits:[0.0D,0.0D],mpside:2},2) pi(1).limited(0)=1-1 ; activate lower boundary pi(1).limits(0)=0.0D ; lower boundary value for b pi(0).limited(0)=1-1 ; activate upper boundary pi(0).limits(0)=0.0D ; lower boundary value for a result_mix=mpfitfun('func_nee_t',double(x_mix),double(y_mix),err,start,yfit=yfit,parinfo=pi,/quiet);,weights=1.0D) result_wlef=result_mix ;======== ;screen the data based on the fitted line, if the data are within the ;10-umol, OK, if not, throw away. plot,x_mix,y_mix,psym=2 oplot,x_mix,yfit,psym=4 legend,['a='+strtrim(string(result_mix(0)),2),'b='+strtrim(string(result_mix(1)),2)],linestyle=[-1,-1],/left rscreen=where(abs(func_nee_t(x_mix,result_mix)-y_mix) ge 10+90, count) ; actually no use here, since 10+90 too big IF count GE 1 THEN y_mix(rscreen)=badval rgood=where(y_mix NE badval,count) IF count GE 1 THEN BEGIN par_screen=x_mix(rgood) nee_screen=y_mix(rgood) ratio_screen=ratio_mix(*,rgood) ENDIF ;========= ; find the relationship from WC and LC towers ;;; ;;result_wc=wc_para(1:3,mon(r(0))-1) ; use the monthly mean result as default values ;;result_lc=lc_para(1:3,mon(r(0))-1) ;;; ;44444444444444444444444444444444444444444444444444444444444444 IF fix(yyyy) GE 1999 THEN BEGIN ; after 1999, we can use WC observtional data directly rwc=where(jday_wc GE day1 AND jday_wc Lt day1+meanday AND ((par_wc le 0 AND par_wc NE badval) OR hour_wc GE 22 OR hour_wc LE 5) AND abs(nee_wc) LE 30 AND nee_wc GE -2 AND nee_flag EQ 0 AND ttt_wc NE badval,count) IF count GE 20 THEN begin xwc=transpose(ttt_wc(rwc)) ywc=nee_wc(rwc) result_wc=mpfitfun('func_nee_t',double(xwc),double(ywc),err,start,yfit=yfit,perror=p_wc,parinfo=pi,/quiet) endif ENDIF ; ;4444444444444444444 IF fix(yyyy) GE 2001 THEN BEGIN ; after 2001, we can use LC observtional data directly rlc=where(jday_lc GE day1 AND jday_lc Lt day1+meanday AND ((par_lc le 0 AND par_lc NE badval) OR hour_lc GE 22 OR hour_lc LE 5) AND abs(nee_lc) LE 30 AND nee_lc GE -2 AND ttt_lc NE badval,count) IF count GE 20 THEN begin xlc=transpose(ttt_lc(rlc)) ylc=nee_lc(rlc) result_lc=mpfitfun('func_nee_t',double(xlc),double(ylc),err,start,yfit=yfit,perror=p_lc,parinfo=pi,/quiet) endif ENDIF ; ;444444444444444444444444444444444444444444444444444444444 result_mix=result_wlef ; add 7/13/2004 start=result_wlef ;start=result_wc IF result_mix(0) LE 0.0 OR result_mix(0) GE 5 $ OR result_mix(1) LE 0 OR result_mix(1) ge 0.2 THEN start=wlef_para(1:2,mon(r(0))-1);wlef_para(1:2,mon(r(0))-1) IF constraint EQ '0' THEN start7=[start,start,start,start,start,start] ; 11/11/2003 added; 11/18, if used wc and lc measurements as constraint, start7 needs to change pi7=replicate({fixed:0.0,limited:[0,0],limits:[0.0D,0.0D],mpside:0+2},2*nk) GOTO,ccc ;-----add constraints------ ; Have not changed yet, because we skip it now FOR k=1,nk DO BEGIN pi7(3*k-1).limited(0)=1 pi7(3*k-1).limits(0)=0.0D pi7(3*k-1).limited(1)=1 pi7(3*k-1).limits(1)=10.0D pi7(3*k-3).limited(0)=1 pi7(3*k-3).limits(0)=-0.12d;-5d*abs(start7(3*k-3)); -0.12d;abs(result_mix(0));-0.1D pi7(3*k-3).limited(1)=1 pi7(3*k-3).limits(1)=0.0D pi7(3*k-2).limited(1)=1 pi7(3*k-2).limits(1)=0.0D pi7(3*k-2).limited(0)=1 pi7(3*k-2).limits(0)=-100D;-5*abs(start7(3*k-2));-100.0;start(1)*2 ENDFOR ccc: ;-- end of constraints ;;;sx=1d*[par,ratio1] sx=1d*[transpose(par_screen),ratio_screen] addconstraint_nee_t,ratio1,start,parratio_con,nee_con help,sx help,parratio_con help,transpose(nee_screen) help,nee_con sxnew=[[sx],[parratio_con]] neenew=[[transpose(nee_screen)],[nee_con]] neenew=transpose(neenew) ; construct error (or say weights) error=fltarr(len(sxnew))+1d*err number_con=len(parratio_con) ;;;adjust=0.4; 10/3/04; 0.4; 1.0D error(len(error)-1-number_con+1:len(error)-1)=adjust*err ; using errors as weight to control print,len(error),number_con result=mpfitfun('func_nee_t6',double(sxnew),double(neenew),double(error),start7,parinfo=pi7,perror=p1,bestnorm=best, covar=cc,status=status,xtol=1.e-8,niter=niter,maxiter=1000,/quiet);;;;,weights=1d/(1d*nee(r))) ;print,len(x),' status=',status ;print,'bestnorm=',best print,'error=',p1 print,'p=',result print,'result_mix',result_mix print,' niter=',niter,' data points',len(rgood) print,'start=',start ;print,'start=',start print,'result_wc=',result_wc IF fix(yyyy) GE 2001 THEN print,'result_lc=',result_lc IF status gt 0 AND niter GE 2 THEN decompose=result ; s='' ; print,'continue?' ; read,s contd: print,day1, ' status=',status month=strmid(strtrim(string(jdatetodate(day1+1,fix(yyyy))),2),4,2) printf,14,month,day1+1,status,decompose(*),format='(2I5,I5,12E15.4)' printf,20,month,day1+1,result_wlef(*),format='(2I5,2E15.4)' IF fix(yyyy) GE 1999 THEN printf,16,month,day1+1,result_wc(*),format='(2I5,2E15.4)' IF fix(yyyy) GE 2001 THEN printf,18,month,day1+1,result_lc(*),format='(2I5,2E15.4)' ENDFOR close,14 close,16 close,18 close,20 print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/decomposedpara'+type+constraint+smeanday+yyyy+'.txt' print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/wlef_para'+smeanday+yyyy+'.txt' IF fix(yyyy) GE 1999 THEN print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/wc_para'+smeanday+yyyy+'.txt' IF fix(yyyy) GE 2001 THEN print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_t/lc_para'+smeanday+yyyy+'.txt' end