# Unverified solution in MATLAB

Jump to: navigation, search
```% Possible solution to Kees' problem in MATLAB
%
% Based on eqns. presented on
% http://websrv.cs.umt.edu/isis/index.php/Kees'_assignment
% Possibly useful as a solution template for those familiar with MATLAB.
% Note possible problem in lines 55-58.

clear all
close all
clc

% Define constants.
g = 9.8;            % m/ s^ 2; gravitational acceleration
rho = 920;          % kg/ m^ 3; density of ice
delta_x = 100;      % m; distance step
x_max = 25* 10^ 3;  % m; length of problem domain
delta_t = 10^ -2;   % yr; time step (needs to be small, about 10^ -2 yr)
t_max = 10^ 1;      % yr; integration time of model
M_0 = 4;            % m; net accumulation at upper end of problem domain
S_M = 4/ 10^ 4;     % 1/ m; net accumulation gradient (note change from
% value given on wiki page, above)
bed_slope = 0.1;    % slope of glacier bed
Glen_n = 3;         % d'less; exponent in Glen's flow law
Glen_A = 10^ -16;   % Pa^ -3/ yr; coefficient in Glen's flow law

% Calculate derived constants.
C = (2/ (Glen_n+ 2))* Glen_A* (rho* g)^ Glen_n;

% Set up initial condition and plotting variables.  Initialize vectors.
x = 0: delta_x: x_max;      % m; distances from upper end of domain for
% each node
t = 0: delta_t: t_max;      % yr; times corresponding to each time step
z_bed = wrev(bed_slope* x); % m; bed height as a fn of x
num_x = numel(x);           % number of nodes
num_t = numel(t);           % number of time steps
D = zeros(1, num_x- 1);     % m^ 2/ yr; diffusivities at staggered points
H = zeros(1, num_x);        % m; thickness of ice at different points
M = zeros(1, num_x);        % m/ time step; net mass balance per time step
num_D = numel(D);           % number of points in offset grid
phi = zeros(1, num_D);      % m^ 2/ yr?; flux of matter between nodes

% Calculate the net mass balance per time step within the problem domain.
M = (M_0- S_M* x)* delta_t; % m/ time step

% Step through time.
for count1 = 1: 1: num_t;
for count2 = 1: 1: num_D;
D(count2) = C* ((H(count2)+ H(count2+ 1))/ 2)^ (Glen_n+ 2)* ...
((H(count2+ 1)+ z_bed(count2+ 1))- ...
(H(count2)+ z_bed(count2)))^ (Glen_n- 1);
phi(count2) = D(count2)* (((H(count2+ 1)+ z_bed(count2+ 1))- ...
(H(count2)+ z_bed(count2)))/ delta_x);
end
for count3 = 2: 1: (num_x- 1);
H(count3) = H(count3)+ (phi(count3)- phi(count3- 1))/ ...
delta_x+ M(count3);
%         H(count3) = H(count3)- (phi(count3)- phi(count3- 1))/ ...
%             delta_x+ M(count3);
end
H(1) = 0;
for count4 = 1: 1: num_x;
if H(count4) < 0;
H(count4) = 0;
end
end
end

% Plot results.
ice_top = z_bed+ H; % m; elevation of ice surface above lowest node
plot(x, ice_top, 'b-')
hold on
plot(x, z_bed, 'k-')
xlabel('Horizontal distance (m)')
ylabel('Elevation (m)')```