Boundary conditions and intensity reconstruction

Once the layer operators are known, the radiation field is fixed by the two external boundary conditions and reconstructed at every internal interface.

Top of atmosphere

The downward diffuse intensity at the top of the atmosphere is set to the user-supplied isotropic emission field:

\[I^{\downarrow}_{\mathrm{top},i} = B_{\mathrm{top}}, \qquad i = 1, \dots, M.\]

In thermal mode this is the Planck function at the topmost level, \(B(T_0)\), unless a separate top-boundary temperature is set (for example top_temperature = 0 for cold space, matching DisORT’s default of no downwelling). In non-thermal mode it is the raw value top_emission.

Lambertian surface

A Lambertian surface with albedo \(\alpha\) and thermal emission \(B_\mathrm{surf}\) is represented as an extra layer with

\[R_{ij} = 2\alpha\,\mu_j w_j\, c, \qquad s^{+}_i = (1-\alpha)\,B_\mathrm{surf},\]

where \(c = 1/(2\sum_j \mu_j w_j)\) is a quadrature-correction factor that enforces energy conservation. The solar surface source is

\[s^{+}_{i,\mathrm{sol}} = \frac{\alpha}{\pi}\,F_\odot\,\mu_0\, \ee^{-\tau^*/\mu_0},\]

with \(\tau^*\) the total (possibly delta-M scaled) optical depth of the atmosphere above the surface.

Diffusion lower boundary condition

For stellar-atmosphere applications where the deepest level is not a physical surface but an optically thick interior, a diffusion-approximation lower boundary can be used instead of a Lambertian surface. The upward intensity at the lower boundary is then set to

(4)\[I^{\uparrow}_{\mathrm{bot},i} = B(T_\mathrm{bot}) + \mu_i \left.\frac{\dd B}{\dd \tau}\right|_\mathrm{bot},\]

with the Planck gradient estimated from the bottom two levels,

\[\left.\frac{\dd B}{\dd \tau}\right|_\mathrm{bot} = \frac{B_\mathrm{bot} - B_\mathrm{bot-1}}{\Delta\tau_\mathrm{last}}.\]

When this boundary condition is active, no surface layer is added to the stack. It is enabled with use_diffusion_lower_bc and requires either use_thermal_emission or explicit planck_levels.

Intensity at internal interfaces

Composites are built from the top down (rtop) and from the bottom up (rbase). At each internal interface, the layer stack above it (rtop) and below it (rbase) are matched by solving for the upward and downward intensities:

\[\begin{split}\vect{I}^{\uparrow} &= (\mat{I} - \mat{R}_{ab,\mathrm{base}} \mat{R}_{ba,\mathrm{top}})^{-1} \bigl[\mat{T}_{ba,\mathrm{base}}\vect{I}_\mathrm{bot} + \mat{R}_{ab,\mathrm{base}}(\mat{T}_{ab,\mathrm{top}}\vect{I}_\mathrm{top} + \vect{s}^{-}_\mathrm{top}) + \vect{s}^{+}_\mathrm{base}\bigr], \\ \vect{I}^{\downarrow} &= (\mat{I} - \mat{R}_{ba,\mathrm{top}} \mat{R}_{ab,\mathrm{base}})^{-1} \bigl[\mat{T}_{ab,\mathrm{top}}\vect{I}_\mathrm{top} + \mat{R}_{ba,\mathrm{top}}(\mat{T}_{ba,\mathrm{base}}\vect{I}_\mathrm{bot} + \vect{s}^{+}_\mathrm{base}) + \vect{s}^{-}_\mathrm{top}\bigr],\end{split}\]

where the source vectors include both the thermal and the solar contributions. The matrices \((\mat{I} - \mat{R}\mat{R})\) are temperature-independent, so their LU factorisations can be cached and reused when computing the temperature Jacobian (Analytic temperature Jacobians).

Fluxes and mean intensity

The upward and downward diffuse hemispheric fluxes at each interface are

\[F^{\uparrow} = 2\pi \sum_{i=1}^{M} w_i\,\mu_i\, I^{\uparrow}_i, \qquad F^{\downarrow} = 2\pi \sum_{i=1}^{M} w_i\,\mu_i\, I^{\downarrow}_i.\]

The direct solar beam is tracked separately:

\[F_\mathrm{direct}(l) = F_\odot\,\mu_0\,\exp\!\Bigl(-\sum_{k=0}^{l-1} \tau^*_k / \mu_0\Bigr),\]

with \(\tau^*_k\) the (possibly delta-M scaled) optical depth of layer \(k\). The total net flux is \(F_\mathrm{net} = F^{\uparrow} - F^{\downarrow} - F_\mathrm{direct}\).

The actinic mean intensity \(J\) follows from the half-range quadrature,

(5)\[J = \frac{1}{2} \sum_{i=1}^{M} w_i\, \bigl(I^{\uparrow}_i + I^{\downarrow}_i\bigr) + \frac{F_\mathrm{direct}}{4\pi\mu_0}.\]

Following the DisORT convention, \(J\) is the full actinic mean intensity, including the collimated direct beam contribution \(F_\mathrm{direct}/(4\pi\mu_0)\). The diffuse-only mean intensity is recovered as mean_intensity - flux_direct / (4*pi*solar_mu).

Flux divergence (heating rate)

The net flux divergence at each level is

(6)\[\frac{\dd F}{\dd\tau} = 4\pi\,(1-\omega)\,(J - B),\]

where \(\omega\) is the unscaled single-scattering albedo. Following DisORT, each interface uses the albedo of the layer immediately above it (the top-of-atmosphere level uses the topmost layer). This quantity is the radiative heating rate per unit optical depth and is the basis of the heating-rate Jacobian (Analytic temperature Jacobians).