not finished adapting yet!
0001 % not finished adapting yet! 0002 0003 function m = pdd_mb(dem,T,precip, minalt,params) 0004 % PURPOSE: calculate mass balance for a DEM based on a temperature index model 0005 % using positive degree days (PDD). Calculation includes separate 0006 % DDF for ice and snow and a simple firn (refreezing) model 0007 % ------------------------------------------------------------------- 0008 % USAGE: m = meltpdd(dem,T,minalt,params) 0009 % where: [dem] is the input topography grid 0010 % [T] is the mean annual air temperature MAAT at sealevel in degC 0011 % [precip] is the annual precipitation in [m], supplied either as a scalar or matrix 0012 % [minalt] (optional) is the minimum altitude to calculate melt for 0013 % (default min(dem(:)) ) 0014 % [params] (optional) is a vector containing the 4 parameters: 0015 % ddf_ice: degree day factor for ice [mm/d degC] 0016 % ddf_snow: degree day factor for snow [mm/d degC] 0017 % wmax: fraction of melting snow that refreezes and forms 0018 % superimposed ice 0019 % lrate: temperature lapse rate in deg/100m 0020 % a_seas: annual temp fluctuation in degrees. 0021 % tsnow: threshold temperature for precipitation to fall as 0022 % snow 0023 % tau_a: length of the year in days to calculate annual fluctuation 0024 % for, should be 365 0025 % [params] default to [8.0,3.0,0.6,0.6,5.0,1.0,365]; 0026 % Values and model taken from 0027 % Payne and Dongelmans and 0028 % SJ Marshall and GKC Clarke 1999: Ice sheet inception: subgrid 0029 % hypsometric parameterisation of mass balance in an ice sheet 0030 % model. Climate Dynamics 15:533-550 0031 % 0032 %% ------------------------------------------------------------------------- 0033 % OUTPUTS: 0034 % m is an array with ablation values for evey grid cell in dem above 0035 % minalt 0036 % ------------------------------------------------------------------- 0037 % NOTE: Seasonal variation is calculated using cos(2*pi*d/tau_a) 0038 % for d 1:tau_a. 0039 % Melt for areas lower than [minalt] is set to 0. 0040 % 0041 % EXAMPLE: m=meltpdd(randn(100,100)*1000,13,2,500,[8.0,3.0,0.6,0.6,5.0,365]) 0042 % 0043 %% Felix Hebeler, Geography Dept., University Zurich, Jan 2007. 0044 0045 if nargin < 2 0046 error('This function needs at least 2 input parameters (DEM,T)'); 0047 end 0048 0049 % default parameter values 0050 ddf_ice = 8.0; %1 DDF for ice in mm/d degC 0051 ddf_snow = 3.0; %2 DDF for snow in mm/d degC 0052 wmax = 0.6; %3 fraction of melting snow that refreezes and forms superimposed ice 0053 lrate = 0.6; %4 temperatur lapse rate in deg/100m 0054 a_seas = 5.0; %5 annual amplitude of temperature fluctuations in degrees 0055 tsnow = 1.0; %6 threshold temperature for precipitation to fall as snow 0056 tau_a = 365; %7 length of the year in days used in PDD 0057 0058 if ~exist('minalt','var') 0059 minalt=min(dem(:)); 0060 end 0061 if ~exist('params','var') 0062 params=[ddf_ice,ddf_snow,lrate,a_seas,tau_a]; 0063 end 0064 params(1)=params(1)/1000; % convert pdd to m/d degC 0065 params(2)=params(2)/1000; % convert pdd to m/d degC 0066 % initialise 0067 T = T - (dem*params(4)/100); % adjust T to be surface temperature for the DEM 0068 m = zeros(size(dem,1),size(dem,2)); % allocate matrix to hold daily melt 0069 0070 %wfrac = wmax .* precip; % depth of superimposed ice 0071 wfrac = params(3) .* precip; % depth of superimposed ice 0072 pablt = % arg. GLIMMER calculates the number PDDs here... 0073 % calc melt 0074 for d=1:params(6) % seasonal variation in melt: 0075 % tempm = (T - (a_seas * cos(2*pi*d/tau_a)))*ddf_ice; 0076 tempm = (T - (params(5)*cos(2*pi*d/params(7))))*params(1); 0077 tempm(tempm<0)=0; 0078 m = m + tempm; 0079 end 0080 clear tempm; 0081 m(dem<minalt)=0; 0082 clear dem;