ttest_iceberg_removal.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
       ---
       ttest_iceberg_removal.py (2582B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # Copyright (C) 2012, 2013, 2014, 2015 Ricarda Winkelmann, Torsten Albrecht,
            4 # Ed Bueler, and Constantine Khroulev
            5 
            6 import numpy as np
            7 import PISMNC
            8 import piktests_utils
            9 
           10 # command line arguments
           11 options = piktests_utils.process_options("test_iceberg_removal.nc",
           12                                          domain_size=1000.0)
           13 p = piktests_utils.Parameters()
           14 
           15 # create arrays which will go in output
           16 thk = np.zeros((options.My, options.Mx))  # sheet/shelf thickness
           17 bed = np.zeros_like(thk)                 # bedrock surface elevation
           18 accum = np.zeros_like(thk)                 # accumulation/ ablation
           19 Ts = np.zeros_like(thk) + p.air_temperature
           20 
           21 bc_mask = np.zeros_like(thk)
           22 ubar = np.zeros_like(thk)
           23 vbar = np.zeros_like(thk)
           24 
           25 dx, dy, x, y = piktests_utils.create_grid(options)
           26 
           27 
           28 def R(x, y):
           29     return np.maximum(np.abs(x), np.abs(y))
           30 
           31 
           32 def C(x, y):
           33     return np.sqrt(x ** 2 + y ** 2)
           34 
           35 
           36 if options.square:
           37     rad = R
           38 else:
           39     rad = C
           40 
           41 # bedrock and ice thickness
           42 for j in range(options.My):
           43     for i in range(options.Mx):
           44         radius = rad(x[i], y[j])
           45 
           46         # grounded
           47         if radius <= p.r_gl and x[i] < 0:
           48             bed[j, i] = 100.0
           49             thk[j, i] = p.H0
           50 
           51         # ice shelf
           52         if radius <= p.r_cf and radius > p.r_gl and options.shelf == True:
           53             thk[j, i] = (4.0 * p.C / p.Q0 * (radius - p.r_gl) + 1 / p.H0 ** 4) ** (-0.25)
           54 
           55         # ocean (shelf and elsewhere)
           56         if radius >= p.r_gl or x[i] >= 0:
           57             accum[j, i] = p.accumulation_rate * p.rho_ice
           58             bed[j, i] = p.topg_min
           59 
           60 # cap ice thickness
           61 thk[thk > p.H0] = p.H0
           62 
           63 # set values and locations of Dirichlet boundary conditions
           64 for j in range(options.My):
           65     for i in range(options.Mx):
           66         radius = rad(x[i], y[j])
           67         width = dx * 3
           68 
           69         if radius <= p.r_gl - width and x[i] < 0:
           70             bc_mask[j, i] = 1.0
           71         elif radius <= p.r_gl and x[i] < 0:
           72             bc_mask[j, i] = 1.0
           73             ubar[j, i] = p.vel_bc * (x[i] / radius)
           74             vbar[j, i] = p.vel_bc * (y[j] / radius)
           75         else:
           76             bc_mask[j, i] = 0.0
           77 
           78 ncfile = PISMNC.PISMDataset(options.output_filename, 'w', format='NETCDF3_CLASSIC')
           79 piktests_utils.prepare_output_file(ncfile, x, y)
           80 
           81 variables = {"thk": thk,
           82              "topg": bed,
           83              "ice_surface_temp": Ts,
           84              "climatic_mass_balance": accum,
           85              "bc_mask": bc_mask,
           86              "u_ssa_bc": ubar,
           87              "v_ssa_bc": vbar}
           88 
           89 piktests_utils.write_data(ncfile, variables)
           90 ncfile.close()
           91 
           92 print("Successfully created %s" % options.output_filename)