disp('hr_bulk_AMMA');
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 P zu  zt zq y1(i,151) 1 0 0 0 ];
%    u       us    ts    ta      qs         qa       rs        rl     rain                   lat

%    1       2     3     4      5           6        7         8       9     10 11 12 13  14    15    16
yflx=cor30a(xx);
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)';
yrhoair=1010*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(:,142);
ycw=y1(:,144);
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(:,146);
ycq=ycqf';%y1(:,148);  %Ophir
ycql=ycqlf'; %Licor

yusib=sqrt(ycs.*zuu^.667./psi_fu(zety));yusib=yusib';
ytsibs=-sqrt(yct.*zuu^.667./psi_ft(zety)).*sign(ywt_sonb');ytsibs=ytsibs';
yqsib_oph=sign(yqsb').*sqrt(ycq.*zuu^.667./psi_ft(zety))*1e-3;yqsib_oph=yqsib_oph';%convert to g/g
yqsib_lic=sign(yqsb').*sqrt(ycql.*zuu^.667./psi_ft(zety))*1e-3;yqsib_lic=yqsib_lic';%convert to g/g
ii=find(yqsib_oph<-1e-3 | yqsib_oph>1e-3);yqsib_oph(ii)=NaN;
ii=find(yqsib_lic<-1e-3 | yqsib_lic>1e-3);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.*y1(:,123)'/1000;
yhlc=yrhoair'*Le.*y3(:,11)'/1000;
ytauc=-yrhoair'.*y1(:,117)';
ytaucx=ytauc;
ytaucy=-yrhoair'.*y1(:,119)';

%*****************   create blended turbulence latent flux
hll=(yhlc+hl_licbar)/2;
ii=find(isnan(yhlc));%case bad hlc
%jj=find(isfinite(yhlib_oph(ii)));
kk=find(isnan(hl_licbar(ii)));
   if ~isempty(kk)
          	hll(ii(kk))=yhlc(ii(kk));%only hlib is good
   end;
kk=find(isfinite(hl_licbar(ii)));
   if ~isempty(kk)
      hll(ii(kk))=hl_licbar(ii(kk))*.5;%both hlib and hl_lic good
   end;
jj=find(isfinite(hl_licbar(ii)));
if ~isempty(jj)
   hll(ii(jj))=hl_licbar(ii(jj));%only hl_lic good
end;

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

%create blended licor and ophir flux

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


%

figure;plot(jda,(yhsib+yhsc)*.5,'o',jda,yhsb,[stdt endt],[0 0]);axis([stdt endt -10 75]);
xlabel('Julian Day (2008)');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('Julian Day (2008)');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 .35]);
xlabel('Julian Day (2008)');ylabel('Stress (N/m^2)');title([jdollars ': line - bulk,  o - turbulence'])

figure;plot(ytaub,(ytauib+ytauc)*.5,'o',[-.05 .35],[-.05 .35],[-.05 .35],[0 0],[0 0],[-.05 .35]);axis([-.05 .35 -.05 .35]);
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)
