/*
This program calculates planar fit coordinate rotation coefficients at monthly seasonal and annual scale
OUTPUT: 
	1. datafile with the mean and standard errors of the estimates at different time steps
	2. graphs of monthly coefficients and their standard errors

07/15/2008, Asko Noormets, anoorme@ncsu.edu
*/

****************************************************************************;
*	EDITABLE PARAMETERS														;
****************************************************************************;

%let path=DIRECTORY_NAME_HERE_like_C:\DATA\;
%let dsn=NAME_OF_ANNUAL_ONLINE_FLUX_FILE_HERE;

****************************************************************************;
*	IMPORTANT!																;
*	please confirm that the variable names in lines 33-50 correspond to		;
*	the locations where the parameters are in your file						;
*																			;
****************************************************************************;
data a1;
infile "&path.&dsn" dsd firstobs=6 lrecl=2000;
input dhms :$22. v1-v79; 
	sasdate=input(scan(dhms,1,' '),yymmdd10.); 
	sastime=input(scan(dhms,2,' '),time12.); 
	sasdt=dhms(sasdate,0,0,sastime); format sasdt datetime22.1; 
	month=month(sasdate); 
/********************************************************/
/********** Data SCREENING criteria *********************/
/********************************************************/
		if v1<-10 or v1>10 then v1bad=1; else v1bad=0; 			* Fc, mg m^-2 s^-1 ;
		if v2<-200 or v2>800 then v2bad=1; else v2bad=0; 		* LE ;
		if v3<-300 or v3>900 then v3bad=1; else v3bad=0; 		* Hs ;
		if v5<0 or v5>6 then v5bad=1; else v5bad=0; 			* ustar ;
		if v27<0 or v27>1200 then v27bad=1; else v27bad=0; 		* co2_ave, mg m^-3 ;
		if v28<0 or v28>40 then v28bad=1; else v28bad=0; 		* h2o_ave, g m^-3 ;
		if v29<-30 or v29>50 then v29bad=1; else v29bad=0; 		* Tsonic ;
		if v30<1 or v30>1.5 then v31bad=1; else v31bad=0; 		* rho;
		if v58<-50 or v58>50 then v34bad=1; else v34bad=0; 		* ta1;
		if v31<80 or v31>110 then v42bad=1; else v42bad=0; 		* pressure;
		if v59<-400 or v59>1000 then v44bad=1; else v44bad=0; 	* Rn;
		if v35<0 or v35>20 then v52bad=1; else v52bad=0; 		* wind speed ;
		if v37<10.5 then v54bad=1; else v54bad=0; 				* battery voltage ;
		if v50>80 then v79bad=1; else v79bad=0; 				* AGC ;
	if (v1bad=1 or v2bad=1 or v3bad=1 or
		v5bad=1 or v27bad=1 or v28bad=1 or v29bad=1  or v31bad=1 or 
		v34bad=1 or v42bad=1 /*or v44bad=1*/ or v52bad=1 or v54bad=1 or v79bad=1
		or v40>500 or v41>0 or v39<17000)
		then gooddata=0; else gooddata=1;						* v40, v41 = CSAT & IRGA warnings ;
	run; 

/* check if wind speed is sensitive to direction, esp. when from behind the tower 100-190 degrees */
proc gplot data=a1; 
plot (v24 v25 v26)*v50 / overlay; 
title "&dsn.; u, v, w, 2005"; symbol value=plus; goptions interpol=none; run;

/* step 2. sort by month, calculate b0, b1, b2 for different time windows & save */
proc sort data=	a1 nodupkey; by month sasdt; run; 
proc reg data=a1 outest=regyear edf outseb; where  gooddata=1; 
title "1 of 3 &dsn.; PLANAR FIT COEFFICIENTS 2007, FULL SEASON";
model v26= v24 v25; run;
proc reg data=a1 outest=regseas edf outseb; where  gooddata=1 and 5<month<9; 
title "2 of 3 &dsn.; PLANAR FIT COEFFICIENTS 2007, JUN-AUG";
model v26= v24 v25; run;
proc reg data=a1 outest=regmonth edf outseb; where  gooddata=1; by month; 
title "3 of 3 &dsn.; PLANAR FIT COEFFICIENTS 2007 BY MONTH";
model v26= v24 v25; run;

data regs; set regyear regseas regmonth;
file "&path.REGPARMS_&dsn.reg";
put _ALL_;
run;

/*plot parameter estimates and error, check for representativeness*/
proc gplot data=regs; where _TYPE_="PARMS";
title "&dsn.; PARAMETERS";
plot (Intercept v24 v25)* month / overlay;
symbol value=dot; goptions interpol=join;
run;
proc gplot data=regs; where _TYPE_="SEB" ;
title "&dsn.; PARAM-ERR";
plot (Intercept v24 v25)* month / overlay;
symbol value=dot; goptions interpol=join;
run;