#!/usr/bin/perl
#
# Updated 14/11/2012
# A simulation of radiative transfer in the atmosphere due to the absorption of IR by CO2
# in the atmosphere. Only the main absorption band 13-17 microns is included. The ground
# temperature is assumed to be 288K and emit black body radiation. An 20% fraction
# of which is assumed to lie within the CO2 absorption band.
# No other greenhouse gases are considered. The atmosphere is divided into 150 x 100m
# layers, each with a single temperature and density. A lapse rate of 6.5 K/km is assumed
# up to 11 km and a fixed temperature above this. Beer-Lambert law is used to calculate
# the transmission and absorption in each layer. Each layer then emits IR in the same band
# according to Kirchoff law. The net upgoing IR flux for each layer is the integral of
# upgoing IR from the layers below minus the downgoing IR from above.
@radup=();
#
# define constants
#
$sigma = 0.0000000567;
#
# frac is the estimated energy fraction of the IR spectrum in the main CO2 band
# this really should be calculated properly in the next iteration.
#
$frac= 0.2;
$TS = 288;
#
# use an average lapse rate for Earth. The dry adiabatic lapse rate is ~10 deg/km
#
$lapse=0.0065;
$CO2 = 0.0003;
# Beer-Lambert parameter for main CO2 band.
$k = 1.48;
#
# emission from the ground
#
$rtup = $frac*$sigma*$TS**4;
$radup[0] = $rtup;
$l=0;
$atmrad=0.0;
#
# average scale height in the atmosphere. This actually depends on absolute temperature
# and so also should be calculated properly for each level/
#
$scale=8600.0;
#
# divide atmosphere into 100m slices
#
$a = 100.0;
for ($i=0;$i<151;$i++) {
$z=$i*100.0;
if ($i < 111)
{
$T = $TS - $z*$lapse;
}
#
# Calculate the partial pressure of CO2 at level $i
#
$P = $CO2*exp(-$z/$scale);
#
# trans is the transmitted radiation applying Beer-Lambert
# $absorp is the absorbed IR in the level $i
#
$trans=exp(-$k*$P*$a);
$absorp=1-$trans;
#
# Using Kirchoff law - an equal flux emitted by level $i to that absorbed )
#
#
$rtup = $radup[$i-1]*$trans + $absorp*$frac*$sigma*$T**4;
#
# now subtract the downward radiation from all the levels above $i.
#
$rdown = 0.0;
for ($j=$i+2;$j<153;$j++) {
$zz=$j*100.0;
$P = $CO2*exp(-$zz/$scale);
$local=exp(-$k*$P*$a);
$absorp=1-$local;
$TU=$TS - $zz*$lapse;
$remit=$absorp*$sigma*$frac*$TU**4;
#
# now integrate down through the layers
#
$rleft=$remit;
for ($m=$j;$m>$i+1;$m--) {
$zk = $m*100.0;
$Pk = $CO2*exp(-$zk/$scale);
$trans=exp(-$k*$Pk*$a);
$rleft=$rleft*$trans;
}
$rdown=$rdown+$rleft;
}
$radup[$i]=$rtup;
$rtup=$rtup-$rdown;
$atmrad = $rtup-$radup[0];
print $i.' '.$rtup.' '.$rdown," \n";
}
#