options nonotes nosource nomprint nosymbolgen nomlogic bufsize=16M; 

/********************************************************************************************;
*																							;
* EC_Processor, v2.2																		;
* Saved on: 	Mar. 26, 2009																;
* Created by: 	© Asko Noormets, anoorme at ncsu dot edu									;
*																							;
********************************************************************************************
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
* FAIR USE POLICY:																			;
*																							;
* Please use this program in a fashion you would like your own work to be used - offering 	; 
* due credit when appropriate. Please cite http://www4.ncsu.edu/~anoorme/ECP/ as the source ;
* when using it in publication. For studies where this program plays a significant role,	;
* I expect that agreement about co-authorship will be reached before submission to a journal; 
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
*																							;
* SYSTEM REQUIREMENTS:																		;
* free space on system disk (or where SAS temporary files are located): 					;
*					~3.8 times more than the largest input file								;
*																							;
* EDITABLE PARAMETERS ARE:																	;
*																							;
* path - directory where the files are located (line 89)								 	;
* type of coordinate rotation (92)															;
* file where planar fit rotation parameters are stored (line 95)							;
* names of datafiles to be processed (from line 480 on)										; 
*																							;
*																							;
********************************************************************************************;
* Planar rotation coefficients are assigned based on information contained in the file name:;
*																							;
* site name 	- 2 symbols, starting in position 1											;
* STARTING DATE																				;
* year			- 2 symbols, starting in position 7											;
* month			- 2 symbols, starting in position 3											;
* day of month	- 2 symbols, starting in position 5											;
* ENDING DATE																				;
* year			- 2 symbols, starting in position 14										;
* month			- 2 symbols, starting in position 10										;
* day of month	- 2 symbols, starting in position 12										;
*																							;
* for example, LP013105-021805.dat															;
* for different file name structure, change the parameters of the "substr" function below	;
********************************************************************************************;
* CHANGES COMPARED TO ECP v2.1  															;
* _1  ITERATIVE DESPIKING 																	;
* _2  PERFORMS EITHER PLANAR FIT OR 2D-ROTATION												;
* _3  MAXIMIZES COVARIANCE WITHIN ± 0.5 sec													;
*																							;
* DIFFERENCES COMPARED TO AMERIFLUX PROCESSING												;
* _4  CALCULATES Hs USING RHO_A (instead of RHO_D like Ameriflux standard)					;
* 																							;
* CHANGES FROM version 2.0 to v2.1 (Oct, 2008)												;
* _1  Units corrected for the WPL term of LE (thanks to H. Guo)								;
* _2  Correction of error in the LE component of the WPL term of Fc (thanks to H. Guo)		;
* _3  Removed screening by flag variables, just replace out-of-range data w/ mean			;
* _4  Deleted extra variables that were not used to speed up run time						;
********************************************************************************************;*/

/* Enter path where datafiles are located;*/
%let path=C:\DATA\CC\CC_TSDATA_TXT_2006\;

/* Do you want to apply 2-D or planar fir rotation?;*/
%let rotation=1; 	/* 1=planar fit, 2=2-axis rotation ;*/

/* Enter the name of the file with planar fit rotation parameters */
%let rotation_parameters=RegParms_CC052506.txt;




/* Read planar rotation parameters from comma-delimited external file;*/
/* have timestamps in the format of: 16feb1995:10:56:00 */
	data rotationparms_by_collection;
	infile "&path.&rotation_parameters." dsd;
	input startdt :datetime. enddt :datetime. intercept Ux Uy;
	format startdt datetime.;
	format enddt datetime.;
	startdate=datepart(startdt); format startdate date.;
	enddate=datepart(enddt); format enddate date.;
	run;
	
