% Code to generate tghe transfer coefficient plots Cd and Ce

%% First run COARE 3.5 for bulk model estimates and to be consitent through the cruises
% output from coare35vn: 42 columns
% y = [usr tau hsb hlb hbb hsbb hlwebb tsr qsr zo zot zoq Cd Ch Ce  L zet dter dqer tkt Urf Trf Qrf RHrf UrfN Rnl Le rhoa UN U10 U10N Cdn_10/1000 Chn_10/1000 Cen_10/1000 RF Qs Evap T10N Q10N RH10 ug wbar];
%      1   2   3   4   5   6    7      8   9   10 11  12  13 14 15  16 17  18   19  20  21  22  23   24   25  26  27  28  29  30  31       32          33          34     35  36  37   38   39  40   41  42
% with additional variables from WarmLayer - set to zero if jwarm = 0
% 43: dt_wrm - warming across entire warm layer deg C
% 44: tk_pwp - warm layer thickness m
% 45: dsea   - dT due to warming at depth of Tsea such that Tsea_true = Tsea + dsea

% constants and other variables for the bulk flux calculations
Rgas = 287.1;
cpa = 1004.67;
be = 0.026;
cpw = 4000;
rhow = 1022;
visw = 1e-6;
tcw = 0.6;
tdk = 273.15;
rhoair = press.*100./(Rgas*(ta+tdk).*(1+0.61*qa/1000));
Le = (2.501-.00237*nanmedian1(tsnk))*1e6;
grav=grv(lat);%9.82;

zref = 10;   % ref height, m
jwarm = 1;   % set to zero to skip warm layer
zi = 600;    % inversion ht default
ts_depth = 0.05;  % bulk water temperature sensor depth _ pay attention to which sensor is selected (snake or tsg)

if jwarm
    y = coare35vnWarm_lb(jdyx,U,zu,ta,zt,RH,zq,press,tsnk,rs,rl,lat,lon,zi,org,ts_depth,zref);
end
ust=y(:,1);
zou=y(:,10);
zoq=y(:,12);
cdn10=y(:,32);
cen10=y(:,34);
hlb=y(:,4);
taub=y(:,2);
hsb=y(:,3);
dter=y(:,18);
dqer=y(:,19)*1e3;  
u10n=y(:,31);
dq10n=qse-y(:,39)-dqer;
dth10n=tsnk-y(:,38)-dter;
ugb=y(:,41);  %1.2.W* in coare35vn_lb
gf=sqrt(1+ugb.^2./U.^2); %gustiness factor


%% Play with the stress and Cd
%Bin average selected data
clear uu
taucovib=(taucx+tauib)/2;  %average cov and ID a la Chris
uu(1:13)=0:12;  uu(14:19)=14:2:24;
ii=find(isfinite(taucx) &  isfinite(taub)&  isfinite(tauib) & (reldir>300 | reldir<60) & ushp<3.5 & abs(tiltx)<15 & Jm<1.5  & J==0 & org<5);

[u10nx, taucx_mean, taucx_med, taucx_SD, taucx_SE, taucx_SDM, taucx_SEM, taucx_nm, taucx_sig] = binave_lb2(u10n(ii),taucx(ii),uu);
[u10nx, taucy_mean, taucy_med, taucy_SD, taucy_SE, taucy_SDM, taucy_SEM, taucy_nm, taucy_sig] = binave_lb2(u10n(ii),taucy(ii),uu);
[u10nx, taub_mean, taub_med, taub_SD, taub_SE, taub_SDM, taub_SEM, taub_nm, taub_sig] = binave_lb2(u10n(ii),taub(ii),uu);
[u10nx, tauib_mean, tauib_med, tauib_SD, tauib_SE, tauib_SDM, tauib_SEM, tauib_nm, tauib_sig] = binave_lb2(u10n(ii),tauib(ii),uu);
[u10nx, taucovib_mean, taucovib_med, taucovib_SD, taucovib_SE, taucovib_SDM, taucovib_SEM, taucovib_nm, taucovib_sig] = binave_lb2(u10n(ii),taucovib(ii),uu);
[u10nx, Cdn10_mean, Cdn10_med, Cdn10_SD, Cdn10_SE, Cdn10_SDM, Cdn10_SEM, Cdn10_nm, Cdn10_sig] = binave_lb2(u10n(ii),cdn10(ii),uu);

