From 8a00c84d870f4e4c9cdb72b3aa797d3821986b8e Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 17:30:59 +0000 Subject: [PATCH 1/6] first outline of new ozone treatment plus failing test --- fair/forcing/ozone.py | 91 +++++++++++++++++++++++++++++++++++++++++ tests/unit/unit_test.py | 4 ++ 2 files changed, 95 insertions(+) create mode 100644 fair/forcing/ozone.py diff --git a/fair/forcing/ozone.py b/fair/forcing/ozone.py new file mode 100644 index 00000000..6d58bbb5 --- /dev/null +++ b/fair/forcing/ozone.py @@ -0,0 +1,91 @@ +from __future__ import division + +import numpy as np +from ..constants import cl_atoms, br_atoms, fracrel + + +def thornhill_skeie( + emissions, + concentrations, + temperature=0, + feedback=-0.037, + beta=np.array([2.33379720e-04, 1.27179106e-03, -6.69347820e-05, + 1.14647701e-04, 5.14366051e-12, 3.78354423e-03]), + emissions_pi=np.zeros(40), + concentrations_pi=np.zeros(31), + ) + """Calculates total ozone forcing from precursor emissions and + concentrations based on AerChemMIP and CMIP6 Historical behaviour + + Skeie et al. (2020) + Thornhill et al. (2021) + + Unlike Stevenson CMIP5 no distinction is made for tropospheric and + stratospheric. + + With this formulation, ozone forcing depends on concentrations of + CH4, N2O, ozone-depleting halogens, and emissions of CO, NVMOC and NOx, + but any combination of emissions and concentrations are allowed. + + Inputs: + emissions: (nt x 40) numpy array in FaIR default units + concentrations: (nt x 31) numpy array of GHGs in FaIR default units + temperature: global mean surface temperature (for feedback) + feedback: temperature feedback on ozone forcing (W/m2/K) - set to zero + or False to turn off + beta: 6-element array of radiative efficiency coefficients in order of + CH4 concentrations, W m-2 ppb-1 + N2O concentrations, W m-2 ppb-1 + ODS concentrations in EESC, W m-2 ppt-1 + CO emissions, W m-2 (Mt yr-1)-1 + NMVOC emissions, W m-2 (Mt yr-1)-1 + NOx emissions, W m-2 (MtN yr-1)-1 + + emissions_pi: pre-industrial/reference emissions + concentrations_pi: pre-industrial/reference concentrations + + Outputs: + ozone ERF time series. + """ + + # we allow 2D output for quick calculation if feedbacks turned off + if emissions.ndim == 1: + nspec = len(emissions) + emissions = emissions.reshape((1, nspec)) + if concentrations.ndim == 1: + nspec = len(concentrations) + concentrations = concentrations.reshape((1, nspec)) + + nt = emissions.shape[0] + + # calculate EESC for halogens + cl = np.array(cl_atoms.aslist) + br = np.array(br_atoms.aslist) + fc = np.array(fracrel.aslist) + + def eesc(c_ods, c_ods_pi): + return ( + np.sum(cl * (c_ods-c_cods_pi) * fc/fc[0]) + + 45 * np.sum(br * (c_ods-c_ods_pi) * fc/fc[0]) + ) * fc[0] + + + c_ch4, c_n2o = concentrations[:, [2, 3]].T + delta_c_ods = eesc(concentrations[:,15:].T, concentrations_pi[15:]) + e_co, e_nmvoc, e_nox = emissions[:,[6, 7, 8]].T + c_ch4_pi, c_n2o_pi = concentrations_pi[2, 3] + e_co_pi, e_nmvoc_pi, e_nox_pi = emissions_pi[6, 7, 8] + + forcing = np.zeros(nt) + + for i in range(nt): + f_ch4 = beta[0] * (c_ch4[i] - c_ch4_pi) + f_n2o = beta[1] * (c_n2o[i] - c_n2o_pi) + f_ods = beta[2] * (delta_c_ods[i]) + f_co = beta[3] * (e_co[i] - e_co_pi) + f_nmvoc = beta[4] * (e_nmvoc[i] - e_nmvoc_pi) + f_nox = beta[5] * (e_nox[i] - e_nox_pi) + ) + forcing[i] = f_ch4 + f_n2o + f_ods + f_co + f_nmvoc + f_nox + feedback * temperature[i] + + return forcing diff --git a/tests/unit/unit_test.py b/tests/unit/unit_test.py index cd23d8c6..366584a9 100755 --- a/tests/unit/unit_test.py +++ b/tests/unit/unit_test.py @@ -523,3 +523,7 @@ def test_inverse_millar_fin(): temperature_function = 'Millar', F_in = rcp85.Forcing.total ) + + +def test_thornhill_skeie(): + raise AssertionError From b8fd4c617c1487cf994f9a8a4df113071d4f9094 Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 17:51:41 +0000 Subject: [PATCH 2/6] update changelog --- CHANGELOG.rst | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 2de0ebbd..d248a5c3 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,12 +4,20 @@ Changelog master ------ +(`#99 `_) Add option to use ozone forcing dependent on N2O concentrations + +v1.6.1 +------ + (`#87 `_) Add Geoffroy temperature and prescribed forcing to inverse (`#84 `_) Fix concentration-driven runs for 45-species FaIR (`#83 `_) Support scmdata >= 0.7.1 +earlier +------- + (`#82 `_) Expand AR6 forcing diagnostics to get aerosol direct forcing by species (`#78 `_) Optimise ``fair.tools.scmdf.scmdf_to_emissions`` a bit From db25b058d64fe67ec46f8fff106e384efe8bed6b Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 19:43:26 +0000 Subject: [PATCH 3/6] get a passing unit test --- fair/forcing/ozone.py | 15 +++++------ fair/forward.py | 56 +++++++++++++++++++++++++++++++++-------- tests/unit/unit_test.py | 10 +++++++- 3 files changed, 62 insertions(+), 19 deletions(-) diff --git a/fair/forcing/ozone.py b/fair/forcing/ozone.py index 6d58bbb5..d988f0b0 100644 --- a/fair/forcing/ozone.py +++ b/fair/forcing/ozone.py @@ -13,7 +13,7 @@ def thornhill_skeie( 1.14647701e-04, 5.14366051e-12, 3.78354423e-03]), emissions_pi=np.zeros(40), concentrations_pi=np.zeros(31), - ) + ): """Calculates total ozone forcing from precursor emissions and concentrations based on AerChemMIP and CMIP6 Historical behaviour @@ -65,27 +65,28 @@ def thornhill_skeie( def eesc(c_ods, c_ods_pi): return ( - np.sum(cl * (c_ods-c_cods_pi) * fc/fc[0]) + + np.sum(cl * (c_ods-c_ods_pi) * fc/fc[0]) + 45 * np.sum(br * (c_ods-c_ods_pi) * fc/fc[0]) ) * fc[0] c_ch4, c_n2o = concentrations[:, [2, 3]].T - delta_c_ods = eesc(concentrations[:,15:].T, concentrations_pi[15:]) +# delta_c_ods = eesc(concentrations[:,15:].T, concentrations_pi[None, 15:]) + c_ods = concentrations[:,15:] e_co, e_nmvoc, e_nox = emissions[:,[6, 7, 8]].T - c_ch4_pi, c_n2o_pi = concentrations_pi[2, 3] - e_co_pi, e_nmvoc_pi, e_nox_pi = emissions_pi[6, 7, 8] + c_ch4_pi, c_n2o_pi = concentrations_pi[[2, 3]] + c_ods_pi = concentrations_pi[15:] + e_co_pi, e_nmvoc_pi, e_nox_pi = emissions_pi[[6, 7, 8]] forcing = np.zeros(nt) for i in range(nt): f_ch4 = beta[0] * (c_ch4[i] - c_ch4_pi) f_n2o = beta[1] * (c_n2o[i] - c_n2o_pi) - f_ods = beta[2] * (delta_c_ods[i]) + f_ods = beta[2] * eesc(c_ods[i], c_ods_pi) f_co = beta[3] * (e_co[i] - e_co_pi) f_nmvoc = beta[4] * (e_nmvoc[i] - e_nmvoc_pi) f_nox = beta[5] * (e_nox[i] - e_nox_pi) - ) forcing[i] = f_ch4 + f_n2o + f_ods + f_co + f_nmvoc + f_nox + feedback * temperature[i] return forcing diff --git a/fair/forward.py b/fair/forward.py index 331c4bc5..945105c4 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -8,8 +8,8 @@ from .constants import molwt, lifetime, radeff from .constants.general import M_ATMOS, ppm_gtc from .defaults import carbon, thermal -from .forcing import ozone_tr, ozone_st, h2o_st, contrails, aerosols, bc_snow,\ - landuse +from .forcing import (ozone, ozone_tr, ozone_st, h2o_st, contrails, aerosols, + bc_snow, landuse) from .gas_cycle.gir import calculate_alpha, step_concentration from .gas_cycle.fair1 import carbon_cycle from .forcing.ghg import co2_log, minor_gases @@ -90,6 +90,7 @@ def fair_scm( ref_isSO2=True, # is Stevens SO2 emissions in units SO2 (T) or S (F) useMultigas=True, tropO3_forcing='stevenson', + ozone_feedback=-0.037, # W m-2 K-1, only for Thornhill-Skeie lifetimes=False, aerosol_forcing="aerocom+ghan", scaleAerosolAR5=True, @@ -441,27 +442,38 @@ def fair_scm( # because SLCFs can still be given as emissions with GHGs as # concentrations if type(emissions) is not bool: - if tropO3_forcing[0].lower()=='s': + if tropO3_forcing[0].lower()=='s': # stevenson F[0,iF_tro3] = ozone_tr.stevenson(emissions[0,:], C[0,1], T=np.sum(T_j[0,:]), feedback=useTropO3TFeedback, fix_pre1850_RCP=fixPre1850RCP, PI=pi_tro3) - elif tropO3_forcing[0].lower()=='c': + elif tropO3_forcing[0].lower()=='c': # cmip6 stevenson F[0,iF_tro3] = 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 tropO3_forcing[0].lower()=='r': + elif tropO3_forcing[0].lower()=='r': # regression F[0,iF_tro3] = ozone_tr.regress(emissions[0,:]-E_pi[:], beta=b_tro3) + elif tropO3_forcing[0].lower()=='t': # thornhill-skeie + F[0,iF_tro3] = ozone.thornhill_skeie( + emissions=emissions[0,:], + concentrations=C[0,:], + temperature=T[0], + ozone_feedback=ozone_feedback, + beta=b_tro3, + emissions_pi=E_pi, + concentrations_pi=C_pi, + ) else: F[0,iF_tro3] = F_tropO3[0] else: F[0,iF_tro3] = F_tropO3[0] - # Stratospheric ozone depends on concentrations of ODSs (index 15-30) - F[0,iF_sto3] = ozone_st.magicc(C[0,15:], C_pi[15:]) + if tropO3_forcing[0].lower()!='t': + # Stratospheric ozone depends on concentrations of ODSs (index 15-30) + F[0,iF_sto3] = ozone_st.magicc(C[0,15:], C_pi[15:]) # Stratospheric water vapour is a function of the methane ERF F[0,iF_ch4h] = h2o_st.linear(F[0,1], ratio=stwv_from_ch4) @@ -711,9 +723,21 @@ def fair_scm( beta=b_tro3) elif tropO3_forcing[0].lower()=='r': F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3) + elif tropO3_forcing[0].lower()=='t': + F[t,iF_tro3] = ozone.thornhill_skeie( + emissions=emissions[t,:], + concentrations=C[t,:], + temperature=T[t-1], + ozone_feedback=ozone_feedback, + beta=b_tro3, + emissions_pi=E_pi, + concentrations_pi=C_pi, + ) else: F[t,iF_tro3] = F_tropO3[t] - F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:]) + + if tropO3_forcing[0].lower()!='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 @@ -824,11 +848,21 @@ def fair_scm( beta=b_tro3) elif tropO3_forcing[0].lower()=='r': F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3) + elif tropO3_forcing[0].lower()=='t': + F[t,iF_tro3] = ozone.thornhill_skeie( + emissions=emissions[t,:], + concentrations=C[t,:], + temperature=T[t-1], + ozone_feedback=ozone_feedback, + beta=b_tro3, + emissions_pi=E_pi, + concentrations_pi=C_pi, + ) else: F[t,iF_tro3] = F_tropO3[t] - else: - F[t,iF_tro3] = F_tropO3[t] - F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:]) + + if tropO3_forcing[0].lower()!='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 diff --git a/tests/unit/unit_test.py b/tests/unit/unit_test.py index 366584a9..02e4def0 100755 --- a/tests/unit/unit_test.py +++ b/tests/unit/unit_test.py @@ -8,6 +8,7 @@ from fair.defaults import carbon from fair.ancil import cmip5_annex2_forcing from fair.forcing.ozone_tr import regress +from fair.forcing.ozone import thornhill_skeie from fair.inverse import inverse_fair_scm import numpy as np import os @@ -526,4 +527,11 @@ def test_inverse_millar_fin(): def test_thornhill_skeie(): - raise AssertionError + forcing = thornhill_skeie( + emissions = rcp85.Emissions.emissions, + concentrations = rcp85.Concentrations.gases, + temperature = np.linspace(0,5,736), + feedback = -0.037, + emissions_pi = rcp85.Emissions.emissions[0,:], + concentrations_pi = rcp85.Concentrations.gases[0,:], + ) From 1c11be8371f61be6d79e61ab6ccd6a61150eb9fa Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 19:47:35 +0000 Subject: [PATCH 4/6] add back accidentally deleted line --- fair/forward.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/fair/forward.py b/fair/forward.py index 945105c4..63e0f296 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -860,6 +860,8 @@ def fair_scm( ) else: F[t,iF_tro3] = F_tropO3[t] + else: + F[t,iF_tro3] = F_tropO3[t] if tropO3_forcing[0].lower()!='t': F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:]) From 7575a0bc3ec61b4e2fdf67279326ef4ee00b93ac Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 21:02:55 +0000 Subject: [PATCH 5/6] add notebook and run diagnostic tests --- fair/forcing/ozone.py | 7 +- fair/forcing/ozone_tr.py | 1 - fair/forward.py | 6 +- notebooks/ozone-comparison.ipynb | 136 +++++++++++++++++++++++++++++++ 4 files changed, 144 insertions(+), 6 deletions(-) create mode 100644 notebooks/ozone-comparison.ipynb diff --git a/fair/forcing/ozone.py b/fair/forcing/ozone.py index d988f0b0..e1860bb8 100644 --- a/fair/forcing/ozone.py +++ b/fair/forcing/ozone.py @@ -70,15 +70,18 @@ def eesc(c_ods, c_ods_pi): ) * fc[0] - c_ch4, c_n2o = concentrations[:, [2, 3]].T + c_ch4, c_n2o = concentrations[:, [1, 2]].T # delta_c_ods = eesc(concentrations[:,15:].T, concentrations_pi[None, 15:]) c_ods = concentrations[:,15:] e_co, e_nmvoc, e_nox = emissions[:,[6, 7, 8]].T - c_ch4_pi, c_n2o_pi = concentrations_pi[[2, 3]] + c_ch4_pi, c_n2o_pi = concentrations_pi[[1, 2]] c_ods_pi = concentrations_pi[15:] e_co_pi, e_nmvoc_pi, e_nox_pi = emissions_pi[[6, 7, 8]] + forcing = np.zeros(nt) + if np.isscalar(temperature): + temperature = np.ones(nt) * temperature for i in range(nt): f_ch4 = beta[0] * (c_ch4[i] - c_ch4_pi) diff --git a/fair/forcing/ozone_tr.py b/fair/forcing/ozone_tr.py index 7c5785c5..10c5c9cc 100644 --- a/fair/forcing/ozone_tr.py +++ b/fair/forcing/ozone_tr.py @@ -87,7 +87,6 @@ def temperature_feedback(T, a=0.03189267, b=1.34966941, c=-0.03214807): F = F_CH4 + F_CO + F_NMVOC + F_NOx + temperature_feedback(T) else: F = F_CH4 + F_CO + F_NMVOC + F_NOx - return F diff --git a/fair/forward.py b/fair/forward.py index 63e0f296..e04f0f7d 100755 --- a/fair/forward.py +++ b/fair/forward.py @@ -461,7 +461,7 @@ def fair_scm( emissions=emissions[0,:], concentrations=C[0,:], temperature=T[0], - ozone_feedback=ozone_feedback, + feedback=ozone_feedback, beta=b_tro3, emissions_pi=E_pi, concentrations_pi=C_pi, @@ -728,7 +728,7 @@ def fair_scm( emissions=emissions[t,:], concentrations=C[t,:], temperature=T[t-1], - ozone_feedback=ozone_feedback, + feedback=ozone_feedback, beta=b_tro3, emissions_pi=E_pi, concentrations_pi=C_pi, @@ -853,7 +853,7 @@ def fair_scm( emissions=emissions[t,:], concentrations=C[t,:], temperature=T[t-1], - ozone_feedback=ozone_feedback, + feedback=ozone_feedback, beta=b_tro3, emissions_pi=E_pi, concentrations_pi=C_pi, diff --git a/notebooks/ozone-comparison.ipynb b/notebooks/ozone-comparison.ipynb new file mode 100644 index 00000000..ed756794 --- /dev/null +++ b/notebooks/ozone-comparison.ipynb @@ -0,0 +1,136 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Comparison of different ozone forcing treatments" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from fair.forward import fair_scm\n", + "from fair.RCPs import rcp85\n", + "import matplotlib.pyplot as pl\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEICAYAAACwDehOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABWFklEQVR4nO3dd3gU1frA8e/Zkt4JvfdeJAFUREBFEBEEReygXgte+xUvNuz3Z79YEC8WiqggRQVEKQIiAkLoEHoPNaT3bDm/P2azbCoJpOf9PM8+O+XMzLubzXln5sycUVprhBBCCABTRQcghBCi8pCkIIQQwk2SghBCCDdJCkIIIdwkKQghhHCTpCCEEMJNkoIoF0oprZRqVdFxXAyl1HCl1HGlVKpS6rJSXO9dSqmlpbW+POvepZTqVxbrFtWbJIUazlXR5bycSqkMj/G7Clmmn1IqprxjrUDvA49prQO01ltKa6Va62+11teX1vryrLuj1nrVxSzrSuBprt/ACaXUh0ops8f8+kqpr5RSp5RSKUqpPUqp15RS/hdaXinVTCm1WCmVoJQ6rZT6VCllKSSOfq7fpOdvdPTFfCZRfJIUajhXRRegtQ4AjgE3eUz7tqLjqySaArsuZkHPyrSK6er6TfQFRgH3AyilwoB1gC9whdY6EBgAhAAtC1j+WuBO4EHX9M+As0B9oJtr/Y8WEcdJz9+o1np66Xw8URhJCqJASilvpdREpdRJ12uia5o/8CvQwGPvrYFSqqdSap1SKtG1B/mpUsqrmNtqoJRaoJSKV0odUEo96DEv0WM7aa690GaueQ+6yse7lm/gsZxWSj2ilNrvWsckpZTymH+/Umq3a491iVKqaSHfQSpgBrYppQ66prdXSq1yrXeXUmqoxzLTlFKTXXvDaUB/pVRjpdR8pVSsUipOKfWpq+wYpdSa4sSslDIrpT5QSp1TSh1WSj3mKl/YXvYRpdR1ruFXlVI/KKVmuPbsdymlIovzt9FaHwD+wqjAAZ4BUoC7tdZHXGWOa62f1FpvL2D5PcCfQCfXpObAD1rrTK31aeA3oGNxYhHlQ5KCKMyLwOUYlUFXoCfwktY6DbiB3HtwJwEH8DQQDlyBsYdY1B6gp1lADNAAuBX4j1LqGgCtdYjHkcxHGBXMCdf8/wNuw9jrPOpaj6chQA+gi6vcQACl1DDgBWAEUNu1zu/zBqW1znJtF4w935ZKKSuwEFgK1AEeB75VSrX1WPRO4C0gEGOvepErvmZAwwLivGDMGHvaN2D8PboDNxexjoIMdW03BFgAfFqchZRS7YA+wAHXpOuA+VprZzGX7+BaPue020TgdqWUn1KqIcZn+q2IVdRRSp1xJcL/5pyiEmVIay0veaG1BjgCXOcaPggM9pg3EDjiGu4HxFxgXU8BP3qMa6BVAeUaYySUQI9p/wdMy1NulCu+2q7xr4B3PeYHADagmcf2rvKY/wMw3jX8K/CAxzwTkA40LeSzuGPHqOBOAyaP+d8Dr7qGpwEzPOZdAcQClgLWOwZYk2c7hcW8AnjYY951rvL51lvA3/JVYLnHvA5ARhF/Ow0kA2mu4e8Bb9e8/cAjF/jb5yyf4PodvZnzfQHtgU2A3VVuGqAKWU89V6wmjCOM1cD/Kvr/pLq/5EhBFKYBxt5tjqOuaQVSSrVRSi1yNR4mA//BOGooznbitdYpebbV0GPdl2Hs2Q7XWscWFJ/WOhWI81wOo/LOkY6ROMBoI/jIdYomEYgHVJ5li4r3uM69p5wrXuC4x3Bj4KjW2l6MdRcVc4M86/Ucvpj1+hR26smlu2vbo4BeQM4eehzGkdmFdNdah2qtW2qtX9JaO5VSJoyjgvmu9YUDocA7Ba1Aa31aax2ttXZqrQ8DzwG3FGPb4hJIUhCFOYlReeZo4poGxh5eXpOBPUBrrXUQxukZVUC5grYTppQKzLOtEwBKqTrAT8A/de4rf3LF5zqtUCtnuQs4jrHXHeLx8tVary1mvI1dFVy+eF08v5/jQJMLVMDFcQpo5DHe+BLXd0Ha8APGKbAJrsnLgeF5Pn9xhWF8V59q49RcHDAVGFzckJA6q8zJFywK8z3wklKqtlIqHKNSmOmadwaopZQK9igfiHHKINV1HnpscTaitT4OrAX+Tynlo5TqAjwAzHRVpHOBma7KKW989ymluimlvDGOTP7WrsbPC/gceF4p1RFAKRWslBpZnHiBvzH2tJ9TSlmVcS/ATRTeTrABo0J/Wynl7/qMvYu5LU8/AE8qpRoqpUKAf1/EOi7W28CDSql6wIdAEDA9p3HeFdOHrr9dobTW54DDwFillMX1OUYD+RqoXevtr5RqqgyNXXH8XGqfShRIkoIozJtAFMY/7A5gs2sa2rii5HvgkOsUTAPgWYwG1hTgC2B2CbZ1B0Yj7EngR+AVrfVyjD3jPsBTKve16k1c818G5mFUui2B24uzMa31jxinLGa5TnXtxGjwLM6y2RhJ4AbgHMYllve6vpOCyjtc5VthXPIbg3FKpqS+wGjc3o7RaLsY47y84yLWVSJa6x0Y5/PHaa3jgSsx2m/+VkqlAL8DSZxvjC7KCGAQRjvLAdd6ns6Z6fr79nGNXoaxw5Dmet8BPFEan0kUTmktD9kRoqpRSt0AfK61zncprRCXQo4UhKgClFK+SqnBrtMuDYFXMI6qhChVcqQgRBWglPID/gDaARnAL8CTWuvkCg1MVDuSFIQQQrjJ6SMhhBBul3rtdIULDw/XzZo1q+gwhBCiStm0adM5rXXtvNOrfFJo1qwZUVFRFR2GEEJUKUqpowVNl9NHQggh3CQpCCGEcJOkIIQQwk2SghBCCDdJCkIIIdwkKQghhHCTpCCEEMKtyt+nIIQQlZ3Dqcm2O8m2O8lyOMi2O7E5NDaH0zVsjNsdTrId5+flzLc7Pcuen3dXr6bUC/Yp1VglKQghqiWtNVl2J5k2B5k217vdGM6poLNdFXSW63V+ev7hrFzjjoLn5QznWd7hLP0+5pSC69rXlaQghKjabA4n6VkO0rLtpGfbSctykJ5tVNhZnhW4zUGm3UmGe57HdJuTTLvDNc/pWu789JzhS2VS4GUx4WU24WUx420xeYyfH/bzs7jHvV2v3GXM54ctJrxd86xmE1azcr27hl3LWs0mLGblHvacZzEpzCaFUsV54m3JSFIQQhSL1prULDvxadnEp2WTkJ5NXKrxnpJpJz3b4VHJe7xnO0jPMt4zsh1kO0peWXuZTfhYTfhYza6Xa9hiJsjHgm+gd+7pVjM+FhPeVjO+eZex5qmkzecr85xxb6ur8jXXvGZXSQpC1FDZdicJ6a4KPi2buDwVfU7ln5MAEtJshVboSoG/lwU/LzP+3q53Lwshfl40DDXj52XB38uMn7fx7pt33GrG21Vhn6/EjXFvixmzqfT3iEXBJCkIUQ1orUnOsBOfnk18WhbxabZcFb1nBZ+TBFKy7IWuL8TPSpifF6H+XjQO86NroxDCArzc02r5G+/GuJUAb0uZnMoQ5U+SghCVUKbNUeBee0Jatqvi93zZSEjPLrQx09tiOl+J+3vRtJYfYf4FVPCuV4ivtUaeNhGGcksKSqmvgSHAWa11pwLmK+AjYDCQDozRWm8ur/iEKA1aa5Iz7SSl20jOtJGcYSMpw0ZKpp3ULDtpWXZSs+2kZrqGsxykZtlIy3K4xo1XerajwPUrBaF+Xu5KvXm4PxFNvQnztxLq50WtAC/j3d+bUH8rYf5e+HnJvp8ovvL8tUwDPgVmFDL/BqC169ULmOx6F6JSSs2ys+5gHNtjEtl1Mpnj8emcSMwotELP4WUxEehtwd/1CvS2EB5g7MEHeFsI8La499w9K/owfy+Cfa1yfl2UqXJLClrr1UqpZkUUGQbM0MZDo9crpUKUUvW11qfKJ0IhLkxrzd+H4/lqzWH+2BtLtsOJ2aRoVTuA5uH+XNU6nAbBvoT4WQn2tRLka7wHeFsI9DGSgFVOzYhKrDIdVzYEjnuMx7imSVIQlcLx+HSen7+DNQfOUcvfi3uuaMp17evSrXEIvl7mig5PiFJRmZJCsSmlHgIeAmjSpEkFRyNqgvWH4nhoRhRODa/c1IE7ejbBxyqJQFQ/lSkpnAAae4w3ck3LR2s9BZgCEBkZWfr3jwvh4c/9sTwwLYomtfyYOqYHjcP8KjokIcpMZTq5uQC4VxkuB5KkPUFUtE1H43loxiZa1PZnzsNXSEIQ1V55XpL6PdAPCFdKxQCvAFYArfXnwGKMy1EPYFySel95xSZEQaJPJjNm6kbqBnkz44GehPp7VXRIQpS58rz66I4LzNfAP8spHCGKlJCWzYMzovD3sjDzH72oE1i6PVEKUVlVpjYFISoFrTXj5m4jNiWLOY9cQaNQOWUkao7K1KYgRKWweMdplu8+y3OD2tK1cUhFhyNEuZKkIISHpAwbry3cRccGQYy5sllFhyNEuZPTR0J4eHXBLuLSsvlydKR0CidqJPnVC+GyeMcpftxygsevaUWXRiEVHY4QFUKSghDA2eRMXvxxB10aBfPP/q0qOhwhKowkBVHjaa0ZP38H6dkOPrytm3RYJ2o0aVMQNd7sjcdZsecsr9zUgVZ1Asp121prlFI4s7JwJCXhTEnBmZ6BMz0dR9w57PEJaLsNHE6caWnY4+NwJCUZTzkzmY0HLDjsKKsVU2AQJh9vtFOD0wlao6wWlJc3mIxEZwrwxxwcgjk4CJOPD8rXF0toKOaQEJSfH8pqlSeo1XCSFESNdjw+nTcWRXNly1qMvqJZqaxTa4397Fns587hSEgk++gRsg8eJPvECZyJSTiSjJczNRVtsxkVtrMYD7NXCnNoKOagIGM7ORW/xYLOzsaRkoLOzASzkSwUoO12dHZ28YM3m7HWrYulXj1MAf6Y/Pwx+flh8vc3EkpAAMrHB3NQMNaGDbE2bOCepszSQWB1IElB1FhOp+bZOdtQSvHurV0wXeTDa7TTSebOnaT99RcZ27aTsWMHjri4XGVMgYF4NW6MOSQEa8MGmAKDMAcGoLy80U4HJh9fY+89KMiohH39sITXwhwairJawWTG5HtxFa/RWQCgNc6UFBzJyTgSk9BZmTjT0nAkJuJITMSZYYzbTp/GfvYsjoREbDEncKan40xLw5mWVmTyMgUHY6lVC3NgICZ/P8y1wrHWr485LBSTr5FYLHVqYwkPx+QfYHx+X185MqlkJCmIGuub9Uf5+3A8797apcR3LWu7nbS//yb1999J+X0F9jNnQCm8WrYgoE8ffDp1wlq/HubgYKyNm2CpU7vCKj/3dpXCHByMOTgYGjcueqECaK3R6ek4MzNxJCRgO3kS28mTRsJIz8CRkIA9Ls51Ciyd7KObST59GhxFPInOZMIUGIg5MBBzUBCm4CDMQcGYgwIxBbmGg4OMMq5ho1ww5sBAlEWqsNIm36iokc4kZ/Lekr1c3aY2IyMaFXs5bbcTP3068dNnYD97FuXrS8BVVxF43bUE9O2LOSSk7IKuYEoplL+/scdfqxberS58lZZ2OIykkZGBMyUF25kzOBIScaam4kxLxZGaijPZOHpxJifjSE4m6+wBHMlJOJOSL3jqy+TnZySIoCDjCCVnOCfBBAa5j8ByJZvgYEze3qX11VQrkhREjfTGomiyHU5eH9qx2Hvw2ceOceKZf5G5cyf+vXtT9+WXCOjTB5OPdJZXGGU2uytp6tYtViLx5G6AdyUMd/JISjYShyuhOJKTcSYlYTt+nMyUFJxJSTjT04uOzdvbOALx98cUEOBqNwlwt52Y/AM8prumuacbydEcUP1OgUlSEDXO6n2xLNp+iqeva0OzcP9iLZO5ezfH/vEg2O00nPhfggYNKuMoBYDJ2xtTnTpQp06Jl9U2G46UlPMJJSeRpKScH05NM45aUlNxpKUap8Pc42lgsxUjSNP5ZOLvbyQNf39Mfr6Y/PxQvr5Gm4qfHyY/33zjJl+Pcn5+xrivr9GWVAEkKYgaJdvuZMLPO2kR7s8j/VoUa5n0qCiOj30UU0AATb6ZgXeL4i0nKpayWrGEhUFY2EWvw5md7U4SzlTX6a7UNFfDe55prlNizrQ0HCkp2M+ecV9e7MzIMK4MK2H8ys8vV6Iw+fqi/I0LEUw+PoQ/OhavUn4ksSQFUaPMXH+UI3HpTL2vB96WC1/Jk7Z+PccffgRrgwY0+epLrA0alEOUorIweXlhusTEkkM7neiM80nCmZGBMy0dZ0a6Md09nnF+WppH2fQ0dHoGjnNx2DJicGZmEHZf6T+LTJKCqDGSMmx8smI/V7UKp1+b2hcsn7lnDzGPPY5XkyY0mT7N2OsU4iIpk8ndUF+Zyf38osaYsvogiRk2xt/Q7oINg7ZTpzj+0MOYAgJo/MUUSQiixpAjBVEjJGfamLH2KIM71adTw+Aiy2qHgxPPjsOZlkbT77/DWq9eOUUpRMWTpCBqhG/XHyMly87Yfi0vWDbuiy/I2LSJBu++g0+bNuUQnRCVh5w+EtWew6mZvvYIV7UKv+BRQsb27cR+8ilBN95I0E03lVOEQlQekhREtffn/lhOJ2dyV6+iL91zpqVx4tlxWOrWod4rE6rVDUlCFJecPhLV3pxNMYT6Wbm2fd0iy53+z3+wxcTQdMZ0d0+kQtQ0cqQgqrXE9GyW7TrDsG4N8bIU/nNP+uUXkubNp9aDD+IXGVmOEQpRuUhSENXawm0nyXY4ubWITu+yjxzh9MsT8L3sMmo/9s9yjE6IykeSgqjW5myKoX39oEIbmJ1ZWcQ89TTKaqXhhx9UWH8zQlQWkhREtXUoNpXtMUnc0r1hoWXOffIJWXv2UP/t/8Nav345RidE5SRJQVRbi7afQikY0qXg/orSt2wh7uuphIwcSWD//uUcnRCVkyQFUS1prVmw7SQ9moVRLzj/8w6cmZmcev4FLPXqUuffz1VAhEJUTuWaFJRSg5RSe5VSB5RS4wuY30QptVIptUUptV0pNbg84xPVx57TKRw4m8rQrgUfJZz7bDLZR47Q4M03MQcElHN0QlRe5ZYUlFJmYBJwA9ABuEMp1SFPsZeAH7TWlwG3A5+VV3yielmw7SRmk+KGTvn7Lcrav5+4r78mePhw/K+8sgKiE6LyKs8jhZ7AAa31Ia11NjALGJanjAZy7hoKBk6WY3yimtBas2j7Sa5qFU6tgNzP4dVOJ6deeRVzQAB1nhtXQREKUXmVZ1JoCBz3GI9xTfP0KnC3UioGWAw8XtCKlFIPKaWilFJRsbGxZRGrqMKiTyVzPD6DwZ3zHyUkzp1LxubN1HnuOSyhoRUQnRCVW2VraL4DmKa1bgQMBr5RSuWLUWs9RWsdqbWOrF37wg9LETXL8uizKAXXtMvdrYX93DnOvv8Bfj16EDz85ooJTohKrjyTwgmgscd4I9c0Tw8APwBordcBPkB4uUQnqo1lu0/TvUkotQNznzo6+9//4szIoN5rr0pnd0IUojyTwkagtVKquVLKC6MheUGeMseAawGUUu0xkoKcHxLFdiopg50nkrkuT+d3WYcOk/TjT4TdeQfeLVpUUHRCVH7llhS01nbgMWAJsBvjKqNdSqnXlVJDXcX+BTyolNoGfA+M0Vrr8opRVH3Lo88AMKBD7qRw7tNPUT4+1HrooYoIS4gqo1y7ztZaL8ZoQPacNsFjOBroXZ4xierlt12naRHuT8va5x+Objt9muTffiNszBgstWpVYHRCVH6VraFZiIt2NiWTdQfjGNKlfq42g8QffgCtCb3zjgqMToiqQZKCqDZ+2X4Kp4abPO5i1tnZJMyZQ8DVV+PVqPDus4UQBkkKotpYuO0k7eoF0rpuoHtayu+/44g9J0cJQhSTJAVRLRyPT2fzsUSGdsvd11HCd99jbdQI/6uuqqDIhKhaJCmIamHhdqNHlJs8usnOOnCA9I0bCb19FMpsrqjQhKhSJCmIKk9rzU9bTtC9SQiNw/zc05N+/hksFoJHjKjA6ISoWiQpiCpvx4kk9p1J5daI8zfMa61J/mUx/ldegSUsrAKjE6JqkaQgqry5m2Lwtpi4scv5x2lmbNmK7eRJgm+8sQIjE6LqkaQgqrRMm4Oft55kYMd6BPta3dOTFy1CeXsTcO11FRidEFWPJAVRpf2++yxJGTZGRp6/B0Hb7SQvWUJAv36YA/yLWFoIkZckBVGlfbfhKPWDfbiy5fnOdNPW/40jLo6gG+VprkKUlCQFUWWt2nuWvw7EcV/vZphN57u1SP7lF0wBAQT07VuB0QlRNUlSEFWSzeHkjUXRNKvlx5grm7unO7OySFm2jMABAzB5exexBiFEQSQpiCpp6l+HORibxos3dsDLcv5nnPrHHzhTUwmSq46EuCjl2nW2EKXhwNkU3l+6jwEd6nJd+zq55iX/shhzrVr4X96rgqITwkVr0E5w2s+/HPbc4047OB2FjNtc7651aIfHfIcx3vp68Cvd+3AkKYgqJTnTxj+/3YK/l5n/DO+cq4tsR2oqqatWEXLLLSiL/LSrHK2Nys6RbVSIjpxXtvHuzBm2u95LUs5+fn5O2ZxlC6ugHbYiKmzPcVvh88vagyskKYia7cnvt3AwNpXp9/fM9wzmlOXL0VlZBA0ZUkHRVSNagz0LslIgO8V4z3llp4E90/XKcr1n5xnPAkdW7nH3Mtm5y3hW6JThgxaVGcxerpfFeDdZjWFTzsvsMex6WbyLnu85zWwtoExB49ZirjOnrMkYV57TzRDY4MKfu4QkKYgqY3tMIiv3xjL+hnb0bhWeb37yL4uxNmyI72Xdyj+4ysLpgOzU3JV4Qa+8FX1BL6etZNs2exsVqMUbLD65383e4BUAfuF5pnudr3TNXkalarZ6VNge42ara5pHpW72Kt6yORWruCBJCqLKmLn+KL5WM3f2apJvnj0ujrS1a6n1wAO5TilVCU6HsfftfqVAlmfFnly8ij4rBWxpxdum1Q+8A42K2jvQeIU0BW+Pce9A8A7KXSZnmbwVv9lLKt1qQpKCqBKS0m0s2HaS4Zc1IsjHmm9+8m+/gcNR/lcdOZ2QcgoSj0JSDGQkQEai8Z6ZBLZ041SJLcPjPQvsGZCd7joVk1G8bSkz+ASBl0cF7VcLQpu5KvOg3BV3TqXurtBd07wCjT1tIQogvwxRJcyOOkamzck9lzctcH7yL4vxbt0an7ZtyjYQhw0Or4b9yyBmI5zeYZwXz8s7CHyCjT1yq49rr9oHfEONd6uvMc/L36jAvfxzD+faW3e9LD5Q1Y6CRJUjSUFUeg6nZvrao/RqHkaHBkH55ttOnCBj82ZqP/102QWRmQQbvoD1kyH9HFh8ocFl0PNBCGth7K0HNzb23H2CZU9cVFnyyxWV3rLoM5xIzODlIe0LnJ+0eDFA2fR15HTCtu9g2StGMmg1ACLvh5bXGEcAQlQzkhREpTdt7WEahvhyXfu6Bc5PXvQLvt264dWoUYHzL9rZ3bDgCYjZAI16wl0/QMOI0t2GEJWMXC4gKrXdp5JZfyiee69oisWc/+eatX8/WXv3lu69CVrD2k/gf1dD/EG4+XO4f4kkBFEjyJGCqNSmrz2Cj9XEqB6NC5yf9MsvYDYTNGhg6WzQng2/joNN06DdELjpI/DPf0+EENWVJAVRaaVk2vhxywlGdG9EiJ9XvvnaZiPpp5/xv+IKLOGlUHEnxcCc+4zTRVc9A9dOkKt9RI0jSUFUWn8dOEeW3cnwyxoWOD9l2TLsp09T75UJl76x/cth/oPGJacjp0HH4Ze+TiGqIEkKotJauSeWQB8L3ZuEFDg/fvoMrE2bXNrDdJxOWP0erPo/qNMBbpsB4a0ufn1CVHHl2tCslBqklNqrlDqglBpfSJnblFLRSqldSqnvyjM+UXlorVm59yxXt6ldYANzxtatZGzbRtjd96AutnuFzGSYfRes+g90GQX/WC4JQdR45XakoJQyA5OAAUAMsFEptUBrHe1RpjXwPNBba52glKpT8NpEdbfrZDJnU7Lo37bgn8C5KV9gCgoiePhFnuZJjYWZI+DMLrjhXej5kLQfCEH5Hin0BA5orQ9prbOBWcCwPGUeBCZprRMAtNZnyzE+UYms2mv86fu2qZ1vXmZ0NKkrVhB2772YA/xLvvLE4zB1EJzbD3fOhl4PS0IQwuWCSUEpNUAp9YVSqptr/KGL3FZD4LjHeIxrmqc2QBul1F9KqfVKqUGFxPSQUipKKRUVGxt7keGIymzl3li6NArO98wEgHOTJ2MKDCTs3ntKvuLYffD1QONI4d6foPWASw9WiGqkOEcK9wPjgLuVUtcA3cowHgvQGugH3AF8oZQKyVtIaz1Fax2ptY6sXTv/nqSo2hLSstlyLIF+BZw6yty7l5Rlywm75x7MQfn7QSpS7F7jCMGRDff9Ak0uL6WIhag+ipMUUrTWiVrrZ4HrgR4Xua0TgOcdSI1c0zzFAAu01jat9WFgH0aSEDXIsugzODUMKKBbi7gvvsTk70/Y6HtLttLkk/DNCOOBLPcvgXqdSylaIaqX4iSFX3IGtNbjgRkXua2NQGulVHOllBdwO7AgT5mfMI4SUEqFY5xOOnSR2xNV1K87T9Eo1JdODXMfCdjOnCX5t98IvmUE5uDg4q/Qng3f3270dHrXHKjVspQjFqL6uGBS0Fr/DO5KGq31JxezIa21HXgMWALsBn7QWu9SSr2ulBrqKrYEiFNKRQMrgXFa67iL2Z6omuJSs1hz4ByDO9fP9wS1hO+/A4eDsLvvLtlK/3gHTm2DEf+D+l1LMVohqp+SXJL6NTD0gqWKoLVeDCzOM22Cx7AGnnG9RA00b3MMNodmZETuHk+dmZkkzppNwDXX4NUk/+M4CxUTBWs+hG53QbtyfiqbEFVQSS5JlWv2RJnSWjNr43EimobSum5grnnJixbhSEwk7N4StCXYs+CnsRDYAAb9XylHK0T1VJKkoMssCiEwGpgPxabl6xFVa0389Bl4t22LX88SXOew9hM4tw9ummg8DU0IcUFypCAqhcPn0vj3vO20qxfIzd1y376S/vffZO3fT9i99+ZrZyhUwlFY/T60Hyr3IghRAiVpU3i+zKIQNZbWmm//PsZ/Fu/GajYx+e4IvCy591Xip8/AHBZG0JAStAn8Nh6USU4bCVFCxU4KWuudSqlI4EWgqWtZZczSXcooPlGNJaZn8/j3W/hz/zn6tA7n7Vu60DDEN1eZ7KNHSV21ivCxj2Dyzn93c4EO/wl7F8N1r0JwKT+iU4hqrqQd4n2LcXfzDsBZ+uGImuT1hdGsPxTHmzd34q5eTQo8NRQ/81uwWAi5/fbirVRrWP4KBDWEXmNLOWIhqr+SJoVYrXXeG86EKLG1B88xf8sJHuvfirsvb1pgGWd6Oknz5hE8+AasdYrZYe7uhXBiEwz9FKw+pRixEDVDSZPCK0qpL4HfgayciVrr+aUalajWsuwOXvpxJ03C/HjsmsKfX5C6Zg3O9HSCh48o3optmbBsAoS3ha53lFK0QtQsJU0K9wHtACvnTx9pQJKCKLb/W7yHQ+fSmH5/T3ys5kLLpf6+AnNwMH6REcVb8bpPIeEw3PMjmOWhgkJcjJL+5/TQWrctk0hEjfDNuiNMW3uEB65qXuCzEnJou52UVasI7NcPZSnGzzQpBv78ANoNgZbXlGLEQtQsJX3IzlqlVIcyiURUe+sPxfHqwmiuaVeHFwa3L7JsetQmnElJBFx3bfFWvuRF0E4Y+J9SiFSImqukRwqXA1uVUocx2hTkklRRbB8u20fdQG8+vuMyzKaib0JL+f13lLc3Ab17X3jFexZD9E9wzUsQWnCjtRCieEqaFAp8EpoQF7LzRBIbDsfz4uD2BHgX/bPTWpPy+3L8e/fG5OdX9Iozk+GXf0GdjnDlk6UYsRA1U4mSgtb6aFkFIqq36WuP4Odl5rY8/RoVJHPHDuwnTxH4z8cuvOJlL0PKKRg1EyxepRCpEDVbSdsUhCix+LRsft52kuGXNSTY13rB8ok//ojy8SHw+gv0WbThC9g0DXo/AY2KeYWSEKJIkhREmfsh6jjZdiejr2x2wbLOtDSSf1lM4HXXYQ4MLLzgttnw63PQ5ga49pXSC1aIGq5YSUEp9ZRSqqdSSi7+FiW2cNtJujUOoU3dIip5l8R583AmJxN2912FF9r8Dfz0CDTtDbd+DabC73UQQpRMcY8UGgETgbNKqT+UUv9RSg1RSoWVXWiiOjgWl86uk8nc2Ln+Bctqu534adPxjYjAt1u3ggstfQkWPAbNr4Y7fwCvCzRECyFKpFhJQWv9rNb6SqAeRhfa8Rh3N+90PU9ZiAIt230GgEGd6l2wbPz0GdhOnqTWAw8UXCD6Z+PBORH3wV3zJCEIUQZKejrIFwgCgl2vkxg9pgpRoLUHztE83J/GYUVX4GkbNnD2v/8l4LprCejfL3+BlNOw8Cmo3w0GvyfdWAhRRor1n6WUmgJ0BFKAv4G1wIda64QyjE1UcXaHk78PxzOsW4Miy6WtX8/xsY/i1bgxDd56K38X2lrDz4+BLR1GfAHmC1/BJIS4OMVtU2gCeAOngRNADJBYRjGJamJbTBKpWXZ6twovtEzyb0s4/vAjeDVqRNNvZmAOLuBZypumwoFlMOB1qN2mDCMWQhTrSEFrPUgZu28dgSuBfwGdlFLxwDqttVwTKPJZe+AcSsEVLWrlm6e1Jn7qNM6++y6+3brRaPJnWEJD868k7qDRr1GL/tDjwXKIWoiarSSP49QYDcuJQJLrNQToCUhSEPn8dfAcHeoHEeqf+05je3w8p195hZRlywkcOJAG775T8KM2nU748WHjdNHNn4FJbqsRoqwVt03hCYwjhCsBG0abwlrga6ShWRQgI9vB5qOJjL4ydwd19oQEjtw6EntsLHXGPUvYmDEocyH3GexfCjEbYdgkCCq6XUIIUTqKe6TQDJgDPK21PlV24YjqIupoPNkOJ1d6tCdkHTjAiaefwXbyJE2//w6/yy4reiXrP4PABtBlVBlHK8qL1hqnduLUThzaUfC7s5Dpnu/OwufnWwfOXNvVGMNaa/dwzrgTj2FX2YKme67TXSbPdM9tObXxTLKCpnuuH02R23LiKuMafrHXizQPbl6qf6Pitik8U6pbFdXe2oNxWEyKns2M+xuzDh3m6L2jccTHE3zrLRdOCGd2weE/jC4s5GqjYtNaY9d2sh3ZZDmyyHZknx92ZucezzvsNIZtDhtZjixj2GkrdD2e021OGw6no+iK3KNyrE5MyoQJE0opY1iZUBjDSin3cN7peZfLNa+w9eVZ3jirX7rkYm9RJtYeOMdlTULw97ZgT0jg2P33g1I0mz0Lny7FePzG+slg8YWIMe5JZ9LO8O2eb2kc2JgRrUZgrmTdW2ityXZmk2nPLLICLW4FXZJK23O65tIrCqvJipfZC2+zN1aTFW+zN15mL/c0L5MXAdYAY77ZitVkxWKyYFImzMpc4LtJmTCbCp6X6910gfke6ynOOkyYMJlc7yWoqHMq56Iq8XyXT1cD5ZoUlFKDgI8AM/Cl1vrtQsrdAszFePxnVDmGKEpBapadHSeS+Gf/VgDETvwIe2wszWbPxrdTx2Ks4Cxs/wEuuwv8wki3pTN111Sm75pOpj0TjWbO3jm80OsFutXpVuAqtNbEZ8ZzOOkwx1KOkZqdmq/CTrOlkZyd7H5l2jPdyysUflY//K3+BFgD8LP6YXfaybBnkGnPNF6OTDLsGeenOTIveU/YpEzuCtjbZFS43mZvd+XrbfYm0CuQWuZauSpqL5NXvoq7sMrcPWzOP5wz36SkUb+mKrekoJQyA5OAARj3OWxUSi3QWkfnKRcIPIlxk5yogjYfTcCpoWfzMOwJCSTOnUvoqFHFSwgAf30EThtc/k8WHlzIh5s+5FzGOQY1G8ST3Z9kZ9xO3tv4Hvf8eg/DWg7j+mbXcyL1BCdSThjvqSeISYkhxZZS4Op9zD54mb3ws/oR5BVEoFcgjQIa4Wvxde/5ObWTDFsGqbZUzqafJc2WhtVkxcfig4/FhxCfEHwtvviYfYx31/ScaTmVrWdlXlCFnLcyt5jk4F1UrPL8BfYEDmitDwEopWYBw4C8fSe9AbwDjCvH2EQp2nA4HrNJ0b1JKCk/zgOHg5CRtxZv4YSjsPFL6HwbO8nkhTUv0KV2Fyb2n0jX2l0BaBTYiKsbXs3/tv+PGbtm8PPBnwHwNnvTIKABDQMa0qV2F5oFNaN5cHOaBjUl2DvYvfdcHQ/5hSgt5ZkUGgLHPcZjgF6eBZRS3YHGWutflFKSFKqoDYfj6dQwGH9vC7E//YR361Z4t2tXvIV/Gw/KDNdO4KvN7xPoFciUAVPwt/rnKuZn9ePpiKe5re1txKbH0jCgIbV8a8lpDyEuUaX5D1JKmYAPMe6WvlDZh5RSUUqpqNjY2LIPThRbps3B1phEejYLJevQYTK2biV42LAL751npcKv42HvYuj3bw6Rxe/HfueOdnfkSwieGgY0pFudbtT2qy0JQYhSUJ7/RScAzwf0NnJNyxEIdAJWKaWOAJcDC5RSkXlXpLWeorWO1FpH1q5duwxDFiW1PSaJbLuTns1rET9jOspqJXjYsMIX0BqiF8CknvD3ZKMri8v/ydc7vsbb7M1d7Yt42I4QotSV5+mjjUBrpVRzjGRwO3BnzkytdRLgvtNJKbUKeFauPqpaNhyOA6B7kCb2x58IvnkYloISty0Djq2Dlf8HMRugbie4dSo06cWp1FP8cugXbmt7G2E+8hwnIcpTuSUFrbVdKfUYsATjktSvtda7lFKvA1Fa6wXlFYsoOxuOJNC2biB64Y/orCzCxow5PzMr1ejtdN1ncCIKtNO4Y3nIRLjsbvdNatOjpwMwpuOYfOsXQpStcr3+TWu9GFicZ9qEQsr2K4+YROmxO5xsOhLPiG71SfhkDv5XXoF3y5bGTIcdZt0Bh1eDTzD0+ZdxdNBmIFh93etIyExg3r55DG4xmPoBF36EpxCidMlF0aLU7D6VQlq2g75JB7GfOkXd8ePPz1zxhpEQut0N/V+A4IYFruPb3d+S6cjkgU6FPJJTCFGmJCmIUvO3qz2h2V9LcNYOJ/Ca/saMY+uNG9IixsBNHxW6/PGU48yInsGApgNoEdKiHCIWQuQl1/CJUrPhcDxdvLOwrV1DyIhbUFZXR3Yr3oTA+nD9W4Uuq7Xm1bWvYlZmnuvxXDlFLITIS5KCKBU2h5N1B+MYmbALnE5CbhlhzDizC478Cb0eBu+AQpefu38uG05v4F+R/6Kef71yiloIkZckBVEqNh9NICXTRpddf+HbvTteTZoYMzZ+CRYf6H5vocseSjrEB1Ef0KteL25pfUs5RSyEKIgkBVEqVu2LpVPiUbxOHCV4+M3GxKwUo7fTjiPAr+D7Dbae3co9i+/B2+zNK1e+Iv0SCVHBJCmIUrFqbyx3xm7B5OdH8ODBxsRtsyA7FXr8AzDaDVKzUwGwOWz8cugXHln+CCHeIcwcPJPGgY0LW70QopzI1Ufikp1JzuTY0TN02R9F0PBhmPz9je4r/v4fNIyARhEAvLPxHWbvnU3/xv3ZcnYL5zLO0Sa0DZOunSTtCEJUEpIUxCX7Y28s/WI2Y7ZlETJypDHxyBqI2w83f26UOf4H3+7+FoDtsdvpUKsDo9qOoneD3pXuCWpC1GSSFMQlW7XvLEOPb8S7fTt8ch6ks3kGeAdDh2HEpsfy8l8v0y6sHd8O/hYvs1fFBiyEKJS0KYhLYnM4ObZhK00TYgi59VajoTgjEXYvgM634rT68NJfL5Fhz+CdPu9IQhCikpOkIC7J5qMJXHFgA9psISingXnHHLBnQvd7mRk9k7Un1zKuxzi5S1mIKkBOH4lLsmrvWfqc3I5vnz5YQkONiZtnQL3O7PXxZeLmifRv3J+RbUZWbKCFsNlsxMTEkJmZWdGhCFEmfHx8aNSoEdacHgYuQJKCuCTbo/YwNCOR4D69jQmnd8Lp7XDDe3y46UMCvQJ57crXKu39BzExMQQGBtKsWbNKG6MQF0trTVxcHDExMTRv3rxYy8jpI3HRziRn4rt7GwB+kT2MiTt+AJOF482uYO3Jtdze7nZCfUIrMMqiZWZmUqtWLUkIolpSSlGrVq0SHQlLUhAX7Y+9sXQ+dxCCgvBu3QqcDtgxD1pey7yY3zEpEyNajajoMC9IEoKozkr6+5akIC7aH3tOccWZPQRedRXKZIJDqyA5BlvnW/jxwI9c3ehq6vrXregwhRAlIElBXJTULDsZK1cSmJVK8JAbjYlbZoJvKCsDAojPjK+0jcvVyeDBg0lMTCyyjNaaF198kTZt2tC+fXs+/vjjEm1jzJgxzJ07N9/0VatWMWTIkHzTt27dyuLFi/NNL4kJEyawfPnyS1rHI488wl9//XVJ6wAYN24cHTt2ZNy4cZe8rn79+hEVZTx2PiCg8F6Di+vIkSN06tTpktfjSRqaxUWZE3WcgXtXo+vVJ6BvX0iPhz2LIOI+5h1cQH3/+vRu0Luiw6z2ilP5Tps2jePHj7Nnzx5MJhNnz54t05i2bt1KVFQUg3MuUb4Ir7/++iXHsX79eiZNmnTJ65kyZQrx8fGYzTXjzntJCqLE7A4nvy1Yw6txh6gz7lmU2Qw75oIjm9Ptb2Dd6id5uOvDVa77itcW7iL6ZHKprrNDgyBeualjkWVmzJjB+++/j1KKLl268M033zBmzBh8fX3ZsmULZ8+e5euvv2bGjBmsW7eOXr16MW3aNACaNWtGVFQUqampDBo0iIiICDZv3kzHjh2ZMWMGfn5+TJ48me+++w6TyTgxUKdOnQLj2Lp1K4888gjp6em0bNmSr7/+mtDQ3BcJ/Pbbbzz11FP4+flx1VVX5VtHdnY2EyZMICMjgzVr1vD888+ze/duAgICePbZZwHo1KkTixYtAuCGG27gqquuYu3atTRs2JCff/4ZX19fxowZw5AhQ7j11ltp1qwZo0ePZuHChdhsNubMmUO7du2IjY3lzjvv5OTJk1xxxRUsW7aMTZs2ER4ezu7du2nTpg1ms5l+/frRq1cvVq5cSWJiIl999RV9+vQhMzOTsWPHEhUVhcVi4cMPP6R///65Ps/QoUNJTU0lIiKC559/nmuuuYZHHnmEY8eOATBx4kR69+5NWloajz/+ODt37sRms/Hqq68ybNgwMjIyuO+++9i2bRvt2rUjIyMj1/qffvppli5dSr169Zg1axa1a9fmiy++YMqUKWRnZ9OqVSu++eYb/Pz8OHPmDI888giHDh0CYPLkyTRo0MC9rkOHDnHLLbcwZcoUevToUeRvrihy+kiU2G+7TtNr6+84vbwJucX1/IMt30C9LixKPYhGM7TF0IoNsorYtWsXb775JitWrGDbtm189NH5x5UmJCSwbt06/vvf/zJ06FCefvppdu3axY4dO9i6dWu+de3du5dHH32U3bt3ExQUxGeffQbAwYMHmT17NpGRkdxwww3s37+/wFjuvfde3nnnHbZv307nzp157bXXcs3PzMzkwQcfZOHChWzatInTp0/nW4eXlxevv/46o0aNYuvWrYwaNarIz79//37++c9/smvXLkJCQpg3b16B5cLDw9m8eTNjx47l/fffB+C1117jmmuuYdeuXdx6663uihrg119/ZdCgQe5xu93Ohg0bmDhxovtzTZo0CaUUO3bs4Pvvv2f06NH5rtJZsGABvr6+7s/y5JNP8vTTT7Nx40bmzZvHP/5h9AD81ltvcc0117BhwwZWrlzJuHHjSEtLY/Lkyfj5+bF7925ee+01Nm3a5F53WloakZGR7Nq1i759+7rjGjFiBBs3bmTbtm20b9+er776CoAnnniCvn37sm3bNnfiz7F3715uueUWpk2bdkkJAeRIQZSQ1pqZS3fwfMxmQocPwxwS4r43QQ96h58P/Ez3Ot1pHFT1usG+0B59WVixYgUjR44kPDwcgLCw88+duOmmm1BK0blzZ+rWrUvnzp0B6NixI0eOHKFbt2651tW4cWN69zZO2d199918/PHHPPvss2RlZeHj40NUVBTz58/n/vvv588//8y1bFJSEomJifTt2xeA0aNHM3Jk7jahPXv20Lx5c1q3bu3expQpUy7p8zdv3tz9OSIiIjhy5EiB5UaMGOEuM3/+fADWrFnDjz/+CMCgQYNyHdUsWbKEqVOnFrh8zjbWrFnD448/DkC7du1o2rQp+/bto0uXLoXGu3z5cqKjo93jycnJpKamsnTpUhYsWOBOWJmZmRw7dozVq1fzxBNPANClS5dc6zaZTO6keffdd7tj3LlzJy+99BKJiYmkpqYycOBAwPitzJgxAwCz2UxwcDAJCQnExsYybNgw5s+fT4cOHQqNvbgkKYgS2XA4nobrluHlsBF2z93GxO2zwGRhe4MOHNk7ifs63VexQVYT3t7egFF55AznjNvt9nzl8156mDPeqFEjd4UzfPhw7rvP+PsMHDiQM2fOEBkZyQcffFAmnwHAYrHgdDrd4557456fy2w25zu9krec2Wwu8LN7Sk9PJzExMdeplZIsXxSn08n69evx8fHJNV1rzbx582jbtu1Frzvn7zVmzBh++uknunbtyrRp01i1alWRywUHB9OkSRPWrFlTKklBTh+JEvli1QGGHlmHT2QkPm3bGs9N2PUztBrAzydW4WP24fqm11d0mFXGNddcw5w5c4iLiwMgPj7+otd17Ngx1q1bB8B3333nPud/8803s3LlSgD++OMP2rRpAxh701u3buXLL78kODiY0NBQ9xHEN9984z5qyNGuXTuOHDnCwYMHAfj+++8LjCMwMJCUlBT3eLNmzdi8eTMAmzdv5vDhwxf9GT317t2bH374AYClS5eSkJAAwMqVK/O1DRSkT58+fPut0Z37vn37OHbs2AUr9euvv55PPvnEPZ5zGm/gwIF88sknaK0B2LJlCwBXX3013333HWAcAWzfvt29rNPpdF/V5fn3SklJoX79+thsNnd8ANdeey2TJ08GwOFwkJSUBBin7H788UdmzJjh3talkKQgii36ZDKpf6yiTlo8te65x5iYcASSjnGu2eUsPLiQQc0HEeB16Zfa1RQdO3bkxRdfpG/fvnTt2pVnnnnmotfVtm1bJk2aRPv27UlISGDs2LEAjB8/nnnz5tG5c2eef/55vvzyywKXnz59OuPGjaNLly5s3bqVCRMm5Jrv4+PDlClTuPHGG+nevXuhDdb9+/cnOjqabt26MXv2bG655Rbi4+Pp2LEjn376qTspXapXXnmFpUuX0qlTJ+bMmUO9evUIDAzM155QmEcffRSn00nnzp0ZNWoU06ZNy3XkUpCPP/6YqKgounTpQocOHfj8c+N5IS+//DI2m40uXbrQsWNHXn75ZQDGjh1Lamoq7du3Z8KECURERLjX5e/vz4YNG+jUqRMrVqxwf99vvPEGvXr1onfv3rRr185d/qOPPmLlypV07tyZiIiIXKex/P39WbRoEf/9739ZsGBB8b/Egmitq/QrIiJCi/IxdmaUnnvVjXpv3/7aabMZEzdN1/qVIP32qn/rrtO76mNJxyo2yBKKjo6u6BBKxeHDh3XHjh0rOoxylZmZqW2u3+HatWt1165dtdZaX3bZZTo7O7sCI6t8CvqdA1G6gDpV2hREsew7k8L+1Rt5PPYgtf79b5TF9dM59AfnAmoz5/hyhrQYUiUbmEXVdOzYMW677TacTideXl588cUXAO5TVeLiSFIQxfLpigOMPPwnyt+fkJG3GhPtWbBvCd+16EJ2Vgz/6PyPig2yBmvWrBk7d+6s6DDKVevWrd3n7kXpkTYFcUGHYlOJXr2R3jFbCR05EnPO7fkHV5BuS2WW/RzXNrmWZsHNKjROIcSlK9ekoJQapJTaq5Q6oJQaX8D8Z5RS0Uqp7Uqp35VSTcszPlGwycv28szmWVhCw6j18EPnZ2yazrxadUhxZMhlqEJUE+WWFJRSZmAScAPQAbhDKZX3ototQKTWugswF3i3vOITBTsWl47f7Gk0TzpFgzdeP/90tbiD2Pb9xoyQUCLqRtClduE3/Aghqo7yPFLoCRzQWh/SWmcDs4BhngW01iu11umu0fVAo3KMTxTgh29/Y+S+FXjfeBOB13hc+71+Mr8FBnLakc79ne6vuACFEKWqPJNCQ+C4x3iMa1phHgB+LWiGUuohpVSUUioqNja2FEMUnmLOJNL5m4/JCgyh6SsvnZ+RkYBt67d8XacBrUJacVXD/B2jiUszceJE0tPTL1xQiFJWKRualVJ3A5HAewXN11pP0VpHaq0ja9euXb7B1SB/v/x/NEk5Q/jrr2MOCjo/Y9N0pvpZOOBM57Fuj2FSlfJnVKVJUhAVpTwvST0BeF7E3sg1LRel1HXAi0BfrXVWOcUmPDhS09jz9ge0W72QfT2uZdigaz1m2jiwaQqfh4YwsNlArm16beErqmp+HQ+nd5TuOut1hhveLrJIWloat912GzExMTgcDkaOHMnJkyfp378/4eHhrFy5kqVLl/LKK6+QlZVFy5YtmTp1KmvWrOGrr75izpw5gPHQm/fff59FixYVWD4gIKDQbqj/+OMPnnzyScDog2f16tUEBATw3HPP8euvv6KU4qWXXmLUqFGsWrWKV199lfDwcHbu3ElERAQzZ86Ux5pWE+W5i7cRaK2Uaq6U8gJuB3Ldj62Uugz4HzBUa122TwIRBcrcu5fdw4bD3Fn82eZKeryX+2En9p3zmeBjJ8Dqzwu9XqigKKuX3377jQYNGrBt2zZ27tzJU089RYMGDVi5ciUrV67k3LlzvPnmmyxfvpzNmzcTGRnJhx9+yHXXXcfff/9NWloaALNnz+b2228vtHyOgrqhfv/995k0aRJbt27lzz//xNfXl/nz57N161a2bdvG8uXLGTduHKdOnQKMvn0mTpxIdHQ0hw4dKpUnnInKodyOFLTWdqXUY8ASwAx8rbXepZR6HeN26wUYp4sCgDmuvY5jWmvpmL8caK1JnDuXU2+8SZLJh5lDnubVCfdQN8jHsxAzN77PDh9v3r38ZcJ8wgpfYVV0gT36stK5c2f+9a9/8e9//5shQ4bQp0+fXPPXr19PdHS0u1vs7OxsrrjiCiwWC4MGDWLhwoXceuut/PLLL7z77rv88ccfBZbPUVA31L179+aZZ57hrrvuYsSIETRq1Ig1a9Zwxx13YDabqVu3Ln379mXjxo0EBQXRs2dPGjUyrgPp1q0bR44cKfChO6LqKdc7mrXWi4HFeaZN8Bi+rjzjEQZnWhqnXnuN5AUL2VanNT/d+DCTHxtAmL9XrnJHoufyqTmdawLbMKjFxT9qUeTWpk0bNm/ezOLFi3nppZe49trcp+S01gwYMKDAXklvv/12Pv30U8LCwoiMjCQwMLDI8lBwN9Ljx4/nxhtvZPHixfTu3ZslS5YUGXPeLq8vpTtqUblIC2ENZ4+N5cg995C0cBEz21/P/FHP8cWT1+dLCE7tZELUe3ijeOm6T+T8cSk6efIkfn5+3H333YwbN47Nmzfn6n768ssv56+//uLAgQOA0Qaxb98+APr27cvmzZv54osvuP322y9YvjAHDx6kc+fO/Pvf/6ZHjx7s2bOHPn36MHv2bBwOB7GxsaxevZqePXuW1dcgKgnp+6gGyz56lGP/eJDMM2d5rdf9+PTpw/R7IvDzyv+z+P6Pl9lCBm/V6U3tILl9pDTt2LGDcePGYTKZsFqtTJ48mXXr1jFo0CB328K0adO44447yMoyrr1488033c8gHjJkCNOmTWP69OkA1K5du9DyhZk4cSIrV67EZDLRsWNHbrjhBry8vFi3bh1du3ZFKcW7775LvXr12LNnT9l/KaLCKO16KERVFRkZqaOioio6jCpF22yk/rmGUxMmkJGRxXMRY2h6VU8+ufMyvC3mfOWPn9jALUvvJ0Jb+ezutSgv3wqIumzs3r2b9u3bV3QYQpSpgn7nSqlNWuvIvGXlSKEG0VoT9/nnxE+fgSMxkbRadXnq8rFE9unGeyO7YjXnP5votGfx6tJHMAOvXP95tUoIQoj8JCnUEM7MTM59+ilxX35FQP/+/N7yCt48G8TIy1vw1s2dMJkKbiOYu/ABNphsvNJkKPUa9SrnqIUQ5U2SQjXnzMwk8YcfOPfFFzhizxF8883MvmYMH604wJ1XNOGtmzsV2mh8avv3fJi4hV7e4dzS761yjlwIUREkKVRjzqwsjt55F5nR0fj16kX4Bx8wLS2Mj5bsZWREI94cVkhCyEola/N0Xt3xKU4vC6/dMFWuNhKihpCkUE1prTn32WQyo6Np8N57BN80hC//PMR7S3Zzc7cGvH1Ll3ynjLQ9m21/vc3Pu79niZeJFG8LL3Z+mIahLSroUwghypskhWpG22yk/b2BuC+/JH39eoKHDSX4piFMX3uEN3/ZzY2d6/P+yK6YPRLCyZQTLFz/LguPr+CoGXx9rFxXrxfDuv6DXvWlHUGImkSSQhWXGR1NysqV4HCSsWsn6Rs2ojMyMIeFUffllwi94w6++/sYryzYxYAOdZl4ezcsZhN2p51lR5cxd9c3bIgzOoHrgZl/tLyVAT2fwt8roII/magoJ0+e5IknnmDu3LkVHYqoAJIUqjBHSgpH77kXp6tDNGuTJoQMH45fr14E9O+HycuLOVHHeeHHHfRvW5tP77yMbEc6s6J/YObeWZxMO0Uju4N/ptu4KfIxGvZ6DEz571MQ5UdrjdYak+niOhuw2+1YLJf2b92gQQNJCDWYJIUqyJmdTdrq1cR/MxNnWhrNfpiNT6dOqDwVybxNMTw3bztXtQrng1Ht+WLzh3wX/Q0pShORkcnzSSlc3ehqTCM/gqD6FfRpKo93NrzDnvjSvVu3XVg7/t3z30WWOXLkCAMHDqRXr15s2rSJ2267jUWLFpGVlcXw4cN57bXXAHjjjTeYOXMmtWvXpnHjxkRERPDss8/Sr18/unXr5u7Arl+/fjzzzDOkpqYSHh7OtGnTqF+/Ph9//DGff/45FouFDh06MGvWrAK7zI6Li2PIkCHs3LmTzMxMxo4dS1RUFBaLhQ8//JD+/fszbdo0FixYQHp6OgcPHmT48OG8+648Pbc6kKRQRTgzM0nfGEXK0qUkL1mCMznZOEX0wvP4dsn9fGStNd+sP8orC3ZxRcswbm3/ByPnPspZnc2AjEzuazqYzv4NoX43aD0A5MqiCrd//36mT59OcnIyc+fOZcOGDWitGTp0KKtXr8bX15d58+axbds2bDYb3bt3JyIiwr18dnY2UVFR2Gw2+vbty88//0zt2rWZPXs2L774Il9//TVvv/02hw8fxtvbm8TEROB8l9m9e/cmNTUVHx+fXHFNmjQJpRQ7duxgz549XH/99e5+lLZu3cqWLVvw9vambdu2PP744zRu3BhRtUlSqKS01qSuWEHqn3+SuX0Hmfv2gd2Oyc+PwAHXETRkCP5XXIHKc6rgUGwqL/64k3WH4ujVJgllfYOXD6TRMdvOB0Gd6dbvaWh6RSFbrdkutEdflpo2bcrll1/Os88+y9KlS7nssssASE1NZf/+/aSkpDBs2DB8fHzw8fHhpptuyrX8qFGjANi7dy87d+5kwIABADgcDurXN44Cu3Tpwl133cXNN9/MzTffDBTcZbanNWvW8PjjjwPQrl07mjZt6k4K1157LcHBwQB06NCBo0ePSlKoBiQpVCJaa2xHj5K+aRMpv68gdcUKTAEB+HTuRK377sMvMgK/Xr0w5dmbA0jLsvPln4eZtOoA3r7n6BW5iui0DdTKdvB6nSsZNugTTFbpoqKy8vf3B4zfwPPPP8/DDz+ca/7EiROLvXzHjh1Zt25dvjK//PILq1evZuHChbz11lvs2LGjwC6z8x4tFEa6z66epOvsUqQdDhJmzSL76NGC52tNwqxZJP/2G9rpBMCRnEz8t99y7P4H2H9VHw4OuoFTL75ExubNhD/xOG3+Xk/TqVOp869nCOjbN19CSMuy8/kfB+n7znLWrp5On6YfYGr0LodT1vNwYjKLWtzJ8Ju+lIRQRQwcOJCvv/6a1NRUAE6cOMHZs2fp3bs3CxcuJDMzk9TUVBYtWlTg8m3btiU2NtadFGw2G7t27cLpdHL8+HH69+/PO++8Q1JSEqmpqQV2me2pT58+fPvttwDs27ePY8eO0bZt2zL8BkRFkyOFUpT862+cfvU1TIGBtP5rDSav3M8kSPvzT06/ajQaerVsiSUsjIytW9E2G96tWxFw9dX4duuGX2QEXi1aFHkXcVKGje83HOOPP5bTz76E28J38kOgFadS3Jml+Ee9ftQa/CjU7Vimn1mUruuvv57du3e7n5QWEBDAzJkz6dGjB0OHDqVLly7UrVuXzp07u0/dePLy8mLu3Lk88cQTJCUlYbfbeeqpp2jTpg133303SUlJaK154oknCAkJ4eWXX87XZXbOIzcBHn30UcaOHUvnzp2xWCxMmzYt1xGCqH6k6+xSdPr110n4znjaVd2XXyLsrrtyzT/54oskzZtPvTdeJ+mnn9E2G37duxM0ZAi+nYpXee8+lcycdfs4sXUZ1+o/IWQbn4aGkGA2cWNYV57o8wYNQpqX+merrqpS19mpqakEBASQnp7O1VdfzZQpU+jevXtFhyWqAOk6u4LY4+LxatECS1gY5yZ/TuCAAVjr1AHg3Of/I2nefLzbtCF05EhCR44s3kqzUkk/sYOt27ZwdP8OmqRsY7xpLxv9LbwXFsZBrzC6h3fls17j6RTeqQw/nahoDz30ENHR0WRmZjJ69GhJCKJMSFIoRY74eMxhodR98QWO3HU3R++6m+DhN+NMSSV+6lQCBwyg7ksvFr0SezYkHCH7xDYyVn1IYOJu/NBcCVyOYmtYSx6r25V1tlgaBzRiYuSzXNPkGumwrgb47rvvKjoEUQNIUigFcVOnkfbnn2QdOYzfZd3xad+eJlP+x+n//IdzH38CJhM+HTpQb8LLWGrXzr1w6lnY9RPsWYSOOwDJp1A48QJOOusyyzIS/6bdadmpDcuz/mTegZ/wI5NnI5/ljnZ34GX2KigkIYS4KJIULpEzPZ2z77zjHjeHhQLgFxlJi/nzcWZmoqxWlNmj+whbJuz7FbZ+jz6wHKUdnPJqxgpbM3ZZOhLj64Oq7QuhZuzmFOIzFnA2+iwO7WBU21E80vURQn1Cy/ujCiFqAEkKlyg7JgaAuhNexhZzgqDrB+Sa776E1OmAY+tg53z0znmozEROW8P5xLsPy338SQ84Cdbj7uUsykIjeyPqetelc+3O1PGtw4g2I2gRLN1YCyHKjiSFS2RzJQXfjh0Ju/PO3DMddjj6F0T/jCN6Aab0WPZbffnapzWrg9qS4hcL6hBeJh+urt+TnvV70jy4Oc2CmtEgoAEWk/x5hBDlS25eK4QjNQ1nVlah8+3nznFuyhecees/YDZjbdLEtaAdDq5AL3gS+3ttSP9mGEui5/CwTy16NGzHLY1q80t4In7hVu5sfzv/G/A/1t65hknXTWJ0x9Fc3ehqmgQ1kYRQQyQmJvLZZ58BsGrVKoYMGVIm23n11Vd5//33C5x35ZVXAkbHfJ06dbqoWMaMGVOinlVztikqH6l5CnHiiSfIOnSIxp9PxqddOwAcSUkkzJpNelQUaevXg82GX48e1H1+HJYz67CvWkTm7l/YY8pgrY8/y4PDOOIbgFYai3LSLTySQS360bthbxoFNrpABKImyEkKjz766CWvy+FwYDaXvOvztWvXXvK2q8I2RfFIUihE5t69OOLiOHLnXQRedy3Y7aStW48jIQGvVi0Ju204Id1CyIj/i92bHmDzTsU6Hz+21w/CbgoCFPV9mnBb4ysY0Lwv3et0x2q2VvTHEkU4/Z//kLW7dLvO9m7fjnovvFDo/PHjx3Pw4EG6deuG1WrF39+fW2+9lZ07dxIREcHMmTNRSvH777/z7LPPYrfb6dGjB5MnT8bb25tmzZoxatQoli1bxnPPPcf48eMZPXo0CxcuxGazMWfOHNq5dmqio6Pp168fx44d46mnnuKJJ54AjLumc7rVKI7x48ezYMECLBYL119/fb4jkJdffpnjx4/z1Vdf8eGHH/LDDz/k6wbcc5vvvfdegWVExZCkUABnVhaOuDhCbh+FIyGR9A0bUVYLfp1b4XVZEDv1ZvZlTyXqhDdbvb3JqhcKGkIsTbm6Tg+GtOlDz/oRBHvn74ZACE9vv/02O3fuZOvWraxatYphw4axa9cuGjRoQO/evfnrr7+IjIxkzJgx/P7777Rp04Z7772XyZMn89RTTwFQq1YtNm/eDBgVdnh4OJs3b+azzz7j/fff58svvwRgz549rFy5kpSUFNq2bcvYsWOxWku2oxIXF8ePP/7Inj17UEq5u+DOMW7cOFJSUpg6dSrLli1j//79+boBv/rqq93lly5desEyonxJUshDa03q6tUA+DSvS2I/fzZ12cX+jMPs8DrGXi+jfyH8gwly1qFjaC8GtbqKG1pdSYhPSMUGLy5JUXv05aVnz57u7qu7devGkSNHCAwMpHnz5rRp0waA0aNHM2nSJHdSyOk2O8eIESMAiIiIYP78+e7pN954I97e3nh7e1OnTh3OnDmTr6vsCwkODsbHx4cHHniAIUOG5Gp3eOONN+jVqxdTpkwBjAq/oG7A8yaFC5UR5atck4JSahDwEWAGvtRav51nvjcwA4gA4oBRWusj5RGbIyWF+B/ncWrGdKwxp8nwhYfiPmYPFrCA1d+PWo46dPGJ5Mrmfbil05XU8Zd7BUTpupjuqHO6zc67jrzLX2xX1wMHDuTMmTNERkby5ZdfsmHDBn7//Xfmzp3Lp59+yooVKwDo0aMHmzZtIj4+nrCwsEK7AfdUnDKifJVbUlBKmYFJwAAgBtiolFqgtY72KPYAkKC1bqWUuh14BxiVf22lI/bQHjbO+hzHxiia7o/Daocj9WHpjSYOtNYEmRvR3+9Krmx9PTe2iyDQR3qHFKUrMDCQlJSUIsu0bduWI0eOcODAAVq1asU333xD3759yylCWLJkiXs4NTWV9PR0Bg8eTO/evWnR4vx9M4MGDWLgwIHceOONLF26lIEDB/Lyyy9z1113ERAQwIkTJ7BardRx9QcGFKuMKF/leaTQEzigtT4EoJSaBQwDPJPCMOBV1/Bc4FOllNJl0JXr7Lt60nFzCs01JATAlo5OzrX3JrDFZdzb8VYiug2+qCs5hCiJWrVq0bt3bzp16oSvry9169bNV8bHx4epU6cycuRId0PzI488UgHR4n4CXGZmJlprPvzww1zzR44cSUpKCkOHDmXx4sXceeed+boB96zwC+sqXJJCxSm3rrOVUrcCg7TW/3CN3wP00lo/5lFmp6tMjGv8oKvMuTzregh4CKBJkyYRRwt5qE1Rvn/ielRcPF59enHFVUOp36w7BNS+8IKiWqlKXWcLcbGqfdfZWuspwBQwnqdwMeu44+OlpRqTEEJUB+V5R/MJwPOp3o1c0woso5SyAMEYDc5CCCHKQXkmhY1Aa6VUc6WUF3A7sCBPmQXAaNfwrcCKsmhPEMKT/MREdVbS33e5JQWttR14DFgC7AZ+0FrvUkq9rpQa6ir2FVBLKXUAeAYYX17xiZrJx8eHuLg4SQyiWtJaExcXh09Ob83FIM9oFjWazWYjJiaGzMzMig5FiDLh4+NDo0aN8t29Xq0amoUoLVarlebNm1d0GEJUGtJ1thBCCDdJCkIIIdwkKQghhHCr8g3NSqlYoOS3NF+6cODcBUtVLImx9FSFOKtCjFA14qwJMTbVWufrxqHKJ4WKopSKKqjlvjKRGEtPVYizKsQIVSPOmhyjnD4SQgjhJklBCCGEmySFizelogMoBomx9FSFOKtCjFA14qyxMUqbghBCCDc5UhBCCOEmSUEIIYSbJAUXpdTXSqmzrqe/5UzrppRar5TaqpSKUkr1dE1XSqmPlVIHlFLblVLdPZYZrZTa73qNLmhbpRxjV6XUOqXUDqXUQqVUkMe8510x7lVKDfSYPsg17YBSqtR7olVKNVZKrVRKRSuldimlnnRND1NKLXN9N8uUUqGu6eX+fRYR40jXuFMpFZlnmXL9PouI8T2l1B7Xd/WjUiqkomK8QJxvuGLcqpRaqpRq4Jpeaf7eHvP/pZTSSqnwioqxqDiVUq8qpU64vsutSqnBHsuU7t9cay0vo13laqA7sNNj2lLgBtfwYGCVx/CvgAIuB/52TQ8DDrneQ13DoWUc40agr2v4fuAN13AHYBvgDTQHDgJm1+sg0ALwcpXpUMrfZX2gu2s4ENjniuddYLxr+njgnYr6PouIsT3QFlgFRHqUL/fvs4gYrwcsrunveHyPFfI3LyLOII8yTwCfV7a/t2u8MUaX/keB8IqK8QLf5avAswWUL/W/uRwpuGitVwPxeScDOXvewcBJ1/AwYIY2rAdClFL1gYHAMq11vNY6AVgGDCrjGNsAq13Dy4BbPGKcpbXO0lofBg4APV2vA1rrQ1rrbGCWq2yp0Vqf0lpvdg2nYDw/o6FrO9NdxaYDN3vEWq7fZ2Exaq13a633FrBIuX+fRcS4VBvPJwFYj/EUwwqJ8QJxJnsU88f4f8qJs1L8vV2z/ws85xFfhcRYjDgLUup/c0kKRXsKeE8pdRx4H3jeNb0hcNyjXIxrWmHTy9Iuzv+xR3L+kaeVIkalVDPgMuBvoK7W+pRr1mmgbmWINU+MhamsMd6PsUdb4TEWFKdS6i3X/89dwITKEKdnjEqpYcAJrfW2PMUq3XcJPOY6lfV1zqnXsohTkkLRxgJPa60bA09jPBmusrkfeFQptQnjcDO7guNxU0oFAPOAp/LsNaKNY98Kvx66qBgri8JiVEq9CNiBbysqNk8Fxam1ftH1//MtxpMXK5RnjBjf3QucT1aVRgHf5WSgJdANOAV8UFbblqRQtNHAfNfwHIxDMoATnN8jB+Pw/UQR08uM1nqP1vp6rXUE8D3GecQKj1EpZcX4UX+rtc75Ds+4DsFxvZ+tyFgLibEwlSpGpdQYYAhwlyvBVliMRcXp4VvOn9qsLN9lS4zz8NuUUkdc29uslKpXUTEWEida6zNaa4fW2gl8QVnWRaXVQFIdXkAzcjfi7gb6uYavBTa5hm8kdyPUBn2+EeowRgNUqGs4rIxjrON6NwEzgPtd4x3J3QB1CKPxyeIabs75BqiOpRyjcsUyMc/098jd0PxuRX2fhcXoMX8VuRuay/37LOJ7HAREA7XzTK+Qv3kRcbb2GH4cmFtZ/96uMkc439BcIf/jRXyX9T2Gn8ZoRyiTv3mpfJDq8MLYyz4F2DDOvz0AXAVscn2hfwMRHn+4SRh75TvyVB73YzT2HADuK4cYn8S4QmEf8Dauu9Rd5V90xbgX11VUrumDXeUPAi+WwXd5Fcapoe3AVtdrMFAL+B3YDyzP+WeqiO+ziBiHu77bLOAMsKSivs8iYjyAcb44Z9rnFfk3LyLOecBO1/SFGI3PlervnafMEc4nhYr6Hy/su/zGFcd2YAG5k0Sp/s2lmwshhBBu0qYghBDCTZKCEEIIN0kKQggh3CQpCCGEcJOkIIQQwk2SghBCCDdJCkIIIdz+H0xpYB8dzxnBAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "F={}\n", + "T={}\n", + "\n", + "beta_ts = np.array([2.33379720e-04, 1.27179106e-03, -6.69347820e-05,\n", + " 1.14647701e-04, 5.14366051e-12, 46/14*1.15151346e-03])\n", + "# the original tuning uses NO2 emissions unit; FaIR is in units of N\n", + "\n", + "beta_default = np.array([2.8249e-4, 1.0695e-4, -9.3604e-4, 99.7831e-4])\n", + "\n", + "for ozone_treatment in ['cmip6-old tuning/no feedback', 'stevenson', 'regression', 'thornhill-skeie']:\n", + " if ozone_treatment=='thornhill-skeie':\n", + " b_tro3 = beta_ts\n", + " else:\n", + " b_tro3 = beta_default\n", + " _,F[ozone_treatment],T[ozone_treatment] = fair_scm(\n", + " emissions = rcp85.Emissions.emissions,\n", + " tropO3_forcing=ozone_treatment,\n", + " diagnostics='AR6',\n", + " efficacy=np.ones(45),\n", + " b_tro3 = b_tro3,\n", + " )\n", + " \n", + " pl.plot(np.arange(1765, 2501), F[ozone_treatment][:,31]+F[ozone_treatment][:,32], label=ozone_treatment)\n", + "pl.ylabel('W m$^{-2}$')\n", + "pl.title('Total ozone forcing in RCP8.5')\n", + "pl.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Temperature change with different ozone treatments')" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEICAYAAABYoZ8gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/d0lEQVR4nO3dd3gVVfrA8e+56Z1AEiChBEKvAUIREFBEBFbBuvaCiliRXXWtq2tfZX+WRVEEFAQLiroioiBIUUBNKAkQCBDSSIdU0nPP748Z8BKTUEwyKe/nee6T6fPOmZn3nnumRGmtEUII0XLYrA5ACCFEw5LEL4QQLYwkfiGEaGEk8QshRAsjiV8IIVoYSfxCCNHCSOIXtVJKaaVUN6vjqAtKqceVUgtqGX+rUuqnP7H8D5RSz5vd5yul9juM66mU2qmUKlBKPaCU8lBKrVRK5SmlPjvXdQpxLpp84ldKFTp87EqpYof+G6yO71wopRKUUhdZHUdzo7V+UWt9B4BSKtT8UnOup3Vt1lr3dBj0CPCj1tpHa/0mcBXQFmijtb66PmKoiVJqnFIqpSHXWdfqu0Li+CXe0Or72IRmkPi11t4nPkAScKnDsGVWx1dVfe7MhlyHOGudgT1V+uO01hVnuyDZv6cnZXQaWutm8wESgIvMbhvwKHAIOAosB1qb40IBDdwGJAM5wExgKBAN5AJzHZZ7K/AzMBfIA/YB4x3G+wELgTTgCPA84FRl3tfMOJ4HwoD1Zn82sAxoZU7/IWAHioFCjJriOCCllm19BvgcWArkA3fUFlM15eYEPG6WVQEQBXQ0x2mzbA6Y5fIWoMxxNW6HQ4wPmWWaB3wKuDuMf8SML9WMWQPdzHFuwByML/MM4B3Ao4b4E4EhZvcN5nL6mv23A185lNNSszvJnK7Q/Jxn7qufzPXmAIeBSbUcb4OA7WaZfQp8Ajxvjju5z8wyqgRKzHV9DJQB5Wb/7eZ004FYc93fA50d1qWBe839cNgc9hdgp7lftgADTlf2gBfGsWV32PbgarbND1gCZJnl+yRgM8ftcpi30IxtnDnuMowvuFxgA9D7LI6HGrenSmybzHUeN9f/1xPlDfwDSMc4j2rMAeZyPjOnzTOXeeKYmWHumzJz+Ssd4n/YjP84xvnVFlhtHgM/AP4Oyx9hbkeuWWbjHMZtAJ7DyA0FwBogoJZjsxuw0Yw1G/j0T+VKq5N1XX44NRnOArYBHTCSyLvAx+a4ULNg38E4GS7GOCm/AoKAECATGGtOfytQAcwGXMwDLY/fv0i+NJfvZc7/K3BXlXnvB5wBD3MnTjDjCjQPuter246qSaSGbX3GPFCnmQe7R20xVVNuDwMxQE9AAQMxmiAwy+kboBXQCSMRXGKOO5Pt+BUIBlpjJLWZ5rhLME66voAnxpeWY+J/DfjanM8HWAm8VEP8S4C/m93zMU70ux3GzXYop6VVjgFnh+XcapbjnRhfhndjfCmpatbpipEQTxwTV5nz/iHxO5zodzj0n4zF7J8KHAR6YxwnTwJbHMZrYK1ZHh4YXzqZwHAz1lvM8nY7g7I/JbZayvR/ZtmHAnGYX1BVppuBURHyBXpgJMQJZpk8Ym6T6xnEVOv2VLPek8eKwzZVAP/GOB49qCUHmPNMN7fPDXgd2Okw7oMT+7LK8bwNI9mfyBHbzdjdMb7gnzanDcH4spmMcU5OMPsDHY6HQ2aZeZj9L9dybH4MPGEuyx0Y/adypRUJur4+nJoMYzm1Vt4e48R0dijYEIfxR4G/OvSvAB50SAinJADzAL7JPAhKcaiNAtdhtOeemDfpNHFPA3ZUtx01naj8MfFvchhXa0zVrH8/MLWWE2y0Q/9y4NGz2I4bHfpfAd4xuxfhkMgxvkS0+VdhJJAwh/HnYdZ0q1nv7cDXDvv9DuATsz8RGOxQTqdL/Acd+j3NadpVs84x1RwTWzj3xL8ah8SKcYIXYdb6zTgudBg/D3iumv049gzK/g/HU5XlOGHUdvs4DLsL2FBlutEYya+H2f8UsLzKNhzh918DtcVU6/bUcFxWTfxlnPoLosYcUM3yWpnL9DP7P6D6xH+DQ/8KYJ5D//38/uvyH8CHVeb/HrjF4Xh40mHcPcB3tRybSzAqNR1q2m9n82nO7WCdgS+VUnaHYZUYSfGEDIfu4mr6vR36j2hzD5gSMWounTFqN2lKqRPjbBhNSCc4dqOUagu8AZyPUeOwYfy8/zMc13EmMTnqiFH7qEm6Q3cRZrmc4XZUnTfY7A4GImuIPxAj6UY5xK8wElJ1NgJzlFLtzWmWA08rpUIxmix21rJtVZ2MV2tdZK7fu5rpgqn+mDhXnYE3lFL/cRimMGqOJ5ZbdR/fopS632GYK7+XL9Rc9qcTgHH8OG5PohmLEZhSHTHK+RatdZw5ONhxHq21XSmV7DhfLTGdyfacTpbWusShv8YcoJRKB14ArsY43k5ME4Dxa74mZ5ozOgNXK6UudRjvAvzo0F/teVWDRzCahn5VSuUA/9FaL6pl+lo158SfDEzXWv9cdYSZEM5WiFJKOZzonTCaIpIxatcBuuYLdbpK/4vmsP5a62NKqWkY1w9qmv44RiI8Eb8TxsFa0zrOJCZHyRjt9bvPYFpHp9uO2qRh/AQ/oaNDdzbGSdRXa33kdAvSWh9UShVh1Lg2aa3zzRN7BvCT1tpe3WxnGGdN0qj+mKjtC7Q2ycALuvYbEqru4xe01i+cw7pOt+3ZGDXjzsBec1gnjNo7SikPjGbR17XWqx3mSwX6n+hRxrdmxxPzncaf2Z4Tqm5XbTngJozmtYswavJ+GJWWEzWNP3t8JGPU+O88h3n/sG6tdTpGEyRKqdHAD0qpTVrrg+cSXJO/q6cW7wAvKKU6AyilApVSU//E8oKAB5RSLkqpqzHaYr/VWqdhXJj5j1LKVyllU0qFKaXG1rIsH4yLNnlKqRCMNnZHGUBXh/44wF0pNUUp5YLR/utW08LPIaYFwHNKqe7KMEAp1aaW+M90O2qzHLhNKdVbKeWJ0UxwIn478B7wmlIqCEApFaKUmljL8jYC95l/wfgp7dhfVRZGLa9rDeNPZytGm/KJY+IKYNg5LguM4/UxpVRfAKWUn3mc1eQ9YKZSari5z7zM48PnDNaVAbRRSvlVN1JrXYmxf15QSvmY59DfMK7DgNFMt09r/UqVWZcDU5RS483j9O8YFZAtZxDT2W5P1XOkOrXlAB8ztqMYlaoXz2H5tVkKXKqUmqiUclJKuZu30XY47ZzVHJtKqasd5s3B+HKorkJzRppz4n8Do0a+RilVgHFRZvifWN4vQHeM2tALwFVa66PmuJsxfpbuxdgpn2O0J9bkX8BgjJ+Uq4Avqox/CXhSKZWrlHpIa52H0Qa4AKP2dBzjDobanE1M/4dx0q7BuCtoIcYFp9M53XbUyKwpvonx0/cgxv4B42QEo430ILBNKZWPccdEz6rLcbAR42TeVEN/1fUXYezHn81yHnGmsZvzlwFXYFwXOIZxwf+Mt7+a5X2JcWHyE3N7dwOTapk+EqMGOBdj/x40YzmTde3DuFgYb257dc0p92McZ/EYdzp9hJHwAa4FLlenPkNzvtZ6P3Aj8F+M8+RSjNury84gprPdnmeAxWb819QwTW05YAlGs9QRjHNkW5V5FwJ9zOV/dbr4q9meZIxfFI9jJPJkjIrRaXNuDcfmUOAXpVShuU2ztNbxZxvXCUrrP/uLpvlTSt2KcWFutNWxNFdKqd4Yyc7tDJunhBDnqDnX+EUjp5S6XCnlppTyx6jtrpSkL0T9k8QvrHQXxu2AhzDutrjb2nCEaBmkqUcIIVoYqfELIUQL0yTu4w8ICNChoaFWhyGEEE1KVFRUtta66jM/TSPxh4aGEhkZefoJhRBCnKSUqvZpcmnqEUKIFqbeEr9SapFSKlMptdthWGul1Fql1AHzr399rV8IIUT16rPG/wHGq3cdPQqs01p3B9aZ/UIIIRpQvSV+rfUmjEfZHU0FFpvdizFe4yuEEKIBNXQbf1vzBWJgvJK0bU0TKqVmKKUilVKRWVlZDROdEEK0AJZd3DVfZVvj02Na6/la6witdURg4B/uRhJCCHGOGjrxZyjjn2Vg/s1s4PULIUSL19CJ/2uM/6WJ+fd/Dbx+IYRo1ErKK9l9JI8vtqfw0upYjpfW/XsL6+0BLqXUxxj/BzNAKZUCPA28DCxXSt2O8S7smt6jLYQQzZrdrknJKWZvWh570wrYn55PXEYhiUePYzcbwV2dbEwLD6F3e986XXe9JX6t9XU1jBpfX+sUQojGqKS8kv3pBcSm5bM3LZ/YtHxi0wooNGvzNgWhAV70aufDld286V12lA55GfikJ9HWZ0ydx9MkXtkghBBNRUFJOTFH8ohJyWNPqpHkD2UVnqzFe7k60bu9L1f1DSDcfoywgnT805OojD5A6aFDVGZnA8b/Vcz39KT1X6/BqUePOo1REr8QQpyjkvJK9qTmE52SS0xKHrtSconPPs6Jt92HtPKgdzsvrmiv6F2cSYdjR3BLiqfs1zjKEhM5MWGhhwdu3bvjPWYMbmFhuHXvhltYGM7t26NsdX8pVhK/EEKcgUq7Ji6jgB1JucQcyWVXch5xGQVUmFX5QB83woO9uSGwjD4FqbTLSIBdcZR+Fof9+HEAygE6dcK9Zw98p0zBrWcP3Hv2xKVjRyPBaw2FGZAZC8lfw/Z9MOFZcPer022RxC+EENUoKClnZ3IuUYk5RCXmsCMp92SbvJ+HC+HtvZjazUbfwlSCMxJR0fsoXb4fXVoKQKmXF269e+E3bdrJBO/WrRs2Ly9jBaWFRoLPWAe790DGHqO/JPf3IDxaw/CZkviFEKKuaW3cYROVmENk4jGiEnPZn56PXRsXXnu29eGmjoqhRVl0TI/Hec9+Sr/4PcmXe3nh3qcP/tddh3vfvrj37YtraGejFm+3Q26ikdh/+w4yYozuY4c5+Qyrmy8E9YG+l0NQbwjsZfz1CgSl6nx7JfELIVocrTVJx4rYeugo2+KPsi3+GOn5JQB4uzkzvJ0b13Uopk9OIm0S4yj/MZrK3FwA7F5e2GpM8pWQtR/StsJ3b0PaLsjYC2UF5poVtO4K7frDwOugbT9o2xdadfpDgrdrO6mFR2jv1R4nm1Odbr8kfiFEs3ci0Z9I8tvij5KWZyT6AC9XJvqXM8ollS4Z8bjGxlJ24IBRUwfsYWF4X3ghHuED8QgPxy0sDOXkZCT5owch9Rf4/h1I3Qnp0VBeZKzU1dtM8NdCu35Gkg/qDa5ef4gt9XgqB3IOcCj3kPHJO8ThvMMUVxSz6vJVdPLtVKflIYlfCNEsZeSXsPlANlsOZrMt/iipDol+km8xo3QyoUficPp5BxXmiyBtPj64DBiA70UX4REejseA/jj5me3reSmQ/Cus/QBSdxi1+XLjoi0untBuAAy+BYIHQXA4tOkGVWrqhWWFHMjcQdyxOOJy4jiQe4ADOQcoLC88OU2QRxBhrcK4svuVhLUKw9e1bh/eAkn8QohmoqS8kl8PH2PzgSw2xWWzP8NoXgnwcGaKz3FGeibTOWU/to07qczJAcApKAjPYcPwHDoUzyGDcQ0LM5psKkqNxL7nQ0j5FZJ/g4JUY0XO7kaSH3SDmeQHQUCPU5K8XdtJyk9iX84+4o7FcSDnAHE5caQeTz05jY+LD939uzOl6xR6+Pegh38PurbqWi+JvipJ/EKIJklrTVxGIZvisth0IItfDx+jtMKOq00xsVUZ97on0i15L04xO7Hn5wPgHBKC59ixeA6NwHPoUOM2SqWgIB2StsKa941En7YLKsuMFbXqBJ1HQsdh0GGo0Xzj5HIyDru2k5ifyN6je09+9h3bd7IW76Sc6OLXhYGBA7m659X08O9B91bdaefVzli3BSTxCyGajOKySn4+mM26fZms35dBRr5xV80gH82jbqkMyN2Pz+7tVKanA+ASEoLnxRPwGjYMz4gIXIKDjQXlpUDCz7DzVePvsUPGcCc3owY//C7oMMxI9j7tTq7fru0k5CWw5+ieU5J8UYXRru/m5EZP/55M6TqFvm360qt1L8JaheHq5NpwhXQGJPELIRq1jPwS1sVmsi42g58PZVNSbqeVM1zjmsXosnjaHozGHrcfAJuvL54jRuA18y68Ro40avRg3E6Z8CP8+jMk/GT0A7j5QefzYMitRq2+3QBw/j1J55XmEZ2ymejsaKKzoonJiqGg3GhCcndyp0frHlwWdhl92vShT5s+dG3VFRebC42dJH4hRKOitWZPaj4/xGawLjaTmCN5APRzLeVxexID0vbgHh2FLi5GubjgPngwXg8+iNeokbj36WPccVOYCfEbIOolSNgMecnGwj38ofMo46Go0FHGnTZm23yFvYJDuYfYlbWLXVm7iM6KJiE/AQCbstG9VXcu6XIJ/QP60y+gH138uuBsa5optGlGLYRoVux2zc6UXFbHpLF6dzopOcU4YecvLjncU3KI0IM74WAcAM7B7fG5fBpeY8bgNWwYNk9PKC+GxC2w7lM4tMF4SAqMRB96Pox8wEj0gb3BfPdNSUUJMZnbicyIJCojiuisaIorigFo7d6aAYEDmNptKgMCBtAvoB+eLp4WlEz9kMQvhLBEpV0TlZjD6t1pfLc7nbS8Erx0OTfYUhmXFUvrPVHo3BxwcsJjUDjef/8b3mPH4ta9u9F8kx4N29+FQz9C0jaoLAUnV+g4HMb/E7peAO0HnqzRF5YVsjNtC1EZUURlRBGTHUOFvQKFomfrnkzrNo2BgQMZEDiADt4dLLvw2hAk8QshGkylXfPL4aOsjknnuz3pZBWU4m8v5SaSGZUWg19MFLqkBCc/P7zGjsF77Fi8R4827qUvLYTDG2Hlf+HAWihIMxYa1AeG3gFhFxjt9OYDUoVlhUQe2cyv6b8SlRHFvmP7sGs7zsqZPgF9uKnPTUS0jSA8KLxBbqFsTCTxCyHqld2uiUzMYeWuVFbvTiO7sIy2lce51Z7I8JQYvPbsgIoKnIOC8LniCnwunoBnRATK2dl4n82+TyDuO+OibGWZ8V6bsAug+0ToNv7kXTflleVEZ0ezLW0b21K3EZMdQ6WuxM3JjQGBA5gxYAZD2g5hQMCAZtVscy4k8Qsh6tyJNvtvdqXxbUwa6fklBFfkc0f5YYYn7cR9327QGpdOnfC99RZ8JkzAvX9/ownnSCSs/xfEfQ/Zxt06tOkOw2ZAj4nQcQQ4u6K15kDuAbbtWcO2tG1EZkRSXFGMTdno16Yf0/tNZ0T7EQwMGoibk5uVxdHoSOIXQtSZbfFHWf5bMpsOZJFdWEa7snxuLotnRPIOPPbvAcCtRw987r0XnwkTcOvRHWWvgMSfYfXDEPsNFKaDzcW4GBtxG3S/GNqEAZBfls/WlB/ZnLKZn478xNGSowCE+oYyNWwqI4JHMLTd0BbXdHO2JPELIf6UjPwSohJzWB6ZzMZ9GYSXZnJPeQoDk2PwjPs92fvOegCfiZfg1rULlJdA/I/wv9dh/7dQnGO876bbRdD7MuhxMbj7nazVb45ZyOYjm9mZuZNKXYmvqy+jgkcxMmQkI9qPoJ1Xu9qDFKeQxC+EOCfb4o/y0up97ErOBa0Zl3eQ5QfX4J1yGAC3Xr1OTfZlRXBgDXz2rPG3rNB4gKrnJOh9KYRdCK6elFSUsDV1K5uPbGbzkc2kHzeewu3VuhfT+03n/A7n0z+gf5O9h74xkJITQpyxw9nHWRWdytq9GexJOsoVx3bz94KDBCTFoXKO4dKhA22eexbvsWNxCQqCsuNGkl/+jPG3vAg8A6D/VUayDx0Dzq7kleaxKXkd65LWsSV1C8UVxXg6e3Je8HnMHDCT0SGjaevV1urNbzYk8QshalVp12xPymHptkS+3pWKT8lxplWm8M+dq/DMTDXehzNmNF4jR+I7aRKqNAfi18OG1b8ne69A4x+P9J1mPDlrcyKzKJP1B79gXdI6ItMjqdAVBHkEcVnYZVzY6UKGth2Ki1Pjf/1BUySJXwhRrezCUt744QDfxqRxtLCUMVmxfJi2hTYJ+0FrXDt3JuideXiPHWs87JQeA8uvh4NrjQVUk+yT8pNYu+cD1ietJzo7GjAuzN7c92bGdxpPv4B+2JTNuo1uISTxCyFOkVdczv92HmHO9/vpmhrHP4ri6HU4Gpejmbh07oTffffiPWoU7n37ohI2wqq/Q8pvxpO07n4w7jHocYnxwjObjfjceNbGLGBt4lr25xi3Z/Zt05cHBj3A+E7j6dqqq8Vb3PJI4hdCoLVmQ1wWS7Yk8NPBbPzzj/Kv+G/pdWgnysPDaMaZeDG+kyejdDnEb4SPnoLDm4x/MdhhqPGahIjpVLr5svfoXjZFz2NtwloO5R1CoRgUNIhHhj7CRZ0uor13e6s3uUWTxC9EC7f7SB7PrtxLzIFUpmbHsCAnlsBDe7C5uhLw97/R+qabsLm7Q1o0fHkn7F8NFcXg3gqm/B8MuhGc3ThafJT50e+wMn4lBWUF2JSNIW2H8Ndef2V8p/EEeQZZvanCJIlfiBYo+VgRyyOT+W53OglpOdyQ9DP/jFuPU3ERrqGh+Nx2K/7XX4+L63H47S3Y/53xn6ncfI1E32sydB7NgYJENsUuJSY7hp+P/Ey5vZxJXSYxtsNYhrUfRmv31lZvqqiGJYlfKTUbuAPQQAxwm9a6xIpYhGhJyivtvLY2jnc3HKBfdjw3HD9IeMIO3HKO4n3ReAJmzMBjwADjgar1z8NvCwFt/Feq8U9DxHTs7r7EHotl2dZnWBm/EoDOvp2Z0nUKt/S9hS5+XazdSHFaDZ74lVIhwANAH611sVJqOXAt8EFDxyJES1FWYeernUd4Z8MhWu+J5KO4VfgeTUe5ueE1YgStb70Fr2FDjTtyvnoH9q2C0nzjXxCOmgW+waQWprI4eh5rE9eSVZyFs82ZO/vfyXW9riPQM9DqTRRnwaqmHmfAQylVDngCqaeZXghxjn46kM0/VkRjS05gRspPROzfilv3brR59BV8xo832u+jP4E3Z0BeknFnTo9L4Lz7yPBty8aUjWz6bRM/H/kZFFzQ8QLGdRzH6JDR0pTTRDV44tdaH1FKzQGSgGJgjdZ6TdXplFIzgBkAnTp1atgghWgGco6X8doPcXy7bhd/j/uG8ISdKHd3/KdPJ3DWA9gKkmD7fNj9uXEPfsgQmPgC9JxEkb2c92Le44MfPqDCXkGIdwg39L6BG/vcKO/FaQasaOrxB6YCXYBc4DOl1I1a66WO02mt5wPzASIiInRDxylEU1VSXsnbGw7x/uZDnHdwGwv3foOb0rS55278b7oJ56LD8NHlxhsxwbjf/vJ3sfe7il8yfmP9b6/wQ9IPZBdnc1nYZUzvN52ufl2b9X+kammsaOq5CDistc4CUEp9AYwEltY6lxDitDLyS7hrSSSev/7E2wk/EpCZhGdEBO1fehFXfzfY+jpsnWs8VTvhWeh3Jdo3hFWHV/HO19NIzE/Ew9mD89qfxy19b2Fw28FWb5KoB1Yk/iRghFLKE6OpZzwQaUEcQjQrUYk5zFq8lRt+/phxiZG4hoYS8Mq/8R0/CrXhRdj+ofEfrPpfBZNeIb48lw3Ja/ju8HfEHoulV+tevDLmFS7oeAHuzu5Wb46oR1a08f+ilPoc2A5UADswm3SEEGcv53gZ/1m7n+82RPNs5Id0yU4k4P77CJg5E3VoHbwzCgozjfvvRz5AtC7i1Q2z2Jm1E4Ce/j15duSzTO02Vd6T00JYcleP1vpp4Gkr1i1Ec7LlYDb3L/2N8/dsZOGBtbjZIOStufiEucOHUyFhMwT2hus+Js23HUv2LuHjfR8T6BnIQxEPMTF0olysbYHkyV0hmqjPo1J46/01vBq1lPY5aXiNHk27xx/F9cAi+GAueAXBxS+Q3W8ar0XPY1X8KgAuDbuUR4Y+go+rj8VbIKwiiV+IJuZ4aQXPr4oldtU6/i9yCZ7engS//RbewweiPr8NDm+EYTPQ45/h2yMbeOnbaykuL+a6XtdxU5+bCPYOtnoThMUk8QvRhEQmHOOhT7Yz/NdveWHfGty6dqHzvHm4HN0Eb8+Akjz0ZXNZ27ot87+/hf05+xkQOIDnRj1HVz95/bEwSOIXoolY+NNh3v90I4/uXE5Y1mF8Jk6k/UN34vTdzXAkCkIiKLrkRZ6JX8HqmFcI9Q3lhdEvMKXLFJxsTlaHLxoRSfxCNAGLtySw6Z1lvL3zM1w93Wk/Zw6+AwJQH18G9gqY9g7JXUby4Ma/cSDnALMGz+K2vrdJwhfVknu3hGjEtNa8+v0+1s77iEeiPsJnYH/C/vcVfoEpqMWXgpsPRbd+wzxyuXzllaQdT+Pti97mjv53SNIXNZIavxCNVFmFnUc+30XlF8t5LOZ/eIaH0/mVJ7F9f6dxm2a3Cewc8wAPbZ5NRlEGF3e+mIeHPiy3Z4rTksQvRCNUVmFnxge/0v3L97ni0Ca8x19IyO3jsC0eDyi47L9826oNT22cTVuvtiy+ZLG8XkGcMUn8QjQydrvmoU93MOyTN7kgZQf+N91E26l9UStug5DBFE6dy/8dXM5nm19lSNshvD7udVq5t7I6bNGESOIXohGptGue+iqGTh++xQUpOwicPZuAwU6w4jYIHsTWCY/z9MYHyCjK4JY+tzBr8CxcnFysDls0MZL4hWgkisoqeOCjHXRc8T5TErbS+o7bCQhNhFVzoftE3u45knkbZhHqG8riSxYTHhRudciiiZLEL0QjkFVQyl2LtnDhNwsZnxxFq2uuJqjzXti6Ah1xB28Hh/JOzHwuC7uMp0Y8JW/PFH+KJH4hLJZVUMqNr//A7Wvm0T/rEIEPPECbDvtQ21dQfuETPGPL4+uY+UzrNo1/jfyXvEFT/GmS+IWwUGlFJbPf28jfVs6hY1E2wf9+GT+PSNjyAbnn3cPsolgiMyK5J/weZg6YKf8FS9QJSfxCWKSi0s6Dy6KY/OVcOhYdpfPb/8UrcwnsWEHe4Ju5rWQfifmJvHz+y0zpOsXqcEUzIolfCAtU2jV//3QHfZa+yaCsA7R/+jG8DrwIyb9QdOET3FOwk8Rjibx90duMaD/C6nBFMyONhUI0MK01T3wRTaclc7koOYqAmbfTKut1SNtF3uXzuKdoL3uO7mXO2DmS9EW9kBq/EA1swaZ4/BbNZXLCNtrcMZ0A7zWQnkrCVe9yX+wCUgtTefn8l7mw04VWhyqaKUn8QjSg6JRcds9byMz4n/C/6UYCQ3ag9v/CbxOf5sEdr+Jsc2bRxEVyj76oV9LUI0QDKSgp599zV3J79ErcRo2ibfcDqP2riBr7IHcdXEobjzYsm7xMkr6od1LjF6IBlFXYmfX+Vm5ZtwAnP186ja1A7f2S5HEPMzt9DSHeIXw46UP83PysDlW0AFLjF6Keaa158vOdjPnkdToUZhF66zCc47/kyOgHuT17E3bsvHnhm5L0RYORxC9EPVu4OZ6gD+YyPCOW9vfeiFf6YtL6XMrtedsoLC/kvQnv0cWvi9VhihZEEr8Q9ehARgG731povHTt5uvxL1tGtn8HptuyyC8t4L0J79G7TW+rwxQtjCR+IepJpV3z6nvfMz3ma1xHjSIocAP245k83qUPWSVHeWfCO/QN6Gt1mKIFksQvRD15Y00sE1YvwubhSecR+agjUSw67ya2HtvDo8MeZUDgAKtDFC2UJH4h6sGmuCyy332XfkcP0+GKPjgf+Z5t59/L3NT1TAqdxJXdr7Q6RNGCSeIXoo5l5Jfw3hvLuSF2Ld7jR+JX9hV7+09jVvoPdPHrwj/P+6e8ZVNYypLEr5RqpZT6XCm1TykVq5Q6z4o4hKhrdrvmn+9v5O6fFmPr1JGQbrvJ8A/hvopEWrm14t0J7+Lt6m11mKKFs6rG/wbwnda6FzAQiLUoDiHq1JIthxm9Yh5+uowuV4Vgz9vH7I6hFFWWMHf8XII8g6wOUYiGT/xKKT9gDLAQQGtdprXObeg4hKhrBzML2fHWIiIy99P+xgtxT13BkoFT2F2YxDMjn6G7f3erQxQCsKbG3wXIAt5XSu1QSi1QSnlVnUgpNUMpFamUiszKymr4KIU4C2UVdl6Y/z23Rq/EdcgA/Ms+Zn/nYbxVuJ/xncYzsfNEq0MU4iQrEr8zMBiYp7UeBBwHHq06kdZ6vtY6QmsdERgY2NAxCnFWXvg6mstWzcfF3ZVOAxMpdPfibz42Wrm14qkRT8nFXNGoWJH4U4AUrfUvZv/nGF8EQjRJX+5IwX3B2/TKSaLT1X1xLorlqZ5DOVKUwZxxc2jj0cbqEIU4RYMnfq11OpCslOppDhoP7G3oOISoC7Fp+ax6fQnT4jfTaup4fEu+Ymnf8aw7tpvZQ2YzKGiQ1SEK8QdWvZb5fmCZUsoViAdusygOIc5ZSXklz8/9hocjP8W5fz/aBf5InHsor5UkMK7DOG7uc7PVIQpRLUsSv9Z6JxBhxbqFqCuLNsRx85p3cfbxIvSCYsqOZvJolwh8K4r516h/Sbu+aLTkyV0hzkFqbjHJ8xfRqSCT0OsG4JK5mTcHTeFAYQrPjnqW1u6trQ5RiBpJ4hfiLGmteWnJRq7euwbXwb3wzvmYdf0mseRoFNf2vJYxHcZYHaIQtZLEL8RZWrg5nsFfvIeb0nTssp34kP48UZpAvzb9eGjoQ1aHJ8RpSeIX4izEpOSx560FjEjfS7tLQil1P84sfy/cnN147YLXcHNyszpEIU5LEr8QZ6ikvJJX3vqa22K+wS1iAP5eG3my+2CSizKYM3YO7bzaWR2iEGfEqts5hWhyFq7bx41rF+Dk50un/nH86BXC+tI0Zg+ZzdB2Q60OT4gzJjV+Ic5Aam4xqQvep1NhJp2ntqesNJGXAwLo7t9d7tcXTY4kfiFOQ2vNK0s3cWXsWtwGdcO7ZC1z+19EemkOTwx/Ameb/HAWTYskfiFO45voNPqtWICbTdOh83Y2dongw/xY/trzrwxpO8Tq8IQ4a5L4hahFZn4JK/+7lFFpu2k7ri1H/TRPuBTTu3VvHh76sNXhCXFOJPELUQOtNc8u3sT0X5fj1KUDrfx/5ZEuvajAzpyxc+TWTdFkSeOkEDX4ZlcqQ1a8i4+9jM7D0ljRvgs7S7N46fyX6OTbyerwhDhnUuMXohol5ZWsfedjRqTvpe3kbhQ5HeFNHzeGtRvGlC5TrA5PiD9FEr8Q1Viwbh+X//oFukM7Wruv47keERTZy3hs2GPy1k3R5EniF6KKxKPHSXtvIe2Kcug8IJVVId1YW5rOfeH30c2/m9XhCfGnSeIXwoHWmv8s/pFr9v2AR//2FPhn85K3M4ODBnNr31utDk+IOiEXd4Vw8P3uNIb8byFOzjaCw3ZzX9e+VFLK86Ofx8nmZHV4QtQJqfELYSouq+Sbtz5mWMY+2o1pw5q2XmytyGH2kNl09OlodXhC1BlJ/EKYFq2L5epfPoeQIJzbbOfVgAD6B/Tn6h5XWx2aEHVKEr8QGE/opi9cRFBxLp0HpvJ6x67k6XL+ed4/pYlHNDuS+IUA3lqxjctj1+HRJ4DYoAJWOJdxU5+b6NW6l9WhCVHn5OKuaPF2H8nD/5P3caOSoC57mNWxF8FuPtw98G6rQxOiXtRa41dK1XhFSyn1l7oPR4iGpbVmwcLVTEj8jdbh7iwNa0t8RQFPjHgCTxdPq8MTol6crqlnrVIqtOpApdR04I16iUiIBvT97nRGfv8heLlTFpbAu95uXNz5YsZ0GGN1aELUm9Ml/r8Ba5RS3U8MUEo9BswGxtZnYELUt7IKO6vnfUz/o/G0H1TEa6GdcHJy5R/D/mF1aELUq1rb+LXW3yqlSoHVSqlpwB3AMGCM1jqnAeITot58+ethpm5bgWrrQ3zXBNY5uXBfv7sI8gyyOjQh6tVp7+rRWq8DbgM2AF2BCyXpi6auvNLOrgXLaF90jOA+Kczp2JUgzyBu7iv/P1c0f7XW+JVSBYAGFOAGjAcylfF6Qq219j3XFSulnIBI4IjWWi4Uiwb1dWQiF+/4Dud27nzXo4Q99uO8NORJPJw9rA5NiHp3uqYen3pc9ywgFjjnLw8hzkWlXbN9wUfcWJyD1/BjvB7QgWFB4fKefdFiWPIAl1KqAzAFWGDF+kXLtjIqifFR3+Ic5Mrcga0pxs4TI56Q9+yLFsOqJ3dfBx4B7DVNoJSaoZSKVEpFZmVlNVhgonkrrahk67tLaV90jLJ+Waz0dObWvrfS1a+r1aEJ0WAaPPGbD35laq2japtOaz1fax2htY4IDAxsoOhEc7dsw36mRK3Epa0r/x3oRytXP27vd7vVYQnRoKyo8Y8CLlNKJQCfABcqpZZaEIdoYXKLykh4dxGBxXkUhGfxs7sz0/vfjrert9WhCdGgGjzxa60f01p30FqHAtcC67XWNzZ0HKLleX/VdqbG/oBHVzfe6OdPgHsbru11rdVhCdHg5O2cokXILCihfMkiPCrLSAnPIMrVxoyBd8ntm6JFsvTtnFrrDRgPhglRr95fsZVJh7bg013zj+5BhPp24qoeV1kdlhCWkBq/aPaSjxXh8ckHOCvNzxGFJNjsPDz0YVxsLlaHJoQlJPGLZm/hsvWMT/gNz56lvNm5NaOCR3F+yPlWhyWEZSTxi2YtJiWPjis+QLna+Hx4BUUKHh76sDysJVo0Sfyi2dJas3T+V4xM34Nr3+Msa+fLld2vIqxVmNWhCWEpSfyi2dqwL4Pz1iwDbxcWjXDC2cmVmQNnWh2WEJaTxC+apUq75vu5y+iZmwwDjrKylQc39rmJQE95ClwISfyiWfpiWzwTt67AuY0zrw/3xdfVh9v63WZ1WEI0CpL4RbNTUl5JzNsLaVeUQ0ZEDlvcnbkn/F58XeUN4EKAxQ9wCVEfPt+4j7/s+g63DornwlvTvVU3/trzr1aHJUSjITV+0ayUVlSSuOB9fMqL+Wl4Eak2zWPDH8PZJnUcIU6Qs0E0K59v3s/Fe9bh1LGCN3v4Min0Eoa2G2p1WEI0KlLjF81GaUUlCfPfx7u8hE9HaWxOrvwt4m9WhyVEoyOJXzQbH63by4Td66joWMHyLh7cOeAu2nm1szosIRodSfyiWcgrKifjnXfxKS/mw/MVAW6tuLGP/JsHIaojiV80C+9/uZXJ+zdQ0q2c1Z3duCv8XnnXvhA1kMQvmrzU3GJcF8/HWdl5+wJnQjzbcmX3K60OS4hGSxK/aPKWLvuBC5KiONa3mG0BLtw7eBYuTvKufSFqIolfNGlJR4sI/Hwx2lXxwjhPBrbpx5SuU6wOS4hGTRK/aNI+WrKaEWl7iAsv5YiXjcdGPIlNyWEtRG3kDBFNVnxWIR2+XIJ2g3+P8mBa2FT6BvS1OiwhGj1J/KLJ+mjB1wzJjGPb0AoqPd24f/Asq0MSokmQxC+apK0Hsuj79QdUesBbI9y4Y8AMede+EGdIEr9ocioq7Xz7f4volZPM/8bYae3lz819b7E6LCGaDEn8oslZvnEfk7d8xvG2dj4d5MKDw/6Bu7O71WEJ0WRI4hdNSn5JOSlz36ZV6XH+c7ETQwP6M1lu3xTirMhrmUWT8uEXW5m070cO9KpkfwcPVpz/Ikopq8MSokmRGr9oMnKLynBePB+bsjPnIlfu6DedLn5drA5LiCZHEr9oMj75cgujE6PYHK7xDQjg9vCZVockRJPU4IlfKdVRKfWjUmqvUmqPUkpuvhanlXO8DJYuRDtplo505qGRT+Pm5GZ1WEI0SVbU+CuAv2ut+wAjgHuVUn0siEM0IW+8v5bzE6L4MVzRq0NPxnYcZ3VIQjRZDZ74tdZpWuvtZncBEAuENHQcoun4ISaVgR+9QYWrnU9HOTF79LNyQVeIP8HSNn6lVCgwCPilmnEzlFKRSqnIrKysBo9NNA6FpRVsfOVteh1L5v0JNkZ0GUrfNvI+HiH+DMsSv1LKG1gBPKi1zq86Xms9X2sdobWOCAyUR/FbqvdWbOWq7f8jPtTOln4uzBr1rNUhCdHkWZL4lVIuGEl/mdb6CytiEI3fgYwCfBb9Fycq+M8kFx4b8TgdfTtaHZYQTZ4Vd/UoYCEQq7X+v4Zev2gatNYs+u/njEqJ4avzFOHdhjCt5zVWhyVEs2BFjX8UcBNwoVJqp/mZbEEcohFbuT2ZCd8vIMdXs+48V5648P/kgq4QdaTBX9mgtf4JkDNY1CjpaBG/vPoWNxYc499XOTH7vEdo49HG6rCEaDbkyV3RqBSVVfDCq59yXfRKdnQD9/AuXNr7OqvDEqJZkZe0iUZDa80z72/ktu/mUuBj570pLiy9eJ408QhRxyTxi0Zj0Y9xjF78Ml6VRfzjKmdmjnyAELmLR4g6J009olE4mFlI7v+9QvdjKbzxFxvDBl3ANQPusDosIZolSfzCcpV2zbzXPmXywS18G6FwGdaTJy98TZp4hKgn0tQjLLdw7R6uXPM26f6wfUJr3vvLUpxtcmgKUV+kxi8s9b+dR6h47UVaHS/i/clO/GfaMjxdPK0OS4hmTapVwhJaa+au3c/x1+cwNWEHX56nmHbJDNr5dbY6NCGaPanxC0u8+kUUbZ5/mKnxP/PtUEXW5DAuj7jf6rCEaBGkxi8a3Dc7k+nx2mN0P5bCu5NsFI/qwpzLPpKLuUI0EKnxiwa15VA2cU/+k57ZKcybZKPjpZfw1hVf4ePqY3VoQrQYUuMXDSLx6HFeXxVD688/4KqD21gVoQifdjV3nP+M1aEJ0eJI4hf1bsuhbJb+ax63Rn+FT0kRP/ZXZF41gL+Pftrq0IRokSTxi3qVV1zO5sef4/5da0gK1sy5wJnWfbvx5qVLpE1fCItI4hf16pOP1/GX6DX83Fvx1RWtuD1iFlN7XYNNyeUlIawiiV/Um+KySoI/fZliV9h9RWdWXP8Vbk5uVoclRIsn1S5Rb775/je6JmWwIVzx5BWLJekL0UhI4hd1amtcJq98+gvllXbsH7yEXYHrxEEEeAVZHZoQwiRNPaLOlJRXEnffXVyatI+Xk55l6v44tvRWXDX5RatDE0I4kBq/qDNfLPuOYUn7ALj0w3/hUgHHxnUiuJW8f0eIxkRq/OJPOZyaw5fzPqX18GF0ee9p8jzBpRI8SyuJDlVMu/JZq0MUQlQhiV+cs5LyCnbMuI1JB/fDZ2BHs+ByF3ofquT8aE3cSD/+GjzM6jCFEFVIU484Z18+/ya9D+4nxwvKneDDi2xcc8crZN46igdnODH+usetDlEIUQ2p8Ysa2e12kncfwNXbi/ZdOwDGe/STU7LZ8PpCBn23mL0dYePdvTl4PIX7wu/hvLBJDO86kSPnH6Gjj/yjdCEaI0n84qSSwiJivt+EvcJOaVERFR+8SfuMLAoVJL0wh8GXTmTVTXfSc+c2hgI5XrDrms68cfmKU16/YFM2SfpCNGKS+Fuw4+kZHPwpkuzsfEr276Xthq/wLq44Of6oDyy50MaN6+2kf/YK367+lp47txHbAVZH2Cjr04Y3r/5Y3rkjRBMjib+F0JWVpG6NJG33PnSlJnfzKoJ37cZVQ7A5za5QReRgG3YPJyorNaE9e/LIlNf54doJ9NuRCaznh3DF8ZtH8tSIR+jUuoeVmySEOEeS+JsprTXx33zH4f8tQ5cU47c/Dp+CCrzM8c4u8P0QRWYPV1q5u+HVxp8BAy/l1cF3o2ynXvM/MiKIfvGZxAVD25l/5Zpx8jplIZoySxK/UuoS4A3ACVigtX7ZijgaG223U5GWhtYaXVaOLi/H5uWJa4cOp0xXmZtL5saNZB2KxTMomLDrb0TZbOTv2UvCunX0nXEnP866jZCNO2ntCqXOENdecWisE84dA6hEE9gxlJvHPkFAm9PX2q+667886XsLfYJ78YgkfSGaPKW1btgVKuUExAETgBTgN+A6rfXemuaJiIjQkZGRZ72u5EO7ifpiLhW//UpwfDF2G5S4w3FfF0r9PHALj2Dqw2/VOL+9tNRIqPH78Q3tjnI795eMlSYlErd0HnmH91OZcwxbbiHORWV0f/cDyhOTOLjkLbwOpOJR/Mf9kX/fDQy/70kAYha8if2NebiW/z4+adIgPJx8CPxmEwAlbjbcS+1sDlf4Xn8Jfn6tGRh2EcEdRpxz/EKIpkcpFaW1jqg63Ioa/zDgoNY6HkAp9QkwFagx8Z+rmLuuoWeKptQZDndxxu5iw7WoAt9j5bQ7WE5m/Hp42Ji2ICOV9c/NwDkxFXtwED57kwjI1jiZefiH4a25cvHP1a5H2+389uqjuHfqwoDr7j5lXNam9eyd8xRtDh7D2Q6erpDrBXnu0P0YZF19IwCuHhDVVZHV3gXtbEPbbGBTTFhXxLGVn3J0zBRiH7yTNqnHORACe8e1xaN9R7oti6Tr6h0AbOkHXoWKgQl2tvZTTPj3B3TpLA9QCSFOZUXiDwGSHfpTgOFVJ1JKzQBmAHTq1OmcVmS/fAqHgFHXPEh4YMgp4z6/MpzgxFIADv26joz77qNHvjnyQCKJ7WDHQFcqXZ3ovbsYW25BtevQdjvrrxlH8O4sitwBh8S/9blZ+C5bg5sn/DjYRsnIwYT0OZ8Ogd0Jb9+DFe/ejWfUQbI7e9L5khu5fsydeLl5nbL85Vf0p//eCg7ddgNOWrPyAieG3PgYj466AYB5aZNx/e4wuwY6c/WTH3M08zAr177O5Vc9L0lfCFGtRntxV2s9H5gPRlPPuSzjL/e+WuM4u6szbuWllBcXkTz7frzKYPvM86koLsLJxYUr73sLD3dPAH4Y0xuncnu1y1n/6K1G0ncDzxLYfM14+v37bfISDuLz0Rp2h4H9vpncO/F+bFUums5+7OvTbsPxYH/Ym0Wl1kTeMpBHZn1yyviZT3xD5NSvmdV7Ii4uHnQO7sfg8EtPu1whRMtlReI/Ajg+3dPBHNag7G4uuJXB149dS5+jmqjrw7nxwfnVTlvmqnAptaO1xp6fj5OfHwB5CYdo8+1vxIQpPGbOpNvD8wiITmX93H/gvTcBPw/wf+hxLr7gpnOOc8CMp/imbDZe4RHcf88HfxivbDaGDph2zssXQrQ8Vryr5zegu1Kqi1LKFbgWOH3Vt45pV1ecNLT59QCJ7eH6J5bWOG2Fiw2Xcs2aOycTN3wEeZkpAGx68R5cKsB+/VVceukDrLm1FwC29ExCEkrZ18/jTyV9gCEDJvDw/N3cU03SF0KIc9HgiV9rXQHcB3wPxALLtdZ7GjoO5e4OQNtjcLRXEDYnpxqnrXCz4VoGnX5KACD5UDS6spKAHUnEhSquuda4xfH+hz+n2BV6ROXgpMF7/OR63w4hhDhblrTxa62/Bb61Yt0nKHePk92hf7ml1mkrXZ0JzP39/sniglzi1n1FqwLYMy745JeGzcmJ4x7gUQZJQXDp9U/VT/BCCPEntNjXMjt5GnfPpLeB8y6pPfHbXU/9fiwtzGX3qvcBCLvohmrnSenjj5uL/HNxIUTj02ITv7OXDwAZ3f1qbeYBqPTyOKW/tDAPt0PJZPnB2ItuPmVcXC834ttC92vvr9uAhRCijrTYxB8cPoYcbwieVnttHyD06hmUO8HOXkZxlRcVEpRaRmqw7Q9fGpe9uYauSz9nzLjr6iVuIYT4sxrtffz1bdiE62HC9Wc07ciLbiDqkxD8flkP+z7DnpWJTxEUBnn/YdpAvyAC/YLqOlwhhKgzLbbGf7aG9B+Hu08rACoOHQTAKbhDLXMIIUTjJIn/LLh7tQLA7Ug2AIF9//CmCSGEaPQk8Z8FD+9WAARkVlDqDANGXmltQEIIcQ4k8Z8Fz1atjb+lkNkaOgWHWRyREEKcPUn8Z8HTL/Bkd55/7beACiFEYyWJ/yz4+f9+t05xwB/v6BFCiKZAEv9Z8PL2P9ltC+lYy5RCCNF4SeI/C07OzqQazfwMuPIea4MRQohzJIn/LBU8fCfrbuvLgIEXWB2KEEKckxb75O65mnb534C/WR2GEEKcM6nxCyFECyOJXwghWhhJ/EII0cJI4hdCiBZGEr8QQrQwkviFEKKFkcQvhBAtjCR+IYRoYZTW2uoYTksplQUkWrDqACDbgvWeDYmx7jSFOJtCjNA04mwJMXbWWgdWHdgkEr9VlFKRWusIq+OojcRYd5pCnE0hRmgacbbkGKWpRwghWhhJ/EII0cJI4q/dfKsDOAMSY91pCnE2hRihacTZYmOUNn4hhGhhpMYvhBAtjCR+IYRoYVpU4ldKLVJKZSqldjsMC1dKbVNK7VRKRSqlhpnDlVLqTaXUQaVUtFJqsMM8tyilDpifWxogxoFKqa1KqRil1EqllK/DuMfMGPcrpSY6DL/EHHZQKfVoXcZoLr+jUupHpdRepdQepdQsc3hrpdRas2zWKqX8zeENXp61xHi12W9XSkVUmadBy7OWGF9VSu0zy+pLpVQrq2I8TZzPmTHuVEqtUUoFm8Mbzf52GP93pZRWSgVYFWNtcSqlnlFKHTHLcqdSarLDPHW7z7XWLeYDjAEGA7sdhq0BJpndk4ENDt2rAQWMAH4xh7cG4s2//ma3fz3H+Bsw1uyeDjxndvcBdgFuQBfgEOBkfg4BXQFXc5o+dVyW7YHBZrcPEGfG8wrwqDn8UeDfVpVnLTH2BnoCG4AIh+kbvDxrifFiwNkc/m+HcrRkn9cSp6/DNA8A7zS2/W32dwS+x3gQNMCqGE9Tls8AD1UzfZ3v8xZV49dabwKOVR0MnKhB+wGpZvdUYIk2bANaKaXaAxOBtVrrY1rrHGAtcEk9x9gD2GR2rwWudIjxE611qdb6MHAQGGZ+Dmqt47XWZcAn5rR1RmudprXebnYXALFAiLmexeZki4FpDrE2aHnWFKPWOlZrvb+aWRq8PGuJcY3WusKcbBvQwaoYTxNnvsNkXhjn04k4G8X+Nke/BjziEJ8lMZ5BnNWp833eohJ/DR4EXlVKJQNzgMfM4SFAssN0KeawmobXpz38vkOvxqi9NJoYlVKhwCDgF6Ct1jrNHJUOtG0MsVaJsSaNNcbpGDVTy2OsLk6l1Avm+XMD8M/GEKdjjEqpqcARrfWuKpM1urIE7jObnRadaCatjzgl8cPdwGytdUdgNrDQ4niqMx24RykVhfHTsMzieE5SSnkDK4AHq9T+0MbvVMvvF64txsaiphiVUk8AFcAyq2JzVF2cWusnzPNnGXCflfHBqTFilN3j/P6F1GhUU5bzgDAgHEgD/lNf65bED7cAX5jdn2H8fAI4wu81azB+ah+pZXi90Vrv01pfrLUeAnyM0a5neYxKKReMA3eZ1vpEGWaYP5cx/2ZaGWsNMdakUcWolLoV+Atwg/klalmMtcXpYBm/N0M2lrIMw2gX36WUSjDXt10p1c6qGGuIE611hta6UmttB96jPnNRXV2waCofIJRTL5zGAuPM7vFAlNk9hVMv/Pyqf7/wcxjjoo+/2d26nmMMMv/agCXAdLO/L6de9InHuODjbHZ34feLPn3rOEZlxvJ6leGvcurF3VesKs+aYnQYv4FTL+42eHnWUo6XAHuBwCrDLdnntcTZ3aH7fuDzxrq/zWkS+P3iriXneC1l2d6hezZGu3697PM62ZCm8sGoLacB5RjtYbcDo4Eos9B+AYY47Jy3MGrXMVUSxHSMCywHgdsaIMZZGFf+44CXMZ+4Nqd/woxxP+bdSebwyeb0h4An6qEsR2M040QDO83PZKANsA44APxw4oSxojxrifFys2xLgQzge6vKs5YYD2K0354Y9o6V+7yWOFcAu83hKzEu+Daq/V1lmgR+T/xWneM1leWHZhzRwNec+kVQp/tcXtkghBAtjLTxCyFECyOJXwghWhhJ/EII0cJI4hdCiBZGEr8QQrQwkviFEKKFkcQvhBAtzP8DI8Mull3z41sAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "for ozone_treatment in ['cmip6-old tuning/no feedback', 'stevenson', 'regression', 'thornhill-skeie']:\n", + " pl.plot(np.arange(1765, 2501), T[ozone_treatment], label=ozone_treatment)\n", + "pl.ylabel('K')\n", + "pl.title('Temperature change with different ozone treatments')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From ce8c2d54878b93da0fd9f6c1459012f65b8d0bab Mon Sep 17 00:00:00 2001 From: Chris Smith Date: Mon, 25 Jan 2021 21:04:46 +0000 Subject: [PATCH 6/6] presumtiously bump version with changelog, necessary for dependency cascade --- CHANGELOG.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d248a5c3..57dd8521 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,6 +4,9 @@ Changelog master ------ +v1.6.2 +------ + (`#99 `_) Add option to use ozone forcing dependent on N2O concentrations v1.6.1