tsia_forward.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 --- tsia_forward.py (3592B) --- 1 #!/usr/bin/env python3 2 # 3 # Copyright (C) 2011, 2012, 2014, 2015, 2016, 2017, 2018, 2019, 2020 David Maxwell and Constantine Khroulev 4 # 5 # This file is part of PISM. 6 # 7 # PISM is free software; you can redistribute it and/or modify it under the 8 # terms of the GNU General Public License as published by the Free Software 9 # Foundation; either version 3 of the License, or (at your option) any later 10 # version. 11 # 12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY 13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15 # details. 16 # 17 # You should have received a copy of the GNU General Public License 18 # along with PISM; if not, write to the Free Software 19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 20 21 import PISM 22 from petsc4py import PETSc 23 import os 24 25 context = PISM.Context() 26 ctx = context.ctx 27 config = context.config 28 29 PISM.set_abort_on_sigint(True) 30 31 usage = """\ 32 sia_forward.py -i IN.nc [-o file.nc] 33 where: 34 -i IN.nc is input file in NetCDF format: contains PISM-written model state 35 notes: 36 * -i is required 37 """ 38 39 PISM.show_usage_check_req_opts(ctx.log(), "sia_forward.py", ["-i"], usage) 40 41 input_filename = config.get_string("input.file") 42 if len(input_filename) == 0: 43 import sys 44 sys.exit(1) 45 46 config.set_string("output.file_name", "sia_" + os.path.basename(input_filename), PISM.CONFIG_DEFAULT) 47 48 output_file = config.get_string("output.file_name") 49 is_regional = PISM.OptionBool("-regional", "Compute SIA using regional model semantics") 50 51 registration = PISM.CELL_CENTER 52 if is_regional: 53 registration = PISM.CELL_CORNER 54 55 input_file = PISM.File(ctx.com(), input_filename, PISM.PISM_NETCDF3, PISM.PISM_READONLY) 56 grid = PISM.IceGrid.FromFile(ctx, input_file, "enthalpy", registration) 57 58 config.set_flag("basal_resistance.pseudo_plastic.enabled", False) 59 60 enthalpyconverter = PISM.EnthalpyConverter(config) 61 62 modeldata = PISM.model.ModelData(grid) 63 modeldata.setPhysics(enthalpyconverter) 64 65 vecs = modeldata.vecs 66 67 vecs.add(PISM.model.createIceSurfaceVec(grid)) 68 vecs.add(PISM.model.createIceThicknessVec(grid)) 69 vecs.add(PISM.model.createBedrockElevationVec(grid)) 70 vecs.add(PISM.model.createEnthalpyVec(grid)) 71 vecs.add(PISM.model.createIceMaskVec(grid)) 72 73 # Read in the PISM state variables that are used directly in the SSA solver 74 for v in [vecs.thk, vecs.topg, vecs.enthalpy]: 75 v.regrid(input_file, critical=True) 76 77 # variables mask and surface are computed from the geometry previously read 78 sea_level = PISM.model.createSeaLevelVec(grid) 79 sea_level.set(0.0) 80 gc = PISM.GeometryCalculator(config) 81 gc.compute(sea_level, vecs.topg, vecs.thk, vecs.mask, vecs.surface_altitude) 82 83 # If running in regional mode, load in regional variables 84 if is_regional: 85 vecs.add(PISM.model.createNoModelMask(grid)) 86 vecs.no_model_mask.regrid(input_file, critical=True) 87 88 if PISM.util.fileHasVariable(input_file, 'usurfstore'): 89 vecs.add(PISM.model.createIceSurfaceStoreVec(grid)) 90 vecs.usurfstore.regrid(input_file, critical=True) 91 else: 92 vecs.add(vecs.surface, 'usurfstore') 93 94 solver = PISM.SIAFD_Regional 95 else: 96 solver = PISM.SIAFD 97 98 PISM.verbPrintf(2, context.com, "* Computing SIA velocities...\n") 99 vel_sia = PISM.sia.computeSIASurfaceVelocities(modeldata, siasolver=solver) 100 101 PISM.verbPrintf(2, context.com, "* Saving results to %s...\n" % output_file) 102 pio = PISM.util.prepare_output(output_file) 103 pio.close() 104 105 # Save time & command line & results 106 PISM.util.writeProvenance(output_file) 107 vel_sia.write(output_file)