pellimeltdm

 Usage: [melt,tmelt,rmelt] = pellimeltdm(dem,lat,cs,MAAT)
 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, MAAT

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