from builtins import range
import numpy as np
from climlab.radiation.greygas import GreyGas
from climlab import constants as const
from climlab.domain import domain, axis, Field
from copy import copy
[docs]
class NbandRadiation(GreyGas):
r'''Process for radiative transfer.
Solves the discretized Schwarschild two-stream equations
with the spectrum divided into N spectral bands.
Every NbandRadiation object has an attribute
``self.band_fraction``
with sum(self.band_fraction) == 1
that gives the fraction of the total beam in each band
Also a dictionary
``self.absorber_vmr``
that gives the volumetric mixing ratio of every absorbing gas
on the same grid as temperature
and a dictionary
``self.absorption_cross_section``
that gives the absorption cross-section per unit mass for each gas
in every spectral band
'''
def __init__(self, absorber_vmr=None, **kwargs):
super(NbandRadiation, self).__init__(**kwargs)
newinput = ['band_fraction',
'absorber_vmr',
'absorption_cross_section',
'cosZen',]
self.declare_input(newinput)
# this should be overridden by daughter classes
self.band_fraction = np.array(1.)
## a dictionary of absorbing gases, in volumetric mixing ratios
# each item should have dimensions of self.Tatm
# Can be passed as input argument
if absorber_vmr is None:
absorber_vmr = {}
self.absorber_vmr = absorber_vmr
# a dictionary of absorption cross-sections in m**2 / kg
# each item should have dimension... (num_channels, 1)
self.absorption_cross_section = {}
self.cosZen = 1. # cosine of the average zenith angle
dp = self.Tatm.domain.lev.delta
self.mass_per_layer = dp * const.mb_to_Pa / const.g
self.albedo_sfc = np.ones_like(self.band_fraction) * self.albedo_sfc
self._define_diagnostics()
def _define_diagnostics(self):
atmaxes = copy(self.Tatm.domain.axes)
atmaxes.update(self.channel_ax)
atmdom = domain.Atmosphere(axes=atmaxes)
sfcaxes = copy(self.Ts.domain.axes)
sfcaxes.update(self.channel_ax)
sfcdom = domain.SlabOcean(axes=sfcaxes)
atmdiag = Field(np.zeros(atmdom.shape), domain=atmdom)
sfcdiag = Field(np.zeros(sfcdom.shape), domain=sfcdom)
self.add_diagnostic('flux_from_sfc', 0. * sfcdiag)
self.add_diagnostic('flux_to_sfc', 0. * sfcdiag)
self.add_diagnostic('flux_to_space', 0. * sfcdiag)
self.add_diagnostic('absorbed', 0. * atmdiag)
self.add_diagnostic('absorbed_total', 0. * sfcdiag)
self.add_diagnostic('emission', 0. * atmdiag)
self.add_diagnostic('emission_sfc', 0. * sfcdiag)
@property
def band_fraction(self):
return self._band_fraction
@band_fraction.setter
def band_fraction(self, value):
self.num_channels = value.size
# abstract axis for channels
ax = axis.Axis(num_points=self.num_channels)
self.channel_ax = {'channel': ax}
dom = domain._Domain(axes=self.channel_ax)
# fraction of the total solar flux in each band:
self._band_fraction = Field(value, domain=dom)
def _compute_optical_path(self):
# this will cause a problem for a model without CO2
tau = np.zeros_like(self.absorber_vmr['CO2']*
self.absorption_cross_section['CO2'])
for gas, vmr in self.absorber_vmr.items():
# convert to mass of absorber per unit total mass
if gas == 'H2O': # H2O is stored as specific humidity, not VMR
q = vmr
else:
q = vmr / (1.+vmr)
try:
# if this gas isn't present in absorption dictionary
# the assumption is that there is no absorption!
kappa = self.absorption_cross_section[gas]
tau += q * kappa
except: pass
tau *= self.mass_per_layer / self.cosZen
return tau
def _compute_absorptivity(self):
# assume that the water vapor etc is current
optical_path = self._compute_optical_path()
# account for finite layer depth
absorptivity = 1. - np.exp(-optical_path)
axes = copy(self.Tatm.domain.axes)
# add these to the dictionary of axes
axes.update(self.channel_ax)
dom = domain.Atmosphere(axes=axes)
self.absorptivity = Field(absorptivity, domain=dom)
def _compute_emission_sfc(self):
# need to split the total emission across the bands
total_emission = super(NbandRadiation, self)._compute_emission_sfc()
return self._split_channels(total_emission)
def _compute_emission(self):
# need to split the total emission across the bands
total_emission = super(NbandRadiation, self)._compute_emission()
band_fraction = self.band_fraction
for n in range(self.Tatm.domain.numdims):
band_fraction = band_fraction[:, np.newaxis]
return total_emission * band_fraction
def _compute_radiative_heating(self):
# need to recompute transmissivities each time because
# water vapor is changing
self._compute_absorptivity()
super(NbandRadiation, self)._compute_radiative_heating()
def _split_channels(self, flux):
split = np.outer(self.band_fraction, flux)
# make sure there's a singleton dimension at the last axis (level)
if np.size(split, axis=-1) != 1:
split = split[..., np.newaxis]
return split
def _join_channels(self, flux):
return np.sum(flux, axis=0)
[docs]
class ThreeBandSW(NbandRadiation):
def __init__(self, emissivity_sfc=0., **kwargs):
r'''A three-band mdoel for shortwave radiation.
The spectral decomposition used here is largely based on the
"Moist Radiative-Convective Model" by Aarnout van Delden, Utrecht University
a.j.vandelden@uu.nl
http://www.staff.science.uu.nl/~delde102/RCM.htm
Three SW channels:
channel 0 is Hartley and Huggins band (UV, 1%, 200 - 340 nm)
channel 1 is Chappuis band (27%, 450 - 800 nm)
channel 2 is remaining radiation (72%)
'''
super(ThreeBandSW, self).__init__(emissivity_sfc=emissivity_sfc, **kwargs)
# fraction of the total solar flux in each band:
self.band_fraction = np.array([0.01, 0.27, 0.72])
if 'CO2' not in self.absorber_vmr:
self.absorber_vmr['CO2'] = 380.E-6 * np.ones_like(self.Tatm)
if 'O3' not in self.absorber_vmr:
self.absorber_vmr['O3'] = np.zeros_like(self.Tatm)
if 'H2O' not in self.absorber_vmr:
self.absorber_vmr['H2O'] = self.q
## absorption cross-sections in m**2 / kg
O3 = np.array([200.E-24, 0.285E-24, 0.]) * const.Rd / const.kBoltzmann
#self.absorption_cross_section['O3'] = np.reshape(O3,
# (self.num_channels, 1))
#H2O = np.array([0.002, 0.002, 0.002])
H2O = np.array([0., 0., 0.001])
for n in range(self.Tatm.domain.numdims):
H2O = H2O[:, np.newaxis]
O3 = O3[:, np.newaxis]
self.absorption_cross_section['O3'] = O3
self.absorption_cross_section['H2O'] = H2O
#self.absorption_cross_section['H2O'] = np.reshape(H2O,
# (self.num_channels, 1))
self.absorption_cross_section['CO2'] = \
np.zeros_like(self.absorption_cross_section['O3'])
self.cosZen = 0.5 # cosine of the average solar zenith angle
self._define_diagnostics()
@property
def emissivity(self):
# This ensures that emissivity is always zero for shortwave classes
return np.zeros_like(self.absorptivity)
[docs]
class FourBandSW(NbandRadiation):
# The most recent RMCM program uses a four-channel SW model
# this is probably a better way to do it... distinguishes between
# visible band with no absorption and near-infrared with weak H2O absorption
# But this needs some tuning and better documentation
def __init__(self, emissivity_sfc=0., **kwargs):
r'''A four-band mdoel for shortwave radiation.
The spectral decomposition used here is largely based on the
"Moist Radiative-Convective Model" by Aarnout van Delden, Utrecht University
a.j.vandelden@uu.nl
http://www.staff.science.uu.nl/~delde102/RCM.htm
Four SW channels:
channel 0 is Hartley and Huggins band (UV, 6%, <340 nm)
channel 1 is part of visible with no O3 absorption (14%, 340 - 500 nm)
channel 2 is Chappuis band (27%, 500 - 700 nm)
channel 3 is near-infrared (53%, > 700 nm)
'''
super(FourBandSW, self).__init__(emissivity_sfc=emissivity_sfc, **kwargs)
# fraction of the total solar flux in each band:
self.band_fraction = np.array([0.06, 0.14, 0.27, 0.53])
if 'CO2' not in self.absorber_vmr:
self.absorber_vmr['CO2'] = 380.E-6 * np.ones_like(self.Tatm)
if 'O3' not in self.absorber_vmr:
self.absorber_vmr['O3'] = np.zeros_like(self.Tatm)
if 'H2O' not in self.absorber_vmr:
self.absorber_vmr['H2O'] = self.q
## absorption cross-sections in m**2 / kg
O3 = np.array([200.E-24, 0., 0.285E-24, 0.]) * const.Rd / const.kBoltzmann
#self.absorption_cross_section['O3'] = np.reshape(O3,
# (self.num_channels, 1))
H2O = np.array([0., 0., 0., 0.0012])
for n in range(self.Tatm.domain.numdims):
H2O = H2O[:, np.newaxis]
O3 = O3[:, np.newaxis]
self.absorption_cross_section['O3'] = O3
self.absorption_cross_section['H2O'] = H2O
#self.absorption_cross_section['H2O'] = np.reshape(H2O,
# (self.num_channels, 1))
self.absorption_cross_section['CO2'] = \
np.zeros_like(self.absorption_cross_section['O3'])
self.cosZen = 0.5 # cosine of the average solar zenith angle
self._define_diagnostics()
@property
def emissivity(self):
# This ensures that emissivity is always zero for shortwave classes
return np.zeros_like(self.absorptivity)
[docs]
class FourBandLW(NbandRadiation):
def __init__(self, **kwargs):
r'''Closely following SPEEDY / MITgcm longwave model
band 0 is window region (between 8.5 and 11 microns)
band 1 is CO2 channel (the band of strong absorption by CO2 around 15 microns)
band 2 is weak H2O channel (aggregation of spectral regions with weak to moderate absorption by H2O)
band 3 is strong H2O channel (aggregation of regions with strong absorption by H2O)
'''
super(FourBandLW, self).__init__(**kwargs)
# SPEEDY uses an approximation to the Planck function
# and the band fraction for every emission is calculated from
# its current temperature
# Here for simoplicity we'll just set an average band_fraction
# and hold it fixed
Tarray = np.linspace(-30, 30) + 273.15
self.band_fraction = np.mean(SPEEDY_band_fraction(Tarray), axis=1)
# defaults from MITgcm/aim:
# these are layer absorptivities per dp = 10^5 Pa
# the water vapor terms are expressed for dq = 1 g/kg
ABLWIN = 0.7
ABLCO2 = 4.0
ABLWV1 = 0.7
ABLWV2 = 50.0
# I'm going to assume that the absorption in window region is by O3.
# not sure if this makes any sense... maybe it should be zero
O3 = np.array([ABLWIN, 0., 0., 0.]) / 1E5 * const.g / 5E-6
self.absorption_cross_section['O3'] = np.reshape(O3,
(self.num_channels, 1))
# the CO2 mixing ratio for which SPEEDY / MITgcm is tuned...
# not clear what this number should be
AIMCO2 = 380E-6
CO2 = np.array([0., ABLCO2, 0., 0.]) / 1E5 * const.g / AIMCO2
#self.absorption_cross_section['CO2'] = np.reshape(CO2,
# (self.num_channels, 1))
# Need to multiply by 1E3 for H2O fields because we use kg/kg for mixing ratio
H2O = np.array([0., 0., ABLWV1, ABLWV2]) / 1E5 * const.g * 1E3
#self.absorption_cross_section['H2O'] = np.reshape(H2O,
# (self.num_channels, 1))
for n in range(self.Tatm.domain.numdims):
CO2 = CO2[:, np.newaxis]
O3 = O3[:, np.newaxis]
H2O = H2O[:, np.newaxis]
self.absorption_cross_section.update({'CO2': CO2, 'H2O': H2O, 'O3': O3})
if 'CO2' not in self.absorber_vmr:
self.absorber_vmr['CO2'] = 380.E-6 * np.ones_like(self.Tatm)
if 'O3' not in self.absorber_vmr:
self.absorber_vmr['O3'] = np.zeros_like(self.Tatm)
if 'H2O' not in self.absorber_vmr:
self.absorber_vmr['H2O'] = self.q
self._define_diagnostics()
[docs]
def SPEEDY_band_fraction(T):
r'''Python / numpy implementation of the formula used by SPEEDY and MITgcm
to partition longwave emissions into 4 spectral bands.
Input: temperature in Kelvin
returns: a four-element array of band fraction
Reproducing here the FORTRAN code from MITgcm/pkg/aim_v23/phy_radiat.F
.. code-block:: fortran
EPS3=0.95 _d 0
DO JTEMP=200,320
FBAND(JTEMP,0)= EPSLW
FBAND(JTEMP,2)= 0.148 _d 0 - 3.0 _d -6 *(JTEMP-247)**2
FBAND(JTEMP,3)=(0.375 _d 0 - 5.5 _d -6 *(JTEMP-282)**2)*EPS3
FBAND(JTEMP,4)= 0.314 _d 0 + 1.0 _d -5 *(JTEMP-315)**2
FBAND(JTEMP,1)= 1. _d 0 -(FBAND(JTEMP,0)+FBAND(JTEMP,2)
& +FBAND(JTEMP,3)+FBAND(JTEMP,4))
ENDDO
DO JB=0,NBAND
DO JTEMP=lwTemp1,199
FBAND(JTEMP,JB)=FBAND(200,JB)
ENDDO
DO JTEMP=321,lwTemp2
FBAND(JTEMP,JB)=FBAND(320,JB)
ENDDO
ENDDO
'''
# EPSLW is the fraction of longwave emission that goes directly to space
# It is set to zero by default in MITgcm code. We won't use it here.
Tarray = np.array(T)
Tarray = np.minimum(Tarray, 320.)
Tarray = np.maximum(Tarray, 200.)
num_band = 4
dims = [num_band]
dims.extend(Tarray.shape)
FBAND = np.zeros(dims)
EPS2=0.95
FBAND[1,:] = 0.148 - 3.0E-6 *(T-247.)**2
FBAND[2,:] = (0.375 - 5.5E-6 *(T-282.)**2)*EPS2
FBAND[3,:] = 0.314 + 1.0E-5 *(T-315.)**2
FBAND[0,:] = 1. - np.sum(FBAND, axis=0)
return FBAND