disp('hr_bulk');
n=length(jda);
uss=zeros(n);
clear yflx yhsb yhlb ytaub yLb yusb ytsb yqsb ydter yhrain
for i=1:n
   xx=[v1(i) 0  yts(i) y1(i,19) yqs(i)  y1(i,21)  y1(i,22) y1(i,23) y1(i,24) zi y1(i,157) zu zt zq y1(i,151) 1 0 0 0];
% %    u       us    ts    ta      qs         qa       rs        rl     rain        press           lat
% 
% %    1          2     3     4      5           6        7         8       9     10 11 12 13  14    15    16
% yflx=cor30a(xx);
   A=coare35vn(v1(i),zu,y1(i,19),zt,rhbar(i),zq,y1(i,157),yts(i),y1(i,22),y1(i,23),y1(i,151),zi,y1(i,24),NaN,NaN);
   yflx=[A(3) A(4) A(2) A(42) A(10) A(11) A(15) A(1) A(8) A(9)/1000 A(17) A(18) A(19) A(34) A(41) A(12) A(13) A(14) A(31)/1000 A(32)/1000 A(33)/1000 A(40) ]; %matches old C3 ouput format.

yhsb(i)=yflx(1);
yhlb(i)=yflx(2);
ytaub(i)=yflx(3);
yLb(i)=yflx(7);
yusb(i)=yflx(8);
ytsb(i)=yflx(9);
yqsb(i)=yflx(10);
ydter(i)=yflx(11);
yhrain(i)=yflx(14);
end;


