
%x=[5 0 25.0 24.5 20.2 17.5 141 419 0 600 980 10 10 10 1.0];
P=1010;
g=9.82;
ts=28.7;
qs=qsea(ts);
t=20;
q=18.5;
u=.01;
xm=[u 0 ts t qs q 141 419 0 1 1010 10 10 10];
ym=coul_lkb(xm);
zo=ym(4);
zot=ym(5);
us=ym(8);
tsr=ym(9);
qsr=ym(10);
L=ym(7);

von=0.4;
tdk=273.16;
grav=9.82;
     %*************  air constants ************
Rgas=287.1;
Le=(2.501-.00237*ts)*1e6;
cpa=1004.67;
rhoa=P*100/(Rgas*(t+tdk)*(1+0.61*q));
visa=1.325e-5*(1+6.542e-3*t+8.301e-6*t*t-4.8e-9*t*t*t);
     %************  cool skin constants  *******
Al=2.1e-5*(ts+3.2)^0.79;
be=0.026;
cpw=4000;
rhow=1022;
visw=1e-6;
tcw=0.6;
bigc=16*grav*cpw*(rhow*visw)^3/(tcw*tcw*rhoa*rhoa);
wetc=0.622*Le*qs/(Rgas*(ts+tdk)^2);
Fb=-grav/(t+tdk)*us*(tsr+.61*qsr*(t+tdk));

rr=zo*us/visa;
pr=250;
zi=11;
Ri=visa*Fb/us^4;
phiko=1/(1+(4*Ri)^.75)^.333;
du=10*visa/us*phiko;
phik=1/(phiuv(du/L)+(4*Ri)^.75)^.333;
duo=10*visa/us;
du=10*visa/us*phik;
dx=du/pr^.5;
D=visa/pr;

for i=1:5
   cc=D*phitv(dx/L)/.4/us;
   dx=(cc+sqrt(cc*cc+4*cc*du))/2;
end;

   
ix=-5:.1:1;
z=10.^ix;
zb=us*z/visa;
qz2=1/von*log(zb/zo);
qz3=qz2-1/von*(psitg(z/L)-psitg(zox/L));
kt0=z*von*us./(1+(du./z));

kt1=z*von*us./(1+(du./z))./phitv(z/L);
zp=z*us/visa;
kt2=von*visa*(zp-zi*atan(zp/zi));
[m,n]=size(ix);
for i=1:n
   if zb(i)<xx
      j=i;
      qz(i)=qz1(i);
   else
      qz(i)=qz2(i);
   end;
end;
aa=D*du;
ba=D;
ca=von*us;
Dx=D/(2*von*us);
Dl=D*sqrt((2*du/Dx-1));
t1=1/2/ca*log((aa+ba*z+ca*z.^2)/aa);
t2=(du-Dx)*2/Dl*(atan((ba+2*ca*z)/Dl)-atan(ba/Dl));
dlx=2*sqrt(D*du/us/von);
dlx=du;

figure;loglog(kt1+D,z,'.',kt0+D,z,'x',kt0+visa,z,[1e-7 10],[du du],'--',[1e-7 10],[dlx dlx],':');
axis([1e-8 1 1e-5 10]);
xlabel('\fontsize{12} D_x+K_t  (m^2/s)');ylabel('\fontsize{12} z (m)');
 set(gca,'fontsize',12);
  
     
