-
Notifications
You must be signed in to change notification settings - Fork 0
/
ism_slidinglaw_linearinit.m
80 lines (57 loc) · 1.68 KB
/
ism_slidinglaw_linearinit.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
function [alpha] = ism_slidinglaw_linearinit(uv,B2,N,vv,aa,pp,gg,oo)
%% Field Equations for SSA velocities
% Inputs:
% v solution struct
% aa prescribed fields, including inputs and boundary conditions
% pp parameters
% gg grid and operators
% oo options
% Outputs:
% aa
%% Linear
if isequal(oo.slidinglaw, 1)
alpha = B2;
else
%% Non Linear
%Setup variables
%Setup variables
u = uv(1:gg.nua);
v = uv(gg.nua+1:end);
U = ((gg.c_uh*u).^2 + (gg.c_vh*v).^2).^(1/2);
%% Basal velocities in the case of hybrid ice sheet model
evFac = 1;
if oo.hybrid
evFac = (1 + (pp.c13*B2).*vv.F2);
end
Ub = sqrt((U./evFac).^2 + pp.U_rp.^2); %regularize, ensure > 0
if isequal(oo.slidinglaw, 2) %6a of Hewitt (2012)
p = pp.p; q = pp.q;
F = pp.c14 .*(N.^p .* Ub.^q);
F = F .* (Ub.^-1);
alpha = B2./F;
elseif isequal(oo.slidinglaw, 3) %6b of Hewitt (2012)
n = pp.n_Glen;
F = pp.c15 .* N .* (Ub./ (pp.c16.*Ub + pp.c17.*N.^n)).^(1/n);
F = F .* (Ub.^-1);
alpha = B2./F;
elseif isequal(oo.slidinglaw, 4) %6b of Hewitt (2012)
n = pp.n_Glen;
F1 = (B2/(pp.c20.*N)).^n;
F2 = pp.c16*Ub;
F3 = pp.c21*(N.^n);
alpha = (Ub./F1 - F2).*(1./F3);
end
%% Reshape
alpha = reshape(gg.S_h'*alpha,gg.nJ,gg.nI);
N = reshape(gg.S_h'*N,gg.nJ,gg.nI);
% Fill regions where N (effective pressure) is too low for an accurate
% estimate
eps = 1.5*pp.N_rp/pp.phi;
mask = N<eps;
mask(gg.next) =0;
alpha(mask) = NaN;
alpha = inpaint_nans(alpha,4);
%% Reshape for output
alpha = gg.S_h*alpha(:);
end
end