function [taux,tauy]=wstr(windu,windv,hgt);
%
% WSTR  calculates neutral stability wind stress
% from Large and Pond algorithm (JPO, 1981)
%
% input wind velocity in m/s and anomometer height in m (hgt),
% get back Pascals (N/m**2)
%
% USAGE: 
% [taux,tauy]=wstr(windu,windv,hgt);  % east and north components
%
%
if hgt==10.0;
  al=0;
else
  al=(1/0.41)*log(hgt/10.0);
end;
%
wind=sqrt(windu^2+windv^2);
%  
% iteratively calculate wind speed at 10m
%
c=0.0;
dl=.03;
%
while abs(c-dl)>0.02;
  c=dl;
  if c<11.0;
    cd=0.001205;
  else;
    cd=(0.49+(.065*c))*0.001;
  end;
  dl=wind/(1.0 + al*sqrt(cd));
end;
%
% Now calculate final drag coef
%
  v10=dl;   
  if dl<11.0;
    cd=0.001205;
  else;
    cd=(0.49+(.065*dl))*0.001;
  end;

t=(.00125*cd*v10*1000)/(1.0 + al*sqrt(cd));
taux=t*windu;
tauy=t*windv;
wstress=sqrt((taux^2)+(tauy^2));
