
%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=6.5;
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.0;!.80;
xx=ruff([pr rr]);
c=9.3*rr^.25;
s=c*pr^.5;
rt=xx/exp(von*s*(1-1/(von*xx*pr)));
zox=rt/rr*zo;
xxu=zo*exp(5*von);
zi=11;
du=11*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;
z1(1)=xxu;z1(2)=z1(1);x1(1)=qz1(1);x1(2)=qz3(n);
semilogy(qz1(1:j),z(1:j),'o',qz2(j:n),z(j:n),'--',qz3(j:n),z(j:n),x1,z1,':');
xlabel('\DeltaX u_{*a}/F_x');ylabel('z (m)');
dxx=diff(z);
dyy=diff(qz);
figure;loglog(kt1+visa,z,kt2+visa,z,dxx./dyy*us,z(2:61))
xlabel('\nu +K_x (m^2/s)');ylabel('z (m)');
   
     
