pellimelt

 Usage: [melt,tmelt,rmelt] = pellimelt(dem,lat,cs,MAAT,M)

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