%macro despike;
	options nomlogic nomprint nosymbolgen;
	%do i=1 %to 5; 
	%let old=%eval(&i.-1);
		%if &i=1 %then %do;
			proc means data=raw noprint; by sasdate halfhr; 
				var u v w co2 h2o ts_ck;
				output out=outlier 
				mean=u30 v30 w30 co230 h2o30 ts_ck30 
				std=ustd vstd wstd co2std h2ostd ts_ckstd; 
			run;
	        data filter; merge raw outlier; by sasdate halfhr;
				if abs(u-u30)>6*ustd then do; u=u30; uflag=1;  end;
		        if abs(v-v30)>6*vstd then do; v=v30; vflag=1;  end;
		        if abs(w-w30)>6*wstd then do; w=w30; wflag=1;  end;
		        if abs(co2-co230)>6*co2std then do; co2=co230; co2flag=1;  end;
		        if abs(h2o-h2o30)>6*h2ostd then do; h2o=h2o30; h2oflag=1;  end;
		        if abs(ts_ck-ts_ck30)>6*ts_ckstd then do; ts_ck=ts_ck30; tsckflag=1;  end;
				drop u30 v30 w30 co230 h2o30 ts_ck30 ustd vstd wstd co2std h2ostd ts_ckstd;
	        run;
			proc means data=filter noprint; by sasdate halfhr;  
				var uflag vflag wflag co2flag h2oflag tsckflag ; 
				output out=countspikes_&i. sum=un&i. vn&i. wn&i. co2n&i. h2on&i. ts_ckn&i.; run;
			data _null_; set countspikes_&i.;
				call symput ("un&i.",un&i.);
				call symput ("vn&i.",vn&i.);
				call symput ("wn&i.",wn&i.);
				call symput ("co2n&i.",co2n&i.);
				call symput ("h2on&i.",h2on&i.);
				call symput ("tsn&i.",ts_ckn&i.);
			run;
		%end;
		%else %if &i ge 2 %then %do ;
			proc means data=filter noprint; by sasdate halfhr; 
				var u v w co2 h2o ts_ck;
				output out=outlier 
				mean=u30 v30 w30 co230 h2o30 ts_ck30 
				std=ustd vstd wstd co2std h2ostd ts_ckstd; 
			run;
			data filter; merge filter outlier; 
				by sasdate halfhr ;
				if abs(u-u30)>6*ustd then do; u=u30; uflag=1;  end;
		        if abs(v-v30)>6*vstd then do; v=v30; vflag=1;  end;
		        if abs(w-w30)>6*wstd then do; w=w30; wflag=1;  end;
		        if abs(co2-co230)>6*co2std then do; co2=co230; co2flag=1;  end;
		        if abs(h2o-h2o30)>6*h2ostd then do; h2o=h2o30; h2oflag=1;  end;
		        if abs(ts_ck-ts_ck30)>6*ts_ckstd then do; ts_ck=ts_ck30; tsckflag=1;  end;
				drop w30 co230 h2o30 ts_ck30 ustd vstd wstd co2std h2ostd ts_ckstd;
	        run;
			proc means data=filter noprint; by sasdate halfhr; 
				var u v w co2 h2o ts_ck uflag vflag wflag co2flag h2oflag tsckflag; 
				output out=countspikes_&i. 
					mean=um vm wm co2m h2om ts_ckm fm1-fm6
					sum=sum1-sum6 un&i. vn&i. wn&i. co2n&i. h2on&i. ts_ckn&i. ; 
			run;
			data countspikes (keep=sasdate halfhr uspike vspike wspike co2spike h2ospike tsspike); 
			merge countspikes_&old. countspikes_&i. (drop=fm1-fm6 sum1-sum6 ts_ckm co2m h2om); by sasdate halfhr;
			if un&i.=un&old. and vn&i.=vn&old. and wn&i.=wn&old. and
				co2n&i.=co2n&old. and h2on&i.=h2on&old. and ts_ckn&i.=ts_ckn&old. then do;
					rename un&i.=uspike vn&i.=vspike wn&i.=wspike co2n&i.=co2spike h2on&i.=h2ospike ts_ckn&i.=tsspike;
				end;				
			run;
			data _null_; set countspikes_&i.;
				call symput ("un&i.",un&i.); 
				call symput ("vn&i.",vn&i.);
				call symput ("wn&i.",wn&i.);
				call symput ("co2n&i.",co2n&i.);
				call symput ("h2on&i.",h2on&i.);
				call symput ("tsn&i.",ts_ckn&i.);
			run;
				%if &&un&i.=&&un&old. and &&vn&i.=&&vn&old. and &&wn&i.=&&wn&old. and
				&&co2n&i.=&&co2n&old. and &&h2on&i.=&&h2on&old. and &&tsn&i.=&&tsn&old. %then %goto endmacro;
		%end;
		%put _user_; run;
	%end;
	%endmacro: 
