disp('CAPRICORN_2016_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 cov ID fluxes
zuu = zu;%-.7;
zety = zuu./yLb;

switch IDflag
    case 'a'
        ycu = y1(:,134);        % ID method a
        ycw = y1(:,136);
        yct = y1(:,138);
    case 'b'
        ycu = y1(:,135);        % ID method b
        ycw = y1(:,137);
        yct = y1(:,139);
end;

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);

ycql = ycqfl; % Licor cq2 - using filtered hourly avg result from hr_ave.m

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);
ytauib = yrhoair.*yusib.*yusib./ygg;

ytsibs = -sqrt(yct.*zuu^.667./psi_ft(zety)).*sign(ywt_sonb);
ytsib = ytsibs-.51*(yta+273.16).*yqsb;
yhsibs = -yrhoair*cpa.*yusib.*ytsibs;
yhsib = -yrhoair*cpa.*yusib.*ytsib;

yqsib_lic = sign(yqsb).*sqrt(ycql.*zuu^.667./psi_ft(zety))*1e-3;
ii = find(yqsib_lic<-3e-4 | yqsib_lic>0);
yqsib_lic(ii) = NaN;
yhlib_lic = -yrhoair*Le.*yusib.*yqsib_lic;

%% covariance fluxes
ywt_son = y1(:,121);
ywt = ywt_son+.51*(yta+tdk).*yusb.*yqsb;

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;     % use streamwise tau as best est of tau total
ytaucy = -yrhoair.*y1(:,119);

%% create blended turbulence latent flux
figure; maxy = max(yhlc)+50;
% yhlc is hr ave heat flux computed from hr ave filtered 10-min wq_lic
subplot(3,1,1); plot(jda,yhlc,'b.-'); title('yhlc'); ylim([0 maxy]);grid
% yhlib_lic is hr ave heat flux from hr ave filtered cq2
subplot(3,1,2); plot(jda,yhlib_lic,'b.-'); title('yhlib_ lic'); ylim([0 maxy]);grid
% hlc_licbar is hr ave of filtered 10-min hl
subplot(3,1,3); plot(jda,hlc_licbar,'b.-'); title('hlc_ licbar'); ylim([0 maxy]);grid

hll = yhlib_lic; % choose hlib as best value

% % now fix mean by filling cases where one of above is bad...
% 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;

%% 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;

%% plots
miny = min([min((yhsib+yhsc)*.5),min(yhsb)])-25; maxy = max([max((yhsib+yhsc)*.5),max(yhsb)])+25;
figure; plot(jda,(yhsib+yhsc)*.5,'ro',jda,yhsb,'g-',[stdt endt],[0 0],'k-'); axis([stdt,endt,miny,maxy]);
xlabel('Day-Of-year'); ylabel('Avg Cov & ID Sensible Heat Flux (W/m^2)'); grid;
title([jdollars,' Hourly Sensible Heat Flux: (cov+id)/2']); legend('(cov+id)/2','bulk','location','best');

miny = min([min(hll+hl_webar),min(yhlb)])-50; maxy = max([max(hll+hl_webar),max(yhlb)])+50;
figure; plot(jda,hll+hl_webar,'o',jda,yhlb,[stdt endt],[0 0],'k-'); axis([stdt,endt,miny,maxy]);
xlabel('Day-Of-year'); ylabel('Webb Corr Latent Heat Flux (W/m^2)'); grid;
title([jdollars ' Hourly Latent Heat Flux: (cov+id)/2']);legend('(cov+id)/2 + webb','bulk','location','best');

miny = min([min((ytauib+ytauc)*.5),min(ytaub)])-0.1; maxy = max([max((ytauib+ytauc)*.5),max(ytaub)])+0.1;
figure; plot(jda,(ytauib+ytauc)*.5,'o',jda,ytaub,[stdt endt],[0 0],'k-');
axis([stdt,endt,miny,maxy]); xlabel('Day-Of-year (2014)');ylabel('Stress (N/m^2)'); grid;
title([jdollars ' Hourly Stress: (cov+id)/2']); legend('(cov+id)/2','bulk','location','best');

