disp('read_profs')
fclose('all');
clear all;
%****
g=9.82;
kon=0.4;
T0=273.16;
sgb=5.67e-8;
%Find the flux_file data 
%f=['D:\cwf\data\sheba\mettwr_1&2.txt'];
f=['D:\matlabstf\cwf\qpro\mettwr_1&4.txt'];

flist=fopen(f,'r');
flx_lab1=fgets(flist);
flx_lab2=fgets(flist);
datzz=fscanf(flist,'%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g%g',[21,2900]);
fclose('all');

%jazz(1)=jd_son(j);%decimal JD
%jazz(2:3)=ws_i(j,1:2);%wind speed (m/s)at 2m, 10m
%jazz(4:5)=wd_i(j,1:2);%wind direction at 2m, 10 m
%jazz(6:7)=t_v_i(i,1:2);%temp(C) at 2 m, 10 m
%jazz(8:9)=q_v_i(i,1:2);%mixing ratio (g/kg) at 2 m, 10 m
%jazz(10)=us_md(j);%ustar (m/s)
%jazz(11)=hs_md(j);%sensible heat flux (W/m2)
%jazz(12)=T_sn_an(i);%General Eastern 12 cm temp (C)
%jazz(13)=Tsnwr(i);%radiative surface temp (C)
%jazz(14)=Tice(i);%snow/ice interface temp (C)
%jazz(15)=Rtdx(i);%longwave down
%jazz(16)=Rtux(i);%longwave up
%jazz(17)=Rsdx(i);%solar down (as if there was any)
%jazz(18)=Rsux(i);%solar up
%jazz(19)=orgx(i);%optical raingauge, mm/h
%jazz(20)=ncrx(i);%ncar weighing raingauge, mm/h
%vectorized print, 20 columns, 1 rows

% number of records in flux file
nflux=2900;
cnt_mx=2900;
%Find the flux_file data 
%e=['D:\matlabstf\cwf\qpro\prof_1&2.txt'];
%e=['D:\cwf\data\sheba\prof_1to4.txt'];
e=['D:\matlabstf\cwf\qpro\profs_1&4.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);
    flx_lab2=fgets(flist);
  datx=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%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%g%g%g%g%g%g%g%g%g%g%g',[81,2900]);
  %load prof datx;
  [a,cnt_mx]=size(datx);
     
          n=1;
          m=0;
     
     		for j=1:cnt_mx;   
         		if datx(77,j)==0 & datx(78,j)==0 & datx(79,j)==0 & datx(80,j)==0 & datx(81,j)==0 ,
                  if datx(8,j)<20 & datx(9,j)<20 & datx(10,j)<20 & datx(11,j)<20 ,
                     if datx(13,j)<110 | datx(13,j)>260,                                                              
                        sgum=median([datx(37,j) datx(38,j) datx(39,j) datx(40,j) datx(41,j)]);
                        sgwm=median([datx(47,j) datx(48,j) datx(49,j) datx(50,j) datx(51,j)]);
								if abs(datx(37,j)-sgum)<.5 & abs(datx(38,j)-sgum)<.5 & abs(datx(39,j)-sgum)<.5 & abs(datx(40,j)-sgum)<.5 & abs(datx(41,j)-sgum)<.5,		
									if abs(datx(47,j)-sgwm)<.1 & abs(datx(48,j)-sgwm)<.1 & abs(datx(49,j)-sgwm)<.1 & abs(datx(50,j)-sgwm)<.1 & abs(datx(51,j)-sgwm)<.1,		
                              data(1:81,n)=datx(1:81,j);
                              datz(1:21,n)=datzz(1:21,j);
                               n=n+1;
                           end;    
                        end;
                      end;
               	end;
          		end;
         end;
   end;
      %    while feof(flist)==0,%go until you run out of 1-hr records
      cnt_my=cnt_mx;
      cnt_mx=n-1;

jd2=datz(2,:);
Tsb=datz(14,:);%barnes ts
Tsi=datz(15,:);%snow-ice interface t
Rld=datz(16,:);
Rlu=datz(17,:);
Tse=((Rlu-.01*Rld)/.99/sgb).^0.25-273.16;
Ts=(Tsb+Tse)/2;
Th10=datz(8,:)+0.1;
U10=datz(4,:)+0.1;


jd(1:cnt_mx)=data(1,1:cnt_mx);
z(1:5,1:cnt_mx)=data(2:6,1:cnt_mx);
ws(1:5,1:cnt_mx)=data(7:11,1:cnt_mx);
wd(1:5,1:cnt_mx)=data(12:16,1:cnt_mx);
Th(1:5,1:cnt_mx)=data(17:21,1:cnt_mx);
Q(1:5,1:cnt_mx)=data(22:26,1:cnt_mx);
us(1:5,1:cnt_mx)=data(27:31,1:cnt_mx);
hs(1:5,1:cnt_mx)=data(32:36,1:cnt_mx);
     
sigu(1:5,1:cnt_mx)=data(37:41,1:cnt_mx);
sigv(1:5,1:cnt_mx)=data(42:46,1:cnt_mx);
sigw(1:5,1:cnt_mx)=data(47:51,1:cnt_mx);
sigT(1:5,1:cnt_mx)=data(52:56,1:cnt_mx);

Cu(1:5,1:cnt_mx)=data(57:61,1:cnt_mx);
Cv(1:5,1:cnt_mx)=data(62:66,1:cnt_mx);
Cw(1:5,1:cnt_mx)=data(67:71,1:cnt_mx);
CT(1:5,1:cnt_mx)=data(72:76,1:cnt_mx);
          
fl(1:5,1:cnt_mx)=data(77:81,1:cnt_mx);

wss=mean(ws);
usm=mean(us);
usn=median(us);

sgum=mean(sigu);
sgvm=mean(sigv);
sgwm=mean(sigw);
Cum=median(z.^.667.*Cu);
Cwm=median(z.^.667.*Cw);
hsm=mean(hs);
Thm=mean(Th);
%
usx=max(.02,(usn));
Linv=-5*kon*g./(Thm+T0).*hsm/1300./usx.^3;
n=0;
for i=1:cnt_mx,
   if U10(i)>3,
      n=n+1;
      usy(n)=usx(i);
      for j=1:5,
         zet(j,n)=z(j,i)/5*Linv(i);
         psimx(j,n)=psiu(zet(j,n));
      end;
      wsn(1:5,n)=ws(1:5,i)+usx(i)/kon.*psimx(1:5,n);
      wsx(1:5,n)=ws(1:5,i);
      zs(1:5,n)=z(1:5,i);
   end;
end;