%mend despike;

%macro planar(dsn);
	data _null_; set rotationparms_by_collection;
	where startdate = mdy(%substr(&dsn.,3,2),%substr(&dsn.,5,2),%substr(&dsn.,7,2)) and
	enddate = mdy(%substr(&dsn.,10,2),%substr(&dsn.,12,2),%substr(&dsn.,14,2));
	call symput('b0',intercept);
	call symput('b1',Ux);
	call symput('b2',Uy);
	run; 

	%put _local_; 
	%let p31=-&b1/(sqrt(&b1*&b1+&b2*&b2+1));
	%let p32=-&b2/(sqrt(&b1*&b1+&b2*&b2+1));
	%let p33=1/(sqrt(&b1*&b1+&b2*&b2+1));
	%let cosa=sqrt(&p32*&p32+&p33*&p33);
	%let sina=&p31;
	%let cosb=&p33/(sqrt(&p32*&p32+&p33*&p33));
	%let sinb=-&p32/(sqrt(&p32*&p32+&p33*&p33));

    %let Cp = 1004.67;              /* heat capacity of air, J kg^-1 K^-1;*/
    %let Lv = 2440;                 /* latent heat of vaporization, J g^-1;*/
    %let Rd = 0.28704;              /* gas constant for dry air, J g^-1 K^-1;*/
    %let Rv = 0.0004615;            /* gas constant of water vapor, 0.0004615 J mg-1 K-1;*/
    %let R = 8.31447e-3;    		/* universal gas constant of dry air, (kPa m^3) / (K mol);*/

	data raw ; 
	infile "&path.&dsn."  dsd  firstobs=5;
	input dhms :$22. u v w co2 h2o ts_ck press diag_csat;  /* UNITS: u,v,w=m s-1, co2=mg m-3, h2o=g m-3, ts=șC, pressure=kPa;*/
	sasdate=input(scan(dhms,1," "),yymmdd10.); sastime=input(scan(dhms,2,' '),time12.);
	sasdt=dhms(sasdate,0,0,sastime); 
	if minute(sasdt) < 30 then halfhr=hour(sasdt)+0.5; else if (30<=minute(sasdt)< 60) then halfhr=hour(sasdt)+1; format halfhr 4.1;
/* define 5-min intervals for stationarity test */
	do j=5 to 60 by 5;
		if minute(sasdt) <j and minute(sasdt) ge (j-5) then fiveminute=j;
	end; drop j;
		if abs(u)>32.768 then u=.;
		if abs(v)>32.768 then v=.;
		if abs(w)>8.192 then w=.;
		if ts_ck<-30 or ts_ck>50 then ts_ck=.;
		if co2<0 or co2>5148 then co2=.;
		if h2o<0 or h2o>42 then h2o=.;
	drop dhms sastime sasdt;
	run;


%despike;


