disp('STRATUS_2017_hr_bulk_C35');

% recomputing bulk values with hourly input data using COARE 3.5
% Note, 'best' met values from y4 are used here
% y4 = [jdy,u_best,ts_best,qs_best,ta_best,rh_best,qa_best,rs_best,rl_best,P_best,rain_best];
%        1    2      3       4       5       6       7       8       9      10        11
% input format: A = coare35vn(u,zu,t,zt,rh,zq,P,ts,Rs,Rl,lat,zi,rain,cp,sigH)
% using 'best' hourly average bulk met
[A,yU10N,yT10N,yQ10N,yu03n,yt03n,yrh03] = coare35vn(v1,zu,y4(:,5),zt,y4(:,6),zq,y4(:,10),y4(:,3),y4(:,8),y4(:,9),y1(:,23),zi,y4(:,11),NaN,NaN);
% A = [hsb hlb tau zo zot zoq L usr tsr qsr/1000 dter dqer tkt RF wbar Cd Ch Ce Cdn_10/1000 Chn_10/1000 Cen_10/1000 ug ];
%       1   2   3   4  5   6  7  8   9    10      11   12  13  14  15  16 17 18     19          20          21      22

clear yflx yhsb yhlb ytaub yLb yusb ytsb yqsb ydter yhrain

yhsb = A(:,1);
yhlb = A(:,2);
ytaub = A(:,3);
yLb = A(:,7);
yusb = A(:,8);
ytsb = A(:,9);
yqsb = A(:,10);
ydter = A(:,11);
ydqer = A(:,12);
yhrain = A(:,14);
ywbar = A(:,15);
yta = y4(:,5);
yqa = y4(:,7);
yU = v1;
yP = y4(:,10);
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 = (yhlc+yhlib_lic+hlc_licbar)/3; % just the mean of the above 3 values...

% 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 = polyval(coef,[minx,maxx]); 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')]);


