function  [shf,lhf,tau] = sfcflux_hf_robust(utrue, tair, rh, tsfcsea, p, z);
% function that calculates surface fluxes and other variables
% 1) Calculates neutral fluxes for all data
% 2) Calculates fluxes for subset of data that are likely to converge
% via the stars routine
% inputs Mx1 vectors of wind, tair, rh or q, tsfcsea & p ,z= measmt. heights
% outputs only shf, lhf, tau
% used for open water or open lead calculations
% based on ustar, tstar, qstar and L 
% Peter Guest 4/3/97 and Scott Tessmer
% modified by Ian Renfrew May 98, Feb 99, august 2002
% error found and corrected august 2003
%
%  Input:
%          utrue   wind speed (m/s) at z(1)
%          tair    air temperature (Centigrade) at z(2)
%          rh      relative humidity (%) at z(3)
%          q       specific humidity (g/kg) at z(3)
%          tsfcsea    surface temperature (Centigrade)
%          p       atmospheric pressure at z(4) (mb)
%          z      measurement level
% a check is added if rh < 1 then it is assumed this variable is q
%
%  stars output:   
%        ustar   friction velocity (m/s)
%        tstar   scaling temperature (k or c)
%        qstar   scaling specific humidity (g/kg)  not (g/g)!!!!
%        L       monin-obukov length scale (m)
%
%  sfcflux output:
% 
%        shf    sensible heat flux
%        lhf    latent heat flux
%        hf     (total) heat flux
%        wstar   free convection scaling length (needs zi)
%        tau     wind stress
%        Cd      drag coeff
%        Ce      Dalton number
%        Cdn     neutral drag coeff
%        Ch      Stanton number
%  plus a bunch of other stuff

%     lv latent heat of vaporization at 2.50e+6 J kg-1
%     cp    specific heat set constant set at 1004 J kg-1 k-1   

% set up 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)
CDn10=1.0e-3;  	% momentum flux transfer parameter at z10=10 (Smith, 1988)
gamma=0.00975; 	% adiabatic lapse rate
k=0.4;         	% Von Karmen's constant 
g=9.81;        	% gravity
Cp = 1004;	% specific heat capacity at constant pressure
Lv = 2.5e6;	% latent heat of vaporization 

% set measurement heights
if (length(z) == 4)
 zu = z(1);
 zt = z(2);
 zq = z(3);
 zp = z(4);
else
 disp('measurement height array is wrong size')
end;
    
% Calculate mixing ratios    
psfc = (p + 0.116*zp);   	% surface pressure (mb) based on standard atms
ess = esat(tsfcsea);         	% saturation vapor pressure wrt flat water
qsat = 622.0.*ess./(psfc - ess);% saturation mixing ratio (g/kg)
qsfc = qsat .* 0.98;       	% for salt water  (different from Smith) 
esa = esat(tair);         	% saturation vapor pressure wrt flat water
q = 622.0.*(rh/100).*esa./( p - esa); 
 
% potential air temperature (K) measured at zt
theta=tair+273.16+gamma.*zt;   
thsfc=tsfcsea+273.16;  
 
% virtual potential temp (K) based on Stull (1988)
thetav=theta.*(1.0+0.61e-3*q);  	% virtual potential air temp
thetavsfc=thsfc.*(1.0+0.61e-3*qsfc); 	% virt pot temp at sfc
rho=psfc./(2.87.*thetavsfc); 		% calc density with ideal gas law 

%-----------------------------------------------------------------------------
% As a first guess:
% calculate surface fluxes using constant exchange coeffs (neutral conditions)
%-----------------------------------------------------------------------------
utrue_save = utrue;
u_min = 0.5;
utrue(find(utrue<u_min)) = u_min;

shf = CHn10*Cp*rho.*utrue.*(tsfcsea-tair);
lhf = CEn10*Lv*1e-3*rho.*utrue.*(qsfc-q);
tau = CDn10*rho.*utrue.*utrue;

% Calculate ustar, tstar, qstar, L
ustar = sqrt(tau./rho);
tstar = -shf./(Cp*rho.*ustar);
qstar = -lhf./(Lv*1e-3*rho.*ustar); 		% qstar is g/kg
tstarv=tstar+(0.61e-3.*theta.*qstar);

%-----------------------------------------------------------------------------
% Recalculate surface fluxes using the stars subroutine (i.e. stability dep.)
% for data that are unstable and windy or stable and windy
%-----------------------------------------------------------------------------
e = find(((tair<tsfcsea) & (utrue>2)) | ((tair>tsfcsea) & (utrue>2)));
disp(['stars subroutine used for ' num2str(length(e)) ' data, from ' num2str(length(tair))])
[ustare tstare qstare Le] = stars_rh(utrue(e),tair(e),tsfcsea(e),rh(e),p(e),z);

% replace points that are convergent only, i.e. a subset of e called e2
% but note ustare(e1)
% error found in this replacement on 15 august 2003 IAR
e1 = find(isnan(ustare)==0);		
disp(['stars subroutine convergences for ' num2str(length(e1)) ' data'])
e2 = e(e1);						
shf(e2)= -rho(e2).*Cp.*ustare(e1).*tstare(e1); 	% sensible heat flux  
lhf(e2)= -rho(e2).*Lv*1e-3.*ustare(e1).*qstare(e1); % latent heat flux 
tau(e2)= rho(e2).*(ustare(e1).^2);   		% wind stress

ustar(e2) = ustare(e1);
tstar(e2) = tstare(e1);
qstar(e2) = qstare(e1);
tstarv=tstar+(0.61e-3.*theta.*qstar);