Description of parradiation

0001 function [melt,tmelt,rmelt] = parradiation_melt(dem,lat,cs,varargin)
0002 % Function to calculate melt using degree day factor and potential solar
0003 % radiation including diffuse and reflected radiation, following Kumar et
0004 % al 1997(JGIS 11(5), pp.475f
0005 % Usage:    [melt,tmelt,rmelt] = radiation_melt(dem,lat,parameters)
0006 % Input:    Mandatory:
0007 %           DEM (elevation grid), lat (latitude grid), cs (cellsize)
0008 %           Optional:   parameters, consisting of string keyword and value
0009 %           in the form {'slope',slopegrid,'MAAT',12,...}
0010 %                         -slope  = slope in radians
0011 %                         -aspect = aspect in radians
0012 %                         -MAAT  = mean annual air temperature in degrees (default:12)
0013 %                         -tau_a = seasonal temperature variation in degC (default:5)
0014 %                         -lapse = atmospheric lapse rate degC/100m (default:0.65)
0015 %                         -tt    = temperature theshold for melt in degC (default:0)
0016 %                         -S0    = solar constant W m^-2 (default: 1367)
0017 %                         -SRF   = shortwave radiation factor m2 mm W-1 h-1 (default: 0.012)
0018 %                         -a_mean = average albedo (default: 0.4)
0019 %                         -tf    = temperature factor mm h degC-1 (default: 0.05)
0020 %
0021 % Felix Hebeler, University of Zurich, August 2007
0022 
0023 display(['Calculating melt for DEM size= ',num2str(size(dem))])
0024 %% default parameters
0025 MAAT = 12;  % mean annual air temperature at sealevel
0026 n=24;       % timesteps per day (hours in day)
0027 tau_a    = 365;     %length of the year in days
0028 a_seas   = 5.0;     % annual amplitude of temperature fluctuations in degrees
0029 lapse  = 0.65;       % atmospheric lapse rate
0030 tt = 0;             % temperature threshold for melt (default 0 deg)
0031 %M      ;           % air mass ratio parameter 1:30
0032 S0 = 1367;          % solar constant W m^-2   default 1367
0033 SRF=0.012;         % shortwave radiation factor m2 mm W-1 h-1  default 0.0098
0034 albedo = 0.4;       % average albedo for model, if not provided
0035 tf = 0.05;          % temperature factor mm h degC-1 default 0.05
0036 %L=lat;              %latitude
0037 dr= pi/180;         % degree to radians conversion factor
0038 asp=0;
0039 slop=0;
0040 %% parse inputs
0041 if isempty(varargin)~=1     % check if any arguments are given
0042     [m1,n1]=size(varargin);
0043     opts={'aspect';'slope';'MAAT';'tau_a';'a_seas';'lapse';'tt';'S0';'albedo';'tf'};
0044     for i=1:n1;             % check which parameters are given
0045         indi=strcmpi(varargin{i},opts);
0046         ind=find(indi==1);
0047         if isempty(ind)~=1
0048             switch ind
0049                 case 1
0050                     asp=varargin{i+1};
0051                 case 2
0052                     slop=varargin{i+1};
0053                 case 3
0054                     MAAT=varargin{i+1};
0055                 case 4
0056                     tau_a=varargin{i+1};
0057                 case 5
0058                     a_seas=varargin{i+1};
0059                 case 6
0060                     lapse=varargin{i+1};                   
0061                 case 7
0062                     tt=varargin{i+1};                    
0063                 case 8
0064                     S0=varargin{i+1};                    
0065                 case 9
0066                     albedo=varargin{i+1};
0067                 case 10
0068                     tf=varargin{i+1};                    
0069             end
0070         end
0071     end
0072 end
0073 
0074 %%  convert factors
0075 if (sum(asp)==0 || sum(slop)==0)
0076     [dslop,dasp]=get_ders(dem,cs);   % calculate slope and aspect in radians using given cellsize cs
0077     if sum(asp)==0
0078         asp=dasp; clear dasp;
0079     end
0080     if sum(slop)==0
0081         slop=dslop; clear dslop;
0082     end
0083 end
0084 a=albedo;
0085 clear albedo;
0086 %[dummy,L]=meshgrid(1:size(dem,2),lat);   % grid latitude
0087 %clear dummy;
0088 lat=lat*dr;                     % convert to radians
0089 fcirc = 360*dr; % 360 degrees in radians
0090 
0091 %% some setup calculations
0092 rmelt=0;
0093 tmelt=0;
0094 T = MAAT - (dem*lapse/100); % get temperature at DEM surface
0095 sinL=sin(lat);
0096 cosL=cos(lat);
0097 tanL=tan(lat);
0098 sinSlop=sin(slop);
0099 cosSlop=cos(slop);
0100 cosSlop2=cosSlop.*cosSlop;
0101 sinSlop2=sinSlop.*sinSlop;
0102 sinAsp=sin(asp);
0103 cosAsp=cos(asp);
0104 term1 = ( sinL.*cosSlop - cosL.*sinSlop.*cosAsp);
0105 term2 = ( cosL.*cosSlop + sinL.*sinSlop.*cosAsp);
0106 term3 = sinSlop.*sinAsp;
0107 %% loop over year
0108 parfor d = 1:tau_a; 
0109     %display(['Calculating melt for day ',num2str(d)])
0110     % get temperature
0111     Td = T + (a_seas * cos(2*pi*d/tau_a) ); % apply seasonal variation to calc daily temp
0112     Td(Td<=tt)=0; % set temperatures below threshold to 0 so they don't contribute negatively to the melt
0113     %% RADIATION PART
0114     % clear sky solar radiation
0115     I0 = S0 * (1 + 0.0344*cos(fcirc*d/tau_a)); % extraterr rad per day
0116     % sun declination dS
0117     dS = 23.45 * dr* sin(fcirc * ( (284+d)/tau_a ) ); %in radians, correct/verified
0118     % angle at sunrise/sunset
0119     hsr = real(acos(-tanL*tan(dS)));  % angle at sunrise
0120     It=round(12*(1+mean(hsr(:))/pi)-12*(1-mean(hsr(:))/pi)); % calc daylength
0121 %  daily loop
0122     I=0;
0123     for t=1:It % loop over sunshine hours
0124         % if accounting for shading should be included, calc hillshade here
0125         % and solar azimuth is needed...
0126         % hourangle of sun hs
0127         hs=hsr-(pi*t/It);               % hs(t)
0128         %solar angle and azimuth
0129         %alpha = asin(sinL*sin(dS)+cosL*cos(dS)*cos(hs));
0130         sinAlpha = sinL.*sin(dS)+cosL.*cos(dS).*cos(hs);
0131         %alpha_s = asin(cos(dS)*sin(hs)/cos(alpha));
0132         % <-- could fit in the hillshade part here to honour relief shading
0133         % correction  using atmospheric transmissivity taub_b
0134         M=sqrt(1229+(614.*sinAlpha).^2)-614.*sinAlpha; % Air mass ratio
0135         tau_b = 0.56 * (exp(-0.65*M) + exp(-0.095*M));
0136         tau_d = 0.271-0.294*tau_b; % radiation diffusion coefficient for diffuse insolation
0137         tau_r = 0.271+0.706*tau_b; % reflectance transmitivity
0138         % correct for local incident angle
0139         cos_i = (sin(dS).*term1) + (cos(dS).*cos(hs).*term2) + (cos(dS).*term3.*sin(hs));
0140         Is = I0 * tau_b; % potential incoming shortwave radiation at surface normal (equator)
0141         % R = potential clear sky solar radiation W m2
0142         R = Is .* cos_i;
0143         R(R<0)=0;  % kick out negative values
0144         Id = I0 .* tau_d .* cosSlop2./ 2.*sinAlpha; %diffuse radiation;
0145         %Id=0;
0146         Ir = I0 .* a .* tau_r .* sinSlop2./ 2.* sinAlpha; % reflectance
0147         R= R + Id + Ir;
0148         R(R<0)=0; 
0149         I=I+R;% melt from solar radiation per day (sunshine hours)
0150      end % end of sun hours in day loop
0151     
0152 %  add up radiation part melt for every day
0153     rtemp=(SRF.*(1-a).*I);  
0154     ttemp=tf.*Td .*n;        % no negative melt can occur because Td>=tt
0155     rtemp(ttemp==0)=0;       % if ttemp=0, Temp is lower tt and no melt occurs
0156     rmelt = rmelt + rtemp;   % sum of melt from radiation part
0157     tmelt = tmelt + ttemp;
0158   
0159 end   % end of days in year loop
0160 melt= tmelt+ rmelt; 
0161 
0162 %%
0163 function [grad,asp] = get_ders(dem,cs)
0164 % calculate slope and aspect (rad) using GRADIENT function
0165 [fx,fy] = gradient(dem,cs,cs); % uses simple, unweighted gradient of immediate neighbours
0166 [asp,grad]=cart2pol(fy,fx); % convert to carthesian coordinates
0167 grad=atan(grad); %steepest slope
0168 asp=asp.*-1+pi; % convert asp 0 facing south
parradiation_melt

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
