function qlw=lwhf(Ts,dlw,dsw)
% LWHF: computes net longwave heat flux following Dickey et al (1994).
% qlw=LWHF(Ts,dlw) computes the net longwave heat flux into the ocean.
% Following Dickey et al (1994), J. Atmos. Oceanic Tech., 11, 1057-1078,
% the incident longwave flux can be corrected for sensor heating due to
% insolation if you are using Epply or Kipp & Zonen CG1 pyrgeometers.
% In this case, use qlw=LWHF(Ts,dlw,dsw). Epply is the default 
% pyrgeometer; change code for the Kipp & Zonen instrument.
%
%   INPUT:  Ts  - sea surface temperature [C]
%           dlw - (measured) downward longwave flux [W/m^2]
%           dsw - (measured) insolation [W/m^2] (needed for Eppley 
%                 or Kipp & Zonen pyrgeometers)
%
%   OUTPUT: qlw - net longwave heat flux [W/m^2]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/19/98: version 1.1 (revised for non-Epply pyrgeometers by RP)
% 4/9/99: version 1.2 (included Kipp & Zonen CG1 pyrgeometers by AA)
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% get constants
as_consts;

% convert degC to degK
ts=Ts+CtoK;

% correct dlw for sensor heating by insolation
if nargin==3,
  % this line is for Epply pyrgeometers
  dlwc=dlw-0.036.*dsw; 

  % this line is for Kipp & Zonen CG1 pyrgeometers  
  % (the offset is specified as 25 W/m^2 at 1000 W/m^2)
  % dlwc=dlw-0.025.*dsw;
else
  dlwc=dlw;
end;

% compute upward gray-body longwave flux
lwup=-emiss_lw.*sigmaSB.*(ts.^4);

% compute net flux
qlw=lwup + emiss_lw.*dlwc;