(* Model from Boeker *)
(* parameters *)
{als,ta,ala,ta2,ala2} ={ 0.11, 0.53,0.30,0.06,0.31}
cvec=2.7

combi=ta*(1-als)+(1-ala)

(* constants *)
sigma=5.67*^-8
T0=(1370/4/sigma)^0.25
gamma=cvec/T0^3/sigma

(* Solve the energy flow equations by iteration *)
{Gtemp,Htemp}={0.5,0.5}
Do[
G=(1-ala-ta2*Htemp^4)^0.25;
Hp4=(combi+gamma*(Gtemp-Htemp))/(1-ala2+ta2);
H=(Hp4)^0.25;
{Gtemp,Htemp}={G,H};
{Git[i],Hit[i]}={G,H},
{i,1,20}]
{Tatmos,Tsurf}={T0*G,T0*H}

(* Check that the iteration converged to a solution *)
thline=AbsoluteThickness[2]
epi={{Blue,Text["H : surface",{10,1.00}]},
{Red,Text["G : atmosphere",{10,0.98}]}}
ListPlot[
{Table[{i,Git[i]},{i,1,20}],
Table[{i,Hit[i]},{i,1,20}]},
PlotRange->{{0,20},Automatic},
PlotStyle->{{Red,thline},{Blue,thline}},
Joined->True,
Epilog->epi]