tpism_python.py - pism-exp-gsw - ice stream and sediment transport experiments
 (HTM) git clone git://src.adamsgaard.dk/pism-exp-gsw
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tpism_python.py (2829B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 # Copyright (C) 2009-2015, 2018 the PISM Authors
            4 
            5 # @package pism_python
            6 # \author the PISM authors
            7 # \brief Creates "from scratch" a boring dataset with the right format
            8 # to use as a PISM bootstrapping file.
            9 # \details Example use of Python for this purpose.
           10 #
           11 # Usage, including a minimal PISM call to bootstrap from this file:
           12 #
           13 # \verbatim $ pism_python.py  # creates foo.nc \endverbatim
           14 # \verbatim $ pismr -i foo.nc -bootstrap -Mx 41 -My 41 -Mz 21 -Lz 4000 -Mbz 5 -Lbz 500 -y 1 \endverbatim
           15 
           16 import sys
           17 import time
           18 import numpy as np
           19 
           20 # try different netCDF modules
           21 try:
           22     from netCDF4 import Dataset as CDF
           23 except:
           24     print("netCDF4 is not installed!")
           25     sys.exit(1)
           26 
           27 # set up the grid:
           28 Lx = 1e6
           29 Ly = 1e6
           30 Mx = 51
           31 My = 71
           32 x = np.linspace(-Lx, Lx, Mx)
           33 y = np.linspace(-Ly, Ly, My)
           34 
           35 # create dummy fields
           36 [xx, yy] = np.meshgrid(x, y)  # if there were "ndgrid" in numpy we would use it
           37 acab = np.zeros((Mx, My))
           38 artm = np.zeros((Mx, My)) + 273.15 + 10.0  # 10 degrees Celsius
           39 topg = 1000.0 + 200.0 * (xx + yy) / max(Lx, Ly)  # change "1000.0" to "0.0" to test
           40 # flotation criterion, etc.
           41 thk = 3000.0 * (1.0 - 3.0 * (xx ** 2 + yy ** 2) / Lx ** 2)
           42 thk[thk < 0.0] = 0.0
           43 
           44 # Output filename
           45 ncfile = 'foo.nc'
           46 
           47 # Write the data:
           48 nc = CDF(ncfile, "w", format='NETCDF3_CLASSIC')  # for netCDF4 module
           49 
           50 # Create dimensions x and y
           51 nc.createDimension("x", size=Mx)
           52 nc.createDimension("y", size=My)
           53 
           54 x_var = nc.createVariable("x", 'f4', dimensions=("x",))
           55 x_var.units = "m"
           56 x_var.long_name = "easting"
           57 x_var.standard_name = "projection_x_coordinate"
           58 x_var[:] = x
           59 
           60 y_var = nc.createVariable("y", 'f4', dimensions=("y",))
           61 y_var.units = "m"
           62 y_var.long_name = "northing"
           63 y_var.standard_name = "projection_y_coordinate"
           64 y_var[:] = y
           65 
           66 fill_value = np.nan
           67 
           68 
           69 def def_var(nc, name, units, fillvalue):
           70     # dimension transpose is standard: "float thk(y, x)" in NetCDF file
           71     var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue)
           72     var.units = units
           73     return var
           74 
           75 
           76 bed_var = def_var(nc, "topg", "m", fill_value)
           77 bed_var.standard_name = "bedrock_altitude"
           78 bed_var[:] = topg
           79 
           80 thk_var = def_var(nc, "thk", "m", fill_value)
           81 thk_var.standard_name = "land_ice_thickness"
           82 thk_var[:] = thk
           83 
           84 acab_var = def_var(nc, "climatic_mass_balance", "m year-1", fill_value)
           85 acab_var.standard_name = "land_ice_surface_specific_mass_balance"
           86 acab_var[:] = acab
           87 
           88 artm_var = def_var(nc, "ice_surface_temp", "K", fill_value)
           89 artm_var[:] = artm
           90 
           91 # set global attributes
           92 nc.Conventions = "CF-1.4"
           93 historysep = ' '
           94 historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
           95 setattr(nc, 'history', historystr)
           96 
           97 nc.close()
           98 print('  PISM-bootable NetCDF file %s written' % ncfile)
           99 print('  for example, run:')
          100 print('    $ pismr -i foo.nc -bootstrap -Mx 41 -My 41 -Mz 21 -Lz 4000 -Mbz 5 -Lbz 500 -y 1')