function [taux,tauy]=stress10(u,v,sensor_height)
%STRESS10 Compute windstress at 10-meter reference level 
% STRESS10 computes the windstress from the windspeed (u,v)
% recorded at sensor_height meters above MSL, as described
% in Large and Pond,  1981, JPO, and Blanton, et. al., GALE...
%
% Instead of iterating until some tolerance is satisfied,
% 4 iterations are made, PERIOD.
%
% STRESS10 assumes neutral stability.  The routine
% STAB10 makes corrections to the neutral stability
% stress estimates from air-sea temperature differences
% and humidity.
%
% Inputs: u,v  -  east/north wind velocity (meters/sec)
%   	   sensor_height - anemometer height (meters)
%
% Outputs:[taux,tauy] - windstress (pascals) at 10 meter ref level.
%
% Call as:[taux,tauy]=stress10(u,v,sensor_height)
%
% Written by: Brian O. Blanton
% Spring '98
%
 
% 3 args in required.
if nargin~=3
   error('3 args to STRESS10 are REQUIRED!')
end

% Set number of iterations
nit=4;

% Von Karman's Constant
kar=.4;

% density of air, kg/m^3
rho=1.22; 

% Compute magnitude
mag=sqrt(u.*u+v.*v);

% We'll need to keep track of NaN and 0 speeds entries.
ikeep=find(~isnan(mag)& mag>eps);
iNaN=find(isnan(mag));
izero=find(mag<eps);

% Make temp arrays of u,v EXCLUDING NaN entries.
uu=u(ikeep);
vv=v(ikeep);

% Compute mag and direction
mag=sqrt(uu.*uu+vv.*vv);
dir=atan2(vv,uu);

% Initial Z0 and u10 calculations
Z0=5e-4*ones(size(uu));
u10=mag.*log(10./Z0)./log(sensor_height./Z0);

for it=1:nit
   Cdn=(.49+.065*u10)*1e-3;
   ustar=sqrt(Cdn.*(u10.*u10));
   ln_z0=log(sensor_height)-kar*mag./ustar;
   Z0=exp(ln_z0);
   u10=mag.*log(10./Z0)./log(sensor_height./Z0);
end
   
%%%%%%%%%% FINAL PART, IF U10<10m/s, 
idx=find(mag<10.15385);
u10(idx)=11.55*mag(idx)./(log(sensor_height./9.67e-5));

% Revert back to u,v from u10 and dir
x10=u10.*cos(dir);
y10=u10.*sin(dir);

%%%%%%%%%%
%%%%%%%%%% IT'S TIME TO CONVERT TO STRESS
%%%%%%%%%%
Cd=NaN*ones(size(u10));
idx=find(mag<10.15385);
Cd(idx)=1.15*1e-3 * ones(size(idx));
idx=find(isnan(Cd)==1);
Cd(idx)=(.49+.065*u10(idx))*1e-3;

tx=rho*Cd.*u10.*x10;
ty=rho*Cd.*u10.*y10;

% Insert stresses into vectors of original length.
% 0-speeds are obviously 0-stress
taux=NaN*ones(size(u));
tauy=NaN*ones(size(u));

taux(ikeep)=tx;
tauy(ikeep)=ty;

taux(izero)=0;
tauy(izero)=0;


%
%        Brian O. Blanton
%        Department of Marine Sciences
%        Ocean Processes Numerical Modeling Laboratory
%        12-7 Venable Hall
%        CB# 3300
%        University of North Carolina
%        Chapel Hill, NC
%                 27599-3300
%
%        919-962-4466
%        blanton@marine.unc.edu
%
%        Spring  1998
%