# Remarks on the Radiative Transfer Equations for Climatology

Claude Bardos111University of Paris Denis Diderot<EMAIL_ADDRESS>François Golse222Ecole Polytechnique, Palaiseau 91128, France
(francois.golse@polytechnique.edu). and Olivier Pironneau333Applied
Mathematics, Jacques-Louis Lions Lab, Sorbonne Université, 75252 Paris cedex
5, France<EMAIL_ADDRESS>https://www.ljll.math.upmc.fr/pironneau/).

###### Abstract

Using theoretical and numerical arguments we discuss some of the commonly
accepted approximations for the radiative transfer equations in climatology.

###### keywords:

Radiative transfer, Electronic Temperature , Integral equation, Numerical
analysis
AMS 3510, 35Q35, 35Q85, 80A21, 80M10

##  1 Introduction

Satellite, atmospheric and terrestrial measurements are numerous to support
the global warming. Various hypothesis are made to explain the disastrous
effect of the GreenHouse Gases (GHG). However a rigorous proof derived from
the fundamental equations of physics is not available and one of the problems
is the complexity of the physical system of planet Earth in its astronomical
setting around the Sun.

In many references such as kaper ,MOD ,CHA  and FOW  it is explained that the
Sun radiates light with a heat flux $Q=1370$Watt/m2, in the frequency range
$(0.5,20)\times 10^{14}$Hz corresponding approximately to a black body at
temperature of 5800K; 70% of this light intensity reaches the ground because
the atmosphere is almost transparent to this spectrum and about 30% is
reflected back by the clouds or the ocean, snow, etc (albedo). The Earth
behaves almost like a black body at temperature $T_{e}=288K$ and as such
radiates rays of frequencies $\nu$ in the infrared spectrum $(0.03,0.6)\times
10^{14}$Hz.

“The absorption coefficient variation with frequencies is such that the
atmosphere is essentially transparent to solar radiation” (FOW , p65).

“Solar radiation that is not absorbed or reflected by the atmosphere (for
example by clouds) reaches the surface of the Earth. The Earth absorbs most of
the energy reaching its surface, a small fraction is reflected444
https://public.wmo.int/en/sun%E2%80%99s-impact-earth.”

Carbone dioxide renders the atmosphere opaque to infrared radiations around
$20$THz . Hence increasing its proportion in air increases the absorption
coefficient in that range. It is believed to be one of the causes of global
warming. However it must be more complex because some experimental
measurements show that increasing opacity decreases the temperature at high
altitude, above the clouds DUF .

In the framework of the radiative transfer equations, one can studies the
conjectures cited above. To support theoretical claims, the numerical solver
developed in CBOP , OP  and FGOP  is used. We will address 4 questions:

  * •

What is the effect of increasing an altitude dependent absorption coefficient?

  * •

Is it true that sunlight crosses the Earth atmosphere unaffected?

  * •

Is 30% of albedo equivalent to 30% reduction of the Sun radiative energy ?

  * •

What is the effect on the temperature of increasing the absorption coefficient
in a specific frequency range?

##  2 The Radiative Transfer Equations for a Stratified Atmosphere

Finding the temperature $T$ in a fluid heated by electromagnetic radiations is
a complex problem because interactions of photons with atoms in the medium
involve rather intricate quantum phenomena. Assuming local thermodynamic
equilibrium leads to a well-defined electronic temperature. In that case, one
can write a kinetic equation for the radiative intensity
$I_{\nu}({\bm{x}},{\bm{\omega}},t)$ at time $t$, at position ${\bm{x}}$ and in
the direction $\omega$ for photons of frequency $\nu$, in terms of the
temperature field $T({\bm{x}},t)$.

$\displaystyle\frac{1}{c}\partial_{t}I_{\nu}+{\bm{\omega}}\cdot\nabla
I_{\nu}+\rho\bar{\kappa}_{\nu}a_{\nu}\left[I_{\nu}-{\tfrac{1}{4\pi}\int_{\SS^{2}}}p({\bm{\omega}},{\bm{\omega}}^{\prime})I_{\nu}({\bm{\omega}}^{\prime})\mathrm{d}\omega^{\prime}\right]$
(1) $\displaystyle=\rho\bar{\kappa}_{\nu}(1-a_{\nu})[B_{\nu}(T)-I_{\nu}]$ .

In this equation $c$ is the speed of light, $\rho$ is the density of the
fluid, $\nabla$ designates the gradient with respect to the position
${\bm{x}}$, while

$B_{\nu}(T)=\frac{2\hbar\nu^{3}}{c^{2}[{\rm e}^{\frac{\hbar\nu}{kT}}-1]}$ (2)

is the Planck function at temperature $T$, with $\hbar$ the Planck constant
and $k$ the Boltzmann constant. Recall the Stefan-Boltzmann identity,

$\int_{0}^{\infty}B_{\nu}(T)\mathrm{d}\nu=\bar{\sigma}T^{4}\,,\qquad\bar{\sigma}=\frac{2\pi^{4}k^{4}}{15c^{2}\hbar^{3}}\,,$
(3)

where $\pi\bar{\sigma}$ is the Stefan-Boltzmann constant.

The intricacy of the interaction of photons with the atoms of the medium is
contained in the mass-absorption $\bar{\kappa}_{\nu}$, which is theoretically
the result of vibration and rotation atomic resonance, but for practical
purpose the fraction of radiative intensity at frequency $\nu$ that is
absorbed by fluid per unit length. The coefficient $a_{\nu}\in(0,1)$ is the
scattering albedo, and
$\frac{1}{4\pi}p({\bm{\omega}},{\bm{\omega}}^{\prime})\mathrm{d}{\bm{\omega}}$
is the probability that an incident ray of light with direction
${\bm{\omega}}^{\prime}$ scatters in the infinitesimal element of solid angle
$\mathrm{d}{\bm{\omega}}$ centered at ${\bm{\omega}}$.

The kinetic equation (1) is coupled to the temperature - or energy
conservation - equations of the fluid . When thermal diffusion is small, the
system decouples and energy balance becomes:

$\int_{0}^{\infty}\rho\bar{\kappa}_{\nu}(1-a_{\nu})\left(\int_{\SS^{2}}I_{\nu}({\bm{\omega}})\mathrm{d}{\bm{\omega}}-4\pi
B_{\nu}(T)\right)\mathrm{d}\nu=0\,.$ (4)

When the wave source is far in the direction $z$ the problem becomes one-
dimensional in $z$ and the radiative intensity scattered in the direction
$\omega$ depends only on $\mu$, the cosine of the angle between $\omega$ and
$Oz$. Consequently the system becomes,

$\displaystyle\mu\partial_{z}I_{\nu}+\rho\bar{\kappa}_{\nu}I_{\nu}=\rho\bar{\kappa}_{\nu}(1-a_{\nu})B_{\nu}(T)+\frac{\rho\bar{\kappa}_{\nu}a_{\nu}}{2}\int_{-1}^{1}p(\mu,\mu^{\prime})I_{\nu}(z,\mu^{\prime})\mathrm{d}\mu^{\prime},$
(5)
$\displaystyle\int_{0}^{\infty}\rho\bar{\kappa}_{\nu}(1-a_{\nu})\left(B_{\nu}(T)-\tfrac{1}{2}\int_{-1}^{1}I_{\nu}\mathrm{d}\mu\right)\mathrm{d}\nu=0,\quad
z\in(0,H),~{}~{}|\mu|<1,~{}~{}\nu\in{\mathbb{R}}^{+}\,.$

For this system to be mathematically well posed, the radiation intensity
$I_{\nu}$ must be given on the domain boundary where radiation enters. For
example,

$~{}I_{\nu}(H,-\mu)=Q^{-}(\mu),~{}I(0,\mu)=Q^{+}(\mu),\quad 0<\mu<1\,.$ (6)

###  2.1 Semi-analytical Solution when the scattering is isotropic

Isotropic scattering is modelled by taking $p\equiv 1$. By introducing an
optical depth (in FOW  and others, the signs are changed),

$\tau=\int_{0}^{z}\rho(\eta)d\eta,\text{ $\rho$ becomes $1$ and $z$ becomes
$\tau\in(0,Z)$, $Z=\int_{0}^{H}\rho(\eta)d\eta$.}$ (7)

Let us denote the exponential integral by

$E_{p}(X):=\int_{0}^{1}e^{-X/\mu}\mu^{p-2}{d}\mu\,.$ (8)

Define

$\displaystyle J_{\nu}(\tau)=\frac{1}{2}\int_{-1}^{1}I_{\nu}d\mu,\quad
S_{\nu}(\tau)=\tfrac{1}{2}\int_{0}^{1}\left(e^{-\frac{\kappa_{\nu}\tau}{\mu}}Q^{+}_{\nu}(\mu)+e^{-\frac{\kappa_{\nu}(Z-\tau)}{\mu}}Q^{-}_{\nu}(\mu)\right){d}\mu.$
(9)

The problem is equivalent to a functional integral equation:

$\left.\begin{aligned}
&J_{\nu}(\tau)=S_{\nu}(\tau)+\frac{\kappa_{\nu}}{2}\int_{0}^{Z}E_{1}(\kappa_{\nu}|\tau-t|)\left((1-a_{\nu})B_{\nu}(T(t))+a_{\nu}J_{\nu}(t)\right)dt\,,\\\
&\int_{0}^{\infty}(1-a_{\nu})\kappa_{\nu}B_{\nu}(T(\tau)){d}\nu=\int_{0}^{\infty}(1-a_{\nu})\kappa_{\nu}J_{\nu}(\tau){d}\nu,\end{aligned}\right\\}$
(10)

Once $J_{\nu},T$ are known, $I_{\nu}$ can be recovered by

