disp('read_stats_day_mod')
fclose('all');
clear all;
jd=input('julian day, ddd','s');
plotit=input('plot plots? 0=n, 1=y');
%plotit=0;
prtit=input('append master data file? 0=n, 1=y');
%prtit=1;
%yd=date;
%jif=datevec(yd);
%xx=[jif(1) jif(2) jif(3) 0 0 0];
%jd=num2str(julday(xx)-1);
hr='00';%input('gmt hr','s');
g=0.0098;%adiabatic lapse rate
eps_i=0.99;%emissivity of snow
sig_sb=5.67e-8;%Stefan Boltzmann constant
blk=1;%10-min index 
stux=zeros(7,50);%basic stats file block size, 7 cols 50 lines
%some empirical offsets determined in halfassed manner by CWF
Tcor=[+.05 -.05 -.05 +.05 +.05];
ucor=[+.04 -.1 -.05  +.00 +.3];
vcor=[-.2 .00 .00 +.00 -.3];
%****
Tpl=zeros(1,6);
zlch=zeros(1,6);
for jam=0:23,%cycle thru 24 hourly stats files
   if jam<10,
      hr=['0'num2str(jam)];
   else
      hr=num2str(jam);
   end;
   
%Find the data 
%e=['C:\SHEBA\data\98' jd hr '.sta'];
%e=['F:\users\data\98' jd hr '.sta'];
e=['E:\98' jd hr '.sta'];
%e=['D:\users\data\97' jd hr '.sta'];
%e=['F:\users\pguest\98' jd hr '.sta'];
e=['F:\DasWin32\data\ascii\98'   jd hr '.sta'];

	flist=fopen(e,'r');
   if flist>0,   %if the file exists, 
      %find instrument heights
      for jj=1:25,
      	bum=fgets(flist);
      	if jj==9,
            bam=bum;
            if jd>304
               z(jam+1,1)=str2num(bum(7:10));
            else   
            	z(jam+1,1)=str2num(bum(7:7));
            end;%if jd
         end;%if jj
      end; %end jj
   	z(jam+1,2)=3.19;%fixed heights levels 2-4
		z(jam+1,3)=5.14;
      z(jam+1,4)=8.94;
      if jd>304
      	z(jam+1,5)=str2num(bum(7:11));
      else
         z(jam+1,5)=str2num(bum(7:8));
      end;%if jd   
      fclose(flist);%close and re-open file to read full data
   	flist=fopen(e,'r');
   
      loop=1;  %loop thru 10-min sections to extract and compute 1-hr aves
      %zero all accumulators
      th_ac=zeros(1,5);u_ac=zeros(1,5);v_ac=zeros(1,5);w_ac=zeros(1,5);
      Rlu_ac=0;Rld_ac=0;Rsu_ac=0;Rsd_ac=0;
      Tcu_ac=0;Tcd_ac=0;Tdd_ac=0;Tdu_ac=0;
      Rhv_ac=zeros(1,5);Tice_ac=0;n_ac=0;Vbr_ac=0;org_ac=0;
      tn_ac=0;ncr_ac=0;
      while feof(flist)==0,%go until you run out of 10-min sections
         dmp=fgets(flist);%read header lines and skip blanks
         dmp1=fgets(flist);
   		dmp2=fgets(flist);
   		dmpx=fgets(flist);
			dat(blk,1:6)=datevec([dmpx(48:55) ' ' dmpx(57:64)]);
         if eval(dmpx(57:58))==12,
            dat(blk,4)=12;
         end;
         jad(blk)=str2num(jd)+dat(blk,4)/24+dat(blk,5)/24/60;
         if loop==1,
            jz(jam+1)=jad(blk);
         end;
   		dmp3=fgets(flist);
			dmp4=fgets(flist);
			dmp5=fgets(flist);
			dmp6=fgets(flist);

         stux=fscanf(flist,'%*s%*s%g%g%g%g%g%g%g',[7,50]);%read main data block
         [aa,bb]=size(stux);%check size
         if bb>49, %only compute on full data blocks
            for i=1:5; %five levels of instruments
               uraw(blk,i)=stux(2,1+4*(i-1));%sonics
					vraw(blk,i)=stux(2,2+4*(i-1));
					wraw(blk,i)=stux(2,3+4*(i-1));
					Traw(blk,i)=stux(2,4+4*(i-1));

					sgu(blk,i)=stux(3,1+4*(i-1));
					sgv(blk,i)=stux(3,2+4*(i-1));
					sgw(blk,i)=stux(3,3+4*(i-1));
					sgT(blk,i)=stux(3,4+4*(i-1));
 
					Tvs(blk,i)=stux(2,26+1*(i-1));%Vaisalas
					Rhv(blk,i)=stux(2,31+1*(i-1));
	
					sgTv(blk,i)=stux(3,26+1*(i-1));
					sgRh(blk,i)=stux(3,31+1*(i-1));
					if i<4,
						Tmn(blk,i)=stux(2,23+1*(i-1));
						sgTm(blk,i)=stux(3,23+1*(i-1));
					end;  %end if  
         	end;%end for i
	         Tan(blk)=stux(2,25)+.5;%Andreas General Eastern T measure, empirical offset
   	      Tsnow(blk)=stux(2,23)-.5;%Tsnow, ice, and PIR empirical offsets .5 C
      	   Tice(blk)=stux(2,24)-.5; %Tice at snow ice interface
         	org(blk)=stux(2,39);%STI optical rain gauge, mm/hr
         	ncr(blk)=stux(2,40);%NCAR weighing rain gauge, mm/min
				Rlu(blk)=stux(2,45);
         	Rld(blk)=stux(2,46);
	         Vbrn(blk)=stux(2,47);%CRREL Barnes Tir
				Rsu(blk)=stux(2,48);
				Rsd(blk)=stux(2,49);

				Tcu(blk)=stux(2,41);%-.5;
				Tcd(blk)=stux(2,42);%-.5;
				Tdu(blk)=stux(2,43);%-.5;
				Tdd(blk)=stux(2,44);%-.5;
         
	         sgRlu(blk)=stux(3,45);
				sgRld(blk)=stux(3,46);
				sgRsu(blk)=stux(3,48);
				sgRsd(blk)=stux(3,49);
         
	         sgTcu(blk)=stux(3,41);
				sgTcd(blk)=stux(3,42);
				sgTdu(blk)=stux(3,43);
         	sgTdd(blk)=stux(3,44);
         	%accumulate for 1-hr aves
      	   th_ac=th_ac+Tvs(blk,1:5);
         	u_ac=u_ac+uraw(blk,1:5);
         	v_ac=v_ac+vraw(blk,1:5);
         	w_ac=w_ac+wraw(blk,1:5);
         	Rhv_ac=Rhv_ac+Rhv(blk,1:5);
         	Tice_ac=Tice_ac+Tice(blk);
         	tn_ac=tn_ac+Tan(blk);
         	org_ac=org_ac+org(blk);
         	ncr_ac=ncr_ac+ncr(blk);
        	 	Rld_ac=Rld_ac+Rld(blk);
         	Rlu_ac=Rlu_ac+Rlu(blk);
         	Rsd_ac=Rsd_ac+Rsd(blk);
         	Rsu_ac=Rsu_ac+Rsu(blk);
         	Tcd_ac=Tcd_ac+Tcd(blk);
         	Tcu_ac=Tcu_ac+Tcu(blk);
         	Tdd_ac=Tdd_ac+Tdd(blk);
         	Tdu_ac=Tdu_ac+Tdu(blk);
         	Vbr_ac=Vbr_ac+Vbrn(blk);
         
	         n_ac=n_ac+1;
   	      loop=loop+1;
      	   blk=blk+1;
         	strk=fgets(flist);
         end;%end if bb
      end;%end while feof
      %compute 1-hr aves
      hlx=th_ac./n_ac;
      Th(jam+1,1:5)=hlx+Tcor(1:5);%
      Th_tn(jam+1)=tn_ac/n_ac;%+g*.12
      ulx(jam+1,1:5)=u_ac./n_ac+ucor;
      vlx(jam+1,1:5)=v_ac./n_ac+vcor;
      wlx(jam+1,1:5)=w_ac./n_ac;
      Rhvx(jam+1,1:5)=Rhv_ac/n_ac;
      Rldx(jam+1)=Rld_ac./n_ac;
      Rlux(jam+1)=Rlu_ac./n_ac;
      Rsdx(jam+1)=Rsd_ac./n_ac;
      Rsux(jam+1)=Rsu_ac./n_ac;
      Tcdx(jam+1)=Tcd_ac./n_ac;
      Tcux(jam+1)=Tcu_ac./n_ac;
      Tddx(jam+1)=Tdd_ac./n_ac;
      Tdux(jam+1)=Tdu_ac./n_ac;
      Ticex(jam+1)=Tice_ac./n_ac;
      orgx(jam+1)=org_ac/n_ac;
      ncrx(jam+1)=ncr_ac/n_ac*60;%convert to mm/hr
      Vbrx(jam+1)=Vbr_ac/n_ac;
   end;%end if flist   
