function [qsw,alb]=swhf(yd,yr,long,lat,dsw)
% SWHF: computes net shortwave heat flux into the ocean and albedo.
% [qsw,alb]=SWHF(yd,yr,long,lat,dsw) computes the net shortwave heat 
% flux into the ocean and albedo, using Payne (1972), J. Atm. Sci., 29,
% 959-970, to estimate the instantaneous albedo given the atmospheric 
% transmittance (the ratio of measured insolation to the no-atmosphere 
% insolation).
%
% INPUT:  yd   - decimal yearday (e.g., 0000Z Jan 1 is 0.0)
%         yr   - year (e.g., 1995)
%         long - longitude  [deg]
%         lat  - latitude  [deg]
%         dsw  - (measured) insolation  [W/m^2]
%
% OUTPUT: qsw - net shortwave heat flux [W/m^2]
%         alb - albedo 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 7/29/99: version 1.1
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute sun altitude and no atm solar radiation 
[sunalt, sorad]=soradna1(yd,yr,long,lat);

% compute atm transmittance  (note: trans=Inf when sorad=0)
trans=inf.*ones(size(sorad));
j=find(sorad>0);
trans(j)=dsw(j)./sorad(j);

% compute albedo  (note: alb=NaN when trans>1 or sunalt<0)
alb=albedo(trans, sunalt);

% compute net shortwave heat flux
%      (note: qsw set equal to 0 when alb=NaN)
qsw=(1-alb).*dsw;
ind=find(isnan(qsw));
qsw(ind)=zeros(size(ind));