{ "cells": [ { "cell_type": "markdown", "id": "9a574759-032d-43f1-a1f4-cdb53eb9692a", "metadata": {}, "source": [ "# 0e - Estimate the response timescale for every glacier model and experiment (for the shifted simulations)\n", "\n", "- response timescale estimated at mass shrinkage rates of -50%, -80% and -90%. The response timescales is only computed if at least 25% of the mass is lost (otherwise the noise is too large)...\n", "\n", "The following file is created, which is later used in other notebooks for the analysis and visualisations:\n", "`'../data/resp_time_shifted_for_deltaT_rgi_reg_roll_volume_21yravg.csv'`" ] }, { "cell_type": "code", "execution_count": 1, "id": "6d79abd8-8da7-4ee9-88a5-8d69ec1a81b8", "metadata": {}, "outputs": [], "source": [ "\n", "DATE = 'Feb12_2024' \n", "fill_option = 'repeat_last_101yrs'\n", "\n", "import xarray as xr\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "\n", "approach = '_via_5yravg'\n", "\n", "\n", "# the non-shifted stuff is just read in for some tests ... \n", "path_merged_runs_scaled_extend = f'../data/GMIP3_reg_glacier_model_data/glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}.nc'\n", "ds_reg_models_extend = xr.open_dataset(path_merged_runs_scaled_extend)\n", "ds_reg_yr_shift= xr.open_dataset(f'../data/GMIP3_reg_glacier_model_data/all_shifted_glacierMIP3_{DATE}_models_all_rgi_regions_sum_scaled_extended_{fill_option}{approach}.nc')\n", "\n", "# only select final submission data\n", "glac_models = ['CISM2', 'GO', 'PyGEM-OGGM_v13', 'GloGEMflow', 'Kraaijenbrink',\n", " 'GLIMB','OGGM_v16', 'GloGEMflow3D']\n", "ds_reg_models_extend = ds_reg_models_extend.sel(model_author = glac_models)\n", "ds_reg_yr_shift = ds_reg_yr_shift.sel(model_author = glac_models)\n", "\n", "pd_avg_exps= pd.read_csv('../data/climate_input_data/temp_ch_ipcc_ar6_isimip3b.csv', index_col=[0])\n", "pd_global_temp_exp_m = pd_avg_exps.set_index(['gcm', 'period_scenario'])" ] }, { "cell_type": "code", "execution_count": 37, "id": "e78d082b-64cc-4e45-8563-a20e58745afb", "metadata": {}, "outputs": [], "source": [ "global_models = ds_reg_models_extend.sel(simulation_year=0).dropna(dim='model_author',how='any').model_author.values" ] }, { "cell_type": "code", "execution_count": 38, "id": "3eea1380-031e-44f9-ae9a-9febcebb2ac5", "metadata": {}, "outputs": [], "source": [ "def resp_time_estimate(model_author='OGGM_v16', perc_change_l=[-50],\n", " rgi_reg='11', min_perc_change=25, roll_volume=21, shift_years=False):\n", " # computes the response timescale estimates for all 80 experiments of a given model_author and rgi_reg \n", " # perc_change : -50 (gives estimates of how fast 50% of total volume change occurs), also -80 is used in the manuscript\n", " # min_perc_change: minimum shrinkage/growing to estimate a \"response time\" (%)\n", " # HERE, the first occurrence of -50% change over the roll_volume avg. is used, \n", " # alternatively we could use the last occurrence, but this would \"mean\" something different of course ... \n", " # we also select the equilibrium volume by averaging over the last 100 years,\n", " # this is necessary due to the large interdecadal variability of some models\n", " # MAYBE could even select there 200 or 300 years (but then we would also need to check how the data is \"extended\"...\n", " if shift_years:\n", " # we shift by maximum by +50 years, so like that there should always be values inside \n", " _ds_reg_models_extend = ds_reg_yr_shift.sel(year_after_2020=slice(0,4950))\n", " else:\n", " _ds_reg_models_extend = ds_reg_models_extend\n", " if rgi_reg == 'All':\n", " if model_author in global_models:\n", " ds_reg_models_vol_all = _ds_reg_models_extend.volume_m3.sum(dim='rgi_reg')\n", " else:\n", " raise Exception(\"not a global model!!!\")\n", " \n", " else:\n", " ds_reg_models_vol_all = _ds_reg_models_extend.volume_m3.sel(rgi_reg=rgi_reg) # need to do the interpolation first !\n", "\n", " ds_reg_models_vol_all = ds_reg_models_vol_all.stack(experiments=['gcm','period_scenario'])\n", " # merge \"gcm\"s and \"period_scenario\" into one coordinate callled \"experiments\"\n", " ds_reg_models_vol_all = ds_reg_models_vol_all.assign_coords(temp_ch_ipcc = ('experiments',\n", " pd_global_temp_exp_m.loc[ds_reg_models_vol_all.experiments, \n", " 'temp_ch_ipcc']))\n", " ds_reg_models_sel_vol_all = ds_reg_models_vol_all.sel(model_author=model_author)\n", " # steady state is the average of the last 100 years!!!\n", " if shift_years:\n", " v_eq = ds_reg_models_sel_vol_all.sel(year_after_2020=slice(4901-50,5000-50)).mean(dim='year_after_2020')\n", " v0 = ds_reg_models_sel_vol_all.sel(year_after_2020=0) \n", " else:\n", " v_eq = ds_reg_models_sel_vol_all.sel(simulation_year=slice(4901,5000)).mean(dim='simulation_year')\n", " v0 = ds_reg_models_sel_vol_all.isel(simulation_year=0)\n", " \n", " pd_category_resp_time = pd.DataFrame(index=ds_reg_models_vol_all.experiments)\n", " if rgi_reg == 'All':\n", " pd_category_resp_time['rgi_reg'] = f'All'\n", " else:\n", " pd_category_resp_time['rgi_reg'] = f'RGI{rgi_reg}'\n", " pd_category_resp_time['model_author'] = model_author\n", " pd_category_resp_time['min_perc_change'] = min_perc_change\n", "\n", " pd_category_resp_time.loc[ds_reg_models_vol_all.experiments, 'temp_ch_ipcc'] = ds_reg_models_vol_all.temp_ch_ipcc.values\n", " # did it grow, shrink or stay the same??? ... this information is later not anymore used\n", " pd_category_resp_time['category'] = r'similar to V$_{0}$'\n", " grow_exps = v_eq.where(v_eq/v0 >=1+min_perc_change/100, drop=True).experiments.values\n", " pd_category_resp_time.loc[grow_exps, 'category'] = 'grow'\n", " shrink_exps = v_eq.where(v_eq/v0 <=1-min_perc_change/100, drop=True).experiments.values\n", " pd_category_resp_time.loc[shrink_exps, 'category'] = 'shrink'\n", " pd_category_resp_time['shift_years'] = shift_years\n", " for perc_change in perc_change_l:\n", " # response time computation\n", " # V(resp_t) = -perc/100 * (V_eq-v0)+v0)\n", " # we need the absolute because we only want to have either shrinking (negative perc) \n", " # or growing (positive perc) experiments\n", " v_resp_t = v0 + (perc_change/100)*np.abs(v0-v_eq)\n", " if shift_years:\n", " v_roll = ds_reg_models_sel_vol_all.rolling(year_after_2020=roll_volume, center=True).mean()\n", " # we take the last year where the anti-condition holds true (so it is basically the first year where the condition holds true)\n", " resp_t = v_roll.where(v_roll<=v_resp_t).idxmax(dim='year_after_2020')\n", " #print(v_roll, resp_t)\n", " else:\n", " v_roll = ds_reg_models_sel_vol_all.rolling(simulation_year=roll_volume, center=True).mean()\n", " # we take the last year where the anti-condition holds true (so it is basically the first year where the condition holds true)\n", " resp_t = v_roll.where(v_roll<=v_resp_t).idxmax(dim='simulation_year')\n", " if perc_change<0:\n", " # if we look at how much it shrinks it should at least shrink by min_perc_change to estimate a \"response time\"\n", " resp_t = resp_t.where(v_eq/v0 <=1-min_perc_change/100, np.NaN)\n", " else:\n", " # if we look at how much it grows it should at least grow by min_perc_change to estimate a \"response time\"\n", " resp_t = resp_t.where(v_eq/v0 >=1+min_perc_change/100, np.NaN)\n", " if perc_change>0:\n", " perc_change_i = f'+{perc_change}'\n", " else:\n", " perc_change_i = perc_change\n", " pd_category_resp_time.loc[resp_t.experiments,f'resp_time_{perc_change_i}%'] = resp_t.values\n", " \n", " return pd_category_resp_time" ] }, { "cell_type": "markdown", "id": "88f5c4e3-4b2b-4df7-a787-5666419f460a", "metadata": {}, "source": [ "just some random response timescale test estimates ... with and without shifting the experiments to start in year 2020..." ] }, { "cell_type": "code", "execution_count": 39, "id": "aeae7971-5833-4477-9b29-10cc0115d65f", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rgi_regmodel_authormin_perc_changetemp_ch_ipcccategoryshift_yearsresp_time_-50%
(gfdl-esm4, 1851-1870_hist)AllGLIMB250.231409similar to V$_{0}$FalseNaN
(gfdl-esm4, 1901-1920_hist)AllGLIMB250.478289shrinkFalse346.0
(gfdl-esm4, 1951-1970_hist)AllGLIMB250.392281shrinkFalse370.0
(gfdl-esm4, 1995-2014_hist)AllGLIMB250.901467shrinkFalse284.0
(gfdl-esm4, 2021-2040_ssp126)AllGLIMB251.493792shrinkFalse237.0
........................
(ukesm1-0-ll, 2061-2080_ssp370)AllGLIMB254.439977shrinkFalse72.0
(ukesm1-0-ll, 2061-2080_ssp585)AllGLIMB255.230543shrinkFalse59.0
(ukesm1-0-ll, 2081-2100_ssp126)AllGLIMB253.038482shrinkFalse99.0
(ukesm1-0-ll, 2081-2100_ssp370)AllGLIMB255.840495shrinkFalse51.0
(ukesm1-0-ll, 2081-2100_ssp585)AllGLIMB256.884361shrinkFalse40.0
\n", "

80 rows × 7 columns

\n", "
" ], "text/plain": [ " rgi_reg model_author min_perc_change \\\n", "(gfdl-esm4, 1851-1870_hist) All GLIMB 25 \n", "(gfdl-esm4, 1901-1920_hist) All GLIMB 25 \n", "(gfdl-esm4, 1951-1970_hist) All GLIMB 25 \n", "(gfdl-esm4, 1995-2014_hist) All GLIMB 25 \n", "(gfdl-esm4, 2021-2040_ssp126) All GLIMB 25 \n", "... ... ... ... \n", "(ukesm1-0-ll, 2061-2080_ssp370) All GLIMB 25 \n", "(ukesm1-0-ll, 2061-2080_ssp585) All GLIMB 25 \n", "(ukesm1-0-ll, 2081-2100_ssp126) All GLIMB 25 \n", "(ukesm1-0-ll, 2081-2100_ssp370) All GLIMB 25 \n", "(ukesm1-0-ll, 2081-2100_ssp585) All GLIMB 25 \n", "\n", " temp_ch_ipcc category \\\n", "(gfdl-esm4, 1851-1870_hist) 0.231409 similar to V$_{0}$ \n", "(gfdl-esm4, 1901-1920_hist) 0.478289 shrink \n", "(gfdl-esm4, 1951-1970_hist) 0.392281 shrink \n", "(gfdl-esm4, 1995-2014_hist) 0.901467 shrink \n", "(gfdl-esm4, 2021-2040_ssp126) 1.493792 shrink \n", "... ... ... \n", "(ukesm1-0-ll, 2061-2080_ssp370) 4.439977 shrink \n", "(ukesm1-0-ll, 2061-2080_ssp585) 5.230543 shrink \n", "(ukesm1-0-ll, 2081-2100_ssp126) 3.038482 shrink \n", "(ukesm1-0-ll, 2081-2100_ssp370) 5.840495 shrink \n", "(ukesm1-0-ll, 2081-2100_ssp585) 6.884361 shrink \n", "\n", " shift_years resp_time_-50% \n", "(gfdl-esm4, 1851-1870_hist) False NaN \n", "(gfdl-esm4, 1901-1920_hist) False 346.0 \n", "(gfdl-esm4, 1951-1970_hist) False 370.0 \n", "(gfdl-esm4, 1995-2014_hist) False 284.0 \n", "(gfdl-esm4, 2021-2040_ssp126) False 237.0 \n", "... ... ... \n", "(ukesm1-0-ll, 2061-2080_ssp370) False 72.0 \n", "(ukesm1-0-ll, 2061-2080_ssp585) False 59.0 \n", "(ukesm1-0-ll, 2081-2100_ssp126) False 99.0 \n", "(ukesm1-0-ll, 2081-2100_ssp370) False 51.0 \n", "(ukesm1-0-ll, 2081-2100_ssp585) False 40.0 \n", "\n", "[80 rows x 7 columns]" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "resp_time_estimate(model_author='GLIMB', rgi_reg='All')" ] }, { "cell_type": "code", "execution_count": 40, "id": "a30d47d0-39a9-4a45-9d84-28692d3b77b7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGxCAYAAACKvAkXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSAUlEQVR4nO3dfVxUZd4/8M+IAoowMyoMcIuJmWIm6JoKma41OOa2PmRleFfqJlsZlmZu/dysbHvAzW2zttItS9u7dWnbO80eVJQE71qxNPEhk5Q0cIXxaZgBFES4fn+cGBhmhnlghjPD+bxfr9mYcz19z5lee76dc53rqIQQAkREREQK0kXuAIiIiIg6GhMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLFYQJEREREisMEiIiIiBSnq9wBBKLGxkacPn0akZGRUKlUcodDREREbhBCoKqqCvHx8ejSpe1rPEyAHDh9+jQSEhLkDoOIiIi8UFZWhr59+7ZZhwmQA5GRkQCkAxgVFSVzNEREROQOi8WChIQE63m8LUyAHGi67RUVFcUEiIiIKMi4M32Fk6CJiIhIcZgAERERkeIwASIiIiLFYQJEREREisMEiIiIiBSHCRAREREpDhMgIiIiUhwmQERERKQ4TICIiIhIcWRNgFavXo3k5GTristpaWnYsmWLtby2thZZWVno3bs3evbsidtvvx1Go7HNPoUQePrppxEXF4fu3bsjPT0dx44d8/euEBERURCRNQHq27cvVqxYgX379mHv3r24+eabMW3aNHz33XcAgEcffRSffPIJPvzwQxQUFOD06dOYMWNGm32+9NJLeO2117BmzRrs2bMHERERmDRpEmpraztilzo1kwkoKQGKioAvvwSOHJG2ERERBRuVEELIHURLvXr1wsqVK3HHHXcgOjoaGzZswB133AEAOHr0KIYMGYLdu3cjNTXVrq0QAvHx8XjsscewZMkSAIDZbIZOp8P69euRkZHhVgwWiwVqtRpms5nvAvvZqVPA8ePA888DeXnN2w0GYO1aICFBvtiIiIgAz87fATMHqKGhATk5OaipqUFaWhr27duH+vp6pKenW+skJSWhX79+2L17t8M+Tpw4gYqKCps2arUaY8aMcdqGXDOZgC1b7JMfAMjNBTIzeSWIiIiCi+xvgz906BDS0tJQW1uLnj17YuPGjbj22mtRVFSE0NBQaDQam/o6nQ4VFRUO+2rartPp3G4DAHV1dairq7N+t1gsXu5N52Q0AnFx9slPk9xcqY5W27FxEREReUv2K0CDBw9GUVER9uzZg/nz52POnDk4cuRIh8aQnZ0NtVpt/STwfo4NsxlwNYXKbO6YWIiIiHxB9gQoNDQUAwcOxMiRI5GdnY2UlBS8+uqriI2NxeXLl1FZWWlT32g0IjY21mFfTdtbPynWVhsAWLp0Kcxms/VTVlbWvp3qZNRqIDzcdR0iIqJgIXsC1FpjYyPq6uowcuRIdOvWDXkt7rsUFxejtLQUaWlpDtsmJiYiNjbWpo3FYsGePXuctgGAsLAw66P4TR9qptMB5eWAXu+43GCQ6hAREQULWecALV26FJMnT0a/fv1QVVWFDRs2ID8/H9u2bYNarca8efOwePFi9OrVC1FRUXj44YeRlpZm8wRYUlISsrOzcdttt0GlUmHRokV4/vnncc011yAxMRFPPfUU4uPjMX36dPl2NMhptcDkycA110jfHT0Fxvk/REQUTGRNgM6cOYPZs2ejvLwcarUaycnJ2LZtGyZOnAgAeOWVV9ClSxfcfvvtqKurw6RJk/Dmm2/a9FFcXAxziwkojz/+OGpqanD//fejsrISN954I7Zu3YpwV/dwqE19+wIREcBf/wpUVQE1NVLSExfH5IeIiIJPwK0DFAi4DhAREVHwCcp1gIiIiIg6ChMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLFYQJEREREisMEiIiIiBSHCRAREREpDhMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLFYQJEREREisMEiIiIiBSHCRAREREpDhMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLFYQJEREREitNV7gCUyGQCjEbAbAY0GiAmBtBq5Y7Kdzr7/hERUfDjFaAOVlYGZGQAQ4YAqalAUpL0vaxM7sh8o7PvHxERdQ5MgDqQyQRkZgK5ubbbc3Ol7SaTPHH5SmffPyIi6jyYAHUgo9E+OWiSmyuVB7POvn9ERNR5MAHqQGZz+8oDXWffPyIi6jyYAHUgtbp95YGus+8fERF1HkyAOpBOBxgMjssMBqk8mHX2/SMios6DCVAH0mqBtWvtkwSDQdoe7I+Kd/b9IyKizkPWBCg7OxujRo1CZGQkYmJiMH36dBQXF1vLT548CZVK5fDz4YcfOu137ty5dvVvueWWjtgllxISgJwc4PvvgcJC6Z85OdL2zqCz7x8REXUOsi6EWFBQgKysLIwaNQpXrlzB73//exgMBhw5cgQRERFISEhAeXm5TZu33noLK1euxOTJk9vs+5ZbbsG6deus38PCwvyyD97Qajv31ZDOvn9ERBT8ZE2Atm7davN9/fr1iImJwb59+zB+/HiEhIQgNjbWps7GjRsxc+ZM9OzZs82+w8LC7NoSERERAQE2B8j883PSvXr1cli+b98+FBUVYd68eS77ys/PR0xMDAYPHoz58+fj/PnzPo2ViIiIgpdKCCHkDgIAGhsbMXXqVFRWVuLLL790WOehhx5Cfn4+jhw50mZfOTk56NGjBxITE1FSUoLf//736NmzJ3bv3o2QkBC7+nV1dairq7N+t1gsSEhIgNlsRlRUVPt2jIiIiDqExWKBWq126/wdMC9DzcrKwuHDh50mP5cuXcKGDRvw1FNPuewrIyPD+vewYcOQnJyMq6++Gvn5+dDr9Xb1s7Oz8eyzz3ofPBEREQWVgLgFtmDBAnz66afYuXMn+vbt67DOv/71L1y8eBGzZ8/2uP8BAwagT58+OH78uMPypUuXwmw2Wz9lfHMnERFRpybrFSAhBB5++GFs3LgR+fn5SExMdFr3nXfewdSpUxEdHe3xOKdOncL58+cRFxfnsDwsLCygnhIjIiIi/5L1ClBWVhbef/99bNiwAZGRkaioqEBFRQUuXbpkU+/48ePYtWsXMjMzHfaTlJSEjRs3AgCqq6vxu9/9DoWFhTh58iTy8vIwbdo0DBw4EJMmTfL7PhEREVHgk/UK0OrVqwEAEyZMsNm+bt06zJ071/r93XffRd++fWFw8p6F4uJi6xNkISEhOHjwIN577z1UVlYiPj4eBoMBzz33nOxXeUwm6Y3oZjOg0QAxMc7Xy2mqazIBPXsCISFA165AdDTX2CEiImqvgHkKLJB4MovcXWVlQGYmkJvbvK3pFRGtV0l2VFevBxYulOq//jpXViYiImrNk/N3QEyC7uxMJvuEBpC+Z2ZK5a7q5uUBr74KDBtm34aIiIg8wwSoAxiN9glNk9xcqdydunl5QGqqfRsiIiLyDBOgDvDz9CS3yl3Vra11rx4RERE5xwSoA6jV7pe7qhse7l49IiIico4JUAfQ6aQJz44YDFK5O3X1eqCw0L4NEREReYYJUAfQaqWnt1onNk1PgbV8rN1Z3aanwA4dsm9DREREnuFj8A744zF4wHYdILVauorjzjpAERHSGkAhIW2vHURERKRkQfkyVCXQat1PXjypS0RERJ7hLTAiIiJSHCZAREREpDhMgIiIiEhxmAARERGR4jABIiIiIsVhAkRERESKwwSIiIiIFIfrAJHbWi7kqNFwUUYiIgpevAJEbikrAzIygCFDgNRUIClJ+l5WJndkREREnmMCRC6ZTEBmJpCba7s9N1fabjLJExcREZG3mACRS0ajffLTJDdXKiciIgomTIDIJbO5feVERESBhgkQuaRWt6+ciIgo0DABIpd0OsBgcFxmMEjlREREwYQJELmk1QJr19onQQaDtJ2PwhMRUbDhOkDkloQEICeneR0gtVq68sPkh4iIghETIHKbVsuEh4iIOgfeAiMiIiLFYQJEREREisMEiIiIiBSHCRAREREpDhMgIiIiUhwmQERERKQ4siZA2dnZGDVqFCIjIxETE4Pp06ejuLjYps6ECROgUqlsPg8++GCb/Qoh8PTTTyMuLg7du3dHeno6jh075s9dISIioiAiawJUUFCArKwsFBYWYvv27aivr4fBYEBNTY1Nvd/+9rcoLy+3fl566aU2+33ppZfw2muvYc2aNdizZw8iIiIwadIk1NbW+nN3OpTJBBw9CuzZAxQXS9+JiIjIPbIuhLh161ab7+vXr0dMTAz27duH8ePHW7f36NEDsbGxbvUphMCqVauwbNkyTJs2DQDwt7/9DTqdDps2bUJGRobvdkAmZWVAZiaQm9u8rem1FAkJ8sVFREQULAJqDpDZbAYA9OrVy2b73//+d/Tp0wfXXXcdli5diosXLzrt48SJE6ioqEB6erp1m1qtxpgxY7B7927/BN6BTCb75AeQvmdm8koQERGROwLmVRiNjY1YtGgRxo4di+uuu866/b//+79x1VVXIT4+HgcPHsQTTzyB4uJifPTRRw77qaioAADoWr2iXKfTWctaq6urQ11dnfW7xWJp7+74jdFon/w0yc2Vyvm6CiIiorYFTAKUlZWFw4cP48svv7TZfv/991v/HjZsGOLi4qDX61FSUoKrr77aJ2NnZ2fj2Wef9Ulf/vbzRTKvy4mIiChAboEtWLAAn376KXbu3Im+ffu2WXfMmDEAgOPHjzssb5orZDQabbYbjUan84iWLl0Ks9ls/ZSVlXm6Cx1GrW5fOREREcmcAAkhsGDBAmzcuBFffPEFEhMTXbYpKioCAMTFxTksT0xMRGxsLPLy8qzbLBYL9uzZg7S0NIdtwsLCEBUVZfMJVDqdNOHZEYNBKiciIqK2yZoAZWVl4f3338eGDRsQGRmJiooKVFRU4NKlSwCAkpISPPfcc9i3bx9OnjyJzZs3Y/bs2Rg/fjySk5Ot/SQlJWHjxo0AAJVKhUWLFuH555/H5s2bcejQIcyePRvx8fGYPn26HLvpU1qt9LRX6ySo6Skwzv8hIiJyTdY5QKtXrwYgLXbY0rp16zB37lyEhoZix44dWLVqFWpqapCQkIDbb78dy5Yts6lfXFxsfYIMAB5//HHU1NTg/vvvR2VlJW688UZs3boV4eHhft+njpCQAOTkSBOezWbptpdOx+SHiIjIXSohhJA7iEBjsVigVqthNpsD+nYYERERNfPk/B0Qk6CJiIiIOhITICIiIlIcJkBERESkOEyAiIiISHGYABEREZHiMAEiIiIixWECRERERIoTMC9DpcBnMjUvvqjRADExXHyRiIiCE68AkVvKyoCMDGDIECA1FUhKkr4H8HtjiYiInGICRC6ZTEBmJpCba7s9N1fabjLJExcREZG3mACRS0ajffLTJDdXKiciIgomTIDIpRbvmfWqnIiIKNAwASKX1Or2lRMREQUaJkDkkk4HGAyOywwGqZyIiCiYMAEil7RaYO1a+yTIYJC281F4IiIKNlwHiNySkADk5DSvA6RWS1d+mPwQEVEwYgJEbtNqmfAQEVHnwFtgREREpDhMgIiIiEhxmAARERGR4jABIiIiIsVhAkRERESKwwSIiIiIFIcJEBERESkOEyAiIiJSHJ8kQDU1NbBYLL7oioiIiMjv2pUAHTlyBNdffz0iIyOh1WoxbNgw7Nu3z1exEREREflFuxKgBx54AAsWLEB1dTXOnz+PGTNmYPbs2b6KjYiIiMgvPEqApk2bhv/85z/W72fPnsXUqVPRo0cPaDQa/OpXv4LRaPR5kERERES+5NHLUO+55x7cfPPNyMrKwsMPP4wFCxZg6NCh+OUvf4n6+np88cUXeOyxx/wVKxEREZFPqIQQwpMGZrMZTzzxBPbv3481a9aga9euyM/PR0NDA8aOHYtRo0b5K9YOY7FYoFarYTabERUVJXc4RERE5AZPzt8ezwFSq9VYs2YNXn75ZcyZMwfr16/HvHnzsGjRIo+Tn+zsbIwaNQqRkZGIiYnB9OnTUVxcbC2/cOECHn74YQwePBjdu3dHv3798Mgjj8BsNrfZ79y5c6FSqWw+t9xyi6e7SkRERJ2UxwnQhQsXsG/fPusTX1FRURgxYgQ+//xzjwcvKChAVlYWCgsLsX37dtTX18NgMKCmpgYAcPr0aZw+fRp/+tOfcPjwYaxfvx5bt27FvHnzXPZ9yy23oLy83Pr5xz/+4XF8RERE1Dl5dAtsw4YNyMzMRFRUFGpra/G3v/0NU6dOxdGjR/Hggw8iJiYGf/nLX6DT6bwK5uzZs4iJiUFBQQHGjx/vsM6HH36Ie+65BzU1Neja1fEUprlz56KyshKbNm3yKg7eAiMiIgo+frsFtnTpUrz77ruoqKhAXl4ennrqKQBAUlIS8vPzMXHiRKSlpXkdeNOtrV69erVZJyoqymny0yQ/Px8xMTEYPHgw5s+fj/Pnz3sdFxEREXUuHl0B6t27N3bs2IERI0agsrISo0aNwrFjx2zqnDlzBjExMR4H0tjYiKlTp6KyshJffvmlwzrnzp3DyJEjcc899+CFF15w2ldOTg569OiBxMRElJSU4Pe//z169uyJ3bt3IyQkxK5+XV0d6urqrN8tFgsSEhJ4BYiIiCiIeHIFyKPH4OfMmYNbb70VEyZMwN69e3Hvvffa1fEm+QGArKwsHD582GnyY7FYcOutt+Laa6/F8uXL2+wrIyPD+vewYcOQnJyMq6++Gvn5+dDr9Xb1s7Oz8eyzz3oVNxEREQUfjx+D/+STT3D06FGkpKTAYDD4JIgFCxbg448/xq5du5CYmGhXXlVVhUmTJqFHjx749NNPER4e7vEY0dHReP755/HAAw/YlfEKEBERUfDz2xUgAJgyZQqmTJnidXAtCSHw8MMPY+PGjcjPz3eY/FgsFkyaNAlhYWHYvHmzV8nPqVOncP78ecTFxTksDwsLQ1hYmMf9EhERUXDyOAE6d+4c3n33XezevRsVFRUAgNjYWNxwww2YO3cuoqOj3e4rKysLGzZswMcff4zIyEhrf2q1Gt27d4fFYoHBYMDFixfx/vvvw2KxWN86Hx0dbZ3Pk5SUhOzsbNx2222orq7Gs88+i9tvvx2xsbEoKSnB448/joEDB2LSpEme7i4RERF1Qh7dAvvmm2+st6LS09Otj7sbjUbk5eXh4sWL2LZtG66//nr3BlepHG5ft24d5s6di/z8fNx0000O65w4cQL9+/e39tPU5tKlS5g+fTr279+PyspKxMfHw2Aw4LnnnnP78fxAeAzeZAKMRsBsBjQaICYG0GplCSXo8NgRESmTJ+dvjxKg1NRUpKSkYM2aNXbJixACDz74IA4ePIjdu3d7F3mAkDsBKisDMjOB3NzmbQYDsHYtkJDQ4eEEFR47IiLl8lsC1L17d+zfvx9JSUkOy48ePYoRI0bg0qVLnkUcYORMgEwmICPD9gTexGAAcnJ4NcMZHjsiImXz20KIsbGx+Prrr52Wf/31116vAk0So9HxCRyQthuNHRtPMOGxIyIid3k0CXrJkiW4//77sW/fPuj1ers5QG+//Tb+9Kc/+SVQpXDxnleX5UrGY0dERO7yKAHKyspCnz598Morr+DNN99EQ0MDACAkJAQjR47E+vXrMXPmTL8EqhRqdfvKlYzHjoiI3OXxY/B33XUX7rrrLtTX1+PcuXMAgD59+qBbt24+D06JdDppvoqzeSy8w+gcjx0REbnLozlALXXr1g1xcXGIi4tj8uNDWq30xFLrRbabnmTiJF7neOyIiMhdHj0F1vpFp0VFRXjllVdw/PhxxMXFYcGCBZgwYYI/4uxQcj8GD9iuZaNWS1cveAJ3D48dEZEy+e1VGHFxcSgvL0dMTAz+/e9/Y8KECbjhhhswduxYFBUVYeLEicjLy8P48ePbtQMknbB50vYOjx0REbni0RWgLl26oKKiAjExMTAYDEhISMA777xjLV+0aBEOHTqEvLw8vwTbUQLhChARERF5xm/rALV0+PBh/Pa3v7XZ9tvf/hYHDx70tksiIiKiDuHxU2BVVVUIDw9HeHi43RvUw8PDcfHiRZ8FR0REROQPHl8BGjRoELRaLU6ePIm9e/falH333XeIj4/3WXBERERE/uDRFaCdO3fafI+Li7P5fuLECdx///3tj4qIiIjIjzyaBK0UnARNREQUfDpkEnSThx56yLoiNBEREVEwaPcVoKioKBQVFWHAgAG+ikl2vALkRMsVBjUaICYmIBbc8SYskwkoLwcuXAAiI4GePYFevQJid4iIyEsdegWId9AUoqwMyMgAhgwBUlOBpCTpe1lZ0IXV1GboUGDcOGD4cOCBB4ADB4BTpzosdCIiklG7EyBSAJMJyMy0f8tobq603WQKmrCctcnLA55/HtiyRbbdISKiDtTuBKiqqsrp7a8VK1agsrKyvUOQ3IxGx69YB6TtRmPHxvMzb8Jqq01eHhAXJ9vuEBFRB/LrFaAXX3wRFy5c8OcQ1BHM5vaV+4k3YblqU1sr2+4QEVEH8msCxPlBnYRa3b5yP/EmLFdtwsNl2x0iIupAnANErul0gMHguMxgkMpl4E1YbbXR66Unw2TaHSIi6kBMgMg1rRZYu9Y+czAYpO0yPTvuTVjO2uj1wLJlwOTJfBSeiEgJPH4ZKilUQgKQk9O84I5aLV0qkTlb8Caspjbl5dITXxER0lpAXAeIiEg5mACR+7TagMwQvAkrQHeFiIg6iF9vgY0bNw7du3f35xBEREREHvM6ASopKcGyZcswa9YsnDlzBgCwZcsWfPfdd9Y6n3/+ud0b44mIiIjk5lUCVFBQgGHDhmHPnj346KOPUF1dDQA4cOAAnnnmGZ8GSERERORrXiVA/+///T88//zz2L59O0JDQ63bb775ZhQWFvosOCIiIiJ/8CoBOnToEG677Ta77TExMTh37ly7gyIiIiLyJ68SII1Gg/Lycrvt+/fvx3/913+1OygiIiIif/IqAcrIyMATTzyBiooKqFQqNDY24quvvsKSJUswe/Zst/vJzs7GqFGjEBkZiZiYGEyfPh3FxcU2dWpra5GVlYXevXujZ8+euP3222F08bZKIQSefvppxMXFoXv37khPT8exY8e82VUiIiLqhLxKgF588UUkJSUhISEB1dXVuPbaazF+/HjccMMNWLZsmdv9FBQUICsrC4WFhdi+fTvq6+thMBhQU1NjrfPoo4/ik08+wYcffoiCggKcPn0aM2bMaLPfl156Ca+99hrWrFmDPXv2ICIiApMmTUJtba03u0tNTCbg6FFgzx6guFj63p56MoTmu4ZERBTURDv89NNP4rPPPhMffPCB+OGHH9rTlRBCiDNnzggAoqCgQAghRGVlpejWrZv48MMPrXW+//57AUDs3r3bYR+NjY0iNjZWrFy50rqtsrJShIWFiX/84x9uxWE2mwUAYTab27E3nUxpqRAGgxBA88dgkLZ7U0+G0HzXkIiIApEn5+92LYTYr18//OpXv8LMmTNxzTXXtDsZM5vNAIBevXoBAPbt24f6+nqkp6db6yQlJaFfv37YvXu3wz5OnDiBiooKmzZqtRpjxoxx2oZcMJmAzEwgN9d2e26utL3pqom79WQIzXcNiYioM/DqVRhCCPzrX//Czp07cebMGTQ2NtqUf/TRRx732djYiEWLFmHs2LG47rrrAAAVFRUIDQ2FRqOxqavT6VBRUeGwn6btulav9G6rTV1dHerq6qzfLRaLx/F3akajfaLQJDdXKtdq3a8nQ2i+a0hERJ2BV1eAFi1ahHvvvRcnTpxAz549oVarbT7eyMrKwuHDh5GTk+NV+/bIzs62iT8hIaHDYwhoP1+Zc1nubj0f8npIGWIlIqLA4dUVoP/5n//BRx99hF/96lc+CWLBggX49NNPsWvXLvTt29e6PTY2FpcvX0ZlZaXNVSCj0YjY2FiHfTVtNxqNNq/hMBqNGD58uMM2S5cuxeLFi63fLRYLk6CWXCW1TeXu1vMhr4eUIVYiIgocXl0BUqvVGDBgQLsHF0JgwYIF2LhxI7744gskJibalI8cORLdunVDXl6edVtxcTFKS0uRlpbmsM/ExETExsbatLFYLNizZ4/TNmFhYYiKirL5UAs6HWAwOC4zGKRyT+rJEJrvGhIRUafgzSzr9evXi4yMDHHx4kVvmlvNnz9fqNVqkZ+fL8rLy62flv0++OCDol+/fuKLL74Qe/fuFWlpaSItLc2mn8GDB4uPPvrI+n3FihVCo9GIjz/+WBw8eFBMmzZNJCYmikuXLrkVF58Cc4BPgRERUYDz5PytEkIIT5OmS5cu4bbbbsNXX32F/v37o1u3bjbl3377rVv9qFQqh9vXrVuHuXPnApAWQnzsscfwj3/8A3V1dZg0aRLefPNNm1tgKpXKpo0QAs888wzeeustVFZW4sYbb8Sbb76JQYMGuRWXxWKBWq2G2Wzm1aCWTCZpcrDZLN0i0ukcTxR2t54MoQVCrERE5B+enL+9SoBmzpyJnTt34o477oBOp7NLZIL9jfBMgIiIiIKPJ+dvryZBf/bZZ9i2bRtuvPFGrwIkIiIikpNXk6ATEhJ4ZYSIiIiCllcJ0Msvv4zHH38cJ0+e9HE4RERERP7n1S2we+65BxcvXsTVV1+NHj162E2CvnDhgk+CIyIiIvIHrxKgVatW+TgMIiIioo7jVQI0Z84cX8dBRERE1GHcToAsFot14rOrl4VygjQREREFMrcTIK1Wi/LycsTExECj0ThcxFAIAZVKhYaGBp8G2Wm0XHRPowFiYlwvuuesjTd9dWK15SZ0OWuEMJvRRauBOSwGJyq1iIpq+9DwMBIRKZPbCdAXX3yBXr16AQB27tzpt4A6rbIyIDMTyM1t3mYwAGvXAs5evOqszZtvAosXA5s3u99XJ3a5pAzdHsxEyI7m46TRG9B14VqMvCkBY8c6PjTe/CRERNQ5eLUSdGlpKRISEuyuAgkhUFZWhn79+vksQDn4fCVokwnIyLA90zYxGICcHPvLDm21SU8HxowBXnjBvb46sdpyE7rNzrBJfppc0RvwamoOlrygtTs03vwkREQU2Dw5f3u1DlBiYiLOnj1rt/3ChQt2b3QnSPdYHJ1pAWm70ehZmx07gNRU9/vqxLqcNTpMfgCga14upqZKx6P1ofHmJyEios7DqwSoaa5Pa9XV1QgPD293UJ2O2ex5uas2tbXejdXJCBf7G1bbXN6yqjc/CRERdR4ePQa/ePFiANLb15966in06NHDWtbQ0IA9e/Zg+PDhPg2wU1CrPS931cZZoumqXSejcrG/deHN5S2revOTEBFR5+FRArR//34A0hWgQ4cOITQ01FoWGhqKlJQULFmyxLcRdgY6nTSxxNmEE53Oszbp6UBhoft9dWKN0To0pBuczgHaXCgdj9aHxpufhIiIOhHhhblz5wqz2eyyXllZmWhoaPBmCFmZzWYBwK19dFtpqRAGgxBA88dgkLZ72qakRIipUz3rqxOrO14qrqTbHqd6vUF8u7lUREQ4PzTe/CRERBS4PDl/e/UUmLuioqJQVFSEAQMG+GsIv/D5U2BNWi46o1ZLlxk8WQeoZRtv+urEWq4DpNKoYQnX4USlFpGRbR8aHkYios7Dk/O3V6/CcJcfc6vgpNV6fnZ11sabvjqx8DgtENd8PPr8/HGFh5GISJm8egqMiIiIKJgxASIiIiLFYQJEREREiuPXBMjRYolEREREcvNrAsRJ0ERERBSI2v0UWFlZGQAgwcHrs48cOYL4+Pj2DkFERETkU15dAbpy5QqeeuopqNVq9O/fH/3794darcayZctQX19vrZeQkICQkBCfBUvUXiYTcPQosGcPcOQIUFICfPMNUFwslRERkTJ4dQXo4YcfxkcffYSXXnoJaWlpAIDdu3dj+fLlOH/+PFavXu3TIIl8oawMyMy0ff2FXg8sXAjcdBMwdiywdi3g4GImERF1Ml6tBK1Wq5GTk4PJkyfbbP/8888xa9YsmIP8Vdp+WwmaZGMyARkZjt/9pdcDqanACy9I7wHLyeHiiEREwciT87dXt8DCwsLQv39/u+2JiYk2L0glChRGo+PkBwDy8qQECJDqGI0dFxcREcnDqwRowYIFeO6551BXV2fdVldXhxdeeAELFizwWXBEvuLqomRtrft1iYgo+Hk1B2j//v3Iy8tD3759kZKSAgA4cOAALl++DL1ejxkzZljrfvTRR76JlKgd1Oq2y8PD3a9LRETBz6sESKPR4Pbbb7fZ5ugxeKJAodNJ83uczQEqLJT+NhikukRE1Ll5NQm6s+Mk6M6prafAZs3iU2BERMHO75OgL126hIsXL1q///TTT1i1ahVync0ydWLXrl2YMmUK4uPjoVKpsGnTJptylUrl8LNy5UqnfS5fvtyuflJSkkdxUeeUkCA94fX999IVn+++A/76VyA2Fti7Vypj8kNEpAxe3QKbNm0aZsyYgQcffBCVlZUYPXo0QkNDce7cOfz5z3/G/Pnz3eqnpqYGKSkpuO+++2zmDTUpLy+3+b5lyxbMmzfP7vZba0OHDsWOHTus37t2bfeC19RJaLX2j7hffbU8sRARkXy8ygy+/fZbvPLKKwCAf/3rX4iNjcX+/fvxv//7v3j66afdToAmT55st5ZQS7GxsTbfP/74Y9x0000YMGBAm/127drVri0RERFRE69ugV28eBGRkZEAgNzcXMyYMQNdunRBamoqfvrpJ58G2MRoNOKzzz7DvHnzXNY9duwY4uPjMWDAANx9990oLS31S0xEREQUnLxKgAYOHIhNmzahrKwM27Ztg8FgAACcOXPGb5OG33vvPURGRjq8VdbSmDFjsH79emzduhWrV6/GiRMnMG7cOFRVVTltU1dXB4vFYvMhIiKizsurBOjpp5/GkiVL0L9/f4wePdr6PrDc3FyMGDHCpwE2effdd3H33XcjvOWCLQ5MnjwZd955J5KTkzFp0iR8/vnnqKysxD//+U+nbbKzs6FWq60fPtJPRETUuXmVAN1xxx0oLS3F3r17sW3bNut2vV5vnRvkS//3f/+H4uJiZGZmetxWo9Fg0KBBOH78uNM6S5cuhdlstn7KysraEy4REREFOK8SIECaoBwZGYnt27fj0qVLAIBRo0b55ZHzd955ByNHjrSuOu2J6upqlJSUIC4uzmmdsLAwREVF2XyIiIio8/IqATp//jz0ej0GDRqEX/3qV9bH1efNm4fHHnvM7X6qq6tRVFSEoqIiAMCJEydQVFRkM2nZYrHgww8/dHr1R6/X4/XXX7d+X7JkCQoKCnDy5En8+9//xm233YaQkBDMmjXLiz0lIiKizsirBOjRRx9Ft27dUFpaih49eli333XXXdi6davb/ezduxcjRoywzhtavHgxRowYgaefftpaJycnB0IIpwlMSUkJzp07Z/1+6tQpzJo1C4MHD8bMmTPRu3dvFBYWIjo62tPd9B+TCTh6FNizBygulr5T0OLPSUQUfLx6FUZsbCy2bduGlJQUREZG4sCBAxgwYAB+/PFHJCcno7q62h+xdhi/vgrD0fsYDAa+gyFI8eckIgocfn8VRk1Njc2VnyYXLlxAWFiYN10qg8lkf7YEpO+Zmbx0EGT4cxIRBS+vEqBx48bhb3/7m/W7SqVCY2MjXnrpJdx0000+C67TMRodv44ckLYbjR0bD7ULf04iouDl1aswVq5ciZtvvhl79+7F5cuX8fjjj+O7777DhQsX8NVXX/k6xs7DbG5fOQUU/pxERMHL4wSovr4ejzzyCD755BNs374dkZGRqK6uxowZM5CVldXm4+aKp1a3r5wCCn9OIqLg5XEC1K1bNxw8eBBarRZPPvmkP2LqvHQ6aYaso/smBoNUTkGDPycRUfDyag7QPffcg3feecfXsXR+Wq30eNDP706zanpsSKuVJy7yCn9OIqLg5dUcoCtXruDdd9/Fjh07MHLkSERERNiU//nPf/ZJcJ1SQgKQkyPNkDWbpfskOh3PlkGKPycRUXDyKgE6fPgwfvGLXwAAfvjhB5sylUrV/qg6O62WZ8hOhD8nEVHw8SoB2rlzp6/jICIiIuowXr8MlYiIiChYMQEiIiIixWECRERERIrDBIiIiIgUhwkQERERKQ4TICIiIlIcrx6DJwUzmZpX/dNogJgY+RfBMZnQUG5E4wUzGiI1uNgzBqpeWqdhmUxAwzkTImqMCKk2Q6XRwBwWgx/OahERAfTsCfTqJf9uERGR//AKELmvrAzIyACGDAFSU4GkJOl7WZmsMYmMDIQMHYJu41IRPjwJUQ9k4PyBMpw6ZV/91Cng/IEyaOdnoPuIIQgdl4puw5KgmZ+B7ufLMHYs8MADwIEDcNieiIg6ByZA5B6TCcjMtH/zZ26utN1kkiUmkZkJVauYuublov/zmSjcYrIJy2QCCreYkPh8JkLy7NsMezUTzy4yIS8PeP55YMsWeXaLiIj8jwkQucdodPzac0DabjR2bDwAYDTaJT9NuublIiXOaBOW0QikxBntkp+WbaamSg3y8oC4OHl2i4iI/I8JELnHbG5fuT+4GDOs1mxTxWyWtrlq06S2Vp7dIiIi/2MCRO5Rq9tX7g8uxqwLV9tUUaulba7aNAkPl2e3iIjI/5gAkXt0OsBgcFxmMEjlHU2ng3AS0xW9AQfKdTZh6XTAgXIdruidt9lcKDXQ64Hycnl2i4iI/I8JELlHqwXWrrVPggwGabscz4xrtVCtXWuXBF3RG3By2VqkTrZ9FF6rBVIna3Fy2Vq7JOiK3oBDC9fimVVa6PXAsmXA5Ml8FJ6IqLNSCSGE3EEEGovFArVaDbPZjKioKLnDCSwt1wFSq6VLJHJnCU3rAJnMaIhQoyZShy6erAOkVsMcrsOxc1p07w5ERnIdICKiYOTJ+ZsLIZJntNrAywy0WoRotQgB0A1AuOvqTf9j3dYHQJ9r/BYhEREFGN4CIyIiIsVhAkRERESKwwSIiIiIFIcJEBERESkOEyAiIiJSHCZAREREpDiyJkC7du3ClClTEB8fD5VKhU2bNtmUz507FyqVyuZzyy23uOz3jTfeQP/+/REeHo4xY8bg66+/9tMeEBERUTCSdR2gmpoapKSk4L777sOMGTMc1rnllluwbt066/ewsLA2+/zggw+wePFirFmzBmPGjMGqVaswadIkFBcXIyYmxqfxUwBpWqCxqkpaxbCuTvpbowFiYprXLmqxkGNDpAaVYTE4WamFRuO8iTdhmM1SP1FRgMUCVFdLYVVVNZf16SO9cb51G2/HJiIiD4gAAUBs3LjRZtucOXPEtGnTPOpn9OjRIisry/q9oaFBxMfHi+zsbLf7MJvNAoAwm80ejU0yKS0VwmAQIiJCiM2bhdDrhQCaPwaDEGVlzfValNXrDaLyUKm49Vb7JqWl3oXRsp/0dCG2bBHi8GHp79Zlx44JMXNm+8cmIiLPzt8BPwcoPz8fMTExGDx4MObPn4/z5887rXv58mXs27cP6enp1m1dunRBeno6du/e3RHhUkczmYDMTCA3F1i0CHj1VSAvz7ZObi6wZUtzvRa65uWi56OZuGm4ya5JZqbUvadhtLRjB1BWBjz6qPR367L584F77rEP15OxiYjIcwGdAN1yyy3429/+hry8PPzxj39EQUEBJk+ejIaGBof1z507h4aGBuhavcJbp9OhoqLC6Th1dXWwWCw2HwoSRmNz1pGaap/8NImLs89OfhayIxdTU41223Nzpe49DcPR0Nu3Oy7bsQOIj7ff7snYRETkuYB+F1hGRob172HDhiE5ORlXX3018vPzodfrfTZOdnY2nn32WZ/1Rx3IbG7+u7bWeb22ygCE1Zodbjc73uxRPRdDw1m+7e7YRETkuYC+AtTagAED0KdPHxw/ftxheZ8+fRASEgJjq/90NhqNiI2Nddrv0qVLYTabrZ+ysjKfxk1+pFY3/x3exmtQ2yoDUBeudrhd7XizR/VcDA1nLyx2d2wiIvJcUCVAp06dwvnz5xEXF+ewPDQ0FCNHjkRei9sgjY2NyMvLQ1pamtN+w8LCEBUVZfOhIKHTAQaD9HdhIeDsymB5eXO9VhrSDdhcqLPbbjBI3XsahqOhJ050XJaeDpw+bb/dk7GJiMhzsiZA1dXVKCoqQlFREQDgxIkTKCoqQmlpKaqrq/G73/0OhYWFOHnyJPLy8jBt2jQMHDgQkyZNsvah1+vx+uuvW78vXrwYb7/9Nt577z18//33mD9/PmpqavCb3/ymo3ePOoJWC6xdK2UMq1YBCxfaJ0EGAzB5cnO9Fq7oDahetRY7i7R2Tdaudf9x9JZhtJSeDiQkAK+8Iv3dumzNGuD99+3D9WRsIiLynEoIIeQaPD8/HzfddJPd9jlz5mD16tWYPn069u/fj8rKSsTHx8NgMOC5556zmeTcv39/zJ07F8uXL7due/3117Fy5UpUVFRg+PDheO211zBmzBi347JYLFCr1TCbzbwaFCycrQOkVkuXUhyuA6RGZZjObh2g1k28CcNslvpRqx2vA6RWA9HR9usAtWdsIiKl8+T8LWsCFKiYABEREQUfT87fQTUHiIiIiMgXmAARERGR4jABIiIiIsVhAkRERESKwwSIiIiIFIcJEBERESkOEyAiIiJSnIB+GSqR25wthKjRADEx9isLmkzSOyhMJqBnT6BHD+mlXG28M87l2Gaz3Xi15SZ0OWuEMJuh0mjQ2CcG4XFad5q2d2giImoDrwBR8CsrAzIygOuvByoqgAceAIYOBVJTgaQkqazlC26b6l93HTBuHDBiBPDQQ8DRo8CPP3o39pAhduNd/vEUus3OQGjKEISNT0VochK6zc7A5ZIyV03bOzQREbnAlaAd4ErQQcRkks76ubnAk09KL0Rt8TJcK4MByMmR/m6q35peD9x1F/DrX0vvqPBk7FaEwQBxxx3ocv/9dmUN6QZcejcHt2dqHYbRFGpbV3LaGNqt9kREnRFXgiblMBqbs4DUVMfJDyDVMRpt67eWlyclPmfPej52K6rcXHRxkkSF7MhFN5PRaRhNoXo5tFvtiYiUjnOAKLiZzc1/19a6X9eZ2lr36rnTXxvxiMq227rqur3lRERKxwSIgpta3fx3eLj7dZ0JD3evnjv9tRGPStN2W1ddt7eciEjpeAuMgptOJ016AaT5P3q943oGg1S3Zf3W9HqgvByIjvZ87FaEwYDG8nKHZQ3pBtRrdU7DaArVy6Hdak9EpHRMgCi4abXA2rXSWX/VKmDhQvskyGCQ6mi1tvVb0uuBZcuA9HT3JkC3HrvVeKq1a3FFPxkN6bZlDekGNPx1LXomaJ01tYbq5dButSciUjo+BeYAnwILQs7WAVKrpcshrtYB6t5dqtvedYBajWezDpBajcZondN1gJyF6uXQRESK48n5mwmQA0yAiIiIgg8fgyciIiJqAxMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLF4aswSHmaFs8xmYCoKCAsDLh0SVpMJzJSWheoVy+XC+q0XINHowFiYgAtHG2078dkAs6cARobgd5dTIi8ZERItRkqjQYNvWNs1gryNYdxc+0gIlIYXgEiZSkrAzIygCFDgIkTgZMngYceApKTgXHjgOHDgQceAA4cAE6dcqub1FRg5Ejg/IEyiJYbk5KkSmVldm3nzgVKSoBuxjL0zspA9xFDEDouFd2GJSF0bgYu/1jmeGAf7n4bIRIRdXpMgEg5TCYgMxPIzZW+L1okvT5j+3bbenl5wPPPA1u2SG1cdAMAzy4yIfH5TKhabgSkSpmZ1n6a2g4bBtScktqE5Nm26ZKbi64PZqK23H7s9nAUt4MQiYgUgQkQKYfRaHv2T00FduxwXDcvT3onmNHoshsAmJpqtEtkrHJzrf00tU1NBYbHO2/TZXsuupy1H7s9HMXtIEQiIkXgHCBSDrPZ9nttbdv1a2vt2zjoBgDCah1sdNCoqW1tLRCGttsIRwO1g6vufDwcEVFAYwJEyqFW234PD2+7fni4fRsH3QBAXbiDjQ4aNbUNDwfq0HYblaOB2sFVdz4ejogooPEWGCmHTgcYDM3fCwulidCO6PVAebnUxkU3ALC5UIcreoNdXQBS5Z/7aWpbWAgUnXbepnGiAY3R9mO3h6O4HYRIRKQITIBIObRaYO3a5ixg1Spg4UL7rECvB5YtAyZPdvh8eOtuAOCZVVqcXLYWonVfBoNU+ed+mtoeOgRE9JXatE6CGg0GXFmz1uePwjuK20GIRESKoBJCCLkG37VrF1auXIl9+/ahvLwcGzduxPTp0wEA9fX1WLZsGT7//HP8+OOPUKvVSE9Px4oVKxAfH++0z+XLl+PZZ5+12TZ48GAcPXrU7bgsFgvUajXMZjOioqK82jcKYC3XAYqMlO5HNa0D1LOntM3DdYDUaukKis06QNaN7q8DBLUajX10HbYOUBshEhEFHU/O37LOAaqpqUFKSgruu+8+zJgxw6bs4sWL+Pbbb/HUU08hJSUFJpMJCxcuxNSpU7F37942+x06dCh2tHi6p2tXTnWiFrRan5zxHXfjXt+2bbU/fzqGj3afiCioyZoZTJ48GZMnT3ZYplarsb3V+iyvv/46Ro8ejdLSUvTr189pv127dkVsbKxPYyUiIqLOI6jmAJnNZqhUKmg0mjbrHTt2DPHx8RgwYADuvvtulJaWdkyAREREFBSC5t5QbW0tnnjiCcyaNavN+3pjxozB+vXrMXjwYJSXl+PZZ5/FuHHjcPjwYURGRjpsU1dXh7q6Out3i8Xi8/iJiIgocARFAlRfX4+ZM2dCCIHVq1e3WbflLbXk5GSMGTMGV111Ff75z39i3rx5DttkZ2fbTZwmIiKizivgb4E1JT8//fQTtm/f7vFTWRqNBoMGDcLx48ed1lm6dCnMZrP1U8Y3QxIREXVqAZ0ANSU/x44dw44dO9C7d2+P+6iurkZJSQni4uKc1gkLC0NUVJTNh4iIiDovWROg6upqFBUVoaioCABw4sQJFBUVobS0FPX19bjjjjuwd+9e/P3vf0dDQwMqKipQUVGBy5cvW/vQ6/V4/fXXrd+XLFmCgoICnDx5Ev/+979x2223ISQkBLNmzero3SMiIqIAJescoL179+Kmm26yfl+8eDEAYM6cOVi+fDk2b94MABg+fLhNu507d2LChAkAgJKSEpw7d85adurUKcyaNQvnz59HdHQ0brzxRhQWFiI6Otq/O0PBreXqgBoNEBPjerEcd9o01ampkRZXrKpqrh8dDYSFuTVudZkJ3UxGiEozVBoN6rUx6Jmg9Wnc3nRFRBS0BNkxm80CgDCbzXKHQh2htFQIg0EIoPljMEjb29OmqU5MjBCHDgmRnt5cNyJCiJ073Rq39lipuJJuW+9KukE0HDsuxNSpPom77nipx10REQUaT87fTIAcYAKkIBcu2CcELTOACxe8a9OyzubNtskPIMSTTwqh17sct6r0gl3y0/RpTE+X+vFB3FfSDeJPT15wuysiokDkyfk7oCdBE/md0Qjk5jouy82Vyr1p07JOfDzQ4tUsAIDUVCAvz+W43UxGhOxwPJZqxw6pHx/EHbIjF1NT7ds464qIKNgFxTpARH5jNnte7mkbRwtr1ta61YeodDGWs368iDus1nG5q90lIgpGTIBI2dRqz8s9beNoWYXwcLf6UGlcjOWsHy/irgt3XO5qd4mIghFvgZGy6XSAweC4zGCQyr1p07LO6dNAerptvcJCQK93OW69VoeGdMdjifR0qR8fxN2QbsDmQvs2zroiIgp6HTAnKehwErTCBO1TYCW+ewqshE+BEVHw8+T8rRJCCLmTsEBjsVigVqthNpu5KrRStFwER62WLnt4sp6OszbO1gFSq6WFdlquA9TGuLbrAKlRr9XZrwPUzri96YqIKJB4cv5mAuQAEyAiIqLg48n5m3OAiIiISHGYABEREZHiMAEiIiIixWECRERERIrDBIiIiIgUhwkQERERKQ4TICIiIlIcvguMSA5Nqw6aTEDPns3be/SQFkx0tqCi2QxoNNIiir5epfDnMRpNZjREamDqFoPj57XQaIC4OKmKOyGUlwNnzzbX69NHat9yF7RaIDJS+tufu0RE5AyvABF1tLIyICMDGDIEuOEGIDkZePRR4ORJYOFC4MAB4NQpx/VTU4GkJOl7WZlfYupyQyq6DUtCr6wMdD9fhtGjpaKDB4Hrr287hJISYPZsICUFGD9e2rXZs4Hjx4EHH5R2Qa8Hioul7f7cJSKitnAlaAe4EjT5jckknelzc+3L9HopGygsBO66C7jjDmm7s/oGA5CT0/7LJm3EdEVvwKupOVjyghYTJwKjRwMvvOA4hPJyKanZscN+iPR04JFHgKlTgSeflHYxL89/u0REysSVoIkCldHoOJkBpIwgNVX6Z1ycVLet+rm5UrkfY+qal4upqdIY27dL4TkL4exZx8kPIG2Pj5f+btpFR3y1S0RErnAOEFFHMpvbLq+tbf6nq7ru9OeDmMJqm8ubwnPU3FUoFovzPjwIh4jIJ5gAEXUktbrt8vDw5n+6qutOfz6IqS68ubwpPEfNXYXSdDXaUR8ehENE5BO8BUbUkXQ6aaKLI3q9NDlGr5cm1Oh0bdc3GKRyP8Z0RW/A5kJpjIkTpfCchRAdLc31cSQ9HTh9Wvq7aRcd8dUuERG5JMiO2WwWAITZbJY7FOqMSkuFMBiEAJo/er0QmzcLceutQuzcKURZWdv1DQZpux9jqtcbxLebS0VEhFRUUCBERETbIRw/LkR6um2o6enS9pkzpe8REdKutq7n610iIuXx5PzNp8Ac4FNg5Hct1wGKiAC6dAEaG6W/Xa0DpFZLl0n8uA7QlZ5qVIbqUHJBC7Xafh2gtkJouQ6QWi1dGWq9DpBGI90Sa1oHyF+7RETK4sn5mwmQA0yAiIiIgg8fgyciIiJqAxMgIiIiUhwmQERERKQ4TICIiIhIcZgAERERkeIwASIiIiLFkTUB2rVrF6ZMmYL4+HioVCps2rTJplwIgaeffhpxcXHo3r070tPTcezYMZf9vvHGG+jfvz/Cw8MxZswYfP31137aAyIiIgpGsr4LrKamBikpKbjvvvswY8YMu/KXXnoJr732Gt577z0kJibiqaeewqRJk3DkyBGEO3mh0AcffIDFixdjzZo1GDNmDFatWoVJkyahuLgYMTEx/t4louBnMgEXLgBVVUB1tbQwY9NKiGfOoPFKA1QNV9BoqUGDuhca+8QgPK7FCoYtV0LUaIA+fYC4OFwpK0eI6SxQaQa0GjRGaRBSexGorJTqxcS0eyXE2nITupw1QpjN6KqJhAgLQ8PZC1BFRQEaNbrWVkNVUw1UVQMaDerVvRHaL85he5VGY79vTsZyVRcALv9Ujq7ms1BVmiE0GlxR90HoVXEO63rat7vHpK2+nNXzZSy+EogxtVZbbkK3c+VQVV6AKjISjT16oj6qF8JjAytOOQTM7+ffRandB0Bs3LjR+r2xsVHExsaKlStXWrdVVlaKsLAw8Y9//MNpP6NHjxZZWVnW7w0NDSI+Pl5kZ2e7HQtfhUGKVVYmvYpDr7d/T8WuXUJ89pld2ZV0g6g7/vM7LBy9C+PXvxaNx46JxlbbG9PTpXdiNL1fo53vwqg7XiqupDt5xUhMjGjcuVOIiRNtYzAYRMOxEqftbfbNxVjO6gohRMOx4w73v+HY8Xb37ekxcdSXs3oNx46LK7+e6pNYfMWXx8df6kpKReNE+38XG3fuFHU/lrls35n5+/fz5PwdsAlQSUmJACD2799vU2/8+PHikUcecdhHXV2dCAkJselHCCFmz54tpk6d6nYsTIBIkS5cEOKtt+yTn5ZJ0F//6rDsSrpBXDlZap/8AEI8+aTdyd8mQXnySdsxLlzwOPRLpy/YJz+tkyAn+9U4caKoP+kgeWqxb5dOX3BrrNZ1hRCi7uRpp/vfmJ4u6k6e9rpvb49Jy77aqteYnm77+3gZi6/48vj4M8aG1slPi38XG956KyDilENH/H6enL8DdhJ0RUUFAEDX6tXQOp3OWtbauXPn0NDQ4FEbAKirq4PFYrH5ECmO0Sjd6srLc1yemwvExzssCtmRiy5mE7Bjh31haipUjrYD0lipqbZjGI0eBg50OWtEyI5c52PExzvdL9X27Qgxm5y2D9mRiy5nm2Nqa6zWdQFIt72c7L9qxw50NZ/1uu+2uNtXW/VUO3bY/j5exuIrvjw+/tLlrBFdtjv/d7FLXFxAxCmHQPv9AjYB6kjZ2dlQq9XWT0JCgtwhEXU8sxmorW27TlvllZWet3FUbja3Xd8B4aqNq/+ocRa7g/5djdW6XFXZdv2W5Z723Z66TeUu+3Ty+3kSi6/48vj4izvHMxDilEOg/X4BmwDFxsYCAIyt/mvQaDRay1rr06cPQkJCPGoDAEuXLoXZbLZ+ysrK2hk9URBSqwEnDxdYtVWu0XjexlG5Wt12fQdUrtq4eqmxs9gd9O9qrNblQtN2/ZblnvbdnrpN5S77dPL7eRKLr/jy+PiLO8czEOKUQ6D9fgGbACUmJiI2NhZ5LS5bWywW7NmzB2lpaQ7bhIaGYuTIkTZtGhsbkZeX57QNAISFhSEqKsrmQ6Q4Op30BJde77jcYABOn3ZY1JBuQKNaC6Sn2xcWFkI42g5IYxUW2o7R6ha2OxqjdWhINzgf4/Rpp/slJk5Eg1rrtH1DugGN0c0xtTVW67oAcEUd7XT/RXo6rqijve67Le721VY9kZ5u+/t4GYuv+PL4+EtjtA6NE53/u9hYXh4Qccoh4H6/ds84aoeqqiqxf/9+sX//fgFA/PnPfxb79+8XP/30kxBCiBUrVgiNRiM+/vhjcfDgQTFt2jSRmJgoLl26ZO3j5ptvFn/5y1+s33NyckRYWJhYv369OHLkiLj//vuFRqMRFRUVbsfFSdCkWN4+BVYSJE+BGWzLGw0G0XDcxVNgJR48BeagrhBtPAX289jt6dvTY+KoL+dPgZU4fgrMi1h8xZfHx1/qSkpFo4FPgTni79/Pk/O3SgghOjblapafn4+bbrrJbvucOXOwfv16CCHwzDPP4K233kJlZSVuvPFGvPnmmxg0aJC1bv/+/TF37lwsX77cuu3111/HypUrUVFRgeHDh+O1117DmDFj3I7LYrFArVbDbDbzahApT8t1gGpqpLV5nK0DFKVFY7TO+TpAajUQHW2/DpBGjUa1tnkdILVauvLjr3WAIiMBraZ5HaDqaiBKjXpNH+frAKnV9vvmZCxXdYHW6wCpcUUd7d46QG707e4xaasvZ/V8GYuvBGJMrVnXATKboOoZgYYekbjCdYAA+Pf38+T8LWsCFKiYABEREQUfT87fATsHiIiIiMhfmAARERGR4jABIiIiIsVhAkRERESKwwSIiIiIFIcJEBERESkOEyAiIiJSHCZAREREpDhMgIiIiEhxmAARERGR4nSVO4BA1PR2EIvFInMkRERE5K6m87Y7b/liAuRAVVUVACAhIUHmSIiIiMhTVVVVUKvVbdbhy1AdaGxsxOnTpxEZGQmVSuW0nsViQUJCAsrKyhT/0lQeC1s8Hs14LJrxWNji8WjGY9GsPcdCCIGqqirEx8ejS5e2Z/nwCpADXbp0Qd++fd2uHxUVpfh/YZvwWNji8WjGY9GMx8IWj0czHotm3h4LV1d+mnASNBERESkOEyAiIiJSHCZA7RAWFoZnnnkGYWFhcociOx4LWzwezXgsmvFY2OLxaMZj0ayjjgUnQRMREZHi8AoQERERKQ4TICIiIlIcJkBERESkOEyAvPTGG2+gf//+CA8Px5gxY/D111/LHZIsdu3ahSlTpiA+Ph4qlQqbNm2SOyTZZGdnY9SoUYiMjERMTAymT5+O4uJiucOSzerVq5GcnGxdyyMtLQ1btmyRO6yAsGLFCqhUKixatEjuUDrc8uXLoVKpbD5JSUlyhyWr//znP7jnnnvQu3dvdO/eHcOGDcPevXvlDqvD9e/f3+7fDZVKhaysLL+MxwTICx988AEWL16MZ555Bt9++y1SUlIwadIknDlzRu7QOlxNTQ1SUlLwxhtvyB2K7AoKCpCVlYXCwkJs374d9fX1MBgMqKmpkTs0WfTt2xcrVqzAvn37sHfvXtx8882YNm0avvvuO7lDk9U333yDv/71r0hOTpY7FNkMHToU5eXl1s+XX34pd0iyMZlMGDt2LLp164YtW7bgyJEjePnll6HVauUOrcN98803Nv9ebN++HQBw5513+mdAQR4bPXq0yMrKsn5vaGgQ8fHxIjs7W8ao5AdAbNy4Ue4wAsaZM2cEAFFQUCB3KAFDq9WKtWvXyh2GbKqqqsQ111wjtm/fLn75y1+KhQsXyh1Sh3vmmWdESkqK3GEEjCeeeELceOONcocRkBYuXCiuvvpq0djY6Jf+eQXIQ5cvX8a+ffuQnp5u3dalSxekp6dj9+7dMkZGgcZsNgMAevXqJXMk8mtoaEBOTg5qamqQlpYmdziyycrKwq233mrz/x9KdOzYMcTHx2PAgAG4++67UVpaKndIstm8eTOuv/563HnnnYiJicGIESPw9ttvyx2W7C5fvoz3338f9913X5vv5GwPJkAeOnfuHBoaGqDT6Wy263Q6VFRUyBQVBZrGxkYsWrQIY8eOxXXXXSd3OLI5dOgQevbsibCwMDz44IPYuHEjrr32WrnDkkVOTg6+/fZbZGdnyx2KrMaMGYP169dj69atWL16NU6cOIFx48ahqqpK7tBk8eOPP2L16tW45pprsG3bNsyfPx+PPPII3nvvPblDk9WmTZtQWVmJuXPn+m0MvgyVyA+ysrJw+PBhRc9tAIDBgwejqKgIZrMZ//rXvzBnzhwUFBQoLgkqKyvDwoULsX37doSHh8sdjqwmT55s/Ts5ORljxozBVVddhX/+85+YN2+ejJHJo7GxEddffz1efPFFAMCIESNw+PBhrFmzBnPmzJE5Ovm88847mDx5MuLj4/02Bq8AeahPnz4ICQmB0Wi02W40GhEbGytTVBRIFixYgE8//RQ7d+5E37595Q5HVqGhoRg4cCBGjhyJ7OxspKSk4NVXX5U7rA63b98+nDlzBr/4xS/QtWtXdO3aFQUFBXjttdfQtWtXNDQ0yB2ibDQaDQYNGoTjx4/LHYos4uLi7P6DYMiQIYq+LfjTTz9hx44dyMzM9Os4TIA8FBoaipEjRyIvL8+6rbGxEXl5eYqe20CAEAILFizAxo0b8cUXXyAxMVHukAJOY2Mj6urq5A6jw+n1ehw6dAhFRUXWz/XXX4+7774bRUVFCAkJkTtE2VRXV6OkpARxcXFyhyKLsWPH2i2X8cMPP+Cqq66SKSL5rVu3DjExMbj11lv9Og5vgXlh8eLFmDNnDq6//nqMHj0aq1atQk1NDX7zm9/IHVqHq66utvkvtxMnTqCoqAi9evVCv379ZIys42VlZWHDhg34+OOPERkZaZ0Tplar0b17d5mj63hLly7F5MmT0a9fP1RVVWHDhg3Iz8/Htm3b5A6tw0VGRtrNBYuIiEDv3r0VN0dsyZIlmDJlCq666iqcPn0azzzzDEJCQjBr1iy5Q5PFo48+ihtuuAEvvvgiZs6cia+//hpvvfUW3nrrLblDk0VjYyPWrVuHOXPmoGtXP6cofnm2TAH+8pe/iH79+onQ0FAxevRoUVhYKHdIsti5c6cAYPeZM2eO3KF1OEfHAYBYt26d3KHJ4r777hNXXXWVCA0NFdHR0UKv14vc3Fy5wwoYSn0M/q677hJxcXEiNDRU/Nd//Ze46667xPHjx+UOS1affPKJuO6660RYWJhISkoSb731ltwhyWbbtm0CgCguLvb7WHwbPBERESkO5wARERGR4jABIiIiIsVhAkRERESKwwSIiIiIFIcJEBERESkOEyAiIiJSHCZAREREpDhMgIiIiEhxmAARkeJNmDABixYt8rq9SqXCpk2bfBYPEfkf3wVGRF6bMGEChg8fjlWrVskdiqzKy8uh1WrlDoOIPMAEiIionWJjY+UOgYg8xFtgROSVuXPnoqCgAK+++ipUKhVUKhVOnjyJw4cPY/LkyejZsyd0Oh3uvfdenDt3ztpuwoQJePjhh7Fo0SJotVrodDq8/fbbqKmpwW9+8xtERkZi4MCB2LJli7VNfn4+VCoVPvvsMyQnJyM8PBypqak4fPiw2/F+9dVXmDBhAnr06AGtVotJkybBZDJZyxsbG/H444+jV69eiI2NxfLly93uu+UtsJMnT0KlUiEnJwc33HADwsPDcd1116GgoMCmzXfffYdf//rXiIqKQmRkJMaNG4eSkhJr+bvvvouhQ4ciLCwMcXFxWLBggdvxEJFrTICIyCuvvvoq0tLS8Nvf/hbl5eUoLy9HZGQkbr75ZowYMQJ79+7F1q1bYTQaMXPmTJu27733Hvr06YOvv/4aDz/8MObPn48777wTN9xwA7799lsYDAbce++9uHjxok273/3ud3j55ZfxzTffIDo6GlOmTEF9fb3LWIuKiqDX63Httddi9+7d+PLLLzFlyhQ0NDTYxBQREYE9e/bgpZdewh/+8Ads377d6+Pzu9/9Do899hj279+PtLQ0TJkyBefPnwcA/Oc//8H48eMRFhaGL774Avv27cN9992HK1euAABWr16NrKws3H///Th06BA2b96MgQMHeh0LETng9/fNE1Gn9ctf/lIsXLjQ+v25554TBoPBpk5ZWZkAIIqLi61tbrzxRmv5lStXREREhLj33nut28rLywUAsXv3biGEEDt37hQARE5OjrXO+fPnRffu3cUHH3zgMs5Zs2aJsWPHtrkfLWMSQohRo0aJJ554wmXfQggBQGzcuFEIIcSJEycEALFixQpreX19vejbt6/44x//KIQQYunSpSIxMVFcvnzZYX/x8fHiySefdGtsIvIO5wARkc8cOHAAO3fuRM+ePe3KSkpKMGjQIABAcnKydXtISAh69+6NYcOGWbfpdDoAwJkzZ2z6SEtLs/7dq1cvDB48GN9//73LuIqKinDnnXe2WadlTAAQFxdnN74nWsbatWtXXH/99dZYi4qKMG7cOHTr1s2u3ZkzZ3D69Gno9XqvxyYi15gAEZHPVFdXY8qUKfjjH/9oVxYXF2f9u/WJX6VS2WxTqVQApHk5vtC9e3eXdRzF5KvxPYnHnViJqP04B4iIvBYaGmozj+YXv/gFvvvuO/Tv3x8DBw60+URERLR7vMLCQuvfJpMJP/zwA4YMGeKyXXJyMvLy8to9vidaxnrlyhXs27fPGmtycjL+7//+z+H8pcjISPTv37/D4yVSGiZAROS1/v37Y8+ePTh58iTOnTuHrKwsXLhwAbNmzcI333yDkpISbNu2Db/5zW9sEiVv/eEPf0BeXh4OHz6MuXPnok+fPpg+fbrLdkuXLsU333yDhx56CAcPHsTRo0exevVqm6fTfO2NN97Axo0bcfToUWRlZcFkMuG+++4DACxYsAAWiwUZGRnYu3cvjh07hv/5n/9BcXExAGD58uV4+eWX8dprr+HYsWP49ttv8Ze//MVvsRIpERMgIvLakiVLEBISgmuvvRbR0dG4fPkyvvrqKzQ0NMBgMGDYsGFYtGgRNBoNunRp///drFixAgsXLsTIkSNRUVGBTz75BKGhoS7bDRo0CLm5uThw4ABGjx6NtLQ0fPzxx+ja1X+zAFasWIEVK1YgJSUFX375JTZv3ow+ffoAAHr37o0vvvgC1dXV+OUvf4mRI0fi7bfftt6GmzNnDlatWoU333wTQ4cOxa9//WscO3bMb7ESKZFKCCHkDoKIqC35+fm46aabYDKZoNFo5A6nTSdPnkRiYiL279+P4cOHyx0OETnBK0BERESkOEyAiCjoNa087ejz4osvtqvvv//97077Hjp0qI/2gIg6Gm+BEVHQ+89//oNLly45LOvVqxd69erldd9VVVUwGo0Oy7p164arrrrK676JSD5MgIiIiEhxeAuMiIiIFIcJEBERESkOEyAiIiJSHCZAREREpDhMgIiIiEhxmAARERGR4jABIiIiIsVhAkRERESK8/8BweqiCUeUky8AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rgi_reg='11'\n", "m = 'PyGEM-OGGM_v13'\n", "perc=-50\n", "_r_no_shift = resp_time_estimate(model_author=m, perc_change_l=[perc],\n", " rgi_reg=rgi_reg, min_perc_change=25, shift_years=False)\n", "_r_w_shift = resp_time_estimate(model_author=m, perc_change_l=[perc],\n", " rgi_reg=rgi_reg, min_perc_change=25, shift_years=True)\n", "sns.scatterplot(data=_r_no_shift.dropna(),\n", " x='temp_ch_ipcc', y=f'resp_time_{perc}%', color='blue')\n", "sns.scatterplot(data=_r_w_shift.dropna(),\n", " x='temp_ch_ipcc', y=f'resp_time_{perc}%', color='red')" ] }, { "cell_type": "markdown", "id": "2dde31b0-55b2-43d6-b3a0-3181773b061c", "metadata": {}, "source": [ "- it makes sense that the response timescale is shorter for 2020 as the glaciers retreated to higher altitudes, where they are likely steeper and thus respond in general faster ... " ] }, { "cell_type": "markdown", "id": "04bde477-6b6a-4736-8283-ab963ad5fc43", "metadata": {}, "source": [ "some previous tests showed (removed core here to not have too much code):\n", "- only when we set a threshold of how much volume loss has to occur, we get response timescales that are not totally noisy \n", " - much less noise, when choosing min_perc_change=25, instead of min_perc_change=10, however,it also results in more models where no response timescale can be estimated, which is not good. " ] }, { "cell_type": "markdown", "id": "7ccc875f-21ec-4f91-94dd-e829561417cb", "metadata": {}, "source": [ "**test: check that the correct volume loss occurs at the computed response timescale**\n" ] }, { "cell_type": "code", "execution_count": 41, "id": "f982829d-709a-44e9-a791-f59af62221e3", "metadata": {}, "outputs": [], "source": [ "rgi_regs = []\n", "for rgi_reg in np.arange(1,20,1):\n", " if rgi_reg < 10:\n", " rgi_reg = '0'+str(rgi_reg)\n", " else:\n", " rgi_reg = str(rgi_reg)\n", " rgi_regs.append(rgi_reg)\n", "\n", "for model_author in ['OGGM_v16', 'GLIMB']: ## just select two random models to test\n", " for rgi_reg in rgi_regs:\n", " pd_test = resp_time_estimate(model_author=model_author, perc_change_l=[-50],\n", " roll_volume=21,\n", " rgi_reg=rgi_reg, min_perc_change=25)\n", " pd_test = pd_test.dropna()\n", " pd_test = pd_test.loc[pd_test.temp_ch_ipcc>1.2]\n", "\n", " for j,ind in enumerate(pd_test.index):\n", " ds_test = ds_reg_models_extend.volume_m3.sel(model_author=model_author, rgi_reg=rgi_reg).stack(experiments=['gcm','period_scenario']).sel(experiments=pd_test.index[j])\n", " ds_init = ds_test[0].values\n", " ds_test = ds_test.rolling(simulation_year=21, center=True).mean()\n", " response_timescale = pd_test.iloc[j]['resp_time_-50%'] \n", " # test works only if the response timescale is long enough, as otherwise issues due to fast changes and rolling averages\n", " if response_timescale>20:\n", " vol_change_at_resp_time_yr = ds_init-ds_test.sel(simulation_year = response_timescale).values\n", " \n", " # total volume change\n", " vol_change_total = ds_init-ds_test[-11]\n", " # is the volume change at response timescale year 50% of the total change? \n", " # it does not have to match perfectly as a lot happens in between specifically for small response timescales\n", " np.testing.assert_allclose(vol_change_at_resp_time_yr, 0.5*vol_change_total, rtol=0.05)\n", " " ] }, { "cell_type": "code", "execution_count": 42, "id": "5cade48b-4f39-4976-a33b-0208c8821891", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(rgi_reg RGI06\n", " model_author GO\n", " min_perc_change 25\n", " temp_ch_ipcc 0.901467\n", " category shrink\n", " shift_years False\n", " resp_time_-50% 10.0\n", " dtype: object,\n", " rgi_reg RGI06\n", " model_author GO\n", " min_perc_change 25\n", " temp_ch_ipcc 6.884361\n", " category shrink\n", " shift_years False\n", " resp_time_-50% 23.0\n", " dtype: object)" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd_test = resp_time_estimate(model_author='GO', perc_change_l=[-50],\n", " rgi_reg='06', min_perc_change=25)\n", "pd_test.dropna().min(), pd_test.dropna().max()" ] }, { "cell_type": "code", "execution_count": null, "id": "1cf05158-6560-42a9-bbe5-989a201954fa", "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 43, "id": "54183372-84c0-4b30-9a0f-6f97f8a70cd5", "metadata": {}, "outputs": [], "source": [ "shift_years = True\n", "rgi_regs_global = ['All', '01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12',\n", " '13', '14', '15', '16', '17', '18', '19']\n", "\n", "if shift_years:\n", " # we shift by maximum by +50 years, so like that there should always be values inside \n", " _ds_reg_models_extend = ds_reg_yr_shift.sel(year_after_2020=slice(0,4950))\n", " xx = 'year_after_2020'\n", " pd_reg_models_count = _ds_reg_models_extend.volume_m3.sel(year_after_2020=0).to_dataframe().groupby(['rgi_reg', 'model_author'])['volume_m3'].count().reset_index()\n", " p_add = '_shifted'\n", " lab_add = 'after 2020 '\n", "\n", "else:\n", " _ds_reg_models_extend = ds_reg_models_extend\n", " xx = 'simulation_year'\n", " pd_reg_models_count = _ds_reg_models_extend.volume_m3.sel(simulation_year=0).to_dataframe().groupby(['rgi_reg', 'model_author'])['volume_m3'].count().reset_index()\n", " p_add = ''\n", " lab_add = ''\n", "\n", "## rolling mean to estimate the response time scale.. \n", "roll_volume=21 \n", "# assume that at least 25% change is necessary to estimate a \"response timescale\"\n", "min_perc_change = 25 # also tested: 1,5,10,20," ] }, { "cell_type": "code", "execution_count": 29, "id": "cefc1ccd-1524-49be-b668-0bfa78113236", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rgi_regmodel_authorvolume_m3
001CISM20
101GLIMB80
201GO0
301GloGEMflow80
401GloGEMflow3D0
............
14719GloGEMflow80
14819GloGEMflow3D0
14919Kraaijenbrink0
15019OGGM_v1680
15119PyGEM-OGGM_v1380
\n", "

152 rows × 3 columns

\n", "
" ], "text/plain": [ " rgi_reg model_author volume_m3\n", "0 01 CISM2 0\n", "1 01 GLIMB 80\n", "2 01 GO 0\n", "3 01 GloGEMflow 80\n", "4 01 GloGEMflow3D 0\n", ".. ... ... ...\n", "147 19 GloGEMflow 80\n", "148 19 GloGEMflow3D 0\n", "149 19 Kraaijenbrink 0\n", "150 19 OGGM_v16 80\n", "151 19 PyGEM-OGGM_v13 80\n", "\n", "[152 rows x 3 columns]" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd_reg_models_count" ] }, { "cell_type": "code", "execution_count": 44, "id": "c5a75de0-6250-42ff-a947-7dfe56559863", "metadata": {}, "outputs": [], "source": [ "\n", "pd_reg_models_count = pd_reg_models_count.loc[pd_reg_models_count.rgi_reg !='All']\n", "pd_reg_models_count = pd_reg_models_count.rename(columns={'volume_m3':'n_experiments'})\n", "# models available for each region...\n", "_regional_model = pd_reg_models_count.loc[pd_reg_models_count.n_experiments > 0 ] # ==80\n", "\n", "## amount of models per region\n", "# at least 1 experiment for that RGI region\n", "reg_n_models = pd_reg_models_count.loc[pd_reg_models_count.n_experiments > 0].groupby('rgi_reg').count()['model_author'] \n", "# models that are globally available \n", "max_exps=pd_reg_models_count.n_experiments.max()\n", "global_models = pd_reg_models_count.groupby('model_author').sum().where(pd_reg_models_count.groupby('model_author').sum().n_experiments == 19 * max_exps).dropna().index\n", "\n", "pd_category_resp_time_l = []\n", "for perc in [-50,-80,-90]: \n", " for j, rgi_reg in enumerate(rgi_regs_global):\n", " if rgi_reg == 'All':\n", " # only take models that exist globally\n", " # at the moment Rounce, Glimb, OGGM,GloGEMflow\n", " n = len(global_models)\n", " mj = 0\n", " for m in glac_models:\n", " if m in global_models:\n", " mj+=1\n", " pd_category_resp_time = resp_time_estimate(model_author=m, perc_change_l=[perc],\n", " rgi_reg=rgi_reg, min_perc_change=min_perc_change,\n", " roll_volume=roll_volume,shift_years=shift_years)\n", " pd_category_resp_time_l.append(pd_category_resp_time)\n", " else:\n", " pass\n", "\n", " else: \n", " models_available = _regional_model.loc[_regional_model.rgi_reg==rgi_reg]['model_author'].values\n", " n = reg_n_models[rgi_reg]\n", " mj=0\n", " for m in glac_models:\n", " if m in models_available:\n", " mj+=1\n", " \n", " pd_category_resp_time = resp_time_estimate(model_author=m, perc_change_l=[perc],\n", " rgi_reg=rgi_reg, min_perc_change=min_perc_change,\n", " roll_volume=roll_volume,shift_years=shift_years)\n", " pd_category_resp_time_l.append(pd_category_resp_time)\n", "\n", "pd_response_time_x_perc_loss = pd.concat(pd_category_resp_time_l)\n", "# just some reformatting (no real averaging...)\n", "pd_response_time_x_perc_loss_comp = pd_response_time_x_perc_loss.copy().drop(columns=['category','shift_years']) # loc[pd_response_time_x_perc_loss.rgi_reg !='All']\n", "pd_response_time_x_perc_loss_comp = pd_response_time_x_perc_loss_comp.reset_index(names=['gcm_period_scenario'])\n", "# Convert the string representations of tuples into actual tuples\n", "tuple_list = [eval(str(item)) for item in list(pd_response_time_x_perc_loss_comp['gcm_period_scenario'].values)]\n", "# Create a new list in 'gcm_period_scenario' format\n", "reformatted_list = [f\"{gcm}_{period_scenario}\" for gcm, period_scenario in tuple_list]\n", "pd_response_time_x_perc_loss_comp['gcm_period_scenario'] = reformatted_list\n", "pd_response_time_x_perc_loss_comp = pd_response_time_x_perc_loss_comp.groupby(['rgi_reg', 'model_author', 'temp_ch_ipcc','gcm_period_scenario']).mean()\n", "ds_response_time_x_perc_loss_comp = pd_response_time_x_perc_loss_comp.to_xarray()\n", "\n", "# this is the difference to the median glacier model \n", "for p in [-50,-80,-90]:\n", " ds_response_time_x_perc_loss_comp[f'diff_resp_time_{p}%'] = ds_response_time_x_perc_loss_comp[f'resp_time_{p}%'] - ds_response_time_x_perc_loss_comp[f'resp_time_{p}%'].median(dim='model_author') \n", "\n", "pd_response_time_x_perc_loss_comp = ds_response_time_x_perc_loss_comp.to_dataframe().reset_index()\n", "pd_response_time_x_perc_loss_comp = pd_response_time_x_perc_loss_comp.dropna()\n", "\n", "# create \"classes\" of global temp. change <--- not anymore used ... \n", "#T0_l, T1_l = [], []\n", "#temp_labels = []\n", "#for t0, t1 in zip([min_t_r, 2, 4],[2, 4, max_t_r]):\n", "# temp_label = f'{t0:.1f}°C-{t1:.1f}°C'\n", "# condi = (pd_response_time_x_perc_loss_comp.temp_ch_ipcc <=t1)&(pd_response_time_x_perc_loss_comp.temp_ch_ipcc >=t0)\n", "# temp_labels_length_l = len(pd_response_time_x_perc_loss_comp.loc[condi].groupby('temp_ch_ipcc').count())\n", "# pd_response_time_x_perc_loss_comp.loc[condi, 'temp_ch_classes'] = temp_label + f' (n={temp_labels_length_l})'\n", "# for the summary plot using the K-means clustering in 3_fitted_glacier_response.... \n", "\n", "pd_response_time_x_perc_loss_comp['min_perc_change'] =min_perc_change\n", "pd_response_time_x_perc_loss_comp.reset_index(drop=True).to_csv(f'../data/resp_time{p_add}_for_deltaT_rgi_reg_roll_volume_{roll_volume}yravg.csv')" ] }, { "cell_type": "code", "execution_count": 49, "id": "48b5d77e-c909-4918-add5-c5276a379785", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
rgi_regmodel_authortemp_ch_ipccgcm_period_scenariomin_perc_changeresp_time_-50%resp_time_-80%resp_time_-90%diff_resp_time_-50%diff_resp_time_-80%diff_resp_time_-90%
0AllGLIMB0.302892ipsl-cm6a-lr_1951-1970_hist25420.01138.01785.088.0208.0282.5
1AllGLIMB0.478289gfdl-esm4_1901-1920_hist25382.01152.01877.087.0260.5397.5
2AllGLIMB0.901467gfdl-esm4_1995-2014_hist25312.01001.01660.0122.0440.0720.0
3AllGLIMB0.904940mpi-esm1-2-hr_1995-2014_hist25329.01032.01763.0148.0511.0888.0
4AllGLIMB0.913313ipsl-cm6a-lr_1995-2014_hist25307.0966.01597.0135.0483.0781.0
....................................
6010RGI19PyGEM-OGGM_v134.677751ipsl-cm6a-lr_2081-2100_ssp37025123.0275.0374.022.025.0-38.0
6011RGI19PyGEM-OGGM_v135.230543ukesm1-0-ll_2061-2080_ssp58525116.0267.0377.030.040.510.0
6012RGI19PyGEM-OGGM_v135.840495ukesm1-0-ll_2081-2100_ssp3702569.0158.0223.017.530.031.0
6013RGI19PyGEM-OGGM_v135.849305ipsl-cm6a-lr_2081-2100_ssp5852593.0204.0277.018.034.537.5
6014RGI19PyGEM-OGGM_v136.884361ukesm1-0-ll_2081-2100_ssp5852558.0132.0184.015.531.540.0
\n", "

6015 rows × 11 columns

\n", "
" ], "text/plain": [ " rgi_reg model_author temp_ch_ipcc gcm_period_scenario \\\n", "0 All GLIMB 0.302892 ipsl-cm6a-lr_1951-1970_hist \n", "1 All GLIMB 0.478289 gfdl-esm4_1901-1920_hist \n", "2 All GLIMB 0.901467 gfdl-esm4_1995-2014_hist \n", "3 All GLIMB 0.904940 mpi-esm1-2-hr_1995-2014_hist \n", "4 All GLIMB 0.913313 ipsl-cm6a-lr_1995-2014_hist \n", "... ... ... ... ... \n", "6010 RGI19 PyGEM-OGGM_v13 4.677751 ipsl-cm6a-lr_2081-2100_ssp370 \n", "6011 RGI19 PyGEM-OGGM_v13 5.230543 ukesm1-0-ll_2061-2080_ssp585 \n", "6012 RGI19 PyGEM-OGGM_v13 5.840495 ukesm1-0-ll_2081-2100_ssp370 \n", "6013 RGI19 PyGEM-OGGM_v13 5.849305 ipsl-cm6a-lr_2081-2100_ssp585 \n", "6014 RGI19 PyGEM-OGGM_v13 6.884361 ukesm1-0-ll_2081-2100_ssp585 \n", "\n", " min_perc_change resp_time_-50% resp_time_-80% resp_time_-90% \\\n", "0 25 420.0 1138.0 1785.0 \n", "1 25 382.0 1152.0 1877.0 \n", "2 25 312.0 1001.0 1660.0 \n", "3 25 329.0 1032.0 1763.0 \n", "4 25 307.0 966.0 1597.0 \n", "... ... ... ... ... \n", "6010 25 123.0 275.0 374.0 \n", "6011 25 116.0 267.0 377.0 \n", "6012 25 69.0 158.0 223.0 \n", "6013 25 93.0 204.0 277.0 \n", "6014 25 58.0 132.0 184.0 \n", "\n", " diff_resp_time_-50% diff_resp_time_-80% diff_resp_time_-90% \n", "0 88.0 208.0 282.5 \n", "1 87.0 260.5 397.5 \n", "2 122.0 440.0 720.0 \n", "3 148.0 511.0 888.0 \n", "4 135.0 483.0 781.0 \n", "... ... ... ... \n", "6010 22.0 25.0 -38.0 \n", "6011 30.0 40.5 10.0 \n", "6012 17.5 30.0 31.0 \n", "6013 18.0 34.5 37.5 \n", "6014 15.5 31.5 40.0 \n", "\n", "[6015 rows x 11 columns]" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pd_response_time_x_perc_loss_comp.reset_index(drop=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "a09c0856-932a-4d69-8ea5-56e5f3800f5e", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:oggm_gmip3_working]", "language": "python", "name": "conda-env-oggm_gmip3_working-py" }, "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.10.14" } }, "nbformat": 4, "nbformat_minor": 5 }