function [cd,u10]=cdnve(sp,z)
% CDNVE: computes neutral drag coefficient following Vera (1983).
% [cd,u10]=CDNVE(sp,z) computes the neutral drag coefficient and wind
% speed at 10m given the wind speed at height z. Uses the expression
% for friction velocity derived by E. Vera (1983) and published as
% eqn. 8 in Large, Morzel, and Crawford (1995), J. Phys. Oceanog., 25,
% 2959-2971. Range of fit to data is 1 to 25 m/s.
%
% INPUT: sp - wind speed [m/s]
% z - measurement height [m]
%
% OUTPUT: cd - neutral drag coefficient at 10m
% u10 - wind speed at 10m [m/s]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/26/98: version 1.1 (modified by RP)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% constants in fit for drag coefficient
A=2.717e-3;
B=0.142e-3;
C=0.0764e-3;
as_consts; % other constants
a=log(z./10)/kappa; % log-layer correction factor
tol=.001; % tolerance for iteration (m/s)
u10o=zeros(size(sp))+.1; % don't start iteration at 0 to prevent blowups.
cd=(A./u10o + B + C*u10o);
u10=sp./(1+a.*sqrt(cd));
ii=abs(u10-u10o)>tol;
while any(ii(:)),
u10o=u10;
cd=(A./u10o + B + C*u10o);
u10=sp./(1+a.*sqrt(cd)); % next iteration
ii=abs(u10-u10o)>tol; % keep going until iteration converges
end;