#This script computes and plots the OLR vs. surface
#temperature and CO2 for the "canonical" case of
#a dry atmosphere with troposphere on the dry adiabat,
#patched to an isothermal stratosphere.  The ccm
#radiation model is called via the functions in ccmradFunctions,
#while the simplified exponential sum code is in miniClimt

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

import ccmradFunctions as ccm
from ClimateUtilities import *
import phys
import math

import miniClimt as homebrew
#import miniClimtFancy as homebrew #For testing temperature scaling
#homebrew.Tstar = 700. #900 looks pretty good

#We need to set the band data for the gas we are using, and
#also say what the greenhouse gas and the background gas is
datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
#Path to the exponential sum tables
EsumPath = datapath + 'Chapter4Data/ExpSumTables/'
#**ToDo: Put in guide to which band data to use. (allwave vs prinband, etc.)
#
homebrew.bandData = homebrew.loadExpSumTable(EsumPath+'CO2TablePrinBand.260K.100mb.air.data')
homebrew.BackgroundGas = phys.air #The transparent background gas
homebrew.GHG =phys.CO2 #The greenhouse gas
#
homebrew.pref = 1.e4 #Reference pressure at which band data is given
#Set the continuum, if there is one
homebrew.Continuum = homebrew.NoContinuum
#-------->Information for homebrew radiation computation is set

#Do this calculation with no stratosphere
ccm.setTstrat(0.)
#These are all done for the default gravity (Earth's)
#If you want another planet, you need to remember to
#set the gravity both in ccm and in homebrew

#First compute the OLR as a function of
#co2, for fixed ground temperature
Tg = 273.
co2List = []
olrccmList = []
olrhomebrewList = []
co2 = .1 #Molar concentration in ppmv
while co2 <= 1.e5:
    olrccmList.append(ccm.OLRDryAir(1.e5,Tg,co2))
    olrhomebrewList.append(homebrew.OLRDryAir(1.e5,Tg,co2))
    co2List.append(co2)
    co2 = co2* 2.**.25
c = Curve()
c.addCurve(co2List,'CO2','CO2 in ppmv')
c.addCurve(olrccmList,'ccmOLR','ccm OLR, W/m**2')
c.addCurve(olrhomebrewList,'homebrewOLR','homebrew OLR, W/m**2')
c.XlogAxis = True
plot(c)

#Now plot OLR vs Tg for 300ppmv
firstTime = True
for co2 in [300.,10000.]:
    TList = [200. + 2.*i for i in range(101)]
    olrccmList = []
    olrhomebrewList = []
    for Tg in TList:
        olrccmList.append(ccm.OLRDryAir(1.e5,Tg,co2))
        olrhomebrewList.append(homebrew.OLRDryAir(1.e5,Tg,co2))
    if firstTime:
        c1 = Curve()
        c1.addCurve(TList,'Tg')
        firstTime = False
    c1.addCurve(olrccmList,'ccmOLR%d'%co2,'ccm OLR%d, W/m**2'%co2)
    c1.addCurve(olrhomebrewList,'homebrewOLR%d'%co2,'homebrew OLR%d, W/m**2'%co2)
plot(c1)

#Compute prad/ps from the olr computations in the curve we just computed
c2 = Curve()
c2.addCurve(c1['Tg'],'Tg')
Tg = c1['Tg']
for name in c1.listVariables():
    if not (name == 'Tg'):
        olr = c1[name]
        Trad = (olr/5.67e-8)**.25
        pradNorm = (Trad/Tg)**(7./2.)
        c2.addCurve(pradNorm,name+'pradNorm')
plot(c2)