disp('make_data_file_1')
fclose('all');
clear all;
%****
% number of records in flux file
nflux=5*24*24;
% number of records in stat file
nstat=23*24;  
cnt_mx=24*24;
plotit=0;
g=0.0098;%adiabatic lapse rate
%Find the flux_file data 
e=['D:\cwf\data\sheba\proc_son_data.txt'];
	flist=fopen(e,'rt');
   if flist>0,   %if the file exists, 
      %read flux data file, find time
      %check whether data is good or not
      % -use good/bad flag
      % -other methods???
      %compute one value of u*,Hs - use median
      %compute u2m,T10m,u10m,T10m - sim theory
      flx_lab1=fgets(flist);
   
      loop=1;  %loop thru 60-min sections 
      n=0;
%    while feof(flist)==0,%go until you run out of 1-hr records
     for icnt=1:cnt_mx;
        n=n+1;
        data=fscanf(flist,'%g%d%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g',[19,5]);
for j=1:5;%loop through 5 sonics
          jd_son(n)=data(1,j);
 	 		 z(n,j)=data(3,j);
          us(n,j)=data(8,j);
          hs(n,j)=data(9,j);
          ws(n,j)=data(4,j);
          wd(n,j)=data(5,j);
          sigu(n,j)=data(11,j);
          sigv(n,j)=data(12,j);
          t(n,j)=data(6,j);
          fl(n,j)=data(19,j);
%test for reasonable values
            if(ws(n,j) > 50);
              fl(n,j)=1;
            end;
            if(t(n,j)<-80 |t(n,j)>30);
              fl(n,j)=1;
            end;
            if(us(n,j)>20);
              fl(n,j)=1;
            end;
            if(abs(hs(n,j)>200));
              fl(n,j)=1;
            end;
         end; % end sonic reading loop
%test for vertical continuity
% not done now

%calculate median u* and hs for each time

%test for flag trip
         ngd=0;          
         for j=1:5,
          if(fl(n,j)<1);
            ngd=ngd+1;
            tst_u(ngd)=us(n,j);
            tst_hs(ngd)=hs(n,j);
            tst_w(ngd)=ws(n,j);
            tst_wd(ngd)=wd(n,j);
				tst_z(ngd)=z(n,j);   
         end;
         end;
         if (ngd>0);
           us_md(n)=median(tst_u);
           hs_md(n)=median(tst_hs);
           ws_md(n)=median(tst_w);
           wd_md(n)=median(tst_wd);
           z_md(n)=median(tst_z);

         else
           us_md(n)=99.99;
           hs_md(n)=99.99;
           ws_md(n)=99.99;
           z_md(n)=99.99;
           wd_md(n)=999;
         end;
         clear tst_u tst_hs tst_w tst_z tst_wd;
%calculate 2-m and 10-m ws,wd and T through interpolation
         ht=[2.5 10];
         ngd=0;
         for j=1:5;
           if fl(n,j) <1;
             ngd=ngd+1;
             ind(ngd)=j;
             zgd(ngd)=z(n,j);
           end;%end if
        	end;%end j
        	for k=1:2;
         	for j=1:ngd-1;
           		if zgd(j)<=ht(k)&zgd(j+1)>=ht(k),
            		ws_i(n,k)=ws(n,ind(j))+(ws(n,ind(j+1))-ws(n,ind(j)))*(ht(k)-z(n,ind(j)))/(z(n,ind(j+1))-z(n,ind(j)));
             		wd_i(n,k)=wd(n,ind(j))+(wd(n,ind(j+1))-wd(n,ind(j)))*(ht(k)-z(n,ind(j)))/(z(n,ind(j+1))-z(n,ind(j)));
             		t_son_i(n,k)=t(n,ind(j))+(t(n,ind(j+1))-t(n,ind(j)))*(ht(k)-z(n,ind(j)))/(z(n,ind(j+1))-z(n,ind(j)));
             		break;
          		else;
            	 	ws_i(n,k)=99.9;
            	 	wd_i(n,k)=999;
            	 	t_son_i(n,k)=99.9;
           		end;%end if
           	end;%end j
          end;%end k  end 2,10-m loop
          if ngd>1,
             if ws_i(n,1)==99.9,
            	if fl(n,2)<1,
                  ws_i(n,1)=ws(n,2);
                  wd_i(n,1)=wd(n,2);
            	end;
         	end;
          	if ws_i(n,2)==99.9,
            	if ws_md(n)<99.9,
                  ws_i(n,2)=ws_md(n)*log(10/0.0004)/log(z_md(n)/0.0004);
                  wd_i(n,2)=wd_md(n);
            	end;
         	end;
         else
            ws_i(n,1:2)=99.99;
            wd_i(n,1:2)=999;
         end;
      end; %end reading/processing sonic data file
