function lwest = blwhf(Ts,Ta,rh,Fc,bulk_f)
% BLWHF: estimates net long-wave heat flux using bulk formulas
% lwest = BLWHF(Ts,Ta,rh,Fc,bulk_f) estimates the net longwave 
% heat flux into the ocean using one of a number of bulk formulas 
% evaluated by Fung et al (1984), Rev. Geophys. Space Physics, 
% 22,177-193. 
%
% INPUT:   Ts - sea surface temperature [C]
%          Ta - air temperature  [C]
%          rh - relative humidity  [%]
%          Fc - cloudiness correction factor F(C) (=1 for clear sky) 
%                  (see cloudcor.m for more information)
%          bulk_f - bulk formulas to be used 
%                  (see Fung et al, 1984) for details). 
%                   Options are:
%                     'brunt'
%                     'berliand'  - probably the best one (default)
%                     'clark'
%                     'hastenrath'
%                     'efimova'
%                     'bunker'
%                     'anderson'
%                     'swinbank' - does not use rh
%
% OUTPUT:  lwest - net downward longwave flux [W/m^2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 8/31/98: version 1.1 (contributed by RP)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load constants and defaults
as_consts; 

if nargin<5,
 bulk_f=bulkf_default;
end;

% compute vapour pressure from relative humidity (using formulas 
% from Gill, 1982)
ew=satvap(Ta);
rw=eps_air*ew./(P_default-ew);
r=(rh/100).*rw;
e_a=r*P_default./(eps_air+r);

% convert to K
ta=Ta+CtoK;
ts=Ts+CtoK;

% signs for flux INTO ocean
switch bulk_f,
  case {'brunt',1}
    lwest =  -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc;

  case {'berliand',2}
    lwest =  -emiss_lw*sigmaSB.*ta.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc ...
             - 4*emiss_lw*sigmaSB.*ta.^3.*(Ts - Ta);

  case {'clark',3}
    lwest =  -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.05*sqrt(e_a)).*Fc ...
             - 4*emiss_lw*sigmaSB.*ts.^3.*(Ts - Ta);

  case {'hastenrath',4}
    q_a = e_a*eps_air./P_default;
    lwest =  -emiss_lw*sigmaSB.*ts.^4.*(0.39 - 0.056*sqrt(q_a)).*Fc ...
             - 4*emiss_lw*sigmaSB.*ts.^3.*(Ts - Ta);

  case {'efimova',5}
    lwest =  -emiss_lw*sigmaSB.*ta.^4.*(0.254 - 0.00495*(e_a)).*Fc ;
    
  case {'bunker',6}
    lwest =  -0.022*emiss_lw*sigmaSB.*ta.^4.*(11.7 - 0.23*(e_a)).*Fc ...
             - 4*emiss_lw*sigmaSB.*ta.^3.*(Ts - Ta);

  case {'anderson',5}
    lwest =  -emiss_lw*sigmaSB.*(ts.^4 - ta.^4.*(0.74+0.0049*(e_a))).*Fc ;
 
  case {'swinbank',5}
    lwest =  -emiss_lw*sigmaSB.*(ts.^4 - 9.36e-6*ta.^6).*Fc ;

   otherwise
    error(['Unrecognized bulk formula specified: ' bulk_f]);
end;