-
Notifications
You must be signed in to change notification settings - Fork 2
/
make_table.py
115 lines (77 loc) · 2.75 KB
/
make_table.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
import numpy as np
import netCDF4
def polyfit_omega(n=6):
'fit an order-n polynomial to the maximum growth rate as a function of delta, the slope parameter.'
delta, mu = np.mgrid[-1.2:2.2:1001j, 0:4.2:1001j]
tmu = np.tanh(mu)
omega = np.sqrt( (1.0+delta)*(mu - tmu)/tmu
-0.25*(delta/tmu + mu)**2 ).real
omega2 = (1.0+delta)*(mu - tmu)/tmu - 0.25*(delta/tmu + mu)**2
omega = np.ma.masked_where(np.isnan(omega), omega)
omega_max = omega.max(axis=1)
idx = np.where(~omega_max.mask)
omega_max = omega_max[idx]
delta = delta[:, 0][idx]
p = np.polyfit(delta, omega_max, n)
return p
omega_poly = polyfit_omega()
deltas = []
Ris = []
Ss = []
M2s = []
fs = []
N2 = 1e-4
alpha = 1e-3
for S in [0.1, 0.2, 0.3, 0.5]:
for Ri in [1, 2, 3, 5, 10]:
Ss.append(S)
Ris.append(Ri)
delta = S*np.sqrt(Ri)
deltas.append(delta)
f = np.sqrt(N2)*alpha/S
fs.append(f)
M2 = np.sqrt(N2*f**2/Ri)
M2s.append(M2)
for delta in [0.1, 0.2, 0.3, 0.5]:
for Ri in [1, 2, 3, 5, 10]:
deltas.append(delta)
Ris.append(Ri)
S = delta/np.sqrt(Ri)
Ss.append(S)
f = np.sqrt(N2)*alpha/S
fs.append(f)
M2 = np.sqrt(N2*f**2/Ri)
M2s.append(M2)
tab = np.vstack((M2s, fs, Ris, deltas, Ss)).T
_, idx = np.unique(tab.sum(axis=-1), return_index=True)
idx.sort()
tab = tab[idx]
# shelfstrat_M2_3.33e-07_N2_1.00e-04_f_5.77e-05
times = []
omegas = []
Tis = []
Tfs = []
for case in tab:
# print('%5.2e, %5.2e, %5.2f, %5.2f, %5.2f' % tuple(case))
filename = 'simulations/shelfstrat_M2_%5.2e_N2_1.00e-04_f_%5.2e/shelfstrat_his.nc' % (case[0], case[1])
print filename
nc = netCDF4.Dataset(filename)
times.append(nc.variables['ocean_time'][-1]/86400)
S = case[4]
delta = case[3]
Ri = case[2]
f = case[1]
omega = np.polyval(omega_poly, delta) # non-dim
omega_dim = 86400.0 * omega * f / np.sqrt(Ri) # rad/days
omegas.append(omega_dim)
Tis.append(np.sqrt(S) / omega_dim)
###################################
diafilename = 'simulations/shelfstrat_M2_%5.2e_N2_1.00e-04_f_%5.2e/shelfstrat_dia.nc' % (case[0], case[1])
nd = netCDF4.Dataset(diafilename)
ubar_bstr = nd.variables['ubar_bstr'][:, 1:51, :]
stress_mag = np.sqrt(ubar_bstr**2)
stress = np.sqrt((stress_mag**2).mean(axis=-1).mean(axis=-1))
ubar = Ri**(-1./2.) * np.sqrt(1e-4) * 50.0 / 2.0
Tfs.append(ubar/stress[1]/86400.0)
tab = np.vstack((tab.T, times, omegas, Tis, Tfs)).T
np.savetxt('table.dat', tab, fmt='%5.2e & %5.2e & %5.2f & %5.2f & %5.2f & %6.1f & %5.3f & %5.3f & %5.3f \\\\')