/* determine critical u* */
/* 1. based on real data */
/* 2. based on binned data */
/* use same infile as that for AQRTa, AQRTaSWC model fitting */
/* CONCLUSION:
	tried 	raw Fc_leuning, raw NEE vs. u* 
				by season, spring-fall (excluding winter)
			binned Fc_leuning, binned NEE vs. binned u* 
				by season, by month, spring-fall
	post-dormancy had higher u*crit - 0.35 (30-35% eliminated)
	growing season and pre-dormancy had lower u*crit - 0.12 (20-35% eliminated)
	winter flux roughly flat, reasonable in March-April, October-November
	===========================================================================
		Jan, Feb, Dec fluxes mostly negative, even at night, not reasonable
	===========================================================================
*/
/* original: 08.11.2005 */
/* modified: 01.30.2007 */
/* Asko */

%let site=cc;
%let year=2007;
%let path=C:\DATA\&site.\;
%let dsn=Flux_&site._&year._Master.txt;
%let timestep=month; 			/*used for developing AQRT models, either monthly or weekly or seasonally*/
%let stability=0<stability<7;	/*range of stability classes to consider*/
*%let method=rfw; 				/* rfw, r1w or wpl */

/* uses Master dataset, which was developed in "Combine_Clean_OO_20042005.sas" */
/* it already includes the quality flags */

data flux1; infile "&path.&dsn." lrecl=4000 dsd;
rename r15=Fcrfw r14=LErfw r1=ustarrfw;
	input sasdt :datetime13. year month day hour
				r1-r21
				ts5 ts20 ta4 rh4 m5 m20 par cs616 vwc gashound hft1 hft2 hft3 
				Fcwpl Lewpl Hs tau ustar v1-v74
				L_obukhov turb_u_model turb_w_model turb_u_meas turb_w_meas itc2_u itc2_w stability
				vpd1 vpd2 co2_stor T_stor gooddata;
		format sasdt datetime13.; sasdate=datepart(sasdt); sastime=timepart(sasdt);
	PARa=v74;
run;

	data flux; merge flux1 flux1(firstobs=2 rename=(Fcrfw=leadFc) keep=Fcrfw);
		lagfc=lag(Fcrfw);
		jumpfc=Fcrfw-lagfc;
	/* choose whether to flag Fc[i] or Fc[i-1]*/
		if 	(abs(jumpfc)>20 and abs(Fcrfw)>abs(lagfc)) or 
			(abs(jumpfc)>20 and abs(Fcrfw)>abs(leadFc)) then do; 
			Fcrfw=.; Fcrfw_bad=1; gooddata=0; end;
	/* calc NEE here rather than in 1st step since Fc has been modified 'til now*/
	/*if Gashound signal varies too much, then magnitude of co2_stor uncertain, --> don't use*/
		if co2_stor ne . then do; NEE=Fcrfw+co2_stor; end;
		if co2_stor=. then do; NEE=Fcrfw; end;
		storrat=co2_stor/NEE; *if abs(storrat)>3 then storrat=.;
		/*This is remnant from prior version*/
		*if season ne "4dor" then year=2004;
	run;


/************************************************************************************/
/* Determine u*rit using binned u* data, since raw data shows no distinct threshold */
/*																					*/
/************************************************************************************/

