diff --git a/fair/constants/preindconc.py b/fair/constants/preindconc.py new file mode 100644 index 00000000..30f5b4dd --- /dev/null +++ b/fair/constants/preindconc.py @@ -0,0 +1,38 @@ +# pre-industrial GHG concentrations in their default units. +# Source: http://www.pik-potsdam.de/~mmalte/rcps. Year 1765 in all cases. +CO2 = 278.05158 # ppm +CH4 = 721.89411 # ppb +N2O = 272.95961 # ppb +CF4 = 35. # ppt +C2F6 = 0. # ppt +C6F14 = 0. # ppt +HFC23 = 0. # ppt +HFC32 = 0. # ppt +HFC43_10 = 0. # ppt +HFC125 = 0. # ppt +HFC134A = 0. # ppt +HFC143A = 0. # ppt +HFC227EA = 0. # ppt +HFC245FA = 0. # ppt +SF6 = 0. # ppt +CFC11 = 0. # ppt +CFC12 = 0. # ppt +CFC113 = 0. # ppt +CFC114 = 0. # ppt +CFC115 = 0. # ppt +CARB_TET = 0. # ppt +MCF = 0. # ppt +HCFC22 = 0. # ppt +HCFC141B = 0. # ppt +HCFC142B = 0. # ppt +HALON1211= 0. # ppt +HALON1202= 0. # ppt +HALON1301= 0. # ppt +HALON2402= 0. # ppt +CH3BR = 5.8 # ppt +CH3CL = 480. # ppt + +aslist = [CO2, CH4, N2O, CF4, C2F6, C6F14, HFC23, HFC32, HFC43_10, HFC125, + HFC134A, HFC143A, HFC227EA, HFC245FA, SF6, CFC11, CFC12, CFC113, + CFC114, CFC115, CARB_TET, MCF, HCFC22, HCFC141B, HCFC142B, + HALON1211, HALON1202, HALON1301, HALON2402, CH3BR, CH3CL] diff --git a/fair/tools/steady.py b/fair/tools/steady.py new file mode 100644 index 00000000..1f1ece4e --- /dev/null +++ b/fair/tools/steady.py @@ -0,0 +1,79 @@ +from __future__ import division + +import numpy as np +from itertools import compress +from ..constants.general import M_ATMOS +from ..constants import molwt as molwt_builtin, preindconc,\ + lifetime as lt_builtin + + +def _lookup(species): + # get list of GHGs from preindconc model + named_gases = dir(preindconc) + + # strip out reserved names and the "aslist" list + include = [i[:2]!='__' and i!='aslist' for i in named_gases] + named_gases = list(compress(named_gases, include)) + + # see if the species appears on the list of builtin gases and return + # properties if so, otherwise raise error + if species in named_gases: + inout = {} + inout['pi'] = preindconc + inout['lt'] = lt_builtin + inout['mw'] = molwt_builtin + exec("C, L, M = pi."+species+", lt."+species+", mw."+species, inout) + return inout['C'], inout['L'], inout['M'] + else: + raise ValueError(species + ' is not in the list of recognised '+ + 'greenhouse gases') + + +def emissions(C=None, lifetime=None, molwt=None, species=None): + """Calculate steady state background emissions from a given lifetime + + Keywords: + C: concentrations, scalar or None + if scalar, steady state concentrations to acheive. + if None, use the pre-industrial concentrations for the specified + gas. + lifetime: scalar or None + if scalar, use the given greenhouse gas atmospheric lifetime + if None, use the default lifetime for the specified gas + molwt: scalar or None + if scalar, use molecular weight of given species + if None, use the default molecular weight for the specified gas + species: string or None + Name of the greenhouse gas concentrations you want to calculate. + See ..constants.lifetime module for used names. + If None, use the values given in C, lifetime and wt. + + Any user-specified values for C, lifetime and molwt overrides the + defaults. They must be specified if species is None. + + Returns: + emissions (scalar) + """ + + # if not using a built-in gas, C, lifetime and molwt must be specified. + if species is None: + if any((C is None, lifetime is None, molwt is None)): + raise ValueError('If species is not given then C, '+ + 'lifetime and molwt must be specified.') + else: + # populate defaults but override if values given + species = species.upper() + C0, lifetime0, molwt0 = _lookup(species) + if C is None: C=C0 + if lifetime is None: lifetime=lifetime0 + if molwt is None: molwt=molwt0 + + # convert units for N2O + if species == 'N2O': + molwt = molwt * molwt_builtin.N2/molwt_builtin.N2O + + # now invert the emissions-concentrations relationship + E = (C * (1.0 - np.exp(-1.0 / lifetime)) * M_ATMOS * molwt / + molwt_builtin.AIR) * 1e-18 + + return E diff --git a/tests/test_fair.py b/tests/test_fair.py index f436e843..fce404dc 100644 --- a/tests/test_fair.py +++ b/tests/test_fair.py @@ -4,12 +4,11 @@ import os import numpy as np from fair.RCPs import rcp3pd, rcp45, rcp6, rcp85 -from fair.tools import magicc +from fair.tools import magicc, steady from fair.ancil import natural, cmip5_annex2_forcing from fair.constants import molwt from fair.forcing.ghg import myhre - def test_no_arguments(): with pytest.raises(ValueError): fair.forward.fair_scm() @@ -360,3 +359,23 @@ def test_contrails_invalid(): contrail_forcing='other') +def test_steady(): + assert np.isclose(steady.emissions(species='CH4'), 209.2492053169677) + assert np.isclose(steady.emissions(species='N2O'), 11.155476818389447) + assert np.isclose(steady.emissions(species='CF4'), 0.010919593149304725) + assert steady.emissions(species='CFC11')==0. + assert np.isclose(steady.emissions(species='CH4', lifetime=12.4), + 159.0269945578832) + assert np.isclose(steady.emissions(species='CH4', molwt=17), + 221.77284852795833) + assert np.isclose(steady.emissions(species='CH4', C=1750.), + 507.2573722823332) + # verify override working + assert steady.emissions(species='CH4', lifetime=12.4, molwt=17, + C=1750.) != steady.emissions(species='CH4') + # not a gas on the list - should report error + with pytest.raises(ValueError): + steady.emissions(species='chocolate') + # and no input equals no output, so again value error + with pytest.raises(ValueError): + steady.emissions()