minx = min(yhsc)-10; maxx = max(yhsc)+10;
figure; plot(yhsb,yhsc,'o',[minx maxx],[minx maxx],'k-'); axis([minx maxx minx maxx]);
xlabel('Bulk Sensible Heat Flux (W/m^2)'); ylabel('Covariance Sesible Heat Flux (W/m^2)'); grid;
title([jdollars,' Hourly Heat Flux, yhsc vs Bulk']);

minx = min(hll+hl_webar)-30; maxx = max(hll+hl_webar)+30;
figure; plot(yhlb,hll+hl_webar,'o',[minx maxx],[minx maxx],'k-'); axis([minx maxx minx maxx]);
xlabel('Bulk Latent Heat Flux (W/m^2)'); ylabel('Turbulent Latent Heat Flux (W/m^2)'); grid;
title([jdollars,' Latent Heat Flux, (hll+hl-webar) Cov & ID Results vs Bulk']);

minx = min(ytauib)-0.02; maxx = max(ytauib)+0.02;
figure; plot(ytaub,ytauib,'o',[minx maxx],[minx maxx],'k-',[minx maxx],[0 0],'k-',[0 0],[minx maxx],'k-');
axis([minx maxx minx maxx]); xlabel('Bulk Stress (N/m^2)'); ylabel('ID Stress (N/m^2)'); grid;
title([jdollars,' Hourly Stress ID vs Bulk']);

%% repeat heat flux QC plots using hourly values

lapse = grav/cpa; yTheta10 = yT10N + lapse*10; % Pot temp @ 10m
yqs_adj = y4(:,4) - ydqer; % qs - dqer
yts_adj = y4(:,3) - ydter; % ts - dter

ywt_U = ywt./yU10N;   % wt_son/u10n
ydTs = yts_adj-yTheta10;
idx = (isfinite(ywt_U) & isfinite(ydTs) & ydTs<3 & ydTs>-2); %wT plot criteria
minx = min(ydTs(idx))-0.5; maxx = max(ydTs(idx))+0.5;
coef = polyfit(ydTs(idx),ywt_U(idx),1); liney = coef(1)*[minx,maxx] + coef(2); xIntcpt = coef(2)/-coef(1);
figure; plot(ydTs(idx),ywt_U(idx),'r.',[minx,maxx],liney,'b'); ylim([-0.01,0.01]); grid;
hold on; plot([0,0],[-0.01,0.01],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wT/U_{10n}'); xlabel('(Ts - \Theta_{10})');
title(['Hourly Ch_{10N}x10^{-3} = ',num2str(coef(1)*1000,'%4.3f'),', y-Int = ',num2str(coef(2),'%3.2e'),', x-Int = ',num2str(xIntcpt,'%4.3f')]);

ywq_U = (y3(:,8) - yqa.*ywbar)./yU10N;   % wq_lic/u10n
ydQs = yqs_adj-yQ10N;
idx = isfinite(ywq_U) & isfinite(ydQs) & ywq_U<4*nanstd1(ywq_U); %wq plot criteria
minx = min(ydQs(idx))-1; maxx = max(ydQs(idx))+1;
coef = polyfit(ydQs(idx),ywq_U(idx),1); liney = polyval(coef,[minx,maxx]); xIntcpt = coef(2)/-coef(1);
figure; plot(ydQs(idx),ywq_U(idx),'r.',[minx,maxx],liney,'b'); ylim([-0.01,0.015]); grid;
hold on; plot([0,0],[-0.015,0.015],'k',[minx,maxx],[0,0],'k'); xlim([minx,maxx]);
ylabel('wq/U_{10N}'); xlabel('(Qs - Q_{10N})');
title(['Hourly Ce_{10N}x10^{-3} = ',num2str(coef(1)*1000,'%4.3f'),', y-Int = ',num2str(coef(2),'%3.2e'),', x-Int = ',num2str(xIntcpt,'%4.3f')]);

