#This script illustrates the runaway greenhouse phenomenon
#for a grey gas. The atmosphere consists of a one-component
#saturated condensible atmosphere, with T(p) determined by the
#Clausius-Clapeyron relation. The surface pressure (and hence
#optical thickness of the atmosphere) increases with T. The
#script generates OLR(Tsurf) for this case. It assumes constant
#specific absorption cross-section kappa, but more general
#cases can be treated by editing the integrand function f
#and the expression for tauInf
#
#In this case, we find the OLR by carrying out the definite
#integral giving the solution to the Schwarzschild equations.
#The integral is carried out by integrating from the top of the
#atmosphere down to the ground, using delTau = tauInf-tau as the
#integration variable.
#
#
#
#**ToDo: Turn the OLR computation into a function that
# returns a Curve object. That will make it easier
# for the students to modify the script to make graphs
# exploring the parameter dependence
#
#**ToDo: Replace the ODE integrator with call to romberg quadrature
#
#
#Data on section of text which this script is associated with
Chapter = '4.**'
Figure = '**'
#
from ClimateUtilities import *
import phys
from math import *
#Temperature profile. The argument is the ratio of
#pressure to the reference pressure
def T(pp0):
pp0 = max(pp0,1.e-20) #Avoids math range error at top of atmosphere
return T0/(1. - RTL*log(pp0))
#Integrand function, for OLR computation
#Function has a single parameter, which is tau0.
#"dummy" is an unused argument, but has to be there for
#the ODE integrator
def f(delTau,dummy,tau0):
delTau = min(delTau,100.)
pp0 = delTau/tau0
return phys.sigma*T(pp0)**4.*exp(-delTau)
#Function to compute OLR. This uses the differential
#equation integrator, but that could be replaced by
#quadrature, since we're just doing a definite integral
def OLR(ps):
tauInf = kappa*ps/g
tau0 = kappa*p0/g
#Determine number of subdivisions
n = int(10*tauInf) #Guarantee 10 subdivisions per unit optical depth
n = max(n,50) #Don't let n get too small when we're optically thin
ddelTau = tauInf/n
#
#Use the differential equation integrator to evaluate the
#definite integral of f(delTau) over the atmosphere
fi = integrator(f,0.,0.,ddelTau)
#Pass the parameter of the function to the integrator
#The integrator doesn't need to know tauInf, since that comes
#in only through the limit of integration and the boundary
#radiation term we add in at the end.
fi.setParams(tau0)
#
#Integration loop
#**ToDo: Could stop integrating when the contribution of the
#deeper atmosphere becomes too small. That would save some time.
for i in range(n-1):
ans = fi.next()
#
#Now we have to add in the contribution from the surface
tauInf = min(tauInf,100.) #Avoids underflow in optically thick case
return ans[1] + phys.sigma*T(ps/p0)**4.*exp(-tauInf)
#--------Main script starts here-------------------------------
#Constants for Clausius-Clapeyron. Thes are globals
T0 = 300. #Reference temperature
p0 = phys.satvpg(T0) #Water vapor saturation pressure
RTL = phys.water.R*T0/phys.water.L_vaporization
#Absorption constants. Globals (used in OLR function)
g = 10.
kappa = .1 #Very roughly appropriate for water vapor
#Note: If you want to set p0 so that kappa p0/g = 1, as in the text,
#and then compute the corresonding T(p0), you can do the following instead
#The answer should be the same as using the fixed reference pressure above.
#p0 = g/kappa
#The next line estimates the temperature at the pressure p0, using
# the simplified form of Clausius-Clapeyron
#T0 = 300./(1.-(phys.water.R*300./phys.water.L_vaporization)*log(p0/phys.satvpg(300.)))
#tau0 = 1.
#RTL= phys.water.R*T0/phys.water.L_vaporization
#List of surface pressures we want to do the computation for
psList = [10.*(1.06)**i for i in range(-100,120)]
TsList = [T(ps/p0) for ps in psList] #Corresponding temperatures
#The following statement computes the OLR for each surface pressure ps.
#By Clausius-Clapeyron, ps also determines Ts. It's more
#convenient to loop over ps than Ts, because we already have
#a function that computes T(p).
OLRList = [OLR(ps) for ps in psList]
#Plot results
c = Curve()
c.addCurve(TsList)
c.addCurve(OLRList)
w = plot(c)