#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.
#
#However, rather than introducing a numerical quadrature ("integration")
#routine, we make use of the differential equation solver in order
#to evaluate the definite integral.  This step could be replaced
#by another routine which evaluates definite integrals, such
#as the trapezoidal rule or more advanced methods.
#
#
#**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
#


#Data on section of text which this script is associated with
Chapter = '4.**'
Figure = '**'
#

from ClimateUtilities import *
import phys
from math import *

#Constants for Clausius-Clapeyron
T0 = 300. #Reference temperature
p0 = phys.satvpg(T0) #Water vapor saturation pressure
RTL = phys.water.R*T0/phys.water.L_vaporization

#Absorption constants
g = 10.
kappa = .1 #Very roughly appropriate for water vapor
tau0 = kappa*p0/g

#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


#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))

def f(delTau,dummy,params):
	delTau = min(delTau,100.)
	pp0 = delTau/tau0
	return phys.sigma*T(pp0)**4.*exp(-delTau)

psList = [10.*(1.06)**i for i in range(1,100)]
OLRList = []
TsList = []
#This loop computes the OLR for each surface pressure ps.
#By Clausius-Clapeyron, ps also determines Ts
for ps in psList:
 tauInf = kappa*ps/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
 #**ToDo: Could stop integrating when the contribution of the
 #deeper atmosphere becomes too small. That would save some time.
 fi = integrator(f,0.,0.,ddelTau)
 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
 OLR = ans[1] + phys.sigma*T(ps/p0)**4.*exp(-tauInf)
 #Save the results for plotting
 TsList.append(T(ps/p0))
 OLRList.append(OLR)
 #Print out what pressure we're at, to keep track of how things are going
 print ps/1.e5

#Plot results
c = Curve()
c.addCurve(TsList)
c.addCurve(OLRList)
w = plot(c)