function psi=psit(zt,L)
% Integrated stablity term for temperature  
% zt=height of temperature 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, 1971) 
  beta = 5;   %  (e.g. Smith, 1988; Dyer, 1974)
% beta = 7;   %  (e.g. Large and Pond, 1982) 
 

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

