tplot_flowline.py - 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
---
tplot_flowline.py (3658B)
---
1 #!/usr/bin/env python3
2 # Copyright (C) 2011, 2014 Torsten Albrecht (torsten.albrecht@pik-potsdam.de)
3
4 import numpy as np
5 from pylab import figure, clf, hold
6 import matplotlib.pyplot as plt
7 try:
8 from netCDF4 import Dataset as NC
9 except:
10 print("netCDF4 is not installed!")
11 sys.exit(1)
12
13 pism_output = "flowline_result_50km.nc"
14 pism_output = "flowline_result_20km.nc"
15 #pism_output = "flowline_result_10km.nc"
16 #pism_output = "flowline_result_5km.nc"
17 #pism_output = "flowline_result_2km.nc"
18 #pism_output = "flowline_result_1km.nc"
19
20 # load the PISM output
21 try:
22 print("Loading PISM output from '%s'..." % (pism_output), end=' ')
23 infile = NC(pism_output, 'r')
24
25 except Exception:
26 print("""ERROR!\nSpecify NetCDF file from PISM run with -p.""")
27 exit(-1)
28
29 printfile = pism_output[:-3] + ".pdf"
30
31 thk = np.squeeze(infile.variables["thk"][:])
32 ubar = np.squeeze(infile.variables["ubar"][:])
33
34 # "typical constant ice parameter" as defined in the paper and in Van der
35 # Veen's "Fundamentals of Glacier Dynamics", 1999
36 U0 = 300 # m/yr
37 H0 = 600 # m
38 spa = 3.1556926e7
39 # spa=3600*365.25*24
40 C = 2.4511e-18
41 #C = (rho_ice * standard_gravity * (1.0 - rho_ice/rho_ocean) / (4 * B0))**3
42 Hcf = 250 # calving thickness
43
44 #grid and time
45 x = np.squeeze(infile.variables["x"][:])
46 dx = x[1] - x[0]
47 t = np.squeeze(infile.variables["time"][:])
48 t = int(np.round(t / spa))
49
50 dist = 50.0
51 idist = int(dist / (dx / 1000.0))
52 H = thk[1, idist:]
53 u = ubar[1, idist:]
54 xd = x[idist:] / 1000.0 - dist
55
56
57 # exact van der veen solution
58 xcf = (Hcf ** -4 - H0 ** -4) * U0 * H0 / (4 * C * spa) * 0.001 # find 250m front position on 5km grid
59 ucf = U0 * H0 * (xcf * 1000.0 * 4 * C / (U0 * H0 / spa) + H0 ** (-4)) ** (0.25) # exact (maximal) velocity at front
60
61 dxan = 1.0 # km
62 xan = np.arange(0, (x[len(x) - 1]) / 1000.0, dxan)
63 Han = np.zeros(len(xan))
64 Uan = np.zeros(len(xan))
65 for k in range(len(xan)):
66 if xan[k] <= xcf:
67 Han[k] = (k * dxan * 1000.0 * 4 * C / (U0 * H0 / spa) + H0 ** (-4)) ** (-0.25)
68 Uan[k] = U0 * H0 * (k * dxan * 1000.0 * 4 * C / (U0 * H0 / spa) + H0 ** (-4)) ** (0.25)
69
70
71 # imshow plot
72 plt.rcParams['font.size'] = 14
73 plt.rcParams['legend.fontsize'] = 14
74 plt.rcParams['contour.negative_linestyle'] = 'solid'
75
76 plt.clf()
77 figure(1, figsize=(10, 9))
78 clf()
79 hold(True)
80 plt.subplots_adjust(left=0.14, bottom=0.14, right=0.86, top=None, wspace=0.08,
81 hspace=0.15)
82 plt.clf()
83 title = 'steady state ice thickness and velocity for ' + str(int(dx / 1000.0)) + ' km after ' + str(t) + ' yrs'
84 plt.suptitle(title)
85 # plot
86 ax1 = plt.subplot(211)
87 ax1.plot(xd, H, '#B26415', xan, Han, '#4E238B', linewidth=2)
88 leg1 = ax1.legend(('H_pism in m', 'H_exact in m'), 'upper right', shadow=True)
89 plt.xlim([0, 400])
90 # plt.xticks([])
91 # plt.xticks(np.arange(0,50,10))
92 plt.ylim([0, 600])
93 plt.grid(True)
94
95 arrowtext2a = 'max(u_exact)=%d m/yr at %d km' % (int(ucf), int(xcf))
96 arrowtext2b = 'max(u_pism)=%d m/yr at %d km' % (int(np.max(u)), int(xd[np.argmax(u)]))
97 ax2 = plt.subplot(212)
98 ax2.plot(xd, u, '#B26415', xan, Uan, '#4E238B', linewidth=2)
99 leg2 = ax2.legend((arrowtext2b, arrowtext2a), 'lower right', shadow=True)
100 plt.xlim([0, 400])
101 # plt.xticks([])
102 plt.ylim([0, 900])
103 plt.yticks(np.arange(0, 1000, 200))
104 # ax2.annotate(arrowtext2a, xy=(int(x[np.argmax(Uexact)]),int(np.max(Uexact))), xycoords='data',
105 # xytext=(40, 20), textcoords='offset points',
106 # arrowprops=dict(arrowstyle="->"))
107 #
108 # ax2.annotate(arrowtext2b, xy=(int(x[np.argmax(u)]),int(np.max(u))), xycoords='data',
109 # xytext=(40, 0), textcoords='offset points', arrowprops=dict(arrowstyle="->"))
110 plt.grid(True)
111 plt.savefig(printfile, format='pdf')
112 plt.show()