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)

_images/1-clearsky-radiance.svg

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"

(Source code)

Clearsky flux from datase

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"

xarr = pyarts.data.xarray_open_dataset("atm.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()

(Source code)

_images/3-clearsky-flux-from-dataset_00.svg

Fig. 1 (svg, pdf)

_images/3-clearsky-flux-from-dataset_01.svg

Fig. 2 (svg, pdf)