3. Disort ========= Clearsky radiance ----------------- .. code-block:: python :name: Clearsky radiance :linenos: import pyarts import numpy as np import matplotlib.pyplot as plt PLOT = False 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, 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_unit = "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 if PLOT: f = ws.frequency_grid - line_f0 plt.semilogy( f, ws.disort_spectral_radiance_field[:, 0, 0, : (NQuad // 2)], label="disort", ) plt.semilogy(f, ws.spectral_radiance[:, 0], "k--", lw=3) plt.semilogy( f, ws.disort_spectral_radiance_field[:, 0, 0, 0], "g:", lw=3, ) plt.semilogy( f, ws.disort_spectral_radiance_field[:, 0, 0, (NQuad // 2) - 1], "m:", lw=3, ) plt.ylabel("Spectral radiance [W sr$^{-1}$ m$^{-2}$ Hz$^{-1}$]") plt.xlabel("Dirac frequency [index count]") plt.title("Downlooking") # %% The last test should be that we are close to the correct values assert np.allclose( ws.disort_spectral_radiance_field[:, 0, 0, NQuad // 2-1] / ws.spectral_radiance[:, 0], 1, rtol=1e-3, ), "Bad results, clearsky calculations are not close between DISORT and ARTS" Clearsky flux ------------- .. code-block:: python :name: Clearsky flux :linenos: import pyarts import numpy as np import matplotlib.pyplot as plt 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_unit = "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( ws.disort_spectral_flux_field[:, :2].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 -------------------------- .. code-block:: python :name: Clearsky flux from dataset :linenos: import pyarts import numpy as np import xarray as xa from dataclasses import dataclass import matplotlib.pyplot as plt if not pyarts.arts.globals.data.is_lgpl: 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() plt.semilogy(pyarts.arts.convert.freq2kaycm(ws.frequency_grid), ws.disort_spectral_flux_field[:, 1]) f, s = pyarts.plots.AtmField.plot(ws.atmospheric_field, alts=np.linspace(0, ws.atmospheric_field.top_of_atmosphere))