3. Disort

Clearsky radiance

 1import pyarts3 as pyarts
 2import numpy as np
 3import matplotlib.pyplot as plt
 4
 5PLOT = False
 6toa = 120e3
 7lat = 0
 8lon = 0
 9NQuad = 40
10
11ws = pyarts.Workspace()
12
13# %% Sampled frequency range
14
15line_f0 = 118750348044.712
16ws.frequency_grid = np.linspace(-20e9, 2e6, 101) + 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
35# %% Checks and settings
36
37ws.spectral_radiance_transform_operatorSet(option="Tb")
38ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
39ws.spectral_radiance_surface_agendaSet(option="Blackbody")
40
41# %% Core Disort calculations
42
43ws.disort_settings_agendaSetup()
44
45ws.disort_spectral_radiance_fieldProfile(
46    longitude=lon,
47    latitude=lat,
48    disort_quadrature_dimension=NQuad,
49    disort_legendre_polynomial_dimension=1,
50    disort_fourier_mode_dimension=1,
51)
52
53# %% Equivalent ARTS calculations
54
55ws.ray_pathGeometricDownlooking(
56    latitude=lat,
57    longitude=lon,
58    max_step=1000.0,
59)
60ws.spectral_radianceClearskyEmission()
61
62# %% Plot results
63
64if PLOT:
65    f = ws.frequency_grid - line_f0
66
67    plt.semilogy(
68        f,
69        ws.disort_spectral_radiance_field.data[:, 0, 0, : (NQuad // 2)],
70        label="disort",
71    )
72    plt.semilogy(f, ws.spectral_radiance[:, 0], "k--", lw=3)
73    plt.semilogy(
74        f,
75        ws.disort_spectral_radiance_field.data[:, 0, 0, (NQuad // 2) - 1],
76        "g:",
77        lw=3,
78    )
79    plt.semilogy(
80        f,
81        ws.disort_spectral_radiance_field.data[:, 0, 0, 0],
82        "m:",
83        lw=3,
84    )
85    plt.ylabel("Spectral radiance [W sr$^{-1}$ m$^{-2}$ Hz$^{-1}$]")
86    plt.xlabel("Dirac frequency [index count]")
87    plt.title("Downlooking")
88
89# %% The last test should be that we are close to the correct values
90
91assert np.allclose(
92    ws.disort_spectral_radiance_field.data[:, 0, 0, 0]
93    / ws.spectral_radiance[:, 0],
94    1,
95    rtol=1e-3,
96), "Bad results, clearsky calculations are not close between DISORT and ARTS"

Clearsky flux

 1import pyarts3 as 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_transform_operatorSet(option="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    np.append(ws.disort_spectral_flux_field.up,
56              ws.disort_spectral_flux_field.down_diffuse, axis=1).flatten()
57    / np.array(
58        [
59            2.65924430e-15,
60            2.66162733e-15,
61            2.75440515e-15,
62            9.57958182e-18,
63            2.08076327e-17,
64            5.39749082e-16,
65            2.93074072e-15,
66            2.93357370e-15,
67            3.03918211e-15,
68            9.99616600e-18,
69            2.34511632e-17,
70            6.12773186e-16,
71            3.19264874e-15,
72            3.19665668e-15,
73            3.33784028e-15,
74            1.03768102e-17,
75            3.05943372e-17,
76            8.07853130e-16,
77            3.35251230e-15,
78            3.36075470e-15,
79            3.65036247e-15,
80            1.07209025e-17,
81            6.78926741e-17,
82            1.56163327e-15,
83            3.37235023e-15,
84            2.86992821e-15,
85            3.97673151e-15,
86            1.52446216e-15,
87            2.80896759e-15,
88            3.94566396e-15,
89        ]
90    ),
91    1,
92), "Mismatch from historical Disort spectral fluxes"

Clearsky flux from dataset

 1import pyarts3 as 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.down_diffuse)
72
73    f, s = pyarts.plots.AtmField.plot(ws.atmospheric_field,
74                            alts=np.linspace(0, ws.atmospheric_field.top_of_atmosphere))