#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