/********************************************************************************************;
* Adjust wind coordinates by planar fit														;
********************************************************************************************;*/
%if &rotation.=1 %then %do;
	data filter; merge filter (keep=sasdate halfhr fiveminute u v w co2 h2o ts_ck press diag_csat) countspikes; 
		by sasdate halfhr ;
	u1=u*&cosa.-(w-&b0.)*&sina.;
	v1=u*&sina.*&sinb.+v*&cosb.+(w-&b0.)*&sinb.*&cosa.;
	w1=u*&sina.*&cosb.-v*&sinb.+(w-&b0.)*&cosb.*&cosa.;
	gamma=atan(v1/u1);
	u2=u1*cos(gamma)+v1*sin(gamma);
	v2=-u1*sin(gamma)+v1*cos(gamma);
	w2=w1;
	drop u v w gamma u1 v1 w1;
	lag1co2=lag(co2); lag2co2=lag2(co2);lag3co2=lag3(co2); lag4co2=lag4(co2);lag5co2=lag5(co2);
	lag1h2o=lag(h2o); lag2h2o=lag2(h2o);lag3h2o=lag3(h2o); lag4h2o=lag4(h2o);lag5h2o=lag5(h2o);
	run;
	data filter; merge filter 
	filter (keep=co2 h2o firstobs=2 rename=(co2=co2lead1 h2o=h2olead1)) 
	filter (keep=co2 h2o firstobs=3 rename=(co2=co2lead2 h2o=h2olead2)) 
	filter (keep=co2 h2o firstobs=4 rename=(co2=co2lead3 h2o=h2olead3)) 
	filter (keep=co2 h2o firstobs=5 rename=(co2=co2lead4 h2o=h2olead4)) 
	filter (keep=co2 h2o firstobs=6 rename=(co2=co2lead5 h2o=h2olead5)) ;
	run;
%end;


/********************************************************************************************;
* Adjust wind coordinates for 2-D rotation													;
********************************************************************************************;*/
%else %if &rotation.=2 %then %do;
		data filter; merge filter (keep=sasdate halfhr fiveminute u30 v30 u v w co2 h2o ts_ck press diag_csat) countspikes; 
			by sasdate halfhr ;
				theta=atan(v30/u30);
				u1=u*cos(theta)+v*sin(theta);
				v1=-u*sin(theta)+v*cos(theta);
				w1=w;
	        run;
			proc means data=filter noprint; by sasdate halfhr;
				var  u1 v1 w1 ; 
				output out=rot2d 
					mean= u1m v1m w1m ; 
			run;
			data filter; merge filter rot2d; by sasdate halfhr;
				psi=atan(w1m/u1m);
				u2=u1*cos(psi)+w1*sin(psi);
				v2=v1;
				w2=-u1*sin(psi)+w1*cos(psi);
				lag1co2=lag(co2); lag2co2=lag2(co2);lag3co2=lag3(co2); lag4co2=lag4(co2);lag5co2=lag5(co2);
				lag1h2o=lag(h2o); lag2h2o=lag2(h2o);lag3h2o=lag3(h2o); lag4h2o=lag4(h2o);lag5h2o=lag5(h2o);
			run;

			data filter; merge filter 
				filter (keep=co2 h2o firstobs=2 rename=(co2=co2lead1 h2o=h2olead1)) 
				filter (keep=co2 h2o firstobs=3 rename=(co2=co2lead2 h2o=h2olead2)) 
				filter (keep=co2 h2o firstobs=4 rename=(co2=co2lead3 h2o=h2olead3)) 
				filter (keep=co2 h2o firstobs=5 rename=(co2=co2lead4 h2o=h2olead4)) 
				filter (keep=co2 h2o firstobs=6 rename=(co2=co2lead5 h2o=h2olead5)) ;
			run;
%end;


