-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathtimeseries_3D.py
186 lines (163 loc) · 6.5 KB
/
timeseries_3D.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
from netCDF4 import Dataset
from numpy import *
from matplotlib.pyplot import *
from os.path import *
from cartesian_grid_3d import *
from rotate_vector_roms import *
# Calculate and plot timeseries of total ocean heat content, average salinity,
# and total kinetic energy during a ROMS simulation.
# Takes 32 GB of memory for Kaitlin's circumpolar quarter-degree grid to 30S.
# Input:
# grid_path = path to ROMS grid file
# file_path = path to ROMS history/averages file
# log_path = path to log file (if it exists, previously calculated values will
# be read from it; regardless, it will be overwritten with all
# calculated values following computation)
def timeseries_3D (grid_path, file_path, log_path):
# Grid parameters
theta_s = 7.0
theta_b = 2.0
hc = 250
N = 31
rho0 = 1000.0 # Reference density (kg/m^3)
Cp = 3974 # Specific heat of polar seawater (J/K/kg)
C2K = 273.15 # Celsius to Kelvin conversion
time = []
# ohc = []
avgsalt = []
# tke = []
# Check if the log file exists
if exists(log_path):
print 'Reading previously calculated values'
f = open(log_path, 'r')
# Skip first line (header for time array)
f.readline()
for line in f:
try:
time.append(float(line))
except(ValueError):
# Reached the header for the next variable
break
# for line in f:
# try:
# ohc.append(float(line))
# except(ValueError):
# break
for line in f:
try:
avgsalt.append(float(line))
except(ValueError):
break
# for line in f:
# tke.append(float(line))
f.close()
print 'Analysing grid'
id = Dataset(grid_path, 'r')
h = id.variables['h'][:-15,1:-1]
zice = id.variables['zice'][:-15,1:-1]
lon = id.variables['lon_rho'][:-15,1:-1]
lat = id.variables['lat_rho'][:-15,1:-1]
mask = id.variables['mask_rho'][:-15,1:-1]
# Keep the overlapping periodic boundary on "angle" for now
angle = id.variables['angle'][:-15,:]
id.close()
id = Dataset(file_path, 'r')
# Read time values and convert from seconds to years
new_time = id.variables['ocean_time'][:]/(60*60*24*365.25)
num_time = size(new_time)
# Concatenate with time values from log file
for t in range(num_time):
time.append(new_time[t])
# Process 10 time indices at a time so we don't use too much memory
start_t = 0
while True:
end_t = min(start_t+10, num_time)
print 'Processing time indices ' + str(start_t+1) + ' to ' + str(end_t)
num_time_curr = end_t-start_t
print 'Calculating time-dependent dV'
# Read time-dependent sea surface height
zeta = id.variables['zeta'][start_t:end_t,:-15,1:-1]
# Calculate time-dependent dz
dz = ma.empty([num_time_curr, N, size(lon,0), size(lon,1)])
for t in range(num_time_curr):
# dx and dy will be recomputed unnecessarily each timestep
# but that's ok
dx, dy, dz_tmp, z = cartesian_grid_3d(lon, lat, h, zice, theta_s, theta_b, hc, N, zeta[t,:,:])
dz[t,:,:,:] = dz_tmp
# Calculate time-dependent dV and mask with land mask
# Here mask, dx, dy are all copied into arrays of dimension
# time x depth x lat x lon
dV = ma.masked_where(tile(mask, (num_time_curr,N,1,1))==0, tile(dx, (num_time_curr,1,1,1))*tile(dy, (num_time_curr,1,1,1))*dz)
print 'Reading data'
# temp = id.variables['temp'][start_t:end_t,:,:-15,1:-1]
salt = id.variables['salt'][start_t:end_t,:,:-15,1:-1]
rho = id.variables['rho'][start_t:end_t,:,:-15,1:-1] + rho0
# Keep overlapping periodic boundary for u and v
# u_xy = id.variables['u'][start_t:end_t,:,:-15,:]
# v_xy = id.variables['v'][start_t:end_t,:,:-15,:]
# print 'Interpolating velocities onto rho-grid'
# # We are actually rotating them at the same time as interpolating
# # which is a bit of unnecessary work (sum of squares won't change with
# # rotation) but not much extra work, and it's conveneint
# u = ma.empty(shape(temp))
# v = ma.empty(shape(temp))
# for t in range(num_time_curr):
# for k in range(N):
# u_tmp, v_tmp = rotate_vector_roms(u_xy[t,k,:,:], v_xy[t,k,:,:], angle)
# u[t,k,:,:] = u_tmp[:,1:-1]
# v[t,k,:,:] = v_tmp[:,1:-1]
print 'Building timeseries'
for t in range(num_time_curr):
# Integrate temp*rho*Cp*dV to get OHC
# ohc.append(sum((temp[t,:,:,:]+C2K)*rho[t,:,:,:]*Cp*dV[t,:,:,:]))
# Average salinity (weighted with rho*dV)
avgsalt.append(sum(salt[t,:,:,:]*rho[t,:,:,:]*dV[t,:,:,:])/sum(rho[t,:,:,:]*dV[t,:,:,:]))
# Integrate 0.5*rho*speed^2*dV to get TKE
# tke.append(sum(0.5*rho[t,:,:,:]*(u[t,:,:,:]**2 + v[t,:,:,:]**2)*dV[t,:,:,:]))
# Get ready for next 10 time indices
if end_t == num_time:
break
start_t = end_t
id.close()
# print 'Plotting ocean heat content'
# clf()
# plot(time, ohc)
# xlabel('Years')
# ylabel('Southern Ocean Heat Content (J)')
# grid(True)
# savefig('ohc.png')
print 'Plotting average salinity'
clf()
plot(time, avgsalt)
xlabel('Years')
ylabel('Southern Ocean Average Salinity (psu)')
grid(True)
savefig('avgsalt.png')
# print 'Plotting total kinetic energy'
# clf()
# plot(time, tke)
# xlabel('Years')
# ylabel('Southern Ocean Total Kinetic Energy (J)')
# grid(True)
# savefig('tke.png')
print 'Saving results to log file'
f = open(log_path, 'w')
f.write('Time (years):\n')
for elm in time:
f.write(str(elm) + '\n')
# f.write('Southern Ocean Heat Content (J):\n')
# for elm in ohc:
# f.write(str(elm) + '\n')
f.write('Southern Ocean Average Salinity (psu):\n')
for elm in avgsalt:
f.write(str(elm) + '\n')
# f.write('Southern Ocean Total Kinetic Energy (J):\n')
# for elm in tke:
# f.write(str(elm) + '\n')
f.close()
# Command-line interface
if __name__ == "__main__":
grid_path = raw_input("Path to ROMS grid file: ")
file_path = raw_input("Path to ROMS history/averages file: ")
log_path = raw_input("Path to logfile to save values and/or read previously calculated values: ")
timeseries_3D(grid_path, file_path, log_path)