#This script computes pure radiative equilibrium by timestepping.
#It can use a variety of different radiation codes, and different
#models of the transmission function.
#This version is different from RadConvEq because it
# (a) incorporates stellar heating and (b) Time-steps
# the surface temperature so that one finds the temperature
# which is in equilibrium with a given solar forcing. With
# solar absorption, one can't do the problem by just holding
# Tg fixed and finding the corresponding atmosphere that's in
# equilibrium with that
#This is an experimental version, originally modified for the AGU2008 talk.
#It has been modified to compute absorption of stellar flux
#impinging at TOA. (mainly for the M-dwarf problem, but useful
#for the discussion of solar heating in the scattering chapter)
#Eventually, this should probably change the zenith angle
#representation, and also incorporate scattering.
#As it stands, the code is mostly useful for absorption of nir.
#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.
#
#Note: The stellar flux past the shortwave cutoff had to be
#added in to the calculation of incoming flux. This needs
#to be put in the main code-tree.
#ToDo: After the data files are moved to the Workbook Dataset
#directory, update the band-data load section so it can find them
#ToDo:
# *Make it possible to use ccm radiation as an option, so a separate
# script isn't needed for ccm computations.
# ->Need radcomp in ccmradFunctions. Return heating in W/kg
# *make it possible to use the fancy version of miniclimt
# *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 = '5'
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 15000 cm**-1. Use
#for calculations with near-IR solar heating
radmodel.bandData = radmodel.loadExpSumTable(EsumPath+'CO2TableExtendedWave.260K.100mb.self.data')
radmodel.BackgroundGas = phys.N2 #The transparent background gas (**Fix self-absorption issue)
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 (Don't use this with solar absorption)
#radmodel.Trans = radmodel.TransOobleck (Don't use this with solar absorption
#Set the continuum, if there is one
radmodel.Continuum = radmodel.CO2Continuum
#radmodel.Continuum = radmodel.NoContinuum
#Define function to do time integration for n steps
def steps(Tg,T,nSteps,dtime):
for i in range(nSteps):
#Tg = T[-1] (No longer needed. Tg computed by iteration)
#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]
#
fluxLW,heatLW = radmodel.radCompLW(p,T,Tg,q)
fluxStellar,heatStellar = radmodel.radCompStellar(p,T,Tg,q)
fluxStellar += -OutBandStellarFlux #Stellar flux past shortwave cutoff
flux = fluxLW+fluxStellar
if doStellarAbs:
heat = heatLW+heatStellar
else:
heat = heatLW
dT = heat*dtime
#Limit the temperature change per step
dT = Numeric.where(dT>5.,5.,dT)
dT = Numeric.where(dT<-5.,-5.,dT)
#Midpoint method time stepping
fluxLW,heatLW = radmodel.radCompLW(p,T+.5*dT,Tg+.5*dT[-1],q)
fluxStellar,heatStellar = radmodel.radCompStellar(p,T+.5*dT,Tg+.5*dT[-1],q)
fluxStellar += -OutBandStellarFlux #Stellar flux past shortwave cutoff
flux = fluxLW+fluxStellar
if doStellarAbs:
heat = heatLW+heatStellar
else:
heat = heatLW
dT = heat*dtime
#Limit the temperature change per step
dT = Numeric.where(dT>5.,5.,dT)
dT = Numeric.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 = Numeric.where(T