Merge branch 'v1.4'

docs/examples.rst

@@ -40,7 +40,9 @@ in which CO2 is the only emitted species.
 In almost every application of FaIR you will probably want to vary the
 ``emissions`` time series going in to ``fair_scm``. In CO2-only mode
 this is a 1D array of CO2 emissions. Setting ``useMultigas=False`` turns
-off the emissions from non-CO2 species.
+off the emissions from non-CO2 species. To include radiative forcing
+from non-CO2 forcers, you can specify a constant or a timeseries to the
+``other_rf`` keyword as in the below example.
 
 The output from FaIR is a 3-tuple of ``(C,F,T)`` arrays. In CO2 mode,
 both ``C`` (representing CO2 concentrations in ppm) and ``F`` (total
@@ -349,7 +351,7 @@ upper limit on the time-integrated airborne fraction.
 
 .. parsed-literal::
 
-    /nfs/see-fs-02_users/mencsm/FAIR/fair/forward.py:233: RuntimeWarning: iirf_h=60.000000, which is less than iirf_max (97.000000)
+    /nfs/see-fs-02_users/mencsm/FAIR/fair/forward.py:234: RuntimeWarning: iirf_h=60.000000, which is less than iirf_max (97.000000)
       % (iirf_h, iirf_max), RuntimeWarning)
 
 
@@ -905,7 +907,7 @@ for illustration). Note this is a completely hypothetical scenario!
 
 .. parsed-literal::
 
-    <matplotlib.text.Text at 0x7f8551658890>
+    <matplotlib.text.Text at 0x7fc03b049810>
 
 
 

fair/inverse.py

