function [u_star,n,small,diff] = charnock(u10,alpha_c,b,max_it)

% n = max_it, did NOT converge
% n < max_it, converged!

% model parameters
% alpha_c 	Charnock constant (0.011 - Liu & Tang; 0.018 ECMWF; 0.013 Zeng etc.
% b		smooth flow constant (0.11 all but NCEP)
% permanent constants
% nu is air kinematic viscosity
g = 9.81;
nu = 1.5e-5;
vonK = 0.4;
eps1 = 1e-8;
eps2 = 1e-6;

% initialize
u_star = 0.04 * u10;

for n = 1:max_it
   z0 = alpha_c* u_star.^2/g + b*nu./u_star;
   u_star_new =  - u10 * vonK ./ log(z0/10);
   diff = abs((u_star_new - u_star)./(u_star+eps1));
   small = max(diff);
   if small < eps2
      break
   else
      u_star = u_star_new;
   end
end

u_star(diff>eps2) = NaN;