Pydisort Documentation

Python bindings for DISORT (Discrete Ordinates Radiative Transfer) Program.

Summary

This module provides a python interface to the C version of the DISORT program. Please consult the DISORT publication [1] for more information on the DISORT program, and the C-DISORT C publication [2] for more information on the C version of the DISORT program.

Small changes have been made to the C-DISORT program to make it compatible with python scripting. The C-DISORT program has been wrapped first in a C++ class (DisortWrapper), and the C++ class has been bound to python using pybind11.

Pydisort features the following benefits over the original C-DISORT program:

• Proper handling of errors rather than abrupt exit of the program. Errors can be caught and and handled in the python script.

• Memory management is handled by the C++ class. The user does not need to worry about memory allocation and deallocation.

• Documentation is automated using sphinx and readthedocs.

• Safety guards are implemented to prevent the user from setting incorrect values for arrays or calling methods in the wrong order.

Note that the underlying calculation engine is still the same as the C-DISORT program. So the speed of pydisort is the same as the origin C-DISORT program.

Examples

• Example 1: Calculate attenuation of radiative flux in a plane-parallel atmosphere

>>> from pydisort import disort, RFLDIR, FLDN, FLUP
>>> ds = disort()
>>> flags = {"onlyfl": True, "usrtau": False, "usrang": False}
>>> ds.set_flags(flags)
>>> ds.set_atmosphere_dimension(nlyr = 4)
>>> ds.seal()
>>> ds.set_optical_thickness([0.1, 0.2, 0.3, 0.4])
>>> ds.fbeam = 3.14159
>>> _, flx = ds.run()
>>> flx = flx[:, [RFLDIR, FLDN, FLUP]]
>>> flx
array([[3.14159   , 0.        , 0.        ],
       [2.84262818, 0.        , 0.        ],
       [2.32734711, 0.        , 0.        ],
       [1.72414115, 0.        , 0.        ],
       [1.15572637, 0.        , 0.        ]])

In the example above, flx has two dimensions. The first dimension is number of atmosphere levels (nlyr + 1 = 5), and the second dimension is number of extracted flux fields (3). RFDLIR, FLDN, and FLUP are the indices representing the three flux fields: direct flux, diffuse flux, and upward flux, respectively. Consult pydisort.run() method for more information on the indices of flux fields. The attenuation of radiative flux is acoording to the Beer-Lambert law, i.e.,

\[F(z) = F(0) \exp(-\tau(z)),\]

where \(F(z)\) is the radiative flux at level \(z\), \(F(0)\) is the radiative flux at the top of the atmosphere, and \(\tau(z)\) is the optical depth from the top of the atmosphere to level \(z\). The default direction of radiative flux is nadir.

• Example 2: Calculate intensity from isotropic scattering in a plane-parallel atmosphere

>>> from pydisort import disort, get_phase_function
>>> ds = disort()
>>> flags = {"planck": False}
>>> ds.set_flags(flags)
>>> ds.set_atmosphere_dimension(nlyr = 1, nstr = 16, nmom = 16)
>>> ds.set_intensity_dimension(nuphi = 1, nutau = 2, numu = 6)
>>> ds.seal()
>>> pmom = get_phase_function(nmom = 16, model = "isotropic")
>>> ds.set_optical_thickness([0.1])
>>> ds.set_single_scattering_albedo([1.0])
>>> ds.set_phase_moments(pmom)
>>> ds.set_user_optical_depth([0.0, 0.1])
>>> ds.set_user_cosine_polar_angle([-1.0, -0.5, -0.1, 0.1, 0.5, 1.0])
>>> ds.set_user_azimuthal_angle([0.0])
>>> ds.umu0 = 1.0
>>> ds.fbeam = 3.14159
>>> rad, flx = ds.run()
>>> rad
array([[[0.        , 0.        , 0.        , 0.18095504, 0.0516168 , 0.02707849],
        [0.02703935, 0.05146774, 0.17839685, 0.        , 0.        , 0.        ]]])

The intensity array rad has three dimensions. The first dimension is the azimuthal angles (1). The second dimension is optical depth (2). The third dimension is the polar angles (6). The result is interpreted as backsattering at the top of the atmosphere (optical depth = 0.0) and forward scattering at the bottom of the atmosphere (optical depth = 0.1).

Troubleshooting

• The most common error is “RuntimeError: DisortWrapper::Run failed”. When this error occurs, please check the error message printed before the error message. The error message printed before the error message usually provides more information on the cause of the error. Once you identify the cause of the error, you can fix the error by calling unseal() method, then setting the correct values, and then calling seal() method again.

