tcreate_warming_climate.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_warming_climate.py (2691B)
       ---
            1 #!/usr/bin/env python3
            2 # Copyright (C) 2017 Andy Aschwanden
            3 
            4 import numpy as np
            5 import time
            6 from netCDF4 import Dataset as NC
            7 from argparse import ArgumentParser
            8 
            9 
           10 # Set up the option parser
           11 parser = ArgumentParser()
           12 parser.description = "Create climate forcing for a warming climate"
           13 parser.add_argument("FILE", nargs='*')
           14 parser.add_argument("-T_max", dest="T_max", type=float,
           15                     help="Maximum temperature", default=1)
           16 parser.add_argument("-t_max", dest="t_max", type=float,
           17                     help="lower time bound for maximum temperature", default=100)
           18 parser.add_argument("-amplitude", dest="amplitude", type=float,
           19                     help="Amplitde of seasonal cycle.", default=12)
           20 
           21 
           22 options = parser.parse_args()
           23 args = options.FILE
           24 start = 0
           25 end = 1000
           26 step = 1./12.
           27 amplitude = options.amplitude
           28 t_max = options.t_max
           29 T_max = options.T_max
           30 bnds_interval_since_refdate = np.linspace(start, end, end * 12 + 1)
           31 time_interval_since_refdate = (bnds_interval_since_refdate[0:-1] +
           32                                np.diff(bnds_interval_since_refdate) / 2)
           33 
           34 infile = args[0]
           35 
           36 nc = NC(infile, 'w')
           37 
           38 
           39 def def_var(nc, name, units):
           40     var = nc.createVariable(name, 'f', dimensions=('time'))
           41     var.units = units
           42     return var
           43 
           44 
           45 # create a new dimension for bounds only if it does not yet exist
           46 time_dim = "time"
           47 if time_dim not in list(nc.dimensions.keys()):
           48     nc.createDimension(time_dim)
           49 
           50 # create a new dimension for bounds only if it does not yet exist
           51 bnds_dim = "nb2"
           52 if bnds_dim not in list(nc.dimensions.keys()):
           53     nc.createDimension(bnds_dim, 2)
           54 
           55 # variable names consistent with PISM
           56 time_var_name = "time"
           57 bnds_var_name = "time_bnds"
           58 
           59 # create time variable
           60 time_var = nc.createVariable(time_var_name, 'd', dimensions=(time_dim))
           61 time_var[:] = time_interval_since_refdate
           62 time_var.bounds = bnds_var_name
           63 time_var.units = 'years since 1-1-1'
           64 time_var.calendar = '365_day'
           65 time_var.standard_name = time_var_name
           66 time_var.axis = "T"
           67 
           68 # create time bounds variable
           69 time_bnds_var = nc.createVariable(bnds_var_name, 'd', dimensions=(time_dim, bnds_dim))
           70 time_bnds_var[:, 0] = bnds_interval_since_refdate[0:-1]
           71 time_bnds_var[:, 1] = bnds_interval_since_refdate[1::]
           72 
           73 var = 'delta_T'
           74 dT_var = def_var(nc, var, "K")
           75 T_0 = 0.
           76 
           77 temp = np.zeros_like(time_interval_since_refdate) + T_max
           78 temp[0:int(t_max/step)] = np.linspace(T_0, T_max, t_max / step)
           79 temp[:] += -np.cos(time_interval_since_refdate * 2 * np.pi) * amplitude
           80 dT_var[:] = temp
           81 
           82 # writing global attributes
           83 script_command = ' '.join([time.ctime(), ':', __file__.split('/')[-1],
           84                            ' '.join([str(x) for x in args])])
           85 nc.history = script_command
           86 nc.Conventions = "CF 1.6"
           87 nc.close()