function A=hfbulktc(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw)
% HFBULKTC: computes sensible and latent heat fluxes and other variables.
% A=HFBULKTC(ur,zr,Ta,zt,rh,zq,Pa,Ts,sal,dlw,dsw,nsw) computes the following:
%
%         Hs      = sensible heat flux INTO ocean [W/m^2]
%         Hl      = latent heat flux INTO ocean [W/m^2]
%         Hl_webb = Webb correction to latent heat flux INTO ocean [W/m^2]
%         stress  = wind stress [N/m^2]
%         U_star  = velocity friction scale [m/s]
%         T_star  = temperature scale [deg C]
%         Q_star  = humidity scale [kg/kg]
%         L       = Monin-Obukhov length [m]
%         zetu    = zr/L
%         CD      = drag coefficient
%         CT      = temperature transfer coefficient (Stanton number)
%         CQ      = moisture transfer coefficient (Dalton number)
%         RI      = bulk Richardson number
%         Dter    = cool-skin temperature difference (optional output) [C]; 
%                   positive if surface is cooler than bulk (presently no 
%                   warm skin permitted by model)
%                    
% Based on the following buoy input data:
%
%           ur     = wind speed [m/s] measured at height zr [m] 
%           Ta     = air temperature [C] measured at height zt [m]
%           rh     = relative humidity [%] measured at height zq [m]
%           Pa     = air pressure [mb]
%           Ts     = sea surface temperature [C]
%           sal    = salinity [psu (PSS-78)]
%                    (optional - only needed for cool-skin)
%           dlw    = downwelling (INTO water) longwave radiation [W/m^2]
%                    (optional - only needed for cool-skin)
%           dsw    = measured insolation [W/m^2]
%                    (optional - only needed for cool-skin)
%           nsw    = net shortwave radiation INTO the water [W/m^2]
%                    (optional - only needed for cool-skin)
%
% where ur, Ta, rh, Pa, Ts, zr, zt, and zq (and optional sal, dlw,
% dsw, and nsw) may be either row or column vectors; and rh, Pa, 
% zr, zt, and zq (and optional sal) may also be fixed scalars.
%
% Output variables are given as column vectors in A:
%
% 1) without cool-skin correction:
%
%   A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI]
%
% 2) with cool-skin correction: 
%
%   A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter];

% Code follows Edson and Fairall TOGA COARE code (version 2.0), modified 
% to include Rogers' weighting factor for unstable conditions.  Code does
% include gustiness, and assumes that the marine boundary layer height is
% known and constant over time for simiplicity. zr/L is limited to 
% be <=3.0 to ensure that the code converges to nonzero stress and heat 
% flux values for strongly stable conditions.  The bulk Richardson number
% is computed between the sea surface and zr as a diagnostic about whether
% turbulent boundary layer theory is applicable.  Code does not include 
% warm layer effects to modify Ts.  See Fairall et al (1996), J. Geophys. 
% Res., 101, 3747-3764, for description of full TOGA COARE code and 
% comparison with data. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/19/98: version 1.1 (rewritten by RP to remove inconsistencies in 
%          virtual and real temperatures, improve loop structure, 
%          correct gustiness component of stress computation) 
% 4/9/99: version 1.2 (rewritten by AA to clarify some variable names
%         and include cool-skin effect and Webb correction to latent 
%         heat flux added to output matrix)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=length(ur);

% change to column vectors
ur = ur(:);
Ta = Ta(:);
rh = rh(:);
Ts = Ts(:);
Pa = Pa(:);
zr = zr(:);
zt = zt(:);
zq = zq(:);

% create vectors for rh, Pa, zr, zt, and zq, if scalars are input
if length(rh)==1 & M>1
  rh=rh*ones(M,1);
end
if length(Pa)==1 & M>1
  Pa=Pa*ones(M,1);
end
if length(zr)==1 & M>1
  zr=zr*ones(M,1);
end
if length(zt)==1 & M>1
  zt=zt*ones(M,1);
end
if length(zq)==1 & M>1
  zq=zq*ones(M,1);
end

if nargin > 8
  % optional cool-skin stuff
  sal = sal(:);
  dlw = dlw(:);
  dsw = dsw(:);
  nsw = nsw(:);
  % create vector for sal if scalar is input
  if length(sal)==1 & M>1
    sal=sal*ones(M,1);
  end
end

% initialize various constants
as_consts;

tol=.001;    % tolerance on Re changes to make sure soln has converged.

onethird=1./3;
o61=1/eps_air-1;   % 0.61 (moisture correction for temperature)