• One common issue that results in “RuntimeError: DisortWrapper::Run failed” is incompatible flags for flux or intensity calculations. For example, usrtau and usrang flags should set to False when onlyfl flag is set to True.

• Another common issue is setting the wrong values for temperature, optical thickness, single scattering albedo, or phase function moments. All these values must be positive. If you set a negative value, you will get “RuntimeError: DisortWrapper::Run failed”.

• The program should not exit unexpectedly. If the program exits unexpectedly, please report the issue to the author (zoey.zyhu@gmail.com).

Important Dimensions

nlyr

number of atmosphere layers

nstr

number of radiation streams

nmom

number of phase function moments in addition to the zeroth moment

Tips

• Number of atmosphere levels is one more than the number of atmosphere layers.

• Temperature is defined on levels, not layers. Other properties such as optical thickness, single scattering albedo, and phase function moments are defined on layers.

• You can use print() method to print some of the DISORT internal states.

• You can chain methods such as set_flags, set_atmosphere_dimension(), set_intensity_dimension(), and seal() together.

• If you want to have more insights into DISORT internal inputs, you can set the print-input flag to True. The DISORT internal inputs will be printed to the standard output when the run() method is called.

References

[1]

Stamnes, K., Tsay, S. C., Wiscombe, W., & Jayaweera, K. (1988). Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media. Applied Optics, 27(12), 2502-2509.

[2]

Buras, R., & Dowling, T. (1996). Discrete-ordinate-method for radiative transfer in planetary atmospheres: Generalization of the doubling and adding method. Journal of Quantitative Spectroscopy and Radiative Transfer, 55(6), 761-779.

class pydisort.disort

Bases: pybind11_object

Wrapper class for the C-DISORT program

property albedo

Surface albedo, defaults to 0. Needed if lamber = True

property btemp

Bottom boundary temperature (K), defaults to 0

dimensions(self: pydisort.disort) → Tuple[int, int, int]

Get the dimensions of the DISORT solver

Parameters:

None –

Returns:

• nlyr (int) – Number of layers

• nstr (int) – Number of streams

• nmom (int) – Number of phase function moments in addition to the zeroth moment

property fbeam

Intensity of top-boundary incident parallel beam illumination, defaults to 0

property fisot

Intensity of top-boundary isotropic illumination, defaults to 0

property fluor

Intensity of bottom-boundary isotropic illumination, defaults to 0

property phi0

Azimuthal angle of incident beam, defaults to 0.

run(self: pydisort.disort) → Tuple[numpy.ndarray[numpy.float64], numpy.ndarray[numpy.float64]]

Run DISORT radiative transfer solver

Parameters:

None –

Returns:

• rad (numpy.ndarray, shape (nuphi, nutau, numu)) – radiation intensity

• flx (numpy.ndarray, shape (nlyr + 1, 8)) – radiation fluxes

Notes

The returned array references the internal memory of the disort object. The first dimension is the number of “levels”, which is one more than the number of layers. The second dimension is the number of flux fields. There are 8 flux fields in total, named as follows:

Index	Name	Description
0	RFLDIR	Direct beam flux
1	FLDN	Diffuse downward flux
2	FLUP	Diffuse upward flux
3	DFDT	Flux divergence, d (net flux) / d (optical depth)
4	UAVG	mean intensity including direct beam
5	UAVGDN	mean diffuse downward intensity
6	UAVGUP	mean diffuse upward intensity
7	UAVGSO	mean direct beam

All quantities are calculated without delta-M scaling.

Examples

>>> from pydisort import disort
>>> ds = pydisort.disort()
>>> ds.seal()
>>> rad, flx = ds.run()

In the example above, rad and flx references the internal memory of the DISORT solver without copying the data. If you want to extract a subset of the flux fields, you can do the following:

>>> from pydisort import disort, RFLDIR, FLDN, FLUP
>>> ds = pydisort.disort()
>>> ds.seal()
>>> _, flx = ds.run()
>>> flx = flx[:, [RFLDIR, FLDN, FLUP]]

However, using a subset of the flux fields will create a copy of the underlying memory. In the example above, flx is the variable that extracts the direct beam, diffuse downward, and diffuse upward fluxes.

seal(self: pydisort.disort) → None

Seal disort object after setting the dimensions

Return type:

None

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.seal()
>>> print(disort)

Notes

Internal memory is allocated after this function is called.

set_accuracy(self: pydisort.disort, arg0: float) → None

Set accuracy of disort convergence

Parameters:

arg0 (float) – Accuracy of disort convergence.

Return type:

None

Notes

Default accuracy is 1e-6

set_atmosphere_dimension(self: pydisort.disort, nlyr: int, nmom: int = 4, nstr: int = 4) → pydisort.disort

