function c = reedcf(yd,lat,dsw);
% REEDCF: computes daily mean cloud cover following Reed (1977).
% c = REEDCF(yd,lat,dsw) computes daily averaged cloud cover c from
% yearday, latitude, and measured insolation following Reed (1977), J. Phys.
% Oceanog., 7, 482-485. Assumes hourly input series are either both
% column or both row vectors of equal length. c is output as a boxcar
% function with c constant over a 24 hr period.
%
% INPUT: yd - yearday (e.g., Jan 10 is yd=10)
% lat - latitude [deg]
% dsw - (measured) insolation [W/m^2]
%
% OUTPUT: c - daily averaged cloud cover
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine if input is either row or column vectors
[m,n] = size(dsw);
% first daily average swr and yd
davesw = binave(dsw,24);
daveyd = binave(yd,24);
% estimate daily averaged cloud factor
c=cloudfac(daveyd,lat,davesw);
% cloudfac always outputs a column vector, if initial input is
% row vector, convert cf to row vector
if m==1
c = c';
end
% convert daily averaged cloudfactor to boxcar function
if n== 1
c = c*ones(1,24);
c = reshape(c',prod(size(c)),1);
elseif m==1
c = ones(24,1)*c;
c = reshape(c,1,prod(size(c)));
end
c(length(dsw)+1:end)=[];
function cf=cloudfac(yd,lat,davesw)
% CLOUDFAC: computes daily mean cloud factor following Reed (1977).
% cf=CLOUDFAC(yd,lat,davesw) computes the daily average cloud factor
% based on Reed (1977), J. Phys. Ocean., 7, 482-485. Assumes year is
% not a leap year. If yd is negative, then yd=yd+365.
%
% INPUT: yd - yearday (e.g., Jan 10th is 10)
% lat - latitude [deg]
% davesw - average measured insolation [W/m^2]
%
% OUTPUT: cf - daily average cloud factor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convert to column vectors
[N,M]=size(davesw);
if N==1
davesw=davesw';
yd=yd';
end
% check if yd is negative
ind=find(yd<0);
yd(ind)=yd(ind)+365;
% truncate to integer yearday for use with formula
yd=fix(yd);
% determine noon solar altitude
nsa=nsunang(yd,lat);
% compute ratio of observed to clear sky insolation
cssw=clskswr(yd,lat);
QdQ=davesw./cssw;
% correct for measurement error
ind1=find(QdQ>1);
QdQ(ind1)=ones(length(ind1),1);
% compute cloud factor
cf=(1-QdQ+.0019.*nsa)./.62;
% truncate limits
ind2=find(cf>1);
cf(ind2)=ones(length(ind2),1);
ind3=find(cf<0.3);
cf(ind3)=zeros(length(ind3),1);
% reconvert if necessary
if N==1
cf=cf';
end
function nsa=nsunang(yd,lat)
% NSUNANG: computes noon solar altitude angle.
% nsa=NSUNANG(yd,lat) computes the noonday solar altitude angle as a
% function of yearday and latitude using Smithsonian table (see
% Reed (1976), NOAA Tech Memo., ERL PMEL-8, 20 pp., or Reed (1977),
% J. Phys. Ocean., 7, 482-485).
%
% INPUT: yd - yearday (e.g., Jan 10 is 10)
% lat - latitude [deg]
%
% OUTPUT: nsa - noonday solar altitude [deg]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3/8/97: version 1.0
% 8/5/99: version 2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ensure that yd is positive
ind=find(yd<0);
yd(ind)=yd(ind)+365;
yd=fix(yd);
% compute noon solar altitude
k=pi./180;
t=2.*pi.*(yd./365);
d=.397+3.630*sin(t)-22.98.*cos(t)+.040.*sin(2.*t)-.388.*cos(2.*t)...
+.075.*sin(3.*t)-.160.*cos(3.*t);
dl=k.*d;
ll=k.*lat;
z=sin(ll).*sin(dl)+cos(ll).*cos(dl);
nsa=asin(z)./k;