tcreate_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
       ---
       tcreate_inputs.py (1943B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 import netCDF4
            4 import numpy as np
            5 
            6 M = 3                           # grid size
            7 L = 1e5                         # domain size
            8 b0 = 0.0                        # bed elevation
            9 H0 = 200.0                      # ice thickness
           10 T_mean_annual = 268.15        # mean annual temperature, Kelvin
           11 T_amplitude = 6             # surface temperature aplitude, Kelvin
           12 summer_peak_day = 365/2
           13 seconds_per_year = 365 * 86400
           14 M0 = 0.0                        # mass balance
           15 n_records = 53
           16 
           17 # use a square grid
           18 Mx = M
           19 My = M
           20 Lx = L
           21 Ly = L
           22 
           23 
           24 def T_surface(time, mean, amplitude, summer_peak_day):
           25     "Surface temperature (cosine yearly cycle)"
           26     day_length = 86400
           27     summer_peak = summer_peak_day * day_length
           28     year_length = float(365 * day_length)
           29 
           30     t = np.mod(time - summer_peak, year_length) / year_length
           31 
           32     return mean + amplitude * np.cos(2 * np.pi * t)
           33 
           34 
           35 f = netCDF4.Dataset("input.nc", "w")
           36 
           37 f.createDimension("x", Mx)
           38 x = f.createVariable("x", np.float64, ("x",))
           39 x.units = "meters"
           40 x[:] = np.linspace(-Lx, Lx, Mx)
           41 
           42 f.createDimension("y", My)
           43 y = f.createVariable("y", np.float64, ("y",))
           44 y.units = "meters"
           45 y[:] = np.linspace(-Ly, Ly, My)
           46 
           47 f.createDimension("time", n_records)
           48 
           49 time = np.linspace(0, 1, n_records) * seconds_per_year
           50 
           51 t = f.createVariable("time", np.float64, ("time",))
           52 t.units = "seconds since 1-1-1"
           53 t[:] = time
           54 
           55 H = f.createVariable("thk", np.float64, ("y", "x"))
           56 H.units = "meters"
           57 H.long_name = "ice thickness"
           58 H[:] = H0
           59 
           60 b = f.createVariable("topg", np.float64, ("y", "x"))
           61 b.units = "meters"
           62 b.long_name = "bed elevation"
           63 b[:] = b0
           64 
           65 T_s = f.createVariable("ice_surface_temp", np.float64, ("y", "x"))
           66 T_s.units = "Kelvin"
           67 T_s[:] = T_mean_annual
           68 
           69 M = f.createVariable("climatic_mass_balance", np.float64, ("y", "x"))
           70 M.units = "kg m-2 s-1"
           71 M[:] = M0
           72 
           73 dT = f.createVariable("delta_T", np.float64, ("time",))
           74 dT.units = "Kelvin"
           75 dT[:] = T_surface(time, 0.0, T_amplitude, summer_peak_day)
           76 
           77 f.close()