-
Notifications
You must be signed in to change notification settings - Fork 11
/
functions.py
182 lines (153 loc) · 6.21 KB
/
functions.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
import math
import numpy as np
from libcloudphxx import common as cm
def mole_frac_to_mix_ratio(mole_frac_g, p, M, T, rhod):
"""
convert mole fraction [1] to mixing ratio [kg/kg dry air]
"""
return mole_frac_g * p * M / cm.R / T / rhod
def mix_ratio_to_mole_frac(mix_r, p, M, T, rhod):
"""
convert mixing ratio [kg/kg dry air] to mole fraction [1]
"""
return mix_r * cm.R * T * rhod / p / M
def rh_to_rv(RH, T, p):
"""
convert relative humidity [%] to water vapor mixing ratio [kg/kg]
"""
return cm.eps * RH * cm.p_vs(T) / (p - RH * cm.p_vs(T))
def rhod_calc(T, p, rv):
"""
calculate dry air density
"""
th_0 = T * (cm.p_1000 / p)**(cm.R_d / cm.c_pd)
return cm.rhod(p, th_0, rv)
def rho_calc(T, p, rv):
"""
calculate total air density (not dry air density)
"""
rho = p / T / (rv/(1.+rv) * cm.R_v + 1./(1.+rv) * cm.R_d)
return rho
def henry_teor(chem, p, T, vol, mixr_g, rhod, conc_H):
"""
calculate theoretical mass of chemical species dissolved into cloud droplets - Henry law
(per kg of dry air)
"""
# read in the molar mass of chemical species dissolved in water ...
if chem in ["SO2", "CO2", "NH3"]:
M_aq = getattr(cm, "M_"+chem+"_H2O")
else:
M_aq = getattr(cm, "M_"+chem)
# ... and the molar mass of chem species in gas phase
M_g = getattr(cm, "M_"+chem)
# calculate the correction due to dissociation
if chem in ["O3", "H2O2"]:
pH_corr = 1 # do nothing for those that don't dissociate
elif chem == "SO2":
K1 = getattr(cm, "K_SO2") * np.exp(getattr(cm, "dKR_SO2") * (1./T - 1./298))
K2 = getattr(cm, "K_HSO3")* np.exp(getattr(cm, "dKR_HSO3") * (1./T - 1./298))
pH_corr = 1. + K1 / conc_H + K1 * K2 / conc_H / conc_H
elif chem == "CO2":
K1 = getattr(cm, "K_CO2") * np.exp(getattr(cm, "dKR_CO2") * (1./T - 1./298))
K2 = getattr(cm, "K_HCO3")* np.exp(getattr(cm, "dKR_HCO3") * (1./T - 1./298))
pH_corr = 1. + K1 / conc_H + K1 * K2 / conc_H / conc_H
elif chem == "HNO3":
K = getattr(cm, "K_HNO3") * np.exp(getattr(cm, "dKR_"+chem) * (1./T - 1./298))
pH_corr = 1. + K / conc_H
elif chem == "NH3":
K = getattr(cm, "K_NH3") * np.exp(getattr(cm, "dKR_"+chem) * (1./T - 1./298))
pH_corr = 1 + K / getattr(cm, "K_H2O") * conc_H
else:
raise Exception("chem should be O3, H2O2, SO3, CO2, HNO3, or NH3")
# correction to Henry constant due to temperature and pH
H = getattr(cm, "H_" +chem)
dHR = getattr(cm, "dHR_"+chem)
henry_T = H * np.exp(dHR * (1./T - 1./298)) * pH_corr
# dissolved = partial prsessure * Henry_const * molar mass * drop volume
partial_prs = mix_ratio_to_mole_frac(mixr_g, p, M_g, T, rhod) * p
return partial_prs * henry_T * M_aq * vol
def dissoc_teor(chem, T):
"""
calculate theoretical dissociation constants (as a function of temperature)
"""
# correction to dissoc constant due to temperature
K = getattr(cm, "K_" +chem)
dKR = getattr(cm, "dKR_"+chem)
return K * np.exp(dKR * (1./T - 1./298))
def log10_size_of_lnr(n_tot, mean_r, lnr, gstdev):
"""
log-normal size distribution (defined as a function of log_10(radius))
"""
return n_tot / math.sqrt(2 * math.pi) / math.log(gstdev, 10)\
* math.exp(-1. * math.pow((lnr - math.log(mean_r, 10)), 2) / 2. / math.pow(math.log(gstdev,10),2))
def diag_n_OH(V, conc_H):
"""
calculate the number of OH moles
"""
return cm.K_H2O * V / conc_H
def diag_n_NH3_H2O(m_N3, T, conc_H):
"""
calculate the number of NH3*H2O moles
"""
return m_N3 / cm.M_NH3_H2O / (1. + dissoc_teor("NH3", T) * conc_H / cm.K_H2O)
def diag_n_NH4(m_N3, T, conc_H):
"""
calculate the number of NH4+ moles
"""
return m_N3 / cm.M_NH3_H2O * conc_H * dissoc_teor("NH3", T) / cm.K_H2O / (1. + dissoc_teor("NH3", T) * conc_H / cm.K_H2O)
def diag_n_HNO3(m_N5, T, conc_H):
"""
calculate the number of HNO3 moles
"""
return m_N5 / cm.M_HNO3 / (dissoc_teor("HNO3", T) / conc_H + 1)
def diag_n_NO3(m_N5, T, conc_H):
"""
calculate the number of NO3- moles
"""
return dissoc_teor("HNO3", T) / conc_H * m_N5 / cm.M_HNO3 / (dissoc_teor("HNO3", T) / conc_H + 1)
def diag_n_CO2_H2O(m_C4, T, conc_H):
"""
calculate the number of CO2*H2O moles
"""
return m_C4 / cm.M_CO2_H2O \
/ (1 + dissoc_teor("CO2", T) / conc_H + dissoc_teor("CO2", T) * dissoc_teor("HCO3", T) / conc_H / conc_H)
def diag_n_HCO3(m_C4, T, conc_H):
"""
calculate the number of HCO3- moles
"""
return m_C4 / cm.M_CO2_H2O * dissoc_teor("CO2", T) / conc_H \
/ (1 + dissoc_teor("CO2", T) / conc_H + dissoc_teor("CO2", T) * dissoc_teor("HCO3", T) / conc_H / conc_H)
def diag_n_CO3(m_C4, T, conc_H):
"""
calculate the number of CO3-- moles
"""
return m_C4 / cm.M_CO2_H2O * dissoc_teor("CO2", T) * dissoc_teor("HCO3", T) / conc_H / conc_H \
/ (1 + dissoc_teor("CO2", T) / conc_H + dissoc_teor("CO2", T) * dissoc_teor("HCO3", T) / conc_H / conc_H)
def diag_n_SO2_H2O(m_S4, T, conc_H):
"""
calculate the number of SO2*H2O moles
"""
return m_S4 / cm.M_SO2_H2O \
/ (1 + dissoc_teor("SO2", T) / conc_H + dissoc_teor("SO2", T) * dissoc_teor("HSO3", T) / conc_H / conc_H)
def diag_n_HSO3(m_S4, T, conc_H):
"""
calculate the number of HSO3- moles
"""
return m_S4 / cm.M_SO2_H2O * dissoc_teor("SO2", T) / conc_H \
/ (1 + dissoc_teor("SO2", T) / conc_H + dissoc_teor("SO2", T) * dissoc_teor("HSO3", T) / conc_H / conc_H)
def diag_n_SO3(m_S4, T, conc_H):
"""
calculate the number of SO3-- moles
"""
return m_S4 / cm.M_SO2_H2O * dissoc_teor("SO2", T) * dissoc_teor("HSO3", T) / conc_H / conc_H \
/ (1 + dissoc_teor("SO2", T) / conc_H + dissoc_teor("SO2", T) * dissoc_teor("HSO3", T) / conc_H / conc_H)
def diag_n_HSO4(m_S6, T, conc_H):
"""
calculate the number of HSO4- moles
"""
return m_S6 / cm.M_H2SO4 * conc_H / (conc_H + dissoc_teor("HSO4", T))
def diag_n_SO4(m_S6, T, conc_H):
"""
calculate the number of SO4-- moles
"""
return m_S6 / cm.M_H2SO4 * dissoc_teor("HSO4", T) / (conc_H + dissoc_teor("HSO4", T))