function [ustar, tstar, qstar, L] = stars_rh(utrue, tair, tsfc, rh, p, z);
% This function determines turbulent scaling parameters based on bulk inputs
% Peter Guest   5/8/95
% Based on "asiem.m" program written by Scott Tessmer
% and used for course cruise.
% Derived from earlier versions by Joe Johnson and me
% This version attempts interations for all winds > 0.5 m/s, 
% although theory is questionable below 3 m/s
% Note suggested humidity changes 11/96 
% Changed 9/8/97 so that Qsfc = 98% saturation instead of 100%.
% NOTE new CHn, CEn constants from Decosmo et al
% 
% input parameters

% utrue      wind speed (m/s)            required  
% tair       air temperature (C)         required for tstar
% tsfc       surface temperature (C)     required for tstar
% rh         relative humidity (%)       recommended  (default = 80 %)
%% q          specific humidity (g/kg)
% p          pressure (mb)               optional   default = 1012 mb      
% z(1)       height of wind speed (m)    optional   default = 10 m
% z(2)       height of air temp (m)      optional   default = z(1)
% z(3)       height of humidity (m)      optional   default = z(2)
% z(4)       height of pressure (m)      optional   default = z(2)
%
% The first five input parameters must be M x 1 or 1 x 1 arrays, where 
% M represents the number of rows (i.e observations).
% A 1 X 1 array represents a constant value for all observations and can be
% used everywhere except utrue (or all can be  1 X 1) 
 
% output parameters 

% ustar      friction velocity (m/s)
% tstar      temperature scaling parameter (C or K)
% qstar      humidity scaling parameter (g/kg)
% L          Obukhov length scale
%
% The output parameters are M x 1 arrays.


% set constants
z10=10.0;      % reference height for Obukhov length scale 
CHn10=1.14e-3;  % heat buoyancy flux transfer parameter at z10=10 (Smith 88, DeCosmo 1996)
CEn10=1.20e-3;  % humidity flux transfer parameter at z10=10 (Smith, 1988,DeCosmo 1996, Fairall et al 2003)
k=0.4;         % Von Karmen's constant 
g=9.81;        % gravity
a=0.011;       % Charnocks's constant from Smith(1988)
b=0.11;        % smooth flow constant from Businger(1973)
nu=1.461e-5;   % dynamic viscosity of air (Smith uses 1.4 e-5) 
gamma=0.00975; % adiabatic lapse rate

% Assign defaults if necessary

if isempty(rh) == 1  % Default humidity
  rh = 80;
end

if isempty(p) == 1  % Default pressure 
  p = 1012;
end
pbad = find(isnan(p));
p(pbad) = 1012  * ones(size(pbad)); 

if length(z) == 0 % Default measurement heights
  zu = 10.;
  zt = 10.; 
  zq = 10.; 
  zp = 10.;
elseif length(z) == 1
  zu = z(1);
  zt = z(1);
  zq = z(1);
  zp = z(1);
elseif length(z) == 2
  zu = z(1);
  zt = z(2);
  zq = z(2);
  zp = z(2);      
elseif length(z) == 3
  zu = z(1);
  zt = z(2);
  zq = z(3);
  zp = z(2);      
elseif length(z) == 4
  zu = z(1);
  zt = z(2);
  zq = z(3);
  zp = z(4);
else
  disp ('Incorrect z array')
end 

% Calculate mixing ratios    
%NOTE q is now an input variable
psfc = (p + 0.116.*zp);   	% surface pressure (mb) based on standard atms
ess = esat(tsfc);         	% saturation vapor pressure wrt flat water
%ess = esati(tsfc);         	%                               ice  
qsat = 622.0.*ess./(psfc - ess);  % saturation mixing ratio (g/kg)
qsfc = qsat;             	% mixing ratio at surface assumed staturated
%qsfc = 0.98.*qsat;          	% mixing ratio at surface assumed 98%
esa = esat(tair);         	% saturation vapor pressure wrt flat water
q = 622.0.*(rh./100).*esa./( p - esa); % mixing ratio at zt  (g/kg)