Set atmosphere dimension for disort

Parameters:

• nlyr (int) – Number of layers

• nmom (int, optional, defaults to 4) – Number of phase function moments

• nstr (int, optional, defaults to 4) – Number of streams (computational polar angles)

Return type:

DisortWrapper object

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.set_atmosphere_dimension(1, 4, 4, 4)
>>> print(disort)

Notes

It is common to set nstr = nmom. Memory was not allocated until disort.seal() is called.

See also

seal

allocate memory for disort

set_flags(self: pydisort.disort, arg0: Dict[str, bool]) → pydisort.disort

Set radiation flags for disort

Parameters:

arg0 (dict) – Dictionary of radiation flags consisting of a key (flag) and a value (True/False).

Return type:

DisortWrapper object

Examples

>>> import pydisort
>>> disort = pydisort.disort()
>>> dict = {'ibcnd': False, 'usrtau': True, 'usrang': True, 'lamber': True, 'plank': True}
>>> disort.set_flags(dict).seal()
>>> rad, flx = disort.run()

Notes

The following flags are supported:

Flag	Default value	Description
‘ibcnd’	False	General or Specific boundary condition
‘usrtau’	True	use user optical depths
‘usrang’	True	use user azimuthal angles
‘lamber’	True	turn on lambertian reflection surface
‘plank’	True	turn on plank source (thermal emission)
‘spher’	False	turn on spherical correction
‘onlyfl’	False	only compute radiative fluxes
‘quiet’	True	turn on disort internal printout
‘intensity_correction’	True	turn on intensity correction
‘old_intensity_correction’	True	turn on old intensity correction
‘general_source’	False	turn on general source
‘output_uum’	False	output azimuthal components of the intensity
‘print-input’	False	print input parameters
‘print-fluxes’	False	print fluxes
‘print-intensity’	False	print intensity
‘print-transmissivity’	False	print transmissivity
‘print-phase-function’	False	print phase function

A General boundary condition is invoked when ‘ibcnd’ is set to False. This allows:

• beam illumination from the top (set fbeam)

• isotropic illumination from the top (set fisot)

• thermal emission from the top (set ttemp and temis)

• interal thermal emission (use set_temperature_on_level)

• reflection at the bottom (set lamber, albedo)

• thermal emission from the bottom (set btemp)

A Special boundary condition is invoked when ‘ibcnd’ is set to True. Special boundary condition only returns albedo and transmissivity of the entire medium.

• current version of pydisort has limited support for this option.

• consult the documentation of DISORT for more details on this option.

set_header(self: pydisort.disort, arg0: str) → None

Set header for disort output

Parameters:

arg0 (str) – Header string

Return type:

None

set_intensity_dimension(self: pydisort.disort, nuphi: int = 1, nutau: int = 1, numu: int = 1) → pydisort.disort

Set intensity dimension for disort

Parameters:

• nuphi (int, optional, defaults to 1) – Number of user azimuthal angles

• nutau (int, optional, defaults to 1) – Number of user optical depths

• numu (int, optional, defaults to 1) – Number of user polar angles

Return type:

DisortWrapper object

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.set_intensity_dimension(1, 4, 4)

Notes

This function should be called when intensity calculations are requested. Memory was not allocated until disort.seal() is called.

See also

set_atmosphere_dimension

set atmosphere dimension

seal

allocate memory for disort

set_optical_thickness(self: pydisort.disort, arg0: List[float]) → None

Set layer optical thickness

Parameters:

arg0 (array of floats) – Optical thickness

Return type:

None

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.set_atmosphere_dimension(4)
>>> disort.seal()
>>> disort.set_optical_thickness([1.0, 2.0, 3.0, 4.0])

Warning

This function must be called after disort.set_atmosphere_dimension() and disort.seal().

Notes

If the size of the input array is inconsistent with the atmosphere dimension, the smaller size is used to set the internal array values.

See also

set_atmosphere_dimension

set atmosphere dimension

seal

allocate memory for disort

set_phase_moments(self: pydisort.disort, arg0: numpy.ndarray[numpy.float64]) → None

Set layer phase moments

Parameters:

arg0 (numpy.ndarray, shape (nlyr, nmom + 1) or (nmom + 1,)) – Phase function moments

Return type:

None

Example

>>> import pydisort
>>> ds = pydisort.disort()
>>> ds.set_atmosphere_dimension(nlyr = 1, nstr = 16, nmom = 16)
>>> ds.seal()
>>> _, _, nmom = ds.dimensions()
>>> pmom = get_phase_function(nmom, "isotropic")
>>> ds.set_phase_moments(pmom)

