wstress

 WSTRESS  computes wind stress using the Large and Pond (JPO, 1981) formulation.

  USAGE: [taux,tauy,u10,v10]=wstress(u,v,z0)

       u,v = eastward, northward components of wind in m/s
       z0  = height of wind sensor in m (assumed to be 10 m if not supplied)
     
       taux,tauy = eastward, northward wind stress components (dynes/cm2)
       u10,v10   = equivalent wind velocity components at 10 m.

0001 function [taux,tauy,u10,v10]=wstress(u,v,z0)
0002 % WSTRESS  computes wind stress using the Large and Pond (JPO, 1981) formulation.
0003 %
0004 %  USAGE: [taux,tauy,u10,v10]=wstress(u,v,z0)
0005 %
0006 %       u,v = eastward, northward components of wind in m/s
0007 %       z0  = height of wind sensor in m (assumed to be 10 m if not supplied)
0008 %
0009 %       taux,tauy = eastward, northward wind stress components (dynes/cm2)
0010 %       u10,v10   = equivalent wind velocity components at 10 m.
0011 %
0012 
0013 % Rich Signell  rsignell@usgs.gov
0014 %
0015 % Version 1.1 2005-04-15
0016 % Bug fixed:  Line 49 was changed from : t=(rho*cd.*d1*1.e4)./(1+a1.*sqrt(cd));
0017 %                                   to : t=(rho*cd.*d1*1.e4)
0018 %             The effect of this bug was to yield slightly incorrect wind stresses for
0019 %             anemometer heights different from 10 m. For example, a 10 m wind at 5 m
0020 %             height yielded 1.806 dyn/cm2 instead of 1.699 dyn/cm2 (6.3% too high), while a
0021 %             10 m wind at 30 m height yielded 1.154 dyn/cm2 instead of 1.261 dyn/cm2 (8.5% too low)
0022 %             Only the wind stresses were affected.  The output 10m winds were correct.
0023 %             Thanks to Charlie Stock for finding this bug.
0024 
0025 if(nargin==2),
0026  z0=10.;
0027 end;
0028 %
0029 nans=ones(length(u),1)*nan;
0030 u10=nans; 
0031 v10=nans;
0032 taux=nans;
0033 tauy=nans;
0034 %
0035 % identify times of zero wind.  These indices will be set to zero stress.
0036 izeros=find((abs(u)+abs(v))==0);
0037 igood=find(finite(u)& ((abs(u)+abs(v))>0));
0038 
0039 u=u(igood);
0040 v=v(igood);
0041 u=u(:);
0042 v=v(:);
0043 w=u+i*v;
0044 v0=abs(w);
0045 k=.41;
0046 rho=1.25e-3;
0047 a1=1/k*log(z0/10.);
0048 d1=1.e35*ones(size(u));
0049 c=zeros(size(u));
0050 while (max(abs(c-d1))>.01),
0051   c=d1;
0052   cd=ones(size(u))*1.205e-3;
0053   ind=find(c>11);
0054   if(~isempty(ind))
0055    cd(ind)=(0.49 + 0.065*c(ind))*1.e-3;
0056   end
0057   d1=v0./(1+a1.*sqrt(cd));
0058 end
0059 %t=(rho*cd.*d1*1.e4)./(1+a1.*sqrt(cd));  % wrong
0060 t=(rho*cd.*d1*1.e4);                     % fix for version 1.1
0061 w10=(d1./v0).*w;
0062 u10(igood)=real(w10);
0063 v10(igood)=imag(w10);
0064 taux(igood)=t.*u10(igood);
0065 tauy(igood)=t.*v10(igood);
0066 % set zero wind periods to zero stress (and U10,V10)
0067 u10(izeros)=zeros(size(izeros));
0068 v10(izeros)=zeros(size(izeros));
0069 taux(izeros)=zeros(size(izeros));
0070 tauy(izeros)=zeros(size(izeros));