#This script computes and plots the OLR as a function
#of temperature with fixed CO2 and various relative
#humidities.  It also computes and plots the OLR as a function of CO2
#for fixed temperature. The atmospheric profile is
#the moist adiabat for water/air patched to an isothermal stratosphere.
#The calculation is done using the ccm radiation code

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

from ClimateUtilities import *
import ccmradFunctions as ccm


ps = 1.e5 #Surface partial pressure of air, in Pa
ccm.Tstrat = 0. # Temperature for isothermal stratosphere.
                # Zero for all-troposphere

#Plot OLR vs T for several different relative humidities
#for fixed CO2
co2 = 300.
ch4 = 1.e-30
def doPlot1(co2,ch4):
    TsList = [200.+i for i in range(120)]
    OLRlist1 = [ccm.OLR(ps,Ts,1.e-6,co2,ch4) for Ts in TsList]
    OLRlist2 = [ccm.OLR(ps,Ts,.1,co2,ch4) for Ts in TsList]
    OLRlist3 = [ccm.OLR(ps,Ts,.5,co2,ch4) for Ts in TsList]
    OLRlist4 = [ccm.OLR(ps,Ts,1.,co2,ch4) for Ts in TsList]

    c = Curve()
    c.addCurve(TsList,'T')
    c.addCurve(OLRlist1,'OLR1','Dry')
    c.addCurve(OLRlist2,'OLR2','RH = .1')
    c.addCurve(OLRlist3,'OLR3','RH=.5')
    c.addCurve(OLRlist4,'OLR4','Saturated')
    c.PlotTitle = 'OLR vs T for CO2 = %d ppmv'%co2
    return c
c = doPlot1(co2,ch4)
#c.dump('OLRvsT.txt')
plot(c)

#Plot OLR vs T for fixed rh, several different CO2
rh = .5
Ts = 280.
def doPlot2(rh,Ts):
    co2List = [10.,100.,1000.,10000.,100000.]
    c2 = Curve()
    c2.PlotTitle = 'OLR vs T for rh = %f'%rh
    TsList = [200.+i for i in range(120)]
    c2.addCurve(TsList,'Ts')
    for co2 in co2List:
        OLRlist = [ccm.OLR(ps,Ts,rh,co2,ch4) for Ts in TsList]
        c2.addCurve(OLRlist,'CO2%fppmv'%co2)
    return c2
c2 = doPlot2(rh,Ts)
plot(c2)

#Plot OLR vs CO2 for fixed Ts, and various rh
#Note: This could be modified to compare different assumptions
#on the same graph: e.g. different rh, and moist adiabat vs. dry adiabat
#
Ts = 280.
rhList = [1.e-6,.1,.5,1.]
ch4 = 1.e-30
def doPlot3(Ts,rhList,ch4):
    co2Start = .1
    co2End = 2.e5
    co2List = [co2Start*2.**(.25*i) for i in range(1000) if co2Start*2.**(.25*i) <= co2End]
    c1 = Curve()
    c1.addCurve(co2List,'co2')
    for rh in rhList:
        OLRlist = [ccm.OLR(ps,Ts,rh,co2,ch4) for co2 in co2List]      
        c1.addCurve(OLRlist,'OLRrh%f'%rh,'rh=%f'%rh)
    c1.XlogAxis = True
    c1.PlotTitle = 'OLR vs CO2 for Ts=%f'%Ts
    return c1
c1 = doPlot3(Ts,rhList,ch4)
#c1.dump('OLRvsco2.txt')
plot(c1)