function sradav = clskswr(yd,lat)

% 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]



% NOTE: The output appears to be very similar to what you would get by

%       averaging soradna.m output over a day and then assuming an 

%       atmospheric transmission of 0.7 (with differences of order 10% 

%       for latitudes below 40N, and increasing to 30% in winter at 60N).

%       In absolute terms the agreement is to within +/-25 W/m^2.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3/8/97: version 1.0

% 8/28/98: version 1.1 (vectorized by RP)

% 8/5/99: version 2.0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



radlat = pi*lat/180;



if length(lat)==1,

 lat=lat+zeros(size(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);



phi = (yd-21)*360/365;

phi = pi*phi/180;



sradav=zeros(size(yd))+zeros(size(lat))+NaN;



ii= lat>=-20 & lat<40;

if any(ii(:)),

  a0 = -15.82 + 326.87*cos(radlat(ii));

  a1 = 9.63 + 192.44*cos(radlat(ii)+pi/2);

  b1 = -3.27 + 108.70*sin(radlat(ii));

  a2 = -0.64 + 7.80*sin(2*(radlat(ii)-pi/4));

  b2 = -0.50 + 14.42*cos(2*(radlat(ii)-5*pi/180));

  sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii));

end



ii= lat>=40 & lat<=60;;

if any(ii(:)),

  l2=lat(ii).^2;

  a0 = 342.61 - 1.97*lat(ii) - 0.018*l2;

  a1 = 52.08 - 5.86*lat(ii) + 0.043*l2;

  b1 = -4.80 + 2.46*lat(ii) -0.017*l2;

  a2 = 1.08 - 0.47*lat(ii) + 0.011*l2;

  b2 = -38.79 + 2.43*lat(ii) - 0.034*l2;

  sradav(ii) = a0 + a1.*cos(phi(ii)) + b1.*sin(phi(ii)) + a2.*cos(2*phi(ii)) + b2.*sin(2*phi(ii));

end



if any(lat>60 | lat<-20)

  warning('Formula only works for latitudes 20S-60N, see help text for further help')

end





