% inertial dissipation flux processing
% assumes dt array of formatted results in in workspace

disp('CAPRICORN_2016_da_id2')

zm = -100:.1:10;
cuxx = psi_fu(zm);
ctxx = psi_ft(zm);
%figure; semilogy(zm,cuxx,zm,ctxx);

zuu = zu;%-.7;
zqq = zq;

zetx = zuu./dt(:,17);
usb = dt(:,18);
tsb = dt(:,19);
qsb = dt(:,20);
hsbb = dt(:,11);
hlbb = dt(:,12);
taubb = dt(:,13);
cdn10bb = dt(:,29);
Lb = dt(:,17);
tsbs = tsb+.51*(ta_best+273.16).*qsb; % sonic T

wtv = -usb.*(tsb+.61.*(ta_best+273.16).*qsb);

ws = (9.83./(ta_best+273.16).*abs(wtv)*600).^.333;  % convective velocity
ug = sqrt(usb.^2+ws.^2);
ugw = sqrt(1.5*usb.^2+ws.^2);
ugu = sqrt(8*usb.^2+1.*ws.^2);
gg = sqrt(1+1.6*ws.^2./u_best.^2);

% base selection criteria: adjust this test to exclude 10-min segments based on
% data quality parameters or experimental conditions. jjj for PSD, kkk for UH
jjj = find((reldir<rdir_hi & reldir>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim &...
            sp2<sp2_lim & rain_best<rain_lim & std_rwdir<sig_rwd_lim &...
            (sqrt(wwj)./ugw < 2+.0015*u_best.^2) & (sqrt(vvj)./ugu<2) & ...
            missing<missing_lim & badSon<badSon_lim);

jplume = ones(length(jdy),1);
ii = find(sqrt(wwj)./ugw<2+.0015*u_best.^2 & sqrt(vvj)./ugu<2);
jplume(ii) = 0; % ship plume contam index based on wind variances, 1=bad, 0=good
jplume(isnan(jplume)) = 1;

jmanuv = 3*ones(length(jdy),1);
ij = find(sig_h<sigh_lim & sig_u<sigu_lim & sp2<sp2_lim);
jmanuv(ij) = sig_u(ij); % ship maneuver index <3 good
jmanuv(isnan(jmanuv)) = 3;

zzl = [-60 -30 -20 -8 -5 -4 -3 -2 -1 -.7 -.5 -.3 -.2 -.1 -.05 0 .05 .1 1];
nn = length(zzl)-1;
clear zetxm csm cnt;
for i=1:nn
   ii = find(zetx(jjj)>zzl(i) & zetx(jjj)<zzl(i+1) & isfinite(cs(jjj)));
   if isfinite(jjj(ii));
      zetxm(i) = median(zetx(jjj(ii)));
      csm(i) = median(cs(jjj(ii))./usb(jjj(ii)).^2*zu^.667);
      cnt(i) = length(jjj(ii));
   else
      zetxm(i) = NaN;
      csm(i) = NaN;
      cnt(i) = NaN;
   end;
end;

figure; semilogy(zetxm,csm,'rO',zm,cuxx,'k-'); grid;
axis([-15 1 1 100]); xlabel('z/L'); ylabel('Psi-fu');
legend('Cs/usb * zu^{2/3}','psi-fu (z/L)','location','best');

figure; plot(zetxm,csm./psi_fu(zetxm),'rO'); grid;
xlabel('z/L'); ylabel('Csm/Psi-fu');

u10n = u_best-usb/.4.*(log(zuu/10)-psiu_30(zuu./Lb));
usib = sqrt(cs.*zuu^.667./psi_fu(zetx)); usib = real(usib);
tsibs = -sqrt(ct.*zuu^.667./psi_ft(zetx)).*sign(wtv); tsibs = real(tsibs);
qsib_lic = sqrt(cql.*zqq^.667./psi_ft(zetx))*1e-3.*sign(qsb); % convert to g/g   %cql+0.0075...

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

figure; plot(u10n(jjj),(usib(jjj)./u10n(jjj)).^2,'.',u10n(jjj),cdn10bb(jjj),'.'); axis([0 20 0 3e-3]);
xlabel('U_{10n}'); ylabel('ID Cd_{10n}'); title([jdollars ' 10-min']); grid;

tsib = tsibs-.51*(ta+273.16).*qsb;
hsib = -rhoair*cpa.*(tsib).*usib;
bfbb = -rhoair*cpa.*(tsb+0.51*(ta+tdk).*qsb).*usb;
bfib = -rhoair*cpa.*(tsibs).*usib;

miny = min(hsib(jjj))-10; maxy = max(hsib(jjj))+10;
figure; plot(hsbb(jjj),hsib(jjj),'.',[miny maxy],[miny maxy]); axis([miny maxy miny maxy]);
xlabel('H_s bulk (W/m^2))'); ylabel('H_s ID (W/m^2))'); title([jdollars ' 10-min']); grid;

hsc = cpa*rhoair.*wt;
hlib_lic = -rhoair*Le.*(qsib_lic).*usib;
ii = find(sgq_lic<5e-3 | sgq_lic>.7 | hlib_lic<-200); hlib_lic(ii)=NaN;

miny = min(hlib_lic(jjj))-30; maxy = max(hlib_lic(jjj))+30;
figure; plot(hlbb(jjj),hlib_lic(jjj)+hl_webb(jjj),'.',[miny maxy],[miny maxy]); axis([miny maxy miny maxy]);
xlabel('Bulk Hl (W/m^2)'); ylabel('Turb Hlib-Lic (W/m^2)'); grid;

miny = min(hlc_lic(jjj))-30; maxy = max(hlc_lic(jjj))+30;
figure; plot(hlbb(jjj),hlc_lic(jjj)+hl_webb(jjj),'.',[miny maxy],[miny maxy]); axis([miny maxy miny maxy]);
xlabel('Bulk Hl (W/m^2)'); ylabel('Turb Hlc-Lic (W/m^2)'); grid;

hlm = hlib_lic;    % choose id latent heat as best value
miny = min(hlm(jjj))-30; maxy = max(hlm(jjj))+30;
figure; plot(hlbb(jjj),hlm(jjj)+hl_webb(jjj),'.',[miny maxy],[miny maxy]); axis([miny maxy miny maxy]);
xlabel('Bulk Hl (W/m^2)'); ylabel('Best Hl (ID) (W/m^2)'); grid;

clear hlmm hlbm nn dqm uu

uu=0:20;
for i=1:21
    ii = find(floor(U)==uu(i));
    ppp = find(ushp(ii)<ushp_lim & (reldir(ii)<rdir_hi & reldir(ii)>rdir_lo) & ...
        isfinite(hlm(ii)) & sig_h(ii)<sigh_lim & sp2(ii)<sp2_lim & ...
        sqrt(wwj(ii))./ug(ii)<2+.0015*u_best(ii).^2 & sqrt(vvj(ii))./ug(ii)<2+0.05*u_best(ii) & ...
        badSon(ii)<badSon_lim);
    if isempty(ppp)
        hlmm(i) = NaN;
        hlbmb(i) = NaN;
        hlbm(i) = NaN;
        nn(i) = 0;
        dqm(i) = NaN;
    else
        hlmm(i) = median(hlm(ii(ppp))+hl_webb(ii(ppp)));
        hlbmb(i) = median(Hlb(ii(ppp))+hl_webb(ii(ppp)));
        hlbm(i) = median(hlbb(ii(ppp)));
        nn(i) = length(ppp);
        dqm(i) = median(qse(ii(ppp))-qa(ii(ppp)))-1e3*mean(dt(ii(ppp),22));%
    end;
end;

miny = min(hlmm)-30; maxy = max(hlmm)+30;
figure; plot(uu,hlbm,'g',uu,hlmm,'b'); xlabel('U (m/s)'); ylabel('hlib+webb (W/m^2)'); grid;
ylim([miny maxy]); title('Wind Bin Hl: hlb median in green, hlib+webb median in blue');

figure; plot(hlbm,hlmm,'o',[0 maxy],[0 maxy]); xlabel('Bulk Hl (W/m^2)'); ylabel('hlib (W/m^2)');
axis([0 maxy 0 maxy]); grid;

ppp = find((reldir<rdir_hi & reldir>rdir_lo) & sig_h<sigh_lim & sig_u<sigu_lim & ...
            sp2<sp2_lim & (sqrt(wwj)./ugw<.7+.0015*u_best.^2) & std_rwdir<sig_rwd_lim &...
            (sqrt(vvj)./ugu<1.4) & isfinite(usib) & badSon<badSon_lim);
kk = find(abs(-wuj(ppp)-usib(ppp).^2)<.02+.006*u_best(ppp));

miny = min([min(-wuj(ppp)),min(usib(ppp).^2)])-0.02; maxy = max([max(-wuj(ppp)),max(usib(ppp).^2)])+0.02;
figure; plot(jdy(ppp),-wuj(ppp),'.',jdy(ppp),usib(ppp).^2,[stdt endt],[0 0]); axis([stdt endt miny maxy]); grid;
xlabel('Day-Of-year (2013)'); ylabel('u_*^2(m^2/s^2)'); title([jdollars 'filtered pass 1: ID and Cov Edson']);

figure; plot(jdy(ppp(kk)),-wuj(ppp(kk)),'.',jdy(ppp),usib(ppp).^2,[stdt endt],[0 0]); axis([stdt endt miny maxy]);
xlabel('Day-Of-year (2013)'); ylabel('u_*^2(m^2/s^2)'); title([jdollars 'filtered pass 2: ID and Cov Edson']); grid;
%------

figure; semilogy(zetx(ppp),cql(ppp)./(1e3*qsb(ppp)).^2.*zuu^.67,'.',zetx(ppp),psi_ft(zetx(ppp)),'.'); axis([-10 2 0.01 5e2])
xlabel('z/L'); ylabel('Dimensionless C-qlic^2'); grid;
figure; semilogy(zetx(ppp),ct(ppp)./(tsbs(ppp)).^2.*zuu^.67,'.',zetx(ppp),psi_ft(zetx(ppp)),'.'); axis([-10 2 0.01 5e2])
xlabel('z/L'); ylabel('Dimensionless Ct^2'); grid;

figure; semilogy(zetx(ppp),yxx(ppp,12)./abs(1e3*qsb(ppp)),'.',zetx(ppp),phi_tt(zetx(ppp)),'.'); axis([-10 2 0.1 5e1])
xlabel('z/L'); ylabel('Dimensionless \sigma qLic'); grid;
figure; semilogy(zetx(ppp),sqrt(tt(ppp))./abs(tsbs(ppp)),'.',zetx(ppp),phi_tt(zetx(ppp)),'.'); axis([-10 2 0.1 5e1])
xlabel('z/L'); ylabel('Dimensionless \sigma t'); grid;

figure; semilogy(zetx(ppp),yxx(ppp,12)./abs(1e3*qsb(ppp))./(sqrt(tt(ppp))./abs(tsbs(ppp))),'.',[-10 2],[.85 .85]);
axis([-10 2 0.1 1e1]); xlabel('z/L'); ylabel('Dimensionless \sigma qLic/ \sigma t'); grid;