@@ -60,6 +60,7 @@ def inverse_carbon_cycle(c1, c_acc0, temp, r0, rc, rt, iirf_max, time_scale_sf,
         time_scale_sf : scale factor for CO2 decay constants
     """
 
+
     iirf = iirf_simple(c_acc0, temp, r0, rc, rt, iirf_max)
     time_scale_sf = root(iirf_interp, time_scale_sf,
       args=(a, tau, iirf_h, iirf))['x']
@@ -88,8 +89,55 @@ def inverse_fair_scm(
     iirf_max      = carbon.iirf_max,
     iirf_h        = carbon.iirf_h,
     C_pi          = 278.,
-    time_scale_sf = 0.16
+    time_scale_sf = 0.16,
+    restart_in    = False,
+    restart_out   = False,
     ):
+
+    """Diagnoses emissions from prescribed concentrations.
+
+    Inputs:
+        C             : concentrations of CO2, ppmv
+        other_rf      : non-CO2 radiative forcing (scalar or numpy array, W/m2)
+        q             : coefficients of slow and fast temperature change.
+                        Overridden if tcrecs is specified.
+        tcrecs        : transient climate response and equilibrium climate
+                        sensitivity array (2-element or (nt, 2))
+        d             : timescales of slow and fast contribution to temperature
+                        change
+        F2x           : radiative forcing from a doubling of CO2 concentrations
+                        (W/m2)
+        tcr_dbl       : timescale over which a 1% compound increase of CO2 acts
+                        (yr)
+        a             : partition fractions for CO2 boxes
+        tau           : time constants for CO2 boxes
+        r0            : pre-industrial time-integrated airborne fraction (yr)
+        rc            : sensitivity of time-integrated airborne fraction to
+                        airborne carbon (yr/GtC)
+        rt            : sensitivity of time-integrated airborne fraction to
+                        temperature (yr/K)
+        iirf_max      : maximum value of time-integrated airborne fraction (yr)
+        iirf_h        : time horizon for time-integrated airborne fraction
+        C_pi          : pre-industrial concentration of CO2, ppmv
+        time_scale_sf : initial guess for scaling factor of CO2 time constants.
+                        Overridden if using a restart.
+        restart_in    : Allows a restart of the carbon cycle from a non-initial
+                        state. A 6-tuple of:
+                          array of accumulated carbon in each atmospheric box,
+                          array of slow and fast temperature contributions,
+                          total accumulated carbon,
+                          emissions in the timestep before restart
+                          time constant scale factor
+                          CO2 concentrations in the timestep before restart
+        restart_out   : if True, return the restart state as an extra output.
+                        See restart_in.
+    Outputs:
+        E             : Timeseries of diagnosed CO2 emissions in GtC
+        F             : Timeseries of total radiative forcing, W/m2
+        T             : Timeseries of temperature anomaly since pre-industrial
+        restart       : if restart_out=True, 6-tuple of carbon cycle state
+                        parameters. See restart_in.
+    """
 
     # Error checking and validation goes here...
 
@@ -99,7 +147,8 @@ def inverse_fair_scm(
     thermal_boxes_shape = (nt, d.shape[0])
 
     # Thermal response
-    q         = calculate_q(tcrecs, d, F2x, tcr_dbl, nt)
+    if type(tcrecs) is np.ndarray:
+        q     = calculate_q(tcrecs, d, F2x, tcr_dbl, nt)
 
     # Allocate intermediate and output arrays
     C_acc     = np.zeros(nt)
@@ -112,10 +161,27 @@ def inverse_fair_scm(
         other_rf = other_rf * np.ones(nt)
 
     # First timestep
-    emissions[0] = root(infer_emissions, 0., args=(C[0], R_i[0],
-        tau, a, C_pi))['x']
-    F[0]         = co2_log(C[0], C_pi, F2x=F2x)# + other_rf[0]
-    T_j[0,:]     = forc_to_temp(T_j[0,:], q[0,:], d, F[0])
+    if restart_in:
+        R_i_minus1    = restart_in[0]
+        T_j_minus1    = restart_in[1]
+        C_acc_minus1  = restart_in[2]
+        E_minus1      = restart_in[3]
+        time_scale_sf = restart_in[4]
+        C_minus1      = restart_in[5]
+        emissions[0], C_acc[0], R_i[0,:], time_scale_sf = (
+            inverse_carbon_cycle(
+                C[0], C_acc_minus1, np.sum(T_j_minus1), r0, rc, rt, iirf_max,
+                time_scale_sf, a, tau, iirf_h, R_i_minus1,
+                C_pi, C_minus1, E_minus1
+            )
+        )
+        F[0]          = co2_log(C[0], C_pi, F2x=F2x) + other_rf[0]
+        T_j[0,:]      = forc_to_temp(T_j_minus1, q[0,:], d, F[0])
+    else:
+        emissions[0]  = root(infer_emissions, 0., args=(C[0], R_i[0,:],
+            tau, a, C_pi))['x']
+        F[0]          = co2_log(C[0], C_pi, F2x=F2x) + other_rf[0]
+        T_j[0,:]      = forc_to_temp(T_j[0,:], q[0,:], d, F[0])
 
     # Second timestep onwards
     for t in range(1,nt):
@@ -131,5 +197,9 @@ def inverse_fair_scm(
 
     # Output temperatures
     T = np.sum(T_j, axis=-1)
-
-    return emissions, F, T
+    if restart_out:
+        restart_out_val = (R_i[-1], T_j[-1], C_acc[-1], emissions[-1],
+            time_scale_sf, C[-1])
+        return emissions, F, T, restart_out_val
+    else:
+        return emissions, F, T

fair/tools/constrain.py

@@ -1,3 +1,5 @@
+from __future__ import division
+
 import numpy as np
 from scipy import stats
 

fair/tools/ensemble.py

@@ -1,3 +1,5 @@
+from __future__ import division
+
 import numpy as np
 import scipy.stats as st
 import os
@@ -83,6 +85,7 @@ def tcrecs_generate(tcrecs_in='cmip5', dist='lognorm', n=1000, correlated=True,
         tcrecs_in = np.loadtxt(filepath, delimiter=',', skiprows=3)
     try:
         assert(type(tcrecs_in) is np.ndarray)
+        assert(tcrecs_in.ndim == 2)
         assert(tcrecs_in.shape[1] == 2)
     except AssertionError:
         raise ValueError('tcrecs_in should "cmip5" or an array of shape (n, 2)')

fair/tools/gwp.py

@@ -1,3 +1,5 @@
+from __future__ import division
+
 import numpy as np
 from ..constants import molwt, radeff
 from ..constants.general import M_ATMOS

fair/tools/magicc.py

@@ -1,3 +1,5 @@
+from __future__ import division
+
 import numpy as np
 from scipy.interpolate import interp1d
 

tests/integration/integration_test.py

@@ -3,6 +3,7 @@
 import fair
 from fair.RCPs import rcp26, rcp3pd, rcp45, rcp6, rcp60, rcp85
 import numpy as np
+import warnings
 from fair.forcing.ghg import myhre
 from fair.constants import molwt
 from copy import deepcopy
@@ -72,23 +73,55 @@ def test_contrails_invalid():
         fair.forward.fair_scm(emissions=rcp45.Emissions.emissions,
             contrail_forcing='other')
 
+
 def test_aerosol_regression_zeros_fix():
     _, F, _ = fair.forward.fair_scm(
         emissions=fair.RCPs.rcp85.Emissions.emissions,
-        b_aero = np.zeros(7)
+        b_aero = np.zeros(7),
+        aerosol_forcing='aerocom'
     )
-    # Index 7 is aerosol forcing
-    assert (F[:,7]==np.zeros(736)).all()
+    # Index 8 is aerosol forcing
+    assert (F[:,8]==np.zeros(736)).all()
 
 
 def test_aerosol_regression_zeros_nofix():
     _, F, _ = fair.forward.fair_scm(
         emissions=fair.RCPs.rcp85.Emissions.emissions,
         b_aero = np.zeros(7),
-        fixPre1850RCP=False
+        fixPre1850RCP=False,
+        aerosol_forcing='aerocom'
     )
     # Index 7 is aerosol forcing
-    assert (F[:,7]==np.zeros(736)).all()
+    assert (F[:,8]==np.zeros(736)).all()
+
+
+def test_aerosol_noscale():
+    _, F1, _ = fair.forward.fair_scm(
+        emissions=fair.RCPs.rcp85.Emissions.emissions,
+        scaleAerosolAR5=False
+    )
+    _, F2, _ = fair.forward.fair_scm(
+        emissions=fair.RCPs.rcp85.Emissions.emissions,
+        scaleAerosolAR5=True
+    )
+    assert (F1[:,8]!=F2[:,8]).any()
+
+
+def test_stevens():
+    # The Stevens aerosol causes an overflow in the root finder routine at
+    # present. This isn't a problem, but switching it off gives nice clean
+    # tests
+    with warnings.catch_warnings():
+        warnings.simplefilter('ignore')
+        _, F1, _ = fair.forward.fair_scm(
+            emissions=fair.RCPs.rcp85.Emissions.emissions,
+            aerosol_forcing='Stevens'
+        )
+    _, F2, _ = fair.forward.fair_scm(
+        emissions=fair.RCPs.rcp85.Emissions.emissions,
+        aerosol_forcing='aerocom+ghan'
+    )
+    assert (F1[:,8]!=F2[:,8]).any()
 
 
 def test_ozone_regression_zero():

tests/reproduction/reproduction_test.py

@@ -4,7 +4,9 @@
 from fair.RCPs import rcp3pd, rcp45, rcp6, rcp85, rcp26, rcp60
 import numpy as np
 import os
-from fair.constants import molwt
+from fair.constants import molwt, radeff, lifetime
+from fair.tools.constrain import hist_temp
+from fair.tools.gwp import gwp
 
 def test_ten_GtC_pulse():
     emissions = np.zeros(250)
@@ -193,3 +195,90 @@ def test_forward_versus_reverse():
     assert np.allclose(E_forward, E_inverse, atol=0.01, rtol=0.01)
     assert np.allclose(F_forward, F_inverse, atol=0.01, rtol=0.01)
     assert np.allclose(T_forward, T_inverse, atol=0.01, rtol=0.01)
+
+
+def test_restart_co2_continuous():
+    """Tests to check that a CO2-only run with a restart produces the same
+    results as a CO2-only run without a restart."""
+
+    C, F, T = fair.forward.fair_scm(
+        emissions   = rcp45.Emissions.co2[:20],
+        useMultigas = False
+        )
+
+    C1, F1, T1, restart = fair.forward.fair_scm(
+        emissions   = rcp45.Emissions.co2[:10],
+        useMultigas = False,
+        restart_out = True
+        )
+
+    C2, F2, T2 = fair.forward.fair_scm(
+        emissions   = rcp45.Emissions.co2[10:20],
+        useMultigas = False,
+        restart_in  = restart
+        )
+
+    assert np.all(C == np.concatenate((C1, C2)))
+    assert np.all(F == np.concatenate((F1, F2)))
+    assert np.all(T == np.concatenate((T1, T2)))
+
+
+def test_inverse_restart():
+    """Tests restarts for inverse FaIR."""
+
+    E, F, T = fair.inverse.inverse_fair_scm(
+        C = rcp85.Concentrations.co2[:20])
+
+    E1, F1, T1, restart = fair.inverse.inverse_fair_scm(
+        C = rcp85.Concentrations.co2[:10], restart_out=True)
+
+    E2, F2, T2 = fair.inverse.inverse_fair_scm(
+        C = rcp85.Concentrations.co2[10:20], restart_in=restart)
+
+    assert np.all(E == np.concatenate((E1, E2)))
+    assert np.all(F == np.concatenate((F1, F2)))
+    assert np.all(T == np.concatenate((T1, T2)))
+
+
+def test_constrain():
+    """Checks that the historical temperature constraining function works"""
+
+    datadir = os.path.join(os.path.dirname(__file__),
+        '../../fair/tools/tempobs/')
+    tempobsdata = np.loadtxt(datadir+'had4_krig_annual_v2_0_0.csv')
+    years   = tempobsdata[:,0]
+    tempobs = tempobsdata[:,1]
+
+    C,F,T = fair.forward.fair_scm(emissions=rcp45.Emissions.emissions)
+    accept1,sm1,im1,so1,io1 = hist_temp(tempobs, T[85:252], years)
+    assert accept1==True
+
+    accept2,sm2,im2,so2,io2 = hist_temp(tempobs, T[85:252], years,
+        inflate=False)
+    assert sm1==sm2
+    assert so1==so2
+    assert accept2==True
+
+    accept3,_,_,_,_ = hist_temp(tempobs, np.zeros(167), years)
+    assert accept3==False
+
+
+def test_gwp():
+    """Checks that GWP calculator produces correct GWPs."""
+
+    # methane uses "perturbation lifetime" for GWP calculations and feedback
+    # factor
+    assert np.round(gwp(100, 12.4, radeff.CH4, molwt.CH4, f=0.65))==28
+
+    # for N2O, I think the IPCC AR5 value is out by one year. Most likely
+    # explanation is that they have rounded off the feedback somewhere.
+    # This is calculated as 1-(1-0.36*(1.65)*radeff.CH4/radeff.N2O). See
+    # eq. 8.SM.20 in the supplement to Chapter 8, AR5
+    assert np.round(
+        gwp(20, lifetime.N2O, radeff.N2O, molwt.N2O, f=-0.071874))==263
+    assert np.round(
+        gwp(100, lifetime.N2O, radeff.N2O, molwt.N2O, f=-0.071874))==264
+
+    # Now check a nice straightforward example
+    assert np.round(
+        gwp(100, lifetime.CFC11, radeff.CFC11, molwt.CFC11), decimals=-1)==4660