  function tau=stresslp(sp,z)
% STRESSLP: computes neutral wind stress following Large and Pond (1981).
% tau = STRESSLP(sp,z) computes the neutral wind stress given the wind
% speed at height z following Large and Pond (1981), J. Phys. Oceanog.,
% 11, 324-336. Assumes z a fixed scalar, and constant air density
% (1.22 kg/m^3). 
%
%   INPUT:   sp - wind speed (m/s)
%            z - measurement height (m)
%
%   OUTPUT:  tau - wind stress (N/m^2)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set constants
k=.4;
a=log(z./10)./k;
rho=1.22;

[N,M]=size(sp);
if N==1
   sp=sp';
end

L=length(sp);
cd=zeros(L,1);
tau=zeros(L,1);
sp10=zeros(L,1);

for m=1:L;
% iterate to find wind speed at 10m
spold=0.;
  for j=1:10;
     if spold < 10.15385
       c=.001*1.15;
     else
       c=.001*(0.49+0.065*spold);
     end
    spnew=sp(m)./(1.+a*sqrt(c));
     if (abs(spnew-spold)) < .001
       break
     else
       spold=spnew;
     end
    sp10(m)=spnew;
  end
end

% compute cd and stress using 10m sp
cd=.001.*1.15.*ones(L,1);
ind=find(sp10>=10.15385);
cd(ind)=.001*(0.49+0.065*sp10(ind));
tau=rho*cd.*(sp10.^2);
if N==1
  tau=tau';
end
end

