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

disp('DYNAMO_2011_da_id2')

zm = -100:.1:10;
cuxx = psi_fu(zm);
ctxx = psi_ft(zm);
% 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;

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 = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<0.6 & sp2<0.8 & rain_best<5 & (sqrt(wwj)./ugw<2+.0015*u_best.^2) & (sqrt(vvj)./ugu<2) & missing<300 & badSon<50); %& isledist>3.5 );

jplume = ones(length(jdy),1);
ii = find(sqrt(wwj)./ugw<2+.0015*U.^2 & sqrt(vvj)./ugu<2);
jplume(ii) = 0; % ship plume contam index based on wind variances, 1=bad, 0=good
jmanuv = 3*ones(length(jdy),1);
ij = find(sig_h<5 & sig_u<0.6 & sp2<0.8);
jmanuv(ij) = sig_u(ij); % ship maneuver index < 3=good


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,'O',zm,cuxx)
axis([-15 1 1 100]);
% axis([-1 1 1 10]);

figure; plot(zetxm,csm./psi_fu(zetxm),'.');
% axis([-30 1 .8 1.2]);

u10n = U-usb/.4.*(log(zuu/10)-psiu_30(zuu./Lb));
usib = sqrt(cs.*zuu^.667./psi_fu(zetx));
tsibs = -sqrt(ct.*zuu^.667./psi_ft(zetx)).*sign(wtv);
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;

tsib = tsibs-.51*(ta+273.16).*qsb;

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

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

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

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;

figure; plot(hlbb(jjj),hlib_lic(jjj)+hl_webb(jjj),'.',[-20 300],[-20 300]); axis([-20 300 -20 300]);  
xlabel('Bulk Hl (W/m^2)'); ylabel('Turb Hlib-Lic (W/m^2)');

figure; plot(hlbb(jjj),hlc_lic(jjj)+hl_webb(jjj),'.',[-20 300],[-20 300]); axis([-20 300 -20 300]);  
xlabel('Bulk Hl (W/m^2)'); ylabel('Turb Hlc-Lic (W/m^2)');

hlm = (hlc_lic+hlib_lic)*.5;	% mean of covariance and id latent heat fluxes
figure; plot(hlbb(jjj),hlm(jjj)+hl_webb(jjj),'.',[-20 300],[-20 300]); axis([-20 300 -20 300]);
xlabel('Bulk Hl (W/m^2)'); ylabel('Mean Turb Hl-Lic + Hl-Webb (W/m^2)');

clear hlmm hlbm nn dqm uu

uu=0:20;
for i=1:21
    ii = find(floor(U)==uu(i));
%   ppp = find(ushp(ii)<3 & jdy(ii)>147  & (reldir(ii) <90 & reldir (ii) >-90)& isfinite(hlm(ii)) & sig_h(ii)<4 & sqrt(wwj(ii))./ug(ii)<.5+.003*u_best(ii).^2 & sqrt(vvj(ii))./ug(ii)<1.1+0.05*u_best(ii));
    ppp = find(ushp(ii)<3  & (reldir(ii) <90 & reldir (ii) >-90) & isfinite(hlm(ii)) & sig_h(ii)<5 & rain_best(ii)<5 & sqrt(wwj(ii))./ug(ii)<2+.0015*u_best(ii).^2 & sqrt(vvj(ii))./ug(ii)<2+0.05*u_best(ii));

    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;

figure; plot(uu,hlbm,'g',uu,hlmm,'b'); xlabel('U (m/s)'); ylabel('Hl (W/m^2)');
title('hlb median in green, [(hlc+hlib)*.5 + hl-webb] median in blue');

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

% ppp=find((reldir<90 & reldir>-90) & sig_h<8 & sig_u<2 & sp2<2  & sqrt(wwj)./ugw<.55+.002*u_best.^2 & sqrt(vvj)./ugu<1.1 & isfinite(usib));
ppp = find((reldir<90 & reldir>-90) & sig_h<5 & sig_u<0.6 & sp2<0.8 & rain_best<5 & sqrt(wwj)./ugw<.7+.0015*u_best.^2 & sqrt(vvj)./ugu<1.4 & isfinite(usib));

kk = find(abs(-wuj(ppp)-usib(ppp).^2)<.05+.01*U(ppp));
%------ test on Edson
% jj2 = find((reldir<90 & reldir>-90) & sig_h<5 &sig_u<0.6 &sp2<0.8  & sqrt(ww)./ugw<.7+.0015*u_best.^2 & sqrt(vv)./ugu<1.4 & isfinite(usib));% & isledist>3.5 );
% kk2 = find(abs(-wuj(ppp)-usib(ppp).^2)<.02+.006*U(ppp));

figure; plot(jdy(ppp),-wuj(ppp),'.',jdy(ppp),usib(ppp).^2,[stdt endt],[0 0]); axis([stdt endt -.05 1]);
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 -.05 1]);
xlabel('Day-Of-year (2013)'); ylabel('u_*^2(m^2/s^2)'); title([jdollars 'filtered pass 2: ID and Cov Edson']);
%------

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

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

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