$\displaystyle I_{\nu}(\tau,\mu)=$ $\displaystyle
e^{-\frac{\kappa_{\nu}\tau}{\mu}}I_{\nu}(0,\mu)\mathbf{1}_{\mu>0}+e^{\frac{\kappa_{\nu}(Z-\tau)}{\mu}}I_{\nu}(Z,\mu)\mathbf{1}_{\mu<0}$
(11) $\displaystyle+{\bf
1}_{\mu>0}\int_{0}^{\tau}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt$
$\displaystyle+{\bf
1}_{\mu<0}\int_{\tau}^{Z}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt.$

###### Remark 1

Notice that two solutions $(T,J_{\nu}),(T^{\prime},J^{\prime}_{\nu})$ with
$J_{\nu}\equiv J^{\prime}_{\nu}$ will have $T\equiv T^{\prime}$ but not
necessarily $I_{\nu}\equiv I_{\nu}^{\prime}$.

##  3 Increasing an Altitude Dependent Absorption Coeficient Decreases the
Temperature in the atmosphere

The following argument goes against a common belief which claims that
increasing the absorption coefficient - which GHG do indeed - implies an
increase of temperature in the atmosphere. Consider the case of an altitude
dependent absorption coefficient which we write as $r(z)\kappa$ with $\kappa$
constant.

###### Proposition 1

If $\kappa$ is independent of $\nu$ but function of the altitude and $T$ is
monotone decreasing with altitude, then, increasing $\kappa$ anywhere will
decrease the temperature.

Proof The proof is a straightforward consequence of the concept of optical
length (7). Here, let $\tau=\int_{0}^{z}{r}(\zeta)d\zeta.$ Observe that
$\partial_{z}=\frac{\partial\tau}{\partial\tau}\partial_{\tau}={r}(z)\partial_{z}$;
consequently ${r}$ disappears from (5), $H$ is replaced by
$Z=\int_{0}^{H}{r}(\zeta)d\zeta$, and $z$ is replaced by $\tau$.

When $\kappa$ and $a$ are constant,
$\displaystyle\bar{I}_{r}=\int_{0}^{\infty}I_{\nu}d\nu$ is solution of the
”grey model”:

$\displaystyle\mu\partial_{\tau}\bar{I}_{1}+\kappa\bar{I}_{1}=\frac{1}{2}\int_{-1}^{1}\left(\kappa(1-a)\bar{I}_{1}+\kappa
a\int_{-1}^{1}p(\mu,\mu^{\prime}){\bar{I}_{1}}(\mu,\mu^{\prime})d\mu^{\prime}\right)d\mu,$
(12) $\displaystyle I_{r}(H,-\mu)=0,\quad\bar{I}_{r}(0,\mu)=Q^{+}(\mu),\quad
0<\mu<1.$

where the Stefan-Boltzmann identity has been used: $\displaystyle\sigma
T^{4}=\tfrac{1}{2}\int_{-1}^{1}{\bar{I}_{r}}\mathrm{d}\mu.$

Let $\bar{I}_{1}(\tau,\mu)$ be the solution of (12). We want to compare
$I_{{r}+\delta{r}}$ with $I_{{r}}$ when $\delta{r}$ has a small support near
altitude $y$. By (7),

$\displaystyle\bar{I}_{{r}+\delta{r}}(z)=\bar{I}_{1}(\int_{0}^{z}({r}(\zeta)+\delta{r}(\zeta))d\zeta)\approx\bar{I}_{1}(\tau(z))+\frac{d}{d\tau}\bar{I}_{1}(\tau(y))\delta{r}(y)$
(13)
$\displaystyle\implies\qquad\frac{1}{2}\int_{-1}^{1}\bar{I}_{{r}+\delta{r}}(z)d\mu-\frac{1}{2}\int_{-1}^{1}\bar{I}_{r}(z)d\mu\approx\frac{d}{d\tau}\left(\frac{1}{2}\int_{-1}^{1}\
\bar{I}_{1}(\tau)d\mu\right)|_{\tau(y)}\delta{r}(y)$

Finally from (3),

$T_{{r}+\delta{r}}(z)-T_{{r}}(z)\approx 4\sigma
T^{3}(y)\frac{\delta{r}(y)}{{r}(y)}\frac{dT_{r}}{dz}(y).$ (14)

$\square$

##  4 Sunlight Crosses the Earth Atmosphere Unaffected?

It is true that the absorption coefficient $\kappa_{\nu}$ is small in the
frequency range of visible light; but does it imply that the light source can
be used as a boundary condition at altitude zero instead of altitude H, the
top of the troposphere?

First eliminate the $z$ dependency of the coefficients by using the optical
length $\tau$ defined in (7). The following result gives an interesting answer
to this question.

###### Proposition 2

Assume $a_{\nu}=0$ and assume that for some $\nu^{*}$,

$\epsilon:=\max_{\nu>\nu^{*}}\kappa_{\nu}<<1\text{ and
}\quad\kappa_{\nu}Q^{0}_{\nu}|_{\nu>\nu^{*}}=O(1).$

Let $I_{\nu},T$ and $I^{\prime}_{\nu},T^{\prime}$ be solutions of (5),(6)
respectively with

$\displaystyle Q^{+}_{\nu}(\mu)=0,\quad Q^{-}_{\nu}(\mu)=\mu Q^{0}_{\nu},\quad
0<\mu<1.$ (15) $\displaystyle{Q^{\prime}}^{+}_{\nu}=\mu{\rm
e}^{-\frac{\kappa_{\nu}Z}{\mu}}Q^{0}_{\nu},\quad{Q^{\prime}}^{-}_{\nu}(\mu)=0,\quad
0<\mu<1.$

Then $T^{\prime}(\tau)=T(\tau)+O(\epsilon)$ and

$I^{\prime}_{\nu}-I_{\nu}=\mu{\rm
e}^{-\frac{\kappa_{\nu}Z}{\mu}}Q^{0}_{\nu}\left({\bf 1}_{\mu>0}{\rm
e}^{-\frac{\kappa_{\nu}}{\mu}\tau}+{\bf 1}_{\mu<0}{\rm
e}^{\frac{\kappa_{\nu}}{\mu}\tau}\right)+O(\epsilon).$

Proof Note that

$\left.\begin{aligned}
&\mu\partial_{\tau}I_{\nu}+\kappa_{\nu}(I_{\nu}-B_{\nu}(T))=0,\quad\nu<\nu^{*},\\\
&I_{\nu}(0)|_{\mu>0}=0,\quad
I_{\nu}(Z)_{\mu<0}=Q^{-}_{\nu}(-\mu),\quad\quad\nu<\nu^{*},\\\ &I_{\nu}\approx
I_{\nu}^{*}(\mu)\text{ independent of $\tau$},~{}\quad\nu>\nu^{*}.\\\
&I^{*}_{\nu}|_{\mu>0}=0,\quad
I^{*}_{\nu}|_{\mu<0}=Q^{-}_{\nu}(-\mu),\quad\quad\nu>\nu^{*}.\\\
&\int_{0}^{\infty}\kappa_{\nu}B_{\nu}(T)d\nu\approx\int_{0}^{\nu^{*}}\frac{\kappa_{\nu}}{2}\int_{-1}^{1}I_{\nu}\,d\mu
d\nu+\int_{\nu^{*}}^{\infty}\frac{\kappa_{\nu}}{2}\int_{-1}^{0}Q^{-}_{\nu}(-\mu)\,d\mu
d\nu\end{aligned}\right\\}$ (16)

Now $I_{\nu}^{\prime}$ is defined by (5) but with

$I^{\prime}_{\nu}(0)|_{\mu>0}={Q^{\prime}}_{\nu}^{+}(\mu),\qquad
I^{\prime}_{\nu}(Z)|_{\mu<0}=0.$

By the same argument, when $\nu>\nu^{*}$,
$I^{\prime}_{\nu}\approx{I^{\prime}}^{*}_{\nu}(\mu)$, independent of $\tau$,
and

$\left.\begin{aligned}
&\mu\partial_{\tau}I^{\prime}_{\nu}+\kappa_{\nu}(I^{\prime}_{\nu}-B_{\nu}(T^{\prime}))=0,\quad\nu<\nu^{*},\\\
&I^{\prime}_{\nu}(0)|_{\mu>0}={Q^{\prime}}_{\nu}^{+}(\mu),\qquad
I^{\prime}_{\nu}(Z)|_{\mu<0}=0,\quad\quad\nu<\nu^{*},\\\
&\int_{0}^{\infty}\kappa_{\nu}B_{\nu}(T^{\prime})d\nu\approx\int_{0}^{\nu^{*}}\frac{\kappa_{\nu}}{2}\int_{-1}^{1}I^{\prime}_{\nu}d\mu
d\nu+\int_{\nu^{*}}^{\infty}\frac{\kappa_{\nu}}{2}\int_{0}^{1}{Q^{\prime}}^{+}_{\nu}(\mu)d\mu
d\nu\end{aligned}\right\\}$ (17)

Let $I^{\prime\prime}_{\nu}=I_{\nu}-I^{\prime}_{\nu}$. It holds