Warning

This function must be called after disort.set_atmosphere_dimension() and disort.seal().

Notes

If the size of the input array is inconsistent with the atmosphere dimension, the smaller size is used to set the internal array values.

See also

dimensions

get DISORT internal dimensions

get_phase_function

get phase function moments

set_single_scattering_albedo

set layer single scattering albedo

set_single_scattering_albedo(self: pydisort.disort, arg0: List[float]) → None

Set layer single scattering albedo

Parameters:

arg0 (array of floats) – Single scattering albedo

Return type:

None

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.set_atmosphere_dimension(4)
>>> disort.seal()
>>> disort.set_single_scattering_albedo([0.1, 0.2, 0.3, 0.4])

Warning

This function must be called after disort.set_atmosphere_dimension() and disort.seal().

Notes

If the size of the input array is inconsistent with the atmosphere dimension, the smaller size is used to set the internal array values.

See also

set_optical_thickness

set layer optical thickness

set_temperature_on_level(self: pydisort.disort, arg0: List[float]) → None

Set temperature (K) on levels (cell interfaces)

Parameters:

arg0 (array of floats) – Temperature on levels

Return type:

None

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.set_atmosphere_dimension(4)
>>> disort.seal()
>>> disort.set_temperature_on_level([100.0, 200.0, 300.0, 400.0, 500.0])

Notes

The number of levels is one more than the number of layers.

Warning

This function must be called after disort.set_atmosphere_dimension() and disort.seal().

Notes

If the size of the input array is inconsistent with the atmosphere dimension, the smaller size is used to set the internal array values.

See also

set_single_scattering_albedo

set layer single scattering albedo

set_user_azimuthal_angle(self: pydisort.disort, arg0: List[float]) → None

Set user azimuthal angle for disort

Parameters:

arg0 (array of floats) – User azimuthal angle in radians

Return type:

None

Warning

Must be called after disort.seal()

set_user_cosine_polar_angle(self: pydisort.disort, arg0: List[float]) → None

Set user cosine polar angle for disort

Parameters:

arg0 (array_like) – User cosine polar angle

Return type:

None

Warning

Must be called after disort.seal()

set_user_optical_depth(self: pydisort.disort, arg0: List[float]) → None

Set user optical depth for disort

Parameters:

tau (array_like) – User optical depth

Return type:

None

set_wavenumber_invcm(self: pydisort.disort, arg0: float) → None

Set a wavenumber for disort

Parameters:

arg0 (float) – Wavenumber

Return type:

None

Notes

This function sets both the minimum and maximum wavenumber to the same value. Internal thermal emission is calculated as the planck function at the specified wavenumber.

See also

set_wavenumber_range_invcm

set a wavenumber range for internal thermal emission

set_wavenumber_range_invcm(self: pydisort.disort, wmin: float, wmax: float) → None

Set a wavenumber range for disort

Parameters:

• wmin (float) – Minimum wavenumber

• wmax (float) – Maximum wavenumber

Return type:

None

Notes

This function sets the range of wavenumbers for internal thermal emission. Internal thermal emission is calculated as the integral of the planck function over the specified wavenumber range.

See also

set_wavenumber_invcm

set a wavenumber for internal thermal emission

property temis

Emissivity of top boundary. Needed if lamber = True, defaults to 0

property ttemp

Top boundary temperature (K), defaults to 0

property umu0

Cosine of incident beam zenith angle (positive), defaults to 1.

unseal(self: pydisort.disort) → None

Unseal disort object to allow resetting the dimensions

Return type:

None

Example

>>> import pydisort
>>> disort = pydisort.disort()
>>> disort.seal()
>>> disort.unseal()
>>> print(disort)

Notes

Internal memory is deallocated after unseal() is called.

pydisort.get_phase_function(nmom: int, model: str, gg: float = 0.0) → List[float]

Get phase function moments based on a phase function model

Parameters:

• nmom (int) – Number of phase function moments

• model (str) – Phase function model.

• gg (float) – Asymmetry factor for henyey-greenstein phase function

Returns:

pmom – Phase function moments, shape (1 + nmom,)

Return type:

List[float]

Notes

The following phase function models are supported:

Model	Description
‘isotropic’	Isotropic phase function, [1, 0, 0, 0, …]
‘rayleigh’	Rayleigh scattering phase function, [1, 0, 0.1, 0, …]
‘henyey_greenstein’	Henyey-Greenstein phase function, [1, gg, gg^2, gg^3, …]
‘haze_garcia_siewert’	Tabulated haze phase function by Garcia/Siewert
‘cloud_garcia_siewert’	Tabulated cloud phase function by Garcia/Siewert