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;

