#------------------------------------------------------------------------------
#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')