3. Disort

Clearsky radiance

 1import pyarts
 2import numpy as np
 3import matplotlib.pyplot as plt
 4
 5PLOT = False
 6toa = 100e3
 7lat = 0
 8lon = 0
 9NQuad = 40
10
11ws = pyarts.Workspace()
12
13# %% Sampled frequency range
14
15line_f0 = 118750348044.712
16ws.frequency_grid = [line_f0]
17ws.frequency_grid = np.linspace(-20e9, 2e6, 101) + line_f0
18
19# %% Species and line absorption
20
21ws.absorption_speciesSet(species=["O2-66"])
22ws.ReadCatalogData()
23
24# %% Use the automatic agenda setter for propagation matrix calculations
25ws.propagation_matrix_agendaAuto()
26
27# %% Grids and planet
28
29ws.surface_fieldPlanet(option="Earth")
30ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
31ws.atmospheric_fieldRead(
32    toa=toa, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
33)
34ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
35
36# %% Checks and settings
37
38ws.spectral_radiance_unit = "Tb"
39ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
40ws.spectral_radiance_surface_agendaSet(option="Blackbody")
41
42# %% Core Disort calculations
43
44ws.disort_settings_agendaSetup()
45
46ws.disort_spectral_radiance_fieldProfile(
47    longitude=lon,
48    latitude=lat,
49    disort_quadrature_dimension=NQuad,
50    disort_legendre_polynomial_dimension=1,
51    disort_fourier_mode_dimension=1,
52)
53
54# %% Equivalent ARTS calculations
55
56ws.ray_pathGeometricDownlooking(
57    latitude=lat,
58    longitude=lon,
59    max_step=1000.0,
60)
61ws.spectral_radianceClearskyEmission()
62
63# %% Plot results
64
65if PLOT:
66    f = ws.frequency_grid - line_f0
67
68    plt.semilogy(
69        f,
70        ws.disort_spectral_radiance_field[:, 0, 0, : (NQuad // 2)],
71        label="disort",
72    )
73    plt.semilogy(f, ws.spectral_radiance[:, 0], "k--", lw=3)
74    plt.semilogy(
75        f,
76        ws.disort_spectral_radiance_field[:, 0, 0, 0],
77        "g:",
78        lw=3,
79    )
80    plt.semilogy(
81        f,
82        ws.disort_spectral_radiance_field[:, 0, 0, (NQuad // 2) - 1],
83        "m:",
84        lw=3,
85    )
86    plt.ylabel("Spectral radiance [W sr$^{-1}$ m$^{-2}$ Hz$^{-1}$]")
87    plt.xlabel("Dirac frequency [index count]")
88    plt.title("Downlooking")
89
90# %% The last test should be that we are close to the correct values
91
92assert np.allclose(
93    ws.disort_spectral_radiance_field[:, 0, 0, NQuad // 2-1]
94    / ws.spectral_radiance[:, 0],
95    1,
96    rtol=1e-3,
97), "Bad results, clearsky calculations are not close between DISORT and ARTS"

Clearsky flux

 1import pyarts
 2import numpy as np
 3import matplotlib.pyplot as plt
 4
 5toa = 100e3
 6lat = 0
 7lon = 0
 8NQuad = 40
 9
10ws = pyarts.Workspace()
11
12# %% Sampled frequency range
13
14line_f0 = 118750348044.712
15ws.frequency_grid = [line_f0]
16ws.frequency_grid = np.linspace(-20e9, 2e6, 5) + line_f0
17
18# %% Species and line absorption
19
20ws.absorption_speciesSet(species=["O2-66"])
21ws.ReadCatalogData()
22
23# %% Use the automatic agenda setter for propagation matrix calculations
24ws.propagation_matrix_agendaAuto()
25
26# %% Grids and planet
27
28ws.surface_fieldPlanet(option="Earth")
29ws.surface_field[pyarts.arts.SurfaceKey("t")] = 295.0
30ws.atmospheric_fieldRead(
31    toa=toa, basename="planets/Earth/afgl/tropical/", missing_is_zero=1
32)
33ws.atmospheric_fieldIGRF(time="2000-03-11 14:39:37")
34# ws.atmospheric_field[pyarts.arts.AtmKey.t] = 300.0
35
36# %% Checks and settings
37
38ws.spectral_radiance_unit = "Tb"
39ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
40ws.spectral_radiance_surface_agendaSet(option="Blackbody")
41
42# %% Core Disort calculations
43
44ws.disort_settings_agendaSetup()
45
46ws.ray_pathGeometricDownlooking(longitude=lon, latitude=lat, max_step=40_000)
47
48ws.disort_spectral_flux_fieldFromAgenda(
49    disort_quadrature_dimension=NQuad,
50    disort_legendre_polynomial_dimension=1,
51    disort_fourier_mode_dimension=1,
52)
53
54assert np.allclose(
55    ws.disort_spectral_flux_field[:, :2].flatten()
56    / np.array(
57        [
58            2.65924430e-15,
59            2.66162733e-15,
60            2.75440515e-15,
61            9.57958182e-18,
62            2.08076327e-17,
63            5.39749082e-16,
64            2.93074072e-15,
65            2.93357370e-15,
66            3.03918211e-15,
67            9.99616600e-18,
68            2.34511632e-17,
69            6.12773186e-16,
70            3.19264874e-15,
71            3.19665668e-15,
72            3.33784028e-15,
73            1.03768102e-17,
74            3.05943372e-17,
75            8.07853130e-16,
76            3.35251230e-15,
77            3.36075470e-15,
78            3.65036247e-15,
79            1.07209025e-17,
80            6.78926741e-17,
81            1.56163327e-15,
82            3.37235023e-15,
83            2.86992821e-15,
84            3.97673151e-15,
85            1.52446216e-15,
86            2.80896759e-15,
87            3.94566396e-15,
88        ]
89    ),
90    1,
91), "Mismatch from historical Disort spectral fluxes"

Clearsky flux from dataset

 1import pyarts
 2import numpy as np
 3import xarray as xa
 4from dataclasses import dataclass
 5import matplotlib.pyplot as plt
 6
 7if not pyarts.arts.globals.data.is_lgpl:
 8    NQuad = 16
 9    max_level_step = 1e3
10    atm_latitude = 0.0
11    atm_longitude = 0.0
12    solar_latitude = 0.0
13    solar_longitude = 0.0
14    surface_temperature = 293.0
15    surface_reflectivity = 0.05
16    cutoff = ["ByLine", 750e9]
17    remove_lines_percentile = 70
18    sunfile = "star/Sun/solar_spectrum_QUIET.xml"
19    planet = "Earth"
20
21    xarr = pyarts.data.xarray_open_dataset("atm.nc")
22
23    ws = pyarts.Workspace()
24
25    ws.frequency_grid = pyarts.arts.convert.kaycm2freq(np.linspace(500, 2500, 1001))
26    ws.atmospheric_field = pyarts.data.to_atmospheric_field(xarr)
27
28    v = pyarts.data.to_absorption_species(ws.atmospheric_field)
29
30    ws.absorption_species = v
31    ws.ReadCatalogData(ignore_missing=True)
32    ws.propagation_matrix_agendaAuto(T_extrapolfac=1e9)
33
34    for band in ws.absorption_bands:
35        ws.absorption_bands[band].cutoff = cutoff[0]
36        ws.absorption_bands[band].cutoff_value = cutoff[1]
37
38    ws.absorption_bands.keep_hitran_s(remove_lines_percentile)
39
40    ws.surface_fieldPlanet(option=planet)
41
42    sun = pyarts.arts.GriddedField2.fromxml(sunfile)
43
44    ws.surface_field["t"] = surface_temperature
45
46    ws.sunFromGrid(
47        sun_spectrum_raw=sun,
48        latitude=solar_latitude,
49        longitude=solar_longitude,
50    )
51
52    ws.disort_quadrature_dimension = NQuad
53    ws.disort_fourier_mode_dimension = 1
54    ws.disort_legendre_polynomial_dimension = 1
55
56    ws.ray_pathGeometricDownlooking(
57        latitude=atm_latitude,
58        longitude=atm_longitude,
59        max_step=max_level_step,
60    )
61
62    ws.disort_settings_agendaSetup(
63        sun_setting="Sun",
64        surface_setting="ThermalLambertian",
65        surface_lambertian_value=surface_reflectivity * np.ones_like(ws.frequency_grid),
66    )
67
68    ws.disort_spectral_flux_fieldFromAgenda()
69
70    plt.semilogy(pyarts.arts.convert.freq2kaycm(ws.frequency_grid),
71                ws.disort_spectral_flux_field[:, 1])
72
73    f, s = pyarts.plots.AtmField.plot(ws.atmospheric_field,
74                            alts=np.linspace(0, ws.atmospheric_field.top_of_atmosphere))