prtit=0;
tdk=273.16;
Rgas=287.1;
grav=9.83;
kon=0.4;%von karman constant
cpa=1004.67;%air specific heat
barpress0=barpress(2);


% ssti=29;th20=24;u20=7.5;rh20=90;
% t=th20-.0098*20;Q=qair_p3([t barpress rh20/100])/1e3;
Le   = (2.501-.00237*ssti)*1e6;
rs=300;
rl=400;
zi=600;
zref_u=1e-4;
zref_t=1e-4;
zref_q=1e-4;


rhoa = barpress0*100./(Rgas*(t+tdk).*(1+0.61*Q));

%%%%%%%%%%%%%%%   compute fluxes with fixed sst
A=[];A=coare31vn_ref3(u20,20,th20-.0098*20,20,rh20,20,barpress0,ssti,rs,rl,lat,zi,zref_u,zref_t,zref_q,-1);
%A=[usr tau hsb hlb hbb hsbb tsr qsr zo zot zoq Cd Ch Ce L zet dter tkt Urf Trf Qrf RHrf];   
%   1    2  3   4   5   6    7   8   9   0  1   2  3  4  5 6   7   8    9   0   1  2
usr=A(:,1);
taub=A(:,2);
hsb=-A(:,3);
hlb=-A(:,4);
Cd=A(:,12);
Ce=A(:,14);
zo=A(:,9);
zot=A(:,10);
zoq=A(:,11);
L=A(:,15);
%%%%%%%%%%%%%%%%%%%%  bulk values at surface
Urf=A(:,19);
Trf=A(:,20);
Qrf=A(:,21);
RHrf=A(:,22);
tsr=A(:,7);
qsr=A(:,8);
fdg=A(:,23);
td=tdew([barpress0' Qrf/.98]);%dew point at inferred surface q
cehf=Ce./sqrt(Cd);

ij=find(zzx>5 & zzx<60);
zk=16.;pht=1-(zk./L'+.005).*(psit_30(zk./L'+.01)-psit_30(zk./L'))/.01;
[pp,s]=polyfit(log(zzx(ij)),tsz(ij),1);tszz=polyval(pp,log(zzx));
[pp1,s]=polyfit(log(zzx(ij)),tsz(ij),2);%tszz=polyval(pp,log(zzx));


% figure;semilogy(tsz,zzx,pp(2)+pp(1)/pht*(log(zzx)-psit_30(zzx/L)),zzx,th20,20,'x',ssti,zot,'o');xlim([th20-1 ssti+1]);ylim([1 250]);
% xlabel('\theta (z)');ylabel('z(m)');
% legend('Measured profile','Linear Fit','Ta','SST');

ath=pp(2);thg=pp(1);
thg1=pp1(2)+pp1(1)*log(zk)*2;%thg=thg1
tss=-thg./pht*.4./cehf'+th20;
[ssti pp(2)-4*pp(1) pp(2)+pp(1)/pht*log(zot) th20-thg./pht*.4./cehf'] 

qssi=qair_p3([ssti barpress0 .98]);

ij=find(zzx>5 & zzx<60 & ~isnan(qsz));
[ppq,s]=polyfit(log(zzx(ij)),qsz(ij),1);qszz=polyval(ppq,log(zzx));


aq=ppq(2);qg=ppq(1);
qss=-qg./pht*.4./cehf'+q20;
[qssi ppq(2)-4*ppq(1) ppq(2)+ppq(1)/pht*log(zot) q20-qg./pht*.4./cehf'] 
tdq=tdew([barpress0' qss/.98])

sst0=ath-4*thg;%pp(2)-4*pp(1);
ssq0=aq-4*qg;%ppq(2)-4*ppq(1);tdq0=tdew([barpress' ssq0/.98]);

sst1=sst0;
for i=1:20
A1=[];A1=coare31vn_ref3(u20,20,th20-.0098*20,20,rh20,20,barpress0,sst1,rs,rl,lat,zi,zref_u,zref_t,zref_q,-1);
%A=[usr tau hsb hlb hbb hsbb tsr qsr zo zot zoq Cd Ch Ce L zet dter tkt Urf Trf Qrf RHrf];   
%   1    2  3   4   5   6    7   8   9   0  1   2  3  4  5 6   7   8    9   0   1  2
usr1=A1(:,1);
taub1=A1(:,2);
hsb1=-A1(:,3);
hlb1=-A1(:,4);
Cd1=A1(:,12);
Ce1=A1(:,14);
zo1=A1(:,9);
zot1=A1(:,10);
zoq1=A1(:,11);
L1=A1(:,15);
%%%%%%%%%%%%%%%%%%%%  bulk values at surface
Urf1=A1(:,19);
Trf1=A1(:,20);
Qrf1=A1(:,21);
RHrf1=A1(:,22);
tsr1=A1(:,7);
qsr1=A1(:,8);
fdg1=A1(:,23);
td1=tdew([barpress0' Qrf1/.98]);%dew point at inferred surface q
cehf1=Ce1./sqrt(Cd1);
zk=16;pht1=1-zk./L1'.*(psit_30(zk./L1'+.01)-psit_30(zk./L1'))/.01;
sst2=th20-thg./pht1*.4./cehf1';
ssq2=q20-qg./pht1*.4./cehf1';
td2=tdew([barpress0' ssq2/.98]);%dew point at inferred surface q
sst1=td2;
sst1=sst2;
end;
pht2=1-(zk./L1'+.01).*(psit_30(zk./L1'+.02)-psit_30(zk./L1'+.01))/.01;
thg2=thg*(-1/zk+(pht2-pht1)/.01/pht1);

[ssti sst0 sst2 td2]
[hsb hsb1 -.4*thg*usr1/pht1*rhoa*cpa  hlb hlb1 -.4*qg*usr1/pht1*rhoa*Le/1e3 ]

figure;semilogy(tsz,zzx,'x',tszz,zzx,sst2+tsr1/.4.*(log(zzx./zot1)-psit_30(zzx/L1)),zzx,'--',th20,20,'x',sst2,1,'o');xlim([th20-1 ssti+1]);ylim([1 500]);xlabel('Theta (C)');ylabel('z (m)');title(['ENRR2016 JD= ' num2str(jd(i)) ' hr= ' num2str(hh(i)) ' lat= ' num2str(lat) ' lon= ' num2str(lon)]);
figure;semilogy(qsz,zzx,'x',qszz,zzx,ppq(2)+ppq(1)/pht*(log(zzx)-psit_30(zzx/L)),zzx,q20,20,'x',ssq2,1,'o');xlim([q20-1 ssq2+1]);ylim([1 500]);xlabel('Spec Hum (g/kg)');ylabel('z (m)');title(['ENRR2016 JD= ' num2str(jd(i)) ' hr= ' num2str(hh(i)) ' lat= ' num2str(lat) ' lon= ' num2str(lon)]);