/********************************************************************************************;
* Maximize covariances 																		;
********************************************************************************************;*/
	proc corr data=filter noprint cov outp=covs30min;
	where halfhr ne . and diag_csat=0;
	by sasdate halfhr; 
		var w2 u2 v2 ts_ck
			co2lead5 co2lead4 co2lead3 co2lead2 co2lead1 co2 lag1co2 lag2co2 lag3co2 lag4co2 lag5co2 
			h2olead5 h2olead4 h2olead3 h2olead2 h2olead1 h2o lag1h2o lag2h2o lag3h2o lag4h2o lag5h2o; 
	run;
		data cov30_vert; set covs30min (keep=sasdate halfhr _NAME_ w2 _TYPE_); where _TYPE_="COV"; drop _TYPE_; run;
		data cov30_lags; set cov30_vert;
			by sasdate halfhr;
			covcomp=abs(w2); 
		if _NAME_="ts_ck"  then classvar="ts_ck";
		else if _NAME_="w2" then classvar="w2";
		else if _NAME_="u2"  then classvar="u2";
		else if _NAME_="v2"  then classvar="v2";
		else if _NAME_ in ("co2lead5","co2lead4","co2lead3","co2lead2","co2lead1","co2","lagco2","lag2co2","lag3co2","lag4co2","lag5co2")
				then classvar="co2";
		else if _NAME_ in ("h2olead5","h2olead4","h2olead3","h2olead2","h2olead1","h2o","lagh2o","lag2h2o","lag3h2o","lag4h2o","lag5h2o")
				then classvar="h2o";
			if _NAME_ in ("co2lead5","h2olead5") then lag=5;
			else if _NAME_ in ("co2lead4","h2olead4") then lag=4;
			else if _NAME_ in ("co2lead3","h2olead3") then lag=3;
			else if _NAME_ in ("co2lead2","h2olead2") then lag=2;
			else if _NAME_ in ("co2lead1","h2olead1") then lag=1;
			else if _NAME_ in ("co2","h2o") then lag=0;
			else if _NAME_ in ("lagco2","lagh2o") then lag=-1;
			else if _NAME_ in ("lag2co2","lag2h2o") then lag=-2;
			else if _NAME_ in ("lag3co2","lag3h2o") then lag=-3;
			else if _NAME_ in ("lag4co2","lag4h2o") then lag=-4;
			else if _NAME_ in ("lag5co2","lag5h2o") then lag=-5;
		run;
		proc sort data=cov30_lags; by sasdate halfhr classvar; run;
		proc means data=cov30_lags noprint; by sasdate halfhr classvar;
			var covcomp;
			output out=cov30max max=covmax;
		run;
		data cov30parm; merge cov30_lags cov30max; by sasdate halfhr classvar;
		if classvar="co2" then do;
			if abs(w2)=covmax then do; cov_w2co2=w2; lagco2=lag; end;
		end;
		if classvar="h2o" then do;
			if abs(w2)=covmax then do; cov_w2h2o=w2; lagh2o=lag; end;
		end;
		if classvar="u2" then cov_w2u2=w2; 
		if classvar="v2" then cov_w2v2=w2; 
		if classvar="w2" then cov_w2w2=w2; 
		if classvar="ts_ck" then cov_w2ts=w2; 
		run;
		data cov_co2; set cov30parm (keep=sasdate halfhr cov_w2co2 lagco2); by sasdate halfhr; where cov_w2co2 ne .; run;
		data cov_h2o; set cov30parm (keep=sasdate halfhr cov_w2h2o lagh2o); by sasdate halfhr; where cov_w2h2o ne .; run;
		data cov_u2; set cov30parm (keep=sasdate halfhr cov_w2u2); by sasdate halfhr; where cov_w2u2 ne .; run;
		data cov_v2; set cov30parm (keep=sasdate halfhr cov_w2v2); by sasdate halfhr; where cov_w2v2 ne .; run;
		data cov_w2; set cov30parm (keep=sasdate halfhr cov_w2w2); by sasdate halfhr; where cov_w2w2 ne .; run;
		data cov_ts; set cov30parm (keep=sasdate halfhr cov_w2ts); by sasdate halfhr; where cov_w2ts ne .; run;
		data cov30lag; merge cov_co2 cov_h2o cov_u2 cov_v2 cov_w2 cov_ts; by sasdate halfhr; run;


	data filter; merge filter cov30lag; by sasdate halfhr;
		if lagco2=-5 then co2=lag5co2;
		else if lagco2=-4 then co2=lag4co2;
		else if lagco2=-3 then co2=lag3co2;
		else if lagco2=-2 then co2=lag2co2;
		else if lagco2=-1 then co2=lag1co2;
		else if lagco2=1 then co2=co2lead1;
		else if lagco2=2 then co2=co2lead2;
		else if lagco2=3 then co2=co2lead3;
		else if lagco2=4 then co2=co2lead4;
		else if lagco2=5 then co2=co2lead5;

		if lagh2o=-5 then h2o=lag5h2o;
		else if lagh2o=-4 then h2o=lag4h2o;
		else if lagh2o=-3 then h2o=lag3h2o;
		else if lagh2o=-2 then h2o=lag2h2o;
		else if lagh2o=-1 then h2o=lag1h2o;
		else if lagh2o=1 then h2o=h2olead1;
		else if lagh2o=2 then h2o=h2olead2;
		else if lagh2o=3 then h2o=h2olead3;
		else if lagh2o=4 then h2o=h2olead4;
		else if lagh2o=5 then h2o=h2olead5;
	drop lagco2 lagh2o 
			lag5co2 lag4co2 lag3co2 lag2co2 lag1co2 co2lead1 co2lead2 co2lead3 co2lead4 co2lead5
			lag5h2o lag4h2o lag3h2o lag2h2o lag1h2o h2olead1 h2olead2 h2olead3 h2olead4 h2olead5;
	run;