proc sort data=flux; by year &timestep.; run;
	proc univariate data=flux noprint; by year &timestep.;
	where PARa<4 and gooddata ne 0 and itc2_w<0.3 and r17<2 and r16>17000 and &stability. and v45<75; 
	var ustarrfw; output out=ustar_p10 
		pctlpts= 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95  pctlpre=p_;
	run;
	data flux_1; merge flux ustar_p10; by year &timestep.;
	if 0 ge ustarrfw and ustarrfw<p_5 then ustarperc=5; 
	else if p_5 ge ustarrfw and ustarrfw<p_10 then ustarperc=10; 
	else if p_10 ge ustarrfw and ustarrfw<p_15 then ustarperc=15; 
	else if p_15 ge ustarrfw and ustarrfw<p_20 then ustarperc=20; 
	else if p_20 ge ustarrfw and ustarrfw<p_25 then ustarperc=25; 
	else if p_25 ge ustarrfw and ustarrfw<p_30 then ustarperc=30;
	else if p_30 ge ustarrfw and ustarrfw<p_35 then ustarperc=35; 
	else if p_35 ge ustarrfw and ustarrfw<p_40 then ustarperc=40; 
	else if p_40 ge ustarrfw and ustarrfw<p_45 then ustarperc=45; 
	else if p_45 ge ustarrfw and ustarrfw<p_50 then ustarperc=50; 
	else if p_50 ge ustarrfw and ustarrfw<p_55 then ustarperc=55; 
	else if p_55 ge ustarrfw and ustarrfw<p_60 then ustarperc=60;
	else if p_60 ge ustarrfw and ustarrfw<p_65 then ustarperc=65; 
	else if p_65 ge ustarrfw and ustarrfw<p_70 then ustarperc=70; 
	else if p_70 ge ustarrfw and ustarrfw<p_75 then ustarperc=75; 
	else if p_75 ge ustarrfw and ustarrfw<p_80 then ustarperc=80; 
	else if p_80 ge ustarrfw and ustarrfw<p_85 then ustarperc=85; 
	else if p_85 ge ustarrfw and ustarrfw<p_90 then ustarperc=90;
	else if p_90 ge ustarrfw and ustarrfw<p_95 then ustarperc=95; 
	else if p_95 ge ustarrfw and ustarrfw<1.5 then ustarperc=100;
	else ustarperc=.;
	run;
	proc sort data=flux_1; by year &timestep. ustarperc; run;
	proc means data=flux_1 mean stderr noprint; by year &timestep. ustarperc; 
	where PARa<4 and gooddata ne 0 and itc2_w<0.3 and r17<2 and r16>17000 and &stability. and v45<75 and (sastime>75600 or sastime<21600); 
	var ustarrfw Fcrfw storrat;
	output out=meanustar mean=ustar Fc storrat stderr=e_ustar e_Fc e_storrat;
	run;
	proc sort data=meanustar; by year &timestep. ustarrfw; run;
		data err; set meanustar;
		Fc_minus=Fc_leun-e_Fc_leun; Fc_plus=Fc_leun+e_Fc_leun; 
		NEE_minus=NEE-e_NEE; NEE_plus=NEE+e_NEE; run;

		goptions ftext=swiss htext=2;
		proc gplot data=err ; by year &timestep.; 
		symbol1 value=dot interpol=join co=black ci=black;
		symbol2 value=point interpol=join co=blue ci=blue repeat=2;
		symbol4 value=dot interpol=join co=black ci=black;
		symbol5 value=point interpol=join co=blue ci=blue repeat=2;
		goptions interpol=join; symbol value=dot;
		plot (Fc_leun Fc_minus Fc_plus)*ustarrfw / overlay /*vaxis=-2 0 2 4 6*/ ; 
		*plot (NEE NEE_minus NEE_plus)*ustarrfw / overlay grid vaxis=-2 0 2 ; 
		run;
	/* check threshold for raw data using growing season data only (that was reasonable in binned analysis)*/
		proc gplot data=flux_1; symbol value=plus i=none; 
		where PARa<4 and gooddata ne 0 and itc2_w<0.3 and r17<2 and r16>17000 and &stability. and v45<75 and (sastime>75600 or sastime<21600);
		plot (Fcrfw )*ustarrfw / overlay /*vaxis=-4 -2 0 2 4 6 8*/; 
		plot (storrat )*ustarrfw / overlay grid vaxis=-2 0 2 ; 
		run;

		data ustarout; set meanustar;
		file "&path.ustarrfw_bin_&dsn.&timestep..txt";
		if _n_=1 then put
		'year month ustarclass ustar Fc storrat err_ustar err_Fc err_storrat';
		put year &timestep. ustarperc ustar Fc storrat e_ustar e_Fc e_storrat; run;
/****************************************************************************************/
/* binned u* calculation over 															*/
/****************************************************************************************/