Merge pull request #34 from OMS-NetZero/equilibriate_natural

Steady state emissions from concentrations tool

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]

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

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()