#Script for illustrating simple models of
#ice-albedo feedback
#
#Data on section of text which this script is associated with
Chapter = '3.**'
Figure = '**'
#

from ClimateUtilities import *
import phys

#Define albedo and OLR functions here
#Radiating pressure pRad, in mb, is a global
def OLR(T):
    #return phys.sigma * (T**4.)*(pRad/1000.)**(4.*2./7.)
    return 112. + 2.177*(T-220.)

#Parameters of albedo function are globals
def albedo(T):
    if T < T1:
        return alpha_ice
    elif (T >= T1)&(T<=T2):
        r = (T-T2)**2/(T2-T1)**2
        return alpha_0 + (alpha_ice - alpha_0)*r
    else:
        return alpha_0


#Main body of script starts here
S = 1.*1370./4. #1.2 times modern give marginal case w/o greenhouse
alpha_ice = .6
alpha_0 = .1
T1 = 260.
T2 = 290.
pRad = 1000.

Tlist = [200.+ 2.*i for i in range(70)]
Slist = []
OLRlist = []
NetFluxlist = []
Glist = []

for T in Tlist:
    Slist.append(S*(1.-albedo(T)))
    OLRlist.append(OLR(T))
    NetFluxlist.append(S*(1.-albedo(T)) - OLR(T))
    Glist.append(OLR(T)/(1.-albedo(T)))

#Plot the results
c1 = Curve()
c1.addCurve(Tlist)
c1.addCurve(Slist,'Sabs','Absorbed Solar')
c1.addCurve(OLRlist,'OLR','Outgoing Longwave')
w1 = plot(c1)

c2 = Curve()
c2.addCurve(Tlist)
c2.addCurve(NetFluxlist,'Net','Net Flux')
c2.addCurve([0. for i in range(len(Glist))],'Zero','Equilibrium')
w2 = plot(c2)

c3 = Curve()
c3.addCurve(Tlist)
c3.addCurve(Glist,'G','Balance Function')
c3.addCurve([S for i in range(len(Glist))],'S','Incident Solar')
c3.switchXY = True #Switches axis so it looks like a hysteresis plot
w3 = plot(c3)

#Make a plot of the energy balance diagram with multiple
#different solar constants
Slist = [.9*S,S,1.1*S,1.7*S]
c4 = Curve()
c4.addCurve(Tlist,'T','Surface Temperature')
for S1 in Slist:
    Solar = [S1*(1.-albedo(T)) for T in Tlist]
    c4.addCurve(Solar,'%f'%(4*S1),'L = %f'%(4*S1))
c4.addCurve(OLRlist,'OLR','Outgoing Longwave')

w4 = plot(c4)