#------------------------------------------------------------------------------ #This script computes radiative-convective equilibrium by timestepping. #It has been modified to use the ccmrad radiation routines # # #Warning: Because of the way the surface temperature stepping is done #to achieve TOA balance , when you turn off the stellar heat computation #using the doStellarAbs flag, it does not affect the radiative flux #(only the heating rate, which is what's used in the computation). #If you want to do a computation with an explicit surface budget #incorporated, you need to modify the way the time-stepping #of surface temperature is done, and also exclude the stellar #absorption term from the flux. # # #ToDo: # *Merge this ccm rad/conv model with the rad/conv model # using the homebrew radiation, so that the radiation model # can be just changed using a switch. # ->Need radcomp in ccmradFunctions. Return heating in W/kg # #Change Log: # # *1/28/2012: Updated the radiation calls so it works # with the current version of climt_lite. Note that # despite the information in the help() documentation, # climt_lite now returns heating rates in K/day, not K/s # There have been a number of changes to climt since the # version used to make the figure in the book, including resolution, # default albedos, units, and possibly there may have been # a misprint in the book concerning the ozone concentration # used. With zero surface albedo, the values of insolation (solin) # and ozone reproduce the results in the book very closely. # # Also updated obsolete calls to Numeric. Replaced with numpy # #------------------------------------------------------------------------------ #Data on section of text which this script is associated with Chapter = '5' Section = 'section:RealGasStellarAbs' Figure = 'fig:OzoneResults' # #Instructions: To reproduce the figure, you need to run this script three #times. # (1) Once with ozone zeroed out but shortwave heating included. # (doStellarAbs = True) This includes the near-IR heating # due to CO2, and a little bit of UV heating from O2 # # (2) Once with ozone set to the specified profile, but with # shortwave heating turned off (doStellarAbs = False). This # also eliminates the CO2 and O2 absorption, but there's no # way to turn off the heating due to just one species. # # (3) Once as in (2), but with the shortwave heatting on (doStellarAbs = True) # #Save the results from each run, and combine in a single graph. #Any one of the above runs also saves the shortwave flux data needed to make #the left panel of the figure. import climt_lite as climt #ccm radiation model (Replace with ccmradFunctions?) from ClimateUtilities import * import math,phys import planets #--------------------Set radmodel options------------------- #---Instantiate the radiation model--- r = climt.radiation() n = r.nlev #Set global constants ps = 1000. #rh = 1.e-30#Relative humidity (not used) #rhbdd = 1.e-30 (bit ysed #dt = 24.*3600. #time step in seconds (not used) #---Set up pressure array (a global)---- ptop = 1. #Top pressure in mb pstart = .995*ps rat = (ptop/pstart)**(1./r.nlev) logLevels = [pstart*rat**i for i in range(r.nlev)] logLevels.reverse() levels = [ptop + i*(pstart-ptop)/(r.nlev-1) for i in range(r.nlev)] #p=numpy.array(levels) p = numpy.array(logLevels) #---Temperature and moisture arrays (initialized) T = numpy.zeros(r.nlev) + 230. q = numpy.zeros(r.nlev) + 1.e-30 #Assume dry atmosphere #Set the co2 value. In the current version of climt_lite # the values are passed in the call to r as keyword arguments. co2 = 300. # #Set the (constant) insolation. #Note this is now passed as a parameter to the radiation object solin = 410. #NOTE: Calculation in the book stated 400., but with new climt #410 needed to give same top-of-atmosphere flux. Possible change #in default surface albedo? #Set the ozone profile o3max = 2.*(3.*16./29.)*2.e-6 #Set o3max to zero if you want to zero out ozone. #NOTE: Book says 2e-6 molar, but climt_lite now says mass mixing ratio is used, not molar # Using corresponding mass mixing ratio to 2.e-6 yields too little absorption, # but the above (4e-6 molar) is needed to yield the same absorption as in the figure # in the book. o3= numpy.zeros(r.nlev) for i in range(n): o3[i] = max(o3max*math.exp(-(p[i]-20.)**2/50.**2),1.e-20) #Define function to do time integration for n steps def steps(Tg,T,nSteps,dtime): for i in range(nSteps): Tg = T[-1] #let it float #Do smoothing if i%10 == 0: for j in range(1,len(T)-1): T[j] = .25*T[j-1] + .5*T[j] + .25*T[j+1] # r(p=p,ps=ps,T=T,Ts=Tg,q=q,co2=co2,o3=o3,solin=solin, aldif=0.,aldir=0.,asdif=0.,asdir=0.)#Set all surface albedos to zero fluxLW,heatLW = -r.lwflx[:,0,0],r.lwhr[:,0,0] fluxStellar,heatStellar = -r.swflx[:,0,0],r.swhr[:,0,0] #Convert heating rates to K/s (for climt_lite) heatLW /= 24.*3600. heatStellar /= 24.*3600. # flux = fluxLW+fluxStellar if doStellarAbs: heat = heatLW+heatStellar else: heat = heatLW dT = heat*dtime #Limit the temperature change per step dT = numpy.where(dT>5.,5.,dT) dT = numpy.where(dT<-5.,-5.,dT) #Midpoint method time stepping r(p=p,ps=ps,T=T+.5*dT,Ts=Tg+.5*dT[-1],q=q,co2=co2,o3=o3,solin=solin, aldif=0.,aldir=0.,asdif=0.,asdir=0.)#Set all surface albedos to zero) fluxLW,heatLW = -r.lwflx[:,0,0],r.lwhr[:,0,0] fluxStellar,heatStellar = -r.swflx[:,0,0],r.swhr[:,0,0] #Convert heating rates to K/s (for climt_lite) heatLW /= 24.*3600. heatStellar /= 24.*3600. # flux = fluxLW+fluxStellar if doStellarAbs: heat = heatLW+heatStellar else: heat = heatLW dT = heat*dtime #Limit the temperature change per step dT = numpy.where(dT>5.,5.,dT) dT = numpy.where(dT<-5.,-5.,dT) T += dT # dTmax = max(abs(dT)) #To keep track of convergence #Using flux[0] in the following is not a mistake. #This is just a trick to relax towards satisfying #the TOA balance. The formulation assumes that #turbulent fluxes keep the low level air temperature #near the ground temperature. #Estimate the surface temperature change needed to #bring the TOA energy budget into balance Trad = (fluxLW[0]/phys.sigma)**.25 dTg = -flux[0]/(4.*phys.sigma*Trad**3) #Relax partway towards that value Tg = Tg + .1*dTg # Tad = Tg*(p/p[-1])**Rcp T = numpy.where(T to the argument list # of the radiation object r, where is the value you # want, or a variable containing the value. Rcp = 2./7. # #Initialize the temperature #c = readTable(outputDir+'TTestTrop80.txt') #p = c['p'] #T = c['T'] #Tg = T[-1] Tg = 295. T = Tg*(p/p[-1])**Rcp # #--------------------------------------------------------- #--------------Initializations Done----------------------- #--------------Now do the time stepping------------------- #--------------------------------------------------------- # #Note that since vertical smoothing is called every 10 timesteps #in step(...), if you let the timestep get to small in the following #loop without changing the frequency of smoothing, you start to #get spurious vertical temperature mixing due to smoothing. #To check convergence of the stratosphere, make sure the net #heating in your converged solution is near zero in the stratosphere. for i in range(0,350): nout = 20*i print dtime if (i%50 == 0) & (i > 200): dtime = .5*dtime Tg,Tad,T,flux,fluxStellar,fluxLW,heat,heatStellar,heatLW = steps(Tg,T,20,dtime) print 'History step',Tg,flux[-1],max(heat),min(heat) history(nout,caseTag) #Compute a flux profile for a model temperature profile #The UV/viz absorption are temperature independent, so we #can just pick any old profile Tg = 300. T[:] = 300. Tad = Tg*(p/p[-1])**Rcp r(p=p,ps=ps,T=T,Ts=Tg,co2=1.e-9,q=q,o3=o3,solin=solin, aldif=0.,aldir=0.,asdif=0.,asdir=0.)#Set all surface albedos to zero, Zero out co2 fluxLW,heatLW = -r.lwflx[:,0,0],r.lwhr[:,0,0] fluxStellar,heatStellar = -r.swflx[:,0,0],r.swhr[:,0,0] #Convert heating rates to K/s (for climt_lite) heatLW /= 24.*3600. heatStellar /= 24.*3600. # flux = fluxLW+fluxStellar if doStellarAbs: heat = heatLW+heatStellar else: heat = heatLW history(0,'FluxPlot300K')