PURPOSE: calculate melt for a DEM based on a simple temperature index model
using positive degree days (PDD).
-------------------------------------------------------------------
USAGE: m = meltpdd(dem,T,minalt,params)
where: [dem] is the input topography grid
[T] is the mean annual air temperature MAAT at sealevel in degC
[minalt] (optional) is the minimum altitude to calculate melt for
(default min(dem(:)) )
[params] (optional) is a vector containing the 4 parameters:
ddf_ice: degree day factor for ice [mm/d degC]
lrate: temperature lapse rate in deg/100m
a_seas: annual temp fluctuation in degrees.
tau_a: length of the year in days to calculate annual fluctuation
for, should be 365
[params] default to [8.0,0.6,5.0,365];
Values and model taken from
SJ Marshall and GKC Clarke 1999: Ice sheet inception: subgrid
hypsometric parameterisation of mass balance in an ice sheet
model. Climate Dynamics 15:533-550
% -------------------------------------------------------------------------
OUTPUTS:
m is an array with melt values for evey grid cell in dem above
minalt
-------------------------------------------------------------------
NOTE: Seasonal variation is calculated using cos(2*pi*d/tau_a)
for d 1:tau_a.
Melt for areas lower than [minalt] is set to 0.
EXAMPLE: m=meltpdd(randn(100,100)*1000,13,500,[8.0,0.6,5.0,365])
% Felix Hebeler, Geography Dept., University Zurich, Jan 2007.0001 function m = parmeltpdd(dem,T,minalt,params) 0002 % PURPOSE: calculate melt for a DEM based on a simple temperature index model 0003 % using positive degree days (PDD). 0004 % ------------------------------------------------------------------- 0005 % USAGE: m = meltpdd(dem,T,minalt,params) 0006 % where: [dem] is the input topography grid 0007 % [T] is the mean annual air temperature MAAT at sealevel in degC 0008 % [minalt] (optional) is the minimum altitude to calculate melt for 0009 % (default min(dem(:)) ) 0010 % [params] (optional) is a vector containing the 4 parameters: 0011 % ddf_ice: degree day factor for ice [mm/d degC] 0012 % lrate: temperature lapse rate in deg/100m 0013 % a_seas: annual temp fluctuation in degrees. 0014 % tau_a: length of the year in days to calculate annual fluctuation 0015 % for, should be 365 0016 % [params] default to [8.0,0.6,5.0,365]; 0017 % Values and model taken from 0018 % SJ Marshall and GKC Clarke 1999: Ice sheet inception: subgrid 0019 % hypsometric parameterisation of mass balance in an ice sheet 0020 % model. Climate Dynamics 15:533-550 0021 % 0022 %% ------------------------------------------------------------------------- 0023 % OUTPUTS: 0024 % m is an array with melt values for evey grid cell in dem above 0025 % minalt 0026 % ------------------------------------------------------------------- 0027 % NOTE: Seasonal variation is calculated using cos(2*pi*d/tau_a) 0028 % for d 1:tau_a. 0029 % Melt for areas lower than [minalt] is set to 0. 0030 % 0031 % EXAMPLE: m=meltpdd(randn(100,100)*1000,13,500,[8.0,0.6,5.0,365]) 0032 % 0033 %% Felix Hebeler, Geography Dept., University Zurich, Jan 2007. 0034 0035 if nargin < 2 0036 error('This function needs at least 2 input parameters (DEM,T)'); 0037 end 0038 0039 % default parameter values 0040 ddf_ice = 8.0; % DDF for ice in mm/d degC 0041 lrate = 0.6; % temperatur lapse rate in deg/100m 0042 a_seas = 5.0; % annual amplitude of temperature fluctuations in degrees 0043 tau_a = 365; %length of the year in days used in PDD 0044 0045 if ~exist('minalt','var') 0046 minalt=min(dem(:)); 0047 end 0048 0049 if ~exist('params','var') 0050 params=[ddf_ice,lrate,a_seas,tau_a]; 0051 end 0052 params(1)=params(1)/1000; % convert pdd to m/d degC 0053 % initialise 0054 T = T - (dem*params(2)/100); % adjust T to be surface temperature for the DEM 0055 %tempm=zeros(size(dem,1),size(dem,2)); % allocate matrix to hold daily melt 0056 m = zeros(size(dem,1),size(dem,2)); % matrix for meltsum 0057 % calc melt 0058 parfor d=1:params(4) % seasonal variation in melt: 0059 % tempm = (T - (a_seas * cos(2*pi*d/tau_a)))*ddf_ice; 0060 tempm = (T - (params(3)*cos(2*pi*d/params(4))))*params(1); 0061 tempm(tempm<0)=0; 0062 m = m + tempm; 0063 end 0064 clear tempm; 0065 m(dem<minalt)=0; 0066 clear dem;