%do plots wiht at least 5 pts per bins
ii=find(taucx_nm>5);
figure(1);plot(taub_med(ii),taucx_med(ii),'o',taub_mean(ii),taucx_mean(ii),'x',taub_med(ii),tauib_med(ii),'d',taub_med(ii),taucy_med(ii),'<',[-.10 1.5],[0 0],'--k',[-.10 1.5],[-0.1 1.5],'k');
set(gca,'fontsize',13);
axis([-.025 1.4 -.075 1.4]);
xlabel('Bulk Streamwise Stress (Nm^{-2})');ylabel('Turb Streamwise Stress (Nm^{-2})');%title([jdollars  '  x-mean  o-median  diamond-ID'  ]);
legend('\tau_x median','\tau_x mean','\tau_{id} median','\tau_y median')

cdn10_ib_med=tauib_med.*Cdn10_med./taub_med;
cdn10_cov_med=taucx_med.*Cdn10_med./taub_med;
cdn10_covib_med=taucovib_med.*Cdn10_med./taub_med;

taucovib_sig(taucovib_sig==0)=NaN;
cdn10_covib_med(cdn10_covib_med==0)=NaN;
Cdn10_mean(Cdn10_mean==0)=NaN;
cdn10_ib_med(cdn10_ib_med==0)=NaN;cdn10_cov_med(cdn10_cov_med==0)=NaN;

figure(2);errorbar(u10nx,cdn10_covib_med,(taucovib_sig./sqrt(taucovib_nm)).*Cdn10_mean./taub_med,'o');
set(gca,'fontsize',13);
xlabel('U_{10n} (ms^{-1})');ylabel('{C_{d10n}}');
title([ 'median: o (<Fcov>/<Fbulk>)*Cd10n, diamond=<Cd>, red = Cd10n']);
hold on;
h=plot(u10nx,Cdn10_mean,'r');set(h,'LineWidth',2);
plot(u10nx,cdn10_ib_med,'x',u10nx,cdn10_cov_med,'d');
hold off
axis([0 25 0.5e-3 3e-3]);
legend('taucovib_med.*Cdn10_med./taub_med','cdn10 bulk mean','tauib_med.*Cdn10_med./taub_med','taucx_med.*Cdn10_med./taub_med','Location','Best')


wbarQ=1.61*hlb/Le./(1+1.61*qa/1e3)./rhoair;
wbarT=hsb./rhoair/cpa./(ta+tdk);%formulation in hlb already includes webb
wbar=wbarQ+wbarT;
     %wbar=1.61*hlb/Le/rhoa+(1+1.61*Q)*hsb/rhoa/cpa/ta;
hl_webb=rhoair.*wbar.*qa*Le/1000;
ii=find(isnan(hlwebb));
hlwebb(ii)=hl_webb(ii);

hll=hlc;
ii=find(isnan(hlc) & isfinite(hlib));
hll(ii)=hlib(ii);
ii=find(isfinite(hlc) & isfinite(hlib));
hll(ii)=(hlib(ii)+hlc(ii))*.5;  %average cov and ID a la Chris
hl_webar=hlwebb';

