tgroundingline.m - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tgroundingline.m (2822B)
---
1 function groundingline(filename)
2 % GROUNDINGLINE Compute grounding line from fields in input file.
3 % Uses fields x,y,thk,topg and applies flotation criterion.
4 % Example:
5 % $ cd pism-dev/examples/eisross
6 % $ ./quickstart.sh # generate ross.nc; run irrelevant
7 % $ cp ../../util/groundingline.m .
8 % $ matlab
9 % >> groundingline('ross.nc')
10
11 FIXME: use Matlab's netcdf
12
13 rhoi = 910.0;
14 rhow = 1028.0;
15 grounded = rhoi * thk + rhow * topg; % if this is positive then
16 % rho_i H > - rho_w b
17 % so ice is grounded
18
19 % show "raw" contour:
20 figure
21 C = contour(x/1000 , y/1000, grounded, [0, 0]); % specify zero contour
22 axis equal, xlabel('x (km)'), ylabel('y (km)')
23
24 % now inspect C, which is 2 x N, to get polygonal grounding line
25 % see "help contourc" for format
26 % note contour C has several closed loops, and we grab the *biggest* one;
27 % the rest are islands and lakes, etc.
28 jpoly = 1;
29 Npoly = C(2,1);
30 j = 1;
31 while j + C(2,j) < size(C,2)
32 j = j + C(2,j) + 1;
33 if Npoly < C(2,j)
34 jpoly = j;
35 Npoly = C(2,j);
36 end
37 end
38
39 % we know what part of C is the main contour
40 xpoly = C(1,jpoly+1:jpoly+Npoly); % note units of xpoly,ypoly are km
41 ypoly = C(2,jpoly+1:jpoly+Npoly);
42
43 figure
44 plot(xpoly, ypoly)
45 axis equal, xlabel('x (km)'), ylabel('y (km)')
46
47 % make polygonal grounding line M times as fine before interpolating;
48 % xpf = "x polygonal fine", so to speak
49 M = 10;
50 finelength = (length(xpoly)-1) * M + 1;
51 xpf = zeros(1,finelength);
52 ypf = zeros(1,finelength);
53 for j=1:length(xpoly)-1
54 xpf((j-1)*M+1:j*M) = linspace(xpoly(j),xpoly(j+1),M);
55 ypf((j-1)*M+1:j*M) = linspace(ypoly(j),ypoly(j+1),M);
56 end
57 xpf(end) = xpoly(end);
58 ypf(end) = ypoly(end);
59 figure
60 plot(xpf, ypf), axis square
61
62 % compute arclength sf (= s fine) along grounding line
63 sf = zeros(1,finelength);
64 sf(2:finelength) = cumsum(sqrt(diff(xpf).^2 + diff(ypf).^2));
65 L = max(sf);
66 disp(['computed length of grounding line is L = ' num2str(L) ' km'])
67
68 % interpolate to get profile; see "help interp2"
69 thkf = interp2(x/1000,y/1000,thk,xpf,ypf,'linear');
70
71 % plot profile
72 figure
73 r = rhoi/rhow;
74 plot(sf, (1-r) * thkf, sf, -r * thkf)
75 xlabel('arclength (km)'), ylabel('profile in (m)')
76 xtloc = linspace(0,L,6);
77 set(gca,'XTick',xtloc)
78 set(gca,'XTickLabel',['0 '; '0.2 L'; '0.4 L'; '0.6 L'; '0.8 L'; 'L '])
79
80 % FIXME profile detail along Siple coast and Ross area
81 figure
82 sfe = sf(sf>0.8*L);
83 thkfe = thkf(sf>0.8*L);
84 plot(sfe, (1-r) * thkfe, sfe, -r * thkfe)
85 xlabel('arclength (km)'), ylabel('profile in (m)')
86 set(gca,'XTick',[0.8*L, L])
87 set(gca,'XTickLabel',['0.8 L'; 'L '])
88
89 % FIXME go back to map-plane curve and label locations
90 figure
91 hold on
92 text(xpf(1),ypf(1),'s = 0 = L')
93 for k=1:4
94 i = max(find(sf< (k/5)*L));
95 text(xpf(i),ypf(i),['s = ' num2str(k/5) 'L'])
96 end
97 hold off
98