Skip to content

Commit

Permalink
Merge pull request #27 from MDS-007-52/main
Browse files Browse the repository at this point in the history 
Updated H2O-H2O, H2O-N2 and N2-N2 continua by Microwave spectorclopy laboratory, IAP RAS
  • Loading branch information
@slarosa
slarosa authored Jul 24, 2024
2 parents 3ca3424 + 8a9a8d6 commit fe96776
Show file tree
Hide file tree
Showing 4 changed files with 143 additions and 21 deletions.
Binary file modified pyrtlib/_lineshape/h2o_lineshape.nc
View file
Binary file not shown.
8 changes: 4 additions & 4 deletions pyrtlib/_lineshape/h2oll.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,14 +78,14 @@
xh = mtx[:, 9]
shs = mtx[:, 10] / 1000.0
xhs = mtx[:, 11]
if H2OAbsModel.model in ['R19', 'R19SD', 'R20', 'R20SD', 'R21SD', 'R22SD', 'R23', 'R23SD', 'R24']:
if H2OAbsModel.model in ['R19', 'R19SD', 'R20', 'R20SD', 'R21SD', 'R22SD', 'R23', 'R23SD', 'R24', 'MWL24']:
aair = mtx[:, 12]
aself = mtx[:, 13]
if H2OAbsModel.model in ['R19SD', 'R20SD', 'R21SD', 'R22SD', 'R23', 'R23SD', 'R24']:
if H2OAbsModel.model in ['R19SD', 'R20SD', 'R21SD', 'R22SD', 'R23', 'R23SD', 'R24', 'MWL24']:
w2 = mtx[:, 14] / 1000.0
if H2OAbsModel.model in ['R19SD', 'R20SD', 'R22SD', 'R23', 'R23SD', 'R24']:
if H2OAbsModel.model in ['R19SD', 'R20SD', 'R22SD', 'R23', 'R23SD', 'R24', 'MWL24']:
w2s = mtx[:, 15] / 1000.0
if H2OAbsModel.model in ['R21SD', 'R22SD', 'R23', 'R23SD', 'R24']:
if H2OAbsModel.model in ['R21SD', 'R22SD', 'R23', 'R23SD', 'R24', 'MWL24']:
xw2 = mtx[:, 15]
w2s = mtx[:, 16] / 1000.0
xw2s = mtx[:, 17]
Expand Down
154 changes: 138 additions & 16 deletions pyrtlib/absorption_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,9 +195,57 @@ def liquid_water_absorption(water: np.ndarray, freq: np.ndarray, temp: np.ndarra
class N2AbsModel(AbsModel):
"""This class contains the absorption model used in pyrtlib.
"""

def n2_absor_mwl24(self, f, p, t):
"""
Calculation of the collision-induced absorption by N2-N2 molecular pairs
Based on the Classic Trajectory Calcs by Vigasin/Chistikov/Finenko
Args:
f (float): frequency (Ghz)
p (float): dry air pressure (mbar) (recalculations to Torr are used inside)
t (float): air temperature (K)
Returns:
float: N2-N2 absorption coefficient in 1/cm
@staticmethod
def n2_absorption(t: np.ndarray, p: np.ndarray, f: np.ndarray) -> np.ndarray:
Convert to Np/km by multiplying by 1.E5
Convert to dry air absorption by multiplying by 0.84 (see [Meshkov, DeLucia 2007 10.1016/j.jqsrt.2007.04.001])
References:
[Serov et al, JQSRT 2024]
"""
t0 = 300.
ti = t0 / t
p_torr = p * 0.7500616827

kbcm = 0.6950354549 # boltzmann in (1/cm)/K

d0 = 108. / (kbcm * t)
d_ref = 108. / (kbcm * 300.)
potential = (np.exp(d0) - (1 + d0)) / (np.exp(d_ref) - (1 + d_ref))

a_e = [1.4359, -0.0005385, 1.6182e-7, 1.8845e-10]
b_e = [0.07988, -0.0005165, 1.131e-6, -5.738e-10]
a0e, a1e, a2e, a3e = [1.754e-18, 1.033, 0.1510, 0.1420]
f1e, f2e, f3e, f4e, f5e = [44., 568., 358., 127., 92.]

m0 = 0.
aa = 0.

t2 = a3e * np.exp(-(f/f5e - 1)**2)
t1_denom = 1. + ((f - f1e)/f2e)**2
t1_num_denom = 1. + ((f - f3e)/f4e)**2
t1_num = 1 + a2e / t1_num_denom
spec_fun = a0e * 0.5 * (1 + a1e * t1_num / t1_denom + t2)

for i in range(4):
m0 += a_e[i] * f**i
aa += b_e[i] * f**i

t_dep = ti ** (m0 + aa * np.log(ti))

return spec_fun * t_dep * potential * p_torr**2 * f**2

def n2_absorption(self, t: np.ndarray, p: np.ndarray, f: np.ndarray) -> np.ndarray:
"""Collision-Induced Power Absorption Coefficient (Neper/km) in air
with modification of 1.34 to account for O2-O2 and O2-N2 collisions, as calculated by [Boissoles-2003]_.
Expand Down Expand Up @@ -227,14 +275,19 @@ def n2_absorption(t: np.ndarray, p: np.ndarray, f: np.ndarray) -> np.ndarray:
l, m, n = 9.95e-14, 3.22, 1
elif N2AbsModel.model == 'R03':
l, m, n = 6.5e-14, 3.6, 1.29
elif N2AbsModel.model == 'MWL24':
l, m, n = 9.95e-14, 3.22, 0.84
elif N2AbsModel.model in ['R98']:
l, m, n = 6.4e-14, 3.55, 1
fdepen = 1
else:
raise ValueError(
'[AbsN2] No model available with this name: {} . Sorry...'.format(N2AbsModel.model))

bf = l * fdepen * p * p * f * f * th ** m

if N2AbsModel.model == 'MWL24':
bf = self.n2_absor_mwl24(f, p, t) * 1.E5
else:
bf = l * fdepen * p * p * f * f * th ** m

abs_n2 = n * bf

Expand Down Expand Up @@ -267,7 +320,7 @@ def set_ll() -> None:
if H2OAbsModel.model not in H2OAbsModel.implemented_models()['WaterVapour']:
raise AbsModelError(
H2OAbsModel.model, f"Model {H2OAbsModel.model} is not available. It is necessary to define water vapour absorption model manually")
H2OAbsModel.h2oll = import_lineshape("h2oll")
H2OAbsModel.h2oll = import_lineshape("h2oll")

def h2o_continuum(self, frq: np.ndarray, vx: np.ndarray, nfreq: int):
nf = 6
Expand All @@ -293,9 +346,67 @@ def h2o_continuum(self, frq: np.ndarray, vx: np.ndarray, nfreq: int):
b = 0.5*p*(1.-p)
b1 = b*(1.-p)
b2 = b*p
cs[i] = -a[j]*b1+a[j+1]*(1.-c+b2)+a[j+2]*(c+b1)-a[j+3]*b2

cs[i] = -a[j]*b1+a[j+1]*(1.-c+b2)+a[j+2]*(c+b1)-a[j+3]*b2
return cs

def h2o_continuum_mwl24(self, frq: np.ndarray, vx: np.ndarray):
"""
H2O self-continuum absorption normalized to the squared water vapor pressure (1/cm)/mbar**2
Based on the known measurements meta-analysis and theoretical calculations
Args:
frq (float): Frequency (GHz)
vx (float): Theta (adim) - (normalised temperature 300/t(K))
Returns:
float: self-continuum term ((1/cm)/mbar**2), multiply by pvap**2 * 1.E5 to get Np/km
Reference:
[Tretyakov, Galanina, Koroleva et al 2024] (not published currently)
[Odintsova 2022 10.1016/j.jms.2022.111603]
[Galanina 2022 10.1016/j.jms.2022.111691]
"""

atm2mbar = 1013.25
t = 300.0 / vx
ti_bd = 268.0 / t
ti_md = 266.0 / t
ti_fw = 296.0 / t
# constants of bound dimers absorption (A0-5) and other components (see Tretyakov et al 2024)
(A0, A1, A2, A3, A4, A5, N_D, C_W, N_BD, N_MD, N_FW, PF_FW) = (5.55E-8, 434.8,
219.5, 4.93E-10,
21.06, 1.34E-14,
0.7, 8.5E-15,
10.0, 2.9, 1.5, 2.2)
# calculations for dimer equilibrium constants
# second virial coef-t parameters
b_svc = (-7.804242E+6, 8.345651E+4, -4.212794E+2, 1.242946, -2.409822E-3,
3.017768E-6, -2.518957E-9, 1.350628E-12, -4.134191E-16, 5.530774E-20)
b_t = b_svc[0]
for i in range(1, len(b_svc)):
b_t += b_svc[i] * t**i
b_t *= -(100./t)**6
r_t = 82.05746 * t # molar gas constant multiplied by temperature
k_bd = 4.7856E-4 * np.exp(1669.8 / t - 5.10485E-3 * t ) # K_BD from paper, in 1/atm
k_md = ((b_t - 30.5) / r_t - k_bd) / atm2mbar # K_MD, in 1/mbar
k_bd = k_bd / atm2mbar # recalc to 1/mbar

# bound dimers contribution

CON_BD = A0 / (A1 * A1 + (frq - A2) * (frq - A2))
CON_BD += A3 / (A4 * A4 + frq * frq) + A5
CON_BD = CON_BD * ti_bd**N_BD

# metastable dimers contribution
S2MON = 1.E-13 * (4.06 + 7.12E-3 * frq) * (1 / atm2mbar)
CON_MD = (S2MON * (1 - N_D) * ti_md**N_MD + N_D * CON_BD / k_bd) * k_md

# far wings contribution (note that it uses not f**2 but some different power)
CON_FW = C_W * frq**PF_FW * ti_fw**N_FW

cs_mwl24 = ((CON_BD + CON_MD) * frq**2 + CON_FW)

return cs_mwl24 * 1.E5

def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray, frq: np.ndarray, amu: Optional[dict] = None) -> Union[
Tuple[np.ndarray, np.ndarray], None]:
Expand Down Expand Up @@ -377,7 +488,7 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray,
rho = ekpa * 10.0 / (rvap * t)
f = frq
# cyh ***********************************************

if rho.any() <= 0.0:
npp = 0
ncpp = 0
Expand All @@ -386,7 +497,7 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray,
pvap = (rho * t) / 216.68
if H2OAbsModel.model in ['R03', 'R16', 'R17', 'R98']:
pvap = (rho * t) / 217.0
if H2OAbsModel.model in ['R22SD', 'R23SD', 'R24']:
if H2OAbsModel.model in ['R22SD', 'R23SD', 'R24', 'MWL24']:
pvap = (constants("Rwatvap")[0] * 1e-05) * rho * t
pda = p - pvap
if H2OAbsModel.model in ['R03', 'R16', 'R98']:
Expand All @@ -399,32 +510,43 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray,
if H2OAbsModel.model in ['R03', 'R98']:
con = (5.43e-10 * pda * ti ** 3 + 1.8e-08 *
pvap * ti ** 7.5) * pvap * f * f
elif H2OAbsModel.model == 'MWL24':
con_self = self.h2o_continuum_mwl24(f, vx) * pvap * pvap
# con_frgn = self.h2oll.cf * pda * pvap * f * f * ti**self.h2oll.xcf # foreign continuum from eariler version
con_frgn = 7.1e-10 * (300./t)**4.4 * f**1.96 * pda * pvap
con = con_self + con_frgn
else:

con = (self.h2oll.cf * pda * ti ** self.h2oll.xcf + self.h2oll.cs * pvap * ti ** self.h2oll.xcs) * \
pvap * f * f
# 2019/03/18 *********************************************************
# add resonances
nlines = len(self.h2oll.fl)



ti = self.h2oll.reftline / t
df = np.zeros((2, 1))

if H2OAbsModel.model.startswith(('R19SD', 'R20SD', 'R21SD', 'R22SD', 'R23SD', 'R24')):
if H2OAbsModel.model.startswith(('R19SD', 'R20SD', 'R21SD', 'R22SD', 'R23SD', 'R24', 'MWL24')):
tiln = np.log(ti)
ti2 = np.exp(2.5 * tiln)
summ = 0.0
if H2OAbsModel.model in ['R23SD', 'R24']:
if self.h2oll.cs > 0:

# npp_cs = np.zeros(1)
con = self.h2oll.cs * ti * self.h2oll.xcs
# for i in len(frq):
npp_cs = con
else:

npp_cs = self.h2o_continuum(frq, vx, 1)
for i in range(0, nlines):
for i in range(0, nlines):
width0 = self.h2oll.w0[i] * pda * ti ** self.h2oll.x[i] + \
self.h2oll.w0s[i] * pvap * ti ** self.h2oll.xs[i]
width2 = self.h2oll.w2[i] * pda + self.h2oll.w2s[i] * pvap
if H2OAbsModel.model in ['R21SD', 'R22SD', 'R23SD', 'R24']:
width2 = self.h2oll.w2[i] * pda + self.h2oll.w2s[i] * pvap
if H2OAbsModel.model in ['R21SD', 'R22SD', 'R23SD', 'R24', 'MWL24']:
if self.h2oll.w2[i] > 0:
width2 = self.h2oll.w2[i] * pda * ti ** self.h2oll.xw2[i] + self.h2oll.w2s[i] * pvap * ti ** \
self.h2oll.xw2s[i]
Expand All @@ -443,7 +565,7 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray,
df[0] = f - self.h2oll.fl[i] - shift
df[1] = f + self.h2oll.fl[i] + shift
base = width0 / (562500.0 + wsq)
if H2OAbsModel.model in ["R21SD", 'R22SD', 'R23SD', 'R24']:
if H2OAbsModel.model in ["R21SD", 'R22SD', 'R23SD', 'R24', 'MWL24']:
delta2 = self.h2oll.d2[i] * pda + self.h2oll.d2s[i] * pvap
res = 0.0
for j in range(0, 2):
Expand All @@ -459,15 +581,15 @@ def h2o_absorption(self, pdrykpa: np.ndarray, vx: np.ndarray, ekpa: np.ndarray,
delta2 = 0.0
xc = complex((width0 - np.dot(1.5, width2)), df[j] + np.dot(1.5, delta2)) / complex(
width2, -delta2)
elif H2OAbsModel.model in ["R21SD", 'R22SD', 'R23SD', 'R24']:
elif H2OAbsModel.model in ["R21SD", 'R22SD', 'R23SD', 'R24', 'MWL24']:
xc = complex(
(width0 - 1.5 * width2), df[j] + 1.5 * delta2) / complex(width2, -delta2)

xrt = np.sqrt(xc)
pxw = 1.77245385090551603 * xrt * \
_dcerror(-np.imag(xrt), np.real(xrt))
sd = 2.0 * (1.0 - pxw) / (
width2 if H2OAbsModel.model not in ['R20SD', 'R21SD', 'R22SD', 'R23SD', 'R24'] else complex(width2, -delta2))
width2 if H2OAbsModel.model not in ['R20SD', 'R21SD', 'R22SD', 'R23SD', 'R24', 'MWL24'] else complex(width2, -delta2))
res += np.real(sd) - base
elif np.abs(df[j]) < 750.0:
res += width0 / (df[j] ** 2 + wsq) - base
Expand Down
2 changes: 1 addition & 1 deletion pyrtlib/rt_equation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -555,7 +555,7 @@ def clearsky_absorption(p: np.ndarray, t: np.ndarray, e: np.ndarray, frq: np.nda
aO2[i] = factor * (npp + ncpp) * db2np
# add N2 term
if N2AbsModel.model not in ['R03', 'R16', 'R17', 'R18', 'R98']:
aN2[i] = N2AbsModel.n2_absorption(
aN2[i] = N2AbsModel().n2_absorption(
t[i], np.dot(pdrykpa, 10), frq)

if not N2AbsModel.model:
Expand Down

0 comments on commit fe96776

Please sign in to comment.