#NOT YET VERIFIED
#---------------------------------------------------------------
#                            DESCRIPTION
# Computes OLR(Tg) for a binary atmosphere on the moist adiabat,
# with a noncondensible gas and a saturated condensible gas. The
# main use of this script is to study the runaway greenhouse, but
# it can also be used instead of the NCAR model to study water
# vapor feedback when there is no additional greenhouse gas besides
# the condensible.  It will work for any condensible, e.g.
# CO2 or CH4, though for cases other than H2O and CO2, miniClimtFancy
# would need to have the appropriate continuum absorption function
# put in. 
#
#
#---------------------------------------------------------------

#Data on section of text which this script is associated with
Chapter = '4'
Section = '{sec:RealGasRunaway}'
Figure = 'fig:PureWVRunaway,fig:WVwithN2Runaway'

import miniClimtFancy as rad
import phys
from ClimateUtilities import *



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


#Initialize the radiation calculation, and over-ride defaults
#For doing the pure water vapor case, we read in the self-broadened
#table and leave rad.SelfRat = 1
rad.bandData = rad.loadExpSumTable(EsumPath+'H2OTable.260K.100mb.self.data')
#
#For doing air/water mixtures, we read in the air-broadened table,
#and use a uniform rad.SelfRat. For high water vapor concentrations,
#this is less accurate than the previous method.
#rad.bandData = rad.loadExpSumTable('H2OTable.260K.100mb.air.data')
#rad.SelfRat = 5.

rad.ForeignContinuum = rad.NoContinuum #No w.v. foreign continuum
rad.BackgroundGas = phys.N2 #The transparent background gas
rad.GHG =phys.H2O #The greenhouse gas
rad.Tstar = 0. #Ignore temperature scaling except for continuum
rad.SelfContinuum = rad.H2OSelfContinuum

def OLR(Tg,psAir=1.e5,n=20):
    #-------------------Initialize arrays
    psat = phys.satvps_function(phys.water)
    ps = psAir + psat(Tg) #Total surface pressure for this temperature
    ptop = .1
    p = rad.setpLog(ps,ptop,n)
    #
    #Set temperature to moist adiabat and mass mixing ratio
    #to saturation
    m = phys.MoistAdiabat()
    m.ptop = ptop
    p,T,molarCon,q = m(psAir,Tg,p)
    #
    #Now compute the OLR and return it
    flux,heat = rad.radCompLW(p,T,Tg,q)
    return flux[0]

#Do the plots
#Copy over as many copies of the following example
#block as you need to cover the cases you want to treat
#(e.g. several background gas pressures or several values of g)
#
#The second argument to the OLR function is the background gas partial pressure
Tg = [270.+5.*i for i in range(27)]
rad.g=10. #Set the surface gravity
OLRList = [OLR(TT,1.e5,100) for TT in Tg]
c = Curve()
c.addCurve(Tg,'Tg')
c.addCurve(OLRList,'OLR')
c.dump('Runaway_g10_1barN2.100.txt')
plot(c)