end;%end jam
[npz,a]=size(z);
%10-min time series
Rtd=Rld+sig_sb*((Tcd+273.16).^4)+4*sig_sb*((Tcd+273.16).^4-(Tdd+273.16).^4);
Rtu=Rlu+sig_sb*((Tcu+273.16).^4)+4*sig_sb*((Tcu+273.16).^4-(Tdu+273.16).^4);
Tsr=((Rtu-(1-eps_i)*Rtd)/sig_sb/eps_i).^0.25-273.16;%surface T from eppleys
%Barnes calculations
Nr=20.265*Vbrn+10.453;
Tbrn=-96.888+5.5859*Nr-0.11727*Nr.*Nr+1.2436e-3*Nr.*Nr.*Nr;%Barnes surface T in Cent, uncorrected for emissivity/reflection
Qdn=min(Rtu,2*Rtd-Rtu);%empirical estimate of Barnes downard based on eppley downward
Rbu=sig_sb*(Tbrn+273.16).^4;
Tsb=((Rbu-(1-eps_i)*Qdn)/sig_sb/eps_i).^0.25-273.16;
%
spd=uraw.*uraw+vraw.*vraw;
spd=spd.^.5;
np=blk-1;


%1-hr time series
Rtdx=Rldx+sig_sb*((Tcdx+273.16).^4)+4*sig_sb*((Tcdx+273.16).^4-(Tddx+273.16).^4);
Rtux=Rlux+sig_sb*((Tcux+273.16).^4)+4*sig_sb*((Tcux+273.16).^4-(Tdux+273.16).^4);
Tsrx=((Rtux-(1-eps_i)*Rtdx)/sig_sb/eps_i).^0.25-273.16;
esx=6.112*exp(Th.*17.502./(Th+241)).*Rhvx/100;
Qx=0.622*esx./(1020-0.378.*esx)*1000;% water vapor mixing ratio
%rhox=102000./((Th+273.16).*(1+0.61*qx/1000)*287);% air density, kg/m^3
Th=Th+g*z;%potential temperature
Th_tn=Th_tn+g*.12;