%Bin average selected data
ii=find(isfinite(hll) & (reldir>270 | reldir<90) & isfinite(dq10n) & org<1 & ushp<9.5 & Jm<1.5  & J==0 & abs(dq10n)>1); 
ii=find(isfinite(hll) & (reldir>285 | reldir<75) & isfinite(dq10n) & org<1 & ushp<9.5 & Jm<1.5  & J==0 & abs(dq10n)>1); 
hlcovib=(hll(ii)+hl_webar(ii)');
[u10nx, hlcovib_mean, hlcovib_med, hlcovib_SD, hlcovib_SE, hlcovib_SDM, hlcovib_SEM, hlcovib_nm, hlcovib_sig] = binave_lb2(u10n(ii),hlcovib,uu);
[u10nx, hlc_mean, hlc_med, hlc_SD, hlc_SE, hlc_SDM, hlc_SEM, nm, hlc_sig] = binave_lb2(u10n(ii),hlc(ii)+hlwebb(ii),uu);
[u10nx, hlb_mean, hlb_med, hlb_SD, hlb_SE, hlb_SDM, hlb_SEM, hlb_nm, hlb_sig] = binave_lb2(u10n(ii),hlb(ii),uu);
[u10nx, hlib_mean, hlib_med, hlib_SD, hlib_SE, hlib_SDM, hlib_SEM, hlib_nm, hlib_sig] = binave_lb2(u10n(ii),hlib(ii)+hlwebb(ii),uu);
[u10nx, gf_mean, gf_med, gf_SD, gf_SE, gf_SDM, gf_SEM, gf_nm, gf_sig] = binave_lb2(u10n(ii),gf(ii),uu);
% [u10nx, dq_mean, dq_med, dq_SD, dq_SE, dq_SDM, dq_SEM, dq_nm, dq_sig] = binave_lb2(u10n(ii),dq10n(ii).*rhoair(ii),uu);
[u10nx, dq_mean, dq_med, dq_SD, dq_SE, dq_SDM, dq_SEM, dq_nm, dq_sig] = binave_lb2(u10n(ii),dq10n(ii),uu);
[u10nx, dt_mean, dt_med, dt_SD, dt_SE, dt_SDM, dt_SEM, dt_nm, dt_sig] = binave_lb2(u10n(ii),dth10n(ii),uu);

vari=(hll(ii)+hl_webar(ii)')./dq10n(ii)./rhoair(ii)./u10n(ii)*1e3./gf(ii)./Le;
[u10nx, Ce_mean, Ce_med, Ce_SD, Ce_SE, Ce_SDM, Ce_SEM, Ce_nm, Ce_sig] = binave_lb2(u10n(ii),vari,uu);
[u10nx, Cen10_mean, Cen10_med, Cen10_SD, Cen10_SE, Cen10_SDM, Cen10_SEM, Cen10_nm, Cen10_sig] = binave_lb2(u10n(ii),cen10(ii),uu);
[u10nx, grav_mean, grav_med, grav_SD, grav_SE, grav_SDM, grav_SEM, grav_nm, grav_sig] = binave_lb2(u10n(ii),grav(ii),uu);
[u10nx, rhoair_mean, rhoair_med, rhoair_SD, rhoair_SE, rhoair_SDM, rhoair_SEM, rhoair_nm, rhoair_sig] = binave_lb2(u10n(ii),rhoair(ii),uu);

%do plots wiht at least 5 pts per bins     
ii=find(hlcovib_nm>5);
figure(3);h=plot(hlb_med(ii),hlcovib_med(ii),'o',hlb_mean(ii),hlcovib_mean(ii),'x',[0 250],[0 250]);axis([0 250 0 250])
set(gca,'fontsize',13);set(h,'LineWidth',2);xlabel('Bulk Hl (Wm^{-2})');ylabel('Turbulent Hl (Wm^{-2})');
title(['x-mean  o-median'  ]);

figure;errorbar(u10nx,hlcovib_med./dq_med/Le./u10nx*1e3./gf_med./rhoair_med,hlcovib_sig./sqrt(hlcovib_nm)./dq_med/Le./u10nx*1e3./gf_med./rhoair_med)
hold on;plot(u10nx(Cen10_mean~=0),Cen10_mean(Cen10_mean~=0),'r','LineWidth',2);
hold on;plot(u10nx(Ce_med~=0),Ce_med(Ce_med~=0),'d','LineWidth',2);
axis([0 21 0.5e-3 2e-3])
hold on;errorbar(u10nx,(hlcovib_med./hlb_med.*Cen10_mean),hlcovib_sig./sqrt(hlcovib_nm)./dq_med/Le./u10nx*1e3./gf_med./rhoair_med,'m+')
legend('Ce from averaged cov+ID/Dq_{10n}...','Ce Bulk averaged','Ce averaged from cov+ID/Dq_{10n}...','Ce from averaged cov+ID/Bulk * Ce Bulk averaged')
set(gca,'fontsize',13);xlabel('U_{10n} (ms^{-1})');ylabel('{C_e}');

