function psi = psim(zu,L)
% Integrated stablity term for momentum  
% zu=height of wind measurement (m) 
% L=Obukhov length scale (m)
% 
% Unstable case based on Paulson(1970)
% Stable case based on Dyer (1974)
% beta = 4.7;    (Bussinger et al) 
  beta = 5;   %  (e.g. Smith, 1988)
% beta = 7;   %  (e.g. Large and Pond, 1982) 
 

% stable case
  psi =-beta*(zu./L);

% unstable case
  jf = find(zu./L < 0); 
    if ~isempty(jf);
    if length(zu) > 1; 
      x = (1-16.*(zu(jf)./L(jf))).^0.25;
    else  
      x = (1-16.*(zu./L(jf))).^0.25;
    end  
    psi(jf) = 2.*log((1+x)/2)+log((1+x.^2)/2)-2.*atan(x)+3.14159/2;     
   end;


