tgenerate_inputs.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
       ---
       tgenerate_inputs.py (3455B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 import sys
            4 import time
            5 import netCDF4
            6 import numpy as np
            7 
            8 
            9 def generate_input(filename):
           10     "Generate a slab-on-a-flat-bed input file for PISM."
           11     # Initialize a square 3x3 grid
           12     Mx = 3
           13     My = Mx
           14     Lx = 50e3   # m
           15     Ly = Lx
           16 
           17     x = np.linspace(-Lx, Lx, Mx)
           18     y = np.linspace(-Ly, Ly, My)
           19     dx = x[1]-x[0]
           20     dy = dx
           21 
           22     topg = np.zeros((My, Mx))
           23     thk = np.zeros((My, Mx)) + 1000.0
           24     artm = np.zeros((My, Mx)) + 273.15 - 30  # -30 degrees Celsius
           25     hflx = np.zeros((My, Mx)) + 0.042       # 42 mW/m2
           26     smb = np.zeros((My, Mx))                # zero
           27 
           28     # Write the data:
           29     nc = netCDF4.Dataset(filename, "w")
           30 
           31     # Create dimensions x and y
           32     nc.createDimension("x", size=Mx)
           33     nc.createDimension("y", size=My)
           34 
           35     x_var = nc.createVariable("x", 'f4', dimensions=("x",))
           36     x_var.units = "m"
           37     x_var.standard_name = "projection_x_coordinate"
           38     x_var[:] = x
           39 
           40     y_var = nc.createVariable("y", 'f4', dimensions=("y",))
           41     y_var.units = "m"
           42     y_var.standard_name = "projection_y_coordinate"
           43     y_var[:] = y
           44 
           45     def def_var(nc, name, units, fillvalue):
           46         var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue)
           47         var.units = units
           48         return var
           49 
           50     bed_var = def_var(nc, "topg", "m", fill_value)
           51     bed_var.standard_name = "bedrock_altitude"
           52     bed_var[:] = topg
           53 
           54     thk_var = def_var(nc, "thk", "m", fill_value)
           55     thk_var.standard_name = "land_ice_thickness"
           56     thk_var[:] = thk
           57 
           58     hflx_var = def_var(nc, "bheatflux", "W m-2", fill_value)
           59     hflx_var.standard_name = "land_ice_basal_heat_flux"
           60     hflx_var[:] = hflx
           61 
           62     smb_var = def_var(nc, "climatic_mass_balance", "kg m-2 s-1", fill_value)
           63     smb_var.standard_name = "land_ice_surface_specific_mass_balance"
           64     smb_var[:] = smb
           65 
           66     artm_var = def_var(nc, "ice_surface_temp", "K", fill_value)
           67     artm_var[:] = artm
           68 
           69     # set global attributes
           70     nc.Conventions = 'CF-1.4'
           71     setattr(nc, 'history', historystr)
           72 
           73     nc.close()
           74 
           75 
           76 def generate_forcing(filename):
           77     "Generate surface temperature forcing."
           78     ts_yr = 0
           79     te_yr = 300000
           80     dt_yr = 10
           81 
           82     Mt = (te_yr - ts_yr) / dt_yr + 1
           83     delta_T = np.zeros(Mt)
           84     time = np.linspace(ts_yr, te_yr, Mt)
           85 
           86     delta_T[np.where((time > 100000))] = 25.0
           87     delta_T[np.where((time >= 150000))] = 0.0
           88 
           89     # Write the data:
           90     nc = netCDF4.Dataset(dTfile, "w")
           91 
           92     # Create dimensions x and y
           93     nc.createDimension("time", size=Mt)
           94 
           95     def def_var(nc, name, units, fillvalue):
           96         var = nc.createVariable(name, 'f', dimensions=("time"), fill_value=fillvalue)
           97         var.units = units
           98         return var
           99 
          100     t_var = def_var(nc, "time", "years since since 1-1-1", fill_value)
          101     t_var[:] = time
          102 
          103     dT_var = def_var(nc, "delta_T", "K", fill_value)
          104     dT_var.standard_name = "land_ice_temperature_at_firn_base"
          105     dT_var.long_name = "Temperature (variation from present)"
          106     dT_var[:] = delta_T
          107 
          108     # set global attributes
          109     nc.Conventions = 'CF-1.4'
          110     setattr(nc, 'history', historystr)
          111 
          112     nc.close()
          113 
          114 
          115 if __name__ == "__main__":
          116     fill_value = np.nan
          117 
          118     prefix = 'box'
          119     ncfile = prefix + '.nc'
          120     dTfile = prefix + '_dT.nc'
          121 
          122     historysep = ' '
          123     historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
          124 
          125     generate_input(ncfile)
          126     print('Wrote PISM input file %s.' % ncfile)
          127 
          128     generate_forcing(dTfile)
          129     print('Wrote PISM surface temperature forcing file %s.' % dTfile)