unphysically large absorption and emission when "dry" layer is present near surface for emission spectra and scattering is turned on
I am encountering what appears to be a bug or incorrect calculation when an Earth-like atmosphere (such as the one in the surface scattering example) has a dry layer beneath a moist layer, and scattering emission is turned on. When this is the case, the absorption bands feature unphysical absorption and emission lines (as in down to -10^150 and up to infinity). Below is a minimal (non-)working example, simply taking the surface scattering example and making the bottom 5 layers devoid of H2O (this does not seem to depend on which H2O line list is used):
os.environ["pRT_input_data_path"] = "/home/adiv/astro/pRT_opacities"
from petitRADTRANS import Radtrans
from petitRADTRANS import nat_cst as nc
import numpy as np
import matplotlib.pyplot as plt
atmosphere = Radtrans(line_species = ['H2O_HITEMP',
'CO2',
'O3'],
rayleigh_species = ['N2', 'O2'],
continuum_opacities = ['N2-N2', 'N2-O2','O2-O2','CO2-CO2'],
wlen_bords_micron = [0.3, 20],
do_scat_emis = True)
pressures = np.logspace(-6, 0, 100)
atmosphere.setup_opa_structure(pressures)
R_pl = nc.r_earth
gravity = nc.G*(nc.m_earth)/R_pl**2
# P-T parameters
kappa_IR = 0.0009
gamma = 0.01
T_int = 250.
T_equ = 220.
temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ)
# Mean molecular weight of Earth's atmosphere
MMW = 28.7 * np.ones_like(temperature)
# Absorber mass fractions
mass_fractions = {}
mass_fractions['N2'] = 0.78 * 28. / MMW * np.ones_like(temperature)
mass_fractions['O2'] = 0.21 * 32. / MMW* np.ones_like(temperature)
mass_fractions['H2O_HITEMP'] = 1.0e-3 * 18. / MMW * np.ones_like(temperature)
mass_fractions['H2O_HITEMP'][-5:] = 0.0
mass_fractions['O3'] = 1e-7 * 48. / MMW * np.ones_like(temperature)
mass_fractions['CO2'] = 0.0004 * 44. / MMW * np.ones_like(temperature)
atmosphere.reflectance = np.ones_like(atmosphere.freq)*0.3
atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW)
plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.freq*atmosphere.flux/1e-6)
plt.xscale('log')
plt.yscale('log')
plt.show()
I've attached a plot of what this looks like.
Is there a plausible fix for this?