%rotate to streamwise
for j=1:npz;
   %for clam=1:5
   um=ulx(j,1:5);
	vm=vlx(j,1:5);
	wm=wlx(j,1:5);
   %rotate to streamwise flow
	s=sqrt(um.^2+vm.^2);
	thet=atan2(vm,um);
	phi=atan2(wm,s);
   
   rot(1,1,1:5)=cos(phi).*cos(thet);
	rot(1,2,1:5)=cos(phi).*sin(thet);
	rot(1,3,1:5)=sin(phi);
	rot(2,1,1:5)=-sin(thet);
	rot(2,2,1:5)=cos(thet);
	rot(2,3,1:5)=0;
	rot(3,1,1:5)=-sin(phi).*cos(thet);
	rot(3,2,1:5)=-sin(phi).*sin(thet);
   rot(3,3,1:5)=cos(phi);
   for clam=1:5
   x=[ulx(j,clam) vlx(j,clam) wlx(j,clam)];
   rotx=rot(1:3,1:3,clam);
   xr=rotx*x';
   ur(j,clam)=xr(1);
   vr(j,clam)=xr(2);
   wr(j,clam)=xr(3);
   end; %end for clam   
end;%end for j
spdh(1:npz,1:5)=ur;  %magnitude of 1-hr mean streamwise wind vector
dirx=atan2(-vlx,ulx)*180/pi;  %wind direction relative to boom point direction (about N)
if (dirx<0) 
  dirx=dirx+360;
end
j_st=floor(jad(1));
j_en=j_st+1;

if prtit==1;
  prt_stats;
end;

if plotit==1;
%  set some plot limits
ull=-10;
vll=-10;
Tll=-25;
hll=70;
uul=10;
vul=10;
Tul=-0;
hul=100;
sul=12;
Rlul=30;
Rlll=-100;
Rtul=300;
Rtll=125;
Rsul=800;
Rsll=-10;

Tmed=median(Th);
Tmed=median(Tmed);
pll=Tmed-5;
pul=Tmed+20;
%plot T profiles
figure;
for jq=1:npz
   Tpl(1)=Th_tn(jq)+.5*jq;
   Tpl(2:6)=Th(jq,1:5)+.5*jq;
   %Tpl(2)=Tpl(2)+.1;
   zlch(1)=.12;
   zlch(2:6)=z(jq,1:5);
   semilogy(Tpl,zlch,'-x');
   hold on;
