commit 548a0043c892eee76690872d37ae27deaf5dd7e5
parent e4dac5432aa0219eae7c6a5b35a6751018a6b55f
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 29 Nov 2021 11:46:44 +0100
add sequential plot of grounding line migration
Diffstat:
2 files changed, 129 insertions(+), 1 deletion(-)
diff --git a/Makefile b/Makefile
@@ -3,7 +3,7 @@ SLSERIES=sealvl.nc
all: \
ABC1_1a_M1_A1-flux.pdf\
- ex_1a_M1_A1_1d.nc\
+ evol_1a_M1_A1_1d.pdf\
#ABC1_1a_M3_A1-flux.pdf\
#ABC1_3a_M1_A1-flux.pdf\
#ABC1_3a_M3_A1-flux.pdf\
@@ -20,6 +20,9 @@ ABC1_3a_M1_A1-flux.pdf: ABC1_3a_M1_A1.nc plot.py
ABC1_3a_M3_A1-flux.pdf: ABC1_3a_M3_A1.nc plot.py
./plot.py ABC1_3a_M3_A1.nc
+evol_1a_M1_A1_1d.pdf: ex_1a_M1_A1_1d.nc
+ ./plot-time-series.py ex_1a_M1_A1_1d.nc
+
ex_1a_M1_A1_1d.nc: ex_ABC1_1a_M1_A1.nc
flowline.py -o $@ --collapse -d y ex_ABC1_1a_M1_A1.nc
diff --git a/plot-time-series.py b/plot-time-series.py
@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+
+from pylab import figure, subplot, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend
+from sys import exit
+
+import MISMIP
+
+import numpy as np
+from optparse import OptionParser
+import os.path
+
+try:
+ from netCDF4 import Dataset as NC
+except:
+ print("netCDF4 is not installed!")
+ sys.exit(1)
+
+def process_options():
+ "Process command-line options and arguments."
+ parser = OptionParser()
+ parser.usage = "%prog <input files> [options]"
+ parser.description = "Plots the ice flux as a function of the distance from the divide."
+ parser.add_option("-o", "--output", dest="output", type="string",
+ help="Output image file name (e.g. -o foo.png)")
+
+ opts, args = parser.parse_args()
+
+ if len(args) == 0:
+ print("ERROR: An input file is requied.")
+ exit(0)
+
+ if len(args) > 1 and opts.output:
+ print("More than one input file given. Ignoring the -o option...\n")
+ opts.output = None
+
+ return args, opts.output, opts
+
+
+def read(filename, name):
+ "Read a variable and extract the middle row."
+ nc = NC(filename)
+
+ try:
+ var = nc.variables[name][:]
+ except:
+ print("ERROR: Variable '%s' not present in '%s'" % (name, filename))
+ exit(1)
+
+ return var
+
+
+def find_grounding_line(x, topg, thk, mask):
+ "Find the modeled grounding line position."
+ # "positive" parts of x, topg, thk, mask
+ topg = topg[x > 0]
+ thk = thk[x > 0]
+ mask = mask[x > 0]
+ x = x[x > 0] # this should go last
+
+ def f(j):
+ "See equation (7) in Pattyn et al, 'Role of transition zones in marine ice sheet dynamics', 2005."
+ z_sl = 0
+ return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk[j])
+
+ for j in range(x.size):
+ if mask[j] == 2 and mask[j + 1] == 3: # j is grounded, j+1 floating
+ nabla_f = (f(j + 1) - f(j)) / (x[j + 1] - x[j])
+
+ # See equation (8) in Pattyn et al
+ return (1.0 - f(j) + nabla_f * x[j]) / nabla_f
+
+ raise Exception("Can't find the grounding line")
+
+
+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)
+
+ # compute the lower surface elevation
+ lsrf = topg.copy()
+ lsrf[mask == 3] = -MISMIP.rho_i() / MISMIP.rho_w() * thk[mask == 3]
+ lsrf = np.ma.array(lsrf, mask=mask == 4)
+
+ # convert x to kilometers
+ x /= 1e3
+
+ figure(1)
+ ax = subplot(111)
+ for i in range(0, thk.shape[0]):
+ # modeled grounding line position
+ #xg_PISM = find_grounding_line(x, topg[i], thk[i], mask[i])
+ #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))
+ plot(x, lsrf[i], color=icecolor)
+ xlabel('distance from the divide, km')
+ ylabel('elevation, m')
+
+ _, _, 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')
+ #vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red')
+
+ legend()
+ print("Saving '%s'...\n" % out_file)
+ savefig(out_file)
+
+if __name__ == "__main__":
+ args, out_file, opts = process_options()
+
+ for in_file in args:
+ plot_profile(in_file, out_file)+
\ No newline at end of file
