; this code is used to decompose the NEE measured at WLEF using ; footprint model ; 11/04/2003, write code ; load data ;11/06 add small tower measurements ; 11/18, changed based on decompose_neeconstraint.pro ; 7/7/2004 copied from decompose_nee1118cons.pro ;7/18/2004, use pref_par as par, yyyy='2003' print,'input yyyy' read,yyyy constraint='0' ; 0--no constraint, otherwise, constraint yr2=strmid(yyyy,2,2) hh=18 ;18 ;14+4-4 hh0=10 ;10 smeanday='07' type='10' meanday=fix(smeanday) ;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 get_data,'/eddy/s3/wang/thesis/ft/decompose_flux/results/vegratio_wlef'+yyyy+'_122.txt'+type,data122,column=24,rowmax=10000L get_data,'/eddy/s3/wang/thesis/ft/decompose_flux/results/vegratio_wlef'+yyyy+'_396.txt'+type,data396,column=24,rowmax=10000L ; load flux from WLEF tower get_data,'/eddy/s3/wang/thesis/ablbudget/data/WLEF'+yyyy+'_flux.txt',flx,column=42,rowmax=10000L ;LST f30=flx(9,*) ; 30m turbulent fluxes f122=flx(10,*) ; 122m f396=flx(11,*) ; 396m ;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,*) ; stoarge 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) ; stor flux is small usually, if it is badval, just ignore it. par_wc=wcreekmet(55,*) ; par from Willow creek .... 0.5 hour data jday_wc=wcreekflux(4,*) hour_wc=wcreekflux(3,*) ;rnee=where((nee_wc LE -2 AND par_wc LE 50) OR nee_wc GE 10,count) rnee=where(abs(nee_wc) GE 50 OR (par_wc LE 100 AND nee_wc GE 10),count) IF count GE 1 THEN nee_wc(rnee)=badval openw,16,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/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,*) eddy_lc=lcreekflux(7,*) nee_lc=lcreekflux(8,*) ; nee from Lost creek tower ;;IF yyyy NE 2003 THEN rr=where(stor_lc LE -990) & nee_lc(rr)=eddy_lc(rr) ; see comments for WC par_lc=lcreekmet(32,*) ; par from Lost creek .... 0.5 hour data jday_lc=lcreekflux(4,*) hour_lc=lcreekflux(3,*) rnee=where((nee_lc LE -2 AND par_lc LE 50) OR abs(nee_lc) GE 50,count) IF yyyy EQ 2001 THEN rnee=where((nee_lc LE -2 AND par_lc LE 50) OR nee_lc GE 10,count) ;rnee=where((par_lc GT 100 AND nee_lc GE 10) OR abs(nee_lc) GE 50,count) IF count GE 1 THEN nee_lc(rnee)=badval openw,18,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/lc_para'+smeanday+yyyy+'.txt' ENDIF ;44444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444 ; load wc and lc, wlef monthly fitted parameters load,'/eddy/s3/wang/thesis/ft/application/nee_par/results/wc_para_month.txt',wc_para load,'/eddy/s3/wang/thesis/ft/application/nee_par/results/lc_para_month.txt',lc_para load,'/eddy/s3/wang/thesis/ft/application/nee_par/results/wlef_para_month.txt',wlef_para data=[[data30],[data122],[data396]] flux=[[f30],[f122],[f396]] ; ;printf,14,met(0:5,k),a50,a95,s50,s95,par1(k),co2(k),ratio,format='(3I5,f6.1,I5,f10.4,6e15.4,12f11.5)' mon=data(1,*) day=data(2,*) hour=data(3,*) jday=data(4,*) par=data(10,*) ; parfit=prepar(3,*) ;preferred par parfit=[[parfit],[parfit],[parfit]] par=parfit ;par=prepar(3,*); preferred par nee=data(11,*) ; NEE ;;;;nee=flux(*); turbulent fluxes ;nee_122=data122(11,*) ;nee_396=data396(11,*) ;rnee=where(abs(nee) GE 20,count) ;IF count GE 1 THEN nee(rnee)=badval rnee=where((nee LE -2 AND par LE 50) OR nee GE 20+30,count) IF count GE 1 THEN nee(rnee)=badval ;rnee=where((nee_122 LE -2 AND par LE 50) OR nee_122 GE 20,count) ;IF count GE 1 THEN nee_122(rnee)=badval ;rnee=where((nee_396 LE -2 AND par LE 50) OR nee_396 GE 20,count) ;IF count GE 1 THEN nee_396(rnee)=badval n=len(hour) nx=col(data) ratio=fltarr(12,n) ;ratio122=fltarr(12,n) ;ratio396=fltarr(12,n) FOR j=0L,n-1 DO BEGIN ratio(*,j)=data(12:nx-1,j) ; ratio122(*,j)=data122(12:nx-1,j) ; ratio396(*,j)=data396(12:nx-1,j) ENDFOR ; combine two data ;ratio=[[ratio],[ratio122],[ratio396]] ;nee=[[nee],[nee_122],[nee_396]] ;par=[[par],[par],[par]] ;jday=[[jday],[jday],[jday]] ;hour=[[hour],[hour],[hour]] ;mon=[[mon],[mon],[mon]] rr=where(par GE 100 AND nee GE 5,count) IF count GE 1 THEN nee(rr)=badval help,ratio help,nee help,par ;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 ; using some-day data to fit openw,14,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/decomposedpara'+type+constraint+smeanday+yyyy+'.txt' openw,20,'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/wlef_para'+smeanday+yyyy+'.txt' FOR day1=0,365-1-meanday*0 DO BEGIN ; decompose=fltarr(3*nk)+badval result_wc=fltarr(3)+badval result_lc=fltarr(3)+badval result_wlef=fltarr(3)+badval result_mix=fltarr(3)+badval r=where(jday GE day1 AND jday Lt day1+meanday AND ratio(0,*) NE badval AND par gt 0 AND nee NE badval AND mon GE 5 AND mon LE 10 AND (hour le hh AND hour GE hh0),count) IF count LE 19 THEN BEGIN status=-1 GOTO, contd ENDIF ;find inital values using mixed measurements start=[-0.05D,-15.0,1.0] x_mix=par(r) y_mix=nee(r) ratio_mix=ratio1(*,r) err=0.2d*abs(median(y_mix)) pi=replicate({fixed:0.0,limited:[0,0],limits:[0.0D,0.0D],mpside:2},3) pi(2).limited(0)=1 ; activate lower boundary pi(2).limits(0)=0.0D ; lower boundary value ;pi(2).limited(1)=1 ;pi(2).limits(1)=10.0D ;pi(1).limited(0)=1 ;pi(1).limits(0)=-40D pi(1).limited(1)=1 ; activate upper boundary pi(1).limits(1)=0.0D ; upper boundary value result_mix=mpfitfun('myfunction1',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 rscreen=where(abs(myfunction1(x_mix,result_mix)-y_mix) ge 10+90, count) 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 gt 0 AND nee_wc NE badval,count); AND (hour_wc le hh AND hour_wc GE hh0),count) IF count GE 20 THEN begin xwc=par_wc(rwc) ywc=nee_wc(rwc) result_wc=mpfitfun('myfunction1',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 gt 0 AND nee_lc NE badval,count); AND (hour_lc le hh AND hour_lc GE hh0),count) IF count GE 20 THEN begin xlc=par_lc(rlc) ylc=nee_lc(rlc) result_lc=mpfitfun('myfunction1',double(xlc),double(ylc),err,start,yfit=yfit,perror=p_lc,parinfo=pi,/quiet) endif ENDIF ; ;444444444444444444444444444444444444444444444444444444444 ;result_mix=result_wc ; add 7/13/2004 start=result_wlef IF result_mix(0) LE -0.2 OR result_mix(0) GE 0 $ OR result_mix(1) GE 0 OR result_mix(1) LE -100 $ OR result_mix(2) GE 20 OR result_mix(2) LE 0 THEN start=wlef_para(1:3,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},3*nk) GOTO,ccc ;-----add constraints------ 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,ratio1,result_mix,parratio_con,nee_con addconstraint,ratio1,start,parratio_con,nee_con ;7/12/2004 use a new constraint ;;;addconstraint1,ratio_screen,start,par_screen,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 error=fltarr(len(sxnew))+1d*err number_con=len(parratio_con) adjust=2.5d ;1.0D error(len(error)-1-number_con+1:len(error)-1)=adjust*err ; using errors as weight to control print,'total numbers=',len(error),' constraint points=',number_con ;result=mpfitfun('myfunction7',[sx(0,r),sx(1,r),sx(2,r),sx(3,r),sx(4,r),sx(5,r),sx(6,r),sx(7,r)],1d*nee(r),err,start7,parinfo=pi7,perror=p1,bestnorm=best, covar=cc,status=status,/quiet,weights=1d/err^2) ;help,sx result=mpfitfun('myfunction6',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,'jday=',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,18E15.4)' printf,20,month,day1+1,result_mix(*),format='(2I5,3E15.4)' IF fix(yyyy) GE 1999 THEN printf,16,month,day1+1,result_wc(*),format='(2I5,3E15.4)' IF fix(yyyy) GE 2001 THEN printf,18,month,day1+1,result_lc(*),format='(2I5,3E15.4)' ENDFOR close,14 close,16 close,18 close,20 print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/decomposedpara'+type+constraint+smeanday+yyyy+'.txt' print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/wlef_para'+smeanday+yyyy+'.txt' IF fix(yyyy) GE 1999 THEN print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/wc_para'+smeanday+yyyy+'.txt' IF fix(yyyy) GE 2001 THEN print,'saving '+'/eddy/s3/wang/thesis/ft/decompose_flux/results/nee_par/lc_para'+smeanday+yyyy+'.txt' end