disp('DYNAMO_2011_hr_bulk_C30');

% recomputing bulk values with hourly input data using COARE 3.5
% output: yhsb, yhlb, ytaub, yLb, yusb, ytsb, yqsb, ydter, yhrain, etc...

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  y4(i,3) y4(i,5) yqs(i)  yq(i)  y4(i,8) y4(i,9) y4(i,11) zi y4(i,10) 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
    [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);
    ywbar(i) = yflx(15);
end;
yhsb = yhsb(:); yhlb = yhlb(:); ytaub = ytaub(:); yLb = yLb(:); yusb = yusb(:);
ytsb = ytsb(:); yqsb = yqsb(:); ydter = ydter(:); yhrain = yhrain(:); ywbar = ywbar(:);

yU = v1;
yP = y4(:,10);
yta = y4(:,5);
yqa = qair_p([yta,y4(:,6),yP]);

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;
ywt_son = y1(:,121);
ywt = ywt_son + 0.51*(yta+tdk).*yusb.*yqsb;
ycu = y1(:,134); % method a
ycw = y1(:,136);
% ycu = y1(:,135); % method b
% ycw = y1(:,137);
ycs = (ycu+.75*ycw)/2;% avg 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(:,138); % method a
% yct = y1(:,139); % method b
ycql = y1(:,141); % Licor cq2 method b

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);
ytsibs = -sqrt(yct.*zuu^.667./psi_ft(zety)).*sign(ywt_sonb);

yqsib_lic = sign(yqsb).*sqrt(ycql.*zuu^.667./psi_ft(zety))*1e-3; %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(:,8)/1000;  % Licor H2O covariance Hl
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;

%% yet another blended licor flux, without yhlib - later saved to 1 hr flux file
% yhlc_bl = yhlc;
% ii = find(isnan(yhlc));
% yhlc_bl(ii) = hlc_licbar(ii);
yhlc_bl = hlc_licbar;   % just take the filtered 10-min covariance result
% 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 200]);
figure; plot(jda,yhsc,'o',jda,yhsb,[stdt endt],[0 0]); axis([stdt endt -20 300]);
xlabel('Day-Of-year (2014)'); ylabel('Avg Cov Sensible Heat Flux (W/m^2)'); grid;
title([jdollars,' Hourly Sensible Heat Flux: line = bulk, o = cov']);

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

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

figure; plot(yhsb,yhsc,'o',[-50 300],[-50 300]); axis([-50 300 -50 300]);
xlabel('Bulk Sensible Heat Flux (W/m^2)'); ylabel('Covariance Sesible Heat Flux (W/m^2)'); grid;
title([jdollars,' Hourly Sensible Heat Flux, Cov vs Bulk']);

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)'); grid;
title([jdollars,' Latent Heat Flux, Blended Webb-corr Cov & ID Results vs Bulk']);

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