#------------------------------------------------------------------------- # #Description: # #This script computes pure radiative equilibrium by timestepping. #It can use a variety of different radiation codes, and different #models of the transmission function. It can also do radiative-convective #equilibrium if the hard convective adjustment step is uncommented. # #This version just does equilibrium with thermal IR radiation. It #does not incorporate solar absorption in the atmosphere. In this #version, the ground temperature Tg is specified, and the atmosphere #in equilibrium with the specified temperature is found. The OLR #is thus computed as a function of ground temperature. #-------------------------------------------------------------------------- #ToDo: # *Make it possible to use ccm radiation. # ->Need radcompLW in ccmradFunctions. Return heating in W/kg # *make it possible to use the fancy version of miniclimt # *Make a script that computes the pure radiative equilibria # in Figure {fig:RealGasRadEq}, and also the logarithmic slopes. # *Handle the contribution of water vapor to surface pressure # in setting up the pressure grid. (Tricky, because # temperature is changing!) That's important when # surface temperatures are much over 300K. The changing # surface pressure makes the time stepping rather tricky. #Data on section of text which this script is associated with Chapter = '4' Section = '**' Figure = '**' # from ClimateUtilities import * import phys import miniClimt as radmodel import planets #Path to the workbook datasets datapath = '/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/' #Path to the exponential sum tables EsumPath = datapath + 'Chapter4Data/ExpSumTables/' #--------------------Read in band data,set radmodel options------------------- # #Data for principal band of CO2. Good up to about 10% CO2 in 1bar air. radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TablePrinBand.260K.100mb.self.data') # #Data for the full CO2 spectrum out to 2300 cm**-1. Use for Venus #and Early Mars calculations #radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableAllWave.260K.100mb.self.data') #radmodel.bandData = radmodel.bandData[:-9] #Eliminate "transparent region" from table # #Data for two-band Oobleck. To do one band, just make #a list having only semiGreyBandParams alone. #radmodel.bandData = \ # [radmodel.semiGreyBandParamsM,radmodel.semiGreyBandParams,radmodel.semiGreyBandParamsP] # radmodel.BackgroundGas = phys.air #The transparent background gas radmodel.GHG =phys.CO2 #The greenhouse gas # radmodel.pref = 1.e4 #Reference pressure at which band data is given radmodel.Trans = radmodel.TransEsums #radmodel.Trans = TransMalkmus #radmodel.Trans = radmodel.TransOobleck #Set the continuum, if there is one #radmodel.Continuum = radmodel.CO2Continuum radmodel.Continuum = radmodel.NoContinuum #Set the gravity radmodel.g = planets.Earth.g #---------------------------------------------------------------------------- #-------------------Initialize arrays #-------Set up the pressure array------------------------------ # We typically use 60 points with extra resolution near the # ground when doing pure radiative equilibrium, where there # is a radiative boundary layer near the ground that needs to # be resolved. For radiative-convective equilibrium, the adiabat # takes over near the ground, and setpLin(...), which gives uniform # resolution, does pretty well with as little as 20 points. If the # tropopause gets too high up, though, you may need extra resolution # in the stratosphere, and in that case using a logarithmic grid # can be advantageous. # # Caution: If you go to surface temperatures much above 300K, # then you need to keep track of the total surface pressure # vs. the surface pressure of dry air. That's not done here. #-------------------------------------------------------------- n = 60 #n = 60 ps = 1.e5 ptop = 100. #Use 100 except for Venus, where we use 1.e4 #p = radmodel.setpExtraGroundRes(ps,ptop,n) p = radmodel.setpLin(ps,ptop,n) #p = radmodel.setpLog(ps,ptop,n) #Puts extra resolution in stratosphere. #---------Initialize the temperature profile, set up the adiabat--------- Tg = 280. # Set the temperature of the ground #Use the following for the dry adiabat Tad = Tg*(p/ps)**phys.air.Rcp #Change gas if desired #Use the following for the moist adiabat #m = phys.MoistAdiabat(phys.water,phys.air) #p1,Tad,molarCon,massCon = m(ps,Tg,p) #First arg. is just the dry air ps! # T = Tg*numpy.ones(len(p),numpy.Float) #Initialize constant T #T[:] = Tad[:] #Alternate initialization. This converges faster #if you are doing a radiative-convective case, with #convective adjustment. #------------------------------------------------------------------------ #------------Set the greenhouse gas concentration------------------------ #q = 10.e-4 #Makes optical depth 10 for Oobleck q = (44./29.)*300.e-6 #Earthlike CO2 #q = 1. #Use for Mars and Venus #------------------------------------------------------------------------ #----------------Set initial time step-------------------------------- dtime = 50. #Timestep in days; 5 days is the usual for Earthlike case #For the radiative convective case, you can get away with #using 50 for the first few hundred time steps, then #reducing to 5 or less later as the solution converges # heatfact = 24.*3600./1000. #Rough conversion to K/day #If we are only interested in the steady #state, and not in accurately reproducing #the detailed time dependence, it doesn't #matter if we use the exact specific heat #here to convert the heating rate in W/kg #into Kelvins/day dtime = dtime*heatfact #---------------------------------------------------------------------- #Define function to do time integration for n steps def steps(T,nSteps,dtime): for i in range(nSteps): #Do smoothing if i%100 == 0 & i>0: for j in range(1,len(T)-1): T[j] = .25*T[j-1] + .5*T[j] + .25*T[j+1] # flux,heat = radmodel.radCompLW(p,T,Tg,q) 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 flux,heat = radmodel.radCompLW(p,T+.5*dT,Tg,q) 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 #Uncomment next line to do hard convective adjustment T = numpy.where(T