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_fieldSetPlanetEllipsoid(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_spectral_radiance_fieldSunlessClearsky(
45    longitude=lon,
46    latitude=lat,
47    disort_quadrature_dimension=NQuad,
48    disort_legendre_polynomial_dimension=1,
49    disort_fourier_mode_dimension=1,
50)
51
52# %% Equivalent ARTS calculations
53
54ws.ray_pathGeometricDownlooking(
55    latitude=lat,
56    longitude=lon,
57    max_step=1000.0,
58)
59ws.spectral_radianceClearskyEmission()
60
61# %% Plot results
62
63if PLOT:
64    f = ws.frequency_grid - line_f0
65
66    plt.semilogy(
67        f,
68        ws.disort_spectral_radiance_field[:, 0, 0, : (NQuad // 2)],
69        label="disort",
70    )
71    plt.semilogy(f, ws.spectral_radiance[:, 0], "k--", lw=3)
72    plt.semilogy(
73        f,
74        ws.disort_spectral_radiance_field[:, 0, 0, 0],
75        "g:",
76        lw=3,
77    )
78    plt.semilogy(
79        f,
80        ws.disort_spectral_radiance_field[:, 0, 0, (NQuad // 2) - 1],
81        "m:",
82        lw=3,
83    )
84    plt.ylabel("Spectral radiance [W sr$^{-1}$ m$^{-2}$ Hz$^{-1}$]")
85    plt.xlabel("Dirac frequency [index count]")
86    plt.title("Downlooking")
87
88# %% The last test should be that we are close to the correct values
89
90assert np.allclose(
91    ws.disort_spectral_radiance_field[:, 0, 0, NQuad // 2-1]
92    / ws.spectral_radiance[:, 0],
93    1,
94    rtol=1e-3,
95), "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_fieldSetPlanetEllipsoid(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_agendaSet(option="SunlessClearsky")
45
46ws.ray_pathGeometricDownlooking(longitude=lon, latitude=lat, max_step=40_000)
47ws.ray_path_atmospheric_pointFromPath()
48ws.ray_path_frequency_gridFromPath()
49ws.ray_path_propagation_matrixFromPath()
50ws.ray_path_pointLowestFromPath()
51
52ws.disort_spectral_flux_fieldFromAgenda(
53    disort_quadrature_dimension=NQuad,
54    disort_legendre_polynomial_dimension=1,
55    disort_fourier_mode_dimension=1,
56)
57
58assert np.allclose(
59    ws.disort_spectral_flux_field[:, :2].flatten()
60    / np.array(
61        [
62            2.65924430e-15,
63            2.66162733e-15,
64            2.75440515e-15,
65            9.57958182e-18,
66            2.08076327e-17,
67            5.39749082e-16,
68            2.93074072e-15,
69            2.93357370e-15,
70            3.03918211e-15,
71            9.99616600e-18,
72            2.34511632e-17,
73            6.12773186e-16,
74            3.19264874e-15,
75            3.19665668e-15,
76            3.33784028e-15,
77            1.03768102e-17,
78            3.05943372e-17,
79            8.07853130e-16,
80            3.35251230e-15,
81            3.36075470e-15,
82            3.65036247e-15,
83            1.07209025e-17,
84            6.78926741e-17,
85            1.56163327e-15,
86            3.37235023e-15,
87            2.86992821e-15,
88            3.97673151e-15,
89            1.52446216e-15,
90            2.80896759e-15,
91            3.94566396e-15,
92        ]
93    ),
94    1,
95), "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_fieldSetPlanetEllipsoid(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    ws.ray_path_atmospheric_pointFromPath()
62
63    ws.ray_path_frequency_gridFromPath()
64    ws.ray_path_propagation_matrixFromPath()
65    ws.disort_settingsInit()
66    ws.disort_settingsOpticalThicknessFromPath()
67    ws.disort_settingsLayerThermalEmissionLinearInTau()
68    ws.disort_settingsSurfaceEmissionByTemperature(ray_path_point=ws.ray_path[0])
69    ws.disort_settingsCosmicMicrowaveBackgroundRadiation()
70    ws.disort_settingsSurfaceLambertian(value=surface_reflectivity)
71    ws.disort_settingsNoSingleScatteringAlbedo()
72    ws.disort_settingsNoFractionalScattering()
73    ws.disort_settingsNoLegendre()
74    ws.disort_settingsSetSun(ray_path_point=ws.ray_path[-1])
75    ws.disort_spectral_flux_fieldCalc()
76
77    plt.semilogy(pyarts.arts.convert.freq2kaycm(ws.frequency_grid),
78                ws.disort_spectral_flux_field[:, 1])
79
80    f, s = pyarts.plots.AtmField.plot(ws.atmospheric_field,
81                            alts=np.linspace(0, ws.atmospheric_field.top_of_atmosphere))