Home > air_sea > reedcf.m

reedcf

PURPOSE ^

REEDCF: computes daily mean cloud cover following Reed (1977).

SYNOPSIS ^

function c = reedcf(yd,lat,dsw);

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 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;

Generated on Thu 26-Apr-2007 12:23:13 by m2html © 2003