cdntc

 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]

0001 function [cd,u10]=cdntc(sp,z,Ta)
0002 % CTDTC: computes the neutral drag coefficient following Smith (1988).
0003 % [cd,u10]=CDNTC(sp,z,Ta) computes the neutral drag coefficient and
0004 % wind speed at 10m given the wind speed and air temperature at height z
0005 % following Smith (1988), J. Geophys. Res., 93, 311-326.
0006 %
0007 % INPUT:   sp - wind speed  [m/s]
0008 %          z - measurement height [m]
0009 %          Ta - air temperature (optional)  [C]
0010 %
0011 % OUTPUT:  cd - neutral drag coefficient at 10m
0012 %          u10 - wind speed at 10m  [m/s]
0013 
0014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0015 % 3/8/97: version 1.0
0016 % 8/26/98: version 1.1 (vectorized by RP)
0017 % 8/5/99: version 2.0
0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019 
0020 % get constants
0021 as_consts; 
0022 
0023 if nargin==2,
0024   Ta=Ta_default;
0025 end;
0026 
0027 % iteration endpoint
0028 tol=.00001;
0029 
0030 visc=viscair(Ta);
0031 
0032 % remove any sp==0 to prevent division by zero
0033 i=find(sp==0);
0034 sp(i)=.1.*ones(length(i),1);
0035 
0036 % initial guess
0037 ustaro=zeros(size(sp));
0038 ustarn=.036.*sp;
0039 
0040 % iterate to find z0 and ustar
0041 ii=abs(ustarn-ustaro)>tol;
0042 while any(ii(:)),
0043   ustaro=ustarn;
0044   z0=Charnock_alpha.*ustaro.^2./g + R_roughness*visc./ustaro;
0045   ustarn=sp.*(kappa./log(z./z0));
0046   ii=abs(ustarn-ustaro)>tol;
0047 end
0048 
0049 sqrcd=kappa./log((10)./z0);
0050 cd=sqrcd.^2;
0051 
0052 u10=ustarn./sqrcd;