#Script for doing OLR for a single well-mixed
#GHG using the homebrew radiation model and choice of temperature
#profiles.  This does not incorporate water vapor effects. This
#script can be used to compare the radiative effects of different
#greenhouse gases in isolation (e.g. CO2 vs methane).
#
#As written, this script illustrates how to compare the
#radiative forcing for two different greenhouse gases
#(CO2 and CH4, in the example)

#The OLR function has been moved to this script, to
#make it easier to customize it to examine effects
#of different assumptions about the lapse rate, etc.
#As written, the OLR is computed for the "canonical"
#case consisting of a mixture of the greenhouse gas
#with a background gas having R/Cp = 2/7.  Note that
#the mixing ratio of the gas can be made non-uniform
#if desired, so this script can be modified
#to do a mixture of a condensible gas in a background,
#for example.

#Data on section of text which this script is associated with
Chapter = '4'
Section = 'section:OLR.CO2vsCH4'
Figure = 'fig:CanonicalOLR.CH4vsCO2'
#

from ClimateUtilities import *
import math,phys
import miniClimtFancy as homebrew

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

#OLR function. Molar concentration in ppmv
def OLR(psBackgroundGas,Tg,GHGmolarcon):
    ptop = 100.
    ps = psBackgroundGas/(1.-1.e-6*GHGmolarcon)
    Mbar = 1.e-6*GHGmolarcon*homebrew.GHG.MolecularWeight +\
    (1.-1.e-6*GHGmolarcon)*homebrew.BackgroundGas.MolecularWeight
    p = homebrew.setpLin(ps,ptop,40)
    #Set up the temperature profile
    T = Tg*(p/ps)**(2./7.) #Use fixed adiabatic exponent
    #Specify a uniform mass concentration profile
    q = 1.e-6*GHGmolarcon*homebrew.GHG.MolecularWeight/Mbar
    q = q*Numeric.ones(40,Numeric.Float)
    #
    #Do the radiation calculation
    flux,heat = homebrew.radCompLW(p,T,Tg,q)
    #Add in contribution of part of spectrum outside the band table
    nu1 = homebrew.bandData[0].nu1
    nu2 = homebrew.bandData[-1].nu2
    m = romberg(homebrew.Planck) #Evaluates definite integral
    flux1 = m([.001,nu1],Tg) + m([nu2,50000.],Tg)
    return flux[0] + flux1

#These calculations all done for the default gravity (Earth's)
#If you want another planet, you need to remember to
#set the gravity in homebrew, using, e.g:  homebrew.g = 2.7

#Set up the list of greenhouse gas concentrations (ppmv)
#ghgList = [10.**(.25*i) for i in range(-4,21)]+[2.5e5,5.e5]
ghgList = [1.,10.,100.,1000.,10000.]


#Gas dependent stuff here
homebrew.Tstar = 900. #Temperature scaling coefficient
#Load band data for your gas
#**ToDo: Put in guide to which band data to use. (allwave vs prinband, etc.)
#
homebrew.bandData = homebrew.loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.air.data')
homebrew.BackgroundGas = phys.air #The transparent background gas
homebrew.GHG =phys.CO2 #The greenhouse gas
gas1Name = homebrew.GHG.name
#
homebrew.pref = 1.e4 #Reference pressure at which band data is given
#Set the continuum, if there is one
homebrew.Continuum = homebrew.NoContinuum

Tg = 280. #Set the ground temperature
olr1 = [ OLR(1.e5,Tg,ghg) for ghg in ghgList]

# Now set up the band data for another gas and do
#it again

#Gas dependent stuff here
homebrew.Tstar = 900. #Temperature scaling coefficient
#Load band data for your gas
#**ToDo: Put in guide to which band data to use. (allwave vs prinband, etc.)
#
homebrew.bandData = homebrew.loadExpSumTable(EsumPath+'CH4TableAllWave.260K.100mb.air.data')
homebrew.BackgroundGas = phys.air #The transparent background gas
homebrew.GHG =phys.CH4 #The greenhouse gas
gas2Name = homebrew.GHG.name
#
homebrew.pref = 1.e4 #Reference pressure at which band data is given
#Set the continuum, if there is one
homebrew.Continuum = homebrew.NoContinuum

#Do the calculation for the new gas
olr2 = [ OLR(1.e5,Tg,ghg) for ghg in ghgList]

#Now plot the results
c = Curve()
c.addCurve(ghgList,'ppmv','Concentration in ppmv')
c.addCurve(olr1,'gas1',gas1Name)
c.addCurve(olr2,'gas2',gas2Name)
c.XLabel = 'Molar concentration in ppmv'
c.YLabel = 'OLR (W/m**2)'
c.XlogAxis = True
plot(c)


#If you uncomment the following, it will create a table interpolation
#function that allows you to quickly do further explorations
#of the dependence of the OLR on the greenhouse gas concentrations.
#If the absorptions are non-overlapping (e.g. co2 and ch4) you can
#add the two OLR reductions to get the joint effects of two gases
#
#These fits were used to make the graphs of OLR reduction as a
#function of ch4/(co2+ch4) ratio
x = [math.log(gas) for gas in ghgList]
dOLR1 = [phys.sigma*Tg**4 - olr for olr in olr1]
dOLR2 = [phys.sigma*Tg**4 - olr for olr in olr2]
dOLR1Fun = interp(x,dOLR1)
dOLR2Fun = interp(x,dOLR2)

#The following function defines the net OLR for the two gases
#It is written assuming that the total molar concentration of
#the two is kept fixed as the second is varied. This is right
#for studying conversion of ch4 to co2, but needs to be changed
#if the stoichiometry of your gases is different.
#**ToDo: give example for converting ethane to ch4 or co2
def net(ch4ppm,totppm):
	return dOLR2Fun(math.log(ch4ppm))+dOLR1Fun(math.log(totppm-ch4ppm))
#The following function makes a graph of net olr reduction
#as a function of the ratio of the second to the first gas, assuming
#fixed total ppmv of the two gases
def netGraph(totppmv,n=100):
	dg = (totppmv-1.)/n
	g2 = [1.+dg*i for i in range(n)]
	dolr = [net(g,totppmv) for g in g2]
	g2norm = [g/totppmv for g in g2]
	c = Curve()
	c.addCurve(g2norm,'RelCon')
	c.addCurve(dolr,'dolr%f'%totppmv)
	return c