
%x=[5 0 25.0 24.5 20.2 17.5 141 419 0 600 980 10 10 10 1.0];
P=1010;
ts=28.7;
qs=qsea(ts);
t=26;
q=18.5;
u=.5;
aax=25;
xm=[u 0 ts t qs q 141 419 0 600 1010 10 10 10];
ym=coul_lkb(xm);
zo=ym(4);
zot=ym(5);
us=ym(8);
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);
     
rr=zo*us/visa;
pr=1.%50;
xp=[pr rr];
yy =ruff(xp);
xx=yy(1);
c=yy(3);%c=9.3*rr^.25;
s=yy(4);%s=c*pr^.5;
rt=xx/exp(von*s*(1-1/exp(xx*pr/s)));
zox=rt/rr*zo;
xxu=zo*exp(5*von);
zi=11;
du=aax*visa/us;
ix=-5:.1:1;
z=10.^ix;
zb=us*z/visa;
qz1=s*(1-exp(-zb*pr/s));
qz2=1/von*log(zb/rt);
qz3=qz2-1/von*(psitg(z/L)-psitg(zox/L));
kt1=z*von*us./(1+(du./z).^2);
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;
D=visa/pr;
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;

z1(1)=xx*visa/us;z1(2)=z1(1);x1(1)=0;x1(2)=100;
semilogy(qz1(1:j),z(1:j),'o',qz2(j:n),z(j:n),x1,z1,'--',(t1+t2)*us,z,'.',[1e-7 100],[dlx dlx],':');
xlabel('\fontsize{12} \DeltaX u_{*a}/F_x');ylabel('\fontsize{12} z (m)');
axis([ 0 30 1e-5 1.0]);
set(gca,'fontsize',12);

dxx=diff(z);
dyy=diff(qz);
figure;loglog(kt1+D,z,'.',kt1+visa,z,dxx./dyy*us,z(2:61),'x',[1e-7 10],[du du],'--',[1e-7 10],[dlx dlx],':');
axis([1e-7 .01 1e-5 .1]);
xlabel('\fontsize{12} D_x+K_t  (m^2/s)');ylabel('\fontsize{12} z (m)');
 set(gca,'fontsize',12);
  
     
