23
23
% Called by: VERIFYSIA, ROUGHICE
24
24
25
25
% physical constants
26
- g = 9.81 ; rho = 910.0 ; secpera = 31556926 ;
27
- A = 1 .0e- 16 / secpera ; Gamma = 2 * A * (rho * g )^3 / 5 ; % see Bueler et al (2005)
28
- H = H0 ;
26
+ g = 9.81 ; rho = 910.0 ; secpera = 31556926 ;
27
+ A = 1 .0e- 16 / secpera ; Gamma = 2 * A * (rho * g )^3 / 5 ; % see Bueler et al (2005)
29
28
30
- dx = 2 * Lx / J ; dy = 2 * Ly / K ;
31
- N = ceil(tf / deltat ); deltat = tf / N ;
32
- j = 2 : J ; k = 2 : K ; % interior indices
33
- nk = 3 : K + 1 ; sk = 1 : K - 1 ; ej = 3 : J + 1 ; wj = 1 : J - 1 ; % north,south,east,west
29
+ N = ceil(tf / deltat ); deltat = tf / N ;
30
+ dx = 2 * Lx / J ; dy = 2 * Ly / K ; j = 2 : J ; k = 2 : K ; % interior indices
31
+ nk = 3 : K + 1 ; sk = 1 : K - 1 ; ej = 3 : J + 1 ; wj = 1 : J - 1 ; % north,south,east,west
34
32
35
33
fprintf(' solving SIA for 0.0 < t < %.3f a\n ' ,tf / secpera )
36
34
t = 0 ; dtlist = [];
35
+ H = H0 ;
37
36
for n= 1 : N
38
37
% staggered grid thicknesses
39
38
Hup = 0.5 * ( H(j ,nk ) + H(j ,k ) ); % up and down
64
63
fprintf(' \n SIA solver done\n ' )
65
64
% to plot the diffusivity for the final state:
66
65
% x = linspace(-Lx+dx,Lx-dx,J-1); y = linspace(-Ly+dy,Ly-dy,K-1);
67
- % figure(99), surf(x,y,0.25*(Dup+Ddn+Drt+Dlt))
66
+ % figure(99), surf(x,y,0.25*(Dup+Ddn+Drt+Dlt))
0 commit comments