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.

The normal usage of pydisort is to create a pydisort.DisortOptions object first and then initialize the pydisort.Disort object with the pydisort.DisortOptions object by:

>>> import pydisort
>>> op = pydisort.DisortOptions().flags("onlyfl,lamber")
>>> op.ds().nlyr = 4
>>> op.ds().nstr = 4
>>> op.ds().nmom = 4
>>> op.ds().nphase = 4
>>> ds = pydisort.Disort(op)

Examples

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

>>> import torch
>>> from pydisort import DisortOptions, Disort
>>> op = DisortOptions().flags("onlyfl,lamber")
>>> op.ds().nlyr = 4
>>> op.ds().nstr = 4
>>> op.ds().nmom = 4
>>> op.ds().nphase = 4
>>> ds = Disort(op)
>>> tau = torch.tensor([0.1, 0.2, 0.3, 0.4]).reshape((4,1))
>>> bc = {"fbeam" : torch.tensor([3.14159]).reshape((1,1))}
>>> flx = ds.forward(tau, bc)
>>> flx
tensor([[[[0.0000, 3.1416],
        [0.0000, 2.8426],
        [0.0000, 2.3273],
        [0.0000, 1.7241],
        [0.0000, 1.1557]]]])

It is important to understand the dimensions of the input and output arrays. The input array tau has two dimensions. In order of appearance, they are:

  1. The layer dimension (nlyr = 4),

  2. The property dimension (nprop = 1).

Since this problem only has optical thickness, the property dimension is 1. If not specified, both the wavelength/wavenumber dimension and the column dimension are assumed to be 1 and are automatically added internally to the input array.

The boundary condition dictionary bc has one key, fbeam, which is the solar beam flux. The key fbeam has two dimensions. In order of appearance, they are:

  1. The wavelength/wavenumber dimension (nwave = 1),

  2. The column dimension (ncol = 1).

In the example above, flx has four dimensions. In order of appearance, they are:

  1. The wavelenth/wavenumber dimension (nwave = 1),

  2. The column dimension (ncol = 1),

  3. The level dimension (nlvl = nlyr + 1 = 5),

  4. The flux field dimension (nflx = 2). The first element is upward flux, and the second element is downward flux.

The attenuation of radiative flux is according to the Beer-Lambert law, i.e., The example code above is in test_attenuation.py.

\[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.

Troubleshooting

  • The most common error is “RuntimeError: DisortImpl::forward”, which indicates that the disort run has failed. This error is mostly due to incorrect input dimensions or values. The error message shall provide more information on the cause of the error.

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

Tip

  • 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.

  • 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 forward() method is called.

References

Contents: