Disort
Clearsky radiance
import os
import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
# Download catalogs
pyarts.data.download()
toa = 120e3
lat = 0
lon = 0
NQuad = 40
ws = pyarts.Workspace()
# %% Sampled frequency range
line_f0 = 118750348044.712
ws.frequency_grid = np.linspace(-20e9, 2e6, 101) + line_f0
# %% Species and line absorption
ws.absorption_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
# %% Use the automatic agenda setter for propagation matrix calculations
ws.propagation_matrix_agendaAuto()
# %% Grids and planet
ws.surface_fieldPlanet(option="Earth")
ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atmospheric_fieldRead(
toa=toa, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
# %% Checks and settings
ws.spectral_radiance_transform_operatorSet(option="Tb")
ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
ws.spectral_radiance_surface_agendaSet(option="Blackbody")
# %% Core Disort calculations
ws.disort_settings_agendaSetup()
ws.disort_spectral_radiance_fieldProfile(
longitude=lon,
latitude=lat,
disort_quadrature_dimension=NQuad,
disort_legendre_polynomial_dimension=1,
disort_fourier_mode_dimension=1,
)
# %% Equivalent ARTS calculations
ws.ray_pathGeometricDownlooking(
latitude=lat,
longitude=lon,
max_step=1000.0,
)
ws.spectral_radianceClearskyEmission()
# %% Plot results
f = ws.frequency_grid - line_f0
fig, ax = plt.subplots()
ax.semilogy(
f,
ws.disort_spectral_radiance_field.data[:, 0, 0, : (NQuad // 2)],
label="disort",
)
ax.semilogy(f, ws.spectral_radiance[:, 0], "k--", lw=3)
ax.semilogy(
f,
ws.disort_spectral_radiance_field.data[:, 0, 0, (NQuad // 2) - 1],
"g:",
lw=3,
)
ax.semilogy(
f,
ws.disort_spectral_radiance_field.data[:, 0, 0, 0],
"m:",
lw=3,
)
ax.set_ylabel("Spectral radiance [W sr$^{-1}$ m$^{-2}$ Hz$^{-1}$]")
ax.set_xlabel("Dirac frequency [index count]")
ax.set_title("Downlooking")
if "ARTS_HEADLESS" not in os.environ:
plt.show()
# %% The last test should be that we are close to the correct values
assert np.allclose(
ws.disort_spectral_radiance_field.data[:, 0, 0, 0] / ws.spectral_radiance[:, 0],
1,
rtol=1e-3,
), "Bad results, clearsky calculations are not close between DISORT and ARTS"
(Source code, svg, pdf)
Clearsky flux
import numpy as np
import pyarts3 as pyarts
# Download catalogs
pyarts.data.download()
toa = 100e3
lat = 0
lon = 0
NQuad = 40
ws = pyarts.Workspace()
# %% Sampled frequency range
line_f0 = 118750348044.712
ws.frequency_grid = [line_f0]
ws.frequency_grid = np.linspace(-20e9, 2e6, 5) + line_f0
# %% Species and line absorption
ws.absorption_speciesSet(species=["O2-66"])
ws.ReadCatalogData()
# %% Use the automatic agenda setter for propagation matrix calculations
ws.propagation_matrix_agendaAuto()
# %% Grids and planet
ws.surface_fieldPlanet(option="Earth")
ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
ws.atmospheric_fieldRead(
toa=toa, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
)
ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
# ws.atmospheric_field[pyarts.arts.AtmKey.t] = 300.0
# %% Checks and settings
ws.spectral_radiance_transform_operatorSet(option="Tb")
ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
ws.spectral_radiance_surface_agendaSet(option="Blackbody")
# %% Core Disort calculations
ws.disort_settings_agendaSetup()
ws.ray_pathGeometricDownlooking(longitude=lon, latitude=lat, max_step=40_000)
ws.disort_spectral_flux_fieldFromAgenda(
disort_quadrature_dimension=NQuad,
disort_legendre_polynomial_dimension=1,
disort_fourier_mode_dimension=1,
)
assert np.allclose(
np.append(
ws.disort_spectral_flux_field.up,
ws.disort_spectral_flux_field.down_diffuse,
axis=1,
).flatten()
/ np.array(
[
2.65924430e-15,
2.66162733e-15,
2.75440515e-15,
9.57958182e-18,
2.08076327e-17,
5.39749082e-16,
2.93074072e-15,
2.93357370e-15,
3.03918211e-15,
9.99616600e-18,
2.34511632e-17,
6.12773186e-16,
3.19264874e-15,
3.19665668e-15,
3.33784028e-15,
1.03768102e-17,
3.05943372e-17,
8.07853130e-16,
3.35251230e-15,
3.36075470e-15,
3.65036247e-15,
1.07209025e-17,
6.78926741e-17,
1.56163327e-15,
3.37235023e-15,
2.86992821e-15,
3.97673151e-15,
1.52446216e-15,
2.80896759e-15,
3.94566396e-15,
]
),
1,
), "Mismatch from historical Disort spectral fluxes"
Clearsky flux from dataset
This example requires the following input file:
import os
import matplotlib.pyplot as plt
import numpy as np
import pyarts3 as pyarts
if pyarts.arts.globals.data.is_lgpl:
print(
"CKDMT models are not available in LGPL mode, compile with -DENABLE_ARTS_LGPL=0"
)
exit(0)
# Download catalogs
pyarts.data.download()
NQuad = 16
max_level_step = 1e3
atm_latitude = 0.0
atm_longitude = 0.0
solar_latitude = 0.0
solar_longitude = 0.0
surface_temperature = 293.0
surface_reflectivity = 0.05
cutoff = ["ByLine", 750e9]
remove_lines_percentile = 70
sunfile = "star/Sun/solar_spectrum_QUIET.xml"
planet = "Earth"
# Input file names are prefixed with the example number
prefix = "3-"
xarr = pyarts.data.xarray_open_dataset(prefix + "atmosphere.nc")
ws = pyarts.Workspace()
ws.frequency_grid = pyarts.arts.convert.kaycm2freq(np.linspace(500, 2500, 1001))
ws.atmospheric_field = pyarts.data.to_atmospheric_field(xarr)
v = pyarts.data.to_absorption_species(ws.atmospheric_field)
ws.absorption_species = v
ws.ReadCatalogData(ignore_missing=True)
ws.propagation_matrix_agendaAuto(T_extrapolfac=1e9)
for band in ws.absorption_bands:
ws.absorption_bands[band].cutoff = cutoff[0]
ws.absorption_bands[band].cutoff_value = cutoff[1]
ws.absorption_bands.keep_hitran_s(remove_lines_percentile)
ws.surface_fieldPlanet(option=planet)
sun = pyarts.arts.GriddedField2.fromxml(sunfile)
ws.surface_field["t"] = surface_temperature
ws.sunFromGrid(
sun_spectrum_raw=sun,
latitude=solar_latitude,
longitude=solar_longitude,
)
ws.disort_quadrature_dimension = NQuad
ws.disort_fourier_mode_dimension = 1
ws.disort_legendre_polynomial_dimension = 1
ws.ray_pathGeometricDownlooking(
latitude=atm_latitude,
longitude=atm_longitude,
max_step=max_level_step,
)
ws.disort_settings_agendaSetup(
sun_setting="Sun",
surface_setting="ThermalLambertian",
surface_lambertian_value=surface_reflectivity * np.ones_like(ws.frequency_grid),
)
ws.disort_spectral_flux_fieldFromAgenda()
fig, ax = plt.subplots()
ax.semilogy(
pyarts.arts.convert.freq2kaycm(ws.frequency_grid),
ws.disort_spectral_flux_field.down_diffuse,
)
f, s = pyarts.plots.AtmField.plot(
ws.atmospheric_field, alts=np.linspace(0, ws.atmospheric_field.top_of_atmosphere)
)
f.suptitle("Atmospheric field")
if "ARTS_HEADLESS" not in os.environ:
plt.show()