end;	  
hold off;
title(['T+.5*hr Profiles, Julian Day ' jd]);
xlabel('Potential T(C)');
ylabel('Height (m)');
axis([-20 -0 .1 20]);
%plot u profiles
figure;
for jq=1:npz
   upl(1:5)=spdh(jq,1:5)+.5*jq;
   semilogy(upl,z(jq,1:5)-.4,'-x');
   hold on;
end;	  
hold off;
title(['U+.5*hr Profiles, Julian Day ' jd]);
xlabel('Wind Speed (m/s)');
ylabel('Height (m)');
%axis([0 20 1 20]);

%plot u,v,w,speed, T, RH, rad time series
figure;
plot(jad(1:np),uraw(1:np,1),jad(1:np),uraw(1:np,2),jad(1:np),uraw(1:np,3),jad(1:np),uraw(1:np,4),jad(1:np),uraw(1:np,5));
xlabel('Jday (1997)');
ylabel('u (m/s)');
axis([j_st j_en ull uul]);
Title('Com#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta'); 
figure;
plot(jad(1:np),vraw(1:np,1),jad(1:np),vraw(1:np,2),jad(1:np),vraw(1:np,3),jad(1:np),vraw(1:np,4),jad(1:np),vraw(1:np,5));
xlabel('Jday (1997)');
ylabel('v (m/s)');
axis([j_st j_en vll vul]);
Title('Com#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta'); 
figure;
plot(jad(1:np),wraw(1:np,1),jad(1:np),wraw(1:np,2),jad(1:np),wraw(1:np,3),jad(1:np),wraw(1:np,4),jad(1:np),wraw(1:np,5));
xlabel('Jday (1997)');
ylabel('w (m/s)');
axis([j_st j_en -1 1]);
Title('Com#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta'); 
figure;
plot(jad(1:np),spd(1:np,1),jad(1:np),spd(1:np,2),jad(1:np),spd(1:np,3),jad(1:np),spd(1:np,4),jad(1:np),spd(1:np,5));
xlabel('Jday (1997)');
ylabel('Speed (m/s)');
axis([j_st j_en 0 sul]);
Title('Com#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta'); 
figure;
plot(jad(1:np),Tvs(1:np,1),jad(1:np),Tvs(1:np,2),jad(1:np),Tvs(1:np,3),jad(1:np),Tvs(1:np,4),jad(1:np),Tvs(1:np,5),jad(1:np),Tsr(1:np),'.',jad(1:np),Tsb(1:np),'x');
xlabel('Jday (1997)');
ylabel('T_{air}, T_{srad} (C)');
axis([j_st j_en Tll Tul]);
Title('lev#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta; Erad=grn dots ;Brad=x'); 
figure;
plot(jad(1:np),Rhv(1:np,1),jad(1:np),Rhv(1:np,2),jad(1:np),Rhv(1:np,3),jad(1:np),Rhv(1:np,4),jad(1:np),Rhv(1:np,5));
xlabel('Jday (1997)');
ylabel('RH (%)');
axis([j_st j_en hll hul]);
Title('lev#: 2=blu; 3=grn; 4=orng; 5=cyan; 6=magenta'); 
figure;
plot(jad(1:np),Rlu(1:np),jad(1:np),Rld(1:np));
xlabel('Jday (1997)');
ylabel('R_{Pile} (W/m^2)');
axis([j_st j_en Rlll Rlul]);
Title('Pyrgeometer: upwell=blu; downwell=grn'); 
figure;
plot(jad(1:np),Rtu(1:np),jad(1:np),Rtd(1:np));
xlabel('Jday (1997)');
ylabel('R_{lu,d} (W/m^2)');
axis([j_st j_en Rtll Rtul]);
Title('Pyrgeometer: upwell=blu; downwell=grn'); 
figure;
plot(jad(1:np),Rsu(1:np),jad(1:np),Rsd(1:np));
xlabel('Jday (1997)');
ylabel('R_{su,d} (W/m^2)');
axis([j_st j_en Rsll Rsul]);
Title('Pyranometer: upwell=blu; downwell=grn'); 


figure;
plot(jad(1:np),org(1:np));
xlabel('Jday (1997)');
ylabel('ORG Precipitation Rate (mm/hr)');
axis([j_st j_en 0 5]);
Title('STI Optical Rain Gauge');
end; 