yta=y1(:,19)';
yqa=y1(:,21);
yU=y1(:,16)';
yP=y1(:,157);  %add pressure from sensor
yrhoair=yP.*100./(Rgas*(yta'+tdk).*(1+0.61*yqa/1000));

ywtvb=-yusb.*(ytsb+.61.*(yta+273.16).*yqsb);
ywt_sonb=-yusb.*(ytsb+.51.*(yta+273.16).*yqsb);
yws=(9.83./(yta+273.16).*abs(ywtvb)*600).^.333;%convective velocity
ygg=sqrt(1+1.44*yws.^2./yU.^2);%hrly average gust factor
yhsbs=yrhoair'*cpa.*ywt_sonb;

%**************   recompute ID fluxes  *************
zuu=zu;%-.7;
zety=zuu./yLb;zety=zety';
ywt_son=y1(:,121)';
ywt=ywt_son+.51*(yta+tdk).*yusb.*yqsb;
% ycu=y1(:,134);
% ycw=y1(:,136);
ycu=y1(:,135);
ycw=y1(:,137);
ycs=(ycu+.75*ycw)/2;%aver cu and cw for smoother structure function
ii=find(isnan(ycu));
ycs(ii)=.75*ycw(ii);
kk=find(isnan(ycw));
ycs(kk)=ycu(kk);
% yct=y1(:,139);
yct=y1(:,138);
ycql=ycqfl';%Licor

 a=psi_fu(zety);ii=find(isnan(a)); %correction of the complex numbers into psi_fu(zety), change them by NaN
 a(ii)=median(a);

yusib=sqrt(ycs.*zuu^.667./a);yusib=yusib';
ytsibs=-sqrt(yct.*zuu^.667./psi_ft(zety)).*sign(ywt_sonb');ytsibs=ytsibs';

yqsib_lic=sign(yqsb').*sqrt(ycql.*zuu^.667./psi_ft(zety))*1e-3;yqsib_lic=yqsib_lic';%convert to g/g

% if funit=='li1'
%     ii=find(yqsib_lic<-3e-5 | yqsib_lic>0);yqsib_lic(ii)=NaN;
% elseif funit=='li2'
%     ii=find(yqsib_lic<-3e-4 | yqsib_lic>0);yqsib_lic(ii)=NaN;    
% elseif funit=='lic'
%     ii=find(yqsib_lic<-3e-4 | yqsib_lic>1e-3);yqsib_lic(ii)=NaN;
% end;
ii=find(yqsib_lic<-3e-4 | yqsib_lic>0);yqsib_lic(ii)=NaN;    

ytsib=ytsibs-.51*(yta+273.16).*yqsb;

%************   

yhsib=-yrhoair'*cpa.*yusib.*ytsib;
yhsibs=-yrhoair'*cpa.*yusib.*ytsibs;
yhlib_lic=-yrhoair'*Le.*yusib.*yqsib_lic;
ytauib=yrhoair'.*yusib.*yusib./ygg;

%**************
yhsc=yrhoair'*cpa.*ywt;
yhscs=yrhoair'*cpa.*ywt_son;
yhlc=yrhoair'*Le.*y3(:,9)'/1000;  %Licor
ytauc=-yrhoair'.*y1(:,117)';
ytaucx=ytauc;
ytaucy=-yrhoair'.*y1(:,119)';

%*****************   create blended turbulence latent flux
figure;
subplot(3,1,1);plot(jda,yhlc);title('yhlc');ylim([0 400]);grid
subplot(3,1,2);plot(jda,yhlib_lic);title('yhlib_ lic');ylim([0 400]);grid
subplot(3,1,3);plot(jda,hlc_licbar);title('hlc_ licbar');ylim([0 400]);grid
hll=(yhlc+yhlib_lic+hlc_licbar)/3;
ii=find(isnan(yhlc));%case bad hlc
jj=find(isfinite(yhlib_lic(ii)));
	kk=find(isnan(hlc_licbar(ii(jj))));
	if ~isempty(kk)
   	hll(ii(jj(kk)))=yhlib_lic(ii(jj(kk)));%only hlib is good
   end;
   kk=find(isfinite(hlc_licbar(ii(jj))));
   if ~isempty(kk)
      hll(ii(jj(kk)))=(yhlib_lic(ii(jj(kk)))+hlc_licbar(ii(jj(kk))))*.5;%both hlib and hl_lic good
   end;
jj=find(isfinite(hlc_licbar(ii)) & isnan(yhlib_lic(ii)));
if ~isempty(jj)
   hll(ii(jj))=hlc_licbar(ii(jj));%only hl_lic good
end;

ii=find(isfinite(yhlc));%case good hlc
jj=find(isfinite(yhlib_lic(ii)));%good hlib
	kk=find(isnan(hlc_licbar(ii(jj))));%bad hl_lic
	if ~isempty(kk)
   	hll(ii(jj(kk)))=(yhlc(ii(jj(kk)))+yhlib_lic(ii(jj(kk))))*.5;%hlc and  hlib good
   end;
jj=find(isnan(yhlib_lic(ii)));%bad hlib
	kk=find(isfinite(hlc_licbar(ii(jj))));%good hl_lic
	if ~isempty(kk)
   	hll(ii(jj(kk)))=(yhlc(ii(jj(kk)))+hlc_licbar(ii(jj(kk))))*.5;%hlc and  hl_lic good
   end;
% hll=hlc_licbar;
%********************  god its a pain in the ass  *********

%create blended licor flux

yhlc_bl=yhlc;
ii=find(isnan(yhlc));
yhlc_bl(ii)=hlc_licbar(ii);
ii=find(isfinite(yhlc) & isfinite(hlc_licbar));
if ~isempty(ii);yhlc_bl(ii)=(yhlc(ii)+hlc_licbar(ii))/2;end;


%

figure;plot(jda,(yhsib+yhsc)*.5,'o',jda,yhsb,[stdt endt],[0 0]);axis([stdt endt -100 300]);
xlabel('Day-Of-year (2013)');ylabel('Sensible Heat Flux (W/m^2)');title([jdollars ': line - bulk,  o - turbulence'])

figure;plot(jda,hll+hl_webar','o',jda,yhlb,[stdt endt],[0 0]);axis([stdt endt -20 300]);
xlabel('Day-Of-year (2013)');ylabel('Latent Heat Flux (W/m^2)');title([jdollars ': line - bulk,  o - turbulence'])

figure;plot(jda,(ytauib+ytauc)*.5,'o',jda,ytaub,[stdt endt],[0 0]);axis([stdt endt -.05 2]);
xlabel('Day-Of-year (2013)');ylabel('Stress (N/m^2)');title([jdollars ': line - bulk,  o - turbulence'])

figure;plot(ytaub,(ytauib+ytauc)*.5,'o',[-.05 2],[-.05 2],[-.05 2],[0 0],[0 0],[-.05 .2]);axis([-.05 2 -.05 2]);
xlabel('Bulk Stress (N/m^2)');ylabel('Turbulent Stress (N/m^2)');title(jdollars)

figure;plot(yhlb,hll+hl_webar','o',[-20 300],[-20 300]);axis([-20 300 -20 300]);
xlabel('Bulk Latent Heat Flux (W/m^2)');ylabel('Turbulent Latent Heat Flux (W/m^2)');title(jdollars)
