tcircular_dirichlet.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_dirichlet.py (2606B)
---
1 #!/usr/bin/env python3
2
3 # Copyright (C) 2012, 2013, 2014, 2015, 2017 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("circular_dirichlet.nc",
12 domain_size=1000.0)
13 p = piktests_utils.Parameters()
14
15 dx, dy, x, y = piktests_utils.create_grid(options)
16
17 xx, yy = np.meshgrid(x, y)
18 radius = np.sqrt(xx ** 2 + yy ** 2)
19 # remove the singularity (does not affect the result):
20 radius[radius == 0] = 1e-16
21
22 # Ice thickness
23 thk = np.zeros((options.My, options.Mx)) # sheet/shelf thickness
24 if options.shelf:
25 thk[radius > p.r_gl] = (4.0 * p.C / p.Q0 * (radius[radius > p.r_gl] - p.r_gl) + 1 / p.H0 ** 4) ** (-0.25)
26 thk[radius >= p.r_cf] = 0.0
27 # cap ice thickness
28 thk[thk > p.H0] = p.H0
29
30 thk[radius <= p.r_gl] = p.H0
31 thk[radius <= p.r_gl - 0.5 * (dx + dy)] = 0.0
32 else:
33 thk[radius <= p.r_gl] = p.H0
34
35 # Ice thickness threshold for "calving at a threshold thickness"
36 thk_threshold = np.zeros((options.My, options.Mx)) # sheet/shelf thickness
37 thk_threshold[:] = 500.0
38 thk_threshold[xx < 0.0] = 200.0
39
40 # Bed topography
41 bed = np.zeros_like(thk) # bedrock surface elevation
42 bed[radius <= p.r_gl] = -p.H0 * (p.rho_ice / p.rho_ocean) + 1.0
43 bed[radius > p.r_gl] = p.topg_min
44
45 # Surface mass balance
46 accum = np.zeros_like(thk) # accumulation/ ablation
47 accum[radius > p.r_gl] = p.accumulation_rate * p.rho_ice # convert to [kg m-2 s-1]
48
49 # Surface temperature (irrelevant)
50 Ts = np.zeros_like(thk) + p.air_temperature
51
52 # Dirichlet B.C locations
53 bc_mask = np.zeros_like(thk)
54 bc_mask[radius <= p.r_gl] = 1
55
56 # SSA velocity Dirichlet B.C.
57 ubar = np.zeros_like(thk)
58 ubar[bc_mask == 1] = p.vel_bc * (xx[radius <= p.r_gl] / radius[radius <= p.r_gl])
59 ubar[bc_mask == 0] = 0
60
61 vbar = np.zeros_like(thk)
62 vbar[bc_mask == 1] = p.vel_bc * (yy[radius <= p.r_gl] / radius[radius <= p.r_gl])
63 vbar[bc_mask == 0] = 0
64
65 ncfile = PISMNC.PISMDataset(options.output_filename, 'w', format='NETCDF3_CLASSIC')
66 piktests_utils.prepare_output_file(ncfile, x, y)
67
68 variables = {"thk": thk,
69 "topg": bed,
70 "ice_surface_temp": Ts,
71 "climatic_mass_balance": accum,
72 "bc_mask": bc_mask,
73 "u_ssa_bc": ubar,
74 "v_ssa_bc": vbar,
75 "calving_threshold": thk_threshold}
76
77 piktests_utils.write_data(ncfile, variables)
78
79 ncfile.variables["calving_threshold"].units = "m"
80 ncfile.close()
81
82 print("Successfully created %s" % options.output_filename)