/********************************************************************************************;
* Calculate covariances for 30-min flux														;
********************************************************************************************;*/
	proc corr data=filter noprint cov outp=covs30min;
	where halfhr ne . and diag_csat=0;
	by sasdate halfhr; 
	var u2 v2 w2 co2 h2o ts_ck press; 
	run;
	data mean30; set covs30min; by sasdate halfhr; where _TYPE_="MEAN"; 
		rename u2=u2m v2=v2m w2=w2m co2=co2m h2o=h2om ts_ck=ts_ckm press=pressm; run;
	data std30; set covs30min (rename=(w2=sigmawrot u2=sigmaurot v2=sigmavrot)); by sasdate halfhr; where _TYPE_="STD"; 
		keep sasdate halfhr sigmawrot sigmaurot sigmavrot; run;
	data cov30; set covs30min; by sasdate halfhr; where _TYPE_="COV" and _NAME_="w2"; 
		rename u2=cov_w2u2 v2=cov_w2v2 w2=cov_w2w2 co2=cov_w2co2 h2o=cov_w2h2o ts_ck=cov_w2ts; drop press; run;
	data n30; set covs30min (rename=(u2=n1)); by sasdate halfhr; where _TYPE_="N"; keep sasdate halfhr n1; run;

/********************************************************************************************;
* Calculate covariances for steady-state test												;
********************************************************************************************;*/
	proc corr data=filter noprint cov outp=covs5min;
	where halfhr ne . and diag_csat=0;
	by sasdate halfhr fiveminute; 
	var u2 v2 w2 co2 h2o ts_ck; 
	run;
	data cov5; set covs5min; by sasdate halfhr fiveminute; where _TYPE_="COV" and _NAME_="w2"; 
		rename u2=cov_w2u25 v2=cov_w2v25 w2=cov_w2w25 co2=cov_w2co25 h2o=cov_w2h2o5 ts_ck=cov_w2ts5 ; run;

		data statnry; merge cov5 cov30; by sasdate halfhr;
		statnry_co2w=abs((cov_w2co25-cov_w2co2)/cov_w2co2);
		statnry_h2ow=abs((cov_w2h2o5-cov_w2h2o)/cov_w2h2o);
		statnry_tsw=abs((cov_w2ts5-cov_w2ts)/cov_w2ts);
		statnry_uw=abs((cov_w2u25-cov_w2u2)/cov_w2u2);
		statnry_vw=abs((cov_w2v25-cov_w2v2)/cov_w2v2);
		run;
		proc means data=statnry noprint; by sasdate halfhr;
		var statnry_co2w statnry_h2ow statnry_tsw statnry_uw statnry_vw;
		output out=statnryout mean=statnry_co2w statnry_h2ow statnry_tsw statnry_uw statnry_vw;
		run;
