From 5fe624f4a485e42c6db51961ea15133220f763a8 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 11 Nov 2019 15:15:22 +0000 Subject: [PATCH 1/7] option added to get split of aerosol ERF into ari and aci --- fair/forcing/aerosols.py | 10 ++--- fair/forward.py | 53 +++++++++++++++---------- tests/reproduction/reproduction_test.py | 28 +++++++++++++ 3 files changed, 65 insertions(+), 26 deletions(-) diff --git a/fair/forcing/aerosols.py b/fair/forcing/aerosols.py index a1d17a04..e3378667 100644 --- a/fair/forcing/aerosols.py +++ b/fair/forcing/aerosols.py @@ -31,8 +31,7 @@ def Stevens(emissions, stevens_params=np.array([0.001875, 0.634, 60.]), ERFari = -alpha * em_SOx ERFaci = -beta * np.log(em_SOx/E_SOx_nat + 1) - F = ERFari + ERFaci - return F + return ERFari, ERFaci def aerocom_direct(emissions, @@ -68,9 +67,9 @@ def aerocom_direct(emissions, F_OC = beta[5] * em_OC F_NH3 = beta[6] * em_NH3 - F = F_SOx+F_CO+F_NMVOC+F_NOx+F_BC+F_OC+F_NH3 + ERFari = F_SOx+F_CO+F_NMVOC+F_NOx+F_BC+F_OC+F_NH3 - return F + return ERFari def ghan_indirect(emissions, fix_pre1850_RCP=True, scale_AR5=False, @@ -145,4 +144,5 @@ def _ERFaci(em, else: scale=1.0 - return (F_pd - F_1765) * scale + ERFaci = (F_pd - F_1765) * scale + return ERFaci diff --git a/fair/forward.py b/fair/forward.py index ce356774..a0d7865e 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -228,6 +228,7 @@ def fair_scm( contrail_forcing='NOx', kerosene_supply=0., landuse_forcing='co2', + ariaci_out=False ): if useStevenson is not None: @@ -295,26 +296,29 @@ def fair_scm( else: raise ValueError( "ghg_forcing should be 'etminan' (default) or 'myhre'") + # aerosol breakdown + ariaci = np.zeros((nt,2)) # Check natural emissions and convert to 2D array if necessary - if type(natural) in [float,int]: - natural = natural * np.ones((nt,2)) - elif type(natural) is np.ndarray: - if natural.ndim==1: - if natural.shape[0]!=2: - raise ValueError( - "natural emissions should be a 2-element or nt x 2 " + - "array") - natural = np.tile(natural, nt).reshape((nt,2)) - elif natural.ndim==2: - if natural.shape[1]!=2 or natural.shape[0]!=nt: - raise ValueError( - "natural emissions should be a 2-element or nt x 2 " + - "array") - else: - raise ValueError( - "natural emissions should be a scalar, 2-element, or nt x 2 " + - "array") + if emissions_driven: # don't check for conc runs + if type(natural) in [float,int]: + natural = natural * np.ones((nt,2)) + elif type(natural) is np.ndarray: + if natural.ndim==1: + if natural.shape[0]!=2: + raise ValueError( + "natural emissions should be a 2-element or nt x 2 " + + "array") + natural = np.tile(natural, nt).reshape((nt,2)) + elif natural.ndim==2: + if natural.shape[1]!=2 or natural.shape[0]!=nt: + raise ValueError( + "natural emissions should be a 2-element or nt x 2 " + + "array") + else: + raise ValueError( + "natural emissions should be a scalar, 2-element, or nt x 2 " + + "array") # check scale factor is correct shape. If 1D inflate to 2D if scale is None: @@ -512,16 +516,20 @@ def fair_scm( # Forcing from aerosols - again no feedback dependence if emissions_driven: if aerosol_forcing.lower()=='stevens': - F[:,8] = aerosols.Stevens(emissions, stevens_params=stevens_params) + ariaci[:,0], ariaci[:,1] = aerosols.Stevens( + emissions, stevens_params=stevens_params) + F[:,8] = np.sum(ariaci, axis=1) elif 'aerocom' in aerosol_forcing.lower(): - F[:,8] = aerosols.aerocom_direct(emissions, beta=b_aero) + ariaci[:,0] = aerosols.aerocom_direct(emissions, beta=b_aero) if 'ghan' in aerosol_forcing.lower(): - F[:,8] = F[:,8] + aerosols.ghan_indirect(emissions, + ariaci[:,1] = aerosols.ghan_indirect(emissions, scale_AR5=scaleAerosolAR5, fix_pre1850_RCP=fixPre1850RCP, ghan_params=ghan_params) + 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") @@ -729,5 +737,8 @@ def fair_scm( E_minus1 = emissions[-1] restart_out_val=(R_i[-1],T_j[-1],C_acc[-1],E_minus1) return C, F, T, restart_out_val + + if ariaci_out: + return C, F, T, ariaci else: return C, F, T diff --git a/tests/reproduction/reproduction_test.py b/tests/reproduction/reproduction_test.py index 6925f477..7f5b768b 100644 --- a/tests/reproduction/reproduction_test.py +++ b/tests/reproduction/reproduction_test.py @@ -349,3 +349,31 @@ def test_direct_o3tr_valueerror(): tropO3_forcing='external', F_tropO3 = np.zeros((4)) ) + + +def test_ariaci_aerocom(): + C,F,T,ariaci = fair.forward.fair_scm( + emissions=rcp45.Emissions.emissions, + ariaci_out=True, + aerosol_forcing='aerocom' + ) + assert np.allclose(F[:,8], np.sum(ariaci,axis=1)) + + +def test_ariaci_aerocom_ghan(): + C,F,T,ariaci = fair.forward.fair_scm( + emissions=rcp45.Emissions.emissions, + ariaci_out=True, + aerosol_forcing='aerocom+ghan' + ) + assert np.allclose(F[:,8], np.sum(ariaci,axis=1)) + + +def test_ariaci_stevens(): + C,F,T,ariaci = fair.forward.fair_scm( + emissions=rcp45.Emissions.emissions, + ariaci_out=True, + aerosol_forcing='stevens' + ) + assert np.allclose(F[:,8], np.sum(ariaci,axis=1)) + From 60cd43088797aa58e90058b3df9975affe98525f Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Sat, 23 Nov 2019 00:33:41 +0000 Subject: [PATCH 2/7] Functionality to run SLCFs emission-based and GHGs concentration-based --- fair/forward.py | 77 +++++++++++++++++++++++++++++++++++++------------ 1 file changed, 58 insertions(+), 19 deletions(-) diff --git a/fair/forward.py b/fair/forward.py index a0d7865e..604b6fbf 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -193,6 +193,7 @@ def fair_scm( iirf_max = carbon.iirf_max, iirf_h = carbon.iirf_h, C_pi=np.array([278., 722., 273., 34.497] + [0.]*25 + [13.0975, 547.996]), + E_pi=np.zeros(40), restart_in=False, restart_out=False, F_tropO3 = 0., @@ -228,12 +229,18 @@ def fair_scm( contrail_forcing='NOx', kerosene_supply=0., landuse_forcing='co2', - ariaci_out=False + ariaci_out=False, + bcsnow_forcing='emissions', ): + # Prevents later errors when SLCFs not specified + if type(emissions) is bool and not emissions_driven: + tropO3_forcing='external' + if useStevenson is not None: warnings.warn('"useStevenson" will be deprecated in v1.6; use '+ - '"stevenson", "regression" or "external"', DeprecationWarning) + 'tropO3_forcing keyword with "stevenson", "regression" or "external"', + DeprecationWarning) # is iirf_h < iirf_max? Don't stop the code, but warn user if iirf_h < iirf_max: @@ -478,18 +485,21 @@ def fair_scm( F[0,3] = np.sum((C[0,3:] - C_pi[3:]) * radeff.aslist[3:] * 0.001) # Tropospheric ozone: - if emissions_driven: + # v1.5 update: don't require emissions driven runs here for GHGs + # 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': - F[0,4] = ozone_tr.stevenson(emissions[0,:], C[0,1], + F[0,4] = ozone_tr.stevenson(emissions[0,:]-E_pi[:], C[0,1], T=np.sum(T_j[0,:]), feedback=useTropO3TFeedback, fix_pre1850_RCP=fixPre1850RCP) elif not useStevenson or tropO3_forcing[0].lower()=='r': - F[0,4] = ozone_tr.regress(emissions[0,:], beta=b_tro3) + F[0,4] = ozone_tr.regress(emissions[0,:]-E_pi[:], beta=b_tro3) else: F[0,4] = F_tropO3[0] else: - F[:,4] = F_tropO3 + F[0,4] = F_tropO3[0] # Stratospheric ozone depends on concentrations of ODSs (index 15-30) F[0,5] = ozone_st.magicc(C[0,15:], C_pi[15:]) @@ -499,9 +509,12 @@ def fair_scm( # Forcing from contrails. No climate feedback so can live outside # of forward model in this version - if emissions_driven: + # v1.5 update: don't require emissions driven runs here for GHGs + # because SLCFs can still be given as emissions with GHGs as + # concentrations + if type(emissions) is not bool: if contrail_forcing.lower()[0]=='n': # from NOx emissions - F[:,7] = contrails.from_aviNOx(emissions, aviNOx_frac) + F[:,7] = contrails.from_aviNOx(emissions-E_pi[0], aviNOx_frac) 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 @@ -514,15 +527,18 @@ def fair_scm( F[:,7] = F_contrails # Forcing from aerosols - again no feedback dependence - if emissions_driven: + # v1.5 update: don't require emissions driven runs here for GHGs + # because SLCFs can still be given as emissions with GHGs as + # concentrations + if type(emissions) is not bool: if aerosol_forcing.lower()=='stevens': ariaci[:,0], ariaci[:,1] = aerosols.Stevens( - emissions, stevens_params=stevens_params) + emissions-E_pi, stevens_params=stevens_params) 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-E_pi, beta=b_aero) if 'ghan' in aerosol_forcing.lower(): - ariaci[:,1] = aerosols.ghan_indirect(emissions, + ariaci[:,1] = aerosols.ghan_indirect(emissions-E_pi, scale_AR5=scaleAerosolAR5, fix_pre1850_RCP=fixPre1850RCP, ghan_params=ghan_params) @@ -535,26 +551,36 @@ def fair_scm( "aerocom, aerocom+ghan or external") else: F[:,8] = F_aerosol + ariaci[:] = np.nan # Black carbon on snow - no feedback dependence - if emissions_driven: - F[:,9] = bc_snow.linear(emissions) + # v1.5 update: don't require emissions driven runs here for GHGs + # because SLCFs can still be given as emissions with GHGs as + # concentrations + if type(emissions) is not bool: + if bcsnow_forcing.lower()[0]=='e': + F[:,9] = bc_snow.linear(emissions) + else: + F[:,9] = F_bcsnow else: F[:,9] = F_bcsnow # Land use change - either use a scaling with cumulative CO2 emissions # or an external time series - if emissions_driven: + # v1.5 update: don't require emissions driven runs here for GHGs + # because SLCFs can still be given as emissions with GHGs as + # concentrations + if type(emissions) is not bool: if landuse_forcing.lower()[0]=='c': F[:,10] = landuse.cumulative(emissions) elif landuse_forcing.lower()[0]=='e': F[:,10] = F_landuse else: raise ValueError( - "landuse_forcing should be one of 'co2' or 'external'") + "landuse_forcing should be one of 'co2' or 'external'") else: F[:,10] = F_landuse - + # Volcanic and solar copied straight to the output arrays F[:,11] = F_volcanic F[:,12] = F_solar @@ -644,13 +670,13 @@ def fair_scm( F[t,3] = np.sum((C[t,3:] - C_pi[3:]) * radeff.aslist[3:] * 0.001) if useStevenson and tropO3_forcing[0].lower()=='s': - F[t,4] = ozone_tr.stevenson(emissions[t,:], + F[t,4] = ozone_tr.stevenson(emissions[t,:]-E_pi, C[t,1], T=T[t-1], feedback=useTropO3TFeedback, fix_pre1850_RCP=fixPre1850RCP) elif not useStevenson or tropO3_forcing[0].lower()=='r': - F[t,4] = ozone_tr.regress(emissions[t,:], beta=b_tro3) + F[t,4] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3) else: F[t,4] = F_tropO3[t] F[t,5] = ozone_st.magicc(C[t,15:], C_pi[15:]) @@ -702,6 +728,19 @@ def fair_scm( 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 type(emissions) is not bool: + if useStevenson and tropO3_forcing[0].lower()=='s': + F[t,4] = ozone_tr.stevenson(emissions[t,:]-E_pi, + C[t,1], + T=T[t-1], + feedback=useTropO3TFeedback, + fix_pre1850_RCP=fixPre1850RCP) + elif not useStevenson or tropO3_forcing[0].lower()=='r': + F[t,4] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3) + else: + F[t,4] = 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) From d4ca2de528860702d8b3eaf3f7b0ea058dcc8b73 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Sat, 23 Nov 2019 10:11:18 +0000 Subject: [PATCH 3/7] More flexibility for Stevenson-style tropospheric ozone relationships tuned on CMIP6 --- fair/forcing/ozone_tr.py | 71 +++++++++++++++++++++++++++++++++++++--- fair/forward.py | 35 ++++++++++++++++---- tests/unit/unit_test.py | 17 ++++++++++ 3 files changed, 111 insertions(+), 12 deletions(-) diff --git a/fair/forcing/ozone_tr.py b/fair/forcing/ozone_tr.py index 77a03e3e..7c5785c5 100644 --- a/fair/forcing/ozone_tr.py +++ b/fair/forcing/ozone_tr.py @@ -33,7 +33,66 @@ def regress(emissions, return F -def stevenson(emissions, C_CH4, T=0, feedback=False, fix_pre1850_RCP=False): +def cmip6_stevenson(emissions, C_CH4, T=0, feedback=False, + PI=np.array([722, 170, 10, 4.29]), + beta=np.array([1.77871043e-04, 5.80173377e-05, 2.09151270e-03, + 1.94458719e-04])): + """Calculates tropospheric ozone forcing from precursor emissions based on + Stevenson et al, 2013 10.5194/acp-13-3063-2013 + + Inputs: + emissions: (nt x 40) numpy array + C_CH4 : (nt) numpy array of methane concentrations, ppb + + Keywords: + T : change in surface temperature since pre-industrial + feedback : True or False - include temperature feedback on ozone + forcing? + PI: : 4-element array of pre-industrial CH4 concentrations, + CO emissions, NMVOC emissions and NOx emissions + beta: : coefficients of how CH4 concentrations, CO emissions, + NMVOC emissions and NOx emissions affect forcing + + Outputs: + tropospheric ozone ERF time series. + """ + + # expand to 2D/1D if not already + if emissions.ndim == 1: + nspec = len(emissions) + emissions = emissions.reshape((1, nspec)) + if np.isscalar(C_CH4): + C_CH4 = np.ones(1)*C_CH4 + + year, em_CO, em_NMVOC, em_NOx = emissions[:,[0, 6, 7, 8]].T + nt = len(year) + F_CH4, F_CO, F_NMVOC, F_NOx = np.zeros((4,nt)) + + for i in range(nt): + F_CH4[i] = beta[0] * (C_CH4[i]-PI[0]) + F_CO[i] = beta[1] * (em_CO[i]-PI[1]) + F_NMVOC[i] = beta[2] * (em_NMVOC[i]-PI[2]) + F_NOx[i] = beta[3] * (em_NOx[i]-PI[3]) + + # Include the effect of climate feedback? We fit a curve to the 2000, 2030 + # and 2100 best estimates of feedback based on middle-of-the-road + # temperature projections. + def temperature_feedback(T, a=0.03189267, b=1.34966941, c=-0.03214807): + if T<=0: + return 0 + else: + return a*np.exp(-b*T)+c + + if feedback: + F = F_CH4 + F_CO + F_NMVOC + F_NOx + temperature_feedback(T) + else: + F = F_CH4 + F_CO + F_NMVOC + F_NOx + + return F + + +def stevenson(emissions, C_CH4, T=0, feedback=False, fix_pre1850_RCP=False, + PI=np.array([722, 170, 10, 4.29])): """Calculates tropospheric ozone forcing from precursor emissions based on Stevenson et al, 2013 10.5194/acp-13-3063-2013 @@ -48,6 +107,8 @@ def stevenson(emissions, C_CH4, T=0, feedback=False, fix_pre1850_RCP=False): fix_pre1850_RCP: Use different relationship for 1750/65 to 1850 based on anthropogenic emissions from Skeie et al (2011) for 1750 (atmos-chem-phys.net/11/11827/2011) + PI: : 4-element array of pre-industrial CH4 concentrations, + CO emissions, NMVOC emissions and NOx emissions Outputs: tropospheric ozone ERF time series. @@ -69,11 +130,11 @@ def stevenson(emissions, C_CH4, T=0, feedback=False, fix_pre1850_RCP=False): for i in range(nt): if year[i]>=1850 or fix_pre1850_RCP==False: - F_CH4[i] = 0.166/960 * (C_CH4[i]-722) - F_CO[i] = 0.058/681.8 * (em_CO[i]-170) - F_NMVOC[i] = 0.035/155.84 * (em_NMVOC[i]-10) + F_CH4[i] = 0.166/960 * (C_CH4[i]-PI[0]) + F_CO[i] = 0.058/681.8 * (em_CO[i]-PI[1]) + F_NMVOC[i] = 0.035/155.84 * (em_NMVOC[i]-PI[2]) F_NOx[i] = 0.119/61.16 * (em_NOx[i] * - molwt.NO / molwt.N - 4.29) + molwt.NO / molwt.N - PI[3]) # The RCP scenarios give a negative forcing prior to ~1780. This is # because the anthropogenic emissions are given to be zero in RCPs but # not zero in the Skeie numbers which are used here. This can be fixed diff --git a/fair/forward.py b/fair/forward.py index 604b6fbf..4c09862c 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -215,6 +215,7 @@ def fair_scm( b_aero = np.array([-6.2227e-3, 0.0, -3.8392e-4, -1.16551e-3, 1.601537e-2, -1.45339e-3, -1.55605e-3]), b_tro3 = np.array([2.8249e-4, 1.0695e-4, -9.3604e-4, 99.7831e-4]), + 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.]), useMultigas=True, @@ -490,10 +491,17 @@ def fair_scm( # concentrations if type(emissions) is not bool: if useStevenson and tropO3_forcing[0].lower()=='s': - F[0,4] = ozone_tr.stevenson(emissions[0,:]-E_pi[:], C[0,1], + F[0,4] = ozone_tr.stevenson(emissions[0,:], C[0,1], T=np.sum(T_j[0,:]), feedback=useTropO3TFeedback, - fix_pre1850_RCP=fixPre1850RCP) + fix_pre1850_RCP=fixPre1850RCP, + PI=pi_tro3) + elif tropO3_forcing[0].lower()=='c': + F[0,4] = ozone_tr.cmip6_stevenson(emissions[0,:], C[0,1], + T=np.sum(T_j[0,:]), + 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[0,4] = ozone_tr.regress(emissions[0,:]-E_pi[:], beta=b_tro3) else: @@ -533,12 +541,12 @@ def fair_scm( if type(emissions) is not bool: if aerosol_forcing.lower()=='stevens': ariaci[:,0], ariaci[:,1] = aerosols.Stevens( - emissions-E_pi, stevens_params=stevens_params) + emissions, stevens_params=stevens_params) F[:,8] = np.sum(ariaci, axis=1) elif 'aerocom' in aerosol_forcing.lower(): - ariaci[:,0] = aerosols.aerocom_direct(emissions-E_pi, beta=b_aero) + ariaci[:,0] = aerosols.aerocom_direct(emissions, beta=b_aero) if 'ghan' in aerosol_forcing.lower(): - ariaci[:,1] = aerosols.ghan_indirect(emissions-E_pi, + ariaci[:,1] = aerosols.ghan_indirect(emissions, scale_AR5=scaleAerosolAR5, fix_pre1850_RCP=fixPre1850RCP, ghan_params=ghan_params) @@ -670,11 +678,18 @@ def fair_scm( F[t,3] = np.sum((C[t,3:] - C_pi[3:]) * radeff.aslist[3:] * 0.001) if useStevenson and tropO3_forcing[0].lower()=='s': - F[t,4] = ozone_tr.stevenson(emissions[t,:]-E_pi, + F[t,4] = ozone_tr.stevenson(emissions[t,:], C[t,1], T=T[t-1], feedback=useTropO3TFeedback, - fix_pre1850_RCP=fixPre1850RCP) + fix_pre1850_RCP=fixPre1850RCP, + PI=pi_tro3) + elif tropO3_forcing[0].lower()=='c': + F[t,4] = 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) else: @@ -735,6 +750,12 @@ def fair_scm( 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], + 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) else: diff --git a/tests/unit/unit_test.py b/tests/unit/unit_test.py index 1f3309ab..a52d59d6 100755 --- a/tests/unit/unit_test.py +++ b/tests/unit/unit_test.py @@ -398,3 +398,20 @@ def test_deprecated_ozone_stevenson(): useStevenson=True ) +def test_cmip6_stevenson(): + C1,F1,T1 = fair.forward.fair_scm( + emissions=rcp85.Emissions.emissions, + tropO3_forcing='cmip6', + E_pi=rcp85.Emissions.emissions[0,:], + C_pi=rcp85.Concentrations.gases[0,:] + ) + + C2,F2,T2 = fair.forward.fair_scm( + emissions=rcp85.Emissions.emissions, + tropO3_forcing='stevenson', + E_pi=rcp85.Emissions.emissions[0,:], + C_pi=rcp85.Concentrations.gases[0,:] + ) + + # check differences + assert np.any(F2[:,4]!=F1[:,4]) From 467b5bd82d5aa2b0170b64a82c0ed50f7ee6c9db Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Sat, 23 Nov 2019 11:33:27 +0000 Subject: [PATCH 4/7] Update Stevens aerosol function to specify a reference emissions level to take forcing relative to --- fair/forcing/aerosols.py | 13 +++++++++---- fair/forward.py | 4 ++-- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/fair/forcing/aerosols.py b/fair/forcing/aerosols.py index e3378667..2a7d2a8c 100644 --- a/fair/forcing/aerosols.py +++ b/fair/forcing/aerosols.py @@ -6,7 +6,7 @@ def Stevens(emissions, stevens_params=np.array([0.001875, 0.634, 60.]), - ref_isSO2=True): + ref_isSO2=True, E_pi=0): """Calculates aerosol forcing based on Stevens (2015) that relates sulphate aerosol forcing to SOx emissions in a logarithmic fashion. @@ -18,8 +18,10 @@ def Stevens(emissions, stevens_params=np.array([0.001875, 0.634, 60.]), 1. scaling parameter for ERFaci (beta) 2. natural emissions of SOx in Mt/yr ref_isSO2: True if E_SOx_nat is in units of SO2 rather than S. + E_pi: pre-industrial/reference emissions of SO2 (or S). Output: - F: aerosol effective radiative forcing + ERFari, ERFaci: aerosol effective radiative forcing due to + aerosol-radiation interactions and aerosol-cloud interactions. """ alpha, beta, E_SOx_nat = stevens_params @@ -28,9 +30,12 @@ def Stevens(emissions, stevens_params=np.array([0.001875, 0.634, 60.]), if ref_isSO2: factor = molwt.SO2/molwt.S em_SOx = emissions[:,5] * factor + em_pi = E_pi * factor - ERFari = -alpha * em_SOx - ERFaci = -beta * np.log(em_SOx/E_SOx_nat + 1) + 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)) ) return ERFari, ERFaci diff --git a/fair/forward.py b/fair/forward.py index 4c09862c..a1fe47fc 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -240,7 +240,7 @@ def fair_scm( if useStevenson is not None: warnings.warn('"useStevenson" will be deprecated in v1.6; use '+ - 'tropO3_forcing keyword with "stevenson", "regression" or "external"', + 'tropO3_forcing keyword with "cmip6", "stevenson", "regression" or "external"', DeprecationWarning) # is iirf_h < iirf_max? Don't stop the code, but warn user @@ -541,7 +541,7 @@ 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) + emissions, stevens_params=stevens_params, E_pi=E_pi[5]) F[:,8] = np.sum(ariaci, axis=1) elif 'aerocom' in aerosol_forcing.lower(): ariaci[:,0] = aerosols.aerocom_direct(emissions, beta=b_aero) From 740bd8315d4952e059a31ddc8ef1016b1cf8d5b1 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Tue, 26 Nov 2019 11:28:27 +0000 Subject: [PATCH 5/7] reference emissions to zero for BC_snow and land use --- fair/forward.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fair/forward.py b/fair/forward.py index a1fe47fc..9b756251 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -567,7 +567,7 @@ def fair_scm( # concentrations if type(emissions) is not bool: if bcsnow_forcing.lower()[0]=='e': - F[:,9] = bc_snow.linear(emissions) + F[:,9] = bc_snow.linear(emissions-E_pi) else: F[:,9] = F_bcsnow else: @@ -580,7 +580,7 @@ def fair_scm( # concentrations if type(emissions) is not bool: if landuse_forcing.lower()[0]=='c': - F[:,10] = landuse.cumulative(emissions) + F[:,10] = landuse.cumulative(emissions-E_pi) elif landuse_forcing.lower()[0]=='e': F[:,10] = F_landuse else: From 970fc1d652692170664fa2e1a4c238f1a11eb9f5 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Fri, 29 Nov 2019 09:39:29 +0000 Subject: [PATCH 6/7] Don't subtract pre-industrial emissions for NOx calculation --- fair/forward.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fair/forward.py b/fair/forward.py index 9b756251..8147f63b 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -522,7 +522,7 @@ def fair_scm( # concentrations if type(emissions) is not bool: if contrail_forcing.lower()[0]=='n': # from NOx emissions - F[:,7] = contrails.from_aviNOx(emissions-E_pi[0], aviNOx_frac) + F[:,7] = contrails.from_aviNOx(emissions, aviNOx_frac) 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 From 936fdedddcecf190958c89da774397b01dbe20e4 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Fri, 29 Nov 2019 10:31:05 +0000 Subject: [PATCH 7/7] add API keyword for scaling contrail forcing to reference year --- fair/forward.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/fair/forward.py b/fair/forward.py index 8147f63b..6ae23c89 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -204,6 +204,7 @@ def fair_scm( F_bcsnow=0., F_landuse=0., aviNOx_frac=0., + F_ref_aviNOx=0.0448, fossilCH4_frac=0., natural=natural.Emissions.emissions, efficacy=np.array([1.]*9 + [3.] + [1.]*3), @@ -522,7 +523,8 @@ def fair_scm( # concentrations if type(emissions) is not bool: if contrail_forcing.lower()[0]=='n': # from NOx emissions - F[:,7] = contrails.from_aviNOx(emissions, aviNOx_frac) + F[:,7] = contrails.from_aviNOx(emissions, aviNOx_frac, + F_ref=F_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