#!/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";
    
    }
#