#A general radiation code to compute radiation fluxes and heating rates #from tabulated temperature/pressure data, for a single well-mixed #greenhouse gas. Note that this code can be used for saturated #pure water-vapor atmopheres, since the mixing ratio is 1 everywhere. #This script doesn't do temperature scaling of band data, however. #The continuum, if present, is computed for a fixed temperature #of 300K, though that can be changed by editing the value used #in the function TransEsums(...) . # #Change history: # 3/17/2009-- Made the cutoff wavenumber for the longwave calc # an adjustable parameter, rather than hardwired. # (needed for high temperature calculations) # #**FOR COMPATIBILITY: Previous calls to radComp in Chapter 4 scripts #need to be changed to radCompLW #This is the simplest "Trainer version" of the code. #A more comprehensive (and slightly more complicated) #version is found in miniClimtFancy.py #**ToDo: # *Speed up optically thin case # *Improve accuracy of optically thick case # # *The path does not currently take into account # self vs foreign broadening. The code assumes either # all foreign (as in dilute CO2 in air) or all self # (as in pure CO2), and whether we are doing self or # foreign is determined just by the exponential sum table. # Should we leave it this way to keep it simple, or generalize # the treatment? # # *TransEsums doesn't check for and handle transparent # bands. This needs to be fixed. (do this in miniClimtFancy, too) import math #import Numeric #Numeric is imported #as part of ClimateUtilities, for compatibity with numpy 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. def CO2Continuum(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. #Temperature dependence #kappa = kappa * ... return 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 #Planck functions in wavenumber space (cm**-1) def Planck(wavenum,T): return 100.*math.pi*phys.c*phys.B(100.*wavenum*phys.c,T) #Routine to read in band data from table 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 functions. #The exponential sum form, TransEsums, is the one we use #for real-gas calculations, but the other two are included #for the sake of comparison. #Semigrey ("Oobleck") model def TransOobleck(p1,p2,q,bandParams): pref = 1.e5 pbar = .5*(p1+p2) #pbar = pref #No pressure broadening path = q*(pbar/pref)*abs(p1-p2)/(g*cosThetaBar) return math.exp(-path*bandParams.kappa) #Malkmus model def TransMalkmus(p1,p2,q,bandParams): #Uncomment following for Strong line variant: #return math.exp(-N*math.sqrt(q*kappa*(.5*(p1+p2)/1.e5)*abs(p1-p2)/g)/Delta) # #With this, singularity is coming because we take the strong line approx #all the way to delP = 0. This give singularity in deriv of sqrt. #Should cut over to weak line (as in Malkmus) for thin layers. #Physical implication: Hard to get a temperature discontinuity #at the surface if lines are strong all the way to near the surface. # #Malkmus model: pref = 1.e5 R = bandParams.R S = bandParams.S Delta = bandParams.nu2-bandParams.nu1 # Rp = R*math.sqrt(((p1+p2)/2.)/pref) path = q*abs(p1-p2)/(g*cosThetaBar) W = 2.*(Rp*Rp/S)*(math.sqrt(1.+ path*S*S/(Rp*Rp)) - 1.) return math.exp(-W/Delta) def TransEsums(p1,p2,q,bandParams): # #Compute the path. # pbar = .5*(p1+p2) path = q*(pbar/pref)*abs(p1-p2)/(g*cosThetaBar) # #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 imax = Numeric.searchsorted(bandParams.kappa*path,15.) wave = (bandParams.nu1+bandParams.nu2)/2. ctrans = math.exp(-path*Continuum(wave,300.)) return ctrans*sum(bandParams.dH[:imax] * Numeric.exp(-path*bandParams.kappa[:imax])) #------------------------------------------------------ #Band parameter data for one-band Malkmus and for #one-band and two-band semigrey models #**ToDo: Move this out of here, into a separate file (e.g. #in PureRadEq, where it is used # malkmusBandParams = Dummy() malkmusBandParams.nu1 = 650. malkmusBandParams.nu2 = 700. malkmusBandParams.R =817.1 malkmusBandParams.S = 8.597e4 semiGreyBandParamsM = Dummy() semiGreyBandParamsM.nu1 = 600. semiGreyBandParamsM.nu2 = 650. semiGreyBandParamsM.kappa = .1 semiGreyBandParams = Dummy() semiGreyBandParams.nu1 = 650. semiGreyBandParams.nu2 = 700. semiGreyBandParams.kappa = 1. semiGreyBandParamsP = Dummy() semiGreyBandParamsP.nu1 = 700. semiGreyBandParamsP.nu2 = 750. semiGreyBandParamsP.kappa = .1 #-------------------------------------------------------- #Compute the thermal flux and heating def radCompBand(p,T,Tg,q,bandParams): n = len(p) ps = p[-1] Delta = bandParams.nu2-bandParams.nu1 wave = (bandParams.nu2+bandParams.nu1)/2. #Temporary arrays for use in flux computation Ip = Numeric.zeros(n,Numeric.Float) Im = Numeric.zeros(n,Numeric.Float) #Upward flux for i in range(n): Sum = 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]) Trans1 = Trans(p[i],(p[j]+p[j+1])/2.,q,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(p[i],ps,q,bandParams) Ip[i] = Sum + (Planck(wave,T[i]) + BddTerm)*Delta #Downward flux for i in range(n): Sum = 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]) Trans1 = Trans(p[i],(p[j]+p[j-1])/2.,q,bandParams) Sum += Trans1*dB*Delta if Trans1 < 1.e-4: break #Add in boundary term BddTerm = (-Planck(wave,T[0]))*Trans(p[i],0.,q,bandParams) Im[i] = Sum + (Planck(wave,T[i]) + BddTerm)*Delta flux = Ip-Im # # #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. # #Heating rate is returned 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 #Compute the stellar flux and heating def stellarRadCompBand(p,T,Tg,q,bandParams): n = len(p) ps = p[-1] Delta = bandParams.nu2-bandParams.nu1 wave = (bandParams.nu2+bandParams.nu1)/2. #Contribution from shortwave flux OrbitFact = (Lstellar/4.)/(phys.sigma*Tstellar**4) swFlux = OrbitFact*Delta*Planck(wave,Tstellar)*Numeric.array([Trans(pp,0.,q,bandParams) for pp in p]) flux = -swFlux # #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. # #Heating rate is returned 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: fluxBandLW,heatBandLW = radCompBand(p,T,Tg,q,band) flux += fluxBandLW heat += heatBandLW 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 the radiation model. ------------------ #The functions that follow are utility functions that make #it easier to use the basic radComp and radCompBand functions #to do radiation calculations. #------------------------------------------------------------- #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 q = 1.e-6*GHGmolarcon*GHG.MolecularWeight/Mbar p = setpLin(ps,ptop,40) T = Tg*(p/ps)**(2./7.) #Ignores effect of GHG on R/cp flux,heat = radCompLW(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) #-------------End of Function definitions --------------- #Set default parameters for radiation computation. These can #be over-ridden in the calling program. For example, if #you imported this module using "import miniClimt as homebrew" #then to change the gravity you would use homebrew.g = 3.5 Trans = TransEsums #Default is to use exponential sums #Trans = TransMalkmus #Trans = TransOobleck #Set the gravity g = 9.8 #Set the photosphere temperature of the star #(for shortwave radiation calculation Tstellar = 3500. #Set the stellar constant at the orbit Lstellar = 4.*700. #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 #--->The following options should be set in the calling program # They customize the radiation calculation for the gases # you are interested in # # #Read in band data #This should be done in the calling program, so you are sure #which gases you are using. Some examples are below. # #Default Data for principal band of CO2. Good up to about 10% CO2 in 1bar air. #bandData = loadExpSumTable('CO2TablePrinBand.260K.100mb.air.data') # #Data for the full CO2 spectrum out to 2300 cm**-1. Use for Venus #and Early Mars calculations #bandData = loadExpSumTable('CO2TableAllWave.260K.100mb.self.data') #bandData = bandData[:-9] #Eliminate "transparent region" from table # #Data for two-band Oobleck. To do one band, just make #a list having only semiGreyBandParams alone. #bandData = [semiGreyBandParamsM,semiGreyBandParams,semiGreyBandParamsP] # #pref = 1.e4 #Default reference pressure at which band data is given #THIS SHOULD BE SET IN THE CALLING PROGRAM, FOR COMPATIBILITY WITH BANDDATA # #Set the continuuum. Use NoContinuum if there is none #THIS SHOULD BE DONE IN THE CALLING PROGRAM #Set the continuum, if there is one #Continuum = NoContinuum #Continuum = CO2Continuum # #Say what the gases are. For the pure GHG case, the background #gas is just ignored, but it has to be specified anyway. These are #used to get needed theromdynamic data. #THIS SHOULD BE DONE IN THE CALLING PROGRAM #BackgroundGas = phys.air #The transparent background gas #GHG =phys.CO2 #The greenhouse gas #