#Script for computing the OLR for a grey-gas atmosphere
#with constant specific absorptivity and arbitrary T(p).
#By modifying the relation between tau1 and pps (which is p/psurf)
#in the derivative function dIpdtau(...), you can also use this
#script for alternate assumptions about absorptivity.
#
#In this script, we solve the Schwartzchild equation by direct
#integration as an ODE, even though the solution can be
#reduced to a definite integral
#
#The parameter called "accuracy" is a global which determines
#how many layers the atmosphere is divided up into in the numerical
#integration of the Schwartzschild equations.   Smaller values
#give more accurate results, but will require more computing time.
#A value of .1 is usually a good place to start, but you should check
#selected results by trying smaller values. 


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

import phys,math
from ClimateUtilities import *


#Specify T(pps) function here. pps is pressure
#divided by surface pressure. Tstrat and Ts are set
#as globals
def T(pps):
    #This if block takes care of round off error, which
    #can make pressure go slightly negative at the top of
    #the atmosphere sometimes
    if pps > 0.:
        Tadiabat = Ts*(pps)**Rcp
    else:
        Tadiabat = 0.
    #Return the adiabatic temperature, or the stratospheric
    #temperature, if it is less
    return max(Tstrat,Tadiabat)

#Derivative function, to hand to integrator. This computes
#the right-hand side of the differential equation for upward
#flux, I-plus. Modify expression for pps if you want to deal
#with non-constant absorptivity.
#
#  We use params to pass tauInf to the derivative function,
#  since using globals can be tricky in this script.
def dIpdtau(tau1,Ip,params):
    #Compute p/ps from tau1, then T from this
    pps = 1. - tau1/params.tauInf
    return -Ip + phys.sigma *T(pps)**4 

#This function uses Runge-Kutta integration to integrate
#the Schwartzschild equations to get the OLR
def OLR(tauInf):
    #Set up integrator
    dtau = min(accuracy,accuracy*tauInf)
    IpStart = phys.sigma*Tg**4
    m = integrator(dIpdtau,0.,IpStart,dtau)
    #The following sets up an object to pass parameters to
    #the derivative function
    params = Dummy() #This creates an empty object where you can
                     #put your parameters. Dummy() is a utility
                     #routine in ClimateUtilities
    params.tauInf = tauInf
    m.setParams(params)
    #
    #Carry out the integration
    tau1 = 0.
    while tau1 + dtau < tauInf+1.e-6:
        tau1,Ip = m.next()
    return Ip

#Set parameters governing the temperature profile
Tstrat = 200.
Rcp = 2./7.
Ts = 290.
Tg = 300. #Ground temperature, could be different from air temperature
#
#Set the accuracy parameter.
#Smaller is more accurate, but will take longer.
accuracy = .1

#Example of usage of the OLR function:
tauInf = 1.
print phys.sigma*Tg**4.,phys.sigma*Tstrat**4.,tauInf,OLR(tauInf)

#Example of plotting the OLR vs. optical thickness of the atmosphere
tauList = [.1*i for i in range(1,100)]
Tstrat = 200.
OLRList1 = [OLR(tau) for tau in tauList]
Tstrat = 250.
OLRList2 = [OLR(tau) for tau in tauList]
#Now check accuracy
Tstrat = 200.
accuracy = .01
OLRList1a = [OLR(tau) for tau in tauList] #Check accuracy

#Put results in curve and plot
c = Curve()
c.addCurve(tauList,'tauInf')
c.addCurve(OLRList1,'OLR200','Tstrat = 200K')
c.addCurve(OLRList1a,'OLR200a','Tstrat=200K, high accuracy')
c.addCurve(OLRList2,'OLR250','Tstrat = 250K')
w = plot(c)