Merge branch 'fix-45'

CHANGELOG.rst

@@ -4,6 +4,8 @@ Changelog
 master
 ------
 
+(`#84 <https://github.com/OMS-NetZero/FAIR/pull/84>`_) Fix concentration-driven runs for 45-species FaIR
+
 (`#83 <https://github.com/OMS-NetZero/FAIR/pull/83>`_) Support scmdata >= 0.7.1
 
 (`#82 <https://github.com/OMS-NetZero/FAIR/pull/82>`_) Expand AR6 forcing diagnostics to get aerosol direct forcing by species

fair/forward.py

@@ -89,7 +89,6 @@ def fair_scm(
     stevens_params = np.array([0.001875, 0.634, 60.]),
     ref_isSO2=True, # is Stevens SO2 emissions in units SO2 (T) or S (F)
     useMultigas=True,
-    useStevenson=True,   # deprecate this switch in v1.7
     tropO3_forcing='stevenson',
     lifetimes=False,
     aerosol_forcing="aerocom+ghan",
@@ -116,11 +115,6 @@ def fair_scm(
     if type(emissions) is bool and not emissions_driven:
         tropO3_forcing='external'
 
-    if useStevenson is not None:
-        warnings.warn('"useStevenson" will be deprecated in the future; use '+
-          'tropO3_forcing keyword with "cmip6", "stevenson", "regression" or "external"',
-          DeprecationWarning)
-
     # is iirf_h < iirf_max? Don't stop the code, but warn user
     if iirf_h < iirf_max:
         warnings.warn('iirf_h=%f, which is less than iirf_max (%f)'
@@ -447,7 +441,7 @@ def fair_scm(
         # because SLCFs can still be given as emissions with GHGs as
         # concentrations
         if type(emissions) is not bool:
-            if useStevenson and tropO3_forcing[0].lower()=='s':
+            if tropO3_forcing[0].lower()=='s':
                 F[0,iF_tro3] = ozone_tr.stevenson(emissions[0,:], C[0,1],
                   T=np.sum(T_j[0,:]), 
                   feedback=useTropO3TFeedback,
@@ -459,7 +453,7 @@ def fair_scm(
                   feedback=useTropO3TFeedback,
                   PI=np.array([C_pi[1],E_pi[6],E_pi[7],E_pi[8]]),
                   beta=b_tro3)
-            elif not useStevenson or tropO3_forcing[0].lower()=='r':
+            elif tropO3_forcing[0].lower()=='r':
                 F[0,iF_tro3] = ozone_tr.regress(emissions[0,:]-E_pi[:], beta=b_tro3)
             else:
                 F[0,iF_tro3] = F_tropO3[0]
@@ -702,7 +696,7 @@ def fair_scm(
                 else:
                     F[t,3] = np.sum(minor_gases(C[t,3:], C_pi[3:]))
 
-                if useStevenson and tropO3_forcing[0].lower()=='s':
+                if tropO3_forcing[0].lower()=='s':
                     F[t,iF_tro3] = ozone_tr.stevenson(emissions[t,:],
                       C[t,1],
                       T=T[t-1], 
@@ -715,7 +709,7 @@ def fair_scm(
                       feedback=useTropO3TFeedback,
                       PI=np.array([C_pi[1],E_pi[6],E_pi[7],E_pi[8]]),
                       beta=b_tro3)
-                elif not useStevenson or tropO3_forcing[0].lower()=='r':
+                elif tropO3_forcing[0].lower()=='r':
                     F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3)
                 else:
                     F[t,iF_tro3] = F_tropO3[t]
@@ -811,29 +805,31 @@ def fair_scm(
 
             if useMultigas:
                 F[t,0:3] = ghg(C[t,0:3], C_pi[0:3], F2x=F2x)
-                F[t,3] = np.sum((C[t,3:] - C_pi[3:]) * radeff.aslist[3:]
-                  * 0.001)
+                if diagnostics=='AR6':
+                    F[t,3:31] = minor_gases(C[t,3:], C_pi[3:])
+                else:
+                    F[t,3] = np.sum(minor_gases(C[t,3:], C_pi[3:]))
                 if type(emissions) is not bool:
-                    if useStevenson and tropO3_forcing[0].lower()=='s':
-                        F[t,4] = ozone_tr.stevenson(emissions[t,:]-E_pi,
+                    if tropO3_forcing[0].lower()=='s':
+                        F[t,iF_tro3] = ozone_tr.stevenson(emissions[t,:]-E_pi,
                           C[t,1],
                           T=T[t-1],
                           feedback=useTropO3TFeedback,
                           fix_pre1850_RCP=fixPre1850RCP)
                     elif tropO3_forcing[0].lower()=='c':
-                        F[t,4] = ozone_tr.cmip6_stevenson(emissions[t,:], C[t,1],
+                        F[t,iF_tro3] = ozone_tr.cmip6_stevenson(emissions[t,:], C[t,1],
                           T=np.sum(T_j[t,:]),
                           feedback=useTropO3TFeedback,
                           PI=np.array([C_pi[1],E_pi[6],E_pi[7],E_pi[8]]),
                           beta=b_tro3)
-                    elif not useStevenson or tropO3_forcing[0].lower()=='r':
-                        F[t,4] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3)
+                    elif tropO3_forcing[0].lower()=='r':
+                        F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3)
                     else:
-                        F[t,4] = F_tropO3[t]
+                        F[t,iF_tro3] = F_tropO3[t]
                 else:
-                    F[t,4] = F_tropO3[t]
-                F[t,5] = ozone_st.magicc(C[t,15:], C_pi[15:])
-                F[t,6] = h2o_st.linear(F[t,1], ratio=stwv_from_ch4)
+                    F[t,iF_tro3] = F_tropO3[t]
+                F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:])
+                F[t,iF_ch4h] = h2o_st.linear(F[t,1], ratio=stwv_from_ch4)
 
                 # multiply by scale factors
                 F[t,:] = F[t,:] * scale[t,:]

tests/integration/integration_test.py

@@ -127,7 +127,7 @@ def test_stevens():
 def test_ozone_regression_zero():
     _, F, _ = fair.forward.fair_scm(
         emissions=fair.RCPs.rcp85.Emissions.emissions,
-        useStevenson=False,
+        tropO3_forcing='regression',
         b_tro3 = np.zeros(4)
     )
     # Index 4 is ozone forcing
@@ -268,6 +268,24 @@ def test_45():
     assert np.allclose(F45.sum(axis=1), F13.sum(axis=1))
 
 
+def test_45_conc_driven():
+    """Test 45-species forcing output when driven with concentrations"""
+    _, F45, _ = fair.forward.fair_scm(
+        emissions = rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C = rcp85.Concentrations.gases,
+        efficacy = np.ones(45),
+        diagnostics = 'AR6'
+    )
+    _, F13, _ = fair.forward.fair_scm(
+        emissions = rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C = rcp85.Concentrations.gases,
+        efficacy = np.ones(13)
+    )
+    assert np.allclose(F45.sum(axis=1), F13.sum(axis=1))
+
+
 def test_geoffroy():
     C, F, T, lambda_eff, ohc, heatflux = fair.forward.fair_scm(
         emissions = rcp85.Emissions.emissions,

tests/reproduction/reproduction_test.py

@@ -2,7 +2,6 @@
 
 import fair
 from fair.RCPs import rcp3pd, rcp45, rcp6, rcp85, rcp26, rcp60
-#from fair.SSPs import historical, ssp119, ssp126, ssp245, ssp370, ssp434, ssp460, ssp534over, ssp585
 import numpy as np
 import os
 from fair.constants import molwt, radeff, lifetime

tests/unit/unit_test.py

@@ -391,13 +391,6 @@ def test_ozone_regression_equivalence():
     assert F1==F2[100]
 
 
-def test_deprecated_ozone_stevenson():
-    with pytest.warns(DeprecationWarning):
-        C,F,T = fair.forward.fair_scm(
-          emissions=rcp85.Emissions.emissions,
-          useStevenson=True
-        )
-
 def test_cmip6_stevenson():
     C1,F1,T1 = fair.forward.fair_scm(
       emissions=rcp85.Emissions.emissions,
@@ -416,12 +409,14 @@ def test_cmip6_stevenson():
     # check differences
     assert np.any(F2[:,4]!=F1[:,4])
 
+
 def test_gir():
     C,F,T = fair.forward.fair_scm(
       emissions=rcp85.Emissions.emissions,
       gir_carbon_cycle=True
     )
 
+
 def test_meinshausen():
     C1,F1,T1 = fair.forward.fair_scm(
       emissions=rcp85.Emissions.emissions,
@@ -432,3 +427,63 @@ def test_meinshausen():
       ghg_forcing="Etminan"
     )
     assert np.any(F1!=F2)
+
+
+def test_ozone_treatments():
+    C1,F1,T1 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        tropO3_forcing='stevenson'
+    )
+    C2,F2,T2 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        tropO3_forcing='cmip6'
+    )
+    C3,F3,T3 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        tropO3_forcing='regress'
+    )
+    C4,F4,T4 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        tropO3_forcing='external',
+        F_tropO3 = 0
+    )
+    assert np.any(F1!=F2)
+    assert np.any(F1!=F3)
+    assert np.any(F1!=F4)
+    assert np.any(F2!=F3)
+    assert np.any(F2!=F4)
+    assert np.any(F3!=F4)
+
+
+def test_ozone_treatments_conc_driven():
+    C1,F1,T1 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C=rcp85.Concentrations.gases,
+        tropO3_forcing='stevenson'
+    )
+    C2,F2,T2 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C=rcp85.Concentrations.gases,
+        tropO3_forcing='cmip6'
+    )
+    C3,F3,T3 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C=rcp85.Concentrations.gases,
+        tropO3_forcing='regress'
+    )
+    C4,F4,T4 = fair.forward.fair_scm(
+        emissions=rcp85.Emissions.emissions,
+        emissions_driven=False,
+        C=rcp85.Concentrations.gases,
+        tropO3_forcing='external',
+        F_tropO3 = 0
+    )
+    assert np.any(F1!=F2)
+    assert np.any(F1!=F3)
+    assert np.any(F1!=F4)
+    assert np.any(F2!=F3)
+    assert np.any(F2!=F4)
+    assert np.any(F3!=F4)