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)