pellimeltd

 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 [melt,tmelt,rmelt] = pellimeltd(dem,lat,cs,MAAT,M)
0002 % Function to calculate melt using degree day factor and potential solar
0003 % radiation including diffuse and reflected radiation
0004 %
0005 % Input: DEM, latitude as vector, cellsize in meters
0006 display(['DEM size= ',num2str(size(dem)),' lat size= ',num2str(size(lat)),' CS= ',num2str(cs)])
0007 %% parameters
0008 %MAAT = 12;          % mean annual air temperature at sealevel
0009 %It = 12 ;           % total hours of daily sunshine
0010 n = 24;             % time steps per day (hours in day, leave this, and only change together with It!)
0011 tau_a    = 365;     %length of the year in days
0012 a_seas   = 5.0;     % annual amplitude of temperature fluctuations in degrees
0013 lapse  = 0.65;       % atmospheric lapse rate
0014 tt = 0;             % temperature threshold for melt (default 0 deg)
0015 %M = 3.5;           % air mass ratio parameter 1:30
0016 S0 = 1367;          % solar constant W m^-2   default 1367
0017 r = 0.20;           % ground reflectance coefficient
0018 SRF=0.012;         % shortwave radiation factor m2 mm W-1 h-1  default 0.0098
0019 a_ice = 0.4;        % albedo ice
0020 a_snow = 0.7;       % albedo snow
0021 a_mean = 0.4;       % average albedo for model
0022 tf = 0.05;          % temperature factor mm h degC-1 default 0.05
0023 %L=lat;              %latitude
0024 dr= 0.0174532925;   % degree to radians conversion factor
0025 
0026 %%  convert factors
0027 [slop,asp]=get_ders(dem,cs);   % calculate slope and aspect in radians using given cellsize cs
0028 a=a_mean;
0029 [dummy,L]=meshgrid(1:size(dem,2),lat);   % grid latitude
0030 clear dummy;
0031 L=L*dr;                     % convert to radians
0032 fcirc = 360*dr; % 360 degrees in radians
0033 
0034 %% some setup calculations
0035 rmelt=0;
0036 tmelt=0;
0037 T = MAAT - (dem*lapse/100); % get temperature at DEM surface
0038 sinL=sin(L);
0039 cosL=cos(L);
0040 tanL=tan(L);
0041 sinSlop=sin(slop);
0042 cosSlop=cos(slop);
0043 cosSlop2=cosSlop.*cosSlop;
0044 sinSlop2=sinSlop.*sinSlop;
0045 sinAsp=sin(asp);
0046 cosAsp=cos(asp);
0047 term1 = ( sinL.*cosSlop - cosL.*sinSlop.*cosAsp);
0048 term2 = ( cosL.*cosSlop + sinL.*sinSlop.*cosAsp);
0049 term3 = sinSlop.*sinAsp;
0050 %% loop over year
0051 for d = 1:tau_a; 
0052     %display(['Calculating melt for day ',num2str(d)])
0053     % get temperature
0054     Td = T + (a_seas * cos(2*pi*d/tau_a) ); % apply seasonal variation to calc daily temp
0055     Td(Td<=tt)=0; % set temperatures below threshold to 0 so they don't contribute negatively to the melt
0056     %% RADIATION PART
0057     % clear sky solar radiation
0058     I0 = S0 * (1 + 0.0344*cos(fcirc*d/tau_a)); % extraterr rad per day
0059     % correction  using atmospheric transmissivity taub_b
0060     % calc M here if intended
0061     tau_b = 0.56 * (exp(-0.65*M) + exp(-0.095*M));
0062     tau_d = 0.271-0.294*tau_b; % radiation diffusion coefficient for diffuse insolation
0063     tau_r = 0.271+0.706*tau_b; % reflectance transmitivity
0064         
0065     Is = I0 * tau_b; % potential incoming shortwave radiation at surface normal (equator)
0066     % sun declination dS
0067     dS = 23.45 * dr* sin(fcirc * ( (284+d)/tau_a ) ); %in radians, correct/verified
0068     % angle at sunrise/sunset
0069     % t = 1:It; % sun hour
0070     hsr = real(acos(-tanL*tan(dS)));  % angle at sunrise
0071         % this only works for latitudes up to 66.5 deg N! Workaround:
0072     % hsr(hsr<-1)=acos(-1);
0073     % hsr(hsr>1)=acos(1);
0074     It=round(12*(1+mean(hsr(:))/pi)-12*(1-mean(hsr(:))/pi)); % calc daylength
0075 %%  daily loop
0076     I=0;
0077     for t=1:It % loop over sunshine hours
0078         % if accounting for shading should be included, calc hillshade here
0079         % hourangle of sun hs
0080         hs=hsr-(pi*t/It);               % hs(t)
0081         %solar angle and azimuth
0082         %alpha = asin(sinL*sin(dS)+cosL*cos(dS)*cos(hs));
0083         sinAlpha = sinL.*sin(dS)+cosL.*cos(dS).*cos(hs);
0084         %alpha_s = asin(cos(dS)*sin(hs)/cos(alpha));
0085         % <-- could fit in the hillshade part here to honour relief shading
0086         % correct for local incident angle
0087         cos_i = (sin(dS).*term1) + (cos(dS).*cos(hs).*term2) + (cos(dS).*term3.*sin(hs));
0088         % R = potential clear sky solar radiation W m2
0089         R = Is * cos_i;
0090         R(R<0)=0;  % kick out negative values
0091         Id = I0 .* tau_d .* cosSlop2./ 2.*sinAlpha; %diffuse radiation;
0092         Ir = I0 .* r .* tau_r .* sinSlop2./ 2.* sinAlpha; % reflectance
0093         R= R + Id + Ir;
0094         R(R<0)=0; 
0095         I=I+R;% melt from solar radiation per day (sunshine hours)
0096      end % end of sun hours in day loop
0097     
0098 %%  add up radiation part melt for every day
0099     rtemp=(SRF*(1-a).*I);  
0100     ttemp=tf.*Td .*n;        % no negative melt can occur because Td>=tt
0101     rtemp(ttemp==0)=0;       % if ttemp=0, Temp is lower tt and no melt occurs
0102     rmelt = rmelt + rtemp;   % sum of melt from radiation part
0103     tmelt = tmelt + ttemp;
0104   
0105 end   % end of days in year loop
0106 melt= tmelt+ rmelt; 
0107 
0108 %%
0109 function [grad,asp] = get_ders(dem,cs)
0110 % calculate slope and aspect (deg) using GRADIENT function
0111 [fx,fy] = gradient(dem,cs,cs); % uses simple, unweighted gradient of immediate neighbours
0112 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0113 grad=atan(grad); %steepest slope
0114 asp=asp.*-1+pi; % convert asp 0 facing south