import climt_lite as climt
from ClimateUtilities import *
import phys
#Set default values of globals
g = 9.8 #Gravity
#ToDo: See warnings in the comments on the OLR function
#
# General questions: Does NCAR radiation us mass mixing ratio
# of water, or specific humidity? This affects high pressure
# and very warm calculations.
# Does ClimT use the same water specification as NCAR, or does
# it do a translation?
#
#
#ToDo: Generalize the specification of the profile, so we
# can easily switch to a dry adiabat or a nonideal gas
# adiabat. It can be handled by using a plug-compatible
# adiabat function.
#
#
# The argument list for the radiation object r(...) has
# gotten unwieldy. Make up an object to simplify this
#
# The setting of the surface albedos to zero has not yet been
# implemented, nor is there a provision to set the albedos.
# The default values are .07. We should set default
# to zero, and allow user to set albedos
#
# The ccmrad package uses mb for pressure and ppmv for greenhouse
# gas concentrations. This ought to be harmonized with the
# miniClimt radiation package, which uses Pascals for pressure
# and mass specific concentrations.
#
# Make a radcomp function analogous to miniClimt, so that
# this package can be used to do radiative-convective models
# in a way analogous to miniClimt, without calling climt directly.
# radcompCCM can return, [flux,heatIR,heatSolar], and can
# take a few extra composition arguments (co2, qH2O,ch4,o3 ...).
# Currently this isn't used in Chapter 4, but it could be useful for
# later chapters.
#
# The call to setTstrat, which is supposed to
# set the globals, doesn't work reliably. Check this out.
# What is more reliable is to "import ccmradFunctions as ccm"
# then do things like ccm.Tstrat = 200.
# This may also be a problem with setGravity
#
# Changing the liquid particle size doesn't seem to affect the OLR
# Is this a bug, or something that ccmrad assumes? The OLR
# does respond correctly to the ice particle size, though,
# which is what was used in the previous version of climt
# to make the plots of particle size effect in Chapter 5.
#
#
#
#Remarks: Note that the ground temperature Ts specified in ClimT and
#NCAR doesn't just affect the radiation through sigma Ts**4. Evidently
#it does an estimate of the boundary layer (subgrid) opacity and
#uses this to compute the net upward flux. The resulting flux gives
#the net cooling of the surface itself correctly. It's probably
#necessary to do this because even just the water in 2 meters
#of air has substantial opacity.
#Changes: Surface pressure now incorporates pressure of CO2.
# This makes the NCAR radiation calculation of OLR accurate
# up to around 2 bars of CO2. The surface radiation at
# CO2 over .2 bar seems strange in the presence of water vapor,
# though, and so the low level fluxes may not be as accurate.
# This needs to be checked against the KPA radiation code.
# (Total pressure computation in GetMoistProfile updated
# 6/29/2009 to be accurate at higher mole fractions of CO2.
# Note that pAir/g is NOT the mass path of air when
# CO2 concentration is high!)
#
# I have not yet incorporated the effect of water vapor
# or CH4 on surface pressure. This is not so important
# for CH4, which rarely gets that massive, but it is
# relevant for water vapor in hot atmospheres.
#
# I have zeroed out the ozone. We could allow the user
# to specify an ozone profile, if desired
#
# I have replaced the moist adiabat computation by
# the climate book version.
#
# 5/2010 I have modified this to work with the new
# version of climt_lite. Not yet fully checked
#General notes on the radiation computation done in
#this script:
#
#Note that unless the vertical resolution is increased
#a great deal, this calculation may not get the OLR right
#for temperatures above 320K, when there is a lot of water
#(hence opacity) in the upper atmosphere.
#-----------------Function definitions-----------------------------
#
#
#These two routines set globals used in the
#other computations: the non-CO2 surface pressure,
#the acceleration of gravity, and the stratospheric
#temperature.
#They are used so that these infrequently changed
#parameters don't have to be included as function arguments.
def setGravity(gval):
global g
g = gval# Gravity in m/s**2
def setTstrat(StratosphereTemperature):
global Tstrat
Tstrat = StratosphereTemperature
#**ToDo: Modify the OLR function to allow the optional
#input of a cloud profile parameter. Do that as a keyword
#argument. (Note that AllRad essentially does all of that,
#and returns shortwave albedo information as well).
#Utility function to set the modified moist adiabat
#temperature and moisture profiles, and set up pressure array
#Note that the co2 concentration (ppmv) is interpreted as
#the mole fraction in the total atmosphere, when the concentration
#is high. (Not the most convenient way to specify things at
#high concentrations. Would be better to specify mass path
#of CO2 and mass path of air.
def getMoistProfile(psAir,Ts,rh,co2=300.):
ps = psAir/(1.-1.e-6*co2) #Corrected computation of ps)
ps = ps/100. #Convert to mb.
p = setpLin(ps)
# p = setpLog(ps,1.,r.nlev) #Top at 1 mb
m = phys.MoistAdiabat()
p1,T,molecon,q = m(100.*ps,Ts,100.*p)#Convert pressure to Pa
q = 1000.*q*rh #convert to g/kg, and set relative humidity
#----
#Simple isothermal Stratosphere model:
T = numpy.where(TTstrat,q,1.e30)
qmin = min(q)
q = numpy.where(q>1000.,qmin,q)
return ps,p,T,q
#Define OLR function. Greenhouse gas amounts
# (e.g. CO2 and CH4) are molar concentrations expressed
#in ppm. (essentially the mixing ratio, for small
#concentratins). Note that the effect of water vapor on surface
#pressure is not currently included, so this won't
#work properly at temperatures much above 310K. The effect
#of added CO2 on surface pressure is included, but not the
#effect of methane; also though the pressure broadening
#from increased surface pressure is taken into account, the
#ccm radiation model does not incorporate the effect of
#the difference between self-broadening and foreign broadening
#of CO2.
#
#psAir is specified in Pa
#
#The cloud water profile can now be passed to
#the function as a keyword argument, cloudWater. .
#If the argument is present, the cloud fraction is set to 1
#at all levels
#
def OLR(psAir,Ts,rh = .5,co2=300.,ch4=1.e-30,**kwargs):
#Check keyword arguments for cloud water
if 'cloudWater' in kwargs.keys():
clwp = kwargs['cloudWater']
#Set up cloud fraction profile
cldf = numpy.zeros(r.nlev,numpy.Float) + 1.
else:
clwp = None
#
#Allow specification of other trace gases as keyword args
#(so far just N2O)
n2o =small
if 'N2O' in kwargs.keys():
n2o = kwargs['N2O']
#
#Set up the profile
ps,p,T,q = getMoistProfile(psAir,Ts,rh,co2)
#
#Note that in the old version of climt, things
#like r.params.value['co2'] = 300. would set the
#value. Now this no longer works; the r.Params is
#now used only for output. The input values must
#be set as arguments to r(...)
#---------
#Do the computation:
if clwp == None:
r(p=p,ps=ps,T=T,Ts=Ts,q=q,co2 = co2,n2o = n2o, ch4=ch4,g=g,o3 = small) #Last arg zeroes out ozone
else:
r(p=p,ps=ps,T=T,Ts=Ts,q=q,co2=co2,ch4=ch4,n2o=n2o,clwp=clwp,cldf=cldf,g=g,o3 = small) #Cloud case
return -r.lwflx[0][0][0] #Extra indices make sure we return a number not array
#This is just like the OLR function, except it returns
#the shortwave albedo as well as the OLR
#**ToDo: Make this optionally return longwave and shortwave
#heating rates, for use in radiative-convective models.
#Also profiles of atmospheric solar absorption? Also allow
#for changing ozone?
#
#ToDo: Put particle size in the call to the radiation
#method. It's not there now, because the old version of climT
#left r_liq out of the kwarg list.
#
def AllRad(psAir,Ts,rh = .5,co2=300.,ch4=1.e-30,**kwargs):
#Check if we want to return profiles of flux and heating
#If doProfile = True, we will return flux and heating profiles
#instead of just the albedo and OLR. This is useful for doing
#radiative-convective models.
#
#(Note: to be useful for radiative-convective models,
#we have to have an option to input T and q as well. )
doProfile = False
if 'doProfile' in kwargs.keys():
doProfile = kwargs['doProfile']
#
#Check for specification of ozone profile
o3 = small
if 'o3' in kwargs.keys():
o3 = kwargs['o3']
#Allow specification of other trace gases as keyword args
#(so far just N2O)
n2o =small
if 'N2O' in kwargs.keys():
n2o = kwargs['N2O']
#Check keyword arguments for cloud water
if 'cloudWater' in kwargs.keys():
clwp = kwargs['cloudWater']
#Set up cloud fraction profile
cldf = numpy.zeros(r.nlev,numpy.float) + 1.
else:
clwp = None
#
#Check for particle size
if 'r_liq' in kwargs.keys():
r_liq = kwargs['r_liq']
else:
r_liq = 10.*numpy.ones(r.nlev,numpy.float)
if 'r_ice' in kwargs.keys():
r_ice = kwargs['r_ice']
else:
r_ice = 30.*numpy.ones(r.nlev,numpy.float)
#Set up the profile
ps,p,T,q = getMoistProfile(psAir,Ts,rh,co2)
#Do the computation:
if clwp == None:
r(p=p,ps=ps,T=T,Ts=Ts,q=q,co2 = co2,n2o = n2o, ch4=ch4\
,solin = S,zen=0.,g=g,o3 = o3) #Last arg zeroes out ozone
else:
r(p=p,ps=ps,T=T,Ts=Ts,q=q,clwp=clwp,cldf=cldf,r_ice = r_ice,co2 = co2,\
n2o = n2o, ch4=ch4,solin = S,zen=0.,g=g,o3 = o3) #Cloud case
if doProfile:
#**Need to append some indices?
return p,r.lwflx,r.swflx,r.lwhr,r.swhr
else:
#Extra indices make sure we return a number not array
return -r.lwflx[0][0][0],1.-r.swflx[0][0][0]/S,\
(r.swflx[0][0][0]-r.swflx[-1][0][0])/S
#"Canonical" dry air OLR function
#psAir is specified in Pascal
#
#Note that the ccm radiation code describes the input greenhouse
#gas concentrations as "volumetric mixing ratios," but the actual
#computes the path by multiplying total pressure by the mixing ratio
#Hence, if we take into account the co2 when computing the pressure,
#we are really treating the value as a molar concentration, not as
#a mixing ratio. The maximum would thus be 1000000 ppmv. The ccm
#radiation code is not really supposed to be valid at concentrations
#where the distinction between mixing ratio and concentration matters,
#but treating things this way allows us to push the envelope of validity
#and at least crudely incorporate co2 self-broadening.
#
#The surface pressure adjustment makes little difference out
#to 10% CO2, but as one goes to 20% and above it starts to become
#important
def OLRDryAir(psAir,Ts,co2 = 300.,ch4 = 1.e-30,**kwargs):
#Allow specification of other trace gases as keyword args
#(so far just N2O)
n2o =small
if 'N2O' in kwargs.keys():
n2o = kwargs['N2O']
#
ps = (psAir/100.)/(1.-1.e-6*co2) #Total surface pressure includes CO2
p = setpLin(ps)
T = Ts*(p/ps)**(2./7.) #Ignores effect of co2 on r/cp
T = numpy.where(TTstrat,q,1.e30)
qmin = min(q)
q = numpy.where(q>1000.,qmin,q)
#----
#Reset the boundarylayer humidity
for i in range(len(q)):
if p[i]>p[-1]-50.:
q[i] = q[i]*rhbdd/rh
r(p=p,ps=ps,T=T,Ts=Tg,q=q,co2 = co2 , ch4=ch4,g=g,o3 = small)
#Extra indices make sure we return a number not array
return (-r.lwflx[0][0][0],-r.lwflx[-1][0][0])
def SurfCoeffs(psAir,Ts,rh = .5,rhbdd=.8,co2=300.,ch4=1.e-30):
olr,IRcool1 = BddRad(psAir,Ts,Ts,rh,rhbdd,co2,ch4)
olr,IRcool2 = BddRad(psAir,Ts,Ts+1.,rh,rhbdd,co2,ch4)
g = IRcool1/(phys.sigma*Ts**4)
gs = (IRcool2-IRcool1)/(4.*phys.sigma*Ts**3)
return g,gs
#------------------------------------------------------------------------
#This function computes the coefficient determining the incremental
#atmospheric transmission caused by nonzero air/ground temperature
#temperature difference Tg-Ts. It is used in energy balance models
#for atmospheres which are not optically thick, in which case the
#OLR depends (somewhat) on Tg-Ts. The effect of Tg-Ts on
#OLR is TopCoeff*(Tg-Ts). Note that this is NOT the same as the
#coefficient (1-a_+) in Eq. (6.6) in the text, which refers to
#the transmission of the net upwelling from the surface not the
#incremental upwelling. TopCoeff is the coefficient giving the linearization
#of OLR as a function of Tg, about the state Tg = Tsa. It is mostly
#the coefficient one wants when doing computations solutions of a
#radiative-convective model for an atmosphere that is not optically
#thick (in which case the OLR depends on Tg-Tsa).
#
#Note: There was a bug in BddRad which didn't properly use the
#ground temperature, and caused TopCoeff to return zero. This
#was fixed 11/7/2011
def TopCoeff(psAir,Ts,rh = .5,rhbdd=.8,co2=300.,ch4=1.e-30):
olr1,IRcool = BddRad(psAir,Ts,Ts,rh,rhbdd,co2,ch4)
olr2,IRcool = BddRad(psAir,Ts,Ts+1.,rh,rhbdd,co2,ch4)
gtop = (olr2-olr1)/(4.*phys.sigma*Ts**3)
return gtop
#Function to get linear fit a and b coefficient,
#for OLR = a + b*(T-T0)
def OLRcoeffs(psAir,T0,T1,rh = .5,co2=300.,ch4=1.e-30):
a = OLR(psAir,T0,rh,co2,ch4)
b = (OLR(psAir,T1,rh,co2,ch4)-OLR(psAir,T0,rh,co2,ch4))/(T1-T0)
return a,b
#CLIMT moist adiabat function no longer needed, since
#we now have our own. But remember to convert q to g/kg!
#Utility function to set up a logarithmic pressure grid
def setpLog(ps,ptop,n):
pfact = (ps/ptop)**(1./(n-1))
p = [ptop]
for i in range(1,n):
p.append( pfact*p[i-1])
return numpy.array(p)
#Utility function to set up a linear pressure grid
# (Fix this so that the lowest level is pbot? But
# NCAR wants pressures on half-levels, staggered from surface
def setpLin(pbot):
return (numpy.arange(r.nlev, dtype ='d')+ 0.5 ) * pbot/r.nlev
#To Do: Incorporate a simple nth order polynomial fit.
#Alternately, could use routine polint from Numerical Recipes
#to generate OLR from a table.
#Initialize globals
Tstrat = 150. #Stratospheric temperature
S = 300. #Incoming solar, for shortwave calculation
#This is now passed as an argument to the call
#(used in AllRad only, to compute albedos)
#---Instantiate the radiation model---
#REMINDER: climt wants moisture mass mixing ratio in g/kg !
# Please keep that in mind if you modify the script to
# set the q profile yourself.
r = climt.radiation()
#
#Set albedos (for shortwave calculation
#**ToDo: These no longer work. Need to set them as arguments
#The default value of the albedos is .07
#r.Params.value['aldir']= 0.
#r.Params.value['aldif']= 0.
#r.Params.value['asdir'] = 0.
#r.Params.value['asdif'] = 0.
#Define array of near-zeros for zeroing out ozone, water vapor, etc
small = numpy.zeros(r.nlev,numpy.float)+1.e-30