#-----------------------------------------------------------
#Script for stability of equilibria
#with ice-albedo feedback. It also provides
#a utility for finding all the equilibria in
#a given temperature interval, and allows for
#exploration of hysteresis phenomena. Two different
#ways are given to make a hysteresis loop: integrating
#the time-dependent equations, or using a (somewhat less
#reliable) Newton-method solver.
#-----------------------------------------------------------
#Data on section of text which this script is associated with
Chapter = '3.**'
Figure = 'None'
#
#ToDo: **It would be nice to be able to take input for
# the initial condition interactively, using a mouse click
#
# **Make L and prad into L(t) and prad(t), so that
# the script can be used to illustrate hysteresis
#
import math
from ClimateUtilities import *
import phys
#Define albedo and OLR functions here
#Radiating pressure pRad, in mb, is passed as a
#parameter. Note that the param argument could be
#used to pass an object containing many values.
def OLR(T,pRad):
return phys.sigma * (T**4.)*(pRad/1000.)**(4.*2./7.)
#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
#Thermal mass function, dE/dT
def M(T):
return 50.*2000.*1000.
#Rate function for use with integrator
# t is not always used, but it has to be
# in the argument list anyway. The
# parameters needed (solar constant and
# radiating pressure) are passed through the
# object params. Note that this object
# can also be used to pass L and pRad as
# functions of time, for hysteresis studies
def NetFlux(t,T,params):
L = params.L
pRad = params.pRad
return ((L/4.)*(1.-albedo(T)) - OLR(T,pRad) )/M(T)
#Function to integrate the time dependent
#temperature equation and return lists of
#time and temperature suitable for plotting.
#The argument is the initial temperature.
#
#If the optional Lfunction argument is present, it
#makes the solar constant L a function of time,
#so that a time-dependent hysteresis loop can
#be generated. tmax (days) is a global
def Tseries(Tinit,Lfunction=None):
# Set up the integrator
dt = 24.*3600. # Time step in seconds
model = integrator(NetFlux, 0.,Tinit,dt)
model.setParams(params)
#Do the integration
tList = []
TList = []
t = 0.
count = 0
while t < tmax*24*3600.:
if not (Lfunction == None):
params.L = Lfunction(t/(24.*3600.))
t,T = model.next()
if count%50 == 0:
tList.append(t/(24.*3600.))
TList.append(T)
count += 1
return tList,TList
#Dummy object for passing parameters. The
#class Dummy is defined in ClimateUtilities
params = Dummy()
#Globals
L = 1500.
params.L = L
alpha_ice = .6
alpha_0 = .28
T1 = 260.
T2 = 295.
pRad = 640.
params.pRad = pRad
tmax = 4000.
#
c = Curve()
first = True
for Tstart in [200. ,220.,240.,260.,265.,270.,280.,300.,310.,320.]:
t,T = Tseries(Tstart)
if first:
c.addCurve(t,'Time')
first = False
c.addCurve(T,'Tstart=%f'%Tstart)
c.Xlabel = 'time'
c.Ylabel = 'Temperature'
c.PlotTitle = 'Time series for various initial temperatures'
w = plot(c)
#------Hysteresis loop in the time-dependent system-----
#This function ramps up L from an initial value up
#to a maximum and then back down to a minimum, and repeates
#with a given period.
#Note t is in days. Period (days) is a global
def Lfun(t):
Lmax = 3000.
Lmin = 500.
return Lmin + .5*(Lmax-Lmin)*(math.sin(2.*math.pi*t/period)+1.)
period = 8000. #Try longer and shorter periods. With longer
#periods you get closer to the steady state for
#the instantaneous value of L. You have to
#go up to 64000 or more before things start to
#look really close to the hysteresis diagram
tmax = 3.*period #Integrate for 3 periods
c1 = Curve()
t,T = Tseries(200.,Lfun)
#This time, we want to plot the
#time series in L-T space, to
#see the effect of the hysteresis
#loop
LList = [Lfun(tt) for tt in t]
c1.addCurve(LList,'L')
c1.addCurve(T,'Temperature')
c1.Xlabel = 'Solar Constant (W/m**2)'
c1.Ylabel = 'Temperature'
c1.PlotTitle = 'Time dependence in bifurcation space'
plot(c1)
#
#-------One way to compute a hysteresis loop ------
#Section to compute a temperature hysteresis loop in which a parameter
#such as pRad is decreased and then increased again to its original
#value. This implementation uses a Newton solver to get a steady
#state for each value of the control parameter
#ToDo: Replace or augment the following with a loop that
#compiles a list of all equilibria in a given interval, and
#marks them as stable or unstable. Doing the graphics to
#show the stable and unstable branches for this data is a little
#tricky, though.
#The current implementation is written in a rather awkward
#and ad-hoc way, and should be cleaned up. This calculation works
#by using Newton's method to find a solution to the nonlinear
#equations determining balance. The problem is that there isn't any
#general way to determine which equilibrium is the right one to
#switch to, so the following code has some tricks in it based
#on what we know about the solutions of the system. It won't really
#be reliable for a general system.
#
#A more reliable approach to hysteresis is to make the hysteresis
#parameter time-dependant in the previous time-evolution code given
#above. If the parameter is made slowly varying enough, then
#the correct stable solutions will be automatically picked out.
#One can also then experiment with making the time-dependence faster,
#to see what happens. Note that in a hysteresis loop, the
#solution can never "pass through" an unstable equilibrium. It always
#goes directly from one equilibrium to a stable one.
#
#Perhaps this hysteresis code should be made into a separate script.
#
# We treat the parameters like pRad and L as function arguments rather
# than globals, since that makes it easier to change their values
# inside some other function. We don't need to change the
# constants of the albedo function very much, so it works well
# enough to leave them as globals.
#
#ToDo: Clean up the business of switching which control parameter
#we are varying. Perhaps just allow for a general p-L variation
#Reset the globals
#Globals
L = 1500.
params.L = L
alpha_ice = .6
alpha_0 = .28
T1 = 260.
T2 = 295.
pRad = 640.
params.pRad = pRad
tmax = 4000.
#
pStart = 1000.
gList1 = [pStart-10.*i for i in range(70)]
gList2 = gList1[:] #Makes a new copy of the list, which we will reverse
#and append to get our list for driving the hysteresis
#loop
gList2.reverse() #Reverses list in place!
gList = gList1 + gList2 #Our list of pRad, which decreases then increases
#Do the same to make a loop of L, if we want to do a hysteresis
#loop over solar constant instead
Lstart = 900.
LList1 = [Lstart + 10.*i for i in range(100)]
LList2 = LList1[:]
LList2.reverse()
LList = LList1+LList2
def f(T,params):
L = params.L
pRad = params.pRad
return (1.-albedo(T))*L/4. - OLR(T,pRad)
root = newtSolve(f)
root.setParams(params)
pRad = gList[0]
T = root(200.)
Thyst = []
cHyst = Curve()
#cHyst.addCurve(gList)
#params.L = 1370.*.7
#for pRad in gList:
# params.pRad = pRad
cHyst.addCurve(LList)
params.pRad = 670.
for params.L in LList:
T = root(T)
for trial in range(1,400):
if T != 'No Convergence':
Thyst.append(T)
if root.deriv(T,params) > 0.:
print 'Warning: branch went unstable at T=%f'%T
break # Get out of the loop. We're done
else:
#If convergence failed, that means the
#solution we were tracking probably failed
#to exist. Try again with different initial
#guesses until we find another solution
#Note that this code does not switch branches
#if the branch goes unstable, and does not guarantee
#that the new branch is stable. If an unstable branch
#is encountered, though, a warning is issued.
guess = Thyst[-1] + .5*trial*(-1)**trial
T = root(guess)
if T == 'No Convergence':
print 'Solution cannot be continued'
break
#Plot hysteresis graph
cHyst.addCurve(Thyst)
plot(cHyst)
#-------Yet another way to make the bifurcation diagram------
#
#ToDo: Find a way to mark unstable branches
#
#This is yet another way to get hysteresis information.
#This time, we find all solutions, stable and unstable,
#for each value of the control parameter. This is
#done by generating initial guesses by a coarse scan
#of temperature space. With this method,we can
#see the unstable branch explicitly. Note we could also
#keep track of which branches are stable, by looking
#at root.deriv(x).
#
#The approach used is brute-force and very inefficient, but
#it's good enough for such a simple problem. For the simple
#control parameter used here, the method in IceAlbedoZeroD.py
#is a lot simpler, but the technique below is useful when
#the dependence of the radiation balance on the control
#parameter is too complicated to use that method. As an example,
#you can try modifying the following code so that you do the
#bifurcation diagram with ice albedo as the control parameter.
#With this method, we don't need to vary the parameter up then
#down. We get all solutions on the first pass.
LList = [Lstart + 5.*i for i in range(200)]
gList = [pStart-5.*i for i in range(140)]
TList = [] #This time, it will be a list of lists,
#because of multiple solutions
params.pRad = 670.
for params.L in LList:
guesses = root.scan([200.,330.],500)
solutions = []
for guess in guesses:
solutions.append(root(guess))
if solutions[-1] == 'No Convergence':
solutions.pop()
TList.append(solutions)
#Now process the solution list into a single two-column
#table suitable for scatter-plotting
Lmerge = []
Tmerge = []
for i in range(len(LList)):
for T in TList[i]:
Lmerge.append(LList[i])
Tmerge.append(T)
c = Curve()
c.addCurve(Lmerge,'L')
c.addCurve(Tmerge,'T')
c.scatter['T'] = True #Do scatter plot, don't connect with lines
plot(c)