%function zo=getzom(x)
%x=[30 29 5 10 .17 -.1 1.5e-5]
tw=x(1);
ta=x(2);
u=x(3);
zref=x(4);
ustr=x(5);
Lin=x(6);
nu=x(7);

betag=1;
betac=1;
betas=0;
wa=20;

w=.8;
w_eql=0.81;
zz=1e4;
s=0.349;
crit=1e-3;
b=0.019;
kv=.4;
g=9.81;
h=1e3;
denwat=1027;
u_sfc=0.55*ustr;
sfcten=(-1.55e-4*(273.16+tw))+0.118;
cmin=sqrt(2*sqrt(g*sfcten/denwat));

for jam=1:3
if ~(betag==0) & wa*ustr<cmin
   betag=0;
   betac=0;
   betas=1;
end;
cmin=sqrt(2*sqrt(g*sfcten/denwat));
zo0=betas*.11*nu/abs(ustr)+betac*b*sfcten/(ustr*ustr*denwat)+.48/wa*ustr*ustr/g;
dum=ustr*ustr*wa*wa;
if dum>cmin*cmin
   wave_len=pi*(dum+(sqrt(dum*dum-cmin^4)))/g;
else
   wave_len=pi*2*dum/g;
end;
period=wave_len/(wa*ustr);
hp=0.062*sqrt(g*ustr*period^3);
hp=min(hp,137*wave_len);
u_sfc=betag*(.5*2*pi*hp/wave_len)^2*ustr*wa+(1-betag)*.55*ustr;
phi=phiu_clay(betag*hp,Lin,hp/zo0);
hzzg=hp*zz/zref;
hzzg=max(hzzg,1e-4);
wa=(log(hzzg+1)+phi)/(0.4*w*w_eql);
wa=min(wa,150);wa=max(wa,1.08);
zz=zref/((betas*.11*nu/abs(ustr))+betac*b*sfcten/(ustr*ustr*denwat)+betag*.33/wa*ustr*ustr/g);
if zref/zz>25;zz=zref/25;end;
if zref/zz<1e-5;zz=zref/1e-5;end;
if betag==0 | wa*ustr>cmin;break;end;
end;
zo=zref/zz;

