{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "daef28b7-0e39-459d-93f5-99cbb96f4c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import datetime\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "\n",
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import gsw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "821de9fd-b284-447b-a377-0bdeb93b86fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path to ITPs data and topography in the Canada basin file\n",
    "# ITP data from https://www2.whoi.edu/site/itp/data/data-products/\n",
    "path='XXX/Arctique/ITP/'\n",
    "file_topo='XXX/Arctique/etopo4.mat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "532cef47-3fcb-401d-a072-d2c7686a18a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_itp(file):\n",
    "    \"\"\"Read an ITP file and returns a xarray dataset with time and z-levels as dimensions/coordinates\n",
    "    and the number of the ITP.\n",
    "    \"\"\"\n",
    "    ds = xr.open_dataset(file,decode_times=False)\n",
    "    ds = ds.assign_coords({'lon':ds['LONGITUDE'],'lat':ds['LATITUDE'],'levels':ds['BIN']})\n",
    "    t0 = ds.REFERENCE_DATE_TIME.values.tolist().decode('utf-8')\n",
    "    ds = ds.where(ds.JULD<700000,drop=True)\n",
    "    ds = ds.assign(time=pd.to_datetime(ds.JULD.values,origin='-'.join([t0[:4],t0[4:6],t0[6:8]]),unit='D'))\n",
    "    ds = ds.set_index(N_PROF='time')\n",
    "    ds = ds.set_index(N_LEVELS='levels').rename({'N_PROF':'time','N_LEVELS':'levels'})\n",
    "    id_itp = file.split('_')[1]\n",
    "\n",
    "    return ds,id_itp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "179b12f6-382d-4088-bb28-a6b48b932ec3",
   "metadata": {},
   "source": [
    "### Geographical area and period to study"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ef11471a-f44f-43a4-b26e-01eff4e55eb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# geographical zone:\n",
    "lonmax=-120\n",
    "lonmin=-160\n",
    "latmin=65\n",
    "latmax=85\n",
    "\n",
    "# time period:\n",
    "year0=2003\n",
    "year1=2007\n",
    "#year0=2008\n",
    "#year1=2013\n",
    "#year0=2014\n",
    "#year1=2019\n",
    "\n",
    "month0=1\n",
    "month1=12\n",
    "\n",
    "# depth range\n",
    "p1 = 300\n",
    "p2 = 700"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98fabb52-8bcb-4789-be12-bedfaeae1347",
   "metadata": {},
   "source": [
    "### Selection and creation of a list of sub dataset for each itp\n",
    "We stock our sub dataset in a dictionnary, list_itp, with the id of the ITPs as keys.\n",
    "We read each file, select the needed time period when possible anse select the data in the good area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "30b366a5-6ba6-452a-baad-8cbbbf52f5ef",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ooops\n",
      "Ooops\n",
      "Ooops\n",
      "Ooops\n",
      "Ooops\n",
      "Ooops\n"
     ]
    }
   ],
   "source": [
    "list_itp={}\n",
    "for file in glob.glob(path+'*L3*.nc'):\n",
    "    ds, id_itp = read_itp(file)\n",
    "    try :\n",
    "        ds_sub = ds.sel(time=slice(datetime.datetime(year0,month0,1),datetime.datetime(year1,month1,30)))\n",
    "    except:\n",
    "        print('Ooops')\n",
    "    else:\n",
    "        ds_sub = ds_sub.where(np.logical_and(ds_sub.lon >= lonmin, ds_sub.lon <= lonmax),drop=True)\n",
    "        ds_sub = ds_sub.where(np.logical_and(ds_sub.lat >= latmin, ds_sub.lat <= latmax),drop=True)\n",
    "        if ds_sub.levels.size==0:\n",
    "            pass\n",
    "            #    print(file)\n",
    "        else:\n",
    "            if ds_sub.time.size==0:\n",
    "                pass\n",
    "            else:\n",
    "                list_itp['%s'%str(id_itp)]= ds_sub\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f9aa66d6-ebfe-4239-bfca-40bb6b535adf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of ITPs :  10\n",
      "--------------------------------------\n",
      "Grid levels step (for each itp):  [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "print('Number of ITPs : ',len(list_itp.keys()))\n",
    "print('--------------------------------------')\n",
    "print('Grid levels step (for each itp): ',[list_itp[itp].levels.diff(dim='levels').mean().values.tolist() for itp in list_itp])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "134631ed-5e30-4c87-8312-54f16a787b8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# from itp\n",
    "DS = {}\n",
    "for itp in list_itp:\n",
    "    ds = list_itp[itp]\n",
    "    Lon, Lat = np.meshgrid(ds.lon.values,ds.levels.values)[0].transpose(),np.meshgrid(ds.lat.values,ds.levels.values)[0].transpose()\n",
    "    SA = xr.DataArray(gsw.SA_from_SP(ds['PSAL_ADJUSTED'].values,ds['PRES_ADJUSTED'].values,Lon,Lat),dims=['time','levels'])\n",
    "    CT = xr.DataArray(gsw.CT_from_t(SA,ds['TEMP_ADJUSTED'].values,ds['PRES_ADJUSTED'].values),dims=['time','levels'])\n",
    "    rho = xr.DataArray(gsw.sigma0(SA,CT),dims=['time','levels'])\n",
    "    DS[itp] = xr.Dataset({'SA': SA,'CT':CT,'dens':rho,'pres':ds['PRES_ADJUSTED']},\n",
    "                         coords={'time':ds.time.values,'levels':ds.levels.values})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26614126-d504-4a6d-b0dd-d0810a3b839e",
   "metadata": {},
   "source": [
    "### Compute steps from ITP data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4b090f90-d770-4db7-9d97-9f199b138d46",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold_rho=5e-4 #threshold used for the definition of steps on the density profiles\n",
    "th=10  #Minimal thickness of steps, under this value we do not count th steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "39236119-9187-443b-8483-4d6136fd488b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimal value for which we consider the gradient as the signature of a step :  0.0005\n",
      "------------------------------------------------------------------------------\n",
      "Minimum thickness of steps (dbar) :  10\n"
     ]
    }
   ],
   "source": [
    "print('Minimal value for which we consider the gradient as the signature of a step : ',threshold_rho)\n",
    "print('------------------------------------------------------------------------------')\n",
    "print('Minimum thickness of steps (dbar) : ',th)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "448fbf3d-6625-4bbc-9eb6-f5e328080a6a",
   "metadata": {},
   "source": [
    "#### Definition of density steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4e4ef048-5241-4b41-8c12-9b9c5ef240f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_test = DS[list(list_itp.keys())[0]]\n",
    "dsi = ds_test.sel(levels=slice(p1,p2))\n",
    "\n",
    "var = dsi.dens.isel(time=114)\n",
    "\n",
    "drho = var.differentiate('levels')\n",
    "rho_th = var.where(drho>threshold_rho,drop=True)\n",
    "rho_g = rho_th.where(rho_th.levels.diff(dim='levels')>th,drop=True)\n",
    "rho_g = xr.concat([rho_th[np.isin(rho_th.levels.values,rho_g.levels.values)],rho_th[np.roll(np.isin(rho_th.levels.values,rho_g.levels.values),-1)]],dim='levels')\n",
    "rho_g = rho_g.sortby('levels')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "5c34f27b-090a-4745-b5dc-0f64929b3847",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x7f910a42d5e0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABNH0lEQVR4nO2dd3gc1dWH37OrZkmWq2yMXHHB2GAc25jeIWAgOIFQAwSnEBNIQhqBkHyBwBcIJCThC8EQSuglIRAHTAdjDMa4YNyL3GXLtuQiWV27e78/ZmY1u9oqaXcl+7zPs8/uzsydOXNn9/7mnHPvHTHGoCiKoiiJ4sm0AYqiKErXQoVDURRFSQoVDkVRFCUpVDgURVGUpFDhUBRFUZJChUNRFEVJChWOMERksIjUiIg307YoipJaRMSIyIhM29HVOOiFQ0Q2ichZzndjzBZjTKExxp9Ju6IhIrki8piIbBaR/SLyuYhMCdvmTBFZLSJ1IvKBiAxxrRMR+b2I7LZf94qIuNZvEpF6WzxrROTtOPYMtY9RZx/zLNc6EZHbRGSLiFSLyAsiUhRjXwNEZKaIbLf/0EPD1t8rIlvtfW0Wkdvs5Se77K21y9a4XoNF5FIR+cS2c3aEY58hIovtfW8QkevinPePRWSHiFSJyOMikuta11tEXrFt2SwiV8bZV6zrdZNtT7VdL38Skaw27ivetY96Le31V9rnUysir4pIb9e6XLsequ16+YlrXV8R+dg+5j4RmSciJ4bZdZeIbLPrc7aIjI1wbiNFpEFEnkn0nLs6InK6fU5VIrIpxnan2r/7u9JimDHmoH4Bm4CzMm1HEvYWALcDQ7GE/wJgPzDUXt8XqAIuAfKA+4BPXeW/B6wBBgIlwEpgelvrA5gH3A90Ay4G9gHF9rpvAquBQUAh8B/gyRj76g98HzgeMM45udYfDhTYn0uAFcBFYdsMtctmhS0/C7gU+B9gdti6bLvOvgcIcAxQAxwdxc5zgJ3AWKAXMBu4x7X+eeBF+5xPsvc9Nsq+4l2v4UBP+3Nv4H3gJ23cV7xrH+tajrV/Z6fY5/Uc8IKr7N3AR3Z9HAHsAM611+XZ185j1+9XgT3ONbKvy3bgMMBr72txhPN72z7GM4mecwK/XwOMyPT/OoZ9k4GrgeuATVG2yQaWAJ8Cd6XFrkxXTIYvytNAAKi3G4qbCWt47EbhLuATe5v/An2AZ4FqYAGuBg4YDbxj/zHWAJem4TyWAhfbn68DPnGtK7DPb7T9/RPgOtf6b4c1LptIUDiAUUAj0N217COnMQL+Bfzcte4EoAHIj7PfLCIIR9g2JcAy4Oaw5SHXL0K579BaOPrbZfJdyxYAV0TZx3PA71zfzwR2uOq7CRgV9ju7J8q+Yl6vsG37AO8Cf2vLvmJd+wSu5e+A51zrhtvn2d3+vg34smv9nbiExbXcA3zFru9+9rJfAC+5thkLNISVuxx4Ceum6ZlEzzmB33BQOIAewFNABbAZ+BXgsdddC8wF/gDsBTYCU9rzv03mhXXjsynKuluAe4F/kCbhOKhDVcaYq4EtwFeMFZ66N8qml2OpfgnWH2Ye8ATWHeAq4DcAIlKAJRrPAf2AK4C/RXK77e3/ZrvukV5LEzkHEemP9adfYS8aC3zhOsdaYL29vNV6+3O4fc+KSIWIvC0iR8c4/FhggzFmf5T9if3C9T0XGBnvvKIhIreISA1QhtVIPNfWfTkYY3ZieQnTRMQrIscDQ7AaCkTkJBHZ5yoSqQ77i0gfrGvhN8asDVsfrGP7+p4UaV8RrpcTIqoGKoGjgYejnEp7rn28axm+7/XYAikivYBDY+zbOY+lWDcOM4FHjTG77FUvACNEZJSIZGN5qm+6yhUBvwV+2oZzTob/wxKPw4BTgWuAaa71x2LdDPbFaqgfc4f63IjIazH+26+1wbaI2GG5b2HVT9o4qIUjCZ4wxqw3xlQBbwDrjTHvGmN8wD+BL9nbXYB1V/CEMcZnjFkMvAx8PdJOjTHfN8b0jPIaF88o+0/2LFb4Z7W9uBDLdXdTBXSPsr4KKHT9Ab6Bddc+BPgAeEtEekYxId6x3gC+Y8fOe2DdWQLkxzu3aBhj7rH3PwHrTj78+G3leawwViPWnfZtxpit9jHnGmN6uraNVIfYdsWrE+zrOzfKviJt/5wxpghLlGZghcki0Z5rn2xZ9/pC1/eI52CfxzigCLgSW5RtyrHqfA2Wt3AJ8GPX+juBx5zrEUbc+ksEsTrDXAbcaozZb4zZBPwR64bRYbMx5u/Gyn8+CQzA8lZbYYy5IMZ/+4JkbIvDA8CvjTE1HbjPuKhwJIb7j1of4bvzxxkCHOu+u8BqiA/paINExIPVcDYBN7pW1WD9Od0UYcWnI60vAmqM47cb87Expt4YU2eMuRsrzn2yfcwV0pJwPjmBYz2O1SDPxvKIPrCXl0loQnsFSWAsPseq+zuSKRsJERmNlZO4BsjBulu9WUTOj1IkUh2Cdd7x6iTevqJub4xZh1WPf2vjvmJd+2TLutfXuL7HO4cGY8zzwC0ub/Y3WHmlQVh5ijuA90UkX0TGY4Vp/hThfOPZlQx9sa79ZteyzVhRBocdrvOosz8W0kGIyC9d/4kZCWz/FaxQ4YsdZUOiqHBYMc6OYivwYdjdRaEx5vpIG4vIDAnt/VOTSGNq3yE+hnW3c7Exptm1egVWOMPZtgArvLYi0nr7c6yG22CHm4wxY+3zKTTGfGSXO0xE3Hd3wf0ZYwLGmN8YY4YaYwbay7cB24wxH7n21ZawAli5kOFtLOvmSGCNMeYt2+Y1wOvAlCjbR6rDncaY3cBaIEtERoatj1bH8a5XOLHOuT3XPua1jLDvw7DCjmuNMXuxvIZkflfZWCEhZ9sXjTFltqf+D6wk+xjgNCwPeIuI7AB+BlwsIosTPOdEqQSasW7+HAZj/V6TRkTeiPHffiNSGWPM71z/iekJHOZMYJJYvdh2YHlMN4nIf9pic1KkI5HSmV9YPRHcCcOhtE6Of8e1/i7gH67vZwGl9ufuWHcpV2P9MbKx7qSO6GCbZ9h2F0ZYV4zlql+Mdff2e0KT39Ox8jIlWHHpFbQkQAcDJ2LdeeUBP8dKFPaJU39/sLf/GqE9cXpj/YkFqxFY7q7rKPvLw8pdGKyeOHn2cg9Wr6Be9v4mYzVWPwwrH3L9XMu99r6nA3Psz9n2uuFYd65n2PseDpQC341i47lYd59jbHveJ7RX1QtYnlaBXZ+xelXFu17foSWJPMa+Xve3cV9Rr30C13IsVmeQk+3zeobQXlX3AB/a9THavjZOr6rjsHqX5WD12PoFlkdwqL3+N1ihq/72db4aqAV6YoU1D3G9/oDV6aI4wXO+lihJZXu9Ozn+DPAK1v94CFaPwO+49jM3WtkUtk8e+7ymYLUteUCOq71x182LWJ5Z71TaZIxR4QCmYiXI92HdzQyljcJhfz8c6261AtiN1aiM70B7h9j2NWA1ds7rG2E2rcYK5cwmtNeXYCX29tivewGx143F6qFVa9v+HjApjj1D7WPUY8Woz3KtG2Uvq7N/9BG7kYbtz4S/7OUerITpHvt81wK/dGwPsyeScFwbYd/u63gplrDtx0q8/56WHjUnY4V03Pv7CVbIshqro0Sua11v4FW7HrcAV4aVrQFOTvB6PWEfpxarx9t92GJqr1/REdc+3rW0119pn08tVtfq3q51uVihyWrb3p+41p2KlcDebx/3Q+AU1/o84EEssakGFmOLToTfx+24elUlcM6/Bp6N83tzhKMXlnhUYEUP/oewXlXRyqbqheVxhf9uZ0fZ9h+4elWF/26x/i9vuL6/AfyyLXY5DYaiKMoBh1gDWH9kjFmVaVsOJFQ4FEVRlKTQ5LiiKIqSFCociqIoSlKocCiKoihJEXWWzQOJvn37mqFDh2baDEVRlC7FokWLKo0xxeHLDwrhGDp0KAsXLsy0GYqiKF0KEdkcabmGqhRFUZSkUOFQFEVRkkKFQ1EURUkKFQ5FURQlKVQ4FEVRlKRQ4VAURVGSQoVDURRFSQoVjhi8t2onf5tdmmkzFEVROhUqHDGYvaaCRz/amGkzFEVROhUqHDHwCAR02nlFUZQQVDhiICIEAiociqIoblQ4YuARQR0ORVGUUFQ4YqChKkVRlNaocMTA4xE0UqUoihKKCkcMRD0ORVGUVqRUOETkXBFZIyKlInJLhPUiIg/Y65eKyIR4ZUXkRRFZYr82iciSVNmvOQ5FUZTWpOxBTiLiBR4EzgbKgAUiMtMYs9K12RRgpP06FngIODZWWWPMZa5j/BGoStU5aI5DURSlNan0OCYDpcaYDcaYJuAFYGrYNlOBp4zFp0BPERmQSFkREeBS4PlUnYBHBL8Kh6IoSgipFI4SYKvre5m9LJFtEil7MrDTGLMu0sFF5DoRWSgiCysqKtpgfkuoyqh4KIqiBEmlcEiEZeEtcLRtEil7BTG8DWPMI8aYScaYScXFrZ61nhAeEXtfbSquKIpyQJKyHAeWlzDI9X0gsD3BbXJilRWRLOAiYGIH2tsKjy1fAWPwRNQyRVGUg49UehwLgJEiMkxEcoDLgZlh28wErrF7Vx0HVBljyhMoexaw2hhTlkL78djKoWM5FEVRWkiZx2GM8YnIjcBbgBd43BizQkSm2+tnALOA84BSoA6YFqusa/eXk8KkuIO4PA5FURTFIpWhKowxs7DEwb1shuuzAW5ItKxr3bUdZ2V0NMehKIrSGh05HgOPehyKoiitUOGIgeNxqHAoiqK0oMIRAxFNjiuKooSjwhEDJ1SlAwAVRVFaUOGIgROq8qvLoSiKEkSFIwY6jkNRFKU1Khwx0FCVoihKa1Q4YuDR5LiiKEorVDhioOM4FEVRWqPCEQPRcRyKoiitUOGIgU45oiiK0hoVjhhoqEpRFKU1Khwx0OS4oihKa1Q4YqDTqiuKorRGhSMGLTkOFQ5FURQHFY4YeD3OlCMZNkRRFKUTkVLhEJFzRWSNiJSKyC0R1ouIPGCvXyoiExIpKyI/sNetEJF7U2W/JscVRVFak7InAIqIF3gQOBsoAxaIyExjzErXZlOAkfbrWOAh4NhYZUXkdGAqMM4Y0ygi/VJ4DoAKh6IoiptUehyTgVJjzAZjTBPwAlaD72Yq8JSx+BToKSID4pS9HrjHGNMIYIzZlaoT0HEciqIorUmlcJQAW13fy+xliWwTq+wo4GQRmS8iH4rIMR1qtQsNVSmKorQmZaEqQCIsC2+Bo20Tq2wW0As4DjgGeElEDjNhXZ9E5DrgOoDBgwcnYXYLOo5DURSlNan0OMqAQa7vA4HtCW4Tq2wZ8G87vPUZEAD6hh/cGPOIMWaSMWZScXFxm05Ax3EoiqK0JpXCsQAYKSLDRCQHuByYGbbNTOAau3fVcUCVMaY8TtlXgTMARGQUkANUpuIEdByHoihKa1IWqjLG+ETkRuAtwAs8boxZISLT7fUzgFnAeUApUAdMi1XW3vXjwOMishxoAr4ZHqbqKDRUpSiK0ppU5jgwxszCEgf3shmuzwa4IdGy9vIm4KqOtTQyweS4KoeiKEoQHTkeA+eZ434NVSmKogRR4YiBjuNQFEVpjQpHDHQch6IoSmtUOGIgmhxXFEVphQpHDNTjUBRFaY0KRwx0HIeiKEprVDhiEBzHoc/jUBRFCaLCEQOdckRRFKU1Khwx0JHjiqIorVHhiIHHrh3NcSiKorSgwhEDr3ociqIorVDhiIEzjkOnHFEURWlBhSMGzjgODVUpiqK0oMIRg5bkuAqHoiiKgwpHDHQch6IoSmtUOGKg4zgURVFao8IRA+d5HKobiqIoLaRUOETkXBFZIyKlInJLhPUiIg/Y65eKyIR4ZUXkdhHZJiJL7Nd5qbJfJzlUFEVpTcqEQ0S8wIPAFGAMcIWIjAnbbAow0n5dBzyUYNk/GWPG269Wj5ftKHTkuKIoSmtS6XFMBkqNMRvs54S/AEwN22Yq8JSx+BToKSIDEiybcjTHoSiK0ppUCkcJsNX1vcxelsg28creaIe2HheRXpEOLiLXichCEVlYUVHRphPQadUVRVFak0rhkAjLwlvgaNvEKvsQMBwYD5QDf4x0cGPMI8aYScaYScXFxQkZHI5OOaIoitKarBTuuwwY5Po+ENie4DY50coaY3Y6C0Xk78BrHWdyKI7H4VflUBRFCZJKj2MBMFJEholIDnA5MDNsm5nANXbvquOAKmNMeayydg7E4WvA8lSdgNi1ozkORVGUFlLmcRhjfCJyI/AW4AUeN8asEJHp9voZwCzgPKAUqAOmxSpr7/peERmPFbraBHwvVefQkuNI1REURVG6HqkMVWF3lZ0VtmyG67MBbki0rL386g42Myo6jkNRFKU1OnI8BjqOQ1EUpTUqHDHIxDiOFxds4f531qbteIqiKMmiwhGDTIzj+MXLy3jgvXVpO56iKEqyqHDEQENViqIorVHhiIEmxxVFUVqjwhEDEUFEPQ5FURQ3Khxx8IgQUOVQFEUJosIRB49oqEpRFMWNCkccRCRtoapPSiuDn3VGXkVROisqHHHwSHoa8U9KK7ny0fnB7z4NjymK0klR4YiDRyQtoaqq+uaQ702+QMqPqSiK0hZUOOLgSVOoqjAvdNqwZr8Kh6IonRMVjjhImpLj3bK9Id/V41AUpbOiwhEHj0haplUP92qa1ONQFKWTosIRh3R1xw1/ymCzX5PjiqJ0TlQ44uD1pCc5Hn4MDVUpitJZSalwiMi5IrJGREpF5JYI60VEHrDXLxWRCUmU/ZmIGBHpm+JzSEtyPFw4NDmuKEpnJWXCISJe4EFgCjAGuEJExoRtNgUYab+uAx5KpKyIDALOBrakyn4Hj5CWKUfCQ1WN6nEoitJJSaXHMRkoNcZsMMY0AS8AU8O2mQo8ZSw+BXqKyIAEyv4JuBnrueMpJV3jONTjUBSlq5BK4SgBtrq+l9nLEtkmalkRuRDYZoz5ItbBReQ6EVkoIgsrKiradgakbxxHuE6ocCiK0llJpXBIhGXhTXC0bSIuF5F84Dbgf+Id3BjziDFmkjFmUnFxcVxjo5GucRzhoSpNjiuK0llJpXCUAYNc3wcC2xPcJtry4cAw4AsR2WQvXywih3So5S7SN45DQ1WKonQNUikcC4CRIjJMRHKAy4GZYdvMBK6xe1cdB1QZY8qjlTXGLDPG9DPGDDXGDMUSmAnGmB2pOol0j+OYcdVEQJPjiqJ0XrLibwIiciKwxBhTKyJXAROAvxhjNkcrY4zxiciNwFuAF3jcGLNCRKbb62cAs4DzgFKgDpgWq2xbT7I9pCvH4YhTXral5ToAUFGUzkpCwoHVTfZoETkaqzfTY8BTwKmxChljZmGJg3vZDNdnA9yQaNkI2wxNwPZ2ka4ch3OM3CxrzirNcSiK0llJNFTlsxv5qViexl+A7qkzq/Ng5TjSEaqy3rvlWMKhOQ5FUToriXoc+0XkVuAq4BR7gF526szqPHg9QiANbbgzyLAlVKXCoShK5yRRj+MyoBH4tp2ILgHuS5lVnQhJ0wBAv30MZ3p1TY4ritJZScjjsMXiftf3LVg5jgOedPeqysvWUJWiKJ2bhDwOEblIRNaJSJWIVIvIfhGpTrVxnYF096ryeoQsj2hyXFGUTkuiOY57ga8YY1al0pjOSLo9Dq8I2V6PehyKonRaEs1x7DwYRQPSOa269e7xCDlZHh3HoShKpyVRj2OhiLwIvIqVJAfAGPPvVBjVmfAIaemO6/Sq8ghke0UfHasoSqclUeEowhrZ/WXXMgMcBMKR3l5VXo9Yx0yHm6MoitIGEu1VNS3VhnRWrEY89cfxBz0OSZtYKYqitIVEe1WNEpH3RGS5/X2ciPwqtaZ1DtI25UigxeOwnnOe8kMqiqK0iUST438HbgWaAYwxS7FmrD3g8XrSM616MFQlYomVKoeiKJ2URIUj3xjzWdgyX0cb0xlJ26NjnVBV0ONQ4VAUpXOSqHBUishw7Cf4icjXgfKUWdWJEGnxBlJJwFjeDaRv0KGiKEpbSLRX1Q3AI8BoEdkGbAS+kTKrOhFpe+a4Mdi6kTaxUhRFaQsxhUNEfuL6Ogv4AMtLqQUuxjV/1YFKOsdxeMRSDm+apnJXFEVpC/E8DueZG4cDxwD/AQS4GpiTQrs6DWkbxxEwoaEqHf+nKEonJWaOwxhzhzHmDqAv1rO9f2aM+SkwERgYb+cicq6IrBGRUhG5JcJ6EZEH7PVLRWRCvLIicqe97RIReVtEDk3mhJNF0jWOwxi8tsehoSpFUToziSbHBwNNru9NwNBYBeyHPT0ITAHGAFeIyJiwzaYAI+3XdViPqI1X9j5jzDhjzHjgNeB/EjyHNtHeSQ5Pu+8DfvD853G3CwQMHtvjsLoAq3BEoskXYPGWvZk2Q1EOahIVjqeBz0TkdhH5DTAfeDJOmclAqTFmgzGmCXgB69GzbqYCTxmLT4GeIjIgVlljjHs69wLsnl6pwnp0bNvLb9pdx3+/2B53O78JDVX5tVtVRO57azUX/e0T1u7cn2lTFOWgJSHhMMb8LzAN2AvsA6YZY+6OU6wE2Or6XmYvS2SbmGVF5H9FZCtWz66IHoeIXCciC0VkYUVFRRxTo+PxpGtadYLJcU8KR46v27mfRp8/NTtPA2t21gCwbW99hi1RlIOXRD0OjDGLjTF/sV/xYy9WEr3VbhLcJmZZY8xtxphBwLPAjVHsfcQYM8kYM6m4uDgBcyOTruS4MQavxzlmasRqydZ9nP2nOTw2d2OH7ztdFORYT0isbTooxp8qSqckYeFoA2XAINf3gUB4zCbaNomUBXgOq1twymhvqCpR/K7uuKkSK0cwahq6bqNbkGt1BKxr6rpek6J0dVIpHAuAkSIyTERysOa2mhm2zUzgGrt31XFAlTGmPFZZERnpKn8hsDqF55C+JwCa0HEcqejJtWSrlVTOt+/auyKO7V1Z/BSlq5PoyPGkMcb4RORG4C3ACzxujFkhItPt9TOwBhWeB5RiPe9jWqyy9q7vEZHDgQCwGZieqnMAO1GdpgGATnI81TPyNnbh55l3y7aEY19dU5wtFUVJFSkTDgBjzCwscXAvm+H6bLCmM0morL08paGpcNI3jiN0ripfCg7qsx9H29SFhcN5pO7euuYMW6IoBy+pDFUdEKR3yhHrc6qex+Gzd9qVPQ6nR9he9TgUJWOocMShoyY5jCc+/rBQVSrGcfjs55h3ZeFwvKV96nEoSsZQ4YhDR43jiNdYB9zJ8RSNHD8QQlVOPe6pVY9DUTKFCkccpIM8jsbmxIUjVVO5t4Squm5XVkf0qurV41CUTKHCEYeOynHUN8durENnx01RqMpOuHdtj8Oqx1R0HlAUJTFUOOLgbcdgPPdzw+MKhyE4yWGqBgAeCMnxJjtPo3N5KUrmUOGIQ3tCVe7xH/VxRjoHAgav3asqFaPV/QET3GdX9jgc251uuYqipB8Vjji05+7ffVfcECevEBKq8nT88zia/S1i0ZVzHI63pB6HomQOFY44WDmOtpUNEY44Hoe/VXK8YxtGty1N/q7vcWiOQ1EyhwpHHDyetj8bIyRUFSfHEWj16NiOFQ6fK7QTr4dXZ0Y9DkXJPCoccWjPvFF+fxLCYUJ7VXV0u9jsukM/MDwOFQ5FyRQqHHFoT6I6meS431iJeHAe5JTCUFUXTo47+Rlj1OtQlEyhwhGH9kyr7g43NSQSqnL1quroUJWTHM/N8nTp7rhu2zXPoSiZQYUjDu1JVLvDKfFCK+5eVd4UjBx3chyFuVld2uNo8gWwHTP1OBQlQ6hwxKFd4zhcBeM1ciFTjqTgOeeOcOXnertsd1xjDI2+APn2Mzk0z6EomUGFIw7OVOdtmXbE3fjHE47Q2XE7PsfhhHUKcrJo9psOD4WlA2fQX779+FifDgJUlIyQUuEQkXNFZI2IlIrILRHWi4g8YK9fKiIT4pUVkftEZLW9/Ssi0jOV5+C1vYC2tLPuO+J4A/r8xgSnHEllqMp59GpX7FnleEoFOY7H0fXOQVEOBFImHCLiBR4EpgBjgCtEZEzYZlOAkfbrOuChBMq+AxxpjBkHrAVuTdU5QMv8UW3xANx39f44d8dWcjx1kxw6IlZg3613xQS5k5vplmOdg+Y4FCUzpNLjmAyUGmM2GGOagBeAqWHbTAWeMhafAj1FZECsssaYt40xPrv8p8DAFJ5DMBHbFuFwexnxPI6AaQmLpSRUZXsYjsfRFfMcjtgFPQ4NVSlKRkilcJQAW13fy+xliWyTSFmAbwFvRDq4iFwnIgtFZGFFRUWSprfgJKzbEhVJJjnuD7hCVZ6On+TQyQ84HkdX7FnV4nFoclxRMkkqhUMiLAv/p0fbJm5ZEbkN8AHPRjq4MeYRY8wkY8yk4uLiBMyNjKc9HkeSvapSGapy9leQ03WFwxmL0s3uVeXXHIeiZISsFO67DBjk+j4Q2J7gNjmxyorIN4ELgDNNKp6x6iLocaRYOEJnx+34UJUz5UhLqKrrNbrNYQl+9TgUJTOk0uNYAIwUkWEikgNcDswM22YmcI3du+o4oMoYUx6rrIicC/wCuNAYU5dC+7GPB7StV1Uy3XEDrl5VqXgeh5MPcMI8zV2wV5Vjc1625jgUJZOkzOMwxvhE5EbgLcALPG6MWSEi0+31M4BZwHlAKVAHTItV1t71X4Fc4B27Uf/UGDM9VefRnnEc7oYtoZHj7lBVh89VFepxHAjCob2qFCUzpDJUhTFmFpY4uJfNcH02wA2JlrWXj+hgM2PiaYfH4W7844WeWk850tFzVdkeR7Y35HtXojnMa9JxHIqSGXTkeBzakxx3t2vxPA5jWkRK7FBVR6ZvnLtzZwxEV/Y4uoWFqnz+AFc9Op9PN+zOmG2KcjChwhGH9gwAdN8Rx5viw3oCoH3Mdng50QhvdLuicDj12S0sVLW7tom5pZV8tnFPxmxTlIMJFY44OI14W27+3WITz+PwBQxerzOOo3X59hKc5NCZcsTX9UJVjs152VYFOedUXd8MwP6G5swYpigHGSoccWjfOI6Wz/E8jkDAkOWa5NAq3/HCcUD1qrI9kGpbMKrrfZELKorSoahwxKE93XHdA9RieRzGGMvj8FiXw0mSd2R+3HcgharCphxxBGN/o3ocipIOVDji0DLlSPs8jljda51du7vjWss70OM4EMZx+EJ7hjkemXocipJeVDji0K5QlV3G65GYs+M6DWCWt2UAoLt8RxAMVWU706p3wRxHeK+qsBxHteY4FCUtqHDEoT09nBwvJcfriSkCjnB4PaHCYTrQKQiGqhyPowtOOeKcQ15OeI7DDlU1qMehKOlAhSMO7ZlW3bkjzsnyxEx0Ow1geKiqIz2O5jCPo0uGqsIGMbbkOJpD3hVFSS0qHHFo6Y7b9gc5xRMOJ4funuQQOjbH4Q8EyPIIOVnWJe+SwhFlHIcTotrf4OvQQZOKokRGhSMOHTHlSI43MY8jPMfR0clxr0eCXX67Yo6jOTiOIzzHYYWomvyBLjnrr6J0NVQ44tCewXiJhqqcdR4JE44ObAOb/YZsrwcRIcfr6Zoehz+ARwh6TeEeB2i4SlHSgQpHHKQdjXhIcjyWcNiilOVJ3chxfyAQ9GiyvRJMNHclmgMBsr2eYEjPET+3WFRrglxRUo4KRxw64kFOOVmxe1U5SV5P2MjxjhSOZtfI9Cyvp2vOjuuzvCbnPJz6rW3yB6ch0WlHFCX1qHDEoeV5HMmXdQtHrJHjwXEcntSFqnz+AFn2yPRsryc4JqIr0ewPkO2VoOfk1Gmjz09RXjbQNR+JqyhdDRWOOLTL47DLZHsl5vOx3QMFrXfafMxo+Pwm2ODmeKVrjuOwQ1WOADqC29gcoKibLRxdUBAVpauhwhEHaceYihaPw0us9izaAMCOHMfR5A+QYytSdlbXTI432aEqxwt08jQNzX6652XZ23S981KUrkZKhUNEzhWRNSJSKiK3RFgvIvKAvX6piEyIV1ZELhGRFSISEJFJqbQf2jeOwx+SHI/hcYSFqqQdx4yG2+PI7qI5DsvjEESsbsUtoapAMFSl3XEVJfWkTDhExAs8CEwBxgBXiMiYsM2mACPt13XAQwmUXQ5cBMxJle1u2jWOI+hxSELdcYOz46bgQU5OmAe6do4jyz6HLK9Vp8YYSzi6df4cx/SnF/H60vJMm6Eo7SaVHsdkoNQYs8EY0wS8AEwN22Yq8JSx+BToKSIDYpU1xqwyxqxJod0hBCc5bEMrHrCf6pfliTcA0BGO0GN25PM4mvwm2OjmeKVLh6rAqlNfwAQFsLOHqgIBw5srdnDDc4szbYqitJtUCkcJsNX1vcxelsg2iZSNiYhcJyILRWRhRUVFMkXD9wO07e7fesaGWLPjJjTJoXU5UjHliM8fINvjDlV1zgY2Fr5AgJzgUxKtsSgNzdZ5tISq/BmzLxZ1zZ3TLkVpC6kUDomwLLwljLZNImVjYox5xBgzyRgzqbi4OJmiIbR0x23bXFUekcSnVQ+fHbdDH+RkQkJVzV3w0bEhoSo7x+EIRVE3y+PorDmOusauMzAxEDCs3F6daTOUTkwqhaMMGOT6PhDYnuA2iZRNC15P+3IcXo/gldgehzNXlSd8dtwODVW1jBzP8krXzHH4DNmuc/AHDI3NTqiqc3fHrelCwvHn99Zx3gMfsW7n/kybonRSUikcC4CRIjJMRHKAy4GZYdvMBK6xe1cdB1QZY8oTLJsW2jOK229s4fDGTo47Ha6CkxymIlTlSo532bmqXOfg5DgcD6PIznE4QtLZqG1sCVXd/K8vuPqx+Rm0JjYfrN4FwD6d90uJQlaqdmyM8YnIjcBbgBd43BizQkSm2+tnALOA84BSoA6YFqssgIh8Dfg/oBh4XUSWGGPOSdV5tOsJgG6PI4HZcVtNctjB3XGzvV07x2GNHG95LrvPHwiGqvKyvWR3Yk/K7XG8tLAsg5bEZ29dEwB1TZqXUSKTMuEAMMbMwhIH97IZrs8GuCHRsvbyV4BXOtbS6LQn3+APGLx2jiO5KUes5R3ZHbfJlR/IzvIE58fqSrjFLzfLQ6OvJTmem+Uhx+vptL2q6pq6TqhqX50+GEuJjY4cj0NwFHdbu+PavapidecNHzkeHMfRgcrh8xtXr6rOe2ceC7f45WV7aWj2h3gcudneTturKlKOo7M+dMqxtUqFQ4mCCkcc2vPo2EafNc1HVoIehzds5HhHTjni8x8AOQ7XtCl52R4amlse3NTZPQ53jsOhM/YAc88uXB0207Ax2ttKsVDhiEN7Ro43+gLkZnvweCSm8IQ/j6M9M/JGwz0AsMtOORISqrK8i8ZgqMpLbnbnFY5IoarwhrkzsLfW/VCsUJtfXbKN8x74iHdW7oxavq7Jx43PLWbmFxnpBKmkCRWOONhj8toUVmhqo8fhTVmvqpbkeGMXHJDW6Gvxmlo8Dn/we47X0+ou/rG5G1mxvSrttobjhH+OGdoruGx/J3zoVL3rd7GzuiEk9Ldim+VtrNsVvZvuJ6W7eW1pOT98/nN2VTekzlAlo6hwxKE9HkeTL0ButhePCMZEz1k4iepWoaqOznHYjW5Brpe6Zn+njbFHo6bRR2Gu1Z8jN9tLg9vjyPaSkxXqcdQ2+rjztZVcMmNeRux1U9voIy/bw+PXHsO0E4cCnV84Xvl8W0jdOb8WiTg+1+KT9buDn3ftb+xw+5TOgQpHHNrTHbfR5yfX/cS6KPto/TyOjh857h4AWJCbhTGhjURnp9kfoMkXoMAWjrwsL40ujyM3y0NuVujkjZt21wKdo1tpbZOfwtwsuudlc95RA4DO+bTC+rC6WlpWFVzm3MjE6iE2b8Pu4O9Xe2UduKS0O+6BgOeaa2D8dwj89rewe3VSZZvGXkm+vxnPnC0w+BT8Z51NtmndiPn7jYPDziHr0kuhuQZPQX846hr8t9wK+9Z3yHn4jv0pOc8+A7//mIL+42HY2dSefyH5zXUdsv9UU+vNg2N+QMHfZ8Cdi8gddjaNvUbS+JdZMPR0cqd+hZxRX6NRBJ7+BQAbex8Ooy5EjIHTTsus/SPOJ7/wUDjtNArzi2Hctez/xW2wZ21G7QqnoecwGP31kGUfXfY9vry3lN0jLoC+R1D5/Mvwu7dbld2blceqST/g/N2reb3PaKp/fivsXZcu05VozJ7d4btUjyMOHvu238Rwz6PR6Mki1/jIMtZdsF8iV7fPXu6xgwHOMQOS/DEj4UcIiCdoR77fuhOs8+R0yP7TQY3XsrXQbw1Oywv4aPBk0eDJCn7PMX4aPd5gmY3drHxCbiDzIaFaTw4Ftu3dfVYIZ783N5MmRaTBk91q2Ya83gBUZBfY7/kRy35aNBiAc/ZYYlGVlZcKE5VOgHoccfA89xzc9wGB226DCQOTKtt0/4fk9CvEO6QXvL4K/6xZkNf6jxn4ZBPMXEHWq69AQQ6e7dXwwEcE7rwLjjyk3efQ3OyHX79J1ne+DaffTcHyHfDMImqeeBIO7dHu/aeD2h374c9zKPj1rTDuUPLeXE3DnA00XjsNzwelZL3/HrlPL6Jpbz38dTYAG19cAp9vo8GbTcM775GX7Y15jFRS/fA8uhvg/2bTvb4Z7nib/T/8MZx8WMZsikT94jJ46QvGHlrEiu3VZHuFqm9+C869l4o/zoaKWiqPPgYeualV2U//s5z8RWWc/NQDcOc7VP3gJjhleMLHvv/tNayvqOX/rvhScNodpXOiHkccWsZxJF+20RcgN8sTjPlGmyHXF9aryunJ1VG9qpz957iS49A5Yv+J4vRKCuY4sr34AobaRj+5WV5EhJwsT0gvoLJ99cHP5VWJ9fCpqmtOyYSEe2qb6FNoe032OVR34uT436+ZxMrfnkOPbtnBgYAVdrK7siZy0nvT7jpG9CukZ342HmndnTcWe2ubeOD9Ul5fVs6/FnXuKVm6Chsra7nwr3NZsGlPh+9bhSMO7ZlwsMkXIMctHNGS4/ZcVeHPHO8w4bATxk5yPD/Harhq0zxja2VNY5un3nBsDfaqyrJ+ulX1zeRmW59zwwYAVtU1Bx/wNLe0Mm633L21TRz927e5/plFbbIxFrtrGuldYAmH1yMU5mYlnDzevq+eN5al58mBTiK8IDeL/JwsivKyqa5vpr7JHxS6yprGiD3yKmsa6VuYi4hQ5BKcRHjdPr/cLA8Pzi6N2gNx/obd/OGttD3HrU0EAobdUcQ1nazZsZ+lZVXBG8aORIUjDu15AmCjzx8qHGH7+MfHGyndVYPTESj8eRwd1Ru3KSgcoR5HpNHMHUUgYELuTI0xTLrrXa5+7LM27c8RjoKcFo8DYHdto6uLbqhw7Ktv4sThfcnyCL9+dTnnPzA3KKKReHjOBgA+WlfZJhuj4fMH2FffTJ/ClpzGwF7d2LonsY4Jlz48j+ufXUxDnF5wHTHdinOMbnb9OgKwsdLqoTZ+UE8amgMRvaXdNU30scWxR5LCsaOqAa9HuO+So9m8u44bnlvMzf/6opVAfesfC/jrB6Vsd3mTYNVx+LiRrXvqkuq5ZozhK/83lzP+OJstu9veaeT2/65g4l3vsqe2qc376AhK7fE2I/oVdvi+VTji0O5xHFne4NxTbuH4Yus+bv/vSm57ZVnQ43CO5QzUi9dQJIozTsR5ep7T+NamcOK9u15fxSTXn2ftzhoAFm3eG7OcP2CY/vQi3l8dOjq5JszjyLO9jMqaFuFwDwA0xrC3rpkhffKDXgdY3UXD2bavnr21Taze0TKdRngj9Maycl5auDW8aELsrWvGGIKNKsDwfoWs21WTUPmyvfVBO6PxzsqdjL/jHXa2c9BdQ3MAr0eCv8Ee3SyPY0OlZesJw/sAtDqOMYbdtY307Z7bUi6JRntPXRO98rP58pj+ZHmEN5bv4KWFZa1GoHezf7sPf7g+JKT44Afrmfy79zjjj7P53tMLmbWsnJPv/YDbXlmesA0VNY0s21bFhopaZq/dFbIuEDA8/OH6uIMad+1v4Kl5mwF4YcEWzn/go1ahos827qE0wWvfHtbtqqGkZ7dgeLcjUeGIQ3vmqmryB6J6HM98av24CnKzgjkIx+MY0KMb2V4J3uW1F0c4sjyOx2H9kFL1VLr6Jj+Pf7wRgNlrrD/gR+taHt8bfrfo5tMNu3lzxQ5++tIXwe87qxtaQlV5TqjK9jhqmoLC4ExyaIyhodka99EzP4dbpowO7v+N5TtCjhcIGE68532mPvgxGytrKbYbvsm/e49V5S1Ccv2zi7n5X0uDd8Afl1ayeIslgou37OU7Ty4I2d6NI55OjgNgZL9Ctu6to6HZz5Kt+2KKgsO2vdG3eW/VTuqb/Xy4JrHHJAcChjeWlbeas6y+2U+3bG9wEKrjcWyosH6Lx9vCEZ4zqqpvptlvguJYlJecx7Gnpole+TnkZXsZ0qel19ZbK3YQCJignY63/OS8zZz+h9ns2m/ZMbfUOu8NFbW8tWIn33/Werb7Oyt3JjzQtXRnS2O+fFtoWHPF9mrufmM11zwe22Oev6FFJP7y7jpWbK/mkhnzgp7u/oZmLn14HlP/Opf/LNnG9c8sCvGCjTHMXVfZIROcrttZw8j+He9tgApHXFqmVU/uQlo/dhOaHLd/DP6A4W17vp/aRh+BgEGkJZ+Sk+VhWN8C1u6I/wS2Zn+Asr2WWx2tQW5qleOwQ1VxkuPWM72tbZ6et4mfvLgkuO4v767jL++uY2NlLY/MWc8al60frm1pvN5bZQnHm64GO1ay7uXFVmI0N8tLTaOPyx/5lGlPLAja6jQcjsexu6YpKISHFOXR0BygoqYx+EyJnvnZXHbMYDbefR6njipmcZjHs6RsHwBb9tSxZU8dF4wbEFz3xMcbmb9hN5tcAv6fJduZ+te5fOPR+Vz0t094ZM56pj2xgHdX7eKNZeUEAiYYHtlb28QHq3cFByL2dnkcI/oVYozVQH31wY+5dMY8Pli9i/vfaRnXcfesVdz/dks8P9L1bfT5qWvy8dlGq04/WLOLu99YxebdtXxSWsn8CB4WwNsrd3D9s4uZtaychmY/P//nFzz4QSl1Tf5g3QL06JZFdYOPDRXW3evQPlaX3J1hwlFZY9V3sdvjcAnHn95Zy9PzNrWyo3RXDat3VLOnrikkB+Swo6qB+95ew7l/nkMgYCjf18CEwT2ZfupwKvY3BhtqJxF//rgBITH9+mY/6yti34Bt3VNHsz9AaYUlHIf3787ybaE3Acvt/NjqHfujdg4AWFlu9UT71onDQqa/cTzuVz7fBlj/vV+9spw3lu/guc+20OQLULG/kVeXbOOqx+bz3GdbACtP926EucE2767l+mcWBa+7m9pGH6vKq1lfUcPIFISpQLvjxiVSqGpTZS3vrtrJJRMH0SPf6l5rjAnepUFLY+32OBzPYv7G3VTVNyNiTcvgs5/b4WZU/+58YTdqGytr+ePba7j36+OCiW2H659ZxLurdvHidcdx2SOfcu/Xx3HppEEh2zgPinL+ULlZ1mj2eMnxH724hNmrd7H09nP49X9WAPDzcw/HHzD86V2rgXPeX1ywlXd/cioiwkfrKijI8fLVL5Xw4oKtvLGsnIWb93LzuYfzp3fWsrK8mqnjS1odb0dVA699YSVJd1Q3BMVmQ2UN+xt8ZHsl6Gnk2jH4Jn8gGKoafUh3ANbuqAk2Qr3s6yMiTBzSiz+9u5aq+mZ6dLOWv7WiRdCMgaMH9uT9n57K32av56WFZa0eunSTSzwH987nd7NaBoWu21XDtf9YwJy1Fdxz0VH8c1FZSGiuT0FLjmOs3Q366/aUHtv21TPtHwsAq9G97JhBwZyLg+OV7G9o5scvLmHq+BKe/2xLcJqPbK8EPaq56ypZYc9k+/drJnHSiL50y2npjuz0XFqwaQ9Fedn80/5+ZElRSLdlJ1exesd+DisuoF9RbvD6gNXTqskfCMbT+9p5nJJe3Xhn1U7qmnx4RPjLe9bYjm8cOySkq+1Z938IWB6YE4u/48IjuX3mCob0yWfF9mreWr6DDZW1zNuwmyZ/gK9NGMglEwfy8Jz1lO6qodHnZ31FDTecPpyfnzOas+//kHW7arhk4kD+/fk2pj+ziJevPyF4zd3894vt/OD5z7nprJGW95qbxemj+/HoRxuorGmke14WizbtZZnLA/n9G6upa/Izf+Nupp04jG7ZXj5Ys4t7vz6OFdurGdGvO1/90qE8/vFGBvTIo7yqgRXbqxhzaBFz1rbkz/bb/7+/f7SBv75fyv4GH6ePLgbgV68uZ2ifAp7/bAuvLyvnzZtOpnJ/ExU1DXx1fAmvfr6dN5bv4I3lO/jp2aPYtq+eG04fQVG3bM66/8NgDzjHQ+xoUupxiMi5IrJGREpF5JYI60VEHrDXLxWRCfHKikhvEXlHRNbZ773C99uReMN6OBljuOnFJdz1+iounvEJdU0+bp+5gqNuf5tXPm9pZNyztjrCUV5Vz8wvtnPl3+eT5REunjCQndUNwScFuhnVvztb99RT2+jjphc+57Wl5cy1k7Z1TT42VdbS0OznXfuO3mkwnp1v3alsrKwNeknBUJUtHCJCfo43ZndcYwyvLy2ntsnPE3bYCeDuWas56fcfAHDsMGtg2OH9u7O+opYlW/cRCBg+XFvB8cP78KMzR5Lt9XD9s4spyPFy5eTBjOjXndXloZ7UjqoGKmsaeXLeJvzGcO/F4wD42+xSAPoX5VHb6AuJ1eZltTRuTqhqlC0c//hkY9AL69Gt5S5/4pBeGAOLNu+hrsnH43M38urn2+hf1NKgj+xfyGHFhfz47FFR6wasyQrf+NHJwe8nj+zL3HWVzLG9rUfnbgyGshzcoaphfQu4+6Kjgr3D3Nz52kpun7mi1fJte+tZvaOaX7+6nHdX7eK3r60MmRvq+lNbxkyscE1//t2nFjL2N2/ym/8sZ1NlLZU1jcy2Q1oLN+1l4eaWu9bl26qDiXGwQk7+gGH1jv1MHtqb3Cwv3bK93P/OWh6bu5Fj/vddTrznfaY/szjkHE8bVUyTL8AlM+aF5CkO++Us3lxezr8Xl/GqffcNluj2ssX++OF9eOvHpzC8XyHbq+rZYHt8LyywckwDe3UjL9vLoF75rK+oYd3OGnwBw+GHFFnHKLa8opNHFfPoNZMo3VXDA++to6HZzwdrdnH/22swxrB8WxU//acVEn1v1S4+Lq1k9IDuXDShBI9H+NUry7n3zTVc+eh8npu/hQmDewLwz0VlvL6snMqaJu57aw2/fW0lH62r5JE5G1i5vYqxhxZxVEkPzhjdjx+fNYqcLA8//9dSnvxkE2t37ue4w3oHz3ton3y27qln1/5G6pv9zFrWciNz3dMLWbJ1HwC/emU5Vz02nx+/+AXzNuxm0Za99C7IoWd+Nn98Zy0vLy7jsofn8etXlwdFA+D4w/q2+h11BCnzOETECzwInA2UAQtEZKYxZqVrsynASPt1LPAQcGycsrcA7xlj7rEF5RbgF6k6j9xs6+58ve3GfrBmF0u27uP8cQOYtaycqx/7jI2VtdQ0+rj138sIBODiiQNp9FuNck6Wh+72oD+nR9EhRXncf+nRLN9exb8W+amqbw7mNxxG9bcawZXl1cFGYPm2Kk4Y0ZcT7n6P6gYfD14Z1FneXWW5s2t37GfR5r1c/NAn/OYrYzj98H7BGLuT8AQrz+EkF53eOB+s3sUTH2/iue8eF9Lj567XVwU/O43A/37tSL4+cSCfbdzDuIE9OeHu93ho9nrOGtOfsr31/Pycw+lnn+dNLy7h9gvH0jM/hyMO6R7S2DU0+znu7vcY0icfn99wysi+nHlEP4BgXH3z7jqe3r05JNST6wqnOB6Hc7f77qpdzLOP0aug5S5z4pBe9MzP5lv/WIjX0/I435+fczh9C3PolZ/DmAFW41PSsxvPf/e44DX87xfbOeuI/jw8Zz13XDiWYX0LEBFm/fBkdlTXs2jz3mBvrAvGDeC1pZbn9PL1J/DAe+tYtq2KXvmhI/WvmDyYr32phLnrKvnOUwtD1oWPZTh2WG9mLS/n367G1t1A9C3M4arjh/DA+6XBZXnZHn479Uhu/tdSBvToxpPzNvPkvM30656LL2A4/6gBzFpeTrbXwxEDiijdtZ9mvwnxTNx36SeMsO5enbEed77m/itbHFJkjRY/xr6pWLG9mpv/tRRoeWqjIzLh9A6rn0OK8oLztXXPzeK/9m9vRLHlmQwvLqB0Vw3vrdqFCBxnH/Ow4kJgJ4f1LeDIkh5cMG4Aj83dyNPzNgcjAX0Kc3n84430zs9h8rDewd/1zecezqj+3Zl2wlAenrMh2KuyuHsu3z9tRPA6/fDMkVwycSBn/HE2uVlejh/ehyc+3gTAiSP6ICI8fu0x1j5fts7/N/bNwCUTB3LJxEH84e01/PjsUfzohSUh5/3UtyYzZ20Fj87dSF2T5WUutD1Xj8Bf3y/li637mPqlEqafMhxfwApzXfbIp8z8YjvTTx1O2d46crI8IdeyI0llqGoyUGqM2QAgIi8AUwH3r20q8JT9CNlPRaSniAwAhsYoOxU4zS7/JDCbFApHXraXiyaU8M+FZfzozFH8c2EZfQpy+PNl4znriH78+EXrjuWKyYOYs7aSn/7zC5Zts+46wBpb4NyZgxVzf+zaSYw9tAcVdqx02776ViNlR9lJrdeXlgdDXJ9v3ceM2euDXSEfm7shOCOs0/OmvtnPX9+3wgJ3/Hcld/x3ZTCEk+2K/RbkZgVj8Rc8MBevR1ht5ylWlVfzvaetsQyv3nAin6yvpLreh0fgb7PXc2RJEd84dggAJ4+0XOvvnz6C+95aw6cbdjN+UE8uPPpQAKYcNYAzjugXDDGNHtCdf3++jfP+8hFPfXsyP7GT4Jvt7o8/O2cUfQpzKcqzYutOYwOhbrfb4yjMbWncrjx2MM/N3xLMifR0eRx52V4unTSIR+Zs4LRRxSzbVsWu/Y2cNKIvRw/qGX7pQ443cYjl2J40MvQObsyhRYw5tIgmn3WNJg/rza/OH8Om3bX4/Ibxg3ry5Lcm0+QLtPIqHZvcx/nzZeNp8geCjS1AQY6Xc8YewnxXPNsRp8P6FnD/ZeM5pCiPft3zOHtMf44YUMQxQ3sxoEc3RvQr5NRRxRQX5jJjznoq9zfx6pJtTBrSi4smlPD6snKWbasKzti7qrw6JFTldCEWgXEDrTq6c+pY/nfWKhqaA5x/1ACG9ytkwuCe9OueR0+78c/2enjpe8fz3acWUlXfzIVHH8ofLjmayppGTrjn/Vb1AKE5ILA8TbBuCG44fTh3/HclQ/rkM6i3lTwf1b87c9ZVsvrdtUwe2pt+9vZnju7Hki37gqGv+y8dz0UTSvh0wx627K5jRXlVsBH/5/Tj2b7PigSM7FfIOWOtmRrOPKI/D8/ZQMDAy9cfz8Qh1n/4G8cO5tn5W/jRmSPxeoSbzxlNQW4WJ4/sG3xOiTOJpcNfLh/Prf9eFvTwRx3SnXPGHsLFEwfS0OxnzIAibjxjBHNLK/H5A5wyqpjxg3vy6FzL03/i2mP4/Zur+eqXSvAHDPfZ41jOHN2PwXZHgsOKC7n8mEH4A4affXlUMLqQKlIpHCWAu/9iGZZXEW+bkjhl+xtjygGMMeUi0i/SwUXkOuA6gMGDB7fxFCy+d+pwXlpYxj1vrObdVTu55vihZHs9fHV8SVA4ThlZzLdPOoyz7v+Qf3yyKVg2N9tDXraXU0cV8+HaCj7+xRnBkEu/7tYPvbyqoZXHMaRPATlZHt5Ybt25HndYbz5aV8lH6yoZN7AHS8uqWLxlH2eM7sf2ffWs3rGfIwYUsaq8mg/CetY4guA+xuGHdGfWsnIen7uxVbfQFxZsYdu+ei4/ZhBHD+zBeLtRXV9Rw99mr+cr4w5tVUffPmkYz823yl1/2vCQfE+uq5E/76gB/G7WalaWV/PjF5cwf8Mezhjdj/dXWyG3KUdafzqn/FlH9Of1ZVYD6faw3AncQld329997SiOPLQHv3xlGWAJtZuffnkUZ4zux7HDetPoC/Dpht0RRSNZzhnbnzdvOpnD+3dHRHjtByeHrM+JEJJyKMjN4vjD+jBvw24OKy5geHFhUDgu+lIJ004cRkmvbvz+zdWcOKIvJwzvw2HFBby2tJzJw3oHrw9Y+YxwnAb4+6eNAKy76oAxISO7Jw3pHUyqupPLp44q5sErJ/ClwT2DNx5XHz+Uc48cwC9eXsotU0YHG/JwJg/rzZ8vH8+0JxZw/WnDycnycGjPbvzl8vHkeD28tqyc15eWc+bofry3elercU5OCPGkEX24YvJg5m/Yw3dPGRZcP+3EYXy6cQ8bK2qYflrL1C2Thvbm+euOC6n7M0b354zR/QGrh9/Vj33GuWMP4Zihvalr8vHDM0cy7YShwd+dE5Y6elDPoGgA3Dn1SH7zlbHBm4DvntJy3JevPwEI/b0DTB1fwoTBvTj5XivEe2RJyzQ/edleZtkhT7fgFOVl8/oPT6JvYS79i/I4fXRLMzeqf3d65mdzzNAWuwDusUO86SCVwhFpspnwrknRtkmkbEyMMY8AjwBMmjSpXX3bhhcXMnlob15eXEbP/Gym27FkEaFnfjb76po5/JDuDOtb0Kqs8yd86KoJ7K1rDonTO3+M7fvqWyW9vR5hRHEhK8urEYF7Lz6aU+77gCyPMOOqicG7tnPHHkJtk487/ruSIw8tCnYJ/fk5hzO4dz6vLd3OWyusO6FsV+P1h68fzfpdNfw2QrjhP0u2k+URbr9wbIgADC8uZN6tZ9C/e+vJ6/Kyvdz79XG88vk2zjqif9S6HNgrn5e+dzyXPjyPj9ZV8v3ThvOTs0cx/ZnFfOO4wcG7XSesdvFE6674prCcg/uuuDA39I96wdEDWFVezYh+ha3mp8rN8nLcYX2C+zjt8Ij3HUkjIoy2Y+xtYcZVE5m5dDtHlfRARPjGsYN5a8UO7r9sfHCbt398Cv2659Etx4s/YPj1BWO47JhB0XcaBadOumV76V2Qw57aJibZD5h6aWFZyHiWnCwP548b0Gofxd1zg6GYWJx+eD/W/e+UEG/X6Rhx1pj+3DplNI2+AMu3V4U0jmDlgbplezn3yEPIy/Yy4+qJIesP6ZHHf244sVXHlHicNKIv91x0FKcebnnL+TlZ/CTs95Xl9bDgtrNCxgCB1fMxJ8o8Wo5XGolBvfP57JdnUt3go6Rnt4TsHBtlHrmzx0T/f6WLVApHGeD+VQ8Ewp8nGW2bnBhld4rIANvbGACEjtRJEb+YMppXPi/jkomDgl0OAZ759rG8sGALQ/tYMe/i7rkhsWcnFp+fk9VKHJw7wbomf6sfKMARA4pYWV5Nz27ZDO6Tzwc/O438HC/9i/K49oSh1Df5uWTSQESEgb3yObKkiMnDerNpdy3ft+/63QnabE/Ln7dbjpefnD2K655exFXHDaZ3fg7HDe/DtCcWsL/Bx1ElPSJOCjigR/Qf/Ykj+nLiiPjJuKNcd1zfPmkYWV4Pj34z9E754asn8tz8rZw2qh+r7zy3lS3ua+AOVYF1t3bnV4+Ma0dnokd+NlcfNyT4/a6vHsldYecwpE/LjYnXI3z7pGG0BxHhqJIebKisoX9RHiePsq6d07W2o8iOEjbJ9noY2MvyVub/8qxW63vm5/D5/5wdd3LKZETD2f7yyfGjEO7fWEfQryiPfm2/t+hUpFI4FgAjRWQYsA24HLgybJuZwI12DuNYoMoWhIoYZWcC3wTusd//k8JzCDJxSK+IdxRHlvTgrpKjgt9f/8FJvLd6F7f+2wqV5Hij/+gLcrPonpvF/kZfcHCem2tPGMrLi8tosHtouT2a2y8cG7KtcxdySVhXXCdZCS3jOBy+PPYQFv/67JDYsvMf/JLtqqeCbjleThrRl37dc0Om4XAzcUjvYIggz9O6Dt2NUWEE0e3qJNsYtpW7vnpkMPZelJfNHReOTVnf/7aQyRmNleik7B9njPGJyI3AW4AXeNwYs0JEptvrZwCzgPOAUqAOmBarrL3re4CXROTbwBbgklSdQ1voV5THBeMGtAhHjNi2tX0u+yt8EROnRw3swR0Xjm3XXDOH9GgRjkh3fuEJyauPG8In63dz4+kj2nzMRHjmO+HpruQZ3DufLXvqKEhRz5GDgfD8xDdPGJoZQ5QuRUpv1Ywxs7DEwb1shuuzAW5ItKy9fDdwZsda2rF0z8tmaJ98Nu2ui9hP303/ojzWV9RGFA5o/x/Z7XFECoeFc9v5Y9p1vHRy1hH9efzjjXpXqihp5sDz8TsJYw/twabddXE9DifPEU042ovTXe/EEX2CxzpQ+OV5ozl9dHFILxVFUVKPCkeKGFtSxOvLyhMKVUHL9O0dzYAe3Zj1w5OD40IOJLK8nuA4EkVR0ocKR4r46vgSdlU3MiRKH3cHp8vd2p01MbdrD2MOPUC6ciiK0ilQ4UgRh/bs1qrnUyTOGZv5PtmKoijJoMKRYXKzvMy4agJNUZ5HriiK0tlQ4egEnHtk65G5iqIonRV9kJOiKIqSFCociqIoSlKocCiKoihJocKhKIqiJIUKh6IoipIUKhyKoihKUqhwKIqiKEmhwqEoiqIkhVgzmx/Y2A+G2tzG4n2Byg40p6NQu5JD7UoOtSs5Oqtd0D7bhhhjWs0kelAIR3sQkYXGmEnxt0wvaldyqF3JoXYlR2e1C1Jjm4aqFEVRlKRQ4VAURVGSQoUjPo9k2oAoqF3JoXYlh9qVHJ3VLkiBbZrjUBRFUZJCPQ5FURQlKVQ4FEVRlKQ4qIVDRPJE5DMR+UJEVojIHfby3iLyjoiss997ucrcKiKlIrJGRM5Js123i8g2EVliv85Lp12uY3lF5HMRec3+ntH6imFXZ6mvTSKyzLZhob0s43UWxa6M15mI9BSRf4nIahFZJSLHd5L6imRXRutLRA53HXuJiFSLyE0pry9jzEH7AgQotD9nA/OB44B7gVvs5bcAv7c/jwG+AHKBYcB6wJtGu24HfhZh+7TY5TreT4DngNfs7xmtrxh2dZb62gT0DVuW8TqLYlfG6wx4EviO/TkH6NlJ6iuSXRmvL9cxvcAOYEiq6+ug9jiMRY39Ndt+GWAq1o8E+/2r9uepwAvGmEZjzEagFJicRruikRa7AERkIHA+8GjY8TNWXzHsikba7IpjQ0brLEnSYpeIFAGnAI8BGGOajDH7yHB9xbArGpm4jmcC640xm0lxfR3UwgHB8MYSYBfwjjFmPtDfGFMOYL/3szcvAba6ipfZy9JlF8CNIrJURB53uZ9pswv4M3AzEHAty3h9RbELMl9fYIn+2yKySESus5d1hjqLZBdkts4OAyqAJ+yw46MiUkDm6yuaXdA5fmMAlwPP259TWl8HvXAYY/zGmPHAQGCyiBwZY3OJtIs02vUQMBwYD5QDf0ynXSJyAbDLGLMo0SIRlqXTrozWl4sTjTETgCnADSJySoxt02lbJLsyXWdZwATgIWPMl4BarFBLNDJtV6bryzqYSA5wIfDPeJtGWJa0XQe9cDjYbuds4Fxgp4gMALDfd9mblQGDXMUGAtvTZZcxZqctKAHg77S4mOmy60TgQhHZBLwAnCEiz5D5+opoVyeoLwCMMdvt913AK7Ydma6ziHZ1gjorA8pcHva/sBrsTNdXRLs6QX05TAEWG2N22t9TWl8HtXCISLGI9LQ/dwPOAlYDM4Fv2pt9E/iP/XkmcLmI5IrIMGAk8Fm67HJ+CDZfA5an0y5jzK3GmIHGmKFYbvH7xpiryHB9RbMr0/UFICIFItLd+Qx82bYj07+xiHZlus6MMTuArSJyuL3oTGAlmf+NRbQr0/Xl4gpawlTO8VNXX6nM8nf2FzAO+BxYinXB/8de3gd4D1hnv/d2lbkNqyfCGmBKmu16GlhmL58JDEinXWE2nkZL76WM1lcMuzJeX1ix8S/s1wrgts5QZzHs6gx1Nh5YaNvwKtAr0/UVw67OUF/5wG6gh2tZSutLpxxRFEVRkuKgDlUpiqIoyaPCoSiKoiSFCoeiKIqSFCociqIoSlKocCiKoihJocKhKB2EiNTE3yrz+1SU9qLCoSiKoiSFCoeipAAR+bmILLAnv3Oep/J7Efm+a5vbReSn0bYP298AEZljP3NhuYicnL6zUZRQVDgUpYMRkS9jTeUwGWu08UR7AsEXgMtcm14K/DPG9m6uBN4y1sSXRwNLUncGihKbrEwboCgHIF+2X5/b3wuBkcaYx0Skn4gcChQDe40xW0Tkh5G2B+a49rkAeFxEsoFXjTFL0nAeihIRFQ5F6XgEuNsY83CEdf8Cvg4cguWBxNseAGPMHNsLOR94WkTuM8Y81cF2K0pCaKhKUTqet4BviUghgIiUiIjzIJ0XsGbw/TqWiMTbHnvZEKxnjvwd6yl0E1J/GooSGfU4FKWDMca8LSJHAPNEBKAGuAqr4V9hT2e+zbQ8oS3q9q7dngb8XESa7fXXpOt8FCUcnR1XURRFSQoNVSmKoihJocKhKIqiJIUKh6IoipIUKhyKoihKUqhwKIqiKEmhwqEoiqIkhQqHoiiKkhT/D9WVfO/03VzjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# density gradient and threshold\n",
    "var.differentiate('levels').plot()\n",
    "plt.hlines(threshold_rho,p1,p2,colors='r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ffdf322b-95c2-4893-8ba5-632a4a41d087",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f9173266f40>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGDCAYAAADEegxVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABNF0lEQVR4nO3dd3hUZdr48e89k95ogQgEKQIiICAgErEEe2d1ZdVV17osr67rrmuBdYu+ylpWd/W1/FzUVazY1oYVkFijCIIoHaWFIqEEEtKT+/fHOYljSJmEzJxJ5v5c11xz5pznnHM/c2bmntOeR1QVY4wxBsDndQDGGGMihyUFY4wxtSwpGGOMqWVJwRhjTC1LCsYYY2pZUjDGGFMrapKCiBwoIkUi4vc6FmNMaImIikh/r+Noi9ptUhCRdSJyQs1rVd2gqimqWuVlXA0RkXgReVxE1otIoYgsEpFT65Q5XkRWiEixiMwTkd4B00RE7hKRHe7jbhGRgOnrRKTETYxFIvJ+E/H0cddR7K7zhIBpIiI3i8gGEdkjIjNFJK2RZXUXkTdEZLP7Ze1TZ/rdIrLRXdZ6EbnZHX90QLx73XmLAh4HisgvROQzN86cetZ9nIh85S77exGZ1ES9/yAiW0Vkt4j8R0TiA6Z1FpFX3VjWi8gvm1hWY9vr9248e9z35V8iEtPCZTW17Rvclu70X7r12Ssir4lI54Bp8e77sMd9X64LmJYuIp+66ywQkVwRGVcnrttFZJP7fuaIyJB66jZAREpF5Jlg69zWich4t067RWRdI+WOdT/3t4ctOFVtlw9gHXCC13E0I95k4BagD06yPgMoBPq409OB3cBEIAH4B/B5wPy/AVYCmUBPYBkwuaXvB5AL/BNIBH4OFABd3WmXACuAXkAK8Dowo5FlZQBXAVmA1tQpYPrBQLI73BNYCpxTp0wfd96YOuNPAH4B/BXIqTMt1n3PfgMIcDhQBAxvIM6TgR+AIUAnIAe4M2D688ALbp2Pcpc9pIFlNbW9DgI6usOdgQ+A61q4rKa2fWPbcoj7OTvGrddzwMyAee8APnbfj0OArcAp7rQEd9v53Pf3Z8DOmm3kbpfNQD/A7y7rq3rq9767jmeCrXMQn18F+nv9vW4kvjHAxcAkYF0DZWKBxcDnwO1hi83rNydEb/jTQDVQ4v4I3EidHxX3C3878Jlb5k2gC/AssAf4koAfL2AQMNv90K8EfhGGeiwBfu4OTwI+C5iW7NZvkPv6M2BSwPQr6vxwrCPIpAAMBMqA1IBxH9f80AAvAzcETDsSKAWSmlhuDPUkhTplegLfADfWGf+T7VfPfFeyb1LIcOdJChj3JXBBA8t4Dvh7wOvjga0B73c5MLDO5+zOBpbV6PaqU7YLMAd4uCXLamzbB7Et/w48FzDtILeeqe7rTcBJAdNvIyBpBIz3AWe673c3d9xNwIsBZYYApXXmOx94EecP0TPB1jmIz3BtUgA6AE8B+cB64M+Az512KfAJcA+wC1gLnLo/39vmPHD+1KxrYNoU4G7gScKYFNrl4SNVvRjYAJypziGjuxsoej5Otu6J82XIBZ7A+ee2HPgbgIgk4ySE54BuwAXAw/XtCrvlH3Z3p+t7LAmmDiKSgfOFXuqOGgJ8HVDHvcB37vh9prvDdeN7VkTyReR9ERneyOqHAN+ramEDyxP3QcDreGBAU/VqiIhMEZEiIA/nB+C5li6rhqr+gPPv/jIR8YtIFtAb50cAETlKRAoCZqnvPcwQkS4426JKVVfVmV77Hrvb96j6llXP9qo5bLMH2A4MB/7dQFX2Z9s3tS3rLvs73OQnIp2AHo0su6YeS3D+FLwBPKaq29xJM4H+IjJQRGJx9jDfDZgvDfhf4I8tqHNzPICTGPoBxwK/Ai4LmH4Ezh+9dJwf4ccDD78FEpFZjXy3Z7Ugtnq5h8oux3l/wqpdJoVmeEJVv1PV3cA7wHeqOkdVK4GXgMPccmfgZPMnVLVSVb8CXgHOrW+hqnqVqnZs4DGsqaDcL9CzOIdkVrijU3B2pwPtBlIbmL4bSAn4cF+I82+7NzAPeE9EOjYQQlPrege40j1W3QHnHyFAUlN1a4iq3ukufyTOP/C662+p53EOLZXh/EO+WVU3uuv8RFU7BpSt7z3Ejaup9wR3+37SwLLqK/+cqqbhJJxHcA5d1Wd/tn1z5w2cnhLwut46uPUYBqQBv8RNuK4tOO/5Spx/+ROBPwRMvw14vGZ71NHk+xcMcS4sOQ+YqqqFqroOuBfnz2CN9ar6qDrnG2cA3XH2Mvehqmc08t0+ozmxNeH/gL+oalErLjMo0Z4UAr+EJfW8rvlS9AaOCPxXgPMje0BrByQiPpwfxXLgtwGTinC+eIHScI4H1zc9DSjSmn1p1U9VtURVi1X1Dpzjyke761wqP568PTqIdf0H58c2B2dPZp47Pk9+enJ4Kc2gjkU47/2tzZm3PiIyCOccwK+AOJx/mTeKyOkNzFLfewhOvZt6T5paVoPlVXU1zvv4cAuX1di2b+68gdOLAl43VYdSVX0emBKwF/o3nPM4vXDOC9wKfCAiSSIyAufQyb/qqW9TcTVHOs62Xx8wbj3O0YEaWwPqUewOptBKRORPAd+JR4IofybO4bsXWiuG5mjPSaE1m3/dCHxY519Biqr+T32FReQR+elVMkXB/FC6/+wex/mX8nNVrQiYvBTnEENN2WScQ15L65vuDjf2o6y4h4BUdYhbnxRV/didr5+IBP4rq12eqlar6t9UtY+qZrrjNwGbVPXjgGW1ZFcfnHMPB7Vw3kBDgZWq+p4b80rgLeDUBsrX9x7+oKo7gFVAjIgMqDO9ofe4qe1VV2N13p9t3+i2rGfZ/XAOBa5S1V04//ab87mKxTlMU1P2BVXNc/ewn8Q5YT0YyMbZc90gIluB64Gfi8hXQdY5WNuBCpw/djUOxPm8NpuIvNPId/ud+uZR1b8HfCcmB7Ga44HR4lzttRVnT+f3IvJ6S2JutnCcuPDigXPGPvDkWx/2PdF8ZcD024EnA16fAKxxh1Nx/l1cjPOhj8X5B3RIK8f8iBt3Sj3TuuLsPv8c51/XXfz0RPJknPMgPXGOAy/lx5OJBwLjcP4xJQA34Jx069LE+3ePW/5sfnrFSmecL6jgfMG/DXyvG1heAs65AsW5YiXBHe/DuXqmk7u8MTg/RL+rM/9Ptl/AeL+77MnAR+5wrDvtIJx/nMe5yz4IWAP8uoEYT8H51zjYjecDfnr10UycPaRk9/1s7OqjprbXlfx4Qnawu73+2cJlNbjtg9iWQ3AurDjardcz/PTqozuBD933Y5C7bWquPhqLcxVWHM6VTTfh/JPv4U7/G87hpAx3O18M7AU64hxqPCDgcQ/OBQxdg6zzpTRwgtadHnii+RngVZzvcW+cK+euDFjOJw3NG8LfJ59br1NxflsSgLiA35vA9+YFnD2qzqGMqTa2cKzEiwcwAedkcwHOv5A+tDApuK8PxvmXmQ/swPnBGNGK8fZ24yvF+SGreVxYJ6YVOIdXcvjp1VGCc5Jsp/u4GxB32hCcK5n2urHPBUY3EU8fdx0lOMeETwiYNtAdV+x+oOu9lLLO8rTuQ3/8crzrxlyE84/8TzWx14mnvqRwaT3LDtyOv8BJWoU4J7Hv4scrT47GOcwSuLzrcA4j7sG56CA+YFpn4DX3fdwA/LLOvEXA0UFuryfc9ezFuTLsH7iJ0p2+tDW2fVPb0p3+S7c+e3EuL+4cMC0e53DhHjfe6wKmHYtzMrjQXe+HwDEB0xOAh3ASyR7gK9yEUs/n4xYCrj4Kos5/AZ5t4vNWkxQ64SSGfJy9/r9S5+qjhuYN1QNnT6nu5zangbJPEnD1Ud3PLc735Z2A1+8Af2ppbDU/GsYY02aIc/Pltaq63OtY2htLCsYYY2q15xPNxhhjmsmSgjHGmFqWFIwxxtSypGCMMaZWg031tgXp6enap0+foMvv3buX5OTk0AUUgazO0cHqHB1aq84LFy7crqpd65vWppNCnz59WLBgQdDlc3JyyM7ODl1AEcjqHB2sztGhteosIusbmmaHj4wxxtSypGCMMaaWJQVjjDG12vQ5BWNMeFVUVJCXl0dpaanXodChQweWL4+uVi6aW+eEhAQyMzOJjY0Neh5LCsaYoOXl5ZGamkqfPn1ooHOysCksLCQ1tVl97rR5zamzqrJjxw7y8vLo27dv0Ouww0fGmKCVlpbSpUsXzxOCaZqI0KVLl2bv1VlSMMY0iyWEtqMl28qSgjGmTfH7/YwYMYIxY8YwfPhw/vnPf1JdXd2q63jkkUd46qmnAHjyySfZvHlzqy6/KStWrGDEiBEcdthhfPfddxx55JEArF+/nqFDh4Z03XZOwRjTpiQmJrJ48WIKCwspKSnhl7/8Jbt37+bWW/e7W+9akyf/2Gvmk08+ydChQ+nRo0erLR+gsrKSmJj6f4Jfe+01JkyYUFunzz77rFXX3ZiQ7SmISIKIzBeRr92O4W91x3cWkdkistp97hQwz1QRWSMiK0Xk5FDFZoxpH7p168b06dN58MEHUVWqqqq44YYbOPzwwxk2bBj//ve/gR/vBD733HMZNGgQF154YU0vZUyZMoXBgwczbNgwrr/+egBuueUW7rnnHl5++WUWLFjAhRdeyIgRI3jrrbc4++yza9c/e/ZszjnnnH3i6tOnDzfddBNjxoxhzJgxrFmzBoBLL72U6667jvHjx3PTTTexePFixo4dy7Bhwzj77LPZtWsXb7/9Nvfddx+PPfYY48ePByAlJWWfdTRU1/0Vyj2FMuA4VS0SkVjgE7dj63OAuap6p4hMAaYAN4nIYOB8nK4jewBzRGSgqlaFMEZjTAvd+uZSlm3e06rLHNwjjb+dOaRZ8/Tr14/q6mq2bdvG66+/TocOHfjyyy8pKytj3LhxnHTSSQAsWrSIpUuX0qNHD8aNG8enn37K4MGDefXVV1mxYgUiQkFBwU+Wfe655/Lggw9yzz33MHr0aFSVP/7xj+Tn59O1a1eeeOIJLrvssnrjSktLY/78+Tz11FP8/ve/Z9asWQCsWrWKOXPm4Pf7GTZsGA888ADHHnssf/3rX7n11lu57777mDx5MikpKbVJqj6PP/54vXVtzpVG9QnZnoI6ityXNZ3dK07fyTPc8TOAn7nDE3A6DC9T1bU4HayPCUVse0ormLdiG5+t2c7C9bv4dtNu1mwrYuPOYvILy9hTWkFZZRUt6pUuNxfuuMN5NsaERc139f333+epp55ixIgRHHHEEezYsYPVq1cDMGbMGDIzM/H5fIwYMYJ169aRlpZGQkICV155Jf/9739JSkpqdD0iwsUXX8wzzzxDQUEBubm5nHrqqfWWveCCC2qfcwN+DyZOnIjf72f37t0UFBRw7LHHAnDJJZfw0UcfBV3nxuq6P0J6TkFE/MBCoD/wkKp+ISIZqroFQFW3iEg3t3hP4POA2fPccXWXOQmYBJCRkUFOTk7Q8RQVFZGTk8P3BVX87+dNX6YlQKwPYv0Q6xNifRDnDjvP7nh33CHrl3HTA1OJqapAY2P5+t572TOkef96WltNnaOJ1Tl0OnToQGFhIQDXZR8YknXULL+pMlVVVRQWFrJ27Vp8Ph+JiYlUVFRw1113ccIJJ/yk/Mcff4zf769ddlVVFUVFRZSUlDB37lxycnJ46aWXuP/++5k1axZlZWXExsbWrmfv3r21806cOJHzzjsPgAkTJlBSUrJPfKpaO09FRUVtzBUVFfh8PgoLCyksLERVa5dbVFREdXU1hYWFP1l/YJ2rq6tryzRU17rvX2lpabM+GyFNCu6hnxEi0hF4VUQaO21e37VT+/xVV9XpwHSA0aNHa3NaDKw5rnh4WSWHjSykrLKa0oqqnzyXVVZTVndcRRWlFdWUVe77XFZZTWFFFWVl1Ry55Ft8lRX4tBqtrGTknj3gcSuO1pJkdAhXnZcvXx4RN4ylpqZSWFhIaWkp119/Pddccw1paWmcfvrpzJgxgzPOOIPY2FhWrVpFz549SUpKIiYmpjb2uLg4EhISEBGqq6s599xzOe644+jfvz+pqanEx8cTHx9PamoqHTt2pLq6unbe1NRUMjMzueeee5g9e3a974eI8NZbbzFlyhSeeeYZjjzySFJTU4mNjSUxMZHU1FRSU1Pp3Lkzixcv5uijj+bVV19l/Pjx+6w/sM4+nw+fz0dqamqDda3btHZCQgKHHXZY0O9tWK4+UtUCEckBTgF+EJHu7l5Cd2CbWywP6BUwWyYQkuvAkuNjOOzATk0XbK5jEqn8+Dkqy8vxxcYhUfbDZEw4lJSUMGLECMrKyoiLi+Piiy/muuuuA+DKK69k3bp1jBw5ElWla9euvPbaaw0uq7CwkAkTJlBaWoqq8q9//WufMpdeeimTJ08mMTGR3NxcEhMTufDCC8nPz2fw4MENLrusrIwjjjiC6upqnn/++XrLzJgxg8mTJ1NcXEy/fv144okngn4fmlvXoKlqSB5AV6CjO5wIfAycAfwDmOKOnwLc7Q4PAb4G4oG+wPeAv7F1jBo1Sptj3rx5zSrfEm//+xW965hf6dZ3Pwj5uoIRjjpHGqtz6Cxbtiws6wnGnj17PFv31VdfrY899liD03v37q35+fmtvt6W1Lm+bQYs0AZ+V0O5p9AdmOGeV/ABL6rqLBHJBV4UkSuADcBENzktFZEXgWVAJXC1tsErj95O7cv8ky/mhpOyvQ7FGBMCo0aNIjk5mXvvvdfrUEIiZElBVZcA+xzIUtUdwPENzDMNmBaqmMIh1idUVClV1UqM35oDMKa9WbhwYZNl1q1bF/pAQsSauWhlpw/rzs695Tw/f4PXoRhjTLNZUmhlxw3qxvDMDrz81SavQzHGmGazpNDKRIQj+6ezdNNuSiva3CkRY0yUs6QQAqpQWa3s2FvudSjGGNMslhRCwOeeX773vZXeBmJMOzRt2jSGDBlCVlYWI0aM4IsvvgDgvvvuo7i42OPofuq0007bpz2lurxomrsxlhRC4IqjnAap/rtoU8vaTzLG1Cs3N5dZs2bx1VdfkZuby5w5c+jVy7nnNRKTwttvv03Hjh0bLWNJIQp0SYkHwO8TCoorPI7GGI+1YiORW7ZsIT09nfh45zuWnp5Ojx49+L//+z82b97M+PHja5ubfv/998nKymLkyJFMnDiRoiKnfc6GmrV+6aWXGDp0KMOHD+eYY47ZZ905OTkcc8wxnH322QwePJjJkyfXdu7z/PPPc+ihhzJ06FBuuumm2nn69OnD9u3bWbduHYcccgi//vWvGTJkCCeddBIlJSX7NM1dUlJSb1PeYdXQXW1t4RGJdzTXGPv3OXr1swvDtr6G2N290SFi72j+7DPVxERVv995/uyz/Vp/YWGhDh8+XAcMGKBXXHGF5uTk1E4LvIs4Pz9fjz76aC0qKlJV1TvvvFNvvfXW2nK33367qqrOmDFDTz/9dFVVHTp0qObl5amq6q5du/ZZ97x58zQ+Pl6/++47rays1BNOOEFfeukl3bRpk/bq1Uu3bdumFRUVOn78eH311Vd/EtPatWvV7/frokWLVFV14sSJ+vTTT6uq6rHHHqtffvmlqqru2LFDBw4cqNXV1fXGEY47mm1PobXl5rL2j3+m96qvSYrzex2NMd7KyYHycqiqcp73syXXlJQUFi5cyPTp00lPT+e8887jySef3Kfc559/zrJlyxg3bhwjRoxgxowZrF+/vnZ6fc1ajxs3jksvvZRHH32Uqqr6rxwcM2YM/fr1w+/3c8EFF/DJJ5/w5Zdfkp2dTdeuXYmJieHCCy+stwnsvn37MmLECMC5K7q+G9ya25R3KFh3nK0pN5eq8cfRq7ycJ/wx/HfYk3y2pie905PpnpaAz2d3OJsok50NcXFOQoiLa5VWg/1+P9nZ2YwaNYrRo0czY8YMLr300p+UUVVOPPHEBhuiC+zQvmb4kUce4YsvvuCtt95ixIgRLF68mC5dujQ4X81rDfK8Yc0hr5o61NfkdkxMDPPnz2fu3LnMnDmTBx98kA8++CCo5bcWSwqtKScHqajAr9VQVcmm197l5h/SAIjxCakJMaQmxJISH0NKQgyp8TGkJjjDKfGx7vQYZ3p8DJkrvmbQyoX4xo+HrCyPK2dMC2Rlwdy5zh5CdvZ+f45XrlyJz+djwIABACxevJjevXsDPzannZ6eztixY7n66qtZs2YN/fv3p7i4mLy8PAYOHAjACy+8wJQpU3jhhRfIcmP67rvvOOKIIzjiiCN488032bhx4z5JYf78+axdu5bevXvzwgsvMGnSJI444giuvfZatm/fTqdOnXj++ee55pprgq5TTdzg9KlQXFzMaaedxtixY+nfv/9+vV8tYUmhNWVn44t3/hX54+K47ObLOKrfUNbvLCZvVzF7SiopKquksLSSwtIKtu4pZfW2mnEVVFT9+I9j5KblPDvzZrS6EuLjnS+WJQbTFmVltdpnt6ioiGuuuYaCggJ8Ph8DBw5k+vTpAEyaNIlTTz2V7t27M2/ePJ588kkuuOACysrKALj99ttrk0J9zVrfcMMNrF69GlXl+OOPZ/jw4fVUJYspU6bwzTff1J509vl83HHHHYwfPx5V5bTTTmPChAlB1ymwae533nmnyaa8Q66hkw1t4RGRJ5o/+0z1739v0Qm1kvJKzS8s1bX5Rbrlpr9qpc+nClrt9zvLbAE76RodIvZEcwi1tOnsljZrPW/evNqT0l5p601nR6f9+FeUEOsnIdZPeko8TDiFivvuRt29Dq97cDPGRAdLCpEqK4vvn3+N1+97jkMvmsCpdujImFbR0mats7Ozo6KbV0sKEazX6cfz9GLIKEzg5Gq1q5eMMSFn9ylEsKS4GC4/qi9rthVx7+yVVFRVex2SMdZ0SxvSkm1lSSHCXTS2N6cMOYCH5n3HA3NXex2OiXIJCQns2LHDEkMboKrs2LGDhISEZs1nh48iWW4uXebN44IeQ3iXGNbuiKzGvkz0yczMJC8vj/z8fK9DobS0tNk/eG1dc+uckJBAZmZms9ZhSSFS5eaixx9PdWkZY/wxjLxgGr+anEVBcTmpCbH47fyC8UBsbCx9+/b1OgzAaaDusMP26Qa+XQtHnS0pRKqcHLSsnBj37uix679h4iM/tjKZGh9DWmKs80iIoYM7nBIfQ3K8n+T4GJLjYti4qYKSb7aQFB9DSryfpDjnbumkOKdMfIxvn1v3jTHRy5JCpMrORuLj0PJyJDaWYydNpN+g4ewpqWB3SQV7St3nkkr2lFawYWcxu0sq2FtWyd7yKqqqfzzm+9g3XzW4Gr9PSIrzkxIfQ1pCLF1S4uiSEk+X5DjSA4ZrnpPjY0iM85MQ4yPGv5+npHJzW635A2NM67CkEKmyshC3zRh/djZHZGVxRJCzqiplldXsLavkg48+Zehhoykur2RvWVVt0nCeKykuq6KorJLi8koKiivYubecbzftZntRGYWllY2uJ9YvJMT4SYjzkxDrI9G9+a52XIzPTSB+EuP8xAeUOXDl15z0uwvxV1YgcXHWjIcxEcKSQiRr4d3RIlJ7d3TXJB+HdE9r0erLKqvYubecHUXlbC8qY0dROcXllZRWVFNSUUVpRZX7XE3pT15Xsbukgm11ypRUVFFe6VxWe1Xu65xUXo5oNdXl5fhyciwpGBMBLCmYBsXH+OneIZHuHRJbbZlV1UpZZRWVn6TiO/0lqsrLKBc/VVlHkdJqazHGtJQlBRNWzjmMGDjuGPhgLl888Qp3lRzAFZ37c7rXwRljLCkYD2Vlcc0He0lLjCX74K5eR2OMwe5oNh7r2SmR4rJKu+/CmAgRsqQgIgkiMl9EvhaRpSJyqzv+FhHZJCKL3cdpAfNMFZE1IrJSRE4OVWwmclww5kA27y7l8+93eB2KMYbQHj4qA45T1SIRiQU+EZF33Gn/UtV7AguLyGDgfGAI0AOYIyIDVbX+HrRNu7DBbbpjWGZHbwMxxgAh3FNwO/gpcl/Guo/GWtGaAMxU1TJVXQusAcaEKj4TGYrLnZwfF2NHMo2JBCH9JoqIX0QWA9uA2ar6hTvptyKyRET+IyKd3HE9gY0Bs+e540w7NiyzAwAfr/K+gTVjDEg4msAVkY7Aq8A1QD6wHWev4Tagu6peLiIPAbmq+ow7z+PA26r6Sp1lTQImAWRkZIyaOXNm0HEUFRWRkhJdV8NHep0/3VTBo9+UM3VMAgd39rfKMiO9zqFgdY4OrVXn8ePHL1TV0fVNC8slqapaICI5wCmB5xJE5FFglvsyD+gVMFsmsLmeZU0HpgOMHj1am9M9Xk5OTlR0pxco0uu8+qPv4ZvlnHnckfTo2Do3yUV6nUPB6hwdwlHnUF591NXdQ0BEEoETgBUi0j2g2NnAt+7wG8D5IhIvIn2BAcD8UMVnIsPJQw4AYNpbyz2OxBgDod1T6A7MEBE/TvJ5UVVnicjTIjIC5/DROuA3AKq6VEReBJYBlcDVduVR+9clJQ6AL9bupKpa7X4FYzwWsqSgqkuAfXqDUNWLG5lnGjAtVDGZyOMTIT0ljh4do6sHLWMilV0HaDyVGOfnt+P7syRvN8u37PE6HGOiniUF47kfCsuI8Qn9uiZ7HYoxUc+SgvHc3rJKYve3FzdjTKuwb6LxXHJ8DCUVVXy1vsDrUIyJepYUjOdi3CuOSirsYjNjvGZJwXiqaN7HHPjv+xm5aTnJ8a1zR7MxpuWskx3jndxc4k89ibPLyznDH8vH2Qcxt/xIuqbG0y01gfSUOGLsXIMxYWVJwXgnJwd/RQU+rYaqChY/9RoPr46rnSwCXZLj6JwcR8ekODonxdEpOZaOSXF0SordZ1x6SjxpCfaRNmZ/2DfIeCc7G198HJSX44+L47KbL+PEQcPZVlhGfmGZ+1zKzr3l7Cqu4Lv8Inatr6CguJzK6vobcoyL8ZEao2Qu/ZSuKXF0TY0nPSWerqnxDOiWytCeaaQuWgA5OZCdDVlZYa2yMZHOkoLxTlYWzJ0LOTlIdjZds7IIpqdmVaWorJKC4gp2FTsJY9fecrYXlZFfVMa3qzfgT4xlU0EpX+ftZkdRGTU5ZOSm5Tz3wp+JrarAFx+PzJ1ricGYAJYUjLeyspr9oywipCbEkpoQS6/OSftMz0n8gezsH/tnqqpWdhSVsXxrIXF3f0BsZQV+rUbLy509BksKxtSys3im3fP7hG5pCRw7sCtZV5wL8XFUig+NjXMOIRljallSMNElK4tH//Yo/zzmIna+8bbtJRhThx0+MlHnw8792XZmJjeemO11KMZEHNtTMFGnoqoan1i/DcbUx5KCiTorthZai6zGNMCSgokq2/aUUl5ZzZbdpV6HYkxEsqRgosrd762kvKqay8b18ToUYyKSJQUTVbokO81onDq0u8eRGBOZLCmYqHJAB6cv6C/W7vQ4EmMikyUFE1XOO7wXB6QlcP+cVV6HYkxEsqRgokpSXAxH9u/C0s172Lm33OtwjIk4lhRM1Dl6QDplldXc/e4KqhtobdWYaGVJwUSd4w/JAGDmlxtZtmWPx9EYE1ksKZiok5YQy1H90wHokBjrcTTGRBZLCiYq/U/2QQC88fVmjyMxJrJYUjBR6Yi+nTmqfzr3z1nN2u17vQ7HmIhhScFEpRi/j3/+YjjxMT6mvbXM63CMiRghTwoi4heRRSIyy33dWURmi8hq97lTQNmpIrJGRFaKyMmhjs1Et25pCUwc3Ys5y7fx8ep8r8MxJiKEY0/hWmB5wOspwFxVHQDMdV8jIoOB84EhwCnAwyLiD0N8JopNzu5H7y5J/OW1b70OxZiIENKkICKZwOnAYwGjJwAz3OEZwM8Cxs9U1TJVXQusAcZgTAh1S03gkqw+rNtRzLebdnsdjjGeC/Wewn3AjUB1wLgMVd0C4D53c8f3BDYGlMtzxxkTUj06Ou0hVdqNbMaErjtOETkD2KaqC0UkO5hZ6hm3z7dURCYBkwAyMjLIyckJOqaioqJmlW8PrM5NW7O9CoAPPltAQUbb7KHWtnN0CEedQ/kNGAecJSKnAQlAmog8A/wgIt1VdYuIdAe2ueXzgF4B82cC+1xErqrTgekAo0eP1uzs7KADysnJoTnl2wOrc9PGVVXz5IoP+DA/gcvOGEMnt3nttsS2c3QIR51DdvhIVaeqaqaq9sE5gfyBql4EvAFc4ha7BHjdHX4DOF9E4kWkLzAAmB+q+IypEev3cetZQ1i+pZBf/DvXGsozUc2L+xTuBE4UkdXAie5rVHUp8CKwDHgXuFpVqzyIz0Sh04d1Z8blY/guv4j/fLLW63CM8UxYkoKq5qjqGe7wDlU9XlUHuM87A8pNU9WDVPVgVX0nHLEZUyProC6M65/O619vQtVOOpvoZHc0GxPgnJE92bizhJlfbmy6sDHtkCUFYwKcOawHRw9I56+vf8vM+RuosstUTZSxpGBMgBi/jwcvGMnwzI5M+e833PzqN16HZExYWVIwpo4OSbG8NDkLgG2FZR5HY0x4WVIwph4iwtCeaeRbUjBRxpKCMfX4dtNuvt20h6yDungdijFhZUnBmHqI2+hKny7J3gZiTJhZUjCmHklxTgswBSV2d7OJLpYUjKlH3/RkMtLirTltE3UsKRhTj+pqZUdROcXl1tKKiS6WFIyph88n/OLwXny2Zgff5xd5HY4xYWNJwZgGXHlUX1ISYjhv+ues/qHQ63CMCQtLCsY0oF/XFF6YNJb8wjKmf/S91+EYExaWFIxpRHJ8DCKwfmex16EYExaWFIxpRMekWIZldmT+2p12JZKJCpYUjGlEUlwMk47uB8CSPEsKpv2zpGBMEzotWcBVuS+y7s3ZXodiTMjFeB2AMREtN5exvz6PMWVlVHw2k6sKSqgYM5aBGSkc2DmJrqnxdEtNoFtqPJ2T44jx2/8s07ZZUjCmMTk5SEU5MVqNr7qKk/NX8MD2Q5m3YhuVdTrg8Ql0To6nW2o8XVPj6ZQUS8ekODokxtIhMZaOSc7DeR1XOxxricREEEsKxjQmOxuJi4PycnxxcUz4/YVMyMqivLKabYWlbCssI7+wrPY5v7C09vXa7XspKC6nsKySxrp8Tk2I4cDOSfTukkTvLsn06ZJEv64pjO7dCalpmc+YMLGkYExjsrJg7lzIyYHsbOc1EBfjI7NTEpmdkppcRFW1UlhaQUFxBQUlFewuqaCguJzdJRXsLq5ge1EZ63cWs2JLIe8v/aF2D+SRi0ZyytDuIaycMfuypGBMU7KyapNBS/h9QsekODomxTVZtrKqmitmLODDVfkMzEht8TqNaSk7mGlMBHn047V8uCofgFlLtvDpmu1sL7Le30z42J6CMREk++CufLhqG9sKy/jn7FW14zsnxzGgWwoDMlLo0yWZnh0T6dkpkR4dE+mS3PQeiDHBsqRgTAQ5pHsaMyc5h6p2F1ewZFMBq34oYvUPhazeVsTrizdTWFr5k3niY3x0ilPGbl3E1eP7M8AOO5n9YEnBmAjVISmWowd05egBXWvHqSq7SyrI21XC5oISNhU4z4tWbeC1xZt54+vNvH71URya2cHDyE1bZknBmDZE5MeT1kN7/vjDn5O8jdL0QUx+ZiHLt+6xpGBazE40G9NOnHBIN9JT4vlk9XavQzFtWMiTgoj4RWSRiMxyX98iIptEZLH7OC2g7FQRWSMiK0Xk5FDHZkx7EuP3UVxeSUlFFeTmwh13OM/GNEM4Dh9dCywH0gLG/UtV7wksJCKDgfOBIUAPYI6IDFRV6yTXmCB1S41nx+wcSq/8M3FVFUh8PDJ37n7dZ2GiS5N7CiKSJSIPicgSEckXkQ0i8raIXC0ijR64FJFM4HTgsSBimQDMVNUyVV0LrAHGBFMJY4zjxclZ3Jy8jdjKCnzV1VSVlpH36tteh2XakEaTgoi8A1wJvAecAnQHBgN/BhKA10XkrEYWcR9wI1BdZ/xv3STzHxHp5I7rCWwMKJPnjjPGBKlbagKjLjkbf0I86vdTGRPLTds7s2DdTq9DM22EaCMtdYlIuqo2etaqoTIicgZwmqpeJSLZwPWqeoaIZADbAQVuA7qr6uUi8hCQq6rPuPM/Drytqq/UWe4kYBJARkbGqJkzZwZd2aKiIlJSUoIu3x5YnaND3TqnLV1Kx8WL2Th4GDcUHsT2EuVn/WM566D2c6ObbeeWGz9+/EJVHV3vRFVt8gEkAz53eCBwFhDbxDx34PzbXwdsBYqBZ+qU6QN86w5PBaYGTHsPyGpsHaNGjdLmmDdvXrPKtwdW5+jQWJ0L9pbrJf/5QvtOmaVbd5eEL6gQs+3ccsACbeB3Ndirjz4CEkSkJzAXuAx4srEZVHWqqmaqah+cE8gfqOpFIhLY7OPZwLfu8BvA+SISLyJ9gQHA/CDjM8Y0oENSLH89YzDVCq8t2uR1OCbCBZsURFWLgXOAB1T1bJxzCy1xt4h8IyJLgPHAHwBUdSnwIrAMeBe4Wu3KI2NaRZeUeERgwfpdXodiIlywl6SKiGQBFwJXNHNeVDUHyHGHL26k3DRgWrDLNcYEp0NiLKMO7MSHq/KpqKq23t5Mg4L9ZFyLc8z/VVVdKiL9gHmhC8sY09qG9uxAeWU1b3692etQTAQL5j4FP3Cmqp6lqncBqOr3qvq7kEdnjGk1mZ0SAbjnvZUUFJd7HI2JVE0mBfe4/qgwxGKMCaHLx/XlxlMOZvPuUl5emOd1OCZCBXv4aJGIvCEiF4vIOTWPkEZmjGlVPp/wm2MOYkiPNKa9vZzdJRVeh2QiULBJoTOwAzgOONN9nBGqoIwxoeH3CVeP748qHPH3OTw0b43XIZkIE9QVRKp6WagDMcaEx2mHdmfWNUcx7a3l3Pv+SgZ3T2P8oG5eh2UiRFB7CiKS4DaA97DbXtF/ROQ/oQ7OGBMaQ3t24OELR3JI9zQmPb2ANdsKvQ7JRIhgDx89DRwAnAx8CGQC9ikypg3rlBzHU5ePISkuhr+8ttTrcEyECDYp9FfVvwB7VXUGTnPYh4YuLGNMOHRJiec3x/Yj9/sd5Kzc5nU4JgIEmxRqLlMoEJGhQAecxuyMMW3chWN60y89mRtfXlLTGKWJYsEmheluvwd/wWm4bhlwV8iiMsaETYekWCYd049thWVM/+h7SwxRLtirj2p6TvsQ6Be6cIwxXjh3VCazlmzhjndWkBwfw0Vje3sdkvFIsFcfdRGRB0TkKxFZKCL3iUiXUAdnjAmPGL+Pa47rD8CTn62zvYUoFuzho5nANuDnwLk4Pae9EKqgjDHhN6ZvZy49sg9rthWx0JrYjlpB39Gsqrep6lr3cTvQMYRxGWPCTES49vgBpKfEc+MrSygqq/Q6JOOBYJPCPBE5X0R87uMXwFuhDMwYE36dkuN48JeHsX5HMTe9YlcjRaNGk4KIFIrIHuA3wHNAmfuYidtjmjGmfRnbrwvXnTiQt5Zs4ePV270Ox4RZo0lBVVNVNc199qlqrPvwqWpauII0xoTXpGP6kZYQwy1vLmXjzmKvwzFh1NSewsjGHuEK0hgTXrF+H/dfcBj5hWX86dVvvA7HhFFT9ync6z4nAKOBrwEBhgFfAEeFLjRjjJfGH9yN8Qd346PV+agqIuJ1SCYMmjp8NF5VxwPrgZGqOlpVRwGHAdYQuzHt3OF9O1NQXMF7S7d6HYoJk2CvPhqkqrX7kKr6LTAiJBEZYyLGBYf3YtABqdw2azkVVdVeh2PCINiksFxEHhORbBE5VkQeBZaHMjBjjPdi/D7OGdmTTQUlbN1d6nU4JgyCTQqXAUuBa4Hf4zSIZ72xGRMFjujrtGjzxKfrvA3EhEWwDeKVAv9yH8aYKDK8V0cuGnsgT3y2lovGHki/rileh2RCqKlLUt8UkTNFJLaeaf1E5H9F5PLQhWeMiQTXHj+QWL+P/5fznTMiNxfuuMN5Nu1KU3sKvwauA+4TkZ1APs7lqX1xrj56UFVfD22IxhivdU2N5+QhB/DSwjxKP/qEf06/Hn9lBcTHIXPmIEce6XWIppU0mhRUdStwI3CjiPQBugMlwCpVtdscjYki084eysgDO5J07yykohyfVlNZWsZDf3mURRfFcPfPh9EtLcHZe8jJgexsyMryOmzTTEGdUwBQ1XXAuuYsXETWAYVAFVCpqqNFpDNOs9t93OX9QlV3ueWnAle45X+nqu81Z33GmNBJS4jlsnF9wfcr9N0ZaHk5xMQSd+Lx5H63g9P+72POq8zjD3de5exFxMUhc+daYmhjgr36aH+MV9URqjrafT0FmKuqA4C57mtEZDBwPjAEOAV4WET8YYjPGNMcWVnI3LnIbbcRM+8D/mfKRTx+yeGM659O6uefQnk5UlVFVVkZq1540+toTTMFvafQiiYA2e7wDCAHuMkdP1NVy4C1IrIGGAPYmSxjIk1W1k/2AI4akM5RA9Kh9xXoR89RXV5OpcTw513pPFxURnpKvIfBmuYIek9BRBJF5OBmLl+B990uPCe54zJUdQuA+9zNHd8T2Bgwb547zhjTVrh7Eb7bbmPjS2+wKHMQf3hhsddRmWYIak9BRM4E7gHigL4iMgL4X1U9q4lZx6nqZhHpBswWkRWNraaecfv08OEml0kAGRkZ5OTkBFEDR1FRUbPKtwdW5+gQcXV29yKOOED5ePV2Hn9tLgd1bN2jwRFX5zAIS51VtckHsBDoACwKGLckmHkDyt8CXA+sBLq747oDK93hqcDUgPLvAVmNLXPUqFHaHPPmzWtW+fbA6hwdIrXOP+wu0SPvmKtH3/WBVldXt+qyI7XOodRadQYWaAO/q8EePqpU1d3NSTYikiwiqTXDwEnAt8AbwCVusUuAmvsc3gDOF5F4EekLDADmN2edxpjI0i0tgWMGpvPDnlKsZ8+2Idik8K2I/BLwi8gAEXkA+KyJeTKAT0Tka5wf97dU9V3gTuBEEVkNnOi+RlWXAi/itKv0LnC1qlY1u0bGmIiS2SmJsspqHvjAWttvC4JNCtfgXCpahtNX826chvEapKrfq+pw9zFEVae543eo6vGqOsB93hkwzzRVPUhVD1bVd1pUI2NMRLny6L4A/GvOKkor7H9epGsyKbj3Cryhqjer6uHu48/qNJJnjDGNio/xM6SH06X7ve+v9Dga05Qmk4J7CKdYRDqEIR5jTDvUo2MiAEN72s9IpAv25rVS4BsRmQ3srRmpqr8LSVTGmHblb2cO5vPvd/D+0h+YMMJuP4pkwSaFt9yHMcY0W2anJDolxfHWN1uYsrOYXp2TvA7JNCDYTnZmhDoQY0z79tczBvOHFxaTfU8OL03OYuSBnbwOydQjqKuPRGStiHxf9xHq4Iwx7ccJgzOYdEw/qqqVpZv3eB2OaUCwh49GBwwnABOBzq0fjjGmPTtqQDr3zl5FabldmhqpgtpTcO8tqHlsUtX7gONCG5oxpr05tGcHYv3C1j12RXukCrZBvJEBL304ew6pIYnIGNNuxfh9HDuwGy8u2MhV2QfRxZrUjjjBHj66N2C4ErfHtFaPxhjT7l009kDmLP+B1xdv5vKj+nodjqkj2KuPxoc6EGNMdNhdUgHAByu2WVKIQMFefXStiKSJ4zER+UpETgp1cMaY9ues4T3IPrgrn363nbJKO+EcaYJtEO9yVd2D0/x1N+Ay3NZNjTGmOUSEXx/dD1V4YO6amv5TTIQINinU9Ip2GvCEqn5N/T2lGWNMk8b1T+fM4T14cN4aLnvyS6/DMQGCTQoLReR9nKTwntt5TnXowjLGtHf3nTeCnh0TyVmZT2WV/ZxEimCTwhXAFOBwVS0GYnEOIRljTIv4fcJNpw4C4N8fWQMJkSLYpJCF05dygYhcBPwZp6MdY4xpsbOG92BozzReW7TJ61CMK9ik8P9w+lQYDtwIrAeeCllUxpiocfygDFZvK2Leym1eh2IIPilUqnOJwATgflW9H7uj2RjTCq4afxADuqXwt9eXUlVtVyJ5LdikUCgiU4GLgbfcLjpjQxeWMSZaxMf4mXRMPzbsLGbxxl1ehxP1gk0K5wFlOPcrbAV6Av8IWVTGmKgyflA30hJieOCDNV6HEvWCbSV1K/AKUNN61Xbg1VAFZYyJLukp8ZwwOINFGwoorbC7nL0UbDMXvwZeBv7tjuoJvBaimIwxUeikwQewu6SCpZvtwkYvBXv46GpgHLAHQFVX4zR3YYwxrWJIjzQAVv1Q5HEk0S3YpFCmquU1L0QkBrDLBIwxraZnx0S6JMcx47N17C2r9DqcqBVsUvhQRP4EJIrIicBLwJuhC8sYE218PuGf541gxdZC7np3hdfhRK1gk8JNQD7wDfAb4G2cu5qNMabVjOnTmeQ4Pzv2ljdd2IREk53siIgPWKKqQ4FHQx+SMSZard2+l73lVWzdbX04e6XJPQVVrQa+FpEDm7twEVknIt+IyGIRWeCOu0VENrnjFovIaQHlp4rIGhFZKSInN3d9xpi2bXCPNCYd04+F63fx6ZrtXocTlYLto7k7sFRE5gN7a0aq6llBzDteVetu3X+p6j2BI0RkMHA+MAToAcwRkYGqahctGxNFrs7uz4Ln32L1H95g3J8uhawsr0OKKsEmhVtDGoVjAjBTVcuAtSKyBhgD5IZh3caYCJE/ex7PzryZ2KpKKmc/xfoX36DnqceREOv3OrSo0OjhIxFJEJHfAxOBQcCnqvphzSOI5SvwvogsFJFJAeN/KyJLROQ/ItLJHdcT2BhQJs8dZ4yJIj0XzyeuqpIYrYbycl659xmG/u09/vji13auIQya2lOYAVQAHwOnAoOBa5ux/HGqullEugGzRWQFTjPct+EkjNuAe4HLqb97z33uhXCTyySAjIwMcnJygg6mqKioWeXbA6tzdGhPdU5L78TwuFiqKyogJpYex4/i2HQ/ry/K482v8zi5dyzjesaQosXtps7BCsd2lsY6zRaRb1T1UHc4BpivqiNbtCKRW4CiwHMJItIHmKWqQ91WWFHVO9xp7wG3qGqDh49Gjx6tCxYsCDqGnJwcsrOzWxJ+m2V1jg7trs65uZCTA9nZtecUNu4s5s53VvDWN1sA6JXq46FLshiW2dGzMMOttbaziCxU1dH1TWvq6qOKmgFVbdYthiKS7PbljIgkAycB34pI94BiZwPfusNvAOeLSLyI9AUGAPObs05jTDuRlQVTp/7kJHOvzkk8dOFIPp1yHH85YzAbC6u57IkvrQ+GVtbU4aPhIrLHHRacO5r3uMOqqmmNzJsBvCoiNet5TlXfFZGnRWQEzqGhdTg3w6GqS0XkRWAZUAlcbVceGWPq6tkxkSuO6st97y9jx95ySiuqSI4P9poZ05RG30lVbfHpflX9Hhhez/iLG5lnGjCtpes0xkSPiw+J5+Gvy/jtc1/xt6576PPNlz853GRaxtKrMaZNGn2An4kxmax7cw4Zz/+JqqpKNC6OivfeJ/HYo70Or80Ktu0jY4yJKD4R/jFxOE/2Kya+uhK/ewnr4qdf8zq0Ns2SgjGmTUs++QR88fGo30+FP5bZ3QZ7HVKbZoePjDFtW1YWzJ2L5OTwQmxvZuzsyE2VVcTH2B3QLWF7CsaYts+9hHXvqDFUVSuLNhR4HVGbZUnBGNNu/CqrN6kJMby0IM/rUNosSwrGmHYjNSGWcQelM2/lNoqsS88WsaRgjGlXfpXVm517y/nji4vZMTsH7rjDaTbDBMVONBtj2pXhvTpyUNdk8t/PIWnmzVRVVeJLiEfmzrUb24JgewrGmHYlOT6GOdcdy78PLKq9f6GqtIyn7pjB7bOWsSSvwOsQI5olBWNMuyMidD3j5Nr7FzQujpWDRvLU5+uZ+EiuJYZGWFIwxrRPNfcv3HYbsfM+YNrdk8idchypCbHcN2e119FFLDunYIxpv7KyfnIeoUtKPMcMSGf2sh/YXVxBh6RYD4OLTLanYIyJKp2S4ygsq+Srjbu8DiUiWVIwxkSVQQekAvDG4s0eRxKZLCkYY6LKxNG9mDgqk3e/3cqe0oqmZ4gylhSMMVHnV1l9KKmo4mVrDmMflhSMMVHn0MwO9OuazOff7/A6lIhjScEYE5UOSEtgSd5uSiusK/hAlhSMMVHpgjEHsnVPKd9s2u11KBHFkoIxJiqlp8QD8Oma7R5HElksKRhjotLoPp0AuG/OavILyzyOJnJYUjDGRKVYv4+DuiYDsHD9Lqqq1eOIIoMlBWNM1Lrz58PI7JTI5GcWcsH0z9lWWOp1SJ6zpGCMiVqH9+lMzvXZ3P3zYXydV8D4f+Tw+CdrvQ7LU5YUjDFRLcbv4xeH9+Kda4/m8L6duW3WMuat3OZ1WJ6xpGCMMUC/rik8ctEoBmak8LvnFvFU7jpUo+88gyUFY4xxJcT6efySwxl4QCp/fX0p89fu9DqksAtpUhCRjiLysoisEJHlIpIlIp1FZLaIrHafOwWUnyoia0RkpYicHMrYjDGmPr06JzFhRA8A1u8o9jia8Av1nsL9wLuqOggYDiwHpgBzVXUAMNd9jYgMBs4HhgCnAA+LiD/E8RljzD6G9uxAXIyPh3PWRN0hpJAlBRFJA44BHgdQ1XJVLQAmADPcYjOAn7nDE4CZqlqmqmuBNcCYUMVnjDENGXlgJ27/2VDW7Shm6eY9XocTVqHcU+gH5ANPiMgiEXlMRJKBDFXdAuA+d3PL9wQ2Bsyf544zxpiw69UpCYBFG6Krh7ZQ9tEcA4wErlHVL0TkftxDRQ2Qesbts98mIpOASQAZGRnk5OQEHVBRUVGzyrcHVufoYHVufXsrlC4Jwm2zlrJo2SrO7BeLSH0/U+ETlu2sqiF5AAcA6wJeHw28BawEurvjugMr3eGpwNSA8u8BWY2tY9SoUdoc8+bNa1b59sDqHB2szqGRt6tYz37oE+190yx94csNIV9fU1qrzsACbeB3NWSHj1R1K7BRRA52Rx0PLAPeAC5xx10CvO4OvwGcLyLxItIXGADMD1V8xhjTlJ4dE/n3xaM5IC2BG19ewmMff+91SCEXysNHANcAz4pIHPA9cBnOeYwXReQKYAMwEUBVl4rIiziJoxK4WlWt9wtjjKe6psYz7/psDvnruyzfUuh1OCEX0qSgqouB0fVMOr6B8tOAaaGMyRhjmisxzs9BXZPZXFDidSghZ3c0G2NMEDokxvLDntJ2f9+CJQVjjAnCRWN78/32vZw3/XM2tOM7nS0pGGNMEM4+rCf/O2EI89fu5KrnFrK7pMLrkELCkoIxxgRBRLh4bG9OOCSDbzftIfsf8yitaH/XwlhSMMaYIIkIj10ymnMO68mu4grW7djrdUitzpKCMcY00w2nHEzHpFjOefgz5q1oXx3yWFIwxphm6t4hkXeuPZq+6clc8/witheVeR1Sq7GkYIwxLdC9QyL3n38YxeWV3DdnldfhtBpLCsYY00L9u6Vw4uAM3lqypd3cv2BJwRhj9sPRA7qyq7iCz79vH113WlIwxpj9MLqP06Pwt5t2exxJ67CkYIwx++HgjFQGd09j9vIfvA6lVVhSMMaY/SAiDOqeynfbirwOpVVYUjDGmP1UXlnN7pIK9pS2/aYvLCkYY8x+OmdkTyqrlXe+2eJ1KPvNkoIxxuynrH7pDOiWwu2zllNS3rbbQ7KkYIwx+ykxzs9vj+tPYVklq7e17d7ZLCkYY0wr6JeeAsAHbbwtJEsKxhjTCg7N7MAxA7ty35zVTP3vkjZ70tmSgjHGtJJ7Jg7jwiMO5Pn5G3lo3hqvw2kRSwrGGNNKuqUmMO3sQzn90O489vFavstve/cuWFIwxphWdsmRfaiqVrYUlHodSrNZUjDGmFbWMSkWgDVt8EokSwrGGNPKBnRLYVTvTjyc812ba1LbkoIxxrQyEeFnh/VkW2EZH67K9zqcZrGkYIwxIXDOYT3p1zWZK2cs4L2lW70OJ2iWFIwxJgSS42N49apxdEuN57VFm7wOJ2iWFIwxJkQ6JMYysncnPlm9nXXb93odTlBCmhREpKOIvCwiK0RkuYhkicgtIrJJRBa7j9MCyk8VkTUislJETg5lbMYYEw43nHwwhWWVTP/4e69DCUqo9xTuB95V1UHAcGC5O/5fqjrCfbwNICKDgfOBIcApwMMi4g9xfMYYE1K9uyRz3KBuPPfFBhZt2OV1OE0KWVIQkTTgGOBxAFUtV9WCRmaZAMxU1TJVXQusAcaEKj5jjAmXG04+GIAv1u70OJKmhXJPoR+QDzwhIotE5DERSXan/VZElojIf0SkkzuuJ7AxYP48d5wxxrRpgw5I5ciDujD9o+8jvqE8CdWNFSIyGvgcGKeqX4jI/cAe4EFgO6DAbUB3Vb1cRB4CclX1GXf+x4G3VfWVOsudBEwCyMjIGDVz5sygYyoqKiIlJWX/K9eGWJ2jg9U58q3bXcUtuaVcOCiOE/vEtmgZrVXn8ePHL1TV0fVOVNWQPIADgHUBr48G3qpTpg/wrTs8FZgaMO09IKuxdYwaNUqbY968ec0q3x5YnaOD1TnyVVZV60FT39Jpby1r8TJaq87AAm3gdzVkh49UdSuwUUQOdkcdDywTke4Bxc4GvnWH3wDOF5F4EekLDADmhyo+Y4wJJ79PyDqoCy8t2Mhn3233OpwGhfrqo2uAZ0VkCTAC+Dtwt4h8444bD/wBQFWXAi8Cy4B3gatVtW13dmqMMQH+duZg0hJj+eWjX/DcFxu8DqdeMaFcuKouBuoet7q4kfLTgGmhjMkYY7zSv1sq7157DFc9u5A/v/YNBx+Qwqjenb0O6yfsjmZjjAmjxDg/D/xyJB2T4rjplW9YvmWP1yH9hCUFY4wJs5T4GG49awhrthUx4aFPKauMnCPllhSMMcYDZw7vwZ9OG0R5ZTUrt0ZOZzyWFIwxxiPnjzkQEZi1ZIvXodSypGCMMR7xi6AK4nUgASwpGGOMR5Li/Aw6IJVnPl/P/AhpF8mSgjHGeEREmHH5GDLSEvjN0wsorfD+hLMlBWOM8VBGWgLXnjCAXcUVEdGfsyUFY4zx2ClDD6B/txR+9/wizw8jWVIwxhiPxcf4efE3WcT5fby+2Nv+nC0pGGNMBOicHMeIAzvy3tKtbNld4lkclhSMMSZC/PWMwewpqWT6R97152xJwRhjIsSAjFQGHpDCJ6u3U1lV7UkMlhSMMSaCXJLVh9XbivjHeys9Wb8lBWOMiSATR/fi9EO78+Rn69hdEv7+nC0pGGNMhLlw7IGUVVZz35xVYV+3JQVjjIkwyXFO/2efrA5/t52WFIwxJsIMy+zAr7J6s3pbUdgvT7WkYIwxEUZEuHhsbwA+XhXevQVLCsYYE4F6d0kmOc7PE5+tC+vlqZYUjDEmAsXF+Ljh5INZvmUPK38IX89slhSMMSZCHTOwKwBL8naHbZ2WFIwxJkJ1To5DBFZs2RO2dcaEbU3GGGOapWNSHBNHZTIjdz3j+qcTF4Z12p6CMcZEsGlnH0qPDglM/+h7VDXk67OkYIwxESzW7+Pq4/qzYP0ulu4IfXedlhSMMSbC/XxkJj07JvLcivKQr8uSgjHGRLiEWD+XjevD5iJlW2FpSNdlScEYY9qAw/t0BuC9b7eGdD0hSwoicrCILA547BGR34tIZxGZLSKr3edOAfNMFZE1IrJSRE4OVWzGGNPWDMvswMBOPu6dvYrtRWUhW0/IkoKqrlTVEao6AhgFFAOvAlOAuao6AJjrvkZEBgPnA0OAU4CHRcQfqviMMaYtEREuHRJPQXEFj+R8F7L1hOvw0fHAd6q6HpgAzHDHzwB+5g5PAGaqapmqrgXWAGPCFJ8xxkS8Hik+DumexmOfrOWfs0PT14KE47pXEfkP8JWqPigiBaraMWDaLlXtJCIPAp+r6jPu+MeBd1T15TrLmgRMAsjIyBg1c+bMoOMoKioiJSVl/yvUhlido4PVOToUFRUh8cm8sLKc7sk+Tu0b26LljB8/fqGqjq5vWsjvaBaROOAsYGpTResZt0/GUtXpwHSA0aNHa3Z2dtCx5OTk0Jzy7YHVOTpYnaNDTZ1PPxFUFZH6fjb3TzgOH52Ks5fwg/v6BxHpDuA+b3PH5wG9AubLBDaHIT5jjGlzQpEQIDxJ4QLg+YDXbwCXuMOXAK8HjD9fROJFpC8wAJgfhviMMca4Qnr4SESSgBOB3wSMvhN4UUSuADYAEwFUdamIvAgsAyqBq1U19Pd0G2OMqRXSpKCqxUCXOuN24FyNVF/5acC0UMZkjDGmYXZHszHGmFqWFIwxxtSypGCMMaaWJQVjjDG1LCkYY4ypZUnBGGNMLUsKxhhjallSMMYYU8uSgjHGmFphaTo7VEQkH1jfjFnSge0hCidSWZ2jg9U5OrRWnXuratf6JrTppNBcIrKgoTbE2yurc3SwOkeHcNTZDh8ZY4ypZUnBGGNMrWhLCtO9DsADVufoYHWODiGvc1SdUzDGGNO4aNtTMMYY04g2mRREpJeIzBOR5SKyVESudce/ICKL3cc6EVncwPx/cOf7VkSeF5EEd/wIEfncnX+BiIwJY7Ua1Qp1vtat71IR+X3A+M4iMltEVrvPncJTo6aFsM7/EJEVIrJERF4VkY5hqVAQQlXngOnXi4iKSHpoaxK8UNZZRK4RkZXutLtDX5vghPCzvf+/Yara5h5Ad2CkO5wKrAIG1ylzL/DXeubtCawFEt3XLwKXusPvA6e6w6cBOV7XtZXqPBT4FkjC6W1vDjDAnXY3MMUdngLc5XVdw1Dnk4AYd/iuaKizO70X8B7OvT3pXtc1DNt5vPs63n3dzeu6hqHO+/0b1ib3FFR1i6p+5Q4XAstxfuwBEBEBfgE838AiYoBEEYnBeWM31ywaSHOHOwSM99x+1vkQ4HNVLVbVSuBD4Gx32gRghjs8A/hZSCrQAqGqs6q+744D+BzIDF0tmieE2xngX8CNOJ/ziBHCOv8PcKeqlrnL3ha6WjRPCOu8/79hXmfMVsi4fYANQFrAuGOABY3Mcy1QBOQDzwaMP8Rd1kZgE85df57XcX/r7NZrFU5/2UlALvCAO62gTtldXtcv1HWuU+5N4CKv6xeG7XwWcL87vI4I2lMIYZ0XA7cCX+D8cB7udf3CUOf9/g2LoQ0TkRTgFeD3qronYNIFNLCX4B4znwD0BQqAl0TkIlV9BuefxR9U9RUR+QXwOHBCCKvQbC2ps6ouF5G7gNk4yfBroLK+spEoVHUWkZvdcc+GIu790Zp1FpEk4Gacw2YRKwTbOQboBIwFDgdeFJF+6v56RoIQ1Hn/f8O8zpL7kV1jcY6PXldnfAzwA5DZwHwTgccDXv8KeNgd3s2Pl+kKsMfrerZGnetZzt+Bq9zhlUB3d7g7sNLreoa6zu7rS3D+YSV5XcdQ1xk4FNiGs4ewzv0B2QAc4HVdQ7mdgXeB7IBp3wFdva5riOu8379hbfKcgnu87XFguar+s87kE4AVqprXwOwbgLEikuQu53ic43ngHH871h0+DljdupG33H7WGRHp5j4fCJzDj/9C3sD5gcR9fr01494foaqziJwC3AScparFoYi9pUJRZ1X9RlW7qWofVe0D5OGc5Nwakko0Uwg/26/hfI8RkYFAHBHSgF4I67z/v2FeZ8sWZtijcE6oLME5brgYOM2d9iQwuU75HsDbAa9vBVbgnMF/mh+vTjgKWIizO/YFMMrrurZinT8Glrl1Oz5gfBdgrvvhmQt09rquYajzGpxjrjXLfMTruoa6znXmWUcEnVMI4XaOA55xv+dfAcd5Xdcw1Hm/f8PsjmZjjDG12uThI2OMMaFhScEYY0wtSwrGGGNqWVIwxhhTy5KCMcaYWpYUjKlDRKrcViaXisjXInKdiLTqd0VEJovIr9zhS0WkR2su35iWsktSjalDRIpUNcUd7gY8B3yqqn8L0fpygOtVdUEolm9Mc9iegjGNUKdlzUnAb8XhF6c/hi/F6Y/hNwAiki0iOSLysjh9NTzr3rWKiNwpIsvc8ve4424Rp2+Dc4HRwLPu3snpIvJqzfpF5EQR+W/4a26iVZtuEM+YcFDV793DR91wGlPcraqHi0g88KmIvO8WPQwYgtPUwKfAOBFZhtOs8SBVVanToY+qviwiv8XdU3ATyb0i0lVV84HLgCfCUU9jwPYUjAmWuM8nAb9ye8T6AqeZkAHutPmqmqeq1TjNFvQB9gClwGMicg7QaFtL6hzPfRq4yE0gWcA7rVkRYxpjewrGNEFE+gFVOC2NCnCNqr5Xp0w2UBYwqgqnd7dKt0vE44Hzgd/iNtLWiCdw+nkoBV7SHzsEMibkLCkY0wgR6Qo8AjzoHv55D/gfEflAVSvc1jc3NTJ/Ck7z3G+LyOc4jfHVVYjTJSMAqrpZRDYDfwZObM36GNMUSwrG7CvRPTwUi9P3wNNATfPGj+EcFvrKPf6fT+NdmKYCr4tIAs5exh/qKfMk8IiIlABZqlqC0/FPV1Vdtr+VMaY57JJUYyKQiDwILFLVx72OxUQXSwrGRBgRWQjsBU5Ut9N5Y8LFkoIxxphadkmqMcaYWpYUjDHG1LKkYIwxppYlBWOMMbUsKRhjjKllScEYY0yt/w+sPFlN4f9ywQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# steps selection\n",
    "plt.figure(figsize=(6,6))\n",
    "var.plot(y='levels',yincrease=False,label='Density profile')\n",
    "rho_g.plot.line('r.',y='levels',yincrease=False,label='Steps points')\n",
    "plt.grid();plt.ylabel('Pressure (dbars)');plt.xlabel('Density')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "26915258-05da-45d3-bb13-a88ff058c599",
   "metadata": {},
   "outputs": [],
   "source": [
    "# statistics : thickness, number of steps, mean thickness\n",
    "nb_steps = xr.DataArray(np.array([rho_g.levels.size//2]),dims='time',coords={'time':var.time.values.reshape(1,)})\n",
    "thickness = rho_g.levels.diff(dim='levels')[::2]\n",
    "mean_thick = xr.DataArray(np.array([thickness.mean()]),dims='time',coords={'time':var.time.values.reshape(1,)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec305a90-4397-4147-b50a-8fbd8865aa14",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}