#------------------------------------------------------------------------
#
#Description:
#
#This script computes and plots the OLR vs. surface
#temperature and pressure for a pure CO2 atmosphere on the dry adiabat,
#with or without the effects of condensation. This version
#is intended for surface pressures up to those of Early Mars
#For the very hot case of Venus, additional treatment of the
#shorter wave part of the spectrum is needed.  The approximate
#computation for Venus described in the book (where the emission
#from the shortwave region is replaced by a blackbody at a specified
#emission pressure) is done in the script SimpleVenusOLR.py

#Note that by just editing the thermodynamics and band data,
#this can be used for any pure condensible atmosphere, notably
#water vapor.  By calculating the plot of OLR vs surface temperature
#or pressure for an all-saturated atmosphere, we get the plot
#that tells us about the generalized runaway greenhouse effect
#(and Kombayashi-Ingersoll Limit), as well as information about
#whether continued buildup will lead to condensation limitation(
#as for water on Earth), or whether the greenhouse heating outruns
#the vapor pressure leading to gaseous accumulation.  On Venus,
#CO2 accumulates as a gas, but on Mars it eventually condenses at the
#surface. 
#
#It could be used for pure Methane as well, if
#the Methane continuum were implemented. At low temperatures
#where Methane condensation is important, all Methane absorption
#is in the continuum, so the radiative transfer doesn't
#require exponential sums at all.  For Titan, the CH4/N2 system
#is more interesting than pure CH4, though, and that would
#be a better basis for a "Methane Runaway" scenario.
#
#**ToDo:  Rationalize the setup of the calculation of the "CO2 runaway"
#graph so the calculation is more transparent.  Perhaps this
#graph should be moved to the RealGasRunaway script, which
#can generate the same curve, and also generalize the calculation
#to allow for a background gas

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

from ClimateUtilities import *
import phys,planets
import math

import miniClimtFancy as radmodel


#Path to the workbook datasets
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
#Path to the exponential sum tables
EsumPath = datapath + 'Chapter4Data/ExpSumTables/'

radmodel.Tstar = 900.

#
#Initialize the radiation model parameters
#
radmodel.BackgroundGas = phys.air #The transparent background gas (not used)
radmodel.GHG =phys.CO2 #The greenhouse gas
radmodel.g = planets.Mars.g
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.self.data')
#Eliminate transparent shortwave bands, as kludge for high-temperature Venus
#radmodel.bandData = radmodel.bandData[:-9] #Uncomment this for "Venus"
radmodel.SelfContinuum = radmodel.CO2SelfContinuum
radmodel.ForeignContinuum = radmodel.CO2ForeignContinuum

#Define adiabat functions
#Set needed thermodynamic constants
stuff = radmodel.GHG
L = stuff.L_sublimation # Latent heat
cp = stuff.cp #Specific heat
R = stuff.R #Gas constant for the atmosphere
Tr = stuff.TriplePointT # Reference temperature
psat_ref = stuff.TriplePointP  # Saturation vapor pressure at
                               #reference temperature

#Function that solves for temperature at which saturation
#occurs, for a given pressure.  This uses the simplified
#form of Clausius-Clapeyron based on the assumption that
#the latent heat is constant. Returns a whole array.
#Note that this assumes that the condensation is always into
#the same phase. For CO2, above 5 bars of pressure the condensation will
#be into liquid rather than solid, and the temperature profile
#should be changed accordingly
def Tsat(p):
    return(1./(1./Tr - (R/L)*Numeric.log(p/psat_ref)))

def OLR(ps,Tg,condense=True,n=20):
    p = radmodel.setpLog(ps,ptop,n)
    q = Numeric.ones(n,Numeric.Float)#1. #Pure greenhouse gas
    T = Tg*(p/ps)**(R/cp)
    #Handle condensation
    if condense:
        T = Numeric.where(T<Tsat(p),Tsat(p),T)
        if Tg < Tsat(ps):
            Tg = Tsat(ps)
    flux,heat = radmodel.radCompLW(p,T,Tg,q)
    #Add in contribution of part of spectrum outside the band table
    nu1 = radmodel.bandData[0].nu1
    nu2 = radmodel.bandData[-1].nu2
    m = romberg(radmodel.Planck) #Evaluates definite integral
    #Comment out the second addition to exclude shortwave
    #flux for doing the simplified Venus problem
    flux1 = m([.001,nu1],Tg) + m([nu2,50000.],Tg)
    return flux[0]+flux1

ptop = 1. #Top of atmosphere pressure,Pa.  A global default

#Do the plots. The default 20 points is OK, but for
#the book we used 40 points for extra accuracy
Tg = 270.
psList = [1000.*10.**(.1*i) for i in range(1000) if 1000.*10.**(.1*i)<=10.e5]
OLR1 = [OLR(ps,Tg,True) for ps in psList] #With condensation
OLR2 = [OLR(ps,Tg,False) for ps in psList] #Without condensation
#No-atmosphere flux, and Early Mars solar abs for reference
NoAtm = [phys.sigma*Tg**4 for ps in psList]
EarlyS = [.7*.65*planets.Mars.L/4. for ps in psList]
#
c = Curve()
c.addCurve(psList,'ps')
c.addCurve(OLR1,'OLRcond')
c.addCurve(OLR2,'OLRnocond')
c.addCurve(NoAtm,'NoAtmosphere')
c.addCurve(EarlyS,'EarlyMarsSolarFlux','Early Mars Solar Absorption')
c.XlogAxis = True
plot(c)

#Plot of critical solar flux to inhibit condensation level
#from reaching surface.  This is the same as the plot
#we use for the runaway greenhouse and finding the
#Kombayashi-Ingersoll limit. It's just the OLR vs
#surface temperature (or pressure) for an atmosphere
#in saturation everywhere.

#Setting Tg to zero forces the entire profile
#to be set to saturation, if condensation is included
psList = [1000.*10.**(.1*i) for i in range(1000) if 1000.*10.**(.1*i)<=50.e5]
c1 = Curve()
c1.addCurve(psList,'p')
for g in [1.,2.,10.,100.]:
    radmodel.g = g
    OLRcrit = [OLR(ps,0.,True) for ps in psList]
    c1.addCurve(OLRcrit,'OLRcrit.g=%f'%g)
c1.XlogAxis = True
plot(c1)

#Do CO2 runaway graph. This amounts to the same graph as the previous
#one, but plotted as a function of surface temperature instead. Since
#the atmosphere is assumed saturated all the way to the surface,
#the surface pressure is a function of the surface temperature
TgList = [120.+ 5.*i for i in range(30)]
psat = phys.satvps_function(phys.CO2)
psList = [psat(TT) for TT in TgList]
c2 = Curve()
c2.addCurve(TgList,'Tg')
for g in [1.,3.,10.,20.,100.]:
    radmodel.g = g
    OLRcrit = [OLR(ps,0.,True) for ps in psList]
    c2.addCurve(OLRcrit,'OLR.g=%f'%g)
plot(c2)
#Put in the surface pressure for later plotting, before dumping
c2.addCurve(psList,'ps')