tplot-till-evolution.py - pism-exp-gsw - ice stream and sediment transport experiments
 (HTM) git clone git://src.adamsgaard.dk/pism-exp-gsw
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tplot-till-evolution.py (3787B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 from pylab import figure, subplots, plot, xlabel, ylabel, title, axis, vlines, savefig, text, tight_layout, cm, legend
            4 from sys import exit
            5 
            6 import MISMIP
            7 
            8 import numpy as np
            9 from optparse import OptionParser
           10 import os.path
           11 
           12 try:
           13     from netCDF4 import Dataset as NC
           14 except:
           15     print("netCDF4 is not installed!")
           16     sys.exit(1)
           17 
           18 def process_options():
           19     "Process command-line options and arguments."
           20     parser = OptionParser()
           21     parser.usage = "%prog <input files> [options]"
           22     parser.description = "Plots the ice flux as a function of the distance from the divide."
           23     parser.add_option("-o", "--output", dest="output", type="string",
           24                       help="Output image file name (e.g. -o foo.png)")
           25 
           26     opts, args = parser.parse_args()
           27 
           28     if len(args) == 0:
           29         print("ERROR: An input file is requied.")
           30         exit(0)
           31 
           32     if len(args) > 1 and opts.output:
           33         print("More than one input file given. Ignoring the -o option...\n")
           34         opts.output = None
           35 
           36     return args, opts.output, opts
           37 
           38 
           39 def read(filename, name):
           40     "Read a variable and extract the middle row."
           41     nc = NC(filename)
           42 
           43     try:
           44         var = nc.variables[name][:]
           45     except:
           46         print("ERROR: Variable '%s' not present in '%s'" % (name, filename))
           47         exit(1)
           48 
           49     return var
           50 
           51 
           52 def find_grounding_line(x, topg, thk, mask):
           53     "Find the modeled grounding line position."
           54     # "positive" parts of x, topg, thk, mask
           55     topg = topg[x > 0]
           56     thk = thk[x > 0]
           57     mask = mask[x > 0]
           58     x = x[x > 0]                        # this should go last
           59 
           60     def f(j):
           61         "See equation (7) in Pattyn et al, 'Role of transition zones in marine ice sheet dynamics', 2005."
           62         z_sl = 0
           63         return (z_sl - topg[j]) * MISMIP.rho_w() / (MISMIP.rho_i() * thk[j])
           64 
           65     for j in range(x.size):
           66         if mask[j] == 2 and mask[j + 1] == 3:  # j is grounded, j+1 floating
           67             nabla_f = (f(j + 1) - f(j)) / (x[j + 1] - x[j])
           68 
           69             # See equation (8) in Pattyn et al
           70             return (1.0 - f(j) + nabla_f * x[j]) / nabla_f
           71 
           72     raise Exception("Can't find the grounding line")
           73 
           74 
           75 def plot_profile(in_file, out_file):
           76 
           77     if out_file is None:
           78         out_file = os.path.splitext(in_file)[0] + "-till-evol.pdf"
           79 
           80     steps = read(in_file, 'thk').shape[0]
           81 
           82     fig, ax = subplots(1, 1, sharex=True, figsize=[6, 3])
           83     for i in range(1, steps):
           84 
           85         mask = read(in_file, 'mask')[i]
           86         uvelbase = read(in_file, 'uvelbase')[i]
           87         utillflux = read(in_file, 'utillflux')[i]
           88         till_deposit = read(in_file, 'tilldeposit')[i]
           89         x = read(in_file, 'x')
           90 
           91         # convert x to kilometers
           92         x /= 1e3
           93 
           94             # modeled grounding line position
           95         #xg_PISM = find_grounding_line(x, lsrf, thk, mask)
           96         #plot(x, np.zeros_like(x), ls='dotted', color='red')
           97         icecolor = cm.cool(i / steps)
           98         #ax[0].plot(x, uvelbase, color=icecolor)
           99         #ax[1].plot(x, utillflux, color=icecolor) #, label='{}'.format(i))
          100         ax.plot(x, till_deposit, color=icecolor)
          101 
          102         ax.set_xlabel('distance from the divide, km')
          103         #ax[0].set_ylabel('$v_{SSA}$, m/a')
          104         #ax[1].set_ylabel('$q_{t,x}$, m$^2$/a')
          105         ax.set_ylabel('$\Delta b$, m')
          106 
          107         #_, _, ymin, ymax = axis(xmin=0, xmax=x.max())
          108         _, _, ymin, ymax = axis(xmin=950, xmax=1150)
          109         #_, _, ymin, ymax = axis(xmin=x.min(), xmax=x.max())
          110 
          111         #vlines(xg / 1e3, ymin, ymax, linestyles='dashed', color='black')
          112         #vlines(xg_PISM / 1e3, ymin, ymax, linestyles='dashed', color='red')
          113 
          114     #legend()
          115     fig.tight_layout()
          116     print("Saving '%s'...\n" % out_file)
          117     savefig(out_file)
          118 
          119 if __name__ == "__main__":
          120     args, out_file, opts = process_options()
          121 
          122     for in_file in args:
          123