n_s=n;
end; %end flist 
      fclose (flist)

%now read in stat file
	
%e='D:\cwf\data\sheba\proc_stat.txt'
	e='D:\matlabstf\cwf\sheba_proc\proc_stat.txt'
flist=fopen(e,'rt');
   if flist>0,   %if the file exists, 
      %read flux data file, find time
      %check whether data is good or not
      % -use good/bad flag
      % -other methods???
      %compute one value of u*,Hs - use median
      %compute u2m,T10m,u10m,T10m - sim theory
      slw_lab=fgets(flist);
      n=0;
%  while feof(flist)==0,%go until you run out of 1-hr records
   for icntv=1:cnt_mx;
      n=n+1;
jaz=fscanf(flist,'%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g',35);
jd_v(n)=jaz(1);
 zv(1:5,n)=jaz(2:6);
 Tice(n)=jaz(17);
 Tsnwr(n)=jaz(18);
 T_sn_an(n)=jaz(19);
 Th(1:5,n)=jaz(20:24);
%convert back to temperature
Th(1:5,n)=Th(1:5,n)-g*zv(1:5,n);
T_sn_an(n)=T_sn_an(n)-.12*g;
 Qx(1:5,n)=jaz(25:29);
 Rtdx(n)=jaz(30);
 Rtux(n)=jaz(31);
 Rsdx(n)=jaz(32);
 Rsux(n)=jaz(33);
 orgx(n)=jaz(34);
 ncrx(n)=jaz(35)*60;
%interpolate temp&q to 2m and 10m
         ht=[2.5 10];
         ngd=0;
         clear  ind
         for j=1:5;
           if Th(j,n) <99;
             ngd=ngd+1;
             ind(ngd)=j;
             zge(ngd)=zv(j,n);
           end;
         end;
         for k=1:2;
         for j=1:ngd-1;
           if zge(j)<=ht(k)&zge(j+1)>=ht(k),
             t_v_i(n,k)=Th(ind(j),n)+(Th(ind(j+1),n)-Th(ind(j),n))*(ht(k)-zv(ind(j),n))/(zv(ind(j+1),n)-zv(ind(j),n));
             q_v_i(n,k)=Qx(ind(j),n)+(Qx(ind(j+1),n)-Qx(ind(j),n))*(ht(k)-zv(ind(j),n))/(zv(ind(j+1),n)-zv(ind(j),n));
          	  break;	   
          else;
             t_v_i(n,k)=99.9;
             q_v_i(n,k)=99.9;
           end;
         end;
         end;%end 2,10-m loop
         if ngd>1,
            if t_v_i(n,1)==99.9,
               if Th(2,n)<99,
                  t_v_i(n,1)=Th(2,n);
                  q_v_i(n,1)=Qx(2,n);
               end;
            end;
         end;   
      end;%end while feof
n_v=n;
 end;%end if flist
fclose(flist)   

%now match times for sonic data and stat data and print to file
nt=min([n_s n_v]);
i=0;
j=0;
for k=1:n_s;
i=i+1;
j=j+1;
if (jd_son(j)==jd_v(i));

prt_mainfile;

elseif (jd_son(j) > jd_v(i));
   i=i+1;
else;
   j=j+1;
end;
end;
   




%1-hr time series

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

Tmed=median(Th);
Tmed=median(Tmed);
pll=Tmed-5;
pul=Tmed+20;
%plot time series
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([-34 -16 .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; radiation=puke grn dots'); 
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),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; 
