-
Notifications
You must be signed in to change notification settings - Fork 1
/
gen_fig08.py
135 lines (121 loc) · 4.07 KB
/
gen_fig08.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
################################################################
# generate results for LoWo21 Figure 8
# raindrop evap for specific known/hypothesized planets
################################################################
import src.fall as fall
from src.planet import Planet
import src.drop_prop as drop_prop
import numpy as np
# set up different planets
# ??? indicates very guessed from a large range of possible values
# EARTH
X = np.zeros(5) # composition
X[2] = 0.8 # f_N2 [mol/mol]
X[3] = 0.2 # f_O2 [mol/mol]
T_surf = 290. # [K]
p_surf = 1.01325e5 # [Pa]
RH_surf = 0.75 # [ ]
R_p = 1. # [R_earth]
M_p = 1. # [M_earth]
Earth = Planet(R_p,T_surf,p_surf,X,'h2o',RH_surf,M_p)
# EARLY MARS
X = np.zeros(5) # composition
X[4] = 1. # f_CO2 [mol/mol]
T_surf = 290. # [K] ???
p_surf = 2e5 # [Pa]
RH_surf = 0.75 # [ ] ???
R_p = 0.532 # [R_earth]
M_p = 0.107 # [M_earth]
Mars = Planet(R_p,T_surf,p_surf,X,'h2o',RH_surf,M_p)
# gas giant compositions from Leconte et al. (2017)
# and cloud levels from Carlson et al. (1988)
# JUPITER
X = np.zeros(5) # composition
X[1] = 0.136 # f_He [mol/mol]
X[0] = 1 - X[1] # f_H2 [mol/mol]
T_LCL = 274. # [K]
p_LCL = 4.85e5 # [Pa]
RH_LCL = 1. # [ ]
R_p = 11.2089 # [R_earth]
M_p = 317.8 # [M_earth]
Jupiter = Planet(R_p,T_LCL,p_LCL,X,'h2o',RH_LCL,M_p)
# SATURN
X = np.zeros(5) # composition
X[1] = 0.12 # f_He [mol/mol]
X[0] = 1 - X[1] # f_H2 [mol/mol]
T_LCL = 284. # [K]
p_LCL = 10.4e5 # [Pa]
RH_LCL = 1. # [ ]
R_p = 9.449 # [R_earth]
M_p = 95.16 # [M_earth]
Saturn = Planet(R_p,T_LCL,p_LCL,X,'h2o',RH_LCL,M_p)
# K2-18b
X = np.zeros(5) # composition
X[1] = 0.10 # f_He [mol/mol] ???
X[0] = 1 - X[1] # f_H2 [mol/mol]
T_LCL = 275 # [K] ???
p_LCL = 1e4 # [Pa] ???
RH_LCL = 1. # [ ]
R_p = 2.610 # [R_earth] Benneke et al. (2019)
M_p = 8.63 # [M_earth] Cloutier et al. (2019) +- 1.35 M_earth
K2_18b = Planet(R_p,T_LCL,p_LCL,X,'h2o',RH_LCL,M_p)
planets = [Earth,Mars,Jupiter,Saturn,K2_18b]
planet_names = ['Earth','Mars','Jupiter','Saturn','K2-18b']
n_z = 25
H_LCL = np.zeros(len(planets)) # [m] atmospheric scale height @ LCL
r_min = np.zeros((n_z,len(planets),2)) # [m]
r_maxs = np.zeros(len(planets)) # [m]
z_maxs = np.zeros(len(planets)) # [m]
zH_maxs = np.zeros(len(planets)) # [H]
zs = np.zeros((n_z,len(planets))) # [m]
zHs = np.zeros((n_z,len(planets))) # [H]
for j,pl in enumerate(planets):
# maximum stable raindrop size
# from Rayleigh-Taylor instability, ℓ_RT = 0.5πa
r_maxs[j] = drop_prop.calc_r_max_RT(pl,np.pi/2.)
# calculate scale height at LCL to normalize against
H_LCL[j] = pl.calc_H(pl.z_LCL)
# set max Δz to consider
# distinguish between planets with surface & no surface
# for planets with no surface, end z set from where r_max
# finishes evaporating
if pl.z_LCL!=0:
z_maxs[j] = pl.z_LCL
else:
sol_rmax = fall.integrate_fall(pl,r_maxs[j])
z_maxs[j] = -fall.calc_last_z(sol_rmax)
# equally space z in meters in logspace
zs[:,j] = np.logspace(1,np.log10(z_maxs[j])-0.01,n_z)
# equally space z in H in logspace
zH_maxs[j] = z_maxs[j]/H_LCL[j]
zHs[:,j] = np.logspace(-3,np.log10(z_maxs[j]/H_LCL[j])-0.01,n_z)
# go through planets
for j,pl in enumerate(planets):
# go through each distance [m]
for i,z in enumerate(zs[:,j]):
# reset planet z to cloud base
pl.z2x4drdz(pl.z_LCL)
# calc end z from where cloud base is
z_end = pl.z_LCL - z
if z_maxs[j]<z:
r_min[i,j,0] = None
else:
# calc r_min for that z
r_min[i,j,0] = fall.calc_smallest_raindrop(pl,z_end=z_end,dr=5e-7)[0]
# go through each distance [H] & do same thing
for i,zH in enumerate(zHs[:,j]):
pl.z2x4drdz(pl.z_LCL)
z = zH*H_LCL[j]
z_end = pl.z_LCL - z
if z_maxs[j]<z:
r_min[i,j,1] = None
else:
r_min[i,j,1] = fall.calc_smallest_raindrop(pl,z_end=z_end,dr=5e-7)[0]
# save results
dir = 'output/fig08/'
np.save(dir+'r_mins_log',r_min)
np.save(dir+'zs_log',zs)
np.save(dir+'zHs_log',zHs)
np.save(dir+'r_maxs',r_maxs)
np.save(dir+'z_maxs',z_maxs)
np.save(dir+'zH_maxs',zH_maxs)