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;