plot_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()