visc=viscair(Ta);                 % viscosity
Qsats=Qsat_coeff*qsat(Ts,Pa);     % saturation specific humidity; the Qsat_coeff
                                  % value is set in routine as_consts.m
Q=(0.01.*rh).*qsat(Ta,Pa);        % specific humidity of air [kg/kg]
T =Ta+CtoK;                       % convert to K
Tv=T.*(1 + o61*Q);                % air virtual temperature
rho=(100*Pa)./(gas_const_R*Tv);   % air density
Dt=(Ta+0.0098.*zt)-Ts;            % adiabatic temperature difference
Dq=Q-Qsats;                       % humidity difference

% compute initial neutral scaling coefficients
S=sqrt(ur.^2 + min_gustiness.^2);
cdnhf=sqrt(cdntc(S,zr,Ta)); % Smith's neutral cd as first guess

z0t=7.5*10^(-5);
ctnhf=kappa./log(zt./z0t);

z0q=z0t;
cqnhf=kappa./log(zq./z0q);

U_star = cdnhf.*S;      % (includes gustiness)
T_star = ctnhf.*Dt;     % 
Q_star = cqnhf.*Dq;     %

Dter   = 0;
Dqer   = 0;
if nargin > 8
% initial cool-skin thickness guess  
  delta = 0.001*ones(size(Ts));
end

Reu=0;Ret=0;Req=0;
% begin iteration loop to compute best U_star, T_star, and Q_star
for iter1=1:80;

    ReuO=Reu; RetO=Ret; ReqO=Req; % Save old values
    
    % Compute Monin-Obukov length (NB - definition given as eqn (7)
    % of Fairall et al (1996) probably wrong, following, e.g.
    % Godfrey and Bellars (1991), JGR, 96, 22043-22048 and original code)
    bs=g*(T_star.*(1 + o61*Q) + o61*T.*Q_star)./Tv; 
    L=(U_star.^2)./(kappa*bs);
    % set upper limit on zr/L = 3.0 to force convergence under 
    % very stable conditions. Assume that zr, zt and zq comparable.
    index_limit   = (L<zr/3 & L>0);
    L(index_limit)=zr(index_limit)/3;
    
    zetu=zr./L;  % nondimensionalized heights
    zett=zt./L;
    zetq=zq./L;

    % surface roughness
    z0=(Charnock_alpha/g).*U_star.^2 + R_roughness.*visc./U_star;

    % compute U_star correction for non-neutral conditions
    cdnhf=kappa./(log(zr./z0)-psiutc(zetu));
    U_star=cdnhf.*S;
  
    Reu=z0.*U_star./visc;   % roughness Reynolds #
    [Ret,Req]=LKB(Reu);  % compute other roughness Reynolds #s

    % compute t and q roughness scales from roughness R#s
    z0t=visc.*Ret./U_star;
    z0q=visc.*Req./U_star;

    % compute new transfer coefficients at measurement heights
    cthf=kappa./(log(zt./z0t)-psittc(zett));
    cqhf=kappa./(log(zq./z0q)-psittc(zetq));

    % compute new values of T_star, Q_star
    T_star=cthf.*(Dt + Dter);
    Q_star=cqhf.*(Dq + Dqer);

    % estimate new gustiness
    Ws=U_star.*(-CVB_depth./(kappa*L)).^onethird;
    wg=min_gustiness*ones(M,1);
    j=find(zetu<0);                 % convection in unstable conditions only
    wg(j)=max(min_gustiness,beta_conv.*Ws(j)); % set minimum gustiness
    S=sqrt(ur.^2 + wg.^2);

    if nargin > 8
    % compute cool-skin parameters
      [delta,Dter,Dqer] = cool_skin(sal,Ts-Dter,rho,cp,Pa, ...
                                    U_star,T_star,Q_star, ...
                                    dlw,dsw,nsw,delta,g,gas_const_R, ...
                                    CtoK,Qsat_coeff);
    end

end % end of iteration loop

