{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "91465e8a-502c-4738-be38-db94bd6a9866",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load_ext watermark # this is so that in the end we can check which module versions we used\n",
    "%load_ext autoreload\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import glob\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import xesmf as xe\n",
    "\n",
    "import cartopy\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as mticker\n",
    "import matplotlib.gridspec as gridspec \n",
    "\n",
    "from matplotlib import cm, ticker\n",
    "from scipy import stats\n",
    "\n",
    "def linregress_1d(x, y):\n",
    "    # Wrapper around scipy linregress to use in apply_ufunc\n",
    "    slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)\n",
    "    return np.array([slope, intercept, r_value, p_value, std_err])\n",
    "\n",
    "\n",
    "def linregress_xr(da, dim=\"time\"):\n",
    "    NT = len(da[dim])\n",
    "    da[dim] = np.arange(0, NT)\n",
    "    stats = xr.apply_ufunc(\n",
    "        linregress_1d,\n",
    "        da[dim],\n",
    "        da,\n",
    "        input_core_dims=[[dim], [dim]],\n",
    "        output_core_dims=[[\"parameter\"]],\n",
    "        vectorize=True,\n",
    "        dask=\"parallelized\",\n",
    "        output_dtypes=[\"float32\"],\n",
    "        output_sizes={\"parameter\": 5},\n",
    "    )\n",
    "    return (NT - 1) * stats[:, :, 0], stats[:, :, 3]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "751ce0df-5772-4979-ad56-ac9eb69f45fd",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b67e0780-ccb5-4e33-bcaa-74c9e909ac9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "target_grid = xr.Dataset({'lat': (['lat'], np.arange(-30.5, 31, 1)),\n",
    "        'lon': (['lon'], np.arange(100.5,291,1)),})\n",
    "encoding={'TEMP':{'dtype':'float32'}, 'lon':{'dtype':'float32'}, 'lat':{'dtype':'float32'}}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "da131548-d7e8-46df-ac8b-c296810453b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "VAR = \"TS\"\n",
    "\n",
    "## Load CESM\n",
    "path = \"/glade/campaign/univ/uclb0039/jingyiz/cesm2.f19_g17.BHIST_flx_adj/combie_data\"\n",
    "ds_cntl = xr.open_dataset(f\"{path}/b.e21.BHIST.f19_g17.cntl.{VAR}_1950-2050.nc\")\n",
    "\n",
    "## Load CESM-FA\n",
    "path = \"/glade/campaign/univ/uclb0039/jingyiz/cesm2.f19_g17.BHIST_flx_adj/combie_data\"\n",
    "ds_fa = xr.open_dataset(f\"{path}/b.e21.BHIST.f19_g17.fa.{VAR}_1950-2050.nc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2264fd4-beda-4174-965f-0b4307229f96",
   "metadata": {},
   "source": [
    "## Calculate trend for each ensemble member"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3478fa3e-4f64-4c4b-b585-ff057f007bc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "y0 = 1970\n",
    "ye = 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "25340c06-52ba-42fb-8128-4de731249a70",
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_t= []\n",
    "lst_p = []\n",
    "for i in range(15): \n",
    "    ds = ds_cntl.isel(ensem_meber = i)\n",
    "    ds = ds.sel(time=slice(str(y0), str(ye))).resample(time='1Y').mean() #resample to annual mean\n",
    "    ds = ds.compute()\n",
    "    # regrid data\n",
    "    regridder=xe.Regridder(ds[VAR], target_grid, 'bilinear')\n",
    "    ds_rg = regridder(ds)\n",
    "    # calculate the trend\n",
    "    t, p = linregress_xr(ds_rg[VAR])\n",
    "    lst_t.append(t.values) \n",
    "    lst_p.append(p.values)\n",
    "\n",
    "ds_trend_cntl = xr.Dataset()\n",
    "ds_trend_cntl['trend'] = (('ensem_meber','lat','lon'), np.array(lst_t))\n",
    "ds_trend_cntl['pvalue']  = (('ensem_meber','lat','lon'), np.array(lst_p))\n",
    "ds_trend_cntl.coords['lat'] = target_grid.lat\n",
    "ds_trend_cntl.coords['lon'] = target_grid.lon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "00c2b26e-cdc1-4a5d-a615-c7cbd05a3b95",
   "metadata": {},
   "outputs": [],
   "source": [
    "lst_t= []\n",
    "lst_p = []\n",
    "for i in range(15): \n",
    "    ds = ds_fa.isel(ensem_meber = i)\n",
    "    ds = ds.sel(time=slice(str(y0), str(ye))).resample(time='1Y').mean()#resample to annual mean\n",
    "    ds = ds.compute()\n",
    "    # regrid data\n",
    "    regridder=xe.Regridder(ds[VAR], target_grid, 'bilinear')\n",
    "    ds_rg = regridder(ds)\n",
    "    # calculate the trend\n",
    "    t, p = linregress_xr(ds_rg[VAR])\n",
    "    lst_t.append(t.values) \n",
    "    lst_p.append(p.values)\n",
    "\n",
    "ds_trend_fa = xr.Dataset()\n",
    "ds_trend_fa['trend'] = (('ensem_meber','lat','lon'), np.array(lst_t))\n",
    "ds_trend_fa['pvalue']  = (('ensem_meber','lat','lon'), np.array(lst_p))\n",
    "ds_trend_fa.coords['lat'] = target_grid.lat\n",
    "ds_trend_fa.coords['lon'] = target_grid.lon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23105233-a58b-427f-89cc-82129e676466",
   "metadata": {},
   "source": [
    "## Calculate SST metrics of all members; Ref=HadISST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c3ba82a3-e07e-48f7-b3ba-278e7afd09af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PC and RMSE\n",
    "xE = 140; xW = 270;\n",
    "yN = 10; yS = -10;\n",
    "\n",
    "# NS gradient index \n",
    "xNS_W = 190; xNS_E = 240;\n",
    "yNS_N = 10; yNS_Eq = 5; \n",
    "\n",
    "# EW gradient index\n",
    "xEW_W1 = 140; xEW_E1 = 170; \n",
    "xEW_W2 = 200; xEW_E2 = 260; \n",
    "yWE = 5;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e489c2f4-ec09-4965-ae70-b3631c03ad1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load observational reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "24d10e31-1c0f-403d-b83c-4bc9dda6fd75",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NS index 1 <xarray.DataArray 'trend' ()>\n",
      "array(0.55352294, dtype=float32)\n",
      "NS index 2 <xarray.DataArray 'trend' ()>\n",
      "array(0.60592264, dtype=float32)\n"
     ]
    }
   ],
   "source": [
    "# Observation Ref: \n",
    "ds_hadi = xr.open_dataset('/glade/work/jingyiz/data/HadISST/HadISST_sst.nc')\n",
    "ds_hadi = ds_hadi.sel(time=slice(str(y0), str(ye))).resample(time='1Y').mean()\n",
    "regridder=xe.Regridder(ds_hadi['sst'], target_grid, 'bilinear')\n",
    "ds_hadi_rg = regridder(ds_hadi)\n",
    "\n",
    "t, p = linregress_xr(ds_hadi_rg['sst'])\n",
    "dsHad = ds_hadi_rg.copy()#xr.Dataset()\n",
    "dsHad['trend'] = (('lat','lon'), t.values)\n",
    "dsHad['pvalue']  = (('lat','lon'), p.values)\n",
    "\n",
    "y = dsHad.trend.sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "print('NS index 1', dsHad.trend.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat']))\n",
    "print('NS index 2', dsHad.trend.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat']))\n",
    "ymean = y.mean(['lat','lon'])\n",
    "ynorm = np.sqrt(((y - ymean)**2).mean(['lat','lon']))\n",
    "vol = np.sqrt((0*y + 1).mean(['lat','lon']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "bf33853e-018d-4e5a-8a09-c7d86e60e70a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HadiSST metrics= 0.04949121177196503 0.22246620059013367 0.22246620059013367 1.0\n",
      "CE: 0.27, EW: 0.50, PCC: 1.00\n",
      "['HadISST', 0.5167531371116638, 0.0, 1.0, 0.2747481167316437, 0.5009777545928955]\n"
     ]
    }
   ],
   "source": [
    "# Test:\n",
    "name = 'HadISST'\n",
    "ds = dsHad.copy(deep=False).sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "x = ds.trend + 0*dsHad.trend\n",
    "xmean = x.mean(['lat','lon'])\n",
    "xnorm = np.sqrt(((x - xmean)**2).mean(['lat','lon']))\n",
    "xycorr = ((x-xmean) * (y-ymean)).mean(['lat','lon'])\n",
    "PCC = (xycorr/xnorm/ynorm) \n",
    "print('HadiSST metrics=',xycorr.item(),xnorm.item(),ynorm.item(),PCC.item())\n",
    "RMS = np.sqrt(((x - y)**2).mean(['lat','lon']))/vol\n",
    "\n",
    "dsS = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat'])\n",
    "dsN = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat'])\n",
    "dsEQ = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_Eq, yNS_Eq)).mean(['lon','lat'])\n",
    "CE_index = (0.5*(dsS + dsN) - dsEQ).trend.item()\n",
    "\n",
    "dsW = ds.sel(lon= slice(xEW_W1,xEW_E1),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "dsE = ds.sel(lon= slice(xEW_W2,xEW_E2),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "EW_index = -(dsE-dsW).trend.item()\n",
    "print(f'CE: {CE_index:.2f}, EW: {EW_index:.2f}, PCC: {PCC:.2f}')\n",
    "\n",
    "data_hadi = [name,float(xmean.values),float(RMS.values),float(PCC.values),CE_index,EW_index]\n",
    "print(data_hadi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3dc8080f-98c2-430d-a91b-f28a21dbd327",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CE: -0.32, EW: -0.88, PCC: -0.62\n",
      "CE: -0.26, EW: -0.77, PCC: -0.67\n",
      "CE: -0.07, EW: -0.95, PCC: -0.47\n",
      "CE: -0.22, EW: -0.06, PCC: -0.03\n",
      "CE: 0.04, EW: 0.28, PCC: 0.37\n",
      "CE: -0.25, EW: -0.42, PCC: -0.45\n",
      "CE: -0.17, EW: 0.02, PCC: 0.03\n",
      "CE: -0.21, EW: 0.03, PCC: 0.05\n",
      "CE: -0.40, EW: 0.10, PCC: 0.05\n",
      "CE: 0.04, EW: 0.35, PCC: 0.41\n",
      "CE: -0.38, EW: -0.67, PCC: -0.66\n",
      "CE: -0.05, EW: -0.06, PCC: 0.04\n",
      "CE: 0.03, EW: 0.43, PCC: 0.38\n",
      "CE: -0.41, EW: -0.94, PCC: -0.54\n",
      "CE: -0.09, EW: 0.02, PCC: 0.15\n",
      "ESM | CE: -0.09, EW: 0.02, PCC: 0.15\n",
      "-0.1804612547159195 -0.2340356667836507 -0.130910899117589\n"
     ]
    }
   ],
   "source": [
    "lst_cntl_CE = []\n",
    "lst_cntl_EW = []\n",
    "lst_cntl_PCC = []\n",
    "name = 'CESM'\n",
    "    \n",
    "for i in range(15):\n",
    "    ds = ds_trend_cntl.isel(ensem_meber=i)\n",
    "    ds = ds.sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "    x = ds.trend + 0*dsHad.trend\n",
    "    xmean = x.mean(['lat','lon'])\n",
    "    xnorm = np.sqrt(((x - xmean)**2).mean(['lat','lon']))\n",
    "    xycorr = ((x-xmean) * (y-ymean)).mean(['lat','lon'])\n",
    "    PCC = (xycorr/xnorm/ynorm) \n",
    "    #print(xycorr,xnorm,ynorm,PCC)\n",
    "    RMS = np.sqrt(((x - y)**2).mean(['lat','lon']))/vol\n",
    "    dsS = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat'])\n",
    "    dsN = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat'])\n",
    "    dsEQ = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_Eq, yNS_Eq)).mean(['lon','lat'])\n",
    "    CE_index = (0.5*(dsS + dsN) - dsEQ).trend.item()\n",
    "    \n",
    "    dsW = ds.sel(lon= slice(xEW_W1,xEW_E1),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "    dsE = ds.sel(lon= slice(xEW_W2,xEW_E2),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "    EW_index = -(dsE-dsW).trend.item()\n",
    "    #print(name,dsS.trend.values)\n",
    "\n",
    "    print(f'CE: {CE_index:.2f}, EW: {EW_index:.2f}, PCC: {PCC:.2f}')\n",
    "    lst_cntl_CE.append(CE_index)\n",
    "    lst_cntl_EW.append(EW_index)\n",
    "    lst_cntl_PCC.append(PCC.item())\n",
    "\n",
    "ds = ds_trend_cntl.mean('ensem_meber').sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "x = ds.trend + 0*dsHad.trend\n",
    "xmean = x.mean(['lat','lon'])\n",
    "xnorm = np.sqrt(((x - xmean)**2).mean(['lat','lon']))\n",
    "xycorr = ((x-xmean) * (y-ymean)).mean(['lat','lon'])\n",
    "PCC_cntl = (xycorr/xnorm/ynorm) \n",
    "#print(xycorr,xnorm,ynorm,PCC)\n",
    "RMS = np.sqrt(((x - y)**2).mean(['lat','lon']))/vol\n",
    "dsS = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat'])\n",
    "dsN = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat'])\n",
    "dsEQ = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_Eq, yNS_Eq)).mean(['lon','lat'])\n",
    "CE_cntl = (0.5*(dsS + dsN) - dsEQ).trend.item()\n",
    "\n",
    "dsW = ds.sel(lon= slice(xEW_W1,xEW_E1),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "dsE = ds.sel(lon= slice(xEW_W2,xEW_E2),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "EW_cntl = -(dsE-dsW).trend.item()\n",
    "print(f'ESM | CE: {CE_index:.2f}, EW: {EW_index:.2f}, PCC: {PCC:.2f}')\n",
    "\n",
    "print(np.mean(np.array(lst_cntl_CE)), np.mean(np.array(lst_cntl_EW)), np.mean(np.array(lst_cntl_PCC)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "6fa5dd01-fa58-4f6b-be9e-0c001ed75a52",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CE: -0.11, EW: -0.26, PCC: -0.11\n",
      "CE: 0.04, EW: 0.41, PCC: 0.45\n",
      "CE: 0.02, EW: -0.37, PCC: -0.13\n",
      "CE: -0.12, EW: 0.08, PCC: 0.11\n",
      "CE: -0.11, EW: -0.26, PCC: -0.11\n",
      "CE: -0.02, EW: -0.27, PCC: -0.38\n",
      "CE: 0.27, EW: 0.57, PCC: 0.46\n",
      "CE: 0.27, EW: 1.32, PCC: 0.59\n",
      "CE: 0.42, EW: 0.34, PCC: 0.52\n",
      "CE: 0.11, EW: -0.06, PCC: 0.18\n",
      "CE: 0.21, EW: 1.08, PCC: 0.67\n",
      "CE: -0.15, EW: -0.07, PCC: 0.05\n",
      "CE: -0.26, EW: -0.59, PCC: -0.56\n",
      "CE: 0.04, EW: -0.25, PCC: -0.01\n",
      "CE: 0.07, EW: 0.86, PCC: 0.56\n",
      "ESM | CE: 0.07, EW: 0.86, PCC: 0.56\n",
      "0.044454526901245114 0.16889124711354572 0.15330835009614627\n"
     ]
    }
   ],
   "source": [
    "lst_fa_CE = []\n",
    "lst_fa_EW = []\n",
    "lst_fa_PCC = []\n",
    "name = 'CESM-FA'\n",
    "    \n",
    "for i in range(15):\n",
    "    ds = ds_trend_fa.isel(ensem_meber=i)\n",
    "    ds = ds.sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "    x = ds.trend + 0*dsHad.trend\n",
    "    xmean = x.mean(['lat','lon'])\n",
    "    xnorm = np.sqrt(((x - xmean)**2).mean(['lat','lon']))\n",
    "    xycorr = ((x-xmean) * (y-ymean)).mean(['lat','lon'])\n",
    "    PCC = (xycorr/xnorm/ynorm) \n",
    "    #print(xycorr,xnorm,ynorm,PCC)\n",
    "    RMS = np.sqrt(((x - y)**2).mean(['lat','lon']))/vol\n",
    "    dsS = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat'])\n",
    "    dsN = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat'])\n",
    "    dsEQ = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_Eq, yNS_Eq)).mean(['lon','lat'])\n",
    "    CE_index = (0.5*(dsS + dsN) - dsEQ).trend.item()\n",
    "    \n",
    "    dsW = ds.sel(lon= slice(xEW_W1,xEW_E1),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "    dsE = ds.sel(lon= slice(xEW_W2,xEW_E2),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "    EW_index = -(dsE-dsW).trend.item()\n",
    "    #print(name,dsS.trend.values)\n",
    "\n",
    "    print(f'CE: {CE_index:.2f}, EW: {EW_index:.2f}, PCC: {PCC:.2f}')\n",
    "    lst_fa_CE.append(CE_index)\n",
    "    lst_fa_EW.append(EW_index)\n",
    "    lst_fa_PCC.append(PCC.item())\n",
    "\n",
    "ds = ds_trend_fa.mean('ensem_meber').sel(lon= slice(xE,xW),lat=slice(yS,yN))\n",
    "x = ds.trend + 0*dsHad.trend\n",
    "xmean = x.mean(['lat','lon'])\n",
    "xnorm = np.sqrt(((x - xmean)**2).mean(['lat','lon']))\n",
    "xycorr = ((x-xmean) * (y-ymean)).mean(['lat','lon'])\n",
    "PCC_fa = (xycorr/xnorm/ynorm) \n",
    "#print(xycorr,xnorm,ynorm,PCC)\n",
    "RMS = np.sqrt(((x - y)**2).mean(['lat','lon']))/vol\n",
    "dsS = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_N,-yNS_Eq)).mean(['lon','lat'])\n",
    "dsN = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice( yNS_Eq, yNS_N)).mean(['lon','lat'])\n",
    "dsEQ = ds.sel(lon=slice(xNS_W,xNS_E),lat=slice(-yNS_Eq, yNS_Eq)).mean(['lon','lat'])\n",
    "CE_fa = (0.5*(dsS + dsN) - dsEQ).trend.item()\n",
    "\n",
    "dsW = ds.sel(lon= slice(xEW_W1,xEW_E1),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "dsE = ds.sel(lon= slice(xEW_W2,xEW_E2),lat=slice(-yWE,yWE)).mean(['lon','lat'])\n",
    "EW_fa = -(dsE-dsW).trend.item()\n",
    "print(f'ESM | CE: {CE_index:.2f}, EW: {EW_index:.2f}, PCC: {PCC:.2f}')\n",
    "\n",
    "print(np.mean(np.array(lst_fa_CE)), np.mean(np.array(lst_fa_EW)), np.mean(np.array(lst_fa_PCC)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9605bfcc-7b57-4c10-a9b9-14e87bd4ca7d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "afccfbcc-415e-4f9e-b607-3ff53b6b5596",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1sAAAFRCAYAAACCOLdKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABxK0lEQVR4nO3deZyN9fvH8dc1Yxi7QgiFQspSGdprWpR8k/YolVTaVNppLy1+raRFSJSKon1PoUULidJCkbKVJUR2c/3+OGemMWY5Z8yZ+8w57+fjcR5z7u1zX/c5zDXXfX/uz23ujoiIiIiIiJSslKADEBERERERSUQqtkRERERERGJAxZaIiIiIiEgMqNgSERERERGJARVbIiIiIiIiMaBiS0REREREJAZUbImIiIiIiMSAii2RYjKzpmb2l5lVDzqWwphZBTP7w8wygo5FRETAzNLMbI6ZHR50LEUxs3Fmdk3QcYiUVSq2pFSYWR0zG2Rmc81so5ktMrN3zaxTrnXmm5nn8xqQa52TzOwLM1tlZmvN7GczG55reWZ4m9VmVilPDC1ytVkrPK+Nmb1oZgvMbL2ZzTaz680skv8b9wJPuPvqHf+E8mdmh5vZG+HPy82sRz7rjMznM/sye7m7bwQeAP4vVnGKiEDi/a43sx4FxNon1zr7mdlWM/s8io+qF7DI3T+JYpuomFk9M3sh/NltNbOR+axT0PGl51rtTuCWeD+xKBKvygUdgCQ+M2sEfA6sAfoBMwkV+kcDQ4Ddcq1+F/BknibWhts5GngZuB04H9gK7AWclM9uVwOnA6NyzbsA+CPP/toCy4BzwsvaA8OANELFVEHH1DC832sLWqeEVAFmAc+GXwWZQOgYsm3Ks/x54EEz28fdfyjZEEVEEvN3fdg6YI888/7J9f4i4AngXDNr4e4/FdEewBXA3RGstyMqAMuBAYSKu4Jsd3zuviHX++/NbB7QHXg8BnGKJDQVW1IangAMyHD3tbnm/2Rmz+dZd427/1lAO52Br9w9d2L8BXgzn3VHAj0JJ2AzSyOUZIcAt2Wv5O4j8mw3z8z2B06l8AR8JvC9u/+RPSN81ekxoAswCGgMfA30dPffCmmrQO7+DvBOuP2Rhay6sZDPDXf/O3zWtRtwS3FiEREpQiL+rg9vnn+sZlYROAs4HKhEqNC7rrDGLNSluxnwVq55jYDfgNOAS4BDgPnAVe7+YRHxFRT0fODKcPunFb5qwfkj7A1C+UPFlkiU1I1QYsrMdgY6Ao/lSb4AuPvKKJr7E9jLzNpEsO5ooL2ZZZ+tO4HQWdNJEWxbDSgqrsOAafnMr0DojG5P4CCgBqGkD4CZHRbuElPY66YIYszrUDNbaqF7AIaZ2S75rPM1cEQx2hYRKVQC/64vymnA7+7+HfAcoatbaUVscxjwq7uvymfZPcCjQBtgKjDGzKpkL4wgf7xbjGOoaGa/m9lCM3vLzPbLZ52vCX3OFYvRvkhS05UtibU9CZ3pjKRbBcA9ZnZHnnld3f0tYDChJDXDzBYCXxHqPjc6n+T+N6EzcT2BmwmdbXwG8MJ2Hj7T2QM4u4g4dwdm5DO/HHC5u88Ot/cg8IyZpbh7FqECbd8i2v67iOV5vQe8QuisaCNCXVM+NrO24fu1si0OLxcRKWmJ+rseoLKZbbNfd88ugC4kVGQBTCbUJe9EYHwh7e0OLClg2SPu/mY4xpuAcwnljM/Cy/ctItb1RSzPazahz24mUBW4CvjczNq4+y+51ltMqMvlrsDcKPchktRUbEmsWZTrPww8nWfeEgB3/xf4X/gM5pHAgcB9QD8za+/uf+XZ7mlghJkNAToQ6pqxZ4GBmjUH3gYGunthiRKgIrAhn/kbswutsOwEVQP4293XA78W0XZU3H1Mrsnvzewb4Hfgf4SKsGzrw3GLiJS0Mv+7Pk9BNdrdLwm/X0c+RY6Z7Umou1+3cNwe7i55IYUXWwXlD4Dvcr1fHP6Z01PB3Us6f3wBfJE9bWZTCJ1IvIJwF8Sw7CJOOUQkSiq2JNZ+IXSGsQXwagTrrygqmbj7XEJn1oab2T3AHOBS4I48q04gdGP1s8DH7r4wnBy3Y2Z7AROBMe7eN4I4lwM75TN/S95wwz9Twvs5DCiqm8e9ee5ViIq7Lw6fDW6aZ9HOhG4QFxEpaYnwu37fXO9zD4DhBcR6IZAK/GGWU2taeD8N3X1BAYe2HMivqx7A5tw7Dbebc8tH3its+fjU3Y8vYp0CuftWM5tG/vkDlENEoqZiS2IqPDDD+0BvM3s0bxcQM6tRQL/1SM0ndNaxSt4F7p4VHlTiNkKjVeXLzPYGPgZecverI9zvt8De0QZLbLoRbsNCQx3XZ/tuKi2B6TvStohIfhLhd300V43MrBxwHqF7dN/Ks/g5QqMo3lXA5t8S+pyyu5dHY98ilkfbjXAbFqruWhPqVphbS2BxPlcVRaQIKrakNFwGTAGmmdmthLpJGKHuIf3YdnjeqmZWN8/26919dbh/fyVCo/P9Tqhr3pWEku8bBez7bkL9//MtYMxsH0LJdyJwb+59FzE60/uE7sUq5+55r2YVKNpuhOEbo7PP0KYAu5nZvoS6JP4RXn4HoS4rSwjdk3UfsJTtzy4fBtwa6b5FRKKUiL/rC/I/oBYwzN1X5NnXGOBSM7u7gGJqIpBOqKiZEc1Oo+1GGM4XEBoMJCs8vcndfwwvvx34ktCVyWqEPufWhK4g5nYYofuDRSRKGo1QYi487Pn+wIeEHqz7HaGkdyJwcZ7VbyNUNOR+ZQ81O5nQcOqjCN2E/T6h4uLEgh4M6e6b3X15IWcPTyfUH/7MfPZbmHcInUE8roj1dlQGobOg3xLqK39n+H32GdOtQCvgdUJdbEYRuuH5IHdfk92ImR0EVAfGxTheEUlSCfq7viAXABPzFlphLxMaBOOYAmJdQeh+2kgG59hR2fnjMEJD6n9L+HEiYTWAoYQ+5w8I9Yo43N2/zl7BQg84PpnQc8lEJErmXuiAPSJSADO7BDjd3Y8OOpaimNnLwLc7ci+YiIiUjPCVtonAnu7+T1HrB8nMLge6uPuxQcciUhbF7MqWmXU0s9lm9quZ5TvggJllmtkMM/vBzCbHKhaRGBkGfGRm1YMOpDBmVoFQ//tHgo5FkpvygkiIu/9A6OHHjYOOJQKbCY1OKBJzZjYi/NzQWQUsNzN7NJxHvgs/xiF7WZE5JggxubJlZqmEujR1ABYSejBft+w+wuF1ahDq290xfO/JLu6+tMSDERGRwCkviIhIUczscEIPJn/W3Vvms7wToeK/E3AAMMjdD4gkxwQlVle22hN6Ovo8d98EjAG65FnnLOAVd/8DQAlVRCShKS+IiEihwvdlFjYqcxdChZi7+5dADTOrR2Q5JhCxKrbqA7mfL7EwPC+3ZsBOZjbJzL4xs3NjFIuIiARPeUFERHZUQbkkkhwTiFgN/Z7fk+Tz9lcsB7QFjiY0ytoXZvalu8/ZrjGzXkAvgMqVK7fda6+9SjhcEZH49s033yx399pBx7EDSiwvKCeIiASbF1KqNXC2bIh6O1+/4gcg94ZD3X1oFE0UlEsiyTGBiFWxtRBomGu6AbA4n3WWu/u/wL9m9gnQhlB/y22Ev4ShABkZGT5t2rSYBC0iEq/M7PegY9hBJZYXlBNERALOC1s3ktbi5Kg32zR9+AZ3z9iBPReUS8oXMD9wsepGOBVoamaNzaw80JXtH0T4OnCYmZUzs0qEbnL7KUbxiIhIsJQXREQSiKWkRv0qAW8A54ZHJTwQWO3uS4gsxwQiJle23H2LmfUm9CDCVGCEu/8Qfi4R7j7E3X8ys/cIPfQwCxju7vkO8ygiImWb8oKISCKxkiqetm3V7EUgE6hlZguB24E0COUJQg/l7gT8CqwDzg8vyzfHlHiAxVDmHmqsLiMikozM7Jsd7HqRkJQTRCRZBZkXUqrU8fRWZ0a93fovByddLovVPVulavPmzSxcuJANG6K/UU+kMOnp6TRo0IC0tLSgQxGRCCknSLxSTpFEYYCllvyVrUSUEMXWwoULqVq1Ko0aNcIsv8FIRKLn7qxYsYKFCxfSuHHjoMMRkQgpJ0g8Uk6RhGJGSgy6ESaiWA2QUao2bNhAzZo1o06qc+ZsN/ChSA4zo2bNmjo7LlLGKCdIPFJOkUQT0AAZZU5CFFtA1El17ty5NG/enLlz58YoIkkEOisuUjYpJ0g8Uk6RhGGmYitCCVNsRWvcuHGkpKQwfvz4Emnvzz//pGvXruyxxx7svffedOrUiTlz5jB//nwqVqzIvvvum/N69tlnARgxYgStWrWidevWtGzZktdffx2AHj16UKlSJdasWZPT/lVXXYWZsXz5chYsWMCRRx5JixYt2GeffRg0aFC+Md1xxx3Ur19/m32vWrWqRI43GlWqVMl3fo8ePRg3blzE7dxxxx2YGb/++mvOvEceeQQzQzfIi8iOSIacAJCamrrNvufPnw+Efpemp6ezevXqfLcrLO7SVFDemDRpEieccELE7cyfPx8z49Zbb82Zt3z5ctLS0ujdu3eJxCoiAglyz1ZxvDLuZa7p2Y1Xxr3MDTfcsENtuTsnn3wy5513HmPGjAFgxowZ/PXXXzRs2JA99tiDGTNmbLPNwoULueeee5g+fTrVq1dn7dq1LFu2LGf5nnvuyeuvv0737t3Jyspi4sSJ1K9fH4By5crx0EMPsf/++7NmzRratm1Lhw4d2HvvvbeL7eqrr+a6667boeOLJ61atWLMmDHccsstQOgPpPyOW0QkGi+NG0+X8y7hpXHjEzonVKxYcbt9A7z44ou0a9eOV199lR49euR7XPnFXZY1adKEt956i/79+wPw8ssvs88++wQclUjZYIClJO01m6gk5af0xx9/MHfePG6/8kJ++fVXFixYsEPtTZw4kbS0NC655JKcefvuuy+HHXZYgdssXbqUqlWr5lz1qVKlyjY3zHbr1o2xY8cCoTN2hxxyCOXKhWrjevXqsf/++wNQtWpVWrRowaJFiyKOd+TIkZxyyil07NiRpk2b5vxhsXXrVnr06EHLli1p1aoVjzzyCBDqXtOxY0fatm3LYYcdxs8//wyEzjBeeumlHHnkkTRp0oTJkyfTs2dPWrRosV2yvvbaa9l///05+uijt/kDIts333zDEUccQdu2bTnuuONYsmRJvrGfdNJJOWd7582bR/Xq1aldu3bO8g8++ICDDjqI/fffn9NPP521a9cCcNddd9GuXTtatmxJr169yH7kQWZmJjfeeCPt27enWbNmfPrppxF/jiKSGP744w/mzZ1Lt8uv5ddff0m6nDB37lzWrl3L3XffzYsvvhj5gYZVqVKFm2++mTZt2nDggQfy119/AaHipWXLlrRp04bDDz8cCOWZ66+/nnbt2tG6dWueeuqpnGM64ogjOOOMM2jWrBl9+/bl+eefp3379rRq1Wqb7p0TJkzgsMMOo1mzZrz11lvbxfPvv//Ss2dP2rVrx3777ZeTM/KqWLEiLVq0yOkZMXbsWM4444yc5cuWLePUU0+lXbt2tGvXjs8//xyAr7/+moMPPpj99tuPgw8+mNmzZwMF51aRxKRuhJFKimJr//32o0qVyjmvZs2a0fWEDlSqmE63zsfStGnTbZbvv99+UbU/a9Ys2rZtW+DyuXPnbtP14tNPP6VNmzbUqVOHxo0bc/755/Pmm29us03Tpk1ZtmwZK1eu5MUXX6Rr1675tj1//ny+/fZbDjjggHyXP/LIIzn7PfLII3Pmz5gxg7Fjx/L9998zduxYFixYwIwZM1i0aBGzZs3i+++/5/zzzwegV69eDB48mG+++YYHH3yQyy67LKedlStX8vHHH/PII4/QuXNnrr76an744Qe+//77nDOg//77L/vvvz/Tp0/niCOO4M4779wmxs2bN3PFFVcwbtw4vvnmG3r27MnNN9+c7/FUq1aNhg0bMmvWLF588UXOPPO/ZzwsX76cu+++mwkTJjB9+nQyMjJ4+OGHAejduzdTp05l1qxZrF+/fpsEvWXLFr7++msGDhy4XWwiknj2229/KlepQuXKoVfTZs04vNNJVKhYiSM6ncyeTZvmLKtcpQr77bd/VO3Hc05Yv359zn5PPvlkIHRVq1u3bhx22GHMnj2bpUuXRhw3hH7HH3jggcycOZPDDz+cYcOGAaGTXO+//z4zZ87kjTfeAODpp5+mevXqTJ06lalTpzJs2DB+++03AGbOnMmgQYP4/vvvee6555gzZw5ff/01F154IYMHD97mGCdPnszbb7/NJZdcst2AE/fccw9HHXUUU6dOZeLEiVx//fX8+++/+R5T165dGTNmDAsXLiQ1NZVdd901Z9lVV13F1VdfzdSpUxk/fjwXXnghAHvttReffPIJ3377LXfddRc33XRTzjb55VaRhKR7tiKWFN0Ihw4bRreuZ7J/iz0YeHMfKldMp1LFdAAG3nwV915zMf+u30CfewYy/ae5DBs+vET3X1DXi/fee4+pU6fy0UcfcfXVV/PNN99wxx135Cw/5ZRTGDNmDF999VXO2b/c1q5dy6mnnsrAgQOpVq1avvsuqBvh0UcfTfXq1QHYe++9+f3339lnn32YN28eV1xxBf/73/849thjWbt2LVOmTOH000/P2Xbjxo057zt37oyZ0apVK+rUqUOrVq0A2GeffZg/fz777rsvKSkpOUVR9+7dOeWUU7aJZfbs2cyaNYsOHToAoTOf9erVy/d44L/k+P777/PRRx/xzDPPAPDll1/y448/csghhwCwadMmDjroICB0pvn+++9n3bp1/P333+yzzz507tw553MGaNu2bc79CyKSuIYNG8oZZ3al4V4tueDG/lSoWIkKFSsCcH7f/px11U1sXL+Op//vVhb8PIvhw4eV6P6DzAn5dSMcM2YMr776KikpKZxyyim8/PLLXH755RHHXb58+Zz7pdq2bcuHH34IwCGHHEKPHj0444wzcn7PfvDBB3z33Xc5912tXr2aX375hfLly9OuXbuc3/177LEHxx57LBDqPj5x4sSc/Z1xxhmkpKTQtGlTmjRpktPbItsHH3zAG2+8wYMPPgiERqf8448/aNGixXaxd+zYkVtvvZU6depsc/IOQlfQfvzxx5zpf/75hzVr1rB69WrOO+88fvnlF8yMzZs356yTX25t2LDhdvsVSQTJWjxFKymKrYyMDKZ/O4Pel13GUedcwdhBd9Gy2R4ApKSk8NvCxZx51W0cePChTP92PFWrVo2q/X322SeqgR6ymRnt27enffv2dOjQgfPPP3+bxNq1a1f2339/zjvvPFLy9IvdvHkzp556KmefffZ2xUskKlSokPM+NTWVLVu2sNNOOzFz5kzef/99Hn/8cV566SUGDhxIjRo1Cuynn91OSkrKNm2mpKSwZcuWAo87N3dnn3324Ysvvogo9s6dO3P99deTkZGxzR8U7k6HDh226wazYcMGLrvsMqZNm0bDhg254447tjkTmh139ucgIoktIyODmTO+5dLLLue2nqdy7UND2b3pXkDod9dfi/7goWt7cfghB/PWjG/LbE5YsGBBzkmlSy65ZJtujdm+++47fvnll5yTXZs2baJJkyb5FlsFSUtLy/m9nvv36JAhQ/jqq694++232XfffZkxYwbuzuDBgznuuOO2aWPSpEnb5ZDc+SX37+a8OSS/nDJ+/HiaN29eZOzly5enbdu2PPTQQ/zwww/bXFHMysriiy++oGK4EM92xRVXcOSRR/Lqq68yf/58MjMzc5bll1tFEpKZHmocoaToRgihfuyjnnuOo487nsee23a0qcdHv0KH4//HyGefjTqpAhx11FFs3Lgxp+sEwNSpU5k8eXKB2yxevJjp06fnTM+YMYPdd999m3V222037rnnnm267UEokVxwwQW0aNGCa665Jup4C7J8+XKysrI49dRT6d+/P9OnT6datWo0btyYl19+OWffM2fOjKrdrKysnD88XnjhBQ499NBtljdv3pxly5blFFubN2/mhx9+KLC9ihUr8n//93/bdTU88MAD+fzzz3NGK1y3bh1z5szJKaxq1arF2rVri/VHkIgklqpVqzL6uWc5oeOxvPvCiG2WvffiM5zYqSPPPTuqTOeEhg0bMmPGDGbMmJFvoQWhLoR33HEH8+fPZ/78+SxevJhFixbx+++/R3XM+Zk7dy4HHHAAd911F7Vq1WLBggUcd9xxPPnkkzlXg+bMmVNgF7+CvPzyy2RlZTF37lzmzZu3XVF13HHHMXjw4Jx7c7/99ttC27v22mv5v//7P2rWrLnN/GOPPZbHHnssZzr7pOPq1atzBicZOXJkVLGLJIrQABnqRhiJpLiyldvPP/7AJacex7IVKxk5/m16nPo/OhzSjqGvfljsNs2MV199lT59+jBgwADS09Np1KgRAwcOBP7r556tZ8+edOnSheuuu47FixeTnp5O7dq1GTJkyHZtX3zxxdvN+/zzz3nuuedo1apVTrv33nsvnTp12m7dRx55hNGjR+dMv/baawUex6JFizj//PPJysoC4L777gPg+eef59JLL+Xuu+9m8+bNdO3alTZt2hT1seSoXLkyP/zwA23btqV69eo5N3lnK1++POPGjePKK69k9erVbNmyhT59+hQ6KlR+9yvUrl2bkSNH0q1bt5yujnfffTfNmjXjoosuolWrVjRq1Ih27dpFHLuIJLZZP/7EQV3OYvXfy/no1bEcffKZtD7ocL56c0yx24znnJDXmDFjePfdd7eZd/LJJzNmzBhuvPHGbebnF/eVV15ZYNvXX389v/zyC+7O0UcfTZs2bWjdujXz589n//33x92pXbt2oXkpP82bN+eII47gr7/+YsiQIaSnp2+z/NZbb6VPnz60bt0ad6dRo0b5DqSRbZ999sk33zz66KNcfvnltG7dmi1btnD44YczZMgQbrjhBs477zwefvhhjjrqqKhiF0kY4Xu2pGiWfeanrMjIyPC8z1T66aef8u2LndeKFSto0rgxzz5wK5fd8RCtWrdm1vff8fjt13Lu9f2Z99tv253ZEon035dILJnZN+6eEXQc8WZHc0Kjxk246r7BDO3fl33btGLmd7O46Jb7GNTvCub/Nk85QUqccoqUlCDzQtpOu3nNo66Neru/XumTdLksaboRArz++uts3ryZy+58mJHPPsd773/AiJHPctmdD7N58+ac0ZJERCTxZeeE4ff05fnnRvHB++/z3KhnGH5PX+UEEZHCmLoRRiqpiq3PPvuUIzMz+XbGzJybgY899li+nTGTzCMO57PP9IwlEZFk8elnn3FEZibfzZixTU74bsYMjjjiCD797LOAIxQRiU+m52xFLKnu2Ro2bDgpKSnbjVy0yy678Pa77+XcqyQiIolv+LBhBeaE9959RzlBRKQQyVo8RSthii133y5h5pVayBCVZlbocklOZe2eRhEJUU6QeKScIglDA2RELCG6Eaanp7NixQr9EpMS5e6sWLFiu5GuRCS+KSdIPFJOkcSiboSRSogrWw0aNGDhwoUsW7Ys6FAkwaSnp9OgQYOgwxCRKCgnSLxSTpGEYeihxhFKiGIrLS2Nxo0bBx2GJLDly5cDoQcji0h8U06Q0qC8IMkse4AMKVpCFFsisXbaaacBMGnSpGADERGRuKC8IElN92xFTMWWSASuvTb6B/eJiEjiUl4QkUio2BKJQOfOnYMOQURE4ojygiQ7XdmKjIotkQj8+eefANStWzfgSEREJB4oL0iyS0kp/PEaEqJiSyQCXbt2BdQ3X0REQpQXJJmZGaZiKyIqtkQi0Ldv36BDEBGROKK8IMmuqAfHF7PNjsAgIBUY7u4D8iy/Hjg7PFkOaAHUdve/zWw+sAbYCmxx94wSD7AYVGyJRKBjx45BhyAiInFEeUGSXUl3IzSzVOBxoAOwEJhqZm+4+4/Z67j7A8AD4fU7A1e7+9+5mjnS3ZeXaGA7SMWWSAQWLFgAQMOGDQOORERE4oHygiQ1IxbdCNsDv7r7PAAzGwN0AX4sYP1uwIslHURJU7ElEoFzzjkHUN98EREJUV6QZGbEpNiqDyzINb0QOCDf/ZtVAjoCvXPNduADM3PgKXcfWtIBFoeKLZEI3HLLLUGHICIicUR5QZKbkVK8e7Zqmdm0XNNDcxVF+TXoBbTTGfg8TxfCQ9x9sZntAnxoZj+7+yfFCbIkqdgSicAxxxwTdAgiIhJHlBckqRW/G+HyQgauWAjk7pfbAFhcwLpdydOF0N0Xh38uNbNXCXVLDLzYSgk6AJGyYN68ecybNy/oMEREJE7ES15wd7KysoIOQ5KQpVjUryJMBZqaWWMzK0+ooHpju/2aVQeOAF7PNa+ymVXNfg8cC8wqoUPdIbqyJRKBnj17AuqbLyIiIUHnhfXr13PdbXfz618rySKF2hVTuP/WG2jQoEEg8UhyMSv50QjdfYuZ9QbeJzT0+wh3/8HMLgkvHxJe9WTgA3f/N9fmdYBXw8PRlwNecPf3SjTAYlKxJRKBO++8M+gQREQkjgSdFy6/4RbWtziWBm13AWDLpg1ceO0tvPPiCFJS1HFJYs9i8M/M3d8B3skzb0ie6ZHAyDzz5gFtSj6iHadiSyQCRxxxRNAhiIhIHAkyL6xatYpF640GO++SM69c+XTK7dmejydO4pijjwosNkkesXiocSJSsSUSgdmzZwPQvHnzgCMREZF4EGReWLVqFSkVq203P61aTRb/tbTU45HkY2Yl3o0wUcXsOrOZdTSz2Wb2q5n1LWS9dma21cxOi1UsIjvq4osv5uKLLw46DJEyTXlBEkmQeWH33XfH/l6A+7ajYq+d/SXHdzg6kJgk+cRggIyEFJMrW2aWCjwOdCA0jONUM3vD3X/MZ73/I3QjnEjcuvfee4MOQaRMU16QRBNkXjAz+vQ8i/tHDKPmgSdQLr0Sy6Z+QKe2e1G7du3A4pLkkqzFU7Ri1Y2wPfBr+GY1zGwM0AX4Mc96VwDjgXYxikOkRBx88MFBhyBS1ikvSEIJOi8ce/SRtN23NSNGv8jaZeu468YL2WOPPQKNSUS2F6tiqz6wINf0QuCA3CuYWX1CQzcehZKqxLlZs0KPamjZsmXAkYiUWcoLklDiIS/UrFmT66/qHdj+JYkZpGiAjIjEqtjK79P3PNMDgRvdfWtRo5mYWS+gF8Buu+1WEvGJRKV371Ay03O2RIqtxPKCcoLEA+UFSWaGuhFGKlbF1kKgYa7pBsDiPOtkAGPCCbUW0MnMtrj7a3kbc/ehwFCAjIyMvMlZJOYeeOCBoEMQKetKLC8oJ0g8UF6Q5Ja8A15EK1bF1lSgqZk1BhYBXYGzcq/g7o2z35vZSOCt/AotkXjQrp16NInsIOUFSSjKC5LUDA39HqGYFFvuvsXMehMaTSoVGOHuP5jZJeHlQwptQCTOzJgxA4B999030DhEyirlBUk0yguS7PRQ48jE7KHG7v4O8E6eefkmU3fvEas4REpCnz59APXNF9kRyguSSJQXJJmF7tkKOoqyIWbFlkgiGThwYNAhiIhIHFFekKSmboQRU7ElEgF1ExERkdyUFyTZaYCMyKjYEonA1KlTAd0QLSIiIcoLktxM92xFSMWWSASuv/56QH3zRUQkRHlBkpmpG2HEVGyJROCxxx4LOgQREYkjyguS7NSNMDIqtkQi0LJly6BDEBGROKK8IMnMDFJVbEVExZZIBKZMmQLAwQcfHHAkIiISD5QXJNmp2IqMii2RCNx0002A+uaLiEiI8oIkM8NUbEVIxZZIBJ566qmgQxARkTiivCBJTd0II6ZiSyQCzZs3DzoEERGJI8oLkswMFVuRSgk6AJGyYPLkyUyePDnoMEREJE4oL4hIJHRlSyQCt99+O6C++SIiEqK8IMnMDMrpylZEVGyJRGDEiBFBhyAiInFEeUGSmboRRk7FlkgEmjRpEnQIIiISR5QXJKmZRiOMlO7ZEonAhAkTmDBhQtBhiIhInFBekGQWurKVEvWryHbNOprZbDP71cz65rM808xWm9mM8Ou2SLcNiq5siUTg7rvvBuCYY44JOBIREYkHyguS7Er6ypaZpQKPAx2AhcBUM3vD3X/Ms+qn7n5CMbctdSq2RCLw3HPPBR2CiIjEEeUFSWYWm+dstQd+dfd5oX3YGKALEEnBtCPbxpSKLZEINGzYMOgQREQkjigvSDIzYnLPVn1gQa7phcAB+ax3kJnNBBYD17n7D1FsW+p0z5ZIBN577z3ee++9oMNIGHfffTeZmZlBhyEiUmzKCyVHOaFsSjWL+gXUMrNpuV69cjWZX/XmeaanA7u7extgMPBaFNsGQsWWSAQGDBjAgAEDgg6jREycOJFd6tRj4sSJEa2fmZmZc29CJPOLo1GjRowePTpn+uWXXyYjI4MaNWpQo0YNWrVqxeDBg3OWz5s3j9NPP526detSpUoVGjZsyMknn8ymTZu45JJLqFKlClWqVKFSpUqYWc50lSpVeP7550skZhFJbomSF5QTlBOKI7sbYbQvYLm7Z+R6Dc3V7EIg9yXjBoSuXuVw93/cfW34/TtAmpnVimTboKgboUgExowZE3QIJWLixImccOJJbKjWjBNOPIm33niNI488MuiwtjFlyhR69uzJyy+/TIcOHdi6dSvff/89v//+e846nTp14thjj2X27NlUq1aNRYsW8dZbb+HuDBkyhCFDhgDw2Wefcdhhh7F27dqgDkdEElQi5AXlBNkRMehGOBVoamaNgUVAV+Cs3CuYWV3gL3d3M2tP6MLRCmBVUdsGRVe2RCJQt25d6tatG3QYOyQ7qW6sczCpu7RiY52DOeHEkyI+m1mY888/n4YNG1K1alX23ntvXnjhhW2Wv/322+y9995UqVKFE044geXLlxfY1hdffEGLFi3o2LEjqamplC9fnrZt23LKKacAsGLFCmbPns0ll1xC9erVMTMaNGjAJZdcQoUKFXb4WEREIlHW80IkOcG9eL2wlBMSnxmUS7GoX4Vx9y1Ab+B94CfgJXf/wcwuMbNLwqudBswK37P1KNDVQ/LdNkaHHxUVWyIRePPNN3nzzTeDDqPYcifVlKr1AEipWq/ECq5DDz2UGTNmsGrVKm677TZ69OjBjz+GBgCaN28ep5xyCjfddBOrVq3iyiuvZNiwYQW2dcghhzB9+nSuuuoq3n33XZYuXbrN8po1a7LPPvtw4YUX8uyzz/Ljjz8W+w8CEZHiKst5oaicMPixJzjzwsvpclEfTrmgN4889mRUv2eTKSd8OuULOp97CYefdRnHnX0xw0YlR7fE7AEyitGNsFDu/o67N3P3Pdz9nvC8Ie4+JPz+MXffx93buPuB7j6lsG3jgYotkQg89NBDPPTQQ0GHUSz5JdVskRZc99xzT05f+ezXZ599lrP8ggsuoGbNmqSmptK1a1dat27NpEmTAHjxxRdp37493bt3p1y5chx77LGcdNJJBe7rwAMPZPLkySxfvpxevXpRt25dMjIy+PTTT3PWmTRpEpmZmQwcOJB9992XOnXq0L9/fxVdIlJqympeiCQnXH3NtZTbI4PWp17C3qddzvfrKvHEsBE56yknhPz2229cP3AUGw/oTsVDz4ZDzmXU1AU8/9K4mO43XsSi2EpEKrZEIjBu3DjGjSubvzzP7HoWG6o12y6pZkupWo8N1ZpxZteCuzbffPPNrFq1apvXoYceCkBWVha33XYbzZs3p3r16tSoUYOZM2eybNkyABYuXEijRo22aa9x48aFxnzIIYfw/PPPs2DBAn7//Xf23HNPTjjhBFatWgVArVq1uPfee5k+fTqrVq3i/vvv56677uKZZ56J8FMREdkxZTUvRJITvHZLXrv/hpx5u7Zqz+SpM3KmlRNCHnnqGSofcApm/xURlfc+jLFvfxzT/caDHRggI+mo2BKJQK1atahVq1bQYRTL2DEvkP7PHLLWLMl3edaaJaT/M4exY17Id3lRXnzxRYYPH8748eNZuXIlq1atok2bNjlnFOvXr8/8+fO32ea3336LuP2GDRty8803888//zBv3rztlleqVIkePXrQunVrZsyYUaxjEBGJVlnNC5HkhNSVP3PyjQ9sM3+rpUbUfjLlhJVr/qVcxcrbzd/oif/ntaFiK1KJ/69BpAS88sorvPLKK0GHUSxHHnkkb73xGhX+mrJdcs1as4QKf03ZoRGo/vnnH8qVK0ft2rXJyspixIgRzJw5M2d5t27d+Oqrr3jxxRfZsmULEyZM4PXXXy+wvddee41nnnmGJUtCsS5fvpyBAwdSq1Yt9tprL1auXEm/fv2YNWsWmzdvZsuWLYwfP55Zs2Zx2GGHFesYRESiVVbzQlE5odziTzmy+6U0avPf82C3bNpIjQqRFVvJlBMOaduKtQvnbDNv66YN1KuWHtP9xgVd2YqYii2RCDz66KM8+uijQYdRbPkl15IotADOO+88DjjgAPbcc0/q16/Pjz/+uE2C22OPPRg3bhx33XUXNWrU4JFHHuHCCy8ssL2aNWvy8ssvs99++1G5cmVatmzJsmXL+PDDD6lUqRLly5dn6dKlnHLKKey8887Url2bu+++m8GDB3P66acX+zhERKJRlvNCYTnh7Tdfp0H5zcyfOpEtmzaybP4cvh8ziNuvvyqitpMpJ1xwztnUXjCFNXO/xd1Z9+d8Nk1+hnv6XR3T/UrZYmXthvKMjAyfNm1a0GFIklm9ejUA1atXDziSHZP7mSrp/8yJy2eqSP7M7Bt3zwg6jnijnCBBSYS8UFBOcHcmTprMex9PpmmT3ene9UwqVqwYdLhxaevWrbz65ttM+GQKezfbgwvOOYvKlbfvWhgLQeaFBs1bee8h0V/Z7XdUs6TLZXqosUgEynIyzS37bOaZXc9irAotEZFiS4S8UFBOMDOOOjKTo47MDDK8MiE1NZXTTjqR0046MehQSl2ydguMlootkQiMHTsWgDPPPDPgSHbckUceydK/8r8xWkREIpMoeUE5QYojezRCKZqKLZEIPPnkk0DZT6oiIlIylBckmWWPRihFU7ElEoF33nkn6BBERCSOKC9IUtOVrYip2BKJQKVKlYIOQURE4ojygiQzw0g1FVuRULElEoHRo0cD0L1794AjERGReKC8UHa4O8uXL6dSpUqlNlJgMkhRsRURFVsiERg+fDigpCoiIiHKC2XDlK++5tYHnmR1ag1St26gRZ1KPPng3RrKfgcZkKpaKyIqtkQi8OGHHwYdgoiIxBHlhfi3Zs0a+vQfjB1wLuXCV2G+W72MK/rewfBB/xdwdGWcQYru2YpISqwaNrOOZjbbzH41s775LD/bzL4Lv6aYWZtYxSKyo9LS0khLSws6DJEyTXlBEonyQvwb9cJLbGySieXq7la+em1mLfibLVu2BBhZ2Re6smVRv5JRTK5smVkq8DjQAVgITDWzN9z9x1yr/QYc4e4rzex4YChwQCziEdlRI0eOBKBHjx6BxiFSVikvSKJRXoh/K1atJrVi/e3mbyWVrVu3Uq6cOnjtCN2zFZlYXdlqD/zq7vPcfRMwBuiSewV3n+LuK8OTXwINYhSLyA4bOXJkTmIVkWJRXpCEorwQ/84942SYN2WbeVlbt1CnYhYVKlQIKKrEkH3PVrSvZBSrkr4+sCDX9EIKPzt5AfBujGIR2WGTJk0KOgSRsk55QRKK8kL8a9y4MWcc3IyXP3+NLQ0y8PWrqLpkKg8/dGfQoZV9ZrpnK0KxKrby+/Q93xXNjiSUVA8tsDGzXkAvgN12260k4hMRkdJVYnlBOUFEItX36ss598zFvPLmu9TdpTFdTrhc99qVAEPdCCMVq2JrIdAw13QDYHHelcysNTAcON7dVxTUmLsPJdR3n4yMjHyTs0gsDRs2DICLLroo4EhEyqwSywvKCRIPlBfKjl133ZXeF18QdBgJJ1m7BUYrVvdsTQWamlljMysPdAXeyL2Cme0GvAKc4+5zYhSHSIkYO3YsY8eODToMkbJMeUESivKCJLPsK1vRvpJRTK5sufsWM+sNvA+kAiPc/QczuyS8fAhwG1ATeCI8JOcWd8+IRTwiO2rChAlBhyBSpikvSKJRXpCkZpCqe7YiEnGxZWZVgBrAKndfW9T67v4O8E6eeUNyvb8QuDDiSEVEJK4oL4iIiBSu0G6EZtbSzAab2TxgNfAHsNrM5prZY2bWqlSiFAnYE088wRNPPBF0GCKBU14QCVFekGSmboSRK7DYMrMXgReAJUB3oBZQPvzzHGAR8LyZjSmFOEUC9eabb/Lmm28GHYZIoJQXRP6jvCDJTs/Zikxh3QhfcPf8fousBKaEX/eZ2QkxiUwkjrz7rh73I4LygkgO5QVJZkZsrlSZWUdgEKF7e4e7+4A8y88GbgxPrgUudfeZ4WXzgTXAVuLont8Ci60CEmp+671VcuGIiEi8Ul4QEREgJgNkmFkq8DjQgdDjQqaa2Rvu/mOu1X4DjnD3lWZ2PKHHgByQa/mR7r68RAPbQUXds9XWzAYVsOxeM2sfm7BE4sugQYMYNCjf/woiSUV5QSREeUGSWeierehfRWgP/Oru89x9EzAG6JJ7BXef4u4rw5NfEnpmY1wr6jlbNxDqFpKfz4FbSzYckfj00Ucf8dFHHwUdhkg8UF4QQXlBJNUs6lcR6gMLck0vDM8ryAVA7v68DnxgZt+YWa9iHVQMFDX0ezsKHob3Q2B4yYYjEp/eeOONolcSSQ7KCyIoL0hyyx6NsBhqmdm0XNND3X1ormbz8nz3b3YkoWLr0FyzD3H3xWa2C/Chmf3s7p8UJ8iSVFSxVQNIK2BZeaBaiUYjIiLxrgbKCyIiyc0gtaj+cflbXsjAFQuBhrmmGwCLt9u1WWtCJ/aOd/cV2fPdfXH451Ize5VQt8TAi62iPqYfgOMKWHYc8FPJhiMSnx588EEefPDBoMMQiQfKCyIoL0hyi9FztqYCTc2ssZmVB7oC21xCNrPdgFeAc9x9Tq75lc2savZ74FhgVskdcfEVdWXrEeApM9sKjHf3rWaWApxMaLSQq2MdoEg8+OKLL4IOQSReKC+IoLwgyS6ie7Ci4u5bzKw38D6hod9HuPsPZnZJePkQ4DagJvCEhfafPcR7HeDV8LxyhB5V8l6JBlhMhRZb7v6Kme1K6FLdaDNbTujhlRuB29z9xVKIUSRw48ePDzoEkbigvCASorwgyWwH7tkqlLu/A7yTZ96QXO8vJJ/7ht19HtCmxAMqAUVd2cLdHzOzkcDBhBLqCuALd/8nxrGJiEgcUl4QEUlyxb9nK+kUWmyZWTvgHmAz0MfdPyiVqETizIABoQeY9+3bN+BIRIKlvCDJZMOGDXz55ZfsvPPOtGrVCst1Jl95QZJZrK5sJaKirmw9AJxEqB/kI8AJsQ5IJB7NmDEj6BBE4oXygiSF119/g/cnTma/Aw/jn5/n8egTQ7jnztupU6cOoLwgolorMkUVWylAFqEx7lNjH45IfBozZkzQIYjEC+UFSXjLly/no8++4PIbb8+Zd3iH47lnwP08+shDgPKCSEq+j8WSvIrqbXkDMA4YBFwT+3BERCTOKS9Iwnvl1dc4/pQzt5lXqVJl0tIrs379+oCiEokfRujKVrSvssTMuuczz8ysXzTtFDUa4ZeExqkXSWr9+/cH4NZbbw04EpFgKS9IMkhNTSFra9Z28x3PuW9LeUEk4d1uZp2BS9x9pZk1AZ4j1Lvjvkgb0TgiIhGYPXs2s2fPDjoMEREpBSefdBLvjN/2KQZr/llN1qYNpKenA8oLIikW/auM2Rf4B/jezPoDXwNvAUdE00iBV7bMbCpwP/C6u2/KZ3l5QjdJX+vuB0SzU5GyZvTo0UGHIBI45QVJFjvvvDOdOx7D4Htvp2Xb9qxe+TcL5s7mnjv/u4dLeUGSWhnsFhgtd//XzG4CDgBuBkYBA9zdo2mnsG6E5wF3AU+a2XRgNrAGqAo0A/YHPgZ6RB29iIiURcoLkjQ6HnccRx91FDNnzqRGjRrsueflQYckEjcMS/gBMszsf8Aw4GXgLOAp4FMzO8fdf4u0nQKLLXf/ETjNzOoCHYBWhB5euRJ4FjjH3ZcW/xBEyo7bbrsNgLvuuivgSESCo7wgySYtLY2MjIx8lykvSLJL9CtbwBDgPHf/EMDMDiN0hWsaUDPSRooa+h13/5PQzWAiSWvBggVBhyASN5QXRJQXYmHt2rU8/dwLzF+wiNNO6MjBBx24zYOkJb6UwXuwotXa3VdmT7h7FtDfzN6OppEiiy0RgWeeeSboEEREJI4oL5Ssn2fP4eKb7iG97f+ouPte3DxmMs1ffo0hDw9QwRWnEv1bCY9A2AHoBtR2985mlgFUi6YdjUYoIiIiEpA5c+YEHUJcuPX+QezcsRdV6u5OaloFau57JHOyavLJp58HHZrkw4AUs6hfZYmZXQE8CcwBDg/PXg/cHU07KrZEItCvXz/69YvqGXYiIpLASiIvzJ07l+bNmzN37twSiqrsWrERUlK37XBVo+UhvPTGOwFFJEVJ9IcaA32AY9x9AKFnawH8DDSPphEVWyIRWLFiBStWrAg6DBERiRMlkRfGjRtHSkoK48ePL6Goyq40tn+I9IZVy6lfb5cAopFIpBTjVcZUBbJvzswe7j0N2O7RJ4WJ6LjNbLiZVcozr56ZvRfNzkTKqqFDhzJ06NCgwxCJG8oLkuxKIi+8NG48Xc67hJfGlWyx5e5M+fIrbrz1Th598in++eefEm0/Fg7fby/W/P5jzrRnZfHv169zac/zAoxKChK6UmVRv8qYT4C+eeZdCUyMppFIB8ioCnwXHlf+CzPrCgwGhkezMxEp2ptvv8NbH3yMm9Fmr6ZcfGFPypXTWDYSd5QXJO5s3bqVF55/np9++gkHjj/+eA4//PAitwvCH3/8wby5c+k7ZAwXdchgwYIFNGzYcIfbdXeu7nsLS9Nq0ah9F35Z/iddL72GB26+ln32blECkcfGTddcSdaDg/j045FszDJ2Lu8MuqUPO+20U9ChSQGSYDTCK4A3zewioKqZzQb+ATpH00hEf8G5+5lmdjbwenhH9YCT3F13LUpSuO666wB48MEHY7qfQU88xbx/IeOs3pgZC3/5gSuu68uTA0P7/e233xj01HDWbthM9UrpXH1ZLxo0aBDTmETyo7wg8eiWW26h0/EdOf20U8nKymLM2LHMnz+fc889t8T3FW1e2G+//Znzy5yczkhbtm6h42lnU6FiJY7odDJ7Nm1Kuex7lgyaNW3Gt99OjzquKV9+xdLytWh66PEA7Fy/ETXOvpp7Bz3B808Njrq90mJm3Hp9HyBUMJbBqyBJJ9G/IndfYmbtgPbAboS6FH4dHgI+YtGcLl8EbACaAD8CuptTksb69etjvo+NGzfy9azZHNnjqpx5DZruw59zf2bWrFlUSE/n+nsf5sBul1E+vSIb1q2l9839eeK+29l1110BmPzpZ4wcM46slDRSszZxac9zade2bcxjl6SlvCBxY9asWezRpDHtwg8hTklJ4axu3bj9jjvZvHkzaWlpJbq/aPPCsGFDOePMrjTcqyUX3NifChUrUaFiRQDO79ufs666iY3r1/H0/93Kgp9nMXz4sGLF9epb79L4wFO3mZeSmsraLWXnjhkVWvHPKJP3YEXN3R34Kvwqlkjv2XoQGEOon2IjYAah7iOnF3fHImXJ448/zuOPPx7TfSxZsoSqdepvN79+8zZ8M+M7Bg0ZzkFn9aZ8eig5p1eqQrszL+GRJ0L3DHw9bRrDXnmPdt2v4sCzL6ft2Vfx4PDn+Xn27JjGLclJeUHizfRvvuHAAw7Ybv7uu+/OX3/9VeL7izYvZGRkMHPGtzTcuRq39TyVpYsXkJIS+jMsJSWFvxb9wa3nn8JuNaszc8a3tC3mibIGu9bln+V/bjc/la3Fak+kIIl4z5aZLTCzP4p6RdNmpEVpC6CNu7/m7pvd/XrgNOD+aA9CRPJXr1491ixdtN38xXNm0abVPqzbvJW0ChW2WVaxajVW/Rs6u/r06DEcdNr5Ob/MUlJSOPjMXjw54tnYBy/JSHlBSsV3333H7bfcxC03Xs/IZ55m8+bN+a7XYu+9mTFz5nbzFy5cSO3atWMdZkSqVq3K6Oee5YSOx/LuCyO2Wfbei89wYqeOPPfsKKpWrVrsfZzf/SzmTXiJrK3/FVdLf51F272aFLtNke1Y6J6taF9lQHfgnPDrUWA10B+4MPxzJTAomgYjKrbc/X/u/leeeZ8AraPZmUhZ1adPH/r06RPTfVSoUIH9mzfhxykfE7pqDUvm/cLWpfPZt00b0sulsGXTxm22Wb92DdUrpQOwlRRSUlO3WV4+PZ0Nm7fENG5JTsoLUhref/dd3n9tLDdc0I07rrqQlvVrcsM1fXJ+R+aWkZHBtG+mM3/+/Jx5EyZ8RP0GDaiQ50RVSdiRvDDrx59oc/ARrP57Oa88/Tir/15O64MOZ9aPP+1wXNWrV2dAv6uZ99oQZr78JLNefpz66/6g37XFi1WkIFaMV7xz98nZL6AH0NHdh7n7B+4+DPgfcH40bUbajbCCmd1jZvPMbHV43rGAxuMUKUHXXtmbg3fbia9fGMxXzz9GlWWzefzh0IWCPpdcyBcvPMHmjRsA2LjuX6aOeZKrLrkQgOqVKrB+7bbD+65esYx6NTWSk5Q85QUpDR++9zZXX3guFSqUB2Dfli04sn0bJk3cfuRlM+P//u//ePf997mz/93cfsedrNuwnssvv7y0wy7UihUr+GbaNFLLpXHtaceyeNbXXHf6cZRLK8+0qVNL5JmOrVruw/NDB/PqsIGMHz6YW66/JqfLoohEbFdgbZ55a4Ht7/kohOV3dmi7lcyeDO9wAPCuu9cws/rAB+6+TzQ73FEZGRk+bdq00tylSNz45ddfefSpp1m3eStV09O45rKL2W233QD4888/ufSGm8k4pSc7163PsoW/M+P1UYwY/DA1atQINnDZYWb2jbtnBB1HtnjJC8oJiWvz5s3ce1tfbr6i1zbz16xdy/BXPuD6vv0CimzHjBgxgssu702NnWrw3KhRdOjQgQ8++IBze/Rg1cpVPPnE45x/flQnziVJBZkX9tt/f5/4yWdRb7dT1cpxlcsKY2YjgcbA3cBCoCHQD/jD3SM+sRjpaIQnAXu6+79mlgXg7ovCiVVESknTPfdk8AP35busbt26PPPoQwx5+hl++HwpTRo24NknH6VKlSqlHKUkiZNQXpAYKleuHGvXbdxu/swfZ9Nsr/h9XlRRPv3sM47IzOS5USPZZZddADj22GP5bsYMzjn3PD797DMVW1ImlIUBL3bQJcAdwBBCJxeXAC8Bd0bTSKTF1qa865pZbaDAa91m1pHQDWSpwHB3H5BnuYWXdwLWAT3cPfoHSoiUguxuKLEekXBH1ahRg77XXh10GJIclBckpsyM1vtn8Pr7H9PluKMAWPH3Ssa+NYHBQ4YGHF3x88LwYcNISUnZ7g/VXXbZhffefYesrKge4SMSCCM2A17sSJ4oattoufsGoG/4VWyRFlsvA6PM7GoAM6sHDCQ07O92zCwVeBzoQOiy21Qze8Pdf8y12vFA0/DrAODJ8E+RuFMx/CwUEcmhvCAxd855PXj9tVe55aEnSU1JoWKVatz34MOUKxfNY0Jjo7h5ITXPQEa5mVmhy0XiSUnXWjuSJyLctjgxNQfaANt0E3L3Eflvsb1If1vdRGg43++BSsAvwDAKvozWHvjV3eeFAx0DdCH00MtsXYBnww8L+9LMaphZPXdfEmnwIqXlwQcfDDoEkXijvCClostJJ9PlpJODDmM7yguS3IyUku9GWOw8Qeh5j0VtGxUzuwm4DZhJ6CpaNgciLrYiHfp9k7v3cfcqQB2gqrtf7e6bCtikPrAg1/RCth+5I5J1ii0zM5ORI0cCoZtsMzMzGT16NADr1q0jMzOTsWPHArB69WoyMzN55ZVXAFi+fDmZmZm8+eabQGjggczMTN577z0AFixYQGZmJhMmTABg3rx5ZGZmMnnyZABmz55NZmYmU6ZMAUJPtc/MzGTq1KkAzJgxg8zMTGbMmAHA1KlTyczMZNasWQBMmTKFzMxMZocfRjt58mQyMzOZN28eABMmTCAzM5MFC0If33vvvUdmZiZ//hl6iOGbb75JZmYmy5cvB+CVV14hMzOT1atXAzB27FgyMzNZty7072b06NFkZmbmPLtk5MiRZGZm5nyWw4YN45hjjsmZfuKJJzj++ONzpgcNGsSJJ56YM/3ggw9y6qn/Pb1+wIABdO3aNWe6f//+dO/ePWf6tttu26Z/er9+/ejV678boq+77rptRpPKO9zu5ZdfznXXXZcz3atXL/r1++/G6fPPP5/bbrstZ7p79+70798/Z7pr164MGPDfleZTTz11myR64oknMmjQf49UOP7443niiSdypo855hiGDRuWM61/e/q3ly2Sf3tlVVnLC/p/qf+X2ZQT9G8vW7z/24trBlaMVxF2JE/EIn/0Adq7+wHufmSu11HRNFLglS0zK+zpd1Wz+xpnV5B5N89nXt5hDyNZJzuWXkAvIGfkNZHSNHbsWH744QeGDg3+PgGRoMRLXlBOkHiQ+49wkWRj7lgEI5rno5aZ5R5Cdqi7Z/9xtSN5IuK6IgrrgZ93sI2Ch34Pjy6VHXz2StkHkrORu2/XudjMDgLucPfjwtP9wuvel2udp4BJ7v5ieHo2kFlUdxEN8ytByD4zdd99+Y8EKBJr8TD0ezzmBeUEyS37b5rSGCVNeUGCFmReaLv//v75Z59GvV3FylUKjHlH8gShboSFbhstMzsXOITQiIR/5V7m7hGPZFPglS13z+liaGbnA8eEd/Y7sDuhPowfFbD5VKCpmTUGFgFdgbPyrPMG0Dvcp/IAYLX65Uu8UjIVUV6Q+LV69Wruv60fW1aGuq+l1qjDjf0HUL169ZjtU3lBkp1FXm9Eqth5wsyWRbBttEaGf16Ya172ycaIR7KJdICM/kBTd18fnv7FzC4G5uQKJIe7bzGz3sD74WBGuPsPZnZJePkQ4B1Cwzb+Suims/PztiMiInFLeUHixu1XX851hzRm5yoNAfh77Tpuv/pyBo4YHXBkIonKoYSLrR3JEwVtu4MhNd7B7YHIi60UQpfnfso1b3cKqerc/R1CH0jueUNyvXfg8rzbicSj7Btmn3nmmYAjEYkbygsSF3755ReaVXZ2rlIpZ97OVSrRvHJoWdOmTWOyX+UFSXrFu2eriCaLnyfy23YHY/m9JNqJtNh6BPjYzJ4hNNJHQ6BHeL5IwmvYsGHQIYjEG+UFiQvLli2jXtUK282vV60CS5cujVmxpbwgSc1L/spWPDCzoe7eK/z+OQoYZMPdz420zYiKLXd/wMy+B04H9gOWAD3d/b1IdyRSFqxevZoVK1aw++67b/NgybvuuivAqETij/KCxIv99tuP259cxXFttp3/xYJV3Nm2bcz2q7wgyS4G92zFg99yvf+1JBqM+BHs4QSqJCoJafPmzfS/+24sJZU6devyy+zZHHXkkXTufELQoYnELeUFiQcVK1akfadTeej9Vzn3oOYAPPvlbDI6nkx6enrA0YkksAQstnKPXujud5ZEmxEVW2ZWnlD3kH2BKnmCivgymki8evjhh+ly8qnsseeeOfMGPvQgLVvuQ+PGjXMeepj9IEqRZKe8IPHktG5ns+DQw3n5udD9Uz3ueCTm3fyUFyS5JWY3wliI9MrWKKAN8CZ5xpkXSQTLV6zYptACOP/Ci3jh2ZH069eP5s2bBxSZSNxSXpC40rBhQ6656bZS25/ygiQ1R8VWhCIttjoCjd19VQxjEQmMpaRsN69ixYps2LgRgFtvvbW0QxKJd8oLktSUF0QkEpEWW38A2w/1I5IgyqWmsnrVKqrXqJEzb/y4l+l8gu7ZEimA8oKISNJyyNKVrUhEWmw9C7xuZoPI013E3T8u8ahEStl1117LTTffwmGHH8HujRsx+eOPqVy5MhkZGQB07doVgDFjxgQZpkg8UV6QpKa8IMkuQUcjxMw6Aye4+8X5LHsKeM3d3420vUiLrd7hn/fmme9Ak0h3JsX37NAnmTF5AmlZWXjVnbji1jv1jI8StPPOO/PkE4/zxRdfsOCP3zm/x3nsuuuuOcv33Xff4IITiU/KC5LUlBck6SVosQVcAxR0A+hzwF1AyRZb7t440gal5A0fPIjqP0+h734NANiweQt3XX4hg8a8SqVKlQKOLnGYGQcffHC+y/r27VvK0YjEN+UFSXbKC5LU3EOvxLS3u39awLLPgX2iaWz7UQEk7vzw+ccc0bhuznR6Wjm6N6vNy6OfDTAqEREREUlanhX9q2yoaGZVC1hWBagYTWMFXtkys5/cvUX4/QJCXUO24+67RbNDiV7a1s3bzdujVg0m/TongGiS06mnngrASy+9xIIFC9h5552pVq1awFGJlC7lBZH/ZOeF8ePHBxyJSDAS9Z4t4FvgNOCZfJadAsyIprHCuhFelOt992galZK1Ma0S7o6Z5cybNO9PjjjzkgCjSi4HHXQQP/30E1dcfR2779GUFcuWUj4Fbr/1ZtLS0oIOT6S0KC+IhB100EFBhyASoIR+qPG9wEtmthMwHlgC1ANOBW4FzoymsQKLLXf/LNf7ycUKVUrE2b2v5v4H7+KiNruzU6V0Pp//F99srczZRx0ddGhJ4/TTT+exoSPoddV1OfPmz/uVBx9+hH433hBgZCKlR3lB5D/XXXdd0SuJJLIELbbc/X0zuwB4CHiQ/3pxLAAudPcPomlP92yVAe0PPoRrBo/gtS21GPT7RsodfQYPDXtmmytdElujX3iRrj0u3GZeoyZ78ufS5QFFJCIiIhIQ90S+Zwt3H+fuuwMtgMMIDZrRyN2j7jcc6dDvErD69evTt3/eEZaltIwa+QwTPp7E0Odf3ma+pajgFRFJRieeeCIAb7zxRsCRiJQ+I3Hv2TKzSsAtQEtgOnCfu28sbnsqtkQi0PmEE/j19wXbzFu75h/Kp+risIhIMjr6aHXllySXlZjFFvAY0I7Qs7ROA2oCVxS3sWIVW2a2k7uvLO5ORcqahx56iIceGciopx7joMOPYtGC3/nqk4+5/757gg5NJC4oL0jQvv76K+658w5q7VSDNf/+yx7N9uK+/7s/Zvu76qqrYta2SPxL6OdsHQ/s7+5LzGww8AmxKrbM7FzgL3d/PzydAbwK7GpmvwInuvvs4u5cpCy59uo+LFq0iMmffEKLRg256Owndd+cJB3lBYlHf/75J/fddQejn3iE9PQKAIx/8x2uvfoqHnpkULDBiSQip0zdgxWlyu6+BMDdF5hZ9R1prKgrW9cC5+aaHgpMIDQyx2XAA8CJOxKASFlw/PHHA/Duu+9yVrduAUcjEijlBYk7t91yM/37XZdTaAGc2rkTr783IWb7zJ0XcnN3nho6lNm//EpqajlSU4xrr+7DLrvsErNYRIKQqPdsAeXM7EhCt6blN427fxxxY0Us3w34HsDMGgKtgGPc/W8z6wv8Gk3kImVV586dgw5BJF4oL0jc+evPJTTbo/F282tUqxqzfRaUFwYOepQmLVrR6fSzAVi7Zg39brmVp554nHLldKu8JIqEfs7WUmBErukVeaYdaBJpY0X9r98ClAc2AAcDP7v73+Fl64CKke5IpCy77LLLgg5BJF4oL0jcyWh/AB998jnHHHFozjx3589lsXs8R355YevWrfyxaDGnnntBzrwqVatyfJdTeOeddznxRJ24kwSSoMWWuzcqyfaKGkptMnCPmbUmdGPYm7mW7QX8WZLBiIhI3FNekLjT76abGTT0Gb6cNh2AlatWc2W/2zn1zK6lGse6deuoWq3advMb79GU+b//XqqxiEh8KOrK1lXAc0Av4Avg/3ItOwd4L0ZxicSVY445BoAJE2LX/1+kjFBekLhTrlw5XnvrHW695WYeeeoZ0tLSuOaGvhxwwAEx22d+eaFKlSqsXLFiu3U/m/gRmYcfFrNYREqdO2RtDTqKMqHQYsvdFwFHFbCsb0wiEolDZ555ZtAhiMQF5QWJV+np6Tzw4EOltr/88oKZcVLnE3jsofs596KLqVKlKpM++oClixfQps1FpRabSGnwxH3OVomK+k5NM3vC3XUDiySViy5SkhQpiPKCJKOC8sLRRx/FHns0YfTzo/h33Toyjzic8+66s5SjE4m10r+yZWY7A2OBRsB84Iy8z3cMD9z0LFAXyAKGuvug8LI7gIuAZeHVb3L3d2Idd3GGxelOaHhfERERUF6QBPTPP/8wevRoFi9Zwl7N9+L000+jQoUKRW8INGrUiFtuvinGEYoEyAmiG2Ff4CN3HxAe/bYvcGOedbYA17r7dDOrCnxjZh+6+4/h5Y+4+4OlGHORA2TkR09xlaSTmZlJZmZm0GGIxCvlBUkoixYt4vobbuSAQw7nmhv6UX/3Rlx5VR/Wr1+fs47ygiQzx/GtW6N+7aAuwKjw+1HASdvF5b7E3aeH368BfgLq7+iOd0RxrmzdW+JRiMS5Hj16BB2CSDxTXpCE8uSQIdxyx11UrlwZgDb77keVKlV55plncoZ8V16QpOZA6d+zVcfdl0CoqDKzQp8UbmaNgP2Ar3LN7m1m5wLTCF0BW5nftiUp6mLL3e+LRSAiJWHOnDk0a9asxNtVUhUpmPKCJJotW7NyCq1se+y5J6+OeylnWnlBklux79mqZWbTck0Pdfeh2RNmNoHQ/VZ53RzNTsysCjAe6OPu/4RnPwn0J1Qq9gceAnpG025xFFpsmdnEcEAFcXc/umRDEimeuXPn0rx5c3799Vf22GOPEm178+bNvPHmm3z15ZeklitHxfR0rrjiCnbaaacS3Y9IvFNekGSwZfMm3B2z/3rIrl27lvLl03KmN2/eDEBaWtp224skPHe8eMXWcnfPKLhZP6agZWb2l5nVC1/VqgcsLWC9NEKF1vPu/kqutv/Ktc4w4K3iHEC0irqyNbqA+fWBK4FKJRuOSPGNGzeOlJQUxo8fzw033FCibbdu3RrM+HraNwCsWLGCvv368djgwUq0kmyUFyThHX/88Tz/7Ci6n9cDAHfn8UcH0uvCC3LW6dChAwCTJk0KIEKROFD63QjfAM4DBoR/vp53BQudIXka+MndH86zrF52N0TgZGBWbMMNKeo5W0/nnjazmkA/QsMmjgXuil1oItF5ZdzLXNOzG6+Me7nEi61GjRtzyqmn5UzXrFmT0884k7feeouTTz65RPclEs+UFyQZHH3UUWxYv4G777iVChXS2bBhA2ef1W2bXhMXXnhhgBGKBK3YV7Z2xADgJTO7APgDOB3AzHYFhrt7J+AQ4BzgezObEd4ue4j3+81sX0K9M+YDF5dG0BHds2Vm1YDrgd6ELrnt7+5zYxmYSDT++OMP5s6bx0cjHmT3I05mwYIFNGzYsETa3rp1K23btqVrt27bzG/ZqhUjn366gK1EEpvygiS6//2vE//7X6cCl3fv3r0UoxGJMwEM/e7uK4Dtuqm7+2KgU/j9ZxQwQq67nxPTAAtQ1D1bFYE+wLXAJOBQd/8h9mGJFG7//fZjzi9zcqa3bNnKhWd2oVLFdLp1PpamTZtSrlxqzvJmTZsx/dtvi7Wv1NRUVq9azbp166hU6b8eUp9/9hnt2hXY7VgkISkviISsW7cOYJu8IJI8PIhuhGVSUVe2fgNSgfsJDZFYx8zq5F7B3T+OUWwiBRo6bBjdup7J/i32YODNfahcMZ1KFdMBGHjzVdx7zcX8u34Dfe4ZyPSf5jJs+PAd2t8nn0zmoPbt+Pqb6VSoUIHvv/uOyRM/ZtCgQSVxOCJlifKCCNCpU+iql+7ZkqTklMRzs5JCUcXWBkIXCi8tYLkDTXLPMLOdCfXbb0SoP+QZecewN7OGwLOEhnbMIjTso/5qlYhlZGQw/dsZ9L7sMo465wrGDrqLls1CfelTUlL4beFizrzqNg48+FCmfzueqlWr7tD+brrpJv744w8efeRhtmzdStM99+Shhx7aZqQqkSShvCACXHppQf8FRJJBsYd+TzrmXtgIvsVo0Ox+4G93H2BmfYGd3P3GPOvUA+q5+3Qzqwp8A5zk7j8W1X5GRoZPmzatqNUkifTu3ZtNKxYzpP9/g2Jccuv9VKhVn8GDBwcYmUjJMbNvChsuN57FMi8oJ4hIsgoyL7Rt1sinPHF71Nuld+hZZnNZcaXEoM0uwKjw+1HASXlXcPcl7j49/H4N8BOhYYNFovbzjz9w7KHtWbZiJQ8MHc2yFSvpcEg7fv6x5G4jWb16NatXry6x9kSSjPKCJBzlBUl2npUV9SsZxaLYqpM9hn345y6FrWxmjYD9gK9iEIskuBUrVjB12jeklUtl/5POZ+LMObQ9+XzKp5Xj66nTWLFiRYnsp0uXLnTp0iVneuvWrZT0VWGRBKa8IAknb14QSS7hboTRvpJQREO/52VmEwj1q8/r5ijbqULoCc993P2fQtbrBfQC2G233aLZhSS4119/nc2bN3PZnQ8z8tnn6NChAx988AHn9ziPzZs388Ybb3D++edH3N6mTZv44osvSE9Pp3379jn3ZF155ZUAfPnll4wZ+xKVKldmw4YNNGnciMsvu0z3bknSK828oJwg8SA7L4iIFCYW92zNBjLdfUm4D/4kd2+ez3pphJ7N8n7eJzwXRv3zJbeePc/nryV/8syoUeyyy38ny5cuXUqPc8+hXv36PP30iIjamjx5Mq+9Mp6jM49g3fr1TPrkM6657jr23HNPABYtWsQjgx6l3y235RRXX34xhYXz59GrV6+SPziRXMr4PVsxywvKCSKSrAK9Z6vp7v75oL5Rb1fxf5eV2VxWXMW6slWEN4DzCD3l+Tzg9bwrWOgv1aeBn6IptETyGjZsOCkpKdtdWdpll114+933yIqwf/D69et5/dVXGHD3XTnzOnU8jptuv5NHBz/G8uXLGTp0KBdfdvk2+zrwoIO598MPSuZgRBKX8oIknOXLlwNQq1atgCMRKX2OJ+09WNGKxT1bA4AOZvYL0CE8jZntambvhNc5BDgHOMrMZoRfBT+mXaQAqampBXbhMzNSU1PzXZbXRx99xIknbPtPMC0tjWZ77sGCBQs47bTTGDlyJDvttHO+MYhIoZQXJOGcdtppnHbaaUGHIRIMR/dsRajEr2y5+wrg6HzmLwY6hd9/BugmF4kbqampbM3n4Xxbt2aRmprKtddey/Tp0/l4woccc+xxOcs3bNiggTJEiqC8IIno2muvDToEkQDpOVuRikU3QpEy56ijjuL6a6/hkIMOyrlStmHDBubNn8+uu+7KrrvuygknnMBNN9/M2n/XcvQxx/LbvLk8/9wobrvlloCjFxGR0ta5c+egQxAJjoPnc5JatqdiSwSoUKEC55zXg+v73cyBB7Tn33//5btZP3Bj334A/PnnnwDce889fP3114x5bhQNGzZk8KBBpKenBxm6iIgEIDsv1K2b3yCcIonOQfdsRUTFlgRq06ZNbNq0iSpVqgQdCu3ataNt27bMmjWL9PR0el16ec6yrl27AjBp0iQOOOAADjjggKDCFBGROJA7L4gkJXUjjIiKLQnExo0bufem69m4eD4VU2AFFeh14y20bN0m0LhSUlJo3br1dvP79o1+eFMREUlcyguS1NxxFVsRUbGVhGbPns2YoU9Rrlw5zr70Mho1alTqMdx70/WcXH0jDRuHHrWzNSuL2267kYdfeIVKlSqVejxF6dixY9AhiIhIHFFekGSnod8jE4uh3yWOPT1oECN6nMPhM6dxwLQvGNjtDF5+dlSpxrB582Y2LPqNhjtXy5mXmpJC1xZ1eO3lsaUaS6QWLFjAggULgg5DRCRurVu3jokTJ/L9998nxSityguS1NzxrVlRv5KRrmwlkdWrV/Pty2M5t/p/90edVT2NESOepkvXbpQvX75U4ti0aRPp+ZT5tSpX5McVy0slhmidc845gPrmi4jkZ/zzzzH5+ZG0q+h8t9V4zCpx99AR1K5dO+jQYkZ5QZKZO0lbPEVLxVYS+frLL2m1eQOw7eh5e25czw8//MB+++1XKnFUrlyZVSnpbM3KIjXlv6rrzdmL6dY/PvvA36Lh3UVE8vXnn3/yxeinuarZf6PyHblpM/dd14eHRz0fYGSxpbwgyc3VjTBCKraSSL369fkqZfuvfEVKKrvsskupxnJx39u5/dbrOKN5XWpVTuetOUuo1uaQQO4fi8QxxxwTdAgiInFp3HOj6FKv2jbzKpdPw35fzObNm0lLSwsosthSXpCkpitbEVOxlURatmzJwFq1+WfdP1QLJ7/lGzfy726NqF+/fqnGsnfLljz8wqu8Pv5lflr6F936943bQgtg3rx5ADRp0iTgSERE4ktqaipZ+dyj5VjOQ+ITkfKCJDsVW5FRsZVkHnzuee7ucxX//jYPM6N6830Y8MjAQGKpWLEiXbufG8i+o9WzZ09AffNFRPI67ZzzuOect7iiWuWceas3bCKlbgPKlUvcPzOUFySZuTtZWzX0eyQS97eg5KtGjRo8OLJ0Rx9MBHfeeWfQIYiIxKXatWvT4eIreWjYE7SukMXfW2BRpZ2456lHgw4tppQXJNnpnq3IqNgSicARRxwRdAgiInHrfyefwrEndGbmzJnUrFmTxo0bBx1SzCkviEgkVGyJRGD27NkANG/ePOBIRETiU1paGhkZGUGHUWqUFySphZ+zJUVTsSUSgYsvvhhQ33wREQlRXpBkV9rFlpntDIwFGgHzgTPcfWU+680H1gBbgS3unhHN9iVNxZZIBO69996gQxARkTiivCDJzD2Q52z1BT5y9wFm1jc8fWMB6x7p7st3YPsSo2JLJAIHH3xw0CGIiEgcUV6QZJdV+t0IuwCZ4fejgElEVyzt6PbFomJLJAKzZs0CQs8qExERUV6QpBbMQ43ruPsSAHdfYma7FLCeAx+YmQNPufvQKLcvUSq2RCLQu3dvQH3zRUQkRHlBklrxB8ioZWbTck0PzVUMYWYTgLr5bHdzFPs4xN0Xh4upD83sZ3f/pDjBlgQVWyIReOCBB4IOQURE4ojygiQzp9jP2VqePWBFvu26H1PQMjP7y8zqha9K1QOWFtDG4vDPpWb2KtAe+ASIaPuSpmIrAaxbt46Hbr+VFXN+xlNTaX3kMfTsfQVmFnRoCaNdu3ZBhyAiInFEeUGSWjBDv78BnAcMCP98Pe8KZlYZSHH3NeH3xwJ3Rbp9LKjYKuPcnavP7sYZvoZdKqUDW5jx7nge+OtPbrhbIyWVlBkzZgCw7777BhqHiIjEB+UFSXYBFFsDgJfM7ALgD+B0ADPbFRju7p2AOsCr4QsO5YAX3P29wraPNRVbZdyXU6bQYs1ydqldPWfevtUrM+2rz1m/fj0VK1YMMLrE0adPH0B980VEJER5QZKaQ1YpD/3u7iuAo/OZvxjoFH4/D2gTzfaxpmKrjJv9/XfsUWH7r7G2b2X58uU0bNgwgKgSz8CBA4MOQURE4ojygiQzJ5BuhGWSiq0y7sAjMhn38vM0qrbt/CWpFahXr14wQSUgdRMREZHclBckqTn41q1BR1EmpAQdgOyYvVq04N9m+zD1739wdzZnZfHq0tUcfNqZlCunWrqkTJ06lalTpwYdhoiIxAnlBUlujmdlRf1KRvprPAHc+8QQXhnzIs+++zYp5ctz2tV30P7AA4MOK6Fcf/31gPrmi4hIiPKCJLVgHmpcJqnYSgBmxqndzuLUbmcFHUrCeuyxx4IOQURE4ojygiQ33bMVKRVbIhFo2bJl0CGIiEgcUV6QZOYOWSq2IqJ7tkQiMGXKFKZMmRJ0GCIiEieUF0QkErqyJRKBm266CVDffBERCVFekOTmSTvgRbRUbEmB3J3ly5dTqVIlKleuHHQ4gXrqqaeCDkFEROKI8oIkNQ2QETEVW5KvaV9+yVN33Eattf+wzlKp0LQZdz7+BBUrVgw6tEA0b9486BBERCSOKC9IUnPwrR50FGWC7tmS7axZs4YnrruaXraFU6tV5pyq6Rwx/xfuvPKKoEMLzOTJk5k8eXLQYYiISJxQXpBk5jhZW7OifiUjXdmS7Yx77lk6pjhmljOvTnoF1sz+iS1btiTlw5Jvv/12QH3zRUQkRHlBkpqDZ+nKViSS769mKdLqlSupm09BVS4rK2mLrREjRgQdgoiIxBHlBUl2WepGGBF1I5TtnHR2dyZu2LzNvC1ZWWzcuRbp6ekBRRWsJk2a0KRJk6DDEBGROKG8IMnMwwNkRPtKRiVebJnZzmb2oZn9Ev65UyHrpprZt2b2VknHIcXXqFEjGp98Ki+s/pff/13HzNVreWLdFq578KGgQwvMhAkTmDBhQqHruDvffvst7777LqtWrSqdwETKAOUFSUSR5AWRhOWOb43+lYxi0R+sL/CRuw8ws77h6RsLWPcq4CegWgziSDjuzgfvvM3EN1+nYuUqnH3p5ey5554x2del19/Aku7n8O4r46lVty5Pn9iFtLS0mOyrLLj77rsBOOaYY/JdvnLlSm659Tb2y8igQYOGPPjwIzTefTcuuOCC0gxTJF4pL0jCKSoviCQ6dSOMTCyKrS5AZvj9KGAS+SRVM2sA/A+4B7gmBnEknNv6XEGTP+dyUYOarNu8huFX9SLzkqvo2LlLTPZXr149el7eOyZtlzXPPfdcocsfeOBBrr2hLzV2Cp2wb3/AgQwfOoSffvqJFi1alEaIIvFMeUESTlF5QSSh6TlbEYvFPVt13H0JQPjnLgWsNxC4AdA3FYFZ339P9QWzyWxYixQzqpRP4+IW9XlrxFDcdWYh1ho2bEjDhg0LXL5569acQitbt7PPYfz48bEOTaQsUF6QhFNUXhBJZA5kZXnUr2RUrCtbZjYBqJvPopsj3P4EYKm7f2NmmRGs3wvoBbDbbrtFHmgC+fSD9zi4dtXt5tfO2sjKlSvZeeedA4gqebz33nsAdOzYMf8V8il4N2/aRFpa+ViGJRI3SjMvKCdIPCgyL4gkMk/ee7CiVaxiy90L7KBsZn+ZWT13X2Jm9YCl+ax2CHCimXUC0oFqZjba3bsXsL+hwFCAjIyMpPxmGzVrzu/fTKRe1crbzF/tqVStun0RJiVrwIABQMFJtXLlSixevJhdd901Z96oZ0bQs8e5pRKfSNBKMy8oJ0g8KCoviCS6ZH1IcbRicc/WG8B5wIDwz9fzruDu/YB+AOEzmNcVVGhJyLGd/kevoU+wz8ZNVK0Quloyc9lqdt2vfVIPXFFaxowZU+jyG2+4gdtuv4Oda9Wibr16fDdjBkccdhi77757KUUoEteUFyThFJUXRBJZaOh3neuKRCyKrQHAS2Z2AfAHcDqAme0KDHf3TjHYZ8JLTU3l/pHP8+Atfdm0aDFbU1JpduBhXHvtdUGHlhTq1s2vd9R/KlasyAP3/x9//vkny5Yto3vXM1UEi/xHeUESTlF5QSShqdiKWIkXW+6+Ajg6n/mLge0SqrtPIjQylRShZs2a3PfksKDDSEpvvvkmAJ07dy50vbp16yoBi+ShvCCJKNK8IJKYXN0IIxSLK1siCeehh0IPdFZSFRERUF4Qkcio2BKJwLhx44IOQURE4ojygiQ1B0/SodyjFYvnbIkknFq1alGrVq2gwxARkTihvCDJzIGsrR71a0eY2c5m9qGZ/RL+uVM+6zQ3sxm5Xv+YWZ/wsjvMbFGuZaVyv7CKLZEIvPLKK7zyyitBhyEiInFCeUGSmju+NSvq1w7qC3zk7k2Bj8LTecLy2e6+r7vvC7QF1gGv5lrlkezl7v7OjgYUCXUjFInAo48+CsApp5wScCQiIhIPlBck2QUwGmEXIDP8fhShgZRuLGT9o4G57v57bMMqnIotkQi8/vp2jwUSEZEkprwgycydHe4WWAx13H1JaP++xMx2KWL9rsCLeeb1NrNzgWnAte6+MgZxbkPFlkgEqlevHnQIIiISR5QXJNl5VrG6BdYys2m5poe6+9DsCTObAOT3DJ2bo9mJmZUHTgT65Zr9JNCf0C1n/YGHgJ7RtFscKrZEIjB27FgAzjzzzIAjERGReKC8IEnNiz3gxXJ3zyi4WT+moGVm9peZ1Qtf1aoHLC1kP8cD0939r1xt57w3s2HAW9GFXjwqtkQi8OSTTwJKqiIiEqK8IEnNA7ln6w3gPGBA+GdhfXm7kacLYXahFp48GZgViyDzUrElEoF33imVAWtERKSMUF6QZOZQEqMLRmsA8JKZXQD8AZwOYGa7AsPdvVN4uhLQAbg4z/b3m9m+hMKfn8/ymFCxJRKBSpUqBR2CiIjEEeUFSWoBDJDh7isIjTCYd/5ioFOu6XVAzXzWOyemARZAxZZIBEaPHg1A9+7dA45ERETigfKCJDcPohthmaRiSyQCw4cPB5RURUQkRHlBkpk7ZLmKrUio2BKJwIcffhh0CCIiEkeUFyTZbVWxFREVWyIRSEtLCzoEERGJI8oLkswcUC/CyKjYEonAyJEjAejRo0egcYiISHxQXpBkpytbkVGxJRIBJVUREclNeUGSma5sRU7FlkgEJk2aFHQIIiISR5QXRCQSKrZERERERCRi7upGGCkVWyIiIiIiEhV1I4yMii0REREREYmY47qyFSEVWyIiIiIiEjENkBE5FVsiIiIiIhIVFVuRUbElIiIiIiIR0wAZkVOxJSIiIiIiUdGVrcio2BIRERERkYiF7tlStRUJFVsiIiIiIhIxDZARORVbIiIiIiISFV3ZioyKLRERERERiVhogIygoygbVGyJiIiIiEhUdGUrMiq2REREREQkYg5kBR1EGaFiS0REREREouC6shUhFVsiIiIiIhIxjUYYORVbIiIiIiISMT1nK3IqtkREREREJHIajTBiKUEHICIiIiIikohKvNgys53N7EMz+yX8c6cC1qthZuPM7Gcz+8nMDirpWEREJHjKCyIiiSW7G2G0rx1hZqeb2Q9mlmVmGYWs19HMZpvZr2bWN9f8iHJRSYvFla2+wEfu3hT4KDydn0HAe+6+F9AG+CkGsYiISPCUF0REEsxWj/61g2YBpwCfFLSCmaUCjwPHA3sD3cxs7/DiSHNRiYpFsdUFGBV+Pwo4Ke8KZlYNOBx4GsDdN7n7qhjEIiIiwVNeEBFJIEFc2XL3n9x9dhGrtQd+dfd57r4JGEMoB0EEuSgWYlFs1XH3JQDhn7vks04TYBnwjJl9a2bDzaxyDGIREZHgKS+IiCSQ7KHfS/nKViTqAwtyTS8Mz4PIclGJK9ZohGY2Aaibz6Kbo9jv/sAV7v6VmQ0idCnv1gL21wvoFZ5ca2ZFVbUAtYDlEcYTzxLhOBLhGEDHEW+S7Th2j3UgO6I084Jygo4jjug44kuyHUdgeWE5m95/it9rFWPTdDOblmt6qLsPzZ4oLJe4++sRtG/5zAt03MRiFVvufkxBy8zsLzOr5+5LzKwesDSf1RYCC939q/D0OArpNxn+EoYWtLyAOKa5e4E3z5UViXAciXAMoOOINzqO+FKaeUE5QccRL3Qc8UXHUXrcvWOM2i0wl0RoIdAw13QDYHH4fSS5qMTFohvhG8B54ffnAdtVoe7+J7DAzJqHZx0N/BiDWEREJHjKCyIiUhqmAk3NrLGZlQe6EspBEEEuioVYFFsDgA5m9gvQITyNme1qZu/kWu8K4Hkz+w7YF7g3BrGIiEjwlBdERGSHmNnJZrYQOAh428zeD8/PySXuvgXoDbxPaETbl9z9h3AT+eaiWCtWN8LCuPsKQmck885fDHTKNT0DiOUl0qi6mMSxRDiORDgG0HHEGx1HGREneSFRPmcdR3zRccQXHUcCc/dXgVfzmZ83l7wDvJPPevnmolgz38FhGEVERERERGR7sehGKCIiIiIikvQSotgys9PN7AczyzKzArugmNl8M/vezGbkGXYyLkRxHB3NbLaZ/WpmpfL062iY2c5m9qGZ/RL+uVMB68Xl91HU52shj4aXf2dm+wcRZ1EiOI5MM1sd/vxnmNltQcRZFDMbYWZLzWxWAcvj/vuI4BjKxHdRligvxBflhfiQCHkhEXICKC8kFXcv8y+gBdAcmARkFLLefKBW0PHuyHEAqcBcQg8ALQ/MBPYOOvY8Md4P9A2/7wv8X1n5PiL5fAn1C36X0LMcDgS+CjruYh5HJvBW0LFGcCyHE3r+0qwClpeF76OoYygT30VZeikvKC+UYOzKC3H0SoScEOFxxP13oVdkr4S4suXuP7l7JA+1jGsRHkd74Fd3n+fum4AxQJfYRxeVLsCo8PtRwEnBhRK1SD7fLsCzHvIlUMNCz2uIJ2Xh30lE3P0T4O9CVon77yOCY5ASprwQd//flReCVxb+nRQpEXICKC8kk4QotqLgwAdm9o2Z9Qo6mGKqDyzINb0wPC+e1HH3JQDhn7sUsF48fh+RfL5l4TuINMaDzGymmb1rZvuUTmglrix8H5FIhO+iLIrH30PRKgv/B5QXgpcseaEsfBeRKuvfhRCDod9jxcwmAHXzWXSzu0f6ULJD3H2xme0CfGhmP4fPLJSaEjgOy2deqQ8pWdhxRNFM4N9HPiL5fOPiOyhCJDFOB3Z397Vm1gl4DWga68BioCx8H0VJlO+iVCkv/NdEPvOUF0qO8kLZUha+i0gkwnchlKFiy92PKYE2Fod/LjWzVwldUi/VX+IlcBwLgYa5phsAi3ewzagVdhxm9peZ1XP3JeFL90sLaCPw7yMfkXy+cfEdFKHIGN39n1zv3zGzJ8yslrsvL6UYS0pZ+D4KlUDfRalSXsgRF/8HlBeC/w6KkCx5oSx8F0VKkO9CSKJuhGZW2cyqZr8HjgXyHQEmzk0FmppZYzMrD3QF3gg4przeAM4Lvz8P2O7MbBx/H5F8vm8A54ZHPDoQWJ3dPSaOFHkcZlbXzCz8vj2h3wcrSj3SHVcWvo9CJdB3UabE8e+haCkvxJbyQtlSFr6LIiXIdyGQMKMRnkzoTMZG4C/g/fD8XYF3wu+bEBp5ZybwA6HuGYHHHu1xhKc7AXMIjSoUj8dRE/gI+CX8c+ey9H3k9/kClwCXhN8b8Hh4+fcUMtJZnB9H7/BnPxP4Ejg46JgLOI4XgSXA5vD/jwvK2vcRwTGUie+iLL2UF4KPPc9xKC/EwSsR8kIi5IQIjyPuvwu9IntZ+AsVERERERGREpQ03QhFRERERERKk4otERERERGRGFCxJSIiIiIiEgMqtkRERERERGJAxZaIiIiIiEgMqNiSpGJm75rZeUWvme+2d5jZ6JKOKdz2i2Z20g628bWZ7VNCIYmIJAXlBRGJJRVbEjgzm29m681sba7XY/msV8/M3Mzq5Jp3cwHz3stvX+5+vLuPis2RFI+ZtQbaEH7IZ/hBjMPMbJKZDcn1UMPa4SfIzzezlWb2jZldb2bp4aYeBO4K5ihEREqO8oLygkiiULEl8aKzu1fJ9eqddwUPPQH+V+DwXLMPB37OZ94nMY22ZF0MPO//PfTuOGCpu2cSelr8sWZWA5gC1AVOAOoB3YEWQPPwdm8AR5pZvdILXUQkZpQXlBdEyjwVW1LWfEI4gZpZKrAfMCjPvIMoIKmGzwpeGH7fw8w+M7MHw2cEfzOz43Ot29jMJpvZGjP7EKiVp60DzWyKma0ys5lmlhmef7CZLTezhuHpNuF19irgmI4HJuea3gjUMbM0Qkl0PXA1oafMn+7us9x9g7v/5O49ge8B3H0D8A1wbASfo4hIolBeUF4QiVsqtqSsyUmqhBLqz8BHeealAV9H2N4BwGxCCfN+4Ons7hnAC4SSVC2gP5DTp9/M6gNvA3cDOwPXAePNrLa7TwGeAkaZWUXgOeAWd/85787NrDLQOBwDAO4+EVgYPta57v4JoUQ5zt235m3D3bNyTf5EqOuJiEiyUF7IQ3lBJH6o2JJ48Vr4LF/266IC1psMtDSznYDDgE/d/RegVq55X7r7pgj3+7u7Dwsnq1GEumHUMbPdgHbAre6+MZzY3sy1XXfgHXd/x92z3P1DYBrQKbz8DqA6oeS+GHi8gP3XCP9ck3umu9/h7ge5+73hWbXC7RRlTa42RUTKMuWFXJQXRMomFVsSL05y9xq5XsMAzOyHXDdHH+bu8wmd3TuU0FnLT8Pbf5FrXjT98v/MfuPu68JvqwC7Aivd/d9c6/6e6/3uwOm5/xAI779euK3NwEigJfBQrn73ea0K/6xaRJwrstsuQtVcbYqIlGXKC4VTXhApA1RsSVxz931y3RydnUA/JZQ8DyJ0c3DueYdSMjdBLwF2CnfnyLZbrvcLgOfy/CFQ2d0HQE53ktuBZ4CHzKxCAcf3LzAXaFZEPBOAk82sqP+zLYCZRawjIlJmKS/kUF4QKQNUbElZ9AlwLrDY3f8Jz/ssPK86obOZO8TdfyfU/eNOMytvZocCnXOtMhrobGbHmVmqmaWbWaaZNQj37R8JPA1cQChB9y9kd+8ARxQR0sNANWC0me1hZilmtruZPWxmrQDCibst8GH0RywiUqYpLygviMQlFVsSL960bZ+n8moh604GdiGUSLPNACoC3+Tq9rGjziJ0o/TfhM5GPpu9wN0XAF2Am4BlhM5oXk/o/9SVQB1C/fodOB8438wOK2A/Q4Gzc92AvR13/xs4mNCIVFOAtcD7wEpCwx4DnAhMcvdI+vCLiMQ75QXlBZEyzwruMiwipcXMXgBecvfXdqCNr4AL3H1WiQUmIiKBUF4QSQwqtkRERERERGJA3QhFRERERERiQMWWiIiIiIhIDKjYEhERERERiQEVWyIiIiIiIjGgYktERERERCQGVGyJiIiIiIjEgIotERERERGRGFCxJSIiIiIiEgP/D+IyHJGzPE3IAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1008x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1,2, figsize=(14,5))#, constrained_layout=True)\n",
    "\n",
    "for ax in axs:\n",
    "    ax.plot([-1.5,1.5],[0,0],':',c='k')\n",
    "    ax.plot( [0,0],[-.5,.5],':',c='k')\n",
    "    ax.set_xlim(-1.6, 1.6) \n",
    "    ax.set_ylim(-.6, .6) \n",
    "    ax.set_xlabel(\"E-W index (℃)\",fontsize=12)\n",
    "    ax.set_ylabel(\"N-S index (℃)\",fontsize=12,labelpad=-8)\n",
    "\n",
    "ax = axs[0]\n",
    "expid = 'CESM2 Ensemble Mean'\n",
    "\n",
    "im = ax.scatter(lst_cntl_EW, lst_cntl_CE, c=lst_cntl_PCC, cmap='RdBu',vmax=1,vmin=-1, edgecolors='k',linewidths=.5)\n",
    "ax.scatter(x=EW_cntl,y=CE_cntl, c=PCC_cntl, label=expid, marker='*', cmap='RdBu',vmax=1,vmin=-1, edgecolors='k',s=160, linewidths=1)\n",
    "ax.scatter(x=[data_hadi[-1]],y=[data_hadi[-2]], c=[data_hadi[-3]],  marker='D', cmap='RdBu',vmax=1,vmin=-1,edgecolors='k',s=50)\n",
    "ax.text(x=data_hadi[-1]+.1,y=data_hadi[-2]-0.02, s='HadISST',fontsize=13)\n",
    "ax.set_title('CESM2 (n=15)',fontsize=14)\n",
    "ax.legend(loc='upper left',fontsize=10)\n",
    "ax.set_box_aspect(.85)\n",
    "\n",
    "ax = axs[1]\n",
    "expid = 'CESM2-FA Ensemble Mean'\n",
    "im = ax.scatter(lst_fa_EW, lst_fa_CE, c=lst_fa_PCC, cmap='RdBu',vmax=1,vmin=-1, edgecolors='k',linewidths=.5)\n",
    "ax.scatter(x=EW_fa,y=CE_fa, c=PCC_fa, label=expid, marker='*', cmap='RdBu',vmax=1,vmin=-1, edgecolors='k',s=160, linewidths=1)\n",
    "ax.scatter(x=[data_hadi[-1]],y=[data_hadi[-2]], c=[data_hadi[-3]],  marker='D', cmap='RdBu',vmax=1,vmin=-1,edgecolors='k',s=50)\n",
    "ax.text(x=data_hadi[-1]+.1,y=data_hadi[-2]-0.02, s='HadISST',fontsize=13)\n",
    "cbar = fig.colorbar(im, ax=axs, orientation='vertical', fraction=0.046, pad=0.04)\n",
    "cbar.set_label('PC index', size=12)\n",
    "ax.set_title('CESM2-FA (n=15)',fontsize=14)\n",
    "ax.legend(loc='upper left',fontsize=10)\n",
    "ax.set_box_aspect(.85)\n",
    "\n",
    "plt.savefig(\"Fig2_trend_indices_1970-2020_ns5deg.pdf\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}