tbuild_tiny.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
---
tbuild_tiny.py (3689B)
---
1 #! /usr/bin/env python3
2 #
3
4 import PISM
5
6 PISM.set_abort_on_sigint(True)
7 context = PISM.Context()
8 config = PISM.Context().config
9
10 config.set_string("output.file_name", "tiny.nc")
11
12 # Default constants that may get overridden later.
13
14 Ly = 25e3 # 25 km
15 Lx = 50e3 # 50 km
16 Lz = 4000
17 My = 13
18 Mx = 23
19 Mz = 21
20
21 H0 = 60. # ice thickness at cliff
22 alpha = 0.008 # constant surface slope
23 Lext = 15e3 # width of strip beyond cliff
24 Lstream_x = 50e3
25 Lstream_y = 30e3
26
27 Hext = 0. # m ice thickeness beyond the cliff
28
29 tauc_hi = 2e6 # Pa
30 tauc_lo = 1e4 # Pa
31 tauc_free_bedrock = config.get_number('basal_yield_stress.ice_free_bedrock')
32
33 EC = PISM.EnthalpyConverter(PISM.Context().config)
34 enth0 = EC.enthalpy(273.15, 0.01, 0) # 0.01 water fraction
35 bed0 = 0
36
37
38 def geometry(x, y):
39 x0 = -Lx + Lext
40 if x < x0:
41 return (0, Hext)
42 return (0, H0 + alpha * (x - x0))
43
44
45 def stream_tauc(x, y):
46 x0 = -Lx + Lext
47 if x < x0:
48 return tauc_free_bedrock
49 if x < x0 + Lstream_x:
50 if abs(y) < Lstream_y / 2:
51 return tauc_lo
52 return tauc_hi
53
54
55 # The main code for a run follows:
56 if __name__ == '__main__':
57
58 config = PISM.Context().config
59
60 config.set_number("grid.Mx", Mx, PISM.CONFIG_DEFAULT)
61 config.set_number("grid.My", My, PISM.CONFIG_DEFAULT)
62
63 # Build the grid.
64 p = PISM.GridParameters(config)
65 p.Mx = int(config.get_number("grid.Mx"))
66 p.My = int(config.get_number("grid.My"))
67 p.Lx = Lx
68 p.Ly = Ly
69 z = PISM.IceGrid.compute_vertical_levels(Lz, Mz, PISM.EQUAL, 4.0)
70 p.z = PISM.DoubleVector(z)
71 p.ownership_ranges_from_options(context.size)
72 p.registration = PISM.CELL_CORNER
73 p.periodicity = PISM.NOT_PERIODIC
74 grid = PISM.IceGrid(context.ctx, p)
75
76 vecs = PISM.model.ModelVecs(grid.variables())
77 vecs.add(PISM.model.createIceSurfaceVec(grid))
78 vecs.add(PISM.model.createIceThicknessVec(grid))
79 vecs.add(PISM.model.createBedrockElevationVec(grid))
80 vecs.add(PISM.model.createYieldStressVec(grid), 'tauc')
81 vecs.add(PISM.model.createEnthalpyVec(grid), 'enthalpy')
82 vecs.add(PISM.model.createIceMaskVec(grid))
83 vecs.add(PISM.model.createNoModelMaskVec(grid), 'no_model_mask')
84 vecs.add(PISM.model.create2dVelocityVec(grid, name='_ssa_bc', desc='SSA Dirichlet BC'))
85 vecs.add(PISM.model.createSeaLevelVec(grid))
86
87 # Set constant coefficients.
88 vecs.enthalpy.set(enth0)
89
90 # Build the continent
91 bed = vecs.bedrock_altitude
92 thickness = vecs.land_ice_thickness
93 sea_level = vecs.sea_level
94
95 with PISM.vec.Access(comm=[bed, thickness, sea_level]):
96 for (i, j) in grid.points():
97 x = grid.x(i)
98 y = grid.y(j)
99 (b, t) = geometry(x, y)
100 bed[i, j] = b
101 thickness[i, j] = t
102 sea_level[i, j] = 0.0
103
104 # Compute mask and surface elevation from geometry variables.
105 gc = PISM.GeometryCalculator(grid.ctx().config())
106 gc.compute(sea_level, bed, thickness, vecs.mask, vecs.surface_altitude)
107
108 tauc = vecs.tauc
109 with PISM.vec.Access(comm=tauc):
110 for (i, j) in grid.points():
111 tauc[i, j] = stream_tauc(grid.x(i), grid.y(j))
112
113 vecs.vel_ssa_bc.set(0.0)
114 no_model_mask = vecs.no_model_mask
115 no_model_mask.set(0)
116 with PISM.vec.Access(comm=[no_model_mask]):
117 for (i, j) in grid.points():
118 if (i == 0) or (i == grid.Mx() - 1) or (j == 0) or (j == grid.My() - 1):
119 no_model_mask[i, j] = 1
120
121 output_filename = config.get_string("output.file_name")
122 pio = PISM.util.prepare_output(output_filename)
123 pio.close()
124 vecs.writeall(output_filename)
125 PISM.util.writeProvenance(output_filename)