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