clskswr

 CLSKSWR: computes clear sky insolation following Seckel&Beaudry (1973).
 sradav = CLSKSWR(yd,lat) computes average clear sky solar insolation
 based on the Seckel and Beaudry (1973) formula presented in Reed (1977),
 J. Phys. Ocean., 7, 482-485. Assumes the year is not a leap year.

 INPUT:   yd  -  yearday (e.g., Jan 10th is 10) 
          lat -  latitude [deg]
  
 OUTPUT:  sradav - clear sky mean daily insolation [W/m^2]

0001 function sradav = clskswr(yd,lat)
0002 % CLSKSWR: computes clear sky insolation following Seckel&Beaudry (1973).
0003 % sradav = CLSKSWR(yd,lat) computes average clear sky solar insolation
0004 % based on the Seckel and Beaudry (1973) formula presented in Reed (1977),
0005 % J. Phys. Ocean., 7, 482-485. Assumes the year is not a leap year.
0006 %
0007 % INPUT:   yd  -  yearday (e.g., Jan 10th is 10)
0008 %          lat -  latitude [deg]
0009 %
0010 % OUTPUT:  sradav - clear sky mean daily insolation [W/m^2]
0011 
0012 % NOTE: The output appears to be very similar to what you would get by
0013 %       averaging soradna.m output over a day and then assuming an
0014 %       atmospheric transmission of 0.7 (with differences of order 10%
0015 %       for latitudes below 40N, and increasing to 30% in winter at 60N).
0016 %       In absolute terms the agreement is to within +/-25 W/m^2.
0017 
0018 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0019 % 3/8/97: version 1.0
0020 % 8/28/98: version 1.1 (vectorized by RP)
0021 % 8/5/99: version 2.0
0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 
0024 radlat = pi*lat/180;
0025 
0026 if length(lat)==1,
0027  lat=lat+zeros(size(yd));
0028 end;
0029 
0030 % check if yd is negative
0031 ind=find(yd<0);
0032 yd(ind)=yd(ind)+365;
0033 % truncate to integer yearday for use with formula
0034 yd=fix(yd);
0035 
0036 phi = (yd-21)*360/365;
0037 phi = pi*phi/180;
0038 
0039 sradav=zeros(size(yd))+zeros(size(lat))+NaN;
0040 
0041 ii= lat>=-20 & lat<40;
0042 if any(ii(:)),
0043   a0 = -15.82 + 326.87*cos(radlat(ii));
0044   a1 = 9.63 + 192.44*cos(radlat(ii)+pi/2);
0045   b1 = -3.27 + 108.70*sin(radlat(ii));
0046   a2 = -0.64 + 7.80*sin(2*(radlat(ii)-pi/4));
0047   b2 = -0.50 + 14.42*cos(2*(radlat(ii)-5*pi/180));
0048   sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii));
0049 end
0050 
0051 ii= lat>=40 & lat<=60;;
0052 if any(ii(:)),
0053   l2=lat(ii).^2;
0054   a0 = 342.61 - 1.97*lat(ii) - 0.018*l2;
0055   a1 = 52.08 - 5.86*lat(ii) + 0.043*l2;
0056   b1 = -4.80 + 2.46*lat(ii) -0.017*l2;
0057   a2 = 1.08 - 0.47*lat(ii) + 0.011*l2;
0058   b2 = -38.79 + 2.43*lat(ii) - 0.034*l2;
0059   sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii));
0060 end
0061 
0062 if any(lat>60 | lat<-20)
0063   warning('Formula only works for latitudes 20S-60N, see help text for further help')
0064 end
0065 
0066