Disort

  1#!/usr/bin/env python
  2# coding: utf-8
  3
  4# # Scattering calculations with DISORT
  5#
  6# This notebook demonstrates how to simulate brightness temperatures of cloudy skies using the DISORT scattering solver.
  7
  8# In[ ]:
  9
 10
 11import pyarts3 as pyarts
 12import numpy as np
 13import matplotlib.pyplot as plt
 14
 15toa = 100e3
 16lat = 0
 17lon = 0
 18NQuad = 40
 19ws = pyarts.Workspace()
 20
 21
 22# In[ ]:
 23
 24
 25pyarts.arts.globals.omp_set_num_threads(1)
 26
 27
 28# In[ ]:
 29
 30
 31ws.frequency_grid = [31.5e9, 165e9, 666e9]
 32
 33# %% Species and line absorption
 34ws.absorption_speciesSet(species=["N2-SelfContStandardType", "O2-PWR98", "H2O-PWR98"])
 35ws.ReadCatalogData()
 36ws.propagation_matrix_agendaAuto()
 37
 38
 39# ## Defining scattering species
 40#
 41# ### Loading scattering data
 42#
 43# We load scattering and scattering meta data from the ARTS .xml format. The scattering and scattering meta data arrays contain scattering data for rain particles (index 0) and ice particle (index 1).
 44
 45# In[ ]:
 46
 47
 48from pyarts3.arts import ParticleHabit
 49from pyarts3.xml import load
 50scat_data_raw = load("scat_data.xml")
 51scat_data_meta = load("scat_meta.xml" )
 52
 53t_grid = scat_data_raw[0][0].T_grid
 54f_grid = scat_data_raw[0][0].f_grid
 55
 56
 57# Next we define a rain habit that holds the scattering data for liquid cloud and rain drops of different sizes.
 58#
 59# > **Note**: We transform the rain habit to the spectral TRO representation that is expected by DISORT already here.
 60
 61# In[ ]:
 62
 63
 64from pyarts3.arts import ParticleHabit
 65rain_habit = ParticleHabit.from_legacy_tro(scat_data_raw[0], scat_data_meta[0])
 66rain_habit = rain_habit.to_tro_spectral(t_grid, f_grid, 39)
 67
 68
 69# The particle habit, which holds the scattering data, needs to be combined with a PSD so that it can be used as a scattering species in ARTS. The PSD needs to be linked to an atmospheric field through a particule property. The particulate property represents the moments of the PSD that vary throughout the atmosphere.
 70#
 71# For the rain example considered here, we use a single moment corresponding to the mass density of the rain.
 72
 73# In[ ]:
 74
 75
 76from pyarts3.arts import MGDSingleMoment, ScatteringSpeciesProperty, ParticulateProperty, ScatteringHabit
 77rain_first_moment = pyarts.arts.ScatteringSpeciesProperty("rain", pyarts.arts.ParticulateProperty("MassDensity"))
 78psd = MGDSingleMoment(rain_first_moment, "Wang16", 270, 300, False)
 79rain = ScatteringHabit(rain_habit, psd)
 80
 81
 82# In[ ]:
 83
 84
 85ws.scattering_species = [rain]
 86
 87
 88# ### Demonstration
 89#
 90# The scattering habit can now be used to calculate the bulk properties of our rain scattering species for any point in the atmosphere.
 91
 92# In[ ]:
 93
 94
 95from pyarts3.arts import AtmPoint
 96point = AtmPoint()
 97point["t"] = 280
 98point[rain_first_moment] = 1e-4
 99
