Merge branch 'patch-aerocom-ghan'

fair/forcing/aerosols.py

@@ -33,16 +33,17 @@ def Stevens(emissions, stevens_params=np.array([0.001875, 0.634, 60.]),
     em_pi  = E_pi * factor
 
     ERFari = -alpha * (em_SOx-em_pi)
-    ERFaci = (
-        (-beta * np.log(em_SOx/E_SOx_nat + 1)) - 
-        (-beta * np.log(em_pi/E_SOx_nat + 1)) )
+#    ERFaci = (
+#        (-beta * np.log(em_SOx/E_SOx_nat + 1)) - 
+#        (-beta * np.log(em_pi/E_SOx_nat + 1)) )
+    ERFaci = (-beta * np.log((em_SOx-em_pi)/E_SOx_nat + 1))
     return ERFari, ERFaci
 
 
 def aerocom_direct(emissions,
         beta = np.array(
           [-6.2227e-3, 0.0, -3.8392e-4, -1.16551e-3, 1.601537e-2, -1.45339e-3,
-           -1.55605e-3])
+           -1.55605e-3]), E_pi=np.zeros(40)
     ):
 
     """Calculates direct aerosol forcing based on linear relationships between
@@ -57,28 +58,30 @@ def aerocom_direct(emissions,
     Keywords:
         beta: 7-element array of forcing efficiencies in W m-2 (Mt yr-1)-1 for
             SOx, CO, NMVOC, NOx, BC, OC, NH3 (in that order)
+        E_pi: pre-industrial emissions (40 element array)
     Outputs:
         Forcing time series
     """
 
     em_SOx, em_CO, em_NMVOC, em_NOx, em_BC, em_OC, em_NH3 = \
         emissions[:,[5, 6, 7, 8, 9, 10, 11]].T
 
-    F_SOx    = beta[0] * em_SOx
-    F_CO     = beta[1] * em_CO
-    F_NMVOC  = beta[2] * em_NMVOC
-    F_NOx    = beta[3] * em_NOx
-    F_BC     = beta[4] * em_BC
-    F_OC     = beta[5] * em_OC
-    F_NH3    = beta[6] * em_NH3
+    F_SOx    = beta[0] * (em_SOx   - E_pi[5])
+    F_CO     = beta[1] * (em_CO    - E_pi[6])
+    F_NMVOC  = beta[2] * (em_NMVOC - E_pi[7])
+    F_NOx    = beta[3] * (em_NOx   - E_pi[8])
+    F_BC     = beta[4] * (em_BC    - E_pi[9])
+    F_OC     = beta[5] * (em_OC    - E_pi[10])
+    F_NH3    = beta[6] * (em_NH3   - E_pi[11])
 
     ERFari = F_SOx+F_CO+F_NMVOC+F_NOx+F_BC+F_OC+F_NH3
 
     return ERFari
 
 
 def ghan_indirect(emissions, fix_pre1850_RCP=True, scale_AR5=False,
-    ghan_params=np.array([-1.95011431, 0.01107147, 0.01387492])):
+    ghan_params=np.array([-1.95011431, 0.01107147, 0.01387492]),
+    E_pi=np.zeros(40)):
     """Estimates the aerosol indirect effect based on the simple model in
     Ghan et al., (2013), doi:10.1002/jgrd.50567.
 

fair/forward.py

@@ -205,6 +205,9 @@ def fair_scm(
     F_landuse=0.,
     aviNOx_frac=0.,
     F_ref_aviNOx=0.0448,
+    E_ref_aviNOx=2.946,
+    F_ref_BC=0.04,
+    E_ref_BC=8.09,
     fossilCH4_frac=0.,
     natural=natural.Emissions.emissions,
     efficacy=np.array([1.]*9 + [3.] + [1.]*3),
@@ -219,6 +222,7 @@ def fair_scm(
     pi_tro3 =np.array([722, 170, 10, 4.29]),
     ghan_params = np.array([-1.95011431, 0.01107147, 0.01387492]),
     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.6
     tropO3_forcing='stevenson',
@@ -231,6 +235,7 @@ def fair_scm(
     contrail_forcing='NOx',
     kerosene_supply=0.,
     landuse_forcing='co2',
+    aCO2land=-0.00113789,
     ariaci_out=False,
     bcsnow_forcing='emissions',
     ):
@@ -524,7 +529,7 @@ def fair_scm(
         if type(emissions) is not bool:
             if contrail_forcing.lower()[0]=='n':   # from NOx emissions
                 F[:,7] = contrails.from_aviNOx(emissions, aviNOx_frac,
-                  F_ref=F_ref_aviNOx)
+                  F_ref=F_ref_aviNOx, E_ref=E_ref_aviNOx)
             elif contrail_forcing.lower()[0]=='f': # from kerosene production
                 F[:,7] = contrails.from_fuel(kerosene_supply)
             elif contrail_forcing.lower()[0]=='e': # external forcing timeseries
@@ -543,22 +548,28 @@ def fair_scm(
         if type(emissions) is not bool:
             if aerosol_forcing.lower()=='stevens':
                 ariaci[:,0], ariaci[:,1] = aerosols.Stevens(
-                  emissions, stevens_params=stevens_params, E_pi=E_pi[5])
+                  emissions, stevens_params=stevens_params, E_pi=E_pi[5], 
+                  ref_isSO2=ref_isSO2)
                 F[:,8] = np.sum(ariaci, axis=1)
             elif 'aerocom' in aerosol_forcing.lower():
-                ariaci[:,0] = aerosols.aerocom_direct(emissions, beta=b_aero)
+                ariaci[:,0] = aerosols.aerocom_direct(emissions, beta=b_aero,
+                  E_pi=E_pi)
                 if 'ghan' in aerosol_forcing.lower():
                     ariaci[:,1] = aerosols.ghan_indirect(emissions,
                       scale_AR5=scaleAerosolAR5,
                       fix_pre1850_RCP=fixPre1850RCP,
                       ghan_params=ghan_params)
+                elif 'stevens' in aerosol_forcing.lower():
+                    _, ariaci[:,1] = aerosols.Stevens(
+                      emissions, stevens_params=stevens_params, E_pi=E_pi[5],
+                      ref_isSO2=ref_isSO2)
                 F[:,8] = np.sum(ariaci, axis=1)
             elif aerosol_forcing.lower()[0] == 'e':
                 F[:,8] = F_aerosol
                 ariaci[:] = np.nan
             else:
                 raise ValueError("aerosol_forcing should be one of 'stevens', " +
-                  "aerocom, aerocom+ghan or external")
+                  "aerocom, aerocom+ghan, aerocom+stevens or external")
         else:
             F[:,8] = F_aerosol
             ariaci[:] = np.nan
@@ -569,7 +580,8 @@ def fair_scm(
         # concentrations
         if type(emissions) is not bool:
            if bcsnow_forcing.lower()[0]=='e':
-               F[:,9] = bc_snow.linear(emissions-E_pi)
+               F[:,9] = bc_snow.linear(emissions-E_pi, F_ref=F_ref_BC,
+                   E_ref=E_ref_BC)
            else:
                F[:,9] = F_bcsnow
         else:
@@ -582,7 +594,7 @@ def fair_scm(
         # concentrations
         if type(emissions) is not bool:
             if landuse_forcing.lower()[0]=='c':
-                F[:,10] = landuse.cumulative(emissions-E_pi)
+                F[:,10] = landuse.cumulative(emissions-E_pi, aCO2land=aCO2land)
             elif landuse_forcing.lower()[0]=='e':
                 F[:,10] = F_landuse
             else: