tcircular_ice_sheet.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
       ---
       tcircular_ice_sheet.py (2533B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # Copyright (C) 2012, 2013, 2014 Ricarda Winkelmann, Torsten Albrecht,
            4 # Ed Bueler, and Constantine Khroulev
            5 
            6 import numpy as np
            7 import scipy.optimize as opt
            8 import PISMNC
            9 import piktests_utils
           10 
           11 # command line arguments
           12 options = piktests_utils.process_options("circular_withshelf.nc",
           13                                          domain_size=3600)
           14 p = piktests_utils.Parameters()
           15 
           16 dx, dy, x, y = piktests_utils.create_grid(options)
           17 
           18 # Compute the calving front radius
           19 x_cf = 1000.0 * options.domain_size / 3  # at 1200 km
           20 y_cf = 1000.0 * options.domain_size / 3  # at 1200 km
           21 r_cf = np.sqrt(x_cf ** 2 + y_cf ** 2)  # calving front position in m
           22 
           23 # create arrays which will go in the output file
           24 thk = np.zeros((options.My, options.Mx))  # sheet/shelf thickness
           25 bed = np.zeros_like(thk)                 # bedrock surface elevation
           26 Ts = np.zeros_like(thk) + p.air_temperature
           27 accum = np.zeros_like(thk) + p.accumulation_rate * p.rho_ice
           28 
           29 
           30 def MISMIP_bed(r):
           31     s = r / 1000.0
           32     n = 729.0
           33     m2 = -2184.80 / (750.0) ** 2
           34     m4 = 1031.72 / (750.0) ** 4
           35     m6 = -151.72 / (750.0) ** 6
           36     return n + m2 * s ** 2 + m4 * s ** 4 + m6 * s ** 6  # in m
           37 
           38 
           39 def MISMIP_thk(r):
           40     thk_cf = 200.0             # ice thickness at calving front
           41     thk_max = 4000.0            # maximal ice thickness in m
           42     a = -(thk_cf - thk_max) / (r_cf) ** 4
           43     b = 2 * (thk_cf - thk_max) / (r_cf) ** 2
           44     c = thk_max
           45     return a * r ** 4 + b * r ** 2 + c
           46 
           47 
           48 # bedrock and ice thickness
           49 for j in range(options.My):
           50     for i in range(options.Mx):
           51         radius = np.sqrt(x[i] ** 2 + y[j] ** 2)  # radius in m
           52 
           53         # set bedrock as in MISMIP experiment
           54         bed[j, i] = MISMIP_bed(radius)
           55 
           56         # set thickness
           57         if radius <= r_cf:
           58             thk[j, i] = MISMIP_thk(radius)
           59 
           60 # clip bed topography
           61 bed[bed < p.topg_min] = p.topg_min
           62 
           63 # Compute the grounding line radius
           64 
           65 
           66 def f(x):
           67     "floatation criterion: rho_ice/rho_ocean * thk + bed = 0"
           68     return (p.rho_ice / p.rho_ocean) * MISMIP_thk(x) + MISMIP_bed(x)
           69 
           70 
           71 r_gl = opt.bisect(f, 0, r_cf)
           72 print("grounding line radius = %.2f km" % (r_gl / 1000.0))
           73 
           74 ncfile = PISMNC.PISMDataset(options.output_filename, 'w', format='NETCDF3_CLASSIC')
           75 
           76 piktests_utils.prepare_output_file(ncfile, x, y, include_vel_bc=False)
           77 
           78 variables = {"thk": thk,
           79              "topg": bed,
           80              "ice_surface_temp": Ts,
           81              "climatic_mass_balance": accum}
           82 
           83 piktests_utils.write_data(ncfile, variables)
           84 
           85 ncfile.close()
           86 print("Successfully created %s" % options.output_filename)