forked from knaughten/roms_tools
-
Notifications
You must be signed in to change notification settings - Fork 0
/
mip_seasonal_cycle.py
213 lines (204 loc) · 8.91 KB
/
mip_seasonal_cycle.py
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
from numpy import *
# Calculate the average amplitude of the seasonal cycle in total basal mass
# loss over 2002-2016, in MetROMS, low-res FESOM, and high-res FESOM, for the
# entire continent as well as the Amery. Print results to the screen.
def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr):
# Year simulations start
year_start = 1992
# Years to avearge over
calc_start = 2002
calc_end = 2016
# Number of output steps per year in FESOM
peryear = 365/5
# Name of each ice shelf
names = ['Total Mass Loss', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf']
num_shelves = len(names)
i_total = 0
i_amery = 16
# Read ROMS logfile
roms_time = []
f = open(roms_logfile, 'r')
# Skip the first line (header for time array)
f.readline()
for line in f:
try:
roms_time.append(float(line))
except(ValueError):
# Reached the header for the next variable
break
# Set up array for mass loss values at each ice shelf
roms_massloss = empty([num_shelves, len(roms_time)])
index = 0
# Loop over ice shelves
while index < num_shelves:
t = 0
for line in f:
try:
roms_massloss[index, t] = float(line)
t += 1
except(ValueError):
# Reached the header for the next ice shelf
break
index +=1
f.close()
# Add start year to ROMS time array
roms_time = array(roms_time) + year_start
# Calculate amplitude for each year
#roms_amplitude_total = []
roms_min_total = []
roms_max_total = []
#roms_amplitude_amery = []
roms_min_amery = []
roms_max_amery = []
for year in range(calc_start, calc_end):
# Find the first index after the beginning of this year
t_start = nonzero(roms_time >= year)[0][0]
# Find the first index after the beginning of next year
tmp2 = nonzero(roms_time >= year+1)[0]
if len(tmp2)==0:
# No such index, but we might not have run out of data
# eg simulation that ends on 31 December 2016
t_end = len(roms_time)
else:
t_end = tmp2[0]
# Total Antarctica
# Get min and max between these bounds
curr_min = amin(roms_massloss[i_total,t_start:t_end])
curr_max = amax(roms_massloss[i_total,t_start:t_end])
roms_min_total.append(curr_min)
roms_max_total.append(curr_max)
# Save amplitude
#roms_amplitude_total.append(curr_max-curr_min)
# Amery
curr_min = amin(roms_massloss[i_amery,t_start:t_end])
curr_max = amax(roms_massloss[i_amery,t_start:t_end])
#roms_amplitude_amery.append(curr_max-curr_min)
roms_min_amery.append(curr_min)
roms_max_amery.append(curr_max)
# Calculate average amplitude
#roms_val_total = mean(array(roms_amplitude_total))
#roms_val_amery = mean(array(roms_amplitude_amery))
#print 'Average amplitude of seasonal cycle in total basal mass loss: '
#print 'ROMS (Total Antarctica): ' + str(roms_val_total)
#print 'ROMS (Amery): ' + str(roms_val_amery)
print 'Average min/max in total basal mass loss: '
print 'ROMS (Total Antarctica): ' + str(mean(array(roms_min_total))) + ' Gt/y min, ' + str(mean(array(roms_max_total))) + ' Gt/y max'
print 'ROMS (Amery): ' + str(mean(array(roms_min_amery))) + ' Gt/y min, ' + str(mean(array(roms_max_amery))) + ' Gt/y max'
# Read FESOM logfiles
# Low-res
f = open(fesom_logfile_lr, 'r')
# Skip the first line (header)
f.readline()
# Read total mass loss
tmp = []
num_time = 0
for line in f:
try:
tmp.append(float(line))
num_time += 1
except(ValueError):
# Reached the header for the first individual ice shelf
break
# Set up array for mass loss values at each ice shelf
fesom_massloss_lr = empty([num_shelves, num_time])
# Save the total values in the first index
fesom_massloss_lr[0,:] = array(tmp)
# Loop over ice shelves
index = 1
while index < num_shelves:
t = 0
for line in f:
try:
fesom_massloss_lr[index,t] = float(line)
t += 1
except(ValueError):
# Reached the header for the next ice shelf
break
index += 1
f.close()
# Calculate amplitude for each year
#fesom_amplitude_total_lr = []
fesom_min_total_lr = []
fesom_max_total_lr = []
#fesom_amplitude_amery_lr = []
fesom_min_amery_lr = []
fesom_max_amery_lr = []
for year in range(calc_start, calc_end):
t_start = peryear*(year-year_start)
t_end = peryear*(year+1-year_start)
# Total Antarctica
# Get min and max between these bounds
curr_min = amin(fesom_massloss_lr[i_total,t_start:t_end])
curr_max = amax(fesom_massloss_lr[i_total,t_start:t_end])
# Save amplitude
#fesom_amplitude_total_lr.append(curr_max-curr_min)
fesom_min_total_lr.append(curr_min)
fesom_max_total_lr.append(curr_max)
# Amery
curr_min = amin(fesom_massloss_lr[i_amery,t_start:t_end])
curr_max = amax(fesom_massloss_lr[i_amery,t_start:t_end])
#fesom_amplitude_amery_lr.append(curr_max-curr_min)
fesom_min_amery_lr.append(curr_min)
fesom_max_amery_lr.append(curr_max)
# Calculate average amplitude
#fesom_val_total_lr = mean(array(fesom_amplitude_total_lr))
#fesom_val_amery_lr = mean(array(fesom_amplitude_amery_lr))
#print 'FESOM low-res (Total Antarctica): ' + str(fesom_val_total_lr)
#print 'FESOM low-res (Amery): ' + str(fesom_val_amery_lr)
print 'FESOM low-res (Total Antarctica): ' + str(mean(array(fesom_min_total_lr))) + ' Gt/y min, ' + str(mean(array(fesom_max_total_lr))) + ' Gt/y max'
print 'FESOM low-res (Amery): ' + str(mean(array(fesom_min_amery_lr))) + ' Gt/y min, ' + str(mean(array(fesom_max_amery_lr))) + ' Gt/y max'
# High-res
f = open(fesom_logfile_hr, 'r')
f.readline()
tmp = []
num_time = 0
for line in f:
try:
tmp.append(float(line))
num_time += 1
except(ValueError):
break
fesom_massloss_hr = empty([num_shelves, num_time])
fesom_massloss_hr[0,:] = array(tmp)
index = 1
while index < num_shelves:
t = 0
for line in f:
try:
fesom_massloss_hr[index,t] = float(line)
t += 1
except(ValueError):
break
index += 1
f.close()
#fesom_amplitude_total_hr = []
fesom_min_total_hr = []
fesom_max_total_hr = []
#fesom_amplitude_amery_hr = []
fesom_min_amery_hr = []
fesom_max_amery_hr = []
for year in range(calc_start, calc_end):
t_start = peryear*(year-year_start)
t_end = peryear*(year+1-year_start)
curr_min = amin(fesom_massloss_hr[i_total,t_start:t_end])
curr_max = amax(fesom_massloss_hr[i_total,t_start:t_end])
fesom_min_total_hr.append(curr_min)
fesom_max_total_hr.append(curr_max)
#fesom_amplitude_total_hr.append(curr_max-curr_min)
curr_min = amin(fesom_massloss_hr[i_amery,t_start:t_end])
curr_max = amax(fesom_massloss_hr[i_amery,t_start:t_end])
#fesom_amplitude_amery_hr.append(curr_max-curr_min)
fesom_min_amery_hr.append(curr_min)
fesom_max_amery_hr.append(curr_max)
#fesom_val_total_hr = mean(array(fesom_amplitude_total_hr))
#fesom_val_amery_hr = mean(array(fesom_amplitude_amery_hr))
#print 'FESOM high-res (Total Antarctica): ' + str(fesom_val_total_hr)
#print 'FESOM high-res (Amery): ' + str(fesom_val_amery_hr)
print 'FESOM high-res (Total Antarctica): ' + str(mean(array(fesom_min_total_hr))) + ' Gt/y min, ' + str(mean(array(fesom_max_total_hr))) + ' Gt/y max'
print 'FESOM high-res (Amery): ' + str(mean(array(fesom_min_amery_hr))) + ' Gt/y min, ' + str(mean(array(fesom_max_amery_hr))) + ' Gt/y max'
# Command-line interface
if __name__ == "__main__":
roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ")
fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ")
fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ")
mip_seasonal_cycle(roms_logfile, fesom_logfile_lr, fesom_logfile_hr)