/********************************************************************************************;
* steady state test over																	;
********************************************************************************************;*/


	data in_4; merge mean30 cov30 std30 n30 statnryout (drop=_TYPE_) countspikes cov30lag; by sasdate halfhr; 
	where sasdate ne . and halfhr ne .;
	hour=floor(halfhr);
	if halfhr-floor(halfhr)=0.5 then minute=30; 
	else if halfhr-floor(halfhr)=0 then minute=0; 
	sasdt=dhms(sasdate,hour,minute,0); format sasdt datetime.;
	sasdatestring=catx(" ",put(sasdate,yymmdd10.),put(hms(hour,minute,0),time.)); 

		h2o_gm3=h2om*18.015/1000; 							/* + convert form "mmol m-3" to "g m-3");*/
	    co2_mgm3=co2m*44;	 								/* convert from "mmol m-3" to "mg m-3";*/
		airconc=pressm/(&R.*(ts_ckm+273.16));				/* ++ molar concentration of air, mol m^-3;*/
		e = (h2om/1000)*&R.*(ts_ckm+273.16);				/* ++ pressure of water vapor, kPa;*/
			chi=e/pressm;									/* chi per WET AIR, not much difference from  ;*/
			Tk_k=(ts_ckm+273.16)/(1+0.32*chi);				/* ++ Schotanus correction, according to Leuning ;*/
			Tk_c=Tk_k-273.16;
		rho_d = (pressm-e)/(Tk_k*&Rd.);						/* as per CSI program, matches, but can't figure out the units ;*/
		rho_a = rho_d+(h2om/1000);							/* origial, density of wet air, kg m^-3;*/
		mixratio_h2o=(&R*Tk_k*h2om*1000)/(pressm*18.01); 	/* ++ mmol mol-1, mix ratio = (R*T*concentr. (mjugram per m^3) /(press*Molarmass); */
		mixratio_co2=(&R*Tk_k*co2m*1000)/(pressm*44); 		/* ++ umol mol-1 ;*/

	tau1 = sqrt((cov_w2u2*cov_w2u2)+(cov_w2v2*cov_w2v2)); 
	ustar = sqrt(tau1); 
	tau = rho_a*tau1; 
	
	Hs=cov_w2ts*rho_a*&Cp; 
	LE=(1+chi)*(cov_w2h2o+(h2om/Tk_k)*cov_w2ts)*&Lv.; 
	Fc=(cov_w2co2+co2m*((LE/(&Lv.*rho_d*1000))+(Hs/(rho_a*&Cp.*Tk_k))))*1000/44; 

	co2WplHnew  =  (co2m*Hs/(rho_a*&Cp.*Tk_k))*1000/44;
	co2WplLEnew = (co2m*LE/(&Lv.*rho_d*1000))*1000/44;
	
	file "&path.&dsn._flux.txt" lrecl=3000 dsd; 
	if _N_=1 then put
	"Time,ustar,tau,ave_CO2,ave_H2O,ave_P,ave_Tc,ave_U,ave_V,ave_W,
		uspike,vspike,wspike,tsspike,co2spike,h2ospike,lagco2,lagh2o,
		Hs,co2wplLE,co2wplH,LE,Fc,Records,Stationarity_CO2W,Stationarity_H2OW,Stationarity_TsW,sigma_u,sigma_v,sigma_w";
	put sasdatestring$ ustar tau mixratio_co2 mixratio_h2o pressm Tk_c u2m v2m w2m
		uspike vspike wspike tsspike co2spike h2ospike lagco2 lagh2o
		Hs co2WplLEnew co2WplHnew LE Fc n1 statnry_co2w statnry_h2ow statnry_tsw sigmaurot sigmavrot sigmawrot;
	run;
	proc datasets nolist; delete raw filter; run;
%mend planar;



/********************************************************************************************;
* Enter files to be processed																;
********************************************************************************************;*/

%planar(ENTER_FILE_NAME_HERE);
%planar(NEXT_FILE_NAME_HERE);
%planar(etc ...);