$\left.\begin{aligned}
&\mu\partial_{\tau}I^{\prime\prime}_{\nu}+\kappa_{\nu}(I^{\prime\prime}_{\nu}-(B_{\nu}(T)-B_{\nu}(T^{\prime}))=0,\quad\quad\nu<\nu^{*},\\\
&\int_{0}^{\infty}\kappa_{\nu}(B_{\nu}(T)-B_{\nu}(T^{\prime}))d\nu=\int_{0}^{\nu^{*}}\frac{\kappa_{\nu}}{2}\int_{-1}^{1}I^{\prime\prime}_{\nu}\,d\mu
d\nu\\\
&+\int_{\nu^{*}}^{\infty}\left(\frac{\kappa_{\nu}}{2}\int_{0}^{1}(Q^{-}_{\nu}(\mu)-{Q^{\prime}}^{+}_{\nu}(\mu))\,d\mu\right)d\nu\\\
&I^{\prime\prime}_{\nu}(0)|_{\mu>0}=-{Q^{\prime}}_{\nu}^{+}(\mu),\qquad
I^{\prime\prime}_{\nu}(Z)|_{\mu<0}=Q_{\nu}^{-}(-\mu),\quad\quad\nu<\nu^{*}.\end{aligned}\right\\}$
(18)

We notice that $T=T^{\prime}$ and $Q_{\nu}^{-}(\mu)={\rm
e}^{-\frac{\kappa_{\nu}}{\mu}Z}{Q^{\prime}}^{+}_{\nu}(\mu)$ implies

$\left.\begin{aligned}
&\mu\partial_{\tau}I^{\prime\prime}_{\nu}+\kappa_{\nu}I^{\prime\prime}_{\nu}=0,\quad\nu<\nu^{*},\\\
&\int_{0}^{\nu^{*}}\left(\frac{\kappa_{\nu}}{2}\int_{-1}^{1}I^{\prime\prime}_{\nu}\,d\mu\right)d\nu=\int_{\nu^{*}}^{\infty}\left(\frac{\kappa_{\nu}}{2}\int_{0}^{1}(1-{\rm
e}^{-\frac{\kappa_{\nu}}{\mu}Z}){Q}^{+}_{\nu}(\mu)\,d\mu\right)d\nu\\\
&I^{\prime\prime}_{\nu}(0)|_{\mu>0}=-{Q^{\prime}}_{\nu}^{+}(\mu),\qquad
I^{\prime\prime}_{\nu}(Z)|_{\mu<0}={\rm
e}^{\frac{\kappa_{\nu}}{\mu}Z}{Q^{\prime}}_{\nu}^{+}(-\mu)\end{aligned}\right\\}$
(19)

Let ${Q^{\prime}}_{\nu}^{+}(\mu)=\mu Q^{0}_{\nu}$. Then

$I^{\prime\prime}_{\nu}(0)|_{\mu>0}=-\mu Q^{0}_{\nu},\qquad
I^{\prime\prime}_{\nu}(Z)|_{\mu<0}=-\mu{\rm
e}^{\frac{\kappa_{\nu}}{\mu}Z}Q^{0}_{\nu}.$

Consequently, from (LABEL:eqI6) and (11), the solution is

$\displaystyle I^{\prime\prime}_{\nu}(\tau,\mu)=-\mu Q^{0}_{\nu}\left({\bf
1}_{\mu>0}{\rm e}^{-\frac{\kappa_{\nu}}{\mu}\tau}+{\bf 1}_{\mu<0}{\rm
e}^{\frac{\kappa_{\nu}}{\mu}\tau}\right).$ (20) $\displaystyle\implies$
$\displaystyle\int_{-1}^{1}I^{\prime\prime}_{\nu}(\tau)d\mu=-Q^{0}_{\nu}E_{3}(\kappa_{\nu}\tau)+Q^{0}_{\nu}E_{3}(\kappa_{\nu}\tau)=0,$

We have satisfied (LABEL:eqI6) with precision $O(\kappa_{\nu}^{2})$ because

$\kappa_{\nu}(1-{\rm e}^{-\frac{\kappa_{\nu}}{\mu}Z})\mu
Q^{0}_{\nu}=\kappa_{\nu}^{2}ZQ^{0}_{\nu}+O(\kappa_{\nu}^{3})\approx
0,\quad\nu>\nu^{*}$

$\square$

###### Remark 2

It is easy to see that applying ${\rm e}^{-\frac{\kappa_{\nu}Z}{\tau}}\mu
Q^{0}_{\nu}$ at $\tau=0$, leads to change $E_{3}(\kappa_{\nu}\tau)$ in (11)
into $E_{3}(\kappa_{\nu}(Z+\tau))$.

The change in temperature is shown on Figure 1. All variables are rescaled as
in FGOPSIAM $\S 7.1$. In particular temperatures are in Kelvin and divided by
4780. The air density is $1.22{\rm e}^{-z}$. Here $Q^{0}=5.66\cdot 10^{-5}$
and $T_{sun}=1.2$.

_Conclusion: For the specific case of sunlight on Earth It is a gross
approximation to assume that sunlight crosses the atmosphere unaffected._