ii= abs(Reu-ReuO)>tol | abs(Ret-RetO)>tol | abs(Req-ReqO)>tol;
if any(ii),
 disp(['Algorithm did not converge for ' int2str(sum(ii)) ' values. Indices are:']);
 disp(find(ii)');
 warning('Not converged!');
end;


% compute latent heat
Le=(2.501-0.00237*(Ts-Dter))*10^6;

% compute fluxes into ocean
Hs=rho.*cp.*U_star.*T_star;
Hl=rho.*Le.*U_star.*Q_star;

% compute transfer coefficients at measurement heights
CD=(U_star./S).^2;
CT=U_star.*T_star./(S.*(Dt + Dter)); % Stanton number
CQ=U_star.*Q_star./(S.*(Dq + Dqer)); % Dalton number

% to compute mean stress, we don't want to include the effects
% of gustiness which average out (in a vector sense).
stress=rho.*CD.*S.*ur;

% compute bulk Richardson number (as a diagnostic) - the "T"
% is probably not quite right - assumes T \ approx. Ts (good enough though)
RI=g.*zr.*((Dt + Dter) + o61*T.*(Dq + Dqer))./(Tv.*S.^2);

% compute Webb correction to latent heat flux into ocean
W = 1.61.*U_star.*Q_star + (1 + 1.61.*Q).*U_star.*T_star./T; % eqn. 21
Hl_webb = rho.*Le.*W.*Q; % eqn. 22, Fairall et al. (1996), JGR, 101, p3751.

% output array
if nargin > 8
  % output additional cool-skin parameters 
  A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI Dter];
else
  % otherwise
  A=[Hs Hl Hl_webb stress U_star T_star Q_star L zetu CD CT CQ RI];
end



function y=psiutc(zet)
% PSIUTC: computes velocity profile function following TOGA/COARE.
% y=PSIUTC(zet) computes the turbulent velocity profile function given 
% zet = (z/L), L the Monin-Obukoff length scale, following Edson and
% Fairall TOGA COARE code (version 2.0) as modified to include Rogers' 
% weighting factor to combine the Dyer and free convection forms for 
% unstable conditions. 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/28/98: version 1.1
% 8/5/99: version 1.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c13=1.0./3.0;
sq3=sqrt(3.0);

% stable conditions
y=-4.7*zet;

% unstable conditions
j=find(zet<0);
zneg=zet(j);

% nearly stable (standard functions)
 x=(1-16.0.*zneg).^0.25;
 y1=2.0.*log((1+x)./2) + log((1+x.^2)./2) -2.*atan(x) + pi/2;

% free convective limit
 x=(1-12.87*zneg).^c13;
 y2=1.5*log((x.^2+x+1)./3) - sq3*atan((2.*x+1)/sq3) + pi/sq3;
	
% weighted sum of the two
 F=1.0./(1+zneg.^2);
 y(j)=F.*y1+(1-F).*y2;




function y=psittc(zet)
% PSITTC: computes potential temperature profile following TOGA/COARE.
% y=PSITTC(zet) computes the turbulent potential temperature profile 
% function given zet = (z/L), L the Monin-Obukoff length scale, following 
% Edson and Fairall TOGA COARE code (version 2.0), as modified to use
% Rogers' weighting factor to combine the Dyer and free convective 
% forms for unstable conditions. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/28/98: version 1.1
% 8/5/99: version 1.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

c13=1.0./3.0;
sq3=sqrt(3.0);

% stable conditions
y=-4.7.*zet;

% unstable conditions
j=find(zet<0);
zneg=zet(j);

% nearly stable (standard functions)
 x=(1-16.0.*zneg).^0.25;
 y1=2.0*log((1+x.^2)./2);

% free convective limit
 x=(1-12.87*zneg).^c13;
 y2=1.5.*log((x.^2+x+1)./3.0) - sq3.*atan((2.*x+1)./sq3) + pi./sq3;

% weighted sum of the two
 F=1.0./(1+zneg.^2);
 y(j)=F.*y1 + (1-F).*y2;





function [Ret,Req]=LKB(Reu);
% LKB: computes rougness Reynolds numbers for temperature and humidity
% [Ret,Req]=LKB(Reu) computes the roughness Reynolds for temperature
% and humidity following Liu, Katsaros and Businger (1979), J. Atmos.
% Sci., 36, 1722-1735.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/28/98: version 1.1
% 8/5/99: version 1.2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   Ret=.177*ones(size(Reu));
   Req=.292*ones(size(Reu));
j=find(Reu>.11 & Reu<=.825);
   Ret(j)=1.376.*Reu(j).^0.929;
   Req(j)=1.808.*Reu(j).^0.826;
j=find(Reu>.825 & Reu<=3);
   Ret(j)=1.026./Reu(j).^0.599;
   Req(j)=1.393./Reu(j).^0.528;
j=find(Reu>3 & Reu<=10);
   Ret(j)=1.625./Reu(j).^1.018;
   Req(j)=1.956./Reu(j).^0.870;
j=find(Reu>10 & Reu<=30);
   Ret(j)=4.661./Reu(j).^1.475;
   Req(j)=4.994./Reu(j).^1.297;
j=find(Reu>30);
   Ret(j)=34.904./Reu(j).^2.067;
   Req(j)=30.790./Reu(j).^1.845;