-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathround_forcing
executable file
·112 lines (91 loc) · 4.26 KB
/
round_forcing
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#!/usr/bin/env python
from __future__ import print_function
import flo_utils as fu
from argparse import ArgumentParser
import sys
import numpy as np
from math import log
def mix_forcing(year, options):
data = np.genfromtxt(options.timeseries)[::-1].transpose()
# time -- sea level -- isolation -- co2
# negative value in insol row/col is 115 ka fraction, positive value is 127 ka fraction
data[1] = 1 + data[1]/130. #
data[1] = data[1] * (data[1] < 1) + 1 * (data[1] >= 1) # cut values above 1
data[1] = data[1] * (data[1] > 0 ) # cut negative values
i000 = 480.394908
i127 = 551.656494
i115 = 442.319239
rl = 1./(i000-i115)
rh = 1./(i127-i000)
data[2] = (data[2] - i000)
data[2] = data[2] * rh * (data[2] > 0) + data[2] * rl * (data[2] < 0)
data[2] = data[2] * (data[2] < 1) + 1 * (data[2] >= 1) # cut values above 1
data[2] = data[2] * (data[2] > -1) - 1 * (data[2] <= -1) # cut values below -1
# co2 (now) = 284.7
co2now = 284.7
data[3] =( np.log(data[3]) - np.log(co2now/2.) ) / np.log((2))
data[3] = data[3] * (data[3] < 1) + 1 * (data[3] >= 1) # cut values above 1
data[3] = data[3] * (data[3] > 0 ) # cut negative values
# solar ueber low neutral high
# co2 level first index
# sea level second index
# 0 = pi 1 = lgm
neutral = [["pmu0121", "pmu0112"], ["pmu0123", "pmu0113"]]
low = [["pmu0111", "pmu0115"], ["pmu0119", "pmu0117"]]
high = [["pmu0110", "pmu0114"], ["pmu0118", "pmu0116"]]
n = np.array([ [ data[1] * (1 - abs(data[2])) * data[3] , (1 - data[1]) * (1 - abs(data[2])) * data[3] ],
[ data[1] * (1 - abs(data[2])) * (1 - data[3]) , (1 - data[1]) * (1 - abs(data[2])) * (1- data[3]) ]]).transpose((2,0,1))
l = (np.array([ [ data[1] * (abs(data[2])) * data[3] , (1 - data[1]) * (abs(data[2])) * data[3] ],
[ data[1] * (abs(data[2])) * (1 - data[3]) , (1 - data[1]) * (abs(data[2])) * (1- data[3]) ]]) * (data[2] < 0)).transpose((2,0,1))
h = (np.array([ [ data[1] * (abs(data[2])) * data[3] , (1 - data[1]) * (abs(data[2])) * data[3] ],
[ data[1] * (abs(data[2])) * (1 - data[3]) , (1 - data[1]) * (abs(data[2])) * (1- data[3]) ]]) * (data[2] > 0)).transpose((2,0,1))
if np.interp(year, data[0], data[2]) > 0 :
factors = h
names = high
else:
factors = l
names = low
pairs = []
for i in (0, 1):
for j in (0, 1):
nv = np.interp(year, data[0], n[:,i,j])
if nv != 0:
pairs.append([nv, neutral[i][j]])
fv = np.interp(year, data[0], factors[:,i,j])
if fv != 0:
pairs.append([fv, names[i][j]])
command = "cdo -a -splityear -settaxis,1000-01-01,00:30:00,1hour "
for p in pairs[:-1]:
command = command + " -add -mulc,%.5f %s/ebm_data_%s.grb "%(p[0], options.ebm_data_path, p[1])
command = command + " -mulc,%.5f %s/ebm_data_%s.grb "%(pairs[-1][0], options.ebm_data_path, pairs[-1][1])
command = command + " %s"%options.atm_file_name
fu.qo ( command )
command = "cdo -r -settaxis,1000-01-01,00:30:00,1month "
for p in pairs[:-1]:
command = command + " -add -mulc,%.5f %s/mpiom_%s_ncea_cdo_filled_xy.nc "%(p[0], options.ebm_data_path, p[1])
command = command + " -mulc,%.5f %s/mpiom_%s_ncea_cdo_filled_xy.nc "%(pairs[-1][0], options.ebm_data_path, pairs[-1][1])
command = command + " %s"%options.oce_file_name
fu.qo ( command )
# command = "ncatted -a units,time,o,c,'days since 0001-01-01 00:00:00' %s"%options.oce_file_name
# fu.qo ( command )
command = "ncks -A %s/1year.nc %s"%(fu.home, options.oce_file_name)
fu.qo ( command )
def parse_args():
'''Parses the command line arguments'''
parser = ArgumentParser()
parser.description = "Mix climate states based on time series data"
parser.add_argument("--year", type=int, required=True)
parser.add_argument("--ebm_data_path", required=True)
parser.add_argument("--atm_file_name", required=True)
parser.add_argument("--oce_file_name", required=True)
parser.add_argument("--timeseries", required=True, help="Timerseries file to use")
parser.add_argument("-v", "--verbose",
help='''Be verbose''', action="store_true")
options = parser.parse_args()
return options
def main():
'''Create forcing file for year specified'''
options = parse_args()
mix_forcing(options.year, options)
if __name__ == "__main__":
main()