However the theory works when $T_{sun}=2.4$ and $\nu^{*}=0.2$, probably
because $\max_{\nu}B_{\nu}(T_{sun}\approx 18$ against 2 with $T_{sun}=1.2$.

$0$$5$$10$$15$$5$$6$$7$$\cdot 10^{-2}$Altitude - kmScaled Temperature
(K)$T_{0}$$T_{1}$$T_{2}$$T_{3}$$T_{4}$$T_{5}$ Figure 1:  Scaled temperatures
$z\to T(z)$ versus altitude, computed with $0.5(0.1+{\bf 1_{\nu<3}})$. The 3
curves $T_{1},T_{2},T_{3}$ corresponds to 1: $Q^{-}_{\nu}=Q^{0}_{\nu}$,
$Q^{+}_{\nu}=0$, 2: $Q^{-}_{\nu}=0$, $I_{\nu}(0,\mu)|_{\mu>0}=\mu{\rm
e}^{-\frac{\kappa_{\nu}Z}{\mu}}Q^{0}_{\nu}$ and to 3: $Q^{-}_{\nu}=0$,
$Q^{+}_{\nu}=Q^{0}$. These can be compared to the grey case $T_{0}$ with
$\kappa=0.5$. Note that $T_{1}$ is very different from $T_{2}$. However if the
Sun temperature is multiplied by 2 and $\kappa_{\nu}=0.5(0.1+{\bf
1_{\nu<0.2}})$, then the conditions in Proposition 2 are met and confirmed by
the numerical simulations: $T_{4}$ is computed like $T_{1}$ and $T_{5}$ is
computed like $T_{2}$.

##  5 Is the 30% Earth Albedo Equivalent to a 30% reduction of the Sun
Radiative Energy?

Note that (11) can be used to find the solution of (5)-a with generalized
albedo (inspired by the model proposed in DUF ):

$I_{\nu}(0,\mu)=\sum_{k}\alpha_{k}I_{\nu}(\tau_{k},-\mu)+Q^{0}_{\nu}\mu,\quad
I_{\nu}(Z,-\mu)=0,\quad 0<\mu<1.$

with $k=1,2,..K$. Indeed by (11)

$\displaystyle I_{\nu}(\tau,\mu)=$ $\displaystyle
e^{-\frac{\kappa_{\nu}\tau}{\mu}}\left(\sum_{k}\alpha_{k}I_{\nu}(\tau_{k},-\mu)+Q^{0}_{\nu}\mu\right)\mathbf{1}_{\mu>0}$
(21) $\displaystyle+{\bf
1}_{\mu>0}\int_{0}^{\tau}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt$
$\displaystyle+{\bf
1}_{\mu<0}\int_{\tau}^{Z}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt$
$\displaystyle=\mathbf{1}_{\mu>0}e^{-\frac{\kappa_{\nu}\tau}{\mu}}\sum_{k}\alpha_{k}\int_{\tau_{k}}^{Z}e^{-\frac{\kappa_{\nu}(t-\tau_{k})}{\mu}}\tfrac{\kappa_{\nu}}{\mu}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt$
(22) $\displaystyle+{\bf
1}_{\mu>0}\left(e^{-\frac{\kappa_{\nu}\tau}{\mu}}Q^{0}_{\nu}\mu+\int_{0}^{\tau}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt\right)$
$\displaystyle+{\bf
1}_{\mu<0}\int_{\tau}^{Z}e^{-\frac{\kappa_{\nu}|\tau-t|}{|\mu|}}\tfrac{\kappa_{\nu}}{|\mu|}\left((1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right)dt$
$\displaystyle\implies$ $\displaystyle
J_{\nu}(\tau)=\tfrac{1}{2}Q^{0}_{\nu}E_{3}(\kappa_{\nu}\tau)+\tfrac{1}{2}\int_{0}^{Z}E_{1}(\kappa_{\nu}|\tau-t|)\left[(1-a_{\nu})B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right]dt\,,$
$\displaystyle+\tfrac{1}{2}\sum_{k}\alpha_{k}\kappa_{\nu}\int_{\tau_{k}}^{Z}E_{1}(\kappa_{\nu}(t+\tau-\tau_{k}))\left[(1-a_{\nu})\kappa_{\nu}B_{\nu}(T(t))+a_{\nu}\kappa_{\nu}J_{\nu}(t)\right]dt.$

###  5.1 Earth Albedo

Consider the case

$I_{\nu}(Z,\mu)|_{\mu<0}=0,\quad I_{\nu}(0,\mu)|_{\mu>0}=Q^{0}_{\nu}+\alpha
I_{\nu}(0,-\mu)\,.$

Assume no scattering, $\tau_{1}=0$, $K=1$. As in FGOP , the following
iterative scheme is considered:

$\displaystyle J_{\nu}^{n+1}(\tau)$
$\displaystyle=\tfrac{1}{2}Q^{0}_{\nu}E_{3}(\kappa_{\nu}\tau)+\frac{\kappa_{\nu}}{2}\int_{0}^{Z}\left(E_{1}(\kappa_{\nu}|\tau-t|)+\alpha
E_{1}(\kappa_{\nu}(t+\tau))\right)B_{\nu}(T^{n}(t))dt\,,$ (23)
$\displaystyle\int_{0}^{\infty}\kappa_{\nu}B_{\nu}({T^{n+1}_{\tau}})\mathrm{d}\nu=\int_{0}^{\infty}\kappa_{\nu}J^{n+1}_{\nu}(\tau)\mathrm{d}\nu\,.$

###  5.2 Numerical Example

Assume $a_{\nu}=0$. Figure 2 compares the numerical solutions of (5) with
$\kappa_{\nu}=0.5[{\bf 1}_{\nu<6}+0.1]$ and

$Q^{-}_{\nu}=0,\quad Q^{+}_{\nu}=\lambda Q^{0}_{\nu},\quad\lambda=1\text{ or
}0.7.$

with the solution of (5) with $\beta=0$ or 1 in

${{Q^{-}_{\nu}}}^{\prime}=0,\quad{{I_{\nu}}}^{\prime}(0,\mu)=0.3{{I_{\nu}}}^{\prime}(0,-\mu)+\mu
Q^{0}_{\nu}[\tfrac{1}{0.7}{\rm
e}^{-\frac{\kappa_{\nu}Z}{\mu}}\beta+1-\beta],\quad 0<\mu<1.$

$0$$2$$4$$6$$8$$10$$12$$6$$6.5$$7$$7.5$$\cdot 10^{-2}$Altitude - kmScaled
Temperature (K)$T,Q^{0}$$T,0.7Q^{0}$$T^{\prime}$$T"$ Figure 2:  Scaled
temperatures $z\to T(z)$ versus altitude. In all cases $Q^{-}_{\nu}=0$. The
doted curves are $T$ when $Q^{+}_{\nu}=\lambda Q^{0}_{\nu}$ with $\lambda=1$
(top dotted curve) and $\lambda=0.7$ (lower dotted curve); the doted line is
the solution $T^{\prime}$ when
$I_{\nu}(0,\mu)|_{\mu>0}=0.3I_{\nu}(0,-\mu)+\mu{\rm
e}^{-\frac{\kappa_{\nu}Z}{\mu}}Q^{0}_{\nu}/0.7$ . The dashed line is when
$I_{\nu}(0,\mu)|_{\mu>0}=0.3I_{\nu}(0,-\mu)+\mu Q^{0}_{\nu}$.

This numerical simulation is an attempt to replace the radiative source at
$\tau=H$ by a radiative source at $\tau=0$ which gives the same result.

Notice alongside that adding 0.3 albedo induces a heating of the atmosphere
(i.e. comparing both computations with $\lambda=1,\beta=1$ above).

##  6 Earth Albedo, General Statement, Accommodation Coefficient

As above one considers a stratified atmosphere (5) with $(z,\mu,\nu)$ in
$(0,H)\times(-1,1)\times{\mathbb{R}}^{+}$; the radiation intensity
$I_{\nu}(z,\mu)$ satisfies on at $z=H$ an incoming condition

$I_{\nu}(H,-\mu)=Q_{\nu}^{-}(\mu),\quad 0<\mu<1,$ (24)

and on Earth, at $z=0$, some “albedo” condition, denoted:

$\displaystyle I_{\nu,}(0,\mu)=\mathcal{A}_{\nu}(I_{\nu}(0,-\mu),\quad
0<\mu<1,$ (25)

with $\mathcal{A}$ representing the albedo effect of the earth, hence being an
operator from the space of outgoing intensities into the space of incoming
intensities.

Below will be given, based on energy estimates, some sufficient conditions on
$\mathcal{A}$ for existence, stability – hence uniqueness – of the solution,
beginning with the grey model and then considering a natural operator which
involves the effective temperature of Earth (see FOW  page 66).

###  6.1 The energy type estimate for the grey problem

For clarity we present the case without scattering ($a_{\nu}=0$) and with a
constant $\kappa_{\nu}$ independent of the frequency $\nu$. By rescaling $z$
with $\kappa$, the problem can be formulated in term of
$\displaystyle\overline{I}(z,\mu)=\int_{0}^{\infty}I_{\nu}(z,\mu,)\mathrm{d}\nu$
as,

$\mu\partial_{z}\overline{I}(z,\mu)+\overline{I}(z,\mu)-\frac{1}{2}\int_{-1}^{1}\overline{I}(z,\mu^{\prime})\mathrm{d}\mu^{\prime}=0\,,$
(26a)
$\overline{I}(H,-\mu)=Q(\mu):=\overline{Q}^{-}_{\nu}(\mu),\quad\overline{I}(0,\mu)=\mathcal{A}({\overline{I}}(0,-\mu)),\quad
0<\mu<1\,.$ (26b)

Note that we have assumed that ${\mathcal{A}}$ commutes with the integration
in $\nu$. First observe that the operator
$\displaystyle{\overline{I}}\mapsto{\overline{I}}-\frac{1}{2}\int_{-1}^{1}{\overline{I}}(\mu^{\prime})\mathrm{d}\mu^{\prime}$
is, in the space $L^{2}(-1,1)$ the orthogonal projection on functions of mean
value $0$, because

$\int_{-1}^{1}({\overline{I}}-\frac{1}{2}\int_{-1}^{1}{\overline{I}}(\mu^{\prime})\mathrm{d}\mu^{\prime}){\overline{I}}(\mu)d\mu=\int_{-1}^{1}({\overline{I}}-\frac{1}{2}\int_{-1}^{1}{\overline{I}}(\mu^{\prime})\mathrm{d}\mu^{\prime})^{2}\mathrm{d}\mu\,.$
(27)

This implies that the left hand side of (27) is nonnegative and equal to $0$
iff
$\displaystyle{\overline{I}}=\frac{1}{2}\int_{-1}^{1}{\overline{I}}(\mu^{\prime})\mathrm{d}\mu^{\prime}$
is independent of $\mu\,.$ Multiplying the equation (26a) by
${\overline{I}}(z,\mu)\,,$ integrating over $(0,H)\times(-1,1)\,,$ using
Green’s formula and inserting in the computations the conditions (26b) one
obtains:

$\displaystyle\int_{0}^{H}\mathrm{d}z\int_{-1}^{1}({\overline{I}}-\frac{1}{2}\int_{-1}^{1}{\overline{I}}(\mu^{\prime})\mathrm{d}\mu^{\prime})^{2}\mathrm{d}\mu+\frac{1}{2}\int_{0}^{1}\mu({\overline{I}}(H,-\mu)^{2}\mathrm{d}\mu$
(28)
$\displaystyle+\frac{1}{2}\int_{0}^{1}\mu\big{(}({\overline{I}}(0,\mu))^{2}-(\mathcal{A}({\overline{I}}(0,\mu))^{2}\big{)}\mathrm{d}\mu=\frac{1}{2}\int_{0}^{1}\mu(Q(\mu))^{2}\mathrm{d}\mu$

This estimate indicates that any property which implies the positivity of

$\int_{0}^{1}\mu\big{(}{\overline{I}}(0,\mu)^{2}-(\mathcal{A}({\overline{I}}(\mu,0))^{2}\big{)}\mathrm{d}\mu$

will lead to a well posed problem with a unique nonnegaative solution; for
instance:

###### Theorem 1

Assume that the operator ${\overline{I}}\mapsto\mathcal{A}({\overline{I}})$
restricted to the positive cone $C^{+}({\overline{I}}(\mu)\geq 0)$ on the
space $(L^{2}(-1,1),\mu,d\mu)$ then for any positive $\mu\mapsto Q(\mu)$ there
is a unique nonnegative solution of the (26a),(26b).

###### Remark 3

Making use of the maximum principle one could prove similar results when the
albedo operator is a contraction in the positive cone of $L^{\infty}(0,1)\,.$

###### Remark 4

The albedo condition which is a relation, for $0<\mu<1\,,$ between
${\overline{I}}(0,\mu)\big{|}_{\mu>0}$ and ${\overline{I}}(0,\mu)|_{\mu<0}$
can also be expressed in term of ${\overline{I}}(0,\mu)|_{\mu>0}$ only.

Indeed, define the operator
${\mathcal{T}}:{\overline{I}}_{+}(\mu)\mapsto{\mathcal{T}}({\overline{I}}_{+}(\mu))=I(0,-\mu)|_{\mu>0}$
(defined for instance in $L^{\infty}(0,1))$ by solving (26a) with

${\overline{I}}(H,-\mu)=Q(\mu),\quad{\overline{I}}(0,\mu)={\overline{I}}_{+}(\mu),\quad
0<\mu<1\,.$ (29)

Then one has:

${\overline{I}}(0,\mu)|_{\mu>0}=\mathcal{A}({\overline{I}}(0,-\mu))\Leftrightarrow{\overline{I}}(0,\mu)|_{\mu>0}=\mathcal{A}(\mathcal{T}({\overline{I}}(0,\mu)|_{\mu>0}))\,.$
(30)

With $Q(\mu)>0$ in $L^{\infty}(0,1)$, $\mathcal{T}$ is an affine operator
which preserves the positivity and the monotonicity. and this leads to the
simple but useful proposition, which , in turn, is the key to the proof of
Theorem 1.

###### Proposition 3

Consider two solutions $\\{I_{i}\\}_{i=1,2}$ of the (26a) with (26b) with the
same $Q(\mu)$ but with two different albedo operators $\mathcal{A}_{i}\,.$
Assume that both $\mathcal{A}_{i}$ are linear contractions and one of them is
a strict linear contraction, which preserve positivity. If

$\forall f\geq 0\in L^{\infty}(0,1)\quad\forall\mu\in(0,1)\,,\quad
A_{2}(f)(\mu)\leq A_{1}(f)(\mu)$ (31)

then one has the same ordering for the corresponding solutions :

$\forall(z,\mu)\in(0,1)\times(0,H)\quad I_{1}(z,\mu)\leq I_{2}(z,\mu)\,.$ (32)

Proof Using the linearity of the operators $\mathcal{A}_{i}$ one observes that
$R=I_{1}-I_{2}$ solves(26a) with $R(H,-\mu)=0$ and, for all $0<\mu<1$,

$R(0,\mu)-\frac{1}{2}(\mathcal{A}_{1}+\mathcal{A}_{2})(\mathcal{T}(R(0,\mu))=\frac{1}{2}(\mathcal{A}_{1}-\mathcal{A}_{2})(\mathcal{T}((I_{1}(0,\mu)+I_{2}(0,\mu)))\,.$
(33)

Moreover, restricted to solutions which satisfy $R(H,-\mu)=0$, $\mu>0$ the
operator $\mathcal{T}$ is a strict monotonicity preserving linear contraction;
the same observation holds for
$R\mapsto\frac{1}{2}(\mathcal{A}_{1}+\mathcal{A}_{2})(\mathcal{T}(R)$. As a
consequence the operator
$R\mapsto(I-\frac{1}{2}(\mathcal{A}_{1}+\mathcal{A}_{2})(\mathcal{T}))$ is
invertible with inverse given by the Neumann series (also preserving the
positivity):

$R=\sum_{k\geq
0}(\frac{1}{2}(\mathcal{A}_{1}+\mathcal{A}_{2})(\mathcal{T})))^{k}\frac{1}{2}(\mathcal{A}_{1}-\mathcal{A}_{2})(\mathcal{T}((I_{1}+I_{2})|_{z=0,\mu>0}\,.$
(34)

$I_{1}$ and $I_{2}$ are positive intensities and for
$\frac{1}{2}(\mathcal{A}_{1}-\mathcal{A}_{2})(\mathcal{T}((I_{1}+I_{2})|_{z=0,\mu>0}$
and the righthandside of (34), (31) is also true. $\square$

###  6.2 The frequency dependent case

Let us come back to (5) with (26b). Let us denote $f_{+}=\max(f,0)$. When not
ambiguous, let us write $I_{\nu}(\mu)$ instead of $I_{\nu}(0,\mu)$.

###### Definition 1

An operator $\mathcal{A}$ defined in $L^{1}_{\mu}(0,1)$ is non-accretive if

$\forall I^{1},I^{2}\in
L^{1}_{\mu}(0,1),~{}\int_{0}^{1}\mu\left((I^{2}-I^{1})_{+}-(\mathcal{A}(I^{2})-\mathcal{A}(I^{1}))_{+}\right)\mathrm{d}\mu\geq
0\,.$ (35)

###### Remark 5

Obviously (35) follows from the same pointwise property:

$\forall(I^{1},I^{2})\quad((\mathcal{A}(I^{2})-\mathcal{A}(I^{1}))_{+}\leq(I^{2}-I^{1})_{+},\quad\mu>0\,.$
(36)

Let us uses the symbol
$\Big{\langle}f\Big{\rangle}=\int_{0}^{\infty}\tfrac{1}{2}\int_{-1}^{1}f_{\nu}(\mu)\mathrm{d}\mu\mathrm{d}\nu$.

Then for two solutions $(I^{2}_{\nu},T^{2})$, $(I^{1}_{\nu},T^{1})$, by (5),
according to Theorem 4.1 in FGOP  one has:

$\displaystyle\begin{aligned}
-\frac{d}{dz}\langle\mu(I^{2}_{\nu}-I^{1}_{\nu})_{+}\rangle=\\\
=\left\langle\rho\bar{\kappa}_{\nu}(1-a_{\nu})((I^{2}_{\nu}-I^{1}_{\nu})-(B_{\nu}(T^{2})-B_{\nu}(T^{1})))\mathbf{1}_{I^{2}_{\nu}>I^{1}_{\nu}}\right\rangle\\\
+\left\langle\rho\bar{\kappa}_{\nu}a_{\nu}\left((I^{2}_{\nu}-I^{1}_{\nu})-\tfrac{1}{2}\int_{-1}^{1}p(\mu,\mu^{\prime})(I^{2}_{\nu}-I^{1}_{\nu})(z,\mu^{\prime})d\mu^{\prime}\right)\mathbf{1}_{I^{2}_{\nu}>I^{1}_{\nu}}\right\rangle\geq
0\end{aligned}$ (37)

and therefore
$\frac{d}{dz}\Big{\langle}\mu(I^{2}_{\nu}(z,\mu)-I^{1}_{\nu}(z,\mu))_{+}\Big{\rangle}\leq
0\,.$ On the other hand

$I_{\nu}^{2}(H,-\mu)\leq
I_{\nu}^{1}(H,-\mu),\quad\mu>0\quad\implies~{}~{}\Big{\langle}\mu(I_{\nu}^{2}(H,\mu)-I^{1}_{\nu}(H,\mu))_{+}\Big{\rangle}=0$

while the non accretivity of $\mathcal{A}$ gives

$\displaystyle\begin{aligned}
\langle\mu(I^{2}_{\nu}(0,\mu)-I^{1}_{\nu}(0,\mu))_{+}\rangle\\\
=\tfrac{1}{2}\int_{0}^{\infty}d\nu\int_{0}^{1}\mu(\mathcal{A}(I^{2}_{\nu}(0,-\mu))-\mathcal{A}(I^{1}(0,-\mu)))_{+}-(I^{2}_{\nu}-I^{1}_{\nu})_{+}(0,-\mu))d\mu\leq
0\end{aligned}$ (38)

Therefore, unless
$\langle\mu(I^{2}_{\nu}(0,\mu)-I^{1}_{\nu}(0,\mu))_{+}\rangle=0$ there is a
contradiction because it is a function which is decreasing, negative at $z=0$
and zero at $z=H$. Then consider the difference of equations (5) with
arguments $(I_{\nu}^{i},T^{i})\,i=1,2$ multiplied by
$\frac{\mu}{\bar{\kappa}_{\nu}}{\bf 1}_{I^{2}_{\nu}>I^{1}_{\nu}}$. Following
the end of the proof of Theorem 4.1 in FGOP  (see also G87 ), one finds by the
same method that existence and uniqueness derives from

$\frac{d}{dz}\left\langle\frac{\mu^{2}}{\bar{\kappa}_{\nu}}(I^{2}_{\nu}-I^{1}_{\nu})_{+}\right\rangle=0\,,\qquad
0<z<H\,.$

$\square$
Eventually one has:

###### Theorem 2

Consider (5) with incoming data $0\leq I_{\nu}(H,-\mu)\leq B_{\nu}(T_{M})$,
depending on a given (but possibly large) temperature $T_{M}$, and with
outgoing data given by an accretive operator, of the form
$I_{\nu}(0,\mu)|_{\mu>0}={\mathcal{A}}(I_{\nu}(0,-\mu))|_{\mu>0}$. The problem
has a unique well defined nonnegative solution.

###  6.3 Examples of Albedo operators

Below are given some examples of albedo operators which combine an
accommodation parameter $\alpha$ and the reflection and the thermalisation
effects of the Earth.

The simplest example would be given (cf. Remark 3) with $0\leq\alpha\leq$ by
the formula, $p>0$:

$\mathcal{A}(I(0,\mu))=\alpha
I(0,-\mu)+(1-\alpha)\int_{0}^{1}\mu^{p}I(0,-\mu^{\prime})\mathrm{d}\mu^{\prime}\,.$
(39)

For the more realistic frequency dependent cases, one should compare the
reflection with the emission from the Earth as global black body under an
effective temperature $T_{e}\simeq 288$ and this leads to try the albedo
operator

$\mathcal{A}(I_{\nu}(0,-\mu))=\alpha
I_{\nu}(0,-\mu)+(1-\alpha)B_{\nu}(T_{e}),~{}~{}\mu>0.$ (40)

This obviously satisfies the hypothesis of non accretivity.

##  7  Calculus of Variations for the $\nu$-dependent Case

The following argument sheds some light on the conditions needed to obtain a
cooling (resp. heating) from a local increase of the absorption coefficient.

###### Conjecture 1

Let $I^{\varepsilon}_{\nu}(\tau,\mu),T^{\varepsilon}(\tau)$ be the solution of
(5) when

$\kappa^{\epsilon}_{\nu}:=\kappa+\epsilon\delta\kappa_{\nu}\,,$ (41)

Let $I^{0}_{\nu}(\tau,\mu),T^{0}(\tau)$ be the solution of (5) with
$\kappa_{\nu}=\kappa$ constant. When
$\delta\kappa_{\nu}\approx{\bm{\delta}_{\nu^{*}}}$, the Dirac mass at
$\nu^{*}$, and $\kappa$ is small, then the sign of
$\displaystyle\frac{dT}{d{\varepsilon}}(\tau)|_{{\varepsilon}=0}$ is governed
by the sign of

$\frac{1}{2}\int_{-1}^{1}I^{0}_{\nu^{*}}(\tau)d\mu-B_{\nu^{*}}(T^{0}(\tau)).$
(42)

###  7.1 Calculus of Variations Support of the Conjecture

It will be convenient to define the Planck function in terms of the quantity
given by the Stefan-Boltzmann law $\Phi:=\sigma T^{4}\,.$ Henceforth we denote

$B_{\nu}(T)=b_{\nu}(\Phi)\,.$

We seek to study how the average radiation intensity or the temperature is
altered if the absorption is modified on various intervals in the frequency
variable. We expect that the simplest situation corresponds to absorption of
the form (41) with $\kappa>0$ independent of $\nu$, while $0<\epsilon\ll 1$.
This is of course an extremely general formulation, but in practice one could
think of

$\delta\kappa_{\nu}:=\kappa\mathbf{1}_{\nu_{1}<\nu<\nu_{2}}\,.$

This corresponds to multiplying the absorption by $(1+\epsilon)$ in the
frequency interval $(\nu_{1},\nu_{2})$ only, and leaving it invariant for all
the other frequencies.

We assume that the scattering is isotropic with constant rate $\lambda:=\kappa
a\geq 0$.

Henceforth, we consider the radiative transfer equation on
$(0,Z)\times(-1,1)\times{\mathbb{R}}$,

$\mu\partial_{\tau}I^{\epsilon}_{\nu}+\kappa^{\epsilon}_{\nu}(I^{\epsilon}_{\nu}-b_{\nu}(\Phi^{\epsilon}))+\lambda(I^{\epsilon}_{\nu}-\widetilde{I}^{\epsilon}_{\nu})=0\,,\qquad\Big{\langle}\kappa^{\epsilon}_{\nu}(\widetilde{I}^{\epsilon}_{\nu}-b_{\nu}(\Phi^{\epsilon}))\Big{\rangle}=0\,,$
(43)

with the boundary conditions

$I^{\epsilon}_{\nu}(0,\mu)=\mu Q^{0}_{\nu},\quad
I^{\epsilon}_{\nu}(Z,-\mu)=0\,,\qquad\qquad 0<\\!\mu\\!<1,$ (44)

and the notations

$\widetilde{\phi}:=\tfrac{1}{2}\int_{-1}^{1}\phi(\mu)d\mu\,,\qquad\Big{\langle}\psi\Big{\rangle}:=\int_{0}^{\infty}\psi(\nu)d\nu\,.$

First we study the case $\epsilon=0$: since $\kappa$ is constant, one can
average in frequency, and set

${\bar{I}^{0}}(\tau,\mu):=\Big{\langle}I^{0}_{\nu}(\tau,\mu)\Big{\rangle},\qquad\widetilde{\bar{I}^{0}}(\tau)=\Big{\langle}b_{\nu}(\Phi^{0})\Big{\rangle}=\Phi^{0}(\tau)\,.$

One finds

$\mu\partial_{\tau}{\bar{I}^{0}}(\tau,\mu)+\kappa({\bar{I}^{0}}(\tau,\mu)-\Phi^{0})+\lambda({\bar{I}^{0}}(\tau,\mu)-\widetilde{\bar{I}^{0}}(\tau))=0\,,$

or, equivalently

$\displaystyle\mu\partial_{\tau}{\bar{I}^{0}}(\tau,\mu)+(\kappa+\lambda)({\bar{I}^{0}}(\tau,\mu)-\widetilde{\bar{I}^{0}}(\tau))=0\,,\qquad\tau>0\,,$
$\displaystyle{\bar{I}^{0}}(0,\mu)=\mu\Big{\langle}Q^{0}_{\nu}\Big{\rangle}\,,\qquad\qquad\qquad\qquad\qquad\quad
0<\\!\mu\\!<1.$

Once $\widetilde{\bar{I}^{0}}$ is known, one recovers $I^{0}_{\nu}$ by using
the semi-analytical formula (11), which in this case is

$\displaystyle I^{0}_{\nu}(\tau,\mu)=$ $\displaystyle\mathbf{1}_{0<\mu<1}\mu
e^{-\frac{(\kappa+\lambda)\tau}{\mu}}Q^{0}_{\nu}+\mathbf{1}_{0<\mu<1}\int_{0}^{\tau}e^{-\frac{(\kappa+\lambda)(\tau-y)}{\mu}}\frac{\kappa
b_{\nu}(\widetilde{\bar{I}^{0}}(y))+\lambda\widetilde{I}^{0}_{\nu}(y)}{\mu}dy$
(45)
$\displaystyle+\mathbf{1}_{-1<\mu<0}\int_{\tau}^{\infty}e^{-\frac{(\kappa+\lambda)(y-\tau)}{|\mu|}}\frac{\kappa
b_{\nu}(\widetilde{\bar{I}^{0}}(y))+\lambda\widetilde{I}^{0}_{\nu}(y)}{|\mu|}dy\,.$

This is precisely the solution of the Milne problem in the case of radiative
transfer in a grey atmosphere.

In order to study the effect of the perturbation of the absorption defined
above, we seek to compute

$I^{\prime}_{\nu}(\tau,\mu)=\frac{\partial
I_{\nu}^{\epsilon}(\tau,\mu)}{\partial\epsilon}\Big{|}_{\epsilon=0}\,,\quad\Phi^{\prime}(\tau)=\frac{\partial\Phi^{\epsilon}(\tau)}{\partial\epsilon}\Big{|}_{\epsilon=0}\,,\quad
u^{\prime}(\tau,\mu)=\frac{\partial
u^{\epsilon}(\tau,\mu)}{\partial\epsilon}\Big{|}_{\epsilon=0}\,.$

Let $\dot{b}_{\nu}(\Phi)=\partial_{\Phi}b_{\nu}(\Phi)$. Let us differentiate
in $\epsilon$ both sides of the equations (43) at $\epsilon=0$: one easily
finds that

$\displaystyle\mu\partial_{\tau}I^{\prime}_{\nu}+\kappa(I^{\prime}_{\nu}-\dot{b}_{\nu}(\Phi^{0})\Phi^{\prime}_{0})+\lambda(I^{\prime}_{\nu}-\widetilde{I}^{\prime}_{\nu})+\delta\kappa_{\nu}(I^{0}_{\nu}-b_{\nu}(\Phi^{0}))=0\,,$
$\displaystyle\kappa\Big{\langle}\widetilde{I}^{\prime}_{\nu}\Big{\rangle}-\kappa\underbrace{\Big{\langle}\dot{b}_{\nu}(\Phi^{0})\Big{\rangle}}_{=1}\Phi^{\prime}+\Big{\langle}\delta\kappa_{\nu}(\widetilde{I}^{0}_{\nu}-b_{\nu}(\Phi^{0}))\Big{\rangle}=0\,,$

and we can recast the second equality as

$\Phi^{\prime}(\tau)=\widetilde{\bar{I^{0}}}^{\prime}(\tau)+\tfrac{1}{\kappa}\Big{\langle}\delta\kappa_{\nu}(\widetilde{I}^{0}_{\nu}(\tau)-b_{\nu}(\Phi^{0}(\tau)))\Big{\rangle}=\widetilde{\bar{I^{0}}}^{\prime}(\tau)+\tfrac{1}{\kappa}\Big{\langle}\delta\kappa_{\nu}(\widetilde{I}^{0}_{\nu}(\tau)-b_{\nu}(\widetilde{\bar{I}^{0}}(\tau)))\Big{\rangle}\,.$
(46)

Averaging in frequency the RT equation, we find that

$\mu\partial_{\tau}{\bar{I^{0}}}^{\prime}(\tau,\mu)+\kappa({\bar{I^{0}}}^{\prime}(\tau,\mu)-\underbrace{\Big{\langle}\dot{b}_{\nu}(\Phi^{0}(\tau))\Big{\rangle}}_{=1}\Phi^{\prime}(\tau))+\Big{\langle}\delta\kappa_{\nu}(I^{0}_{\nu}(\tau,\mu)-b_{\nu}(\Phi^{0}(\tau)))\Big{\rangle}=0\,,$

and we can further eliminate $\Phi^{\prime}(\tau)$ between the last two
equations, to find,

$\mu\partial_{\tau}{\bar{I^{0}}}^{\prime}(\tau,\mu)+(\kappa+\lambda)({\bar{I^{0}}}^{\prime}(\tau,\mu)-\widetilde{\bar{I^{0}}}^{\prime}(\tau))=-\Big{\langle}\delta\kappa_{\nu}(I^{0}_{\nu}(\tau,\mu)-\widetilde{I}^{0}_{\nu}(\tau))\Big{\rangle}\,.$
(47)

If $\delta\kappa$ is small everywhere but large near $\nu^{*}$ and if
$\kappa+\lambda\ll 1$ so that the righthandside dominates the seccond term on
the left, then (47) becomes

$\mu\partial_{\tau}{\bar{I^{0}}}^{\prime}(\tau,\mu)=-\delta\kappa(I^{0}_{\nu^{*}}(\tau,\mu)-\widetilde{I}^{0}_{\nu^{*}}(\tau))\,.$

Together with (46) it gives the sign of $T^{\prime}(z)$ by,

$\Phi^{\prime}(\tau)=-\delta\kappa\Big{|}\int^{\tau}\widetilde{\frac{I^{0}_{\nu^{*}}(t,\mu)}{\mu}}dt\Big{|}+\frac{\delta\kappa}{\kappa}(\widetilde{I}^{0}_{\nu^{*}}(\tau)-b_{\nu^{*}}(\widetilde{\bar{I}^{0}}(\tau)))\,.$
(48)

when $\kappa\ll 1$ the second terms dominate the first one and it is possible
to obtain a positive sign when
$\widetilde{I}^{0}_{\nu^{*}}(\tau)>b_{\nu^{*}}(\widetilde{\bar{I}^{0}}(\tau))$.

###  7.2 Numerical Simulations

We study the sign of (42). As it is a complex surface of $\nu$ and $\tau$, we
look at the frequencies and altitudes for which (42) changes sign when one of
the two variable varies. Figure 3 shows the frequencies and altitudes where
there is a change of sign, for the grey case $\kappa=0.5$.

$0$$1$$2$$3$$4$$5$$6$$0.45$$0.5$$0.55$Altitude - kmFrequency$\times
10^{14}$Switch sign $(\nu,z)$ points

$0$$1$$2$$3$$4$$5$$6$$0.35$$0.4$$0.45$$0.5$$0.55$Altitude - kmFrequency$\times
10^{14}$Switch sign $(\nu,z)$ points

Figure 3: Switch sign points: left when $\nu$ is increased at fixed $z$. Right
: when altitude is increase at fixed $\nu$.

Consider now $\kappa_{\nu}=0.5+\delta\kappa{\bf 1_{(\nu^{1},\nu^{2})}}$. We
use an AD (automatic differentiation) extension of the computer program for
(45) and display the sensitivity of $T$ with respect to $\delta\kappa$ for 3
ranges of support $(\nu^{1},\nu^{2})$.

On Figure 4 it is seen that when $\nu^{1}=0.2$, $\nu^{2}=0.3$ the sensitivity
of $T$ is negative, meaning that increasing $\kappa_{\nu}$ in the range
$(\nu^{1},\nu^{2})$ triggers a decrease of temperature at all altitudes.
Conversely with $\nu^{1}=0.6$, $\nu^{2}=0.8$ the sensitivity is positive above
800m meaning that an increase of $\kappa_{\nu}$ in that range triggers a
temperature decrease below 800m and an increase above 800m. The same is true
for $T^{\prime}_{1}$ but at 1500m instead of 800m.

$0$$1$$2$$3$$4$$5$$-4$$-3$$-2$$-1$$0$$1$$\cdot 10^{-3}$Altitude -
kmTemperature Sensitivities
(K)$T^{\prime}_{0}$$T^{\prime}_{1}$$T^{\prime}_{2}$ Figure 4:  Scaled
temperatures sensitivities $z\to T^{\prime}(z)$ versus altitude, computed with
$\kappa_{\nu}=0.5+\delta\kappa{\bf 1_{(\nu^{1},\nu^{2})}}$ and $\delta\kappa$,
at $\delta\kappa=0$ when case 0: $\nu^{1}=0.2$, $\nu^{2}=0.3$. Case 1:
$\nu^{1}=0.3$, $\nu^{2}=0.4$. Case 2 is with $\nu^{1}=0.6$, $\nu^{2}=0.8$ .
Notice the change of signs in the sensitivity $T^{\prime}_{2}$.

## References

  * [1] C. Bardos and O Pironneau.  Radiative Transfer for the Greenhouse Effect.  submitted to SeMA J. Springer, 2021.
  * [2] S. Chandrasekhar.  Radiative Transfer.  Clarendon Press, Oxford, 1950.
  * [3] J. Dufresne, V. Eymet, C. Crevoisier, and J. Grandpeix.  Greenhouse effect: The relative contributions of emission height and total absorption.  Journal of Climate, American Meteorological Society, 33(9):3827–3844, 2020.
  * [4] Andrew Fowler.  Mathematical geoscience.  Berlin: Springer, 2011.
  * [5] F. Golse.  The Milne Problem for the Radiative Transfer Equations (with Frequency Dependence).  Trans. Am. Math. Soc., 303:125–143, 1987.
  * [6] F. Golse and O. Pironneau.  Stratified radiative transfer in a fluid.  Nov. 2021.
  * [7] F. Golse and O. Pironneau.  Stratified radiative transfer for multidimensional fluids.  Compte-Rendus de l’Académie des Sciences (Mécanique), to appear 2022\.
  * [8] H. Kaper and H. Hengler.  Mathematics & Climate.  SIAM publications, Philadelphia., 2013.
  * [9] Michael F. Modest.  Radiative Heat transfer.  Academic Press, 2013.
  * [10] O. Pironneau.  A fast and accurate numerical method for radiative transfer in the atmosphere.  Compte Rendus de l’académie des sciences (Math), December 2021.

⬇

1

2// main.cpp : solves the radiative trnsfer equations

3// for any absorption coefficient kappa

4// does automatic differentiation if _NOAUTODIF is undefined

5// Created by Olivier Pironneau on 10/10/2021.

6//

7

8#include <iostream>

9#include <fstream>

10#include <cmath>

11#include <string>

12#include <time.h>

13using namespace std;

14#define sqr(x) (x*x)

15

16#define _NOAUTODIF

17

18#if defined(_NOAUTODIF)

19 #define ddouble double

20 #define ffabs(x) (fabs(x))

21 const double dknu=0;

22#else

23 #include ”ddouble.h”

24 #define ffabs(x) (fabs(x.val[0]))

25 const ddouble dknu(0.,1.); // means dknu=0.1 and AD computes d_T/d_dknu

26#endif

27const bool greycase=true; // see beginning of main() for parameters. If
false use below

28const string kappa012=”kappa0f”; // defines 3 runs with kappa vs nu in
columns 2,3,4

29const bool verbose=false;

30const int Ntau=60; // nb points in tau

31const int kmax=16; // nb fixed point iterations

32const double Z=1-exp(-12.); // max tau after change of var

33const double SBsun =2.03e-5*(2*sqrt(2.))/0.7; // scaled sunlight power

34const double Tsun = 1.209; // scaled sun temperature

35const int jmaxmax=600; // max of max nb of points for integration in nu
range

36const int newton = 50; // to compute the temperature from int
k*Botlzmann=int k*Imean

37const double epsdycho=0.01, epsnewton=1.e-12; // precision for dychotomy
before Newton

38const double tmin=1.e-10; // min t in ExpInt(t)

39const double dtt = 0.005; //min integration step size in ExpInt integrals

40const double kappamin=0.001; // if kappa read is too small max it with
kappamin

41const int nt = 5; // min nb of integration step in anal formula

42const double ealb=0.3; // and Earth albedo

43const double alpha=1; // proportion of light affected at altitude 0 vs Z

44const double ais=0., ars=0.; // max isotropic and Rayleigh scattering

45const double tm1 = Z*0.6, tm2=Z*0.9; // altitudes for scattering

46double nu[jmaxmax],// uneven discretization of [numin,numax]

47 ai[jmaxmax], ar[jmaxmax], aux[jmaxmax];// isotropic and Rayleigh scattering

48ddouble Inut[Ntau], // mu integral of I_nu

49 Snut[Ntau], // mu integral of mu^2 I_nu

50 G[jmaxmax][Ntau], // J(nu,tau)

51 S[jmaxmax][Ntau], // K(nu,tau)

52 T[Ntau], kappanu[jmaxmax]; // T0 with kappanu0[]

53string basedir(”/Users/pironneau/Dropbox/afaire/Golse-Bardos/greenhouse4/”);

54string myresulttemperature(”temperature”);

55string myresultmeanintensity(”imean”);

56string mykappafile(basedir+kappa012+”.txt”); // has nu, kappa0, kappa1,
kappa2

57int jmax;

58

59struct ddouble2{ ddouble d1,d2;};

60

61// don’t use if kappa >18/Z

62ddouble expint_E1(const ddouble t=1){

63 // if your compiler has it or if you can link to gsl you may adapt this
function

64 // it computes E1(t)*B when t<18

65 double abst=ffabs(t);

66 const int K = 9+(abst-1)*4;; // precision in the exponential integral
function E1

67 const double gaNtaua =0.577215664901533; // special integration for log(t)

68 if(abst<tmin) return 0.;// because integral_0^t(logx)dx ~0

69 if(abst>18) {cout << ”value of E_1 is incorrect with ”<<t<<”>18”<<endl;
return 0;}

70 ddouble ak=(t<0)?-t: t, soNtaue=-gaNtaua - log(abst)+ak;

71 for(int k=2;k<K;k++){

72 ak *= -abst*(k-1)/sqr(k);

73 soNtaue += ak;

74 }

75 return soNtaue;

76}

77ddouble expint_E2(const ddouble t=1){

78 ddouble aux = exp(-t) - t*expint_E1(t);

79 return aux;

80}

81ddouble expint_E3(const ddouble t=1){

82 return (exp(-t) - t*expint_E2(t))/2;

83}

84ddouble expint_E4(const ddouble t=1){

85 return (exp(-t) - t*expint_E3(t))/3;

86}

87ddouble expint_E5(const ddouble t=1){

88 return (exp(-t) - t*expint_E4(t))/4;

89}

90

91int readkappa(string mykappafile, int which=0){

92 // on each line of kappafile: nu kappa0 & optional kappa1 kappa2

93 ifstream kappafile(mykappafile);

94 int j=-1;

95 double kappaux=0, dummy, nuj=0.01;

96 while((j++<=jmaxmax)&&(!kappafile.eof())){

97 if(which==0) kappafile >> nuj >>kappaux >> dummy >> dummy >> dummy;

98 else if(which==1) kappafile >> nuj >> dummy >>kappaux >> dummy >> dummy;

99 else if(which==2) kappafile >> nuj >> dummy >> dummy >>kappaux >> dummy;

100 kappanu[j] = fmax(kappaux,kappamin);

101 nu[j]=nuj;

102 }

103 kappafile.close();

104 jmax=j;

105 return 0;

106}

107

108double Bsun(const double nu){ return SBsun*sqr(nu)*nu/(exp(nu/Tsun) -1);}
// Boltzmann

109

110ddouble BB(const double nu, const ddouble T){

111 if(T<1.e-7) return 0;

112 if(nu<1e-10) return T*sqr(nu);

113 return sqr(nu)*nu/(exp(nu/T) -1);} // Boltzmann

114ddouble dBB(const double nu, const ddouble T){

115 if(T<1.e-7) return 0;

116 if(nu<1.e-10) return sqr(nu);

117 ddouble a = exp(nu/T); return a*sqr(nu*nu/(a -1)/T);

118}

119

120ddouble2 intBS(const int jnu,const double tau,const double tmin,const
double tmax){

121 // returns the convolution t-integral for the mean in mu of kappa*I_nu and
mu^2 kappa*I_nu

122 ddouble Imean=0, Imu2mean=0;

123 ddouble kappa=kappanu[jnu];

124 const double dt=fmin(dtt,nt/(tmax-tmin));

125 for(double t=tmin;t<tmax;t+=dt){

126 int it = int((Ntau-1)*t/Z); // parabolic length for frequency step

127 double
arnu4=ar[it]*sqr(nu[jnu]-0.8)*sqr(nu[jnu]-1.2)*(nu[jnu]>0.8)*(nu[jnu]<1.2)*40;

128 ddouble H0 = kappa*( BB(nu[jnu],T[it])*(1-arnu4) +

129 (ai[it]+1.125*arnu4)*G[jnu][it] -1.125*arnu4*S[jnu][it]) ;

130 ddouble H2 = -0.375*arnu4*kappa*(G[jnu][it] -3*S[jnu][it]);

131 if(kappa*(t-tau)!=0){

132 Imean +=
dt*H0*(expint_E1(kappa*fabs(tau-t))+ealb*expint_E1(kappa*(tau+t)))/2;

133 Imu2mean +=
dt*H0*(expint_E3(kappa*fabs(tau-t))+ealb*expint_E3(kappa*(tau+t)))/2;

134 if(H2 !=0){

135 Imean +=
dt*H2*(expint_E3(kappa*fabs(tau-t))+ealb*expint_E3(kappa*(tau+t)))/2;

136 Imu2mean +=
dt*H2*(expint_E5(kappa*fabs(tau-t))+ealb*expint_E5(kappa*(tau+t)))/2;

137 }

138 }

139 }

140 ddouble2 mean; mean.d1=Imean; mean.d2 = Imu2mean;

141 return mean;

142}

143

144ddouble thefunc(const ddouble rhs, const ddouble T0){

145 ddouble myeq=-rhs;

146 for(int j=1; j<jmax;j++){

147 myeq += kappanu[j]*BB(nu[j],T0)*(nu[j]-nu[j-1]);

148 }

149 return myeq;

150}

151ddouble getTbydycho(const int i,const ddouble Tstart){

152 ddouble T0=Tstart;

153 if(Tstart<0.1) T0=0.1;

154 ddouble Taux, myeq0=1,myeq1=-1, T1=T0, rhs=0;

155 for(int j=1; j<jmax;j++){

156 rhs+=kappanu[j]*G[j][i]*(nu[j]-nu[j-1]);

157 }

158 myeq0 =-rhs;myeq1=-rhs;

159 while (myeq0>0){

160 T0=T0/2;

161 myeq0=thefunc(rhs,T0);

162 if(verbose) cout<<T0<<” down ”<<myeq0<<endl;

163 }

164 while (myeq1<0){

165 T1=2*T1; myeq1=thefunc(rhs,T1);

166 if(verbose) cout<<T1<<” up ”<<myeq1<<endl;

167 }

168 while (T1-T0 > epsdycho){

169 Taux=(T1+T0)/2;

170 myeq0=thefunc(rhs,Taux);

171 if(myeq0>0) T1=Taux; else T0=Taux;

172 if(verbose) cout<<T0<<” middle ”<<T1<<” ”<<myeq0 << endl;

173 }

174 return (T1+T0)/2;

175}

176

177int genT(){

178 for(int i=0;i<Ntau;i++){

179 T[i] = getTbydycho(i,T[i]);

180 ddouble presfunc = 1;

181 int inewton =0;

182 while(inewton++<newton && ffabs(presfunc) > epsnewton){

183 ddouble T0 = T[i];

184 ddouble left=0, rhs=0, deriv=0;

185 double nu1=0;

186 for(int j=1; j<jmax;j++){

187 double dnu=nu[j]-nu[j-1]; // variable integral increment

188 nu1=(nu[j]+nu[j-1])/2;

189 rhs+=kappanu[j]*G[j][i]*dnu;

190 left += kappanu[j]*BB(nu1,T0)*dnu;

191 deriv += kappanu[j]*dBB(nu1,T0)*dnu;

192 }

193 presfunc = rhs-left;

194 if(ffabs(deriv)>1e-10) T[i] = T0+presfunc/deriv;

195 if(verbose) cout<<T[i]<<” newton ”<<presfunc<<endl;

196 }

197 if(inewton>=newton)cout << ”Newton precision doubtful ”<<endl;

198 }

199 return 0;

200}

201

202int getT(){ // return temperature by Kirchhof’s law when kappa is constant

203 const double pi = 4*atan(1.);

204 for(int i=0;i<jmax;i++){

205 ddouble rhs=0;

206 for(int j=1; j<jmax;j++)

207 rhs+=G[j][i]*(nu[j]-nu[j-1]);

208 T[i]=sqrt(sqrt(15*rhs))/pi;

209 }

210 return 0;

211}

212

213int getISnu(const int j){ // returns mean I_nu & mean mu^2 I_nu

214 double nu1 = nu[j];

215 ddouble kappa=kappanu[j];

216 if(kappa<0.01) kappa=0.01;

217 for(int i=0;i<Ntau;i++){

218 double x=i*Z/(Ntau-1);

219 ddouble2 aux = intBS(j,x,0,Z);

220 Inut[i] = aux.d1 + Bsun(nu1)*(alpha*expint_E3(kappa*(x+Z)) +
(1-alpha)*expint_E3(kappa*(Z-x))

221 )/2;// +ealb*alpha*expint_E3(kappa*(Z+x)))/2;

222 Snut[i] = aux.d2 + Bsun(nu1)*(alpha*expint_E5(kappa*x)
+(1-alpha)*expint_E5(kappa*(Z-x))

223 )/2;// +ealb*alpha*(expint_E5(kappa*(Z+x))))/2;

224 }

225 return 0;

226}

227

228int multiBlock(double initT)

229{

230 for(int i=0;i<Ntau;i++){ T[i]=initT; // initialize

231 for(int j=0; j<=jmax;j++){ G[j][i]=0; S[j][i]=0; }

232 }

233 for(int k=0;k<kmax; k++){ // fixed point loop: first update F

234 for(int i=0;i<Ntau;i++){

235 Inut[i]=0; Snut[i]=0;

236 }

237 for(int j=0; j<jmax;j++){

238 getISnu(j);

239 for(int i=0;i<Ntau;i++){

240 G[j][i]=Inut[i];

241 S[j][i]=Snut[i];

242 }

243 }

244 if( greycase )

245#if defined(_NOAUTODIF)

246 genT(); // replace by getT() to speed-up but know what you do!

247#else

248 genT();// Then update T. AD need genT() and false with getT() because
kappa varies

249#endif

250 else

251 genT();

252 ddouble normG=0, normS=0;

253 for(int j=0; j<jmax;j++)

254 for(int i=0;i<Ntau;i++){ normG+=ffabs(G[j][i]); normS+=ffabs(S[j][i]); }

255 cout << ”k= ”<<k <<” ”<<T[2]<<” ”<<normG<<” ”<<normS<<endl;

256// ofstream tempfile(basedir+”tempHistScatLow”+std::to_string(k)+”.txt”);
//study convergence

257// for(int i=1;i<Ntau;i++) tempfile<<-log(1-i*Z/(Ntau-1))<<”\t” <<T[i]
<<endl;

258 }

259 return 0;

260}

261

262int main(int argc, const char * argv[]) {

263 for(int it=0;it<Ntau;it++){ // defines scattering values vs tau

264 double tm = it*Z/Ntau;

265 ai[it] = ais*fmax(tm-tm1,0.)*fmax(tm2-tm,0.)*4/sqr(tm2-tm1);

266 ar[it] = ars*fmax(tm-tm2,0.)/(Z-tm2);

267 }

268 for(int which=0;which<1;which++){ // change to which<1 for a single
computation

269 if(greycase){

270 const double kappa0=0.5, numin = 0.05, numax = 15;

271 jmax = 400;

272 for(int j=0;j<jmax;j++){

273 nu[j]= numax/(1+(jmax-j)*(numax-numin)/numin/jmax); // uniform in
wavelength

274 // kappanu[j]=kappa0;//+dknu*(nu[j]>0.2)*(nu[j]<0.3);

275 kappanu[j]=kappa0*(nu[j]<6) + 0.1;

276 }

277 } else{

278 readkappa(mykappafile,which);

279// for(int j=0;j<jmax;j++) if((nu[j]<3./4)*(nu[j]>3./7))
kappanu[j]=fmax(0.5,1.5*kappanu[j]);

280 }

281 ofstream resultfile =

282 ofstream(basedir+myresulttemperature+to_string(which+10*(ais>0))+”.txt”);

283 cout<<”\n iterations \t [T] near earth [T] far ||G|| and ||S||\n”;

284 double t0 = clock();

285 multiBlock(0.);//,”tempmin”);

286 printf( ” Time CPU = %10.6f\n”, (clock() - t0)/CLOCKS_PER_SEC);

287 cout<<”\n tau\t [T]:”<<endl;

288 for(int i=1;i<Ntau;i++){

289 cout << -log(1-i*Z/(Ntau-1))<<”\t”<<T[i] <<endl;

290 resultfile<< -log(1-i*Z/(Ntau-1))<<”\t” // altitude

291 <<T[i]<<endl; // T(kappa) Milne by multigroup

292 }

293 ofstream imeanz0 =

294
ofstream(basedir+myresultmeanintensity+to_string(10*(ais>0)+which)+”0.txt”);

295 ofstream imeanzZ =

296
ofstream(basedir+myresultmeanintensity+to_string(10*(ais>0)+which)+”Z.txt”);

297 for(int j=10; j<jmax;j++){

298 imeanz0<< 3/nu[j] <<” ”<<1e5*G[j][1]<<endl; // mean intensity near the
ground

299 imeanzZ<< 3/nu[j] <<” ”<<1e5*G[j][Ntau-1]<<endl; // mean intensity at max
altitude

300 }

301 resultfile.close();

302 imeanz0.close();

303 }

304 return 0;

305}