#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.