#Description:
#======================================================================
#Script to solve the two-stream equations with scattering,
#by numerical integration as an ODE problem.
#
#This version is designed to work in the optically thick
#case. It computes the solution in the vicinity of
#optical thickness tau0 by integrating in both directions
#from tau0 and patching the solutions together. If
#one of the integrations hits a boundary before the
#exponential blowup criterion is exceeded, the correct
#boundary condition is applied there. If the blowup criterion
#is satisfied before hitting a boundary, then an arbitrary
#zeroing-out of the blowup mode is used as b.c, since the
#details don't affect the solution near tau0.
#
#Important: Note that the heating obtained from I+ - I-
#is only the diffuse radiation part of the heating. To get
#the full heating, one must also add in the heating due to
#direct beam attenuation. If you leave that out, you get
#a big spurious net cooling in the conversion layer at the
#top, where direct beam is being converted to diffuse radiation.
#
#
#Usage:
# First run the script to define the functions. It will
# print out the value of tauInf, which tells you the allowable
# range of tau for the local computation.
#
# Then do
# ctau,cp = doCalc(30.)
#
# where the argument of doCalc is the value of tau at the
# center of the region where you want the local flux calculation
# to take place (replace the value above with the value you want.)
# ctau is a curve giving fluxes as a function of optical depth.
# cp gives fluxes as a function of pressure. To plot them, do
# plot(ctau)
# plot(cp)
#
# For the figure included in the text, this calculation was done
# for several different tau covering the full domain, and the
# results for each tau were saved out into a file. The results were
# then composited onto a single graph (and prettied up) using a
# commercial plotting program.
#
#-----------------------------------------------------------------------
#
#Data on section of text which this script is associated with
Chapter = '5'
Section = '**'
Figure = '**'
#
#ToDo:
# *Double-check the direct-beam source term; put in
# code to deal with the optically thick case, where
# the conversion layer becomes too thin to resolve.
#
#
# *There is a problem with upward integration when
# atmosphere is too optically thin; needs fixing
#
# *Minor improvement:
# When the integration doesn't hit the boundary, set
# the boundary flux to the optically thick limit, rather
# than to zero. (Add this as an option; not mandatory)
#
from ClimateUtilities import *
import phys
import math
#Utility to invert a 2x2 matrix
def invert2(A):
det = A[0][0]*A[1][1] - A[0][1]*A[1][0]
return numpy.array( [ [A[1][1]/det,-A[0][1]/det],\
[-A[1][0]/det,A[0][0]/det] ])
#Planck functions in wavenumber space (cm**-1)
def Planck(wavenum,T):
return 100.*math.pi*phys.c*phys.B(100.*wavenum*phys.c,T)
#Function returning the derivatives of the
#upward and downward flux, for use with the
#integrator. I = [Iplus,Iminus].
#
def dIdtau(tau,I,source):
Iplus = I[0]
Iminus = I[1]
omg0val = omg0(pfun(tau))
Tval = T(pfun(tau))
gamma1 = gamma + gammaPrime*(1.-omg0val) #Assumes ghat=0
gamma2 = gamma - gammaPrime*(1.-omg0val)
gammaB = gamma1-gamma2
gammaPlus = gammaMinus = .5*omg0val #Assumes ghat = 0
Fsun = Lsun*Planck(wave,Tsun)/(phys.sigma*Tsun**4)
DirectBeam = Fsun*math.exp(-(tauInf-tau)/cosz)
dIpdtau = -gamma1*Iplus + gamma2*Iminus +\
source*(gammaB*Planck(wave,Tval) + gammaPlus*DirectBeam)
dImdtau = gamma1*Iminus - gamma2*Iplus -\
source*(gammaB*Planck(wave,Tval) + gammaMinus*DirectBeam)
return numpy.array([dIpdtau,dImdtau])
#This section is where we define the structure of the
#atmosphere, and generate tau(p) and omg0(p). The
#atmosphere we consider consists of a well-mixed
#gaseous absorber plus a scatterer with scattering
#cross section per unit mass chibar, having mass concentration
#q(p). The scattering particles do not absorb. If the
#scattering is Rayleigh scattering, the absorption is treated
#as part of the gaseous absorption coefficient.
#
g = 10.#Acceleration of gravity
pref = 1.e4 #Reference pressure for pressure broadening
alpha_g = 0. #Surface albedo
e_g = 1.-alpha_g #Surface emissivity
Tg = 270. #Surface temperature
Lsun = 1370. #Solar constant (for direct beam term)
Tsun = 6000. #Photosphere temperature of star
cosz = .5 #Cosine of zenith angle (for direct beam term)
ptop = 100.
pbot = 2.e5
dp = 100.
kappa0 = 1.5e-4#gaseous absorption coeff at p = pref
chibar = 1. #Scattering cross section per unit mass
qmax = .01 #Maximum cloud particle mass mixing ratio
pCloud = .5e5 #Pressure of cloud center
dpCloud = 1.e4 #Geometric thickness of cloud
def q(p):
return qmax*math.exp(-(p-pCloud)**2/(dpCloud**2))
def dtaudp(p,dummy):
return -(kappa0*(p/pref) + q(p)*chibar)/g
def omg0(p):
return q(p)*chibar/(kappa0*(p/pref)+ q(p)*chibar)
mtau = integrator(dtaudp,pbot,0.,-dp)
#Make up table of tau(p)
taupList = [0.]
pList = [pbot]
p = pbot
while p > ptop+1.e-4:
p,tau = mtau.next()
pList.append(p)
taupList.append(tau)
#Interpolation functions
pfun = interp(taupList,pList)
taufun = interp(pList,taupList)
#Temperature profile (specified as a function here,
#but you could make it an interpolation from a table
Rcp = 2./7.
Ts = 270.
def T(p):
return Ts*(p/pbot)**Rcp
wave = 600.
#Coefficients for the two-stream approx
gamma = 1.
gammaPrime = 1.
tauInf = taupList[-1]
print "tauInf = ",tauInf
#tau0 = taufun(p0) #Level in whose vicinity we want the solution
#p0 is the level where we want to calculate the profiles.
#Find the corresponding tau0 using tau0 = taufun(p0)
def doCalc(tau0):
big = math.exp(10.) #Criterion for when the exponentially
#growing term dominates and swamps accuracy.
E1 = numpy.array([1.,0.])
E2 = numpy.array([0.,1.])
Ez = numpy.array([0.,0.])
#Upward integration from tau1
#
#Loop to construct a particular solution and two independent
#homogeneous solutions.
# ToDo: To improve accuracy, use the
# optically thick solution as the initial condition
# for the fluxes (work when we are not too close
# to the pure scattering case; what's the best
# initialization to use in general?)
#
dtau = .1
#Initial conditions
mHom1 = integrator(dIdtau,tau0,E1,dtau)
mHom1.setParams(0.) #Homogeneous. No source
mHom2 = integrator(dIdtau,tau0,E2,dtau)
mHom2.setParams(0.) #Homogeneous. No source
IStartPart = Ez.copy()
mPart = integrator(dIdtau,tau0,IStartPart,dtau)
mPart.setParams(1.)
tau = tau0
tauListP = [tau0]
IListHom1P = [ E1 ]
IListHom2P = [ E2 ]
IListPartP = [ IStartPart ]
test = 0.
while (tau <= tauInf - dtau + 1.e-6) & (test < big):
tau,I = mHom1.next()
tauListP.append(tau)
IListHom1P.append(I.copy())
test1 = max(abs(I[0]),abs(I[1]))
#
tau,I = mHom2.next()
IListHom2P.append(I.copy())
test2 = max(abs(I[0]),abs(I[1]))
#
tau,I = mPart.next()
IListPartP.append(I.copy())
test3 = max(abs(I[0]),abs(I[1]))
test = max(test1,test2,test3)
nP = len(tauListP)
#Superpose the homogeneous solutions to
#produce one that has Ip = 0 at the top and
#one that has Im = 0 at the top
coeffA1 = IListHom2P[-1][0]
coeffA2 = -IListHom1P[-1][0]
coeffB1 = IListHom2P[-1][1]
coeffB2 = -IListHom1P[-1][1]
for i in range(nP):
IListHom1P[i],IListHom2P[i] = coeffA1*IListHom1P[i]+coeffA2*IListHom2P[i]\
,coeffB1*IListHom1P[i]+coeffB2*IListHom2P[i]
#Normalize
normA = IListHom1P[-1][1]
normB = IListHom2P[-1][0]
IListHom1P = [I/normA for I in IListHom1P]
IListHom2P = [I/normB for I in IListHom2P]
#Find correct superposition of particular and
#homogeneous solution to satisfy the upper b.c.
for i in range(nP):
IListPartP[i] += -IListHom1P[i]*IListPartP[-1][1]
#Upward integration is complete. Now we need
#to do it all over again for the downward integration,
#and then make the solutions continuous at tau0
#Loop to construct a particular solution and two independent
#homogeneous solutions.
# ToDo: To improve accuracy, use the
# optically thick solution as the initial condition
# for the fluxes (work when we are not too close
# to the pure scattering case; what's the best
# initialization to use in general?)
#
dtau = -.1
#Initial conditions
mHom1 = integrator(dIdtau,tau0,E1,dtau)
mHom1.setParams(0.) #Homogeneous. No source
mHom2 = integrator(dIdtau,tau0,E2,dtau)
mHom2.setParams(0.) #Homogeneous. No source
IStartPart = Ez.copy()
mPart = integrator(dIdtau,tau0,IStartPart,dtau)
mPart.setParams(1.)
tau = tau0
tauListM = [tau0]
IListHom1M = [ E1 ]
IListHom2M = [ E2 ]
IListPartM = [ IStartPart ]
test = 0.
while (tau >= -dtau) & (test < big):
tau,I = mHom1.next()
tauListM.append(tau)
IListHom1M.append(I.copy())
test1 = max(abs(I[0]),abs(I[1]))
#
tau,I = mHom2.next()
IListHom2M.append(I.copy())
test2 = max(abs(I[0]),abs(I[1]))
#
tau,I = mPart.next()
IListPartM.append(I.copy())
test3 = max(abs(I[0]),abs(I[1]))
test = max(test1,test2,test3)
nM = len(tauListM)
#Superpose the homogeneous solutions to
#produce one that has Ip = 0 at the bottom and
#one that has Im = 0 at the bottom
coeffA1 = IListHom2M[-1][0] - alpha_g*IListHom2M[-1][1]
coeffA2 = -(IListHom1M[-1][0] - alpha_g*IListHom1M[-1][1] )
coeffB1 = IListHom2M[-1][1]
coeffB2 = -IListHom1M[-1][1]
for i in range(nM):
IListHom1M[i],IListHom2M[i] = coeffA1*IListHom1M[i]+coeffA2*IListHom2M[i]\
,coeffB1*IListHom1M[i]+coeffB2*IListHom2M[i]
#Normalize
normA = IListHom1M[-1][1]
normB = IListHom2M[-1][0] - alpha_g*IListHom2M[-1][1]
IListHom1M = [I/normA for I in IListHom1M]
IListHom2M = [I/normB for I in IListHom2M]
#Find correct superposition of particular and
#homogeneous solution to satisfy the lower b.c.
#
#(Check this boundary condition
Fsun = Lsun*Planck(wave,Tsun)/(phys.sigma*Tsun**4)
DirectBeam = Fsun*math.exp(-tauInf/cosz) * cosz
if tauListM[-1] > abs(dtau):
hitBdd = 0. #Set to 1. if you hit the lower boundary
else:
hitBdd = 1.
for i in range(nM):
IListPartM[i] += \
-IListHom2M[i]*(IListPartM[-1][0] - alpha_g*IListPartM[-1][1])\
+hitBdd*IListHom2M[i]*(e_g*Planck(wave,Tg) + alpha_g*DirectBeam)
#Downward integration complete
#
#Now add in correct homogeneous
#solutions to make the solution continuous
#at tau0
A = numpy.array( [ [IListHom2P[0][0],IListHom1M[0][0]],\
[IListHom2P[0][1],IListHom1M[0][1]] ])
diff = IListPartP[0]- IListPartM[0]
coeffs = -numpy.dot(invert2(A),diff)
for i in range(nP):
IListPartP[i] += coeffs[0]*IListHom2P[i]
for i in range(nM):
IListPartM[i] += -coeffs[1]*IListHom1M[i]
#Suture together lower and upper solutions for output
tauListM.reverse()
tauList = tauListM + tauListP[1:] #Don't repeat tau0
IListPartM.reverse()
IList = IListPartM + IListPartP[1:] #Don't repeat tau0 results in list
#Compute heating, and transform back to p-space here if desired
p = [pfun(tau) for tau in tauList]
#Plot vs tau
ctau = Curve()
ctau.addCurve(tauList,'tau')
ctau.addCurve([I[0] for I in IList],'Iplus')
ctau.addCurve([I[1] for I in IList],'Iminus')
ctau.Xlabel = 'tau'
ctau.switchXY = True
#Plot vs p
cp = Curve()
cp.addCurve(p,'p')
cp.addCurve([I[0] for I in IList],'Iplus')
cp.addCurve([I[1] for I in IList],'Iminus')
cp.Xlabel = 'Pressure'
cp.switchXY = True
cp.reverseY = True
return ctau,cp