100
101# In[ ]:
102
103
104bulk_props = rain.get_bulk_scattering_properties_tro_spectral(point, f_grid, 1)
105
106
107# ## Grids and Planet
108
109# In[ ]:
110
111
112from pyarts3.xml import load
113
114p_grid = load("p_grid.xml")
115t_field = load("t_field.xml")
116z_field = load("z_field.xml")
117vmr_field = load("vmr_field.xml")
118pbf_field = load("particle_bulkprop_field.xml")
119pbf_names = load("particle_bulkprop_names.xml")
120ws.surface_fieldPlanet(option="Earth")
121ws.surface_field[pyarts.arts.SurfaceKey("t")] = t_field[0, 0, 0]
122
123
124# In[ ]:
125
126
127from pyarts3.arts import GriddedField3, Tensor3, Vector
128
129lat_grid = np.array([0.0])
130lon_grid = np.array([0.0])
131z_grid = z_field[..., 0, 0]
132
133pressure = GriddedField3("p", p_grid[..., None, None], ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
134temperature = GriddedField3("t", t_field, ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
135n2 = GriddedField3("N2", vmr_field[0], ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
136o2 = GriddedField3("O2", vmr_field[1], ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
137h2o = GriddedField3("H2O", vmr_field[2], ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
138rwc = GriddedField3("RWC", 1.0 * pbf_field[0], ["alt", "lon", "lat"], (z_grid, lon_grid, lat_grid))
139
140
141# In[ ]:
142
143
144ws.atmospheric_field["p"] = pressure
145ws.atmospheric_field["t"] = temperature
146ws.atmospheric_field["N2"] = n2
147ws.atmospheric_field["O2"] = o2
148ws.atmospheric_field["H2O"] = h2o
149ws.atmospheric_field[rain_first_moment] = rwc
150ws.atmospheric_field.top_of_atmosphere = 12.0e3
151
152
153# In[ ]:
154
155
156ws.propagation_matrix_scattering_spectral_agenda
157
158
159# ## Checks and settings
160
161# In[ ]:
162
163
164ws.spectral_radiance_transform_operatorSet(option="Tb")
165ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
166ws.spectral_radiance_surface_agendaSet(option="Blackbody")
167
168ws.disort_settings_agendaSetup(scattering_setting="ScatteringSpecies")
169ws.disort_quadrature_dimension = 40
170ws.disort_fourier_mode_dimension = 1
171ws.propagation_matrix_scattering_spectral_agendaSet()
172ws.disort_legendre_polynomial_dimension = 40
173ws.spectral_radiance_transform_operatorSet(option="Tb")
174ws.spectral_radiance_space_agendaSet(option="UniformCosmicBackground")
175ws.spectral_radiance_surface_agendaSet(option="Blackbody")
176
177
178# In[ ]:
179
180
181def calculate_tbs_disort():
182    ws.disort_settings_agendaSetup(scattering_setting="ScatteringSpecies")
183    ws.disort_spectral_radiance_fieldProfile(
184        longitude=lon,
185        latitude=lat,
186        disort_quadrature_dimension=NQuad,
187        disort_legendre_polynomial_dimension=40,
188        disort_fourier_mode_dimension=1,
189        max_step=100
190    )
191    disort_stokes = [[ws.disort_spectral_radiance_field.data[f_ind, 0, 0, 0], 0.0, 0.0, 0.0] for f_ind in range(3)]
192    ws.spectral_radiance = disort_stokes
193    ws.spectral_radianceApplyForwardUnit(ray_path_point=ws.ray_path[0])
194    return ws.spectral_radiance.value.copy()[:, 0]
195
196
197
198# In[ ]:
199
200
201tbs_cloudy = calculate_tbs_disort()
202
203
204# ## Cloudy-sky brightness temperatures
205#
206# The cloudy-sky brightness temperature are obviously off. I checked the single-scattering albedo and extinction matrix between ARTS 2.6 and the new implementation and they agree. Currently, I am assuming there is an issue in DISORT.
207
208# ARTS 2.6 results: ``271.694859567588, 272.601957925916, 251.643215266136``
209
210# In[ ]:
211
212
213tbs_cloudy
214
215
216# ## Clear-sky brightness temperatures
217#
218# The clearsky brightness temperatures, however, agree well with the ARTS 2.6 results.
219
220# ARTS 2.6 results: ``298.566120236439 283.35611518369 251.643322551348``
221
222# In[ ]:
223
224
225ws.atmospheric_field[rain_first_moment] = 0.0
226tbs_clear = calculate_tbs_disort()
227tbs_clear