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 cover0001 function c = reedcf(yd,lat,dsw); 0002 % REEDCF: computes daily mean cloud cover following Reed (1977). 0003 % c = REEDCF(yd,lat,dsw) computes daily averaged cloud cover c from 0004 % yearday, latitude, and measured insolation following Reed (1977), J. Phys. 0005 % Oceanog., 7, 482-485. Assumes hourly input series are either both 0006 % column or both row vectors of equal length. c is output as a boxcar 0007 % function with c constant over a 24 hr period. 0008 % 0009 % INPUT: yd - yearday (e.g., Jan 10 is yd=10) 0010 % lat - latitude [deg] 0011 % dsw - (measured) insolation [W/m^2] 0012 % 0013 % OUTPUT: c - daily averaged cloud cover 0014 0015 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0016 % 3/8/97: version 1.0 0017 % 8/5/99: version 2.0 0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0019 0020 % determine if input is either row or column vectors 0021 [m,n] = size(dsw); 0022 0023 % first daily average swr and yd 0024 davesw = binave(dsw,24); 0025 daveyd = binave(yd,24); 0026 0027 % estimate daily averaged cloud factor 0028 c=cloudfac(daveyd,lat,davesw); 0029 0030 % cloudfac always outputs a column vector, if initial input is 0031 % row vector, convert cf to row vector 0032 if m==1 0033 c = c'; 0034 end 0035 0036 % convert daily averaged cloudfactor to boxcar function 0037 if n== 1 0038 c = c*ones(1,24); 0039 c = reshape(c',prod(size(c)),1); 0040 0041 elseif m==1 0042 c = ones(24,1)*c; 0043 c = reshape(c,1,prod(size(c))); 0044 end 0045 c(length(dsw)+1:end)=[]; 0046 0047 0048 function cf=cloudfac(yd,lat,davesw) 0049 % CLOUDFAC: computes daily mean cloud factor following Reed (1977). 0050 % cf=CLOUDFAC(yd,lat,davesw) computes the daily average cloud factor 0051 % based on Reed (1977), J. Phys. Ocean., 7, 482-485. Assumes year is 0052 % not a leap year. If yd is negative, then yd=yd+365. 0053 % 0054 % INPUT: yd - yearday (e.g., Jan 10th is 10) 0055 % lat - latitude [deg] 0056 % davesw - average measured insolation [W/m^2] 0057 % 0058 % OUTPUT: cf - daily average cloud factor 0059 0060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0061 % 3/8/97: version 1.0 0062 % 8/5/99: version 2.0 0063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0064 0065 % convert to column vectors 0066 [N,M]=size(davesw); 0067 if N==1 0068 davesw=davesw'; 0069 yd=yd'; 0070 end 0071 0072 % check if yd is negative 0073 ind=find(yd<0); 0074 yd(ind)=yd(ind)+365; 0075 % truncate to integer yearday for use with formula 0076 yd=fix(yd); 0077 0078 % determine noon solar altitude 0079 nsa=nsunang(yd,lat); 0080 0081 % compute ratio of observed to clear sky insolation 0082 cssw=clskswr(yd,lat); 0083 QdQ=davesw./cssw; 0084 0085 % correct for measurement error 0086 0087 ind1=find(QdQ>1); 0088 QdQ(ind1)=ones(length(ind1),1); 0089 0090 % compute cloud factor 0091 cf=(1-QdQ+.0019.*nsa)./.62; 0092 0093 % truncate limits 0094 0095 ind2=find(cf>1); 0096 cf(ind2)=ones(length(ind2),1); 0097 ind3=find(cf<0.3); 0098 cf(ind3)=zeros(length(ind3),1); 0099 0100 % reconvert if necessary 0101 if N==1 0102 cf=cf'; 0103 end 0104 0105 0106 function nsa=nsunang(yd,lat) 0107 % NSUNANG: computes noon solar altitude angle. 0108 % nsa=NSUNANG(yd,lat) computes the noonday solar altitude angle as a 0109 % function of yearday and latitude using Smithsonian table (see 0110 % Reed (1976), NOAA Tech Memo., ERL PMEL-8, 20 pp., or Reed (1977), 0111 % J. Phys. Ocean., 7, 482-485). 0112 % 0113 % INPUT: yd - yearday (e.g., Jan 10 is 10) 0114 % lat - latitude [deg] 0115 % 0116 % OUTPUT: nsa - noonday solar altitude [deg] 0117 0118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0119 % 3/8/97: version 1.0 0120 % 8/5/99: version 2.0 0121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0122 0123 % ensure that yd is positive 0124 ind=find(yd<0); 0125 yd(ind)=yd(ind)+365; 0126 yd=fix(yd); 0127 0128 % compute noon solar altitude 0129 k=pi./180; 0130 t=2.*pi.*(yd./365); 0131 d=.397+3.630*sin(t)-22.98.*cos(t)+.040.*sin(2.*t)-.388.*cos(2.*t)... 0132 +.075.*sin(3.*t)-.160.*cos(3.*t); 0133 dl=k.*d; 0134 ll=k.*lat; 0135 z=sin(ll).*sin(dl)+cos(ll).*cos(dl); 0136 nsa=asin(z)./k;