#Homebrew radiation code
#This version includes computation of
#the heating and flux due to stellar absorption. The
#thermal infrared calculation has not changed.
#Changes:
# 3/16/2009-- Replaced hardwired 260K ref temperature with Tref
#
# 3/17/2009-- Made the cutoff wavenumber for the longwave calc
# an adjustable parameter, rather than hardwired.
# (needed for high temperature calculations)
#Notes:
# --The temperature scaling for the continua assumes all continua
# are defined at a reference temperature of 300K; this introduces
# a slight error for the H2O continuum, which is defined relative
# to 296K, but it's small compared to other issues with the
# continuum. It should be fixed, eventually, though.
#
#Note regarding the stellar absorption: For simplicity
#the zenith angle is left at the same mean angle used
#in the thermal calculation. It ought to be separated
#out as a separate variable, so that one can average
#over time if desired. In other words, instead of using
#the same cosThetaBar for thermal and stellar, this should
#be made into an argument of the stellar transmission function.
#
#Note also that this code does
#not include the effects of Rayleigh scattering.
#This is an elaborated version of miniClimt.py,
#which incorporates temperature weighting and allows
#for an inhomogeneous path, and for distinguishing
#self from foreign broadening (needed for continuum).
#It has a very approximate method for dealing with the
#distinction beween the self and foreign line broadening,
#incorporated through a single band-independent absorption
#ratio, called SelfRat. This could easily be improved upon
#by reading in exponential sum tables for self and air broadening,
#and interpolating between the two (or using a band-dependent
#selfrat; if doing that, one might as well put in a band-dependent
#temperature scaling at the same time. For now, we've tried to keep
#things simple.
#
#Make sure to set selfrat consistently with the band table you
#are using. If you read in a band table that already was computed
#using self-broadening, then SelfRat should be set to 1 .
#
#This version doesn't include the semi-grey and Malkmus
#transmissions, which were in the simpler version
#only for pedagogical purposes.
#**ToDo:
# *Speed up optically thin case
# *Improve accuracy of optically thick case
# *Put in temperature weighting (only here,not in
# the simple version
# *Write an alternate version based on
# ODE integration (as prelude to scattering)
# *Separate the functions out from the work loop,
# so the radiation code can be used for multiple purposes
# *The continuum as written assumes a self-continuum. This
# is correct for water vapor, but not CO2. Need to find
# a way to generalize that. Perhaps define separate self
# and foreign continuum functions
#
import math
import string
#Things written for ClimateBook
import phys
from ClimateUtilities import *
#Mean slant path
cosThetaBar = .5
#Continuum functions (set continuum to desired function)
def NoContinuum(wave,T):
return 0.
#
#CO2 foreign continuum
def CO2ForeignContinuum(wave,T):
if (wave >= 25.) and (wave <= 550.):
kappa = math.exp(
-8.853 + 0.028534 *wave -0.00043194 *wave**2 + \
1.4349e-6* wave**3 -1.5539e-9* wave**4)
elif (wave >= 1050.) and (wave <= 1900.):
kappa = math.exp(
-537.09 +1.0886*wave -0.0007566 *wave**2 +1.8863e-7*wave**3\
-8.2635e-12 *wave**4)
else:
kappa = 0.
return (300./T)**1.7*kappa
#CO2 self continuum
def CO2SelfContinuum(wave,T):
if (wave >= 25.) and (wave <= 550.):
kappa = math.exp(
-8.853 + 0.028534 *wave -0.00043194 *wave**2 + \
1.4349e-6* wave**3 -1.5539e-9* wave**4)
elif (wave >= 1050.) and (wave <= 1900.):
kappa = math.exp(
-537.09 +1.0886*wave -0.0007566 *wave**2 +1.8863e-7*wave**3\
-8.2635e-12 *wave**4)
else:
kappa = 0.
return 1.3*(300./T)**1.7*kappa
#Water vapor self-continuum
def H2OSelfContinuum(wave,T):
if (wave>= 500.) and (wave<= 1400.):
kappaC = math.exp(12.167-0.050898*wave +8.3207e-05*wave*wave
-7.0748e-08*wave**3 + 2.3261e-11*wave**4)
elif (wave>=2100.) and (wave <= 3000.):
x = wave - 2500.
kappaC = math.exp(-6.0055 + -0.0021363*x + 6.4723e-07*x**2 +
-1.493e-08*x**3 + 2.5621e-11*x**4 + 7.328e-14*x**5)
else:
kappaC = 0.
#Temperature dependence:
return kappaC*(296./T)**4.25
H2OForeignContinuum = NoContinuum
ForeignContinuum = NoContinuum #Default
SelfContinuum = NoContinuum
#Planck functions in wavenumber space (cm**-1)
def Planck(wavenum,T):
return 100.*math.pi*phys.c*phys.B(100.*wavenum*phys.c,T)
#Temperature derivative of Planck function
#Not currently used, but might be useful for implementing the
#optically thick limit to improve accuracy
#def dPlanck(wavenum,T):
# nu = 100.*wavenum*phys.c
# u = min(phys.h*nu/(phys.k*T),500.) #To prevent overflow
# dB = (2.*phys.h*nu**3/phys.c**2)*(phys.h*nu/phys.k)*(1/T**2)*math.exp(u)/(math.exp(u)-1.)**2
# return math.pi*100.*phys.c*dB
#Routine to read in band data from table
#by wavenumber.
def loadExpSumTable(fileName):
c = readTable(fileName)
bandList = []
headers = c.listVariables()
binHeaders = [h for h in headers if 'dH' in h]
logKappaHeaders = [h for h in headers if 'logKappa' in h]
nHeaders = len(binHeaders)
for i in range(nHeaders):
h = logKappaHeaders[i]
h1 = binHeaders[i]
waves = [string.atof(h.split('.')[1]),string.atof(h.split('.')[2])]
bandParams = Dummy()
bandParams.nu1 = string.atof(h.split('.')[1])
bandParams.nu2 = string.atof(h.split('.')[2])
bandParams.bandType = 'ExponentialSum'
bandParams.kappa = Numeric.exp(c[h])
bandParams.dH = c[h1]
#
bandList.append(bandParams)
return bandList
#Band-averaged transmission function
#To keep track of self-collision vs. foreign collisions,
#we need two different pressure-weighted paths. "path"
#is the foreign-collision path, whereas "pathq" is the
#self-collision path. If we were getting really fancy, we
#might need more paths, since the temperature dependence
#of the continuum is different from that of the line absorption
def TransEsums(path,pathq,bandParams):
#Speed this up by dropping small terms from the sum.
#It's easy to find which are smalll, since kappa is
#sorted in increasing order.
#
#Handle the continuum (Reference Temperature curently hard-wired to 300K)
#Temperature scaling is currently handled in the path calculation
#Fix this, to generalize (use separate self and foreign contin fn)
nu = (bandParams.nu1 + bandParams.nu2)/2.
continOpacity = \
ForeignContinuum(nu,300.)*path + SelfContinuum(nu,300.)*pathq
contin = math.exp(-continOpacity)
#
#The following statement cuts off the exponential sum
#when the argument of the exponential gets too large. It
#prevents underflow and also speeds things up in the optically
#thick case
pathTot = path + SelfRat*pathq #**This is where we handle the enhancement of
#self-collisions. Make this band-dependent
#for greater accuracy.
imax = Numeric.searchsorted(bandParams.kappa*pathTot,15.)
return contin*sum(bandParams.dH[:imax] * Numeric.exp(-pathTot*bandParams.kappa[:imax]))
#Note: q is now an array
#In this version, we compute the path in the transmission
#loop, which is more efficient than doing the whole integral
#from scratch in the transmission function routine.
#Note:All the path computations need to be cleaned up
def radCompBand(p,T,Tg,q,bandParams):
n = len(p)
Delta = bandParams.nu2-bandParams.nu1
wave = (bandParams.nu2+bandParams.nu1)/2.
#
#Mass ratio for computing molar concentration
# (used to compute proportion of self-collisions)
massrat = GHG.MolecularWeight/BackgroundGas.MolecularWeight
#Temporary arrays for flux computation
Ip = Numeric.zeros(n,Numeric.Float)
Im = Numeric.zeros(n,Numeric.Float)
#-------------------------------------
#Upward flux
for i in range(n):
Sum = 0.
path = 0.
pathq=0.
for j in range(i,n-1):
#Note: we can use dB directly instead of
# dB/dT * dT/dp * dp
dB = Planck(wave,T[j+1])-Planck(wave,T[j])
moleCon = q[j]/(q[j] + (1-q[j])*massrat)
#Temperature scaling------------
#As currently written we
#have one Tstar for all bands. Further, if there
#is a continuum, the scaling is taken from the continuum,
#whereas it would be better to introduce a separate
#path for the continuum
if ForeignContinuum(wave,300.) > 1.e-6:
wTF = ForeignContinuum(wave,T[j])/ForeignContinuum(wave,300.)
else:
wTF = math.exp(-Tstar*(1./T[j] - 1./Tref))
if SelfContinuum(wave,300.) > 1.e-6:
wTS = SelfContinuum(wave,T[j])/SelfContinuum(wave,300.)
else:
wTS = math.exp(-Tstar*(1./T[j] - 1./Tref))
#---------------------------------
#**Move this to midpoint for accuracy
path = path + wTF*((1-moleCon)*p[j]/pref)*q[j]*abs(p[j+1]-p[j])/(g*cosThetaBar)
#Compute the self-broadened pressure path
pathq = pathq + wTS*moleCon*(p[j]/pref)*q[j]*abs(p[j+1]-p[j])/(g*cosThetaBar)
Trans1 = Trans(path,pathq,bandParams)
#This will be inaccurate for j=i when an individual
#layer becomes optically thick
Sum += Trans1*dB*Delta
if Trans1 < 1.e-4:
break
#Add in boundary term
BddTerm = (Planck(wave,Tg)-Planck(wave,T[-1]))*Trans(path,pathq,bandParams)
Ip[i] = Sum + (Planck(wave,T[i]) + BddTerm)*Delta
#Downward flux
for i in range(n):
Sum = 0.
path = 0.
pathq= 0.
for j in range(i,0,-1):
#Note: we can use dB directly instead of
# dB/dT * dT/dp * dp
dB = Planck(wave,T[j-1])- Planck(wave,T[j])
moleCon = q[j]/(q[j] + (1-q[j])*massrat)
#Temperature scaling------------
#As currently written we
#have one Tstar for all bands. Further, if there
#is a continuum, the scaling is taken from the continuum,
#whereas it would be better to introduce a separate
#path for the continuum
if ForeignContinuum(wave,300.) > 1.e-6:
wTF = ForeignContinuum(wave,T[j])/ForeignContinuum(wave,300.)
else:
wTF = math.exp(-Tstar*(1./T[j] - 1./Tref))
if SelfContinuum(wave,300.) > 1.e-6:
wTS = SelfContinuum(wave,T[j])/SelfContinuum(wave,300.)
else:
wTS = math.exp(-Tstar*(1./T[j] - 1./Tref))
#---------------------------------
#**Fix this; move to midpoint
path = path + wTF*((1-moleCon)*p[j]/pref)*q[j]*abs(p[j-1]-p[j])/(g*cosThetaBar)
#Compute the self-broadened pressure path
pathq = pathq + wTS*moleCon*(p[j]/pref)*q[j]*abs(p[j-1]-p[j])/(g*cosThetaBar)
Trans1 = Trans(path,pathq,bandParams)
Sum += Trans1*dB*Delta
if Trans1 < 1.e-4:
break
#Add in boundary term
BddTerm = (-Planck(wave,T[0]))*Trans(path,pathq,bandParams)
Im[i] = Sum + (Planck(wave,T[i]) + BddTerm)*Delta
flux = Ip-Im
heat = Numeric.zeros(n,Numeric.Float)
#
#Compute heat by taking the gradient of the flux.
#The second order centered difference below allows
#the computation to be done accurately even if
#the pressure level spacing is non-uniform.
#
#Returns heating rate as W/kg.
#Divide by specific heat to get K/s
#
delPlus = p[2:] - p[1:-1] #Should wind up as a [1:-1] array
delMinus = p[1:-1]-p[:-2] #likewise
A = (delMinus/delPlus)/(delPlus+delMinus)
B = 1./delMinus - 1./delPlus
C= (delPlus/delMinus)/(delPlus+delMinus)
heat[1:-1] = A*flux[2:] + B*flux[1:-1] - C*flux[:-2]
heat = g*heat #Convert to W/kg
heat[0]=heat[1]
heat[-1] = heat[-2]
return flux,heat
def stellarRadCompBand(p,T,Tg,q,bandParams):
n = len(p)
Delta = bandParams.nu2-bandParams.nu1
wave = (bandParams.nu2+bandParams.nu1)/2.
#
#Mass ratio for computing molar concentration
# (used to compute proportion of self-collisions)
massrat = GHG.MolecularWeight/BackgroundGas.MolecularWeight
#Now compute the transmission function between each
#level and the top.
OrbitFact = (Lstellar/4.)/(phys.sigma*Tstellar**4)
swFlux = OrbitFact*Delta*Planck(wave,Tstellar)
flux = Numeric.zeros(n,Numeric.Float)
flux[0] = -swFlux #Downward flux is negative
#
path = 0.
pathq= 0.
for j in range(1,n):
moleCon = q[j]/(q[j] + (1-q[j])*massrat)
#Temperature scaling------------
#As currently written we
#have one Tstar for all bands. Further, if there
#is a continuum, the scaling is taken from the continuum,
#whereas it would be better to introduce a separate
#path for the continuum
if ForeignContinuum(wave,300.) > 1.e-6:
wTF = ForeignContinuum(wave,T[j])/ForeignContinuum(wave,300.)
else:
wTF = math.exp(-Tstar*(1./T[j] - 1./Tref))
if SelfContinuum(wave,300.) > 1.e-6:
wTS = SelfContinuum(wave,T[j])/SelfContinuum(wave,300.)
else:
wTS = math.exp(-Tstar*(1./T[j] - 1./Tref))
#---------------------------------
#**Fix this; move to midpoint
path = path + wTF*((1-moleCon)*p[j]/pref)*q[j]*abs(p[j-1]-p[j])/(g*cosThetaBar)
#Compute the self-broadened pressure path
pathq = pathq + wTS*moleCon*(p[j]/pref)*q[j]*abs(p[j-1]-p[j])/(g*cosThetaBar)
Trans1 = Trans(path,pathq,bandParams)
flux[j] = -swFlux*Trans1
#Flux computed. Now compute the heating
heat = Numeric.zeros(n,Numeric.Float)
#
#Compute heat by taking the gradient of the flux.
#The second order centered difference below allows
#the computation to be done accurately even if
#the pressure level spacing is non-uniform.
#
#Returns heating rate as W/kg.
#Divide by specific heat to get K/s
#
delPlus = p[2:] - p[1:-1] #Should wind up as a [1:-1] array
delMinus = p[1:-1]-p[:-2] #likewise
A = (delMinus/delPlus)/(delPlus+delMinus)
B = 1./delMinus - 1./delPlus
C= (delPlus/delMinus)/(delPlus+delMinus)
heat[1:-1] = A*flux[2:] + B*flux[1:-1] - C*flux[:-2]
heat = g*heat #Convert to W/kg
heat[0]=heat[1]
heat[-1] = heat[-2]
return flux,heat
#Sum of flux and heating over all bands
def radCompLW(p,T,Tg,q):
n = len(p)
flux = Numeric.zeros(n,Numeric.Float)
heat = Numeric.zeros(n,Numeric.Float)
for band in bandData:
if band.nu1 < LWCutoff:
fluxBand,heatBand = radCompBand(p,T,Tg,q,band)
flux += fluxBand
heat += heatBand
return flux,heat
def radCompStellar(p,T,Tg,q):
n = len(p)
flux = Numeric.zeros(n,Numeric.Float)
heat = Numeric.zeros(n,Numeric.Float)
for band in bandData:
#Shortwave (stellar) contribution
fluxBand,heatBand = stellarRadCompBand(p,T,Tg,q,band)
flux += fluxBand
heat += heatBand
return flux,heat
#----------------End of radiation model------------------
#
#Misc. utility functions
#
#OLR function for dry air adiabat
#Although our basic radiation calculation uses mass
#concentrations, this one takes molar concentrations
#of the greenhouse gas (in ppmv) as input, to make the
#units more familiar. It also increases surface pressure
#to account for the additional mass of greenhouse gas.
#The surface pressure adjustment makes little difference out
#to 10% CO2, but as one goes to 20% and above it starts to become
#important
#
#Do we want to put in an isothermal stratosphere, to make
#this function more parallel to the ccmrad function?
def OLRDryAir(psAir,Tg,GHGmolarcon):
ptop = 100.
ps = psAir/(1.-1.e-6*GHGmolarcon)
Mbar = 1.e-6*GHGmolarcon*GHG.MolecularWeight + (1.-1.e-6*GHGmolarcon)*BackgroundGas.MolecularWeight
p = setpLin(ps,ptop,40)
T = Tg*(p/ps)**(2./7.) #Ignores effect of GHG on R/cp
q = 1.e-6*GHGmolarcon*GHG.MolecularWeight/Mbar
q = q*Numeric.ones(40,Numeric.Float)
flux,heat = radComp(p,T,Tg,q)
#Add in contribution of part of spectrum outside the band table
nu1 = bandData[0].nu1
nu2 = bandData[-1].nu2
m = romberg(Planck) #Evaluates definite integral
flux1 = m([.001,nu1],Tg) + m([nu2,50000.],Tg)
return flux[0] + flux1
#Functions for setting up pressure arrays
#
#Linear pressure array
def setpLin(ps,ptop,n):
p = [ptop + i*(ps-ptop)/(n-1) for i in range(n)]
return Numeric.array(p)
#Log pressure array (more resolution in stratosphere)
def setpLog(ps,ptop,n):
rat = (ps/ptop)**(1./(n-1))
p = [ptop*rat**i for i in range(n)]
return Numeric.array(p)
#Linear array with extra resolution near ground
def setpExtraGroundRes(ps,ptop,n):
n1 = 2*n/3
n2 = n-n1
p = [ptop + ((.9*ps-ptop)*i)/(n1-1) for i in range(n1)]
p1 = [.9*ps + .1*ps*i/(n2) for i in range(1,n2+1)]
p = p + p1 #Concatenation, not addition!
return Numeric.array(p)
#-----------------------------------------------------------------------
#Remember to set gravity for your planet! It can be over-ridded
g = 9.8
#Read in band data (Moved to the calling program)
#
#Default data
BackgroundGas = phys.air #The transparent background gas
GHG =phys.CO2 #The greenhouse gas
#
pref = 1.e4 #Reference pressure at which band data is given
Tref = 260. #Reference temperature at which band data is given
#Used for in-band scaling. Continuum has separate
#temperature scaling.
#**ToDo: Eliminate these defaults and force user to choose
Trans = TransEsums
SelfRat = 1. #Set mean ratio of self to foreign broadening for the
#gas in use
#Set the continuum, if there is one
#SelfContinuum = CO2SelfContinuum
#ForeignContinuum = CO2ForeignContinuum
#To speed up the calculation, when doing both
#stellar absorption and longwave flux together,
#you can set a cutoff wavenumber for the longwave
#flux. That can be helpful, since the atmosphere
#is usually a lot cooler than the photosphere and
#therefore doesn't emit much at the shorter
#waves that are important for the incoming stellar
#absorption. Use with caution, though: if the
#atmosphere gets very hot, then you need the shorter
#wave emission
LWCutoff = 1.e6 #Wavenumber cutoff in 1/cm for longwave calc
#Set the photosphere temperature of the star
#(for shortwave radiation calculation
Tstellar = 3500.
#Set the stellar constant at the orbit
Lstellar = 4.*700.