#Description: #============================================================== #Script to solve the two-stream equations with scattering, #by numerical integration as an ODE problem. This version #uses exponential sums to handle real-gas absorption properties, #and sums fluxes over bands. # #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. # #The function doCalc(tau) is essentially the same as in the #script TwoStreamScatter.py. The script here calls doCalc(...) #repetitively for each absorption coefficient in the exponential #sum table, and then makes the weighted sums of the results. #It does this for each wavenumber in the table, and plots #the result as a spectrum. This script is set up to #just compute the OLR spectrum, so it only calls doCalc for #a tau near the top of the atmosphere, since that is sufficient #to compute the OLR. It could easily be modified to compute #profiles of the fluxes, if those are needed. # #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. # #------------------------------------------------------------ #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. # # # *Fix problem with upward integration when # atmosphere is too optically thin; # # *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,planets import math #Path to the workbook datasets datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/' #Path to the exponential sum tables EsumPath = datapath + 'Chapter4Data/ExpSumTables/' #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] ]) #------------Utility functions copied over from miniClimt---------- #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 = numpy.exp(c[h]) bandParams.dH = c[h1] # bandList.append(bandParams) return bandList #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 #----------------End copies from miniClimt------------------------- #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]. # # #ToDo: Need params to pass kappa0 and wave, not just source def dIdtau(tau,I,params): source = params.source kappa0 = params.kappa0 wave = params.wave Iplus = I[0] Iminus = I[1] omg0val = omg0(pfun(tau),kappa0) #Rewrite this so we call omg0Tau(tau); get rid of kappa0? 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]) #Coefficients for the two-stream approx gamma = 1. gammaPrime = 1. #Yes, I know we can integrate this analytically, #but it's providede here as a function in case we #want to incorporate temperature scaling or a more #complicated form of pressure scaling def dtaudpGas(p,dummy): return -(p/pref)/g def dtaudpParticle(p,dummy): return -q(p)/g def omg0(p,kappa0): return q(p)*chibar/(kappa0*(p/pref)+ q(p)*chibar) #Sets up arrays needed to do tau(p) and p(tau) #interpolation. We implement this mapping as #a sum of optical thickness from gaseous absorption #and particulate absorption/scattering so the integral #doesn't have to be done all over again for each #kappa0. For now, the cross-section chibar is #assumed independent of p, but that could be changed. #If chibar is a function of both p and wavenumber in a #non-separable way, then some additional recomputation would #have to be done inside the doCalc function, for each wavenumber #(though not for each kappa0). def setupTauP(): m1 = integrator(dtaudpGas,pbot,0.,-dp) m2 = integrator(dtaudpParticle,pbot,0.,-dp) pList = [pbot] taupGasList = [0.] taupParticleList = [0.] p = pbot while p > ptop+1.e-4: p,tauGas = m1.next() pList.append(p) taupGasList.append(tauGas) p,tauParticle = m2.next() taupParticleList.append(tauParticle) return numpy.array(pList),numpy.array(taupGasList),numpy.array(taupParticleList) def pTauFunctions(kappa0,chibar): taupList = kappa0*taupGas + chibar*taupParticle #Interpolation functions pfun1 = interp(taupList,pList) taufun1 = interp(pList,taupList) tauInf1 = taupList[-1] #print "tauInf = ",tauInf1 return tauInf1,pfun1,taufun1 #p0 is the level where we want to calculate the profiles. #Find the corresponding tau0 using tau0 = taufun(p0) #kappa0 is the absorption coefficient at the reference pressure p0 # #Optimization note: If an entire profile of flux were needed, #as for radiative-convective equilibrium calculations, it would #be more efficient to modify doCalc to loop over the full #range of p0 needed to cover the atmosphere, and interpolate #each chunk of data to a standard pressure grid, then return #the whole profile. That allows data from each subintegration #to be re-used. The main speedup would be for the more optically #thin spectral regions, where it takes only a few subintegrations #to cover the whole atmosphere. def doCalc(p0,kappa0=1.5e-4,wave=600.): global tauInf,pfun,taufun tauInf,pfun,taufun = pTauFunctions(kappa0,chibar) tau0 = taufun(p0) paramsHom1 = Dummy() paramsHom2 = Dummy() paramsPart = Dummy() paramsHom1.kappa0 = paramsHom2.kappa0 = paramsPart.kappa0 = kappa0 paramsHom1.wave = paramsHom2.wave = paramsPart.wave = wave 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) paramsHom1.source = 0. mHom1.setParams(paramsHom1) #Homogeneous. No source mHom2 = integrator(dIdtau,tau0,E2,dtau) paramsHom2.source = 0. mHom2.setParams(paramsHom2) #Homogeneous. No source IStartPart = Ez.copy() mPart = integrator(dIdtau,tau0,IStartPart,dtau) paramsPart.source = 1. mPart.setParams(paramsPart) 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) paramsHom1.source = 0. mHom1.setParams(paramsHom1) #Homogeneous. No source mHom2 = integrator(dIdtau,tau0,E2,dtau) paramsHom2.source = 0. mHom2.setParams(paramsHom2) #Homogeneous. No source IStartPart = Ez.copy() mPart = integrator(dIdtau,tau0,IStartPart,dtau) paramsPart.source = 1. mPart.setParams(paramsPart) 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] #Compute the heating here, and return, if desired #Return the fluxes at p0 return IListPartP[0][0],IListPartP[0][1] ## #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 #Optimization note: It would be more efficient to hand #the integrator in doCalc a list of kappa0 for which the #calculation is to be done, since then vector arithmetic #instead of looping could be used. That would also make the #calculation easy to run in parallel on beowulf clusters. def radCompBand(p0,bandParams): Delta = bandParams.nu2-bandParams.nu1 wave = (bandParams.nu2+bandParams.nu1)/2. n = len(bandParams.kappa) Iplus = 0. Iminus = 0. for i in range(n): #Skip contribution from transparent gas #(Because of the way the exp sum tables are constructed, #this still allows the continuum to be included) if bandParams.dH[i] < 1.e-20: break kappa0 = bandParams.kappa[i] + Continuum(wave,300.) #Add continuum Ip,Im = doCalc(p0,kappa0,wave) Iplus += Ip*bandParams.dH[i] Iminus += Im*bandParams.dH[i] return Iplus*Delta,Iminus*Delta #Load in the exponential sum table bandData = loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.self.data') pref = 1.e4 #Reference pressure for pressure broadening Continuum = CO2Continuum #Set the continuum ptop = 100. pbot = 2.e5 dp = 100. g = planets.Mars.g#Acceleration of gravity alpha_g = 0. #Surface albedo (This should be wavenumber-dependent) e_g = 1.-alpha_g #Surface emissivity Tg = 275. #Surface temperature Lsun = 0. #1370. #Solar constant (for direct beam term) (Zeroed out for thermal IR calc) Tsun = 6000. #Photosphere temperature of star cosz = .5 #Cosine of zenith angle (for direct beam term) #--------------Define atmospheric structure. # (this could also be done by interpolation from a list of values) # #**ToDo: Rearrange setting of atmospheric profile so # that it's easier to turn this into an importable radiation # module like miniClimt. That means changing the computation # of setupTauP so it is redone for each new profile, and also # making things like qmax arguments rather than globals def q(p): return qmax*math.exp(-(p-pCloud)**2/(dpCloud**2)) #Temperature profile (specified as a function here, #but you could make it an interpolation from a table Rcp = phys.CO2.Rcp Ts = Tg #Set needed thermodynamic constants stuff = phys.CO2 L = stuff.L_sublimation # Latent heat cp = stuff.cp #Specific heat R = stuff.R #Gas constant for the atmosphere Tr = stuff.TriplePointT # Reference temperature psat_ref = stuff.TriplePointP # Saturation vapor pressure at #reference temperature #Function that solves for temperature at which saturation #occurs, for a given pressure. This uses the simplified #form of Clausius-Clapeyron based on the assumption that #the latent heat is constant. def Tsat(p): return(1./(1./Tr - (R/L)*math.log(p/psat_ref))) def T(p): Tdry = Ts*(p/pbot)**Rcp return max(Tdry,Tsat(p)) #Includes CO2 condensation effect #-------------- chibar = 100. #Scattering cross section per unit mass qmax = 1.e-5#1.e-5 pCloud = .5e5 #Pressure of cloud center dpCloud = 1.e4 #Geometric thickness of cloud #Set up tau(p) tables. Redo if you change qmax, etc pList,taupGas,taupParticle = setupTauP() #Sum over bands OLR = 0. nuList = [] OLRSpecList = [] for band in bandData: Ip,Im = radCompBand(100.,band) OLR += Ip nu = .5*(band.nu1+band.nu2) Delta = band.nu2-band.nu1 nuList.append(nu) OLRSpecList.append(Ip/Delta) print nu,Ip,Im,OLR c = Curve() c.addCurve(nuList,'wavenumber') c.addCurve(OLRSpecList,'OLR') c.addCurve([Planck(wave,Tg) for wave in nuList],'Planck') plot(c)