solarradiation_shaded

 Usage: srad = solarradiation(dem,lat,cs)
 Function to calculate melt using degree day factor and potential solar
 radiation including diffuse and reflected radiation

 Input: DEM, latitude as vector, cellsize in meters

0001 function srad = solarradiation_shaded(dem,lat,cs)
0002 % Usage: srad = solarradiation(dem,lat,cs)
0003 % Function to calculate melt using degree day factor and potential solar
0004 % radiation including diffuse and reflected radiation
0005 %
0006 % Input: DEM, latitude as vector, cellsize in meters
0007 
0008 %display(['DEM size= ',num2str(size(dem)),' lat size= ',num2str(size(lat)),' CS= ',num2str(cs)])
0009 %% parameters
0010 %It ;               % total hours of daily sunshine (calculated inline)
0011 n = 24;             % time steps per day (hours in day)
0012 tau_a    = 365;     %length of the year in days
0013 %M ;           % air mass ratio parameter (calculated inline)
0014 S0 = 1367;          % solar constant W m^-2   default 1367
0015 r = 0.20;           % ground reflectance coefficient (more sensible to give as input)
0016 %L=lat;              %latitude
0017 dr= 0.0174532925;   % degree to radians conversion factor
0018 
0019 %%  convert factors
0020 [slop,asp]=get_ders(dem,cs);   % calculate slope and aspect in radians using given cellsize cs
0021 [dummy,L]=meshgrid(1:size(dem,2),lat);   % grid latitude
0022 clear dummy;
0023 L=L*dr;                     % convert to radians
0024 fcirc = 360*dr; % 360 degrees in radians
0025 
0026 %% some setup calculations
0027 srad=0;
0028 sinL=sin(L);
0029 cosL=cos(L);
0030 tanL=tan(L);
0031 sinSlop=sin(slop);
0032 cosSlop=cos(slop);
0033 cosSlop2=cosSlop.*cosSlop;
0034 sinSlop2=sinSlop.*sinSlop;
0035 sinAsp=sin(asp);
0036 cosAsp=cos(asp);
0037 term1 = ( sinL.*cosSlop - cosL.*sinSlop.*cosAsp);
0038 term2 = ( cosL.*cosSlop + sinL.*sinSlop.*cosAsp);
0039 term3 = sinSlop.*sinAsp;
0040 %% loop over year
0041 for d = 1:tau_a; 
0042     %display(['Calculating melt for day ',num2str(d)])
0043     % clear sky solar radiation
0044     I0 = S0 * (1 + 0.0344*cos(fcirc*d/tau_a)); % extraterr rad per day
0045     % sun declination dS
0046     dS = 23.45 * dr* sin(fcirc * ( (284+d)/tau_a ) ); %in radians, correct/verified
0047     % angle at sunrise/sunset
0048     % t = 1:It; % sun hour
0049     hsr = real(acos(-tanL*tan(dS)));  % angle at sunrise
0050     % this only works for latitudes up to 66.5 deg N! Workaround:
0051     % hsr(hsr<-1)=acos(-1);
0052     % hsr(hsr>1)=acos(1);
0053     It=round(12*(1+mean(hsr(:))/pi)-12*(1-mean(hsr(:))/pi)); % calc daylength in hours
0054 %%  daily loop
0055     I=0; % reset daily solar radiation
0056     for t=1:It % loop over sunshine hours
0057             
0058         % hourangle of sun hs
0059         hs=hsr-(pi*t/It);               % hs(t)
0060         %solar angle and azimuth
0061         alpha = asin(sinL.*sin(dS)+cosL.*cos(dS).*cos(hs)); % solar altitude angle
0062         sinAlpha = sin(alpha); % sinus of alpha
0063        % alpha_s = asin(cos(dS).*sin(hs)./cos(alpha));  % solar azimuth angle
0064         % calculate hillshade for relief shading
0065        % h = ( (cos(alpha).*cosSlop ) + ( sinAlpha.*sinSlop.*cos(alpha_s-asp)) ); %
0066        % h(h>0)=1; % set h to one if not shaded
0067        % h(h<0)=0; % set h to 0 is shaded.
0068         % correction  using atmospheric transmissivity taub_b
0069         M=sqrt(1229+((614.*sinAlpha)).^2)-614.*sinAlpha; % Air mass ratio
0070         tau_b = 0.56 * (exp(-0.65*M) + exp(-0.095*M));
0071         tau_d = 0.271-0.294*tau_b; % radiation diffusion coefficient for diffuse insolation
0072         tau_r = 0.271+0.706*tau_b; % reflectance transmitivity
0073         % correct for local incident angle
0074         cos_i = (sin(dS).*term1) + (cos(dS).*cos(hs).*term2) + (cos(dS).*term3.*sin(hs));
0075         Is = I0 * tau_b; % potential incoming shortwave radiation at surface normal (equator)
0076         % R = potential clear sky solar radiation W m2
0077         R = Is .* cos_i;% * h;
0078         R(R<0)=0;  % kick out negative values
0079         Id = I0 .* tau_d .* cosSlop2./ 2.*sinAlpha; %diffuse radiation;
0080         Ir = I0 .* r .* tau_r .* sinSlop2./ 2.* sinAlpha; % reflectance
0081         R= R + Id + Ir;
0082         R(R<0)=0; 
0083         I=I+R;% solar radiation per day (sunshine hours)
0084      end % end of sun hours in day loop
0085 %%  add up radiation part melt for every day
0086     srad = srad + I;
0087 end   % end of days in year loop
0088 
0089 
0090 %%
0091 function [grad,asp] = get_ders(dem,cs)
0092 % calculate slope and aspect (deg) using GRADIENT function
0093 [fx,fy] = gradient(dem,cs,cs); % uses simple, unweighted gradient of immediate neighbours
0094 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0095 grad=atan(grad); %steepest slope
0096 asp=asp.*-1+pi; % convert asp 0 facing south