% potential air temperature (K) measured at zt
 theta=tair+273.16+gamma.*zt;   %

% virtual potential temp (K) based on Stull (1988)
 thetav=theta.*(1.0+0.61e-3*q);  % virtual potential air temp
 thetavsfc=(tsfc+273.16).*(1.0+0.61e-3*qsfc); % virt pot temp at sfc    

% initialization for iteration loop ustar... (these values will change)
  checker=1;
  ctr=0;
  L=1e10*ones(size(utrue));
  zo=1e-4*ones(size(utrue));
  ustar=0.035*ones(size(utrue));  % first guess 
  nanm=NaN*ones(size(utrue));   % NAN array to fill in NAN points

% Values below 2 m/s are invalid for these formula
% find indices of winds below threshold then replaces them with nan
  utruesave = utrue;
  ool=find(utrue<2);
  disp(['no. of data below wind threshold ' num2str(length(ool))])
   if ~isempty(ool) == 1;
     utrue(ool)=nanm(ool); 
   end;

% Start of interation loop             
% exit if all ustars converge within 0.05% or after 40 iterations 
   
while (checker>0.0005) & (ctr<40);  
 ctr=ctr+1;     %;

  ustar1=ustar;  % save old ustar for next check 
 
  zot=z10*exp(-k*k./(CHn10*log(z10./zo))); 	% Using constant CHn10 
  zoq=z10*exp(-k*k./(CEn10*log(z10./zo))); 	% Using constant CEn10 

  ustar=utrue*k./(log(zu./zo)-psim(zu,L));
  tstar=(theta-(tsfc+273.16))*k./(log(zt./zot)-psit(zt,L));
  qstar=(q-qsfc)*k./(log(zq./zoq)-psit(zq,L));
  tstarv=tstar+(0.61e-3.*theta.*qstar);
 
  zc=a.*ustar.*ustar./g;                 	% charnocks formula        
  zs=b.*nu./ustar;                       	% smooth formula
  zo=zc+zs;                             	% Smith's (1988) zo 
 
  okl = find(abs(tstarv) >= 1.0e-30); 	% to find L: find ok data 
  bigl = find(abs(tstarv) < 1.0e-30); 	% prevent underflow
  zerol = find(abs(tstarv) == 0);     	% prevent division by zero
  L(bigl) = 1.0e20.*ones(size(bigl)).*sign(tstarv(bigl));
  L(zerol) = 1.0e20.*ones(size(zerol));  
  L(okl)=(ustar(okl).^2).*thetav(okl)./(g*k*tstarv(okl));     % most refs
  L(find((L<z10)&(L>0))) = z10;		% set small stable L to 1 (IAR aug02)

  checker=isnan(ustar); % find NAN points to avoid checking
  j=find(checker~=1);
  checker=max(abs((ustar(j)-ustar1(j))./ustar(j))); % convergence check

end % while iteration loop

 % finding indices where ustar does not converge    
 noncon=find(abs((abs(ustar)-abs(ustar1))./ustar)>0.0005);
 lnoncon = length(noncon);
 if (lnoncon > 0) 
  disp('non-converging:')
  disp(lnoncon)
  disp('iterations')
  disp(ctr)
 end;
 % setting ustar, tstar, qstar, & L indices that are not valid 
 X=[ool; noncon]; 	% May have to change dimension of entering matrix

 % choice of setting values to be small or NAN
  if ~isempty(X) == 1
%     utrue = utruesave; 
%     ustar(X)=0.001 .* utrue(X);  
%     tstar(X)=0.001 .* (theta(X) - tsfc(X) - 273.16);
%     qstar(X)=0.001 .* (q(X) - qsfc(X));
%     tstarv(X)=tstar(X)+(0.61e-3.*theta(X).*qstar(X));
%     L(X)= (ustar(X).^2).*thetavsfc(X)./(g*k*tstarv(X));   % 
     ustar(X) = nanm(X);
     tstar(X) = nanm(X);
     qstar(X) = nanm(X);
     L(X) = nanm(X);
  end;





