pism-exp-gsw

commit 58a7957458f2651dfd59dc888d30cb0d804fa062
parent 5fdab4dccadedc71c81a54b3e472a05e241658e7
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  6 Dec 2021 16:21:53 +0100

fix plotting of ice geometry

Diffstat:
M
plot-time-series.py
 | 
46
+++++++++++++++++++++++++---------------------

1 file changed, 25 insertions(+), 21 deletions(-)
diff --git a/plot-time-series.py b/plot-time-series.py
@@ -77,37 +77,41 @@ def plot_profile(in_file, out_file):
     if out_file is None:
         out_file = os.path.splitext(in_file)[0] + "-profile.pdf"
 
-    mask = read(in_file, 'mask')
-    usurf = read(in_file, 'usurf')
-    topg = read(in_file, 'topg')
-    thk = read(in_file, 'thk')
-    x = read(in_file, 'x')
-
-    # mask out ice-free areas
-    usurf = np.ma.array(usurf, mask=mask == 4)
-
-    # convert x to kilometers
-    x /= 1e3
+    steps = read(in_file, 'thk').shape[0]
 
     figure(1)
     ax = subplot(111)
-    for i in range(1, thk.shape[0]):
+    for i in range(1, steps):
+
+        mask = read(in_file, 'mask')[i]
+        usurf = read(in_file, 'usurf')[i]
+        topg = read(in_file, 'topg')[i]
+        lsrf = topg.copy()
+        thk = read(in_file, 'thk')[i]
+        sea_lvl = read(in_file, 'sea_level')[i]
+        x = read(in_file, 'x')
+
+        # mask out ice-free areas
+        usurf = np.ma.array(usurf, mask=mask == 4)
+
+        # convert x to kilometers
+        x /= 1e3
+
         # compute the lower surface elevation
-        lsrf = topg[i].copy()
-        thk_ = thk[i].copy()
-        lsrf[mask[i] == 3] = -MISMIP.rho_i() / MISMIP.rho_w() * thk_[mask[i] == 3]
-        lsrf = np.ma.array(lsrf, mask=mask[i] == 4)
+        lsrf[mask == 3] = usurf[mask == 3] - thk[mask == 3]
+        lsrf = np.ma.array(lsrf, mask=mask == 4)
     	# modeled grounding line position
-        #xg_PISM = find_grounding_line(x, lsrf, thk_, mask[i])
+        #xg_PISM = find_grounding_line(x, lsrf, thk, mask)
         #plot(x, np.zeros_like(x), ls='dotted', color='red')
-        icecolor = cm.cool(i / thk.shape[0])
-        plot(x, topg[i], color='black')
-        plot(x, usurf[i], color=icecolor, label='{}'.format(i))
+        icecolor = cm.cool(i / steps)
+        plot(x, topg, color='black')
+        plot(x, usurf, color=icecolor, label='{}'.format(i))
         plot(x, lsrf, color=icecolor)
+
         xlabel('distance from the divide, km')
         ylabel('elevation, m')
 
-        _, _, ymin, ymax = axis(xmin=0, xmax=x.max())
+        #_, _, ymin, ymax = axis(xmin=0, xmax=x.max())
         #_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max())
 
         #vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')