function [cd,u10]=cdntc(sp,z,Ta)
% CTDTC: computes the neutral drag coefficient following Smith (1988).
% [cd,u10]=CDNTC(sp,z,Ta) computes the neutral drag coefficient and 
% wind speed at 10m given the wind speed and air temperature at height z 
% following Smith (1988), J. Geophys. Res., 93, 311-326. 
%
% INPUT:   sp - wind speed  [m/s]
%          z - measurement height [m]
%          Ta - air temperature (optional)  [C] 
%
% 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 (vectorized by RP)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get constants
as_consts; 

if nargin==2,
  Ta=Ta_default;
end;

% iteration endpoint
tol=.00001;

visc=viscair(Ta);

% remove any sp==0 to prevent division by zero
i=find(sp==0);
sp(i)=.1.*ones(length(i),1);

% initial guess
ustaro=zeros(size(sp));
ustarn=.036.*sp;

% iterate to find z0 and ustar
ii=abs(ustarn-ustaro)>tol;
while any(ii(:)),
  ustaro=ustarn;
  z0=Charnock_alpha.*ustaro.^2./g + R_roughness*visc./ustaro;
  ustarn=sp.*(kappa./log(z./z0));
  ii=abs(ustarn-ustaro)>tol;
end

sqrcd=kappa./log((10)./z0);
cd=sqrcd.^2;

u10=ustarn./sqrcd;