{ "cells": [ { "cell_type": "code", "execution_count": 3, "id": "45c9fdb9-a8e4-487e-9207-e05a6c0b909d", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import xarray as xr\n", "import cmocean.cm as cmo\n", "import numpy as np\n", "import my_functions as my\n", "import gsw\n", "from matplotlib.dates import date2num\n", "import matplotlib.dates as mdates\n", "import my_plot_params" ] }, { "cell_type": "code", "execution_count": 4, "id": "4366ef40-295a-4f98-b8fc-fa74218ca5bf", "metadata": {}, "outputs": [], "source": [ "dat_saz = xr.open_dataset('../data/dat_saz_6H.nc')\n", "dat_pfz = xr.open_dataset('../data/dat_pfz_6H.nc')\n", "dat_miz = xr.open_dataset('../data/dat_miz_6H.nc')" ] }, { "cell_type": "code", "execution_count": 219, "id": "602a6967-53da-4a6f-905c-b943a45c8d2a", "metadata": {}, "outputs": [], "source": [ "mld_05 = my.calc_mld(dat_saz.density.values, dat_saz.depth.values, den_lim=0.03)" ] }, { "cell_type": "code", "execution_count": 5, "id": "9cf0fd6f-ebc3-4696-98e1-9ee1cb1e627b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "good: 65.04065040650407\n", "ok: 22.22222222222222\n", "bad: 12.737127371273713\n", "good: 86.58892128279884\n", "ok: 13.119533527696793\n", "bad: 0.2915451895043732\n", "good: 85.96059113300493\n", "ok: 11.576354679802956\n", "bad: 2.4630541871921183\n" ] } ], "source": [ "for i, data in enumerate([dat_saz, dat_pfz, dat_miz]):\n", " \n", " data['mld_03'] = (('time'), my.calc_mld(data.density.values, data.depth.values, den_lim=0.03))\n", "\n", " mld = data['mld_03'].astype(int).values\n", " mld_15 = (data['mld_03'].values*1.5).astype(int)\n", " \n", " QI = np.ndarray(mld.shape)\n", " \n", " for i, val in enumerate(mld):\n", " \n", " den_mld = data['density'].isel(time=i)[:val].values\n", " den_mld_mean = data['density'].isel(time=i)[:val].mean().values\n", " rmse_mld = np.std(den_mld[den_mld>0]-den_mld_mean)\n", " \n", " val_15 = mld_15[i]\n", " den_mld_15 = data['density'].isel(time=i)[:val_15].values\n", " den_mld_mean_15 = data['density'].isel(time=i)[:val_15].mean().values\n", " rmse_mld_15 = np.std(den_mld_15[den_mld_15>0]-den_mld_mean_15)\n", " \n", " QI[i] = 1 - (rmse_mld/rmse_mld_15) \n", " \n", " good = QI>0.8\n", " ok = ((QI<0.8) & (QI>0.5))\n", " bad = QI<0.5\n", " \n", " good = data.sel(time=good)\n", " ok = data.sel(time=ok)\n", " bad = data.sel(time=bad)\n", " \n", " print('good: '+str(100*len(good.time)/len(data.time)))\n", " print('ok: ' +str(100*len(ok.time) /len(data.time)))\n", " print('bad: ' +str(100*len(bad.time) /len(data.time)))" ] }, { "cell_type": "code", "execution_count": 10, "id": "7d4049a1-08db-4dd1-8afc-c786ae0bb09d", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvIAAAG8CAYAAAC8IUHNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABoIElEQVR4nO3dd3gc1dn38e+9q96t5t4rtsEUA8YQeiAQAiSBFAglBUjPm04CT5InPU96TwhJCCEEAikQSmihN2OqDe69W7JsS5asft4/ZlYWiyxL8q5mZvX7XNde8s62e7x7du49c59zzDmHiIiIiIhESyzoAEREREREpP+UyIuIiIiIRJASeRERERGRCFIiLyIiIiISQUrkRUREREQiSIm8iIiIiEgERTqRN7MrzGyFme01s6fN7LigYxKRnqm9ikSH2qtINEQ2kTezS4HfADcB7wR2AfeZ2cQg4xKRN1J7FYkOtVeR6LAoLghlZgasAe51zn3E35YNLAPucs59Msj4RGQftVeR6FB7FYmWqPbITwHGA3cmNjjn2oC7gbcEFZSI9EjtVSQ61F5FIiSqifw0/+/KpO2rgclmFh/keERk/9ReRaJD7VUkQqKayJf4fxuStjfg7VPh4IYjIr1QexWJDrVXkQjJCjqAATL/b3KBf2J7Z9cGsyuBK5OfoLCw8KgZM2akJzqRCHr++edrnXNVaXjqPrdXUJsV6as0tVm1V5E0SNcxNqqJ/G7/bzGwrdv2IrwvmcbEBufcdcB1yU8wd+5ct3DhwnTGKBIpZrYuTU/d5/YKarMifZWmNqv2KpIG6TrGRrW0ZoX/d1LS9knAMhfFqXhEMpfaq0h0qL2KREiUE/kNwPmJDf70WG8FHgooJhHpmdqrSHSovYpESCRLa5xzzsy+C/zCzHYCTwIfByqBHwcanIi8jtqrSHSovYpESyQTeQDn3K/MLB/4FPBp4CXgTOfc6kADE5E3UHsViQ61V5HoiGwiD+Cc+yHww6DjEJEDU3sViQ61V5FoiGqNvIiIiIjIkKZEXkREREQkgpTIi4iIiIhEkBJ5EREREZEIUiIvIiIiIhJBSuRFRERERCJIibyIiIiISAQpkU+RptZ2fvPoKm59bn3QoYhIH6zf0cTvHlvN/a9uDToUEelFXWMrv39iDbct3BB0KCKho0Q+RXY2tfHde5fy0wdXBB2KiPTBsm0NfOueJdz6nJIDkTCr3dPCN+56jese08KyIsmUyKfI3tYOAPJy4gFHIiJ90dHZCUA8ZgFHIiK9yYl7qUpLe2fAkYiEjxL5FEkk8vnZSuRFoqC90wGQHdfXoEiY5WZ7bbRVibzIG+gIliJ725TIi0RJh5/Iq0deJNxys7zjauI4KyL7KJFPka5EXqU1IpHQ3uEl8llK5EVCrSw/m5ysGLv3ttHQ3BZ0OCKhokQ+Rbpq5NUjLxIJ6pEXiYZYzJhYUQjA2tqmgKMRCRcl8inS3KZEXiRK2vzBrllxJfIiYTepykvkV9XsCTgSkXBRIp8irR1eUpCjgXMikaAeeZHomD26FIAFa+sCjkQkXJR1pkhnV1IQcCAi0if7auTVaEXC7vgplQA8saI24EhEwkVHsBTpcOrdE4mSRI+8BruKhN+ho0spyctifV0T63Y0Bh2OSGgokU+RRI98zJQUiERBYh75uGrkRUIvHjNOml4NwJ0vbQ44GpHwUCKfIqq3FYmWxMqu6pEXiYYLjhoDwN+e39DVeSYy1CmRTxG/3FY98iIR0daR+PGtr0GRKDhhSiWjSvPYULeXp1fvCDockVDQESxFnF8jrzxeJBq6zqKp0YpEQjxmvOvosQD8+pFVAUcjEg5K5FMk0RPvdLZPJBISNfKaR14kOi6fP4Hi3CyeWFnL06vUKy+iRD5FEmW2ncrkRSIh0VZVIy8SHWUFOVxx4iQAfnj/sq6z4SJDVVbQAWSKuD+BfLsG4IhEQnuHBqiLRNH7j5/AH59cw8J1O3l0eQ0n+7PZSAZ48GTY/ui+629+Eqrm77u+5iZ4+pJ91ydeBsfdsO9xFcfCmc/0/hp3TIDGda/fFsuFvCqoPglmXg1lsw9uPwaReuRTpCA7DsDe1o6AIxGRvtCsNSLRVJyXzYdPmgzAD+9frl75TNY9qQfY/lh6XqezBZo2wtq/wH3HwNYH0/M6aaBEPkUKc71EvrGlPeBIRKQv9s0jr69Bkai59LgJVBXnsmjTbu5/bVvQ4Ui6JCfyNSlM5EefC+dvgPPWwlkvwpzvQDwPOvbCU++D1t2pe6000hEsRQpzvSqlPUrkRSIhUSOvDnmR6MnPifPxU6YA8LOHVqhXPtMUeLMTUfMkdPp5VfN2qF/m3z7u4F8jng8FY6BwPAw7HGZdDUf+yH+tbbD6hoN/jUGgRD5FygtzANixpzXgSESkL7IT41o6lACIRNG7jx5LRWEOr26uZ8GauqDDkVSqnAeWBe17oO4Fb1uirCZ/FBRNTM/rTrwcYtnev7fcm57XSDEl8ilSVZwLQM2eloAjEZG+yM3yvv5a2jWuRSSK8rLjXDxvPAC/f2JNwNFISsULvF5y2Fdek0jkq05I3+tm5UPxVO/fuxal73VSSIl8ilQU5hIzqGtspbW9M+hwROQAcrO8cS0tbWqvIlF1ybzx5MRjPLBkGxt3NgUdjqRS1fHe3zck8sen93WzS72/rTvT+zopokQ+ReIxY0RJHgCbd+0NOBoROZB9PfJK5EWiqqo4lzNmDcc5+PfLW4IOR1IpkbDXPAEtdbDb7yFPZ488AP7AKReNs7VK5FNoXEUBAOvq1CsgEna52SqtEckE5x8+GoA7XtoUcCSSUpV+It+2G5b/ElwnZBVB2Zz0vm5bvfc3pyy9r5MiSuRTaEJFIQDrdjQGHImIHEiev/ZDs0prRCLtxGlVlOZns3RrA8u3NQQdjqRKwSgo9Ae1LvuJ97dyHsTi6XvNznbYs9L7d+mh6XudFFIin0KTq4oA9EUiEgFF/pSx9c1tAUciIgcjJyvG6YcMB+Cx5TUBRyMplSivafVnJapMc338xjugo9n79+hz0vtaKZIVdACZ5JCRJQAs2aJEXiTsElPG1jVqyliRqDthagV/f2EjT6ys5UNvmhR0OJIqVcfD2pv2Xa8+QH18227Y/J+ebys/EvKq913v2Out5grQ3gS1T8OLn/eu51XD5A8MPO5BpEQ+hWaMLAZg6ZZ6OjsdMa00k1kePPn1q8y9+Umomr/v+pqb4OlL9l2feBkcd8O+x1UcC2c+0/tr3DEBGte9flssF/KqoPokmHk1lM0+uP0QwJtpCrT2g0gmmD+5EoAFa+pobe8kJ0sFBxmh+ww1FoeKeb3fv34pPHJWz7e96Z8w9vx91zfd6V2SxfPhuD9Ddkm/ww2CPukpVFmUy8jSPBpbO1hZsyfocCTdkpeO3p7CpaO762zxeg3W/gXuOwa2Ppie1xlihhV6i36oR14k+oaX5DG+ooCm1g7W1GqcWsYonQXZZd6/y+ZAdlF6XieW460mO/FSeMvzMPKM9LxOGqhHPsWOnlDOnS9v5tk1dUwbXhx0OJJO2x+FWV/ad70mhYn86HPh6F9601+17vROFS7+X+9U4FPvg3OWQU5p6l5vCEr0yNc1teKcw0xn0ESibGp1Eet2NLGqZg/TR+j4G0mnP/L66xaDC3uYzz35fsnXe3Pe2v7FFHLqkU+xoyeWA/CclovOXAVjvb81T3oj3AGat0P9Mv/2cQf/GvF8KBgDheO91e1mXQ1H/sh/rW2w+oaDf40hLj8nTn52nNb2ThpbNQVlxnjwZLjZ9l1qnnr97Wtuev3tT1/++sfdd4BT9+CVwHV/jpsNbsmDf431fmjvWpzafZI+mVzt9dau3K4z4jJ0KJFPsWMm+In82jqccwFHI2lROQ8sC9r3QN0L3rZEWU3+KCiamJ7XnXg5xLxyELbcm57XGGKqir1e+a27mwOORNJGJXBDRmLmOCXyMpQokU+xqdVFlOZns2V3Mxt3aoXXjBQv8HrJoYelo9O44lxWPhRP9f69a1H6XmcIGe8v4ra+TjW1GSs5kU91Cdz5G7xT9We9CHO+A/G8fSVwrbtT91pyQFPUIy9DkBL5FIvFjKMnDAO80fOSoRIj6d+QyKd5jttsvy6+tYeaQem3RCK/tlarMWcclcANOYke+dW1e+js1BlxGRqUyKfBMX6dvBL5DJZI2GuegJY62O33kKezRx4Af0CmU013KiRWY9YsUxlIJXBDTml+NlXFuTS3dbJpl86Iy9CgRD4NjplYAXh18pKhEqvLte2G5b8E1wlZRd70WOnUVu/9zSlL7+sMEbNHe2c4Xtm4K9hAJPVUAjckTfdni1u6VQszytCgRD4NZo0qoSAnzuraRrbXaxBdRioYBYV+j96yn3h/K+dBLJ6+1+xshz0rvX+XHpq+1xlCZo8uxQyWbmmguU1nOTKOSuCGnMS0k8u21gccicjgUCKfBtnxGEeN9+vk1SufuRLJQKv/HlemOTnYeAd0+D8MR5+T3tcaIopys5haXUR7p2PRJg1MzDgqgRtyEom8euRlqFAinybHqk4+8yX36lUfIDlo2+0t7NTTpXn76+/bsdebyq5pI9Qvh9V/guc+4t2WVw2TP5C6/RjijpvklcI9vqI24Egk5VQCN+TM6OqRVyIvQ4NWdk2TRJ38s6uVyGes7om8xaHiAAvJ1C+FR87q+bY3/RPGnr/v+qY7vUuyeD4c92fILul3uNKzN02t4k9Pr+PxFTV85s3Tgg5HUilRAte4RiVwQ8TU6mLMYHVtIy3tHeRmpfG9FgkBJfJpctiYUnKyYizb1sDOxlaGFeYEHZKkWuksyC6Dtl1eD192UXpeJ5YDecNh+Ckw82ooPSQ9rzNEHTe5guy48fKGXdQ1tlKutppZqo73EnmVwIXPzZbyp8wHxmdfx9rWUaz801Rm5a9J+WuIvMFFwU13qkQ+TfKy4xwxtoxn19Tx3No6zpg1IuiQ5GCd/sjrr1sMLuxhMFvy/ZKv9+a8tf2LSQ5aYW4W8ydX8ujyGu5dvIWLjx0fdEiSSlXHw9qb9l3vawlcT8qP9ErbEhIlcADtTVD7NLz4ee+6SuACMzN/NWtbR/Hq3slK5CXjhTKRN7M48CngCmAcsA74FfBL55wzMwO+DFwFVAJPAp9wzi0NKOQeHTuxnGfX1LFgjRJ5yVyZ0F7fNmcUjy6v4c6XNiuRzzQqgXudULXXNPViHvrIKu75z1IWj/0x7zpvdlpeQyQswjrY9X+AbwM3AecCfwN+AvhdHXwFuBb4AfAeoBR4yMxKBz3SXnTVyWvAq2S2yLfXM2YNJycrxoK1dWzZrYVkMkqiBA7SXwJXMBYmXgpveR5GnpGe1zl4kW+vB3Kovz6EZqKSoSB0PfJmFgM+A3zfOfctf/NDZlYFfM7Mfg18Dviac+5n/mMex+tV+CDwowDC7tGR48vIihmvbt5NQ3MbxXnZQYckklKZ0l5L8rI5/ZBq7lm0lVuf28D/O12DXiNLJXD7lSnt9UBmj/bOhCzZUk97RydZ8bD2WYocvDB+ukuBG4F/JG1fBlQBpwJFQNf5TOfcTuBR4C2DFGOfFORkceiYUjodPL9OC4NIRsqY9pooqbllwQbaOzoDjkYkLTKmvfamrCCHseX5NLd1srJmT9DhiKRV6BJ559xO59zHnXMvJt30NmAjMMa/virp9tVA6LrRjvHnk1d5jWSiTGqv8ydXMKmykK31zTy0dPuBHyASMZnUXg+kq7xmo8prJLOFLpHviZl9CDgd+D+gBGhxzrUm3a3Bvy35sVea2cLkS01NTfoDRwtDydBzMO3Vf3wgbdbMuOjYcQDc9My6tL6WSFhEtb0eyGw/kV+sOnnJcKFP5M3sYuA3wO3AL/DWvu5pqLsBbzgf7py7zjk3N/lSVVWV1rgT5k4oxwxe2biLva1arlsy28G2Vwi2zV5w1Bhys2I8vqKWtbWNaX89kSBFvb325rDRZQC8okReMlyoE3kz+zTwZ+Au4GLnnAN2A7lmljxytMi/LVRK8rKZObKEtg7HixtUJy+ZKxPaa1lBDuccNgqAmxesDzgakfTJhPbam+QBryKZKrSJvJl9G2+E/J+BC7qd6luB1zswMekhk/AG7ITO0RO88poX1+8KNhCRNMmk9poor/nXi5vo7AxutT6RdMmk9ro/GvAqQ0UoE3kz+xTwJeCnwOXOufZuNz8FNAPnd7v/MOAk4KFBDLPPZo7yegZe21IfcCQiqZdp7fXIcWWMGZbP9oYWnl+vs2iSWTKtvfZGA15lKAjjPPIjge8Bi4BbgGO9hea6LAR+DnzTzDqB5cA1QD1wfUqDudkOfJ8+mNk0GfgpS5YtgJuPSslzivRLmlZQDFV7TREz4+xDR3LdY6u5Z9GWrjNqIlGXie21N7NHl3LPoq0s3rSbC+eODTockbQIXSIPnAnkAocCT/dwexXe8tGdeAtXFOH1IlzmnAvlz+4peevJop21LaPY25lLfqwl6JBEUiXj2ivAWbNHcN1jq3ngtW189W2zgg5HJFUysr3uz+xRXo/80q0NAUcikj6hS+SdczcAN/Thrlf7l/RJUS9mHjD5x4+xbFsDy07ayuFjy1LyvCJBC1V7TaHDxpRRmp/Nxp172VDXxNjygqBDEjlomdpe92dUWR4A2xvUeSaZK5Q18pnokJHFALy2WXXyImEXjxnzJnklNU+tqg04GhEZiKpiP5Gvbw44EpH0CV2PfKY6ZGQJ/3ppM8u2KpEXSZsUjWsBOG7nOdzHh3nmoZ/y7hU/StnzivRJmsa1DCUleVnkZsVobO2gqbWdghylPJJ51CM/SCZWFgKwrq4p4EhEpC/mFCwH4LXmSQFHIiIDYWYU53nJ+56W9gPcWySa9PN0kCRqbDcokRdJnxT2Yk5vbce+eh+rWifS8q52crPiKXtuERkcedleu21u1aJQkpnUIz9IEon8xp178RbQE5EwK8jJYmJFIe2djpXbtaCMSBQV5HiJfFObeuQlMymRHyRFuVkMK8impb2TGo2gF4mECX5J3MadewOOREQGItFvFrPUjZ8RCRMl8oMo0Su/XuU1IpEwotSb9WLrbs16IRJFbR1eSU12XOmOZCZ9sgfR2GF+nfxOJfIiUTDKT+Q371aPvEgUtbYnEnn1yEtmUiI/iMYMywdg8y717olEQbU/D3VtQ2vAkYjIQLR2eLU1OVlKdyQz6ZM9iBKn6beod08kEsoLcwCoa9S4FpEoSpTW5Ki0RjKUPtmDaKTqbUUiZVgikW9qCzgSERmIfaU1SnckM+mTPYhGlqq0RiRK1CMvEm2JHvks1chLhlIiP4i6euTrlciLREEikd/ZqB55kahxztHe6dXIZ8eU7khm0id7EFUU5ZIVM+oaW2lu6wg6HBE5gJK8LLJixp6Wdlra1WZFoqTDT+JjBrGYeuQlMymRH0TxmDG8xOuV36ZeeZHQM7OuOnn1yotES6I3Pkv18ZLB9OkeZInyGtXJi0RDeYGXyO9QnbxIpHQl8uqNlwymRH6Qda0UWa8pKEWiQHXyItHUnhjoqkReMpgS+UE2qsybuWaLpqAUiYSumWuatCiUSJR0DXRVaY1kMH26B9mIEs0lLxIlXYn8HpXWiERJu7+qa1w98pLBlMgPMtXIi0SLFoUSiab2Ti0GJZlPn+5Bphp5kWip0KJQIpGkHnkZCpTID7JEjbxKa0SiQdNPikST8/8qj5dMpkR+kFUW5RKPGbV7WrXAjEgEJHrkNf2kiIiEjRL5QRaPGcOLcwH1yotEwbDEPPJ7NGuNSBS5A99FJLKUyAdgTHkBABvqVCcvEnaj/XK4Tbv24pxSApGoyM3yUpyWts6AIxFJHyXyAZhQ4SXya3Y0BhyJiBxIaUE2wwqyaWrtYHuDymtEoqIwNwuAPS3tAUcikj5K5AMwobIQgHW1SuRFoiDRZlfXqM2KREWRn8g3trbT2amzaZKZlMgHYEKFlxSs3dEUcCQi0heHjCwB4IX1OwOORET6Kh4zKgpzcA62NWhMmmQmJfIB2JfIq3dPJAqOn1wJwJMrawOORET6Y6J/Nm2NzqZJhlIiH4Dxfo38+h1NdOh0n0joHTe5AjNYuHYne1s1baxIVEyq8hL5VSpllQylRD4AhblZjCjJo7Wjk/V1Kq8RCbvywhwOG1NGa0cn97+2NehwRKSPpg0vBmDxxt0BRyKSHkrkAzJzlFdz+9rm+oAjEZG+uODI0QD8beGGgCMRkb46ZmI5AM+u2RFwJCLpoUQ+IDP9wXOvbVEvgUgUnHv4aHKzYjy5cgcbdCZNJBJmjiyhKDeLtTua2LJba7dI5lEiH5BZfo/8q+qRF4mE0vxszj50JADXP7464GhEpC+y4rGuXvmHl9YEHI1I6imRD4hKa0Si56qTJgHw1wUb1LsnEhFnzR4BwL9f3hxwJCKpp0Q+IGOHFVCcm8X2hhZqtFqkSCTMGFHCWw8bSWtHJ796eFXQ4YhIH5wxawQ58RjPrNnB9nrNJy+ZRYl8QGIx4xC/V37xJtXJi0TF/zttKmbw1wXrWb6tIehwROQASvOzOWl6Fc7B7S9sDDockZRSIh+gI8cNA2DhurqAIxGRvpo6vJj3HjOO9k7HNf9cpKXfRSLg4mPHAXDjU+to6+gMOBqR1FEiH6CjJ3iJ/HNrtOy7SJR88cwZVBbl8Nzandz2vKajFAm7E6dWMbmqkK31zdyzaEvQ4YikjBL5AM0d742kf2njLlratVqkSFSUFmTzP+fMBOBbdy9R3a1IyMVixgdOmAjA9Y+vwTmdSZPMoEQ+QKUF2UwfXkxre6fq5EUi5tw5ozhpWhX1ze38zx2LlRiIhNw7jxxDRWEOizbt5qlVWiBKMoMS+YDNTZTXrFV5jUiUmBnffsehFOVmcd+r27hn0dagQxKRXuRlx7t65X/1yMqAoxFJDSXyATt6gldes3CtBryKRM3osny+dPYMAL5yx2LqGlsDjkhEevO+eeMpys3iyZU7eHnDrqDDETloSuQD1r1HXrNfiETPe48ex3GTKtjR2Mr//vvVoMMRkV6U5md3zWDzm0e1FoREnxL5gI0uy2dkaR6797axsmZP0OGISD/FYsZ333ko+dlx7nhpMw+8ti3okESkFx88YSI58Rj/eXUrK7fruCvRpkQ+YGbWVV6zYI3Ka0SiaHxFIZ8/czoA1/xzEfXNbQFHJCL7U12SxzuPGoNzcN1j6pWXaFMiHwKJ+eSVyItE12XzJ3DkuDK2N7Tw4weWBx2OiPTiqhMnETP454ub2LJ7b9DhiAyYEvkQOHZSBQDPrtmhKexEIioeM75x/mxiBn96ai2vba4POiQR2Y8JlYWcfehI2joc1z++JuhwRAZMiXwITK0uorwwh231Lazb0RR0OCIyQLNGlXLpcRPodN4sNhrALhJeHz5pMgB/XbCenZpxSiJKiXwImBnHTvTq5J9ZrUUqRKLsM2dMo7Iol4XrdvL3FzYGHY6I7Mfs0aWcOK2KptYObl6wPuhwRAYk1Im8meWa2RIzu6HbNjOza8xsvZk1mdkDZjYjwDBTIpHIP6s6eYmwodRm96ckL5tr3urt3nfvXcruJg18lXBSe4UP+QtE3fj0WlrbOwOORqT/Qp3IA18Fkr9AvgJcC/wAeA9QCjxkZqWDHFtKddXJr1advETakGmzvTn/8NEcM7GcHY2t/OD+ZUGHI7I/Q769vmlqJVOri9hW38I9i7YEHY5Iv4U2kTezI4BPArXdthUDnwO+5pz7mXPuTuBMoBj4YCCBpsj04cWUFWSzeXczG+o0gl6iZ6i12d6YGd84bzbxmPGXZ9exdKsGvkq4qL16zIwP+L3yf3hyjTrSJHJCmcibWRbwB+D7wKZuN80DioA7ExucczuBR4G3DGaMqRaLGcf488k/s0Z18hItQ7HNHsj0EcVcMm88nQ6+dfcSJQgSGmqvr/f2I0YzrCCbVzbuZuG6nUGHI9IvoUzkgS8COcB3krZP8/8mr+CwutttkbWvvEZ18hI5Q7LNHsinTptKSV4Wj6+o5ZHlNUGHI5Kg9tpNXnaci48dD8DvNRWlREzoEnl/UM01wIecc8nzQZUALT1sb/Bvi7R5k7we+cdX1NChaeskIoZymz2QYYU5fPK0qQB849+v0dLeEXBEMtSpvfbs0uPGkx037n9tKyu37wk6HJE+C1Uib2Yx4PfA751zT/d0F6CnDNeAHoebm9mVZrYw+VJTE77esZkjSxhXXsD2hhaeWlV74AeIBGyot9m+uPS4CUyuKmR1bSO/e2x10OHIEKb2un/VJXlcOHcsnQ5+8qBWZpboCFUiD3wCGA98xcyy/Do+8GbEygJ2A7lmlp30uCL/tjdwzl3nnJubfKmqqkrbTgyUmfGOI0cD8M8XNh3g3iKhMKTbbF/kZMX4xnmzAfj5f1eyoU6Lvklg1F578fFTppATj3HXK1tYskUD1CUawpbIvx0YDdQBbf5lDnBpt+sGTEx63CQgI+Z4e/sRXiJ/z+ItbN3dHHA0Igc05NtsX8yfUsl5h4+ipb2TL/1jkVZ8laCovfZiVFk+Fx07DoBv3v2aBqhLJIQtkb8KODrpshy4y//3LUAzcH7iAWY2DDgJeGiQY02L8RWFnDlrOM1tnXzn3iVBhyNyIEO+zfbVtW+dSUVhDk+srOV3j6vERgKh9noAHz91CsMKsnly5Q5uelarvUr4ZR34LoPHOfeGX/xmthfY4Zxb6F//OfBNM+vE+wK6BqgHrh/MWNPp2rfO5JFlNdzx0mYuOmZc12w2ImGjNtt3VcW5fP/Cw/jADQv5/n3LmDepgjljy4IOS4YQtdcDqyzK5VtvP5SP/uUFvn33Ek6cWsn4isKgwxLZr7D1yPfFl4Ef4S1acTNe3d7pzrke6/eiaGx5AR85eTIAX/7nIva0tAcckchByfg221enzhjO5fMn0N7p+MRfX2RnY/LkICKBG/Lt9exDR3Le4aPY29bBJ295ieY2zTYl4RX6RN45d7hz7vJu19udc1c750Y454qcc2c455YGGGJafPikyUypLmJVTSOfvvUl1dRKZAzVNttXV581g9mjS1hf18SHb3qe1vYeJwMRGRRqrz3733NnMbosn5c37OLzt7+ienkJrdAn8kNVXnac3106l5K8LB54bRs/1nRYIhkh0bari3N5dk0dX71zsZIEkZApK8jh95fPpSg3i3+/vJkfP7gi6JBEeqREPsQmVhbyi4uOJGbetHV/W7gh6JBEJAVGlubzu0vnkpsV468LNvBbzS8vEjozRpTw84uOIGbws4dWcOtzGvwq4aNEPuROnFbF/5wzE4Av3P4KNz2zLuCIRCQV5owt44fvmgPAd+9dqrYtEkKnTK/mq2+bBcAX/76IWxYomZdwUSIfAe8/fiJfPnsGANf+azHXa+o6kYxwzmGj+Mb53mJR/3PHYv754saAIxKRZJfNn9B1DL76H4v4y7P60S3hoUQ+Iq48cTJfP8/rFfjm3Uv40f3LNABWJANcMm88XzprBs7B5257hbte2Rx0SCKS5MoTJ3PtWw8B4Jp/LuYPT6wJOCIRjxL5CLn0uAn83zsPwwx+9t+VXHHjQnbvbQs6LBE5SFedNJlPnjqFDn9ayr/q9L1I6HzoTZO6Sl2/ftdrfO3OV+lQh5oETIl8xLzr6LH88fKjKc3P5qGl2zn3F0+wZEt90GGJyEH69Jun8dk3T8M5+NI/FvHrR1YFHZKIJPngCRP50bvmkB03bnhqLVfcuFBrvUiglMhH0MnTq7nrEycwa1QJ63Y08fZfPcm/XtwUdFgichDMjE+cNpVvnD8bM/jef5bynXuWaGpKkZB5x5FjuOmDx1JWkM1/l27ngl8/xaZde4MOS4YoJfIRNba8gL9/ZD4XHDWG5rZO/t+tL/G1O1/V4jIiEXfJvPH85N2HkxUzfvvYar7491do71C7FgmTYydV8K+PHs+kykKWbm3gvF88yfPrdgYdlgxBSuQjLC87zvcvOIxvvX1212m+9/7uGbbVNwcdmogchPMOH831l80lLzvG3xZu5GM3v6Bl4kVCZkJlIf/46HyOn1JB7Z4W3nvdM/z9ec08JYNLiXzEmRkXHzuev111HCNL83h+3U7O+fkTLFhTF3RoInIQTp5ezV8+dCwleVnc9+o2PnDDc6rFFQmZsoIcbnj/MVx63HhaOzr57G0v8517l2gQrAwaJfIZ4ohxw/j3J07guEkV1DS0cPH1z2hOapGIO2p8ObdedRxVxbk8tWoHF//uGXY2tgYdloh0kx2P8fXzZvON82cTjxm/fXQ1V964kIZmzSon6adEPoNUFuXy5w8ewweOn0hbh+PTt77MTx5crsFyIhF2yMgSbv/wcYwtz+fljbt593VPs13lcyKhc8m88fz5g8dQVuDNKvfOXz/F+h1NQYclGU6JfIbJisf4yttm8r/nziJm8JMHV/DZ217WIFiRCBtfUchtV81nanURy7ft4YLfPM2GOiUIImEzf3Ild3zseKb4bfW8Xz7BM6t3BB2WZDAl8hnqsvkTuO6SueRnx/nHC5u44saFGiwnEmEjSvO49arjOHR0KevrmrjgN0+xYltD0GGJSJLxFd4g2FOmV7GzqY1Lfv8sd7ykKaIlPZTIZ7DTZw7nb1cdR3lhDo8ur+FDf1rI3lYl8yJRVV6Yw81XHMsxE8vZVt/Cu377NIs27g46LBFJUpKXzfWXHd1V6vqpW17id4+tDjosyUBK5DPcoWNKueXKeVQW5fLEylref8MCGjXzhUhkFedlc+MHjunq7Xvv757hWZ26FwmdeMz4yttmcs3ZhwDwrXuW8I27XqNTM9pICimRHwKmDS/mlivnUV2cyzOr6/jADc+pzEYkwvKy4/z2krmcc9hI9rS0c+kfFvDY8pqgwxKRHlxx4iR++p7DyY4bv39iDZ+69SVa2nUMltRQIj9ETKku4tarjqO6OJdn19Tx2b+9rF4BkQjLyYrx0/ccwXuOHktLeycfunEhjyzbHnRYItKD8w4fzQ3vP4ai3Cz+/fJm3v/H5zQ9paSEEvkhZGJlYdcXyd2LtvCNu1/T1JQiERaPGd9++6G8b944Wts7ufLG5/nv0m1BhyUiPTh+SiW3XjWva12Id/32GU0lKwdNifwQM3NUCb+95Ciy48Yfn1zLDU+tDTokETkIsZjxjfNmc/n8CbR2dHLVn5/ngdeUzIuE0axRpfzjI/OZVFnIki31vP1XT7G2tjHosCTClMgPQcdPqeQHF84B4Ft3L+HF9TsDjkhEDoaZ8dW3zeyaIeMjNz3PfxZvDTosEenB2PICbv/IfA4fW8amXXu58LdPs2yrppKVgelzIm9mo83sfWZ2rZl938yuNrOLzGx4OgOU9Djv8NG8//gJtHc6Pn7zi+xq0rLvIlFmZvzPOYdw5YmTaO90fOzmF7j7lS1BhyUiPSgvzOEvHzqW+ZMrqGlo4d3XaSpZGZgDJvJmdp6ZPQ6sB24EPg28B7gWuAnYZGaPm9nb0hqppNyXzjqEOWNK2bRrL5+77WXVy4tEnJnxpbNm8JGTJ9PR6fjkLS/y75c3Bx2WiPSgMDeLP1x+NKfOqGZXUxsX/e4Znl+nM+TSP/tN5M1sqpk9BvwSeAE4BShxzlU458Y654qAYcDbgOeA35vZ02Y2bTACl4OXkxXjFxcdSUleFg8u2c7vn1gTdEgicpDMjC+cOZ1PnDqFjk7H/7v1Je57VWU2ImGUlx3nN+87irceNpKGlnYu/8MCFm9Sz7z0XW898ncB1wPjnXOfcs495pzb0/0Ozrndzrl7nXOfAUYDvwfuTl+4kmpjywv4vl8v/917l/KC6uVFIs/M+Mybp3X1zH/85hc0NaVISOVkxfjpuw/n7ENH0NDSziW/f1Y189JnvSXyc5xzNzrn+rRqgXOuzTl3PXBoakKTwXLmrBF88ISJXr38X15gZ6Pq5UWiLtEz//7jJ9DW4bjqz8/z1KraoMMSkR5kxWP85N1HcOqManY2tXHx9c9qNhvpk/0m8s65AU1uOtDHSbC++JYZzBlbxubdzXz0Ly/Q1tEZdEgicpDMjK+cM5OLjh1Hiz/PvE7bi4RTTlaMX118JPMnV1C7p4UP/uk56rVolBxAn2atMbNKM/u1mb1gZqt7uKxKd6CSXjlZMX598ZFUFefy9OodfOWOxRr8KpIBzIxvnjebc+eMYk9LO5f/cQHrdqinTySM8rLjXHfpXKYPL2ZVTSOfvuUlrcIuverr9JN/BC4HFgN39HC5Mx3ByeAaVZbP7y6dS25WjL8u2MAfnlwbdEgikgKxmPGDC+dwwpRKave0cukfFlDT0BJ0WCLSg6LcLH536VzKCrJ5aOl2fvjAsqBDkhDrayJ/MvAJ59ylzrlP93RJY4wyiA4fW9Y1+PWbd7+mqetEMkROVozfXHIUh44uZd2OJj5040Ka2/o0BEpEBtm4igJ+edGRxGPGLx9exZMrNb5FetbXRH4ToE/REHHunFF8/szpOAefvvUlzXYhkiGK/HmrxwzL5+UNu/js317WaXuRkDp+SiX/77SpAHzh9ldoUL289KCvifyXgG+a2fFmlpfOgCQcPnryZK54kzeTzYdvep7n1tYFHZKIpEBVcS5/uPxoinOzuHvRFp22Fwmxj5w8mUNHews3fvuepUGHIyHU10R+CZAPPAY0mllH8iV9IUoQzIwvn30I7547lua2Tj7wx+d4fp2SeZFMMG14Mb+8eN9p+zte2hR0SCLSg6x4jB9cOIeceIy/LljPQnWqSZK+JvJ/wkvkvwd8Zj8XyTBmxrffcShvmzPKX6RiAc+u3hF0WCKSAidOq+Jrb5sJwBf//gqvbta0lCJhNH1EMVedNAmA//vPMs0oJ6/T10T+MOAq59yXnXM/7emSziAlOPGY8eN3zeEdR4ymqbWDy/64QINuRDLE++aN58KjxtDc1smHb3pei8GJhNQVJ06irCCbBWvreGR5TdDhSIj0NZFfARSkMxAJr6x4jO9fOId3zfUO+O+/4Tnue3Vr0GGJyEEyM75x/mwOG1PKhrq9fOpWzVktEkYledl87OQpAPzo/uXqlZcufU3kPw9828wuMLPxZlaefElnkBK8eMz47jsO45J542lt7+QjNz3PX55dF3RYInKQ8rLj/OZ9RzGsIJvHltfw28dWBx2SiPTgkuPGM6wgm0WbdvPq5vqgw5GQ6GsifxMwGrgVWA3U9HCRDBeLGV8/bxafPn0anQ6u+edifvSAegZEom5UWT4/fJe3fsQP7l+mge0iIZSXHee8w0cDcNvCDQFHI2HR10T+c8BVwAd6ucgQYGZ86vSpfOcdhxIz+NlDK/jSPxbR3tEZdGgichBOnTGcK0+cREen45N/fYndezVntUjYXDh3DAD/emmzFnQTALL6cifn3J/SHYhEy3uPGUdlUS4fv/kFbnluA7V7Wvj5e48kPycedGgiMkCfO2M6z67ewcsbd/Ode5bw3XceFnRIItLNrFGlzBpVwqub63lwyTbOOWxU0CFJwPbbI29mD5nZEf15MjM72swePviwJArePHM4N19xLGUF2Ty4ZDsXXf+MZr0QibCcrH1zVt/y3AYeX6GqSZGwufAor1f+toUbA45EwqC30pofAf8ys/vM7HIzK+vpTmY2yr/9MeAO4MdpiFNC6qjx5dz+4eMYXZbPi+t38c7fPMWW3XuDDktEBmjq8GI+dbq3LPzVf19EY0t7wBGJSHfnHT6anHiMx1bU6Hgr+0/knXN3A7OBJ4HvArVmtsrMHjGze8zsSTNbB2wAfgg8BBzinLtzMAKX8JhSXczfPzKfGSOKWV3TyHuve4atu5uDDktEBujKEycxa1QJm3bt5fdPrAk6HBHpZlhhDm+eNRzn4IYn1wYdjgSs18GuzrkG59zXgXHAucDfge1AHFgH3AKcBYxwzv2vc05LAw5RI0rzuOXKecwaVcLaHU2893fPsK1eybxIFGXHY/zPOd6qr799dBU79rQEHJGIdHfVid5Kr39+Zh11Kmkd0vo0a41zrtU5d49z7gvOuXc55850zl3knPuic+5+55ymNxDKCnL4y4eOZebIEtbUej3z2xuUzItE0bxJFZw6o5rG1g5+/t+VQYcjIt0cNqaMk6ZV0dTawR901mxI6+v0kyJ9kkjmDxlZwuraRj5ww3OqsRWJqC+8ZTpmcPOz66lVr7xIqHziVG+l1z88uUblrEOYEnlJuWGFOdz0wWMYX1HA4k31fPzmFzTPvEgEzRhRwqnTq2nt6OTW57QAjUiYzJ1QzpmzhtPU2sG37lkSdDgSECXykhYVRbn88fKjGVaQzcPLavjKna9qBViRCLps/gQAbnpmnX6Qi4TM/5wzk7zsGP9+eTNPr9oRdDgSACXykjaTqoq4/rK55GTFuPnZ9fzl2fVBhyQi/XTClEomVRayZXczjy7XvPIiYTJmWAEfO9krsfnyPxfR1KpS1qEmtIm8mZ1mZs+a2V4zW2dm/2tmcf82M7NrzGy9mTWZ2QNmNiPomOWNjhpfzvcv8FaH/Pq/X2PRRk1slInUXjNXLGacd/hoAO5/dVvA0UgqqL1mlitOnMSMEcWsqW3kG3e9FnQ4Msj6nMib2VvN7Admdr2Z/SHp8vtUBmVmxwP3AkuAtwK/AL4IXOvf5Sv+v38AvAcoBR4ys9JUxiGpcd7ho3nfvHG0dnTy0ZufZ3eTJjnKJGqvme+MWcMBeGjpNjo7VSIXZWqvmScvO85P3nM4OVkx/rpgA/e9ujXokGQQ9SmRN7MvA/8G3g8cDRzRwyWVvgvc75y73Dn3X+fc94GfAKeYWTHwOeBrzrmf+QtQnQkUAx9McRySIte+dSazR5ewoW4v1/xrUdDhSGqpvWa4GSOKGTMsn9o9rby8cVfQ4cjBUXvNQDNGlPDFt3gnTr7491fYUNcUcEQyWPraI/8x4Dqg2jk3xzl3RNLlyFQFZGZVwPH+63Vxzl3tnDsZmAcUAXd2u20n8CjwllTFIamVlx3nVxcdRUFOnLte2cL96jHICGqvQ4OZcdykCgBeUXlcZKm9Zrb3z5/AKdOr2NXUxhU3LlS9/BDR10S+DLjFOdeRxlgSDgUMaDSzf5tZs5ltN7OvmVkMmObfb1XS41Z3u01CaFxFAZ87YzoA/3PHYnbvVYlNBlB7HSJmj/YqKxZvUiIfYWqvGSwWM37yniOYWFnI0q0NfP62VzRb3BDQ10T+EeBNaYyjuyr/743AUuAs4Fd4NXufB0qAFudc8prEDf5tEmKXzZ/AEePK2Fbfwv/9Z2nQ4cjBU3sdImaP9t6uxZvrA45EDoLaa4Yrzc/md5ceRVFuFncv2sKvH03+TSaZJmt/N5jZO7pd/S/wTTMbATwNvKH4yjn3jxTFlO3/vc8593n/3w+bWSXel813gZ5+YhrwhkmOzexK4Mrk7ePGjUtNtNIv8ZjxvXcexlk/fZy/LljP5fMnMHV4cdBhycCltL2C2mxYTan22una2kacc5hZwBHJAKi9DgFTqov58bsP54obF/L9+5ZxyIgSTplRHXRYkia99cjf3u3yfSAX+AjeL/nbky63pTCmPf7f/yRtfwCvdm8XkGtm2Um3FwFvOOfrnLvOOTc3+VJVVZV8Vxkk04YXc9Ex4+h08J171SsfcSltr6A2G1YleVkU5Waxt62DXZp5KqrUXoeIN88czmfePA3n4JO3vMjqmj0HfpBEUm+J/MR+XCalMKaV/t+cpO2JL5Y2vN6BiUm3TwKWpTAOSaNPnT6Votws/rt0O0+tqg06HBk4tdchwswYXZYPwKZdewOORgZI7XUI+fgpUzhz1nAamtu58s/Ps7d1MIY5ymDbbyLvnFuXuACXAW3dt3W7DeCzKYzpNWATcGHS9rcCm4FbgGbg/MQNZjYMOAl4KIVxSBpVFuVy1Yne77+fPrgi4GjkIKi9DiGjyvIAJfIRpvY6hMRixg/fdThTq4tYuX0PP3lwedAhSRrsN5E3s3L/UgF8FZjZbVvXBTgDuCJVATnnOoEvA+ea2a/9Fei+g/dj4uvOuXrg53g1+58zs3PxThPWA9enKg5Jv8uOn0BxXhbPrqlj4dq6oMORAVB7HVpG+T3ym5XIR5La69BTlJvF9y+cgxlc/8QazTqVgXorrfkLUANs96/f519PvvwGeDiVQTnnbgQuAk4A7gYuAD7snPutf5cvAz/CW7jiZrzavdOdc/qERkhJXjaXz58AwC8eXtn7nSW01F6HjtHDlMhHndrr0HP42DIunz+Bjk7H1f94hQ6tzpxR9jtrDfAh4HS8erk/AN/kjXPLduANjkn5KTfn3F+Bv+7ntnbgav8iEfb+4ydy/eNreGRZDatq9jC5qijokGQA1F6HhpGlXmnNlt3NAUciB0Ptdej53BnTuW/xVhZvqueeRVt425xRQYckKdJbjfwm59yfnHM3AO8HfuZf7365yTl3l3NO3TMyIOWFOZzrf6H89dn1AUcjIr2pKMwFoK4xeZpxEQmzwtwsPnbqFAB++fBKLRSVQfq0IJRz7k/AbjO7xMyuM7PbzexXZnaBvxqcyIBddKw33/DtL2ykuU2j6kXCqrzQm+xEibxI9Fxw1BiGl+SydGsDDy3ZfuAHSCT0KQk3s/F4o91vAE4BRuENcv0bsMAf1S4yIIeNKWX26BJ2NbXpy0UkxCqKvER+hxJ5kcjJzYpzxZu82eJufGbdAe4tUdHX3vRf+H8Pc85Ndc7Nd85NAeYApcCP0xKdDAlmxvmHjwbgnkVbAo5GRPYn0SO/s7FVp+ZFIugdR44hHjOeWlnLbi3slhH6msifBHzROfdq943OuUXANcC5qQ5MhpazDh0JwH+XbteiFSIhlZsVpyg3i/ZOR/3e9qDDEZF+Ki/M4diJ5bR3Oh5aui3ocCQF+prI1+Et0bw/msJADsrosnwOH1vG3rYOHltRE3Q4IrIfiV75HY0tAUciIgNx1uwRACplzRB9TeSvAb5rZm/uvtHM5gLfQdNUSQqcMr0agKdX7Qg4EhHZn9L8bAAamtUjLxJFR40vB2DJ1vqAI5FU6GsifzVej/x/zKzOzF41s63As8AE4BdmVu9ftGiEDMj8KRUAPLWqNuBIRGR/inK95UeUyItE06SqQmIG63Y00dKuUtao621BqO5uT2sUIsCcMWXkZ8dZvm0PtXtaqCzKDTokEUlSnJdI5DVQTiSK8rLjjK8oZE1tI2tqG5kxoiTokOQg9CmRd879b7oDEcnJijFnbCnPrK5j0cbdnDKjOuiQRCRJcZ5Ka0SibsywfNbUNrJlV7MS+Yjr82JOZlZqZtea2cNmtsTMZpnZF83szHQGKEPL7FGlACzepAotkTBK9MjXq0deJLKqir0z3tsbNFdJ1PV1QagJwCLgM0A9MA3IBQ4D7jKzs9IVoAwts0f7ifxmJfIiYVTiJ/J7WtQjLxJV1cV5ANQ0aPapqOtrjfxPgS3AaXhTTbYCOOcuNrNs4CvAvWmJUIaU2aO9U3yLN2k0vUgYFeVpsKtI1FV39cgrkY+6vpbWnAp82zm3B0hezu+3wOyURiVD1sTKIvKz42zatZedWgZeJHT21cirtEYkqqpL/ES+Xol81PU1kW8F8vdzWzmgT4KkRDxmzBzl9covUp28SOgUq0deJPKqilQjnyn6msjfDXzTzKZ22+bMrBz4EnBfyiOTIeuQkcUALN/WEHAkIpIs0SOvGnmR6Kou8Wvk96gfNur6msh/Fq/X/VXgFX/b74FVQCnw+dSHJkPVxMoiANbuaAw4EhFJllgQql498iKR1VUjX9+Cc8kV0xIlfZ1HvsbMjgIuA04GNgG7gT8Bf3DOqetUUmZSZSEAa2qVyIuETYkWhBKJvMLcLApz4jS2dlDf3E5pfnbQIckA9XXWGpxzzXgDW3+bvnBEYIKfyK+tbQo4EhFJpgWhRDJDVXEujTuaqGloUSIfYQdM5M2sFLgQmA8M9zdvBJ4A/unPZCOSMmOG5ZMVMzbv3ktzWwd52fGgQxIRX2Kw6x4l8iKRVl2cx9odTWxvaGZKdVHQ4cgA9Vojb2bvBtYA1wGXAEf5l/fjldWsMbN3pjtIGVqy4zHGlhfgHKzboV55kTApyIkTM9jb1kFbR2fQ4YjIAFX5U1BqUaho228ib2anAjcDL+MtBJXrnBvhnBsBFANn4g18/auZzR2MYGXomFBRAMA6DXgVCRUzU3mNSAboPuBVoqu3HvnPAQ87505xzj3snOvqenHOtTjnHnDOnQY8hmatkRQbUepNjbVNPQUioVNemANAXaPap0hUVfmJvKagjLbeEvmjgd/04Tl+BxyTmnBEPFXF/hy39VqsQiRsEon8jj1afVkkqqr94+x2HWcjrbdEvgzY0ofn2AiMTEk0Ir7hfu3eNp3yEwmdfT3ySuRFoqqrtEZnviOtt0Q+DvRlouB2QPMWSUp19RRo+WiR0Kks8nvklciLRFZ1iRL5THCglV213JcEQj0FIuGlHnmR6Ksq0qw1meBA88j/0Mx2HeA+ZakJRWQf9RSIhFd5odc+lciLRNewghyyYsbuvW1asyXCeuuRfwzowJtqsrdLh39fkZSpLMrFDGr3tNCuuapFQqXC75Gv1WwXIpEVixnDSxIDXtWWo2q/PfLOuZMHMQ6R18mOx6gozKF2Tys7Glu7vmxEJHgVRSqtEckEI0rz2LRrL1t272Wcv36LRMuBauRFAlNVrJ4CkTBSjbxIZkis2bJVU1BGlhJ5Ca3EYhU6fS8SLhV+jbxmrRGJtlF+Ir95lxL5qFIiL6GlOlyRcBpW6M04vLOxFec0uZlIVI0ozQdg6+69AUciA6VEXkJLp+9Fwik3K05xXhbtnY76ve1BhyMiAzTS75Hfsls98lGlRF5CSwPqRMKrMjEH9R4lACJRpRr56FMiL6G1r7RGibxI2CQWbdumwegikTXKL61Rj3x0KZGX0KroWnRGiYJI2CR68rapJ08ksqqKc4nHjNo9LbS2a82WKFIiL6FVrtIakdBKrO2gHnmR6IrHjOriXJzTj/KoUiIvoaXSGpHw2ldao4O/SJSpTj7alMhLaFUUJUprlMiLhM2+Hnkd/EWiTDPXRJsSeQmtwpw4OVkx9rZ10NSqKe5EwkSJvEhmGKm55CNNibyElplR6ZfX7FB5jUiojFCNvEhGUI98tCmRl1DTgFeRcKou8Urftjc009mp1V1FoipRI79llxL5KFIiL6FW7k9BuUNTUIqESl52nNL8bNo6HDub9ENbJKq6euRVJhdJSuQl1FRaIxJew0u0KJRI1I1QjXykKZGXUCtPJPIqrREJHQ14FYm+6uJcYgbbG1po69CiUFGjRF5CTVNQioRXIpHX/NMi0ZUdj1HlLwq1vUFn16JGibyEWoVKa0RCa5RfW7t5l07Ji0RZYgrKLWrLkaNEXkJtX2mNeglEwmZMeQEAG+qaAo5ERA7GqDL/R7mmoIwcJfISahWaflIktMYO8xP5nerFE4myUX6PvM6uRY8SeQm1isT0kyqtEQmdseXewV898iLRNqpMiXxUhTKRN7O4mX3BzFaa2R4ze9bMTu12u5nZNWa23syazOwBM5sRZMySHoke+R2NLTinRWfCSO116BpRkkc8ZmxvaKG5rSPocKQP1F6lJ12lNVoUKnJCmcgDnwe+DfwBOB9YBfzHzI7wb/8KcC3wA+A9QCnwkJmVDn6okk4FOXHysmM0t3XS1KpEIaTUXoeorHisKwHYqPKaqFB7lTdQj3x0hTWRvwy42Tn3befcg8AlwFbgg2ZWDHwO+Jpz7mfOuTuBM4Fi4IOBRSxpYWZd5TWqkw8ttdchLFEnv3GnymsiQu1V3iAxa81mLQoVOWFN5HOB+sQV51wHsBsoB+YBRcCd3W7fCTwKvGVww5TBkCivqd2jmWtCSu11CNOA18hRe5U3qCjMIScrxq6mNppa24MOR/ohrIn8L4FLzOw0Mys1s08Bs4BbgGn+fVYlPWZ1t9skgySmoFSPfGipvQ5hiQGvGzXgNSrUXuUNYjHrti6E6uSjJCvoAPbj18CpwIPdtl3rnLvTzL4EtDjnkrO6BqAk+YnM7ErgyuTt48aNS2G4kk6auSb0UtZeQW02asYm5pJXaU1UqL1Kj0aW5rN2RxObd+1lSnVR0OFIH4UukTczA+4DZgIfBZYApwNfNbNdgAE9TV9iQGfyRufcdcB1ydvnzp2rKVAiojJRWqNFoUIn1e0V1GajZkyitKZOpTVhp/YqvUkMeN2iOvlICV0iDxwPnAC8yzl3m7/tETPLAv4P+DKQa2bZzrm2bo8rwqvzkwzTVVqjHvkwUnsd4sYO8+eSV498FKi9yn6N9meg2qTSmkgJY438WP/vM0nbnwAK8HoLDJiYdPskYFl6Q5MgVBT5pTWqkQ8jtdchrqo4l1x/kFxDc9uBHyBBUnuV/RqpKSgjKYyJ/HL/7/FJ248F2oF/AM14898CYGbDgJOAhwYhPhlkFYWatSbE1F6HODNjjN8rr7nkQ0/tVfZLpTXRFLrSGufc82Z2N/ArMyvHq+E7Gfgi8FPn3EYz+znwTTPrxPtiugZvOq3rAwpb0igx/aRmrQkftVcBb8DrqppGNtQ1ccjIHsdESgiovUpvRmt110gKXSLvuxD4Jt4XSDmwAvgk8Fv/9i/jDbz5HF7t3lPAZc451fBloK7SGtXIh5Xa6xCnueQjRe1VetS1KNSuvTjn8MZGS9iFMpF3zu0FPutferq9Hbjav0iGq+g2j7y+XMJH7VW65pLXgNfQU3uV/SnMzaI0P5vde9uoa2zt6kSTcAtjjbzI6+RlxynMidPa0UlDi1acEwkbTUEpkhlGalGoyFEiL5FQ7tfJq7xGJHwSpTXqkReJttGJmWs04DUylMhLJCRWd63TolAioZOYtWaTpq0TibRRmoIycpTISyR0re6qHnmR0CkryCY/O05Dczv1mkteJLJGds1co0Q+KpTISyQkVndVaY1I+JhZVwKwRbW1IpHVVVqjdhwZSuQlEvZNQanSGpEwGq1T8iKRl2jHG9WOI0OJvERClZ/Ia3VXkXAaVapBciJRN7Y8MQOVBq5HhRJ5iYTKYi+Rr1EiLxJKGiQnEn1VRbnkZsWoa2xlj6Z7jgQl8hIJXT3yDaqRFwmjkVreXSTyYjHrmoVKvfLRoEReIqGq2Bvsqh55kXBK1NZqCkqRaBun8ppIUSIvkVDZ1SOvRF4kjFRaI5IZEnXy65XIR4ISeYmE0vxssuNGQ0s7zW0dQYcjIkkSS7tvq2+mo9MFHI2IDFSiR37jTv0ojwIl8hIJZtbVK1+jXnmR0MnLjlNZlENbh9PsUiIRNmaYeuSjRIm8REZVsaagFAmzkaWqkxeJOtXIR4sSeYkM9ciLhNsoLe8uEnljy/1Za3Y24ZzK5MJOibxExr5FoTQFpUgYacCrSPQV52UzrCCb5rZOzRQXAUrkJTIqE1NQqkdeJJRGdyXymkteJMq0wmt0KJGXyOiaglI9BCKhlKiRV4+8SLSNHZZI5NWWw06JvERGhZ/I1zWqtEYkjLpq5Hfr4C8SZZpLPjqUyEtkVBR6pTU7GtUjLxJGKq0RyQxdA16VyIeeEnmJjHI/kVePvEg4VRblkh036hpbtXCbSISNU498ZCiRl8hQIi8SbrGYMaJUU1CKRF2iRl6ru4afEnmJjGEFXiK/s6mNTi0BLxJKo0pVXiMSdSP98S5b65tp7+gMOBrpjRJ5iYycrBjFeVl0dDrqm9uCDkdEejBac8mLRF5uVpzhJbl0dDq27NaP8jBTIi+Rsm/Aq8prRMIosSjUJiXyIpE2xi+vUVsONyXyEimqkxcJt0Qiv0VTUIpE2phhXltWnXy4KZGXSCkv9OaS37FHibxIGCVqa1UjLxJtiTK5jTs1c02YKZGXSKlQj7xIqCVmu1hT2xhwJCJyMMZo5ppIUCIvkTKsK5HXolAiYTShooCcrBibdu3VoHSRCNtXWqMe+TBTIi+Rsq9HXgmCSBhlxWNMH14MwNItDQFHIyIDNXqYBq5HgRJ5iZRy9ciLhN4hI71EfsmW+oAjEZGBGl2Wj5k33qVNc8mHlhJ5iZTyIk0/KRJ2M0aUALB0qxJ5kajKy44zqjSfjk7H+jqV14SVEnmJFA12FQm/Q0Z6ifxrKq0RibTJ1UUArNq+J+BIZH+UyEukaB55kfBLlNYs21qvU/IiETapshCA1ZqFKrSUyEukVPjzyNc1tuKcCzgaEelJWUEOkyoLaW7r5LXNKq8RiapEj/zqGvXIh5USeYmU/Jw4edkxWto7aWrtCDocEdmPoyeUA/Dc2rqAIxGRgZrs98ivqlGPfFgpkZfI6d4rLyLhNHfCMECJvEiUddXI1+zRWfCQUiIvkZOok9fMNSLhleiRX7h2pxIAkYiqLs6lOC+LXU1t1O7RMTeMlMhL5GgueZHwG19RQFVxLjsaWzVQTiSizIxp/gJvy7dpFqowUiIvkZOYgnKHegdEQsvMONovr1mwRuU1IlGlRD7clMhL5FQVezXy2xvUIy8SZvMmVQDw5MragCMRkYGaPtyrk1ciH05K5CVyRpbmAbBl996AIxGR3rxpahUAT6yspaNTdfIiUZTokV+2VYl8GCmRl8gZUZoPwJZdzQFHIiK9mVBRwJhh+exqauPVzbuDDkdEBmDaCC+RX7FNM9eEkRJ5iZxRZYkeeSXyImFmZl298o+vUHmNSBRVFuVSXphDQ0u7jrshpEReImdkokdepTUioXfi1EoAHlteE3AkIjJQ01QnH1pK5CVyKgpzyI4bO5va2KvVXUVCbf7kSmIGL6zfyZ6W9qDDEZEB0Mw14aVEXiInFjNGaMCrSCSUFmQzZ2wZbR2OZ1fvCDocERmAfQNe9wQciSRTIi+RlCiv2ap6PZHQU528SLRNTwx43a4e+bBRIi+RlJiCcrMSeZHQ66qTX6E6eZEomla9b+aaTk0lGypK5CWSEj3ym3eptEYk7A4fW0ZxbharaxrZuLMp6HBEpJ9KC7IZXpLL3rYONu7UcTdMAk3kzexcM2tI2mZmdo2ZrTezJjN7wMxmJN0n18x+bGZbzazBzG43s1GDG70EaXxFAQDrdigpGCxqrzJQWfEY86d4q7yqvGbwqM1KKnXVyWvAa6gElsib2XzgJsCSbvoKcC3wA+A9QCnwkJmVdrvPb4BLgauB9wNzgHvMLJ7uuCUc9iXyjQFHMjSovcrB2lcnr/KawaA2K6mmmWvCadATef+X/heAh4H2pNuKgc8BX3PO/cw5dydwJlAMfNC/z2S8L5iPOuducM7dDpwNHAacN3h7IkGaUFEIwFr1yKeV2qukyol+Iv/Eilo6VGObNmqzki6aSz6cguiRPwv4EvB54OdJt80DioA7ExucczuBR4G3+JtO9f/e1e0+K4BXu91HMtyIkjxysmLU7mnR3NTppfYqKTGuooDxFQXUN7fzysZdQYeTydRmJS2m+ANeV27XFJRhEkQi/xww0Tn3MyC5W2aa/3dV0vbV3W6bBmx1ziXXVHS/j2S4WMwYX67ymkGg9iop86auVV5VJ59GarOSFlOqvB75VTWauSZMBj2Rd85tcs7t2s/NJUCLc641aXuDf1viPj2d1+l+HxkCxvvlNetVXpM2aq+SSqqTTz+1WUmX0oJsKotyaW7rZJNmjAuNrKADSGK8sQchsb2zH/fZt9HsSuDK5O3jxo0beJQSCokBr6qTD0zK2yuozWay+ZMriMeMFzfsor65jZK87KBDGmp0jJWDMqW6kNo9Layq2cNY/6y4BCts88jvBnLNLPnbvci/LXGf4h4e2/0+XZxz1znn5iZfqqqqUhq4DL4JmrkmaClvr6A2m8mK87I5clwZHZ2Op1ftCDqcoUjHWDkoU6q98hrVyYdH2BL5FXi/+icmbZ8ELOt2nxFmlt/LfWQISJTWrKlVIh8QtVfpt0R5zWPLVV4TALVZOSjd6+QlHMKWyD8FNAPnJzaY2TDgJOAhf9NDQBx4W7f7TAVmdbuPDAETK5XIB0ztVfrtxGmJOnkNeA2A2qwclMnqkQ+dUNXIO+f2mNnPgW+aWSewHLgGqAeu9++zysxuA37nL2CxE/gO8Arwr0ACl0CMLssnLzvG9oYW1dsGQO1VBuLQ0aWU5mezvq6JdTsau86sSfqpzcrBUmlN+IStRx7gy8CP8BatuBmvJu9051z32rz3A7cC38P78nkZONs51zHIsUqAYjFjUqV/mk9fKkFRe5V+iceME6YkpqFUeU0A1GZlwEaU5FGUm8XOpjZ27GkJOhwh4ETeOfc151xR0rZ259zVzrkRzrki59wZzrmlSfdpdM5d6Zwrd86VOecucM5tHtzoJQx0mm/wqL1KqhzvJ/LPrKkLOJLMpjYrqWZmXcfdZVrhNRTC2CMv0meTq7zT8qtqVCcvEhVzJwwD4IV1OwOORET6a+ZIbzmB1zbXBxyJgBJ5iTjV64lEz5SqIkrystiyu1kLy4hEzKxRfiK/RYl8GCiRl0ib7E+FtVpTYYlERixmHDne65V/Xr3yIpEyc5R65MNEibxE2sTKQsxgXV0Tre09LhQqIiE0N5HIr1WdvEiUzBhRjJl3JrylXeOfg6ZEXiItLzvO2GEFdHQ6rfAqEiGJHvmF6pEXiZSCnCwmVhbS3ulYsU1nw4OmRF4iT3XyItFz+Ngy4jFjyZZ6mlrbgw5HRPqha8Cr6uQDp0ReIm/fzDVK5EWioiAni6nVRXQ61dqKRI3q5MNDibxEnnrkRaJpzpgyAF7euLv3O4pIqGgKyvBQIi+Rl5i5ZqV65EUi5dAxpQAs2rgr2EBEpF9mdpuCsrPTBRzN0KZEXiJvanUx4PXId+gLRSQyEj3yr2xSj7xIlFQX5zGyNI89Le0qaw2YEnmJvNKCbEaU5NHc1smGuqagwxGRPpo+opiceIzVNY3UN7cFHY6I9MMR48oAeHHDrkDjGOqUyEtGmDbC65Vftq0h4EhEpK9ysmIcMtJru4tUJy8SKUeM9aaQfXH9rmADGeKUyEtGmD7cq5NfvlWJvEiUHDW+HIBnV+8IOBIR6Y+uHvn1WgsiSErkJSNMG64eeZEomj+5AoCnVimRF4mS2aNLyYoZy7c10NiitSCCokReMsKMEd4I+uVK5EUi5ZhJ5cQMXtqwS8mASITkZcc5ZGQJnQ5e1sxTgVEiLxlhSnURZrC6ppHW9s6gwxGRPirJy+bwsWW0dzoeX1EbdDgi0g9Hjffq5BeuVXlNUJTIS0bIz4kzvryA9k7HmtrGoMMRkX44feZwAB5csi3gSESkP46Z6I9xWaPSuKAokZeMoTp5kWh68yFeIv/fpdu1FoRIhCQS+efX7dTZ8IAokZeMMd2fglIz14hEy5TqIiZUFFDX2KrZa0QipLIolynVRTS3dbJIC7sFQom8ZAz1yItEk5lx7pxRANzx0uaAoxGR/jhW5TWBUiIvGaOrR16JvEjknHu4l8jfs3gLzW0dAUcjIn3VVSe/ui7gSIYmJfKSMSZUFJIdN9bXNdHUqmnsRKJkSnUxs0aV0NDcrkGvIhEyb5K3FsRza+tUJx8AJfKSMXKyYkyqLMI5WLl9T9DhiEg/vWvuWABuWbAh4EhEpK+Gl+QxbXgRTa0dPL9O01AONiXyklGm+eU1yzTgVSRyzj98NLlZMZ5YWcv6HU1BhyMifXTi1CoAHltRE3AkQ48Secko04cXAaqTF4mi0oJszj50JAB/W6heeZGoOHGan8gvVyI/2JTIS0bZN3ONSmtEoug9R3vlNbc9v4H2DtXbikTBMRPLyc2K8ermemoaWoIOZ0hRIi8ZRXPJi0TbMRPLmVRVyLb6Fh5ept49kSjIy45zrD/o9XGV1wwqJfKSUcYOKyA/O87W+mZ2NrYGHY6I9JOZdfXK37JgfcDRiEhfnTi1EkA/wAeZEnnJKLGYcejoUgBeWK/R8yJR9I4jx5AdNx5etp0tu/cGHY6I9MEZM0cA8N8l27QWxCBSIi8Z56gJwwBYqGmwRCKpsiiXN88cTqeD2xZuDDocEemDcRUFHDq6lMbWDh7VoNdBo0ReMs7c8X4iv1arzIlE1XuOHgfArc9toLPTBRyNiPRFYtapu1/ZEnAkQ4cSeck4cyeUEzN4acMu9rRohVeRKDphSiVjhuWzaddenlhZG3Q4ItIHb/UT+YdUXjNolMhLxinNz+bwsWW0dTieXb0j6HBEZABiMePdiZVen9OgV5Eo6F5e84gGvQ4KJfKSkd40VYtTiETdhXPHEjN44LVt1O7R3NQiUfDWw/zymkUqrxkMSuQlI5003UvkH1yyHedUXysSRSNK8zhlejVtHY6/P69BryJRcM5hIzGD+1/dyu6mtqDDyXhK5CUjHT6mjOriXDbt2sviTfVBhyMiA/SeY/YNetWPcpHwGzOsgBOmVNLS3sk/X9QP8HRTIi8ZKRYzzpg1HID/vKrTeyJRdcr0KoaX5LK6tpEFazQTlUgUJGadukU/wNNOibxkrLNme3V6d768WdPXiURUVjzGhUclBr1uCDgaEemLN88cTkVhDku3NvDShl1Bh5PRlMhLxpo3qYJRpXlsqNvLs+rJE4msdx/tJfL3LNqimluRCMjJivHOo8YAcMsC/QBPJyXykrHiMev6IrnteX2RiETV2PIC3jRVNbciUZL4Af7vVzbT0Kwf4OmiRF4y2gV+In/voq1aHEokwhJJwd9f2BRwJCLSF5OrijhmYjlNrR3c+fLmoMPJWErkJaONryjkmInl7G3r4O5X9EUiElWnHzKcotwsFm3azZraxqDDEZE+uPhYb9DrTc+s16DXNFEiLxnvwkR5zUKdkheJqrzseNdMVHe+pB/lIlHwltkjKC/MYcmWel7UoNe0UCIvGe/sQ0dSmBNn4bqdLN/WEHQ4IjJA584ZBcCdL29S755IBORmxblwrteZ9pdn1gccTWZSIi8ZrzA3i3cc6X2R/PHJNQFHIyIDdfyUSsoLc1hV08hrW7TQm0gUXOQv6nbXK5vZ1dQacDSZR4m8DAmXHz8BgH+8sIm6Rn2RiERRdjzG2YeOANDgOZGIGF9R2DXr1O3Pq8Q11ZTIy5AwuaqIU6ZX0dLeyV8X6PSeSFSdO2c0AHe9vEXlNSIR8b554wH46wINek01JfIyZHzwhEkA/OmptbS2dwYcjYgMxNzxw6guzmXTrr0s3qTyGpEoOG1GNZVFuayqadRKrymmRF6GjOOnVDBteBHbG1q4Z9GWoMMRkQGIxYwzZ3nlNfe9ujXgaESkL7LiMd5+hDdYXeU1qaVEXoYMM+MDx08E4PonVuv0nkhEJRL5/yiRF4mMxErrd768mea2joCjyRxK5GVIOf+I0VQW5bJ4Uz2PLq8JOhwRGYBjJ5VTmp/Nyu17WLl9T9DhiEgfzBhRwqGjS2lobueB17YFHU7GCDSRN7NzzawhaVu+mX3LzFaa2R4ze9HM3p10n1wz+7GZbTWzBjO73cxGDW70EkV52XGueJPXK//z/65Ur3w/qL1KWGTHY5x2SDWg8preqM1K2Fzg98qrvCZ1AkvkzWw+cBNgSTf9GvgY8BPgfOBx4BYze1e3+/wGuBS4Gng/MAe4x8zi6Y1aMsHF88ZTVpDN8+t28szquqDDiQS1Vwmbt6hOvldqsxJG584ZRXbceHxFDVt3NwcdTkYY9ETe/6X/BeBhoD3ptirgMuCzzrlfOOcedM59ErgH+Jx/n8l4XzAfdc7d4Jy7HTgbOAw4bxB3RSKqKDerq1b+Fw+vCDiacFN7lbA6cVoV+dlxXtm4m0279gYdTmiozUqYDSvM4dQZ1XQ6b4EoOXhB9MifBXwJ+Dzw86TbivF6Au5P2r4MmOj/+1T/712JG51zK4BXgbekOljJTJfNn0BxbhZPrtzB8+t2Bh1OmKm9SijlZcc5eXoVAPerV747tVkJtcRaEP9+RbPHpUIQifxzwETn3M+A1xUoO+dWO+c+4pzbkNjmn8o7C1jqb5oGbHXONSY972r/NpEDKs3P5tL53gIVv/iveuV7ofYqofWW2f7sNYuVyHejNiuhduqMagpz4ry8YRfrdzQFHU7kZQ32CzrnNvXzIf8LzADO9a+XAA093K8BGJu80cyuBK7s4f57zGwHUNvPeKKuEu3z69wA3PCBQYtlsAzkfR6fvGGw2yscsM0u62c8qZBJbSYj92Ud8LcPBxvMQRro+xJ4mx1Ae82kz2BfDKX97fe+jv9emiJJv5QcY1Nh0BP5/jCzLwLXAD90zv07sZmkXoZu29+wXKdz7jrguv08/0Ln3NwUhRsJ2uehIYh9TkV7hd7bbBAy6fOjfQmnoPYl3cfY/bxmxrxvfTGU9lf7GoxQJvJmZsAPgU8Dv8Kr9UvYjVfnl6zIv01EBpHaq0i0qM2KZI7QLQhlZjHgRrwvmG875z7mXj/Z9wpghJnlJz10Et6AHREZJGqvItGiNiuSWUKXyOP1ErwPb3qsa3q4/SEgDrwtscHMpgKz/NtEZPCovYpEi9qsSAYJVWmNmR0JfAp4AHjKzOZ1u7nDOfecc26Vmd0G/M7MSoGdwHeAV4B/DXbMIkOV2qtItKjNimSeUCXyeKPmDXizf+muEa9GD7yV5n4MfA/vrMKDwCedcx39fL3QDKgbRNrnoWEw9nmw22tQMunzo30Jp8HalzC02Ux63/piKO2v9jUA9vrSOBERERERiYIw1siLiIiIiMgBKJEXEREREYki51xoL3j1fA1J2wxvAYv1QBPeoJ0ZSfcpx5sbdx1QDzwNnNbt9hvwFrzo6fJwL/HEgc8AS/DqCV8DPo5fouTfZ+5+nvcHEd3niv085vb+xBeVfQZO7uUxDhgf5vfZv89Y4G9ADbAN+DNQ3YeYrsCbem6v/7zH9Te+MFzC9HnqIbZcvO+PG6K6L8BpwLP+52Qd3sqg8ajtC973+ReAlcAef59ODfJ9cWlsv1FtP6TxuBvCfU3r8TZM+0uaj7Xp2td0ttGBvrf9buSDdQHm+/+Je5K2f9X/T/ik/0YtADYBpd3+I/4LbAQuB84A/gp0JP7TgMnAvKTL1f6H4+JeYvoa0Oz/R5/mX28HvtDtPh/AOygkP/+4iO7zqf59zkh67NS+xhelfcZbnjz5MSfjLcV8PxAL+fucjTe7xErgAv+yGniGXpIt4FL/eb4KnA3c68c4MRXv81D+3kiK49v+/W+I4r4AxwOteAfuU/EWEmoGvhrBfbka7/v7y8DpwM3+vh2Rie03qu2HNB13Q7qvaTvehm1/SeOxNs37Grpj7EEdNNNxweux+gLQAtR1fyPwVptrAL7Ybdsw/z/jM/71o/0PT/dfizH/P/5vvXyg1gN/6iWumP8630ja/ktge7frPwGeyYR99u/3/4Ctvdx+wPiits89PO4neL+8q8L+PgPH+vc5tdt9zvC3zd1PXAasBX7dbVs23pfTzw7mfR6sSxQ+T8AReAekGnpJ5MO8L8DjwF1J274LPBLBfVkC3Njtetx/3C+C2hfS1H6j2n5Iw3E3rPvq3+//keLjbZj3t4fH/YSDONYOxr4SwmNsGGvkzwK+hNfT8/Ok2+bhTY91Z2KDc24n8CjwFn9TJ3A98GS3+3Ti/XqauJ/X/BLeB+/z+7kdoBRvNbx/JG1fBlSZWaF//TC8N70/wrrPcOD96Ut8PQnzPncxs5l4p3Gvdc7VdLsprO9zrv+3vttz7/D/lu8nrinA+KTXbgPu7vbaA32fB0uoP09mlgX8Afg+Xg9L5PbFzKrweuRfN+2ac+5q59zJUdoXXy7d2onzplbczf7bSZTb74GE9X1Kx3E3rPsK6Tnehnl/u6ToWBvlNjrwY2x/fi0NxgUYDZT5//4ar/9F9TG8Xz05SY/5KbC2l+cswTtd8+cebhuBdyrjSwOM9wFgQ7frNXinTF7CO027ErgsqvsMLASeAJ7CO725Ee8Xrx1kfKHd56TH/ROvJjOWtD2U7zPe2hAv4p2aHOdfHsTrHSnYz3O81X/tqUnbP413KjA+0PgG6xL2zxNeWcCrQI7/mbkhavvCvtP+pwH/xvs+2O7HGIvSvvj3/Sxe4n4aXsL4KbwD+blB7Qtpar9Rbz89PH7Ax90w7ytpON6GeX+THnfQx9rB2FdCeIwN24JQOOd6660qAVqcc61J2xv82/bnl3hf1D/q4bYP49Xb/bo/cQKY2Yfwais/6V8fBVQCU/F+Fe4E3gvcYGbOOXdjT88T1n02sxgwE2+A0efwPqhn463ylwd8faDxhXWfuzOziXh1alc671d5Ynto32fnXLuZXYH3xbfOv08dcLJzrqmX1068VvJrx4DCg4hvUIT582RmM/Dre51zrWbW6/1DvC9V/t8b8erJfwScBFyLd+D+XvIDQrwv+Pc5Fe8gnHCtc+7Onu4c8fZbTy9C/j69zsEed8O6r+k63oZ1f7tL1bE24m10wMfY0CXyB2B4v1h62t75ho3eEfMXwPvwVqV7sYfbP4RXx7WrX4GYXQz8Brjdfw2AXXinQF5xzm3xtz3ofxi/incA7K8g99mAc4D1zrmV/raHzawI+KKZ/V9/4+ujsLzPV+B9cdyUtH0XIX2fzexw4BHgBbzEyuEdFO4zsxO7vY/Jr0EPr5/Y3tnf+EImsM+Tf3D+PfB759zTA4r+jTEH1Tay/b/3OecSp8wfNrNK4Foz+4Hr38qfQb4vBtyHlzh9FK9e/nTgq2a2yzn3y37sR8r2JY3t92CE5ft4MI67Q+14G5b3djCOtWFvowN+b8NYI9+b3UCumWUnbS/yb+tiZjnALXhf0lc755LrpcAb2DDav1+fmdmn8aYbugtvRLYDcM41Oefu6/aBS/gPMMlvkP0V2D475zqcc//t4YP5H6AAr+6rz/H1QyjeZ+B84F/OuZbuG0P+Pn8cr0fnbOfc3c65e/B6dVqAr/Ty2uANtkl+7U7/+dLxPg+WID9Pn8CrjfyKmWX5tfL+S3X9uz+C3Jc9/t//JG1/wH/9CX14ju6C3JfjgROADzvnfu2ce8Q5dy1er9v/DaANh739HoxQfB8P0nF3qB1vQ/HeMjjH2rC30QG/t1FL5Ffg/TqZmLR9Et7gFwDMLB+4B29aoI84595wytf3Frw5QJ/cz+1vYGbfxvuy/zNwQffTIGY2zcw+bGa5SQ/Lxzv1PJAv1MD22cxGmdmV5g1y6y7f/1vb1/j6KQzv8zjgEN44yCrs7/NYYLFzLpF04ZxrBp7H633c32snXusNr+0fMNPxPg+WID9Pb8c7sNUBbf5lDt5UZG1mNqHPe+EJcl8SCUZO0vbEgaen3qTeBLkvY/2/zyRtfwIvaZrQh+foLuzt92CE4ft4sI67Q+14G4b3drCOtWFvowN+b6OWyCcGgJyf2GBmw/DqNB/qdr+/+Nve65z7TS/PdwywoK9fdGb2Kbw6rZ8Clzvn2pPuMhqvLuzsbo8x4B3A4wP8Qg1yn3OB3+KdWuruncBy59zWfsTXH4G+z90eA94iMcnC/D4vBw7r3lPh9y4cAazZz2uvADYkvXY23gCdxGun430eLEF+nq7C66XqflmO16t4NLC5z3vhCXJfXsObcefCpO1vxduPtX14ju6C3Jfl/t/jk7Yfi1ffu7EPz9Fd2NvvwRhKx92hdrwdSsfasLfRgb+3rh+jigf7QtKoY3/b/+GNXP4c3uCIZ/G+dEv929+O1zP0J964gMBhSc+1Fvh+L68/D5js/3uk/5/8Sg/POw9vvEEcb57lbXiLGJwF/Mt/XI/zi4Z5n/3rN+OdUv8U3lypvyNpZocDxRe1fe4WU81+7h/a9xlvYY7dwGP+c5yDVwv8uth6eJ8/6r+v38L70rwHb4DcpFS9z4N1CePnKen2l+j7yq6h2he8MwkO7+B6Gt5AvE7gqgjuy114Z0o+CpyCt0JtK31YnTld+0Ia228U2w9pPu6GaV/962k93oZtf7vFlPJjbbr2lRAeY/vVwAf7sp83IgtvAZKt/gf+frotYUvvSwUvTnquJuCaXl7f4R9w8Vb52t/zOqDSv1853mCcjXinfp4E3hTFffav5+OtRrnG/6C+CLy9P/FFbZ/9bb8CVvTymNC+z3jz7t7rP0ctXj3h4X3Y58+yb2nop3jj8tEH9T4P1iWMn6ek21/q7faw7wverBGL8L4PVuDNNBG5fcH7bvsh3lmGvXjJ4ofBm+ovqH0hTe03iu2HNB93w7Sv3T6TaTvehm1//W1pOdamc18J2TE2MTepiIiIiIhESNRq5EVEREREBCXyIiIiIiKRpEReRERERCSClMiLiIiIiESQEnkRERERkQhSIi+R5i8OISIRoTYrEh1qr+GnRF72y8weMTPX7dJuZrVmdq+ZnRZAPDeY2eJu168AvjHYcYiEldqsSHSovUoqKJGXA3kSOM6/nIK3Olku8ICZvXeQY/kGcFG369cAZYMcg0jYqc2KRIfaqxyUrKADkNDb5Zx7pvsGM7sd+C/wazP7j3Nu52AE4pxbNRivIxJxarMi0aH2KgdFPfLSb865TuDrQClwIYCZVZvZjWZWZ2Z7zOxOM5uYeIyZfc3MFprZe81suZk1m9lzZja/230Kzex6M9tiZnvN7AUze0e327tO+5nZWmA88DH/lOSh/t8LusdqZheZWauZVaTz/0QkzNRmRaJD7VX6Q4m8DNSjQAcw38zygYeBE4BPAJcAI4DHzGxYt8dMw/ty+hrwTiAfuM3MEmeGfgCcCnwSeCvwmn/7IT28/tuBrcDtwHHOuUXAS0DyqciLgbudczsOZmdFMoDarEh0qL1Kn6i0RgbEOddhZjuA4cClwHRgtnNuKYCZPQSsw/vS+br/sGLgdOfcAv8+ceAOYA7wPHAi8IBz7jb/9ieAbfTwOXXOvWhmLcC2bqcl/wR818xKnXO7zawSOAN4d8r/A0QiRm1WJDrUXqWv1CMvqXAKsAJYaWZZ/q//JuBxoPvI+3ZgYbfrG/2/hf7fp4Ar/FOGVwKVzrnP+j0BfXEzEMfrSQB4D9AA3N3fHRLJcGqzItGh9ir7pUReBsTM8oByYBNQAcwA2pIubwNGdntYi1/7l5D4d+Jz+Em8UfOzgd8CG8zsdjMr6UtMzrntwH/Yd+rvYuBvzrmW/u2dSOZRmxWJDrVX6Ssl8jJQb8I7HfcEsBt4GTi6h8s7+/qEzrm9zrmvOucm4X1pfQU4B/heP+K6ETjVzGYD84A/9+OxIplMbVYkOtRepU9UIy/9ZmYGXA3UAf/Am2f2zcBa51xtt/vcBCwGDnjazq/lexm43jn3E+fcMuBbZnY6MG4/D+voYdudwB7gV8Aa59yT/dg1kYykNisSHWqv0h9K5OVAysxsnv/vLGAM8CHgJOAi51y9mf0B75TdA2b2Hbwvnyvxegre1pcX8Qf2PAt81cyagaV4v/bfBFy1n4ftAo4ysxOBx52nxcxu9R/z9f08TiSTqc2KRIfaqxwUJfJyIMcDT/v/bgO2AAuA451zzwL4XzQnAt8HfoO3Kt1i4Dzn3D39eK1PAo14q8lV443I/6xz7vf7uf+3/df7D960W4mBPffifcnc1I/XFskUarMi0aH2KgfFnHNBxyCSUmb2K+Aw59wJQcciIgemNisSHWqv4aIeeckYZvYB4Ai8U47vCTgcETkAtVmR6FB7DScl8pJJ5uItnPFz59ztQQcjIgekNisSHWqvIaTSGhERERGRCNI88iIiIiIiEaREXkREREQkgpTIi4iIiIhEkBJ5EREREZEIUiIvIiIiIhJBSuRFRERERCLo/wPbyoL6uqv6uwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig,ax = plt.subplots(1,3,figsize=[12, 7])\n", "\n", "prof=[15,19,20]\n", "\n", "dat = bad\n", "\n", "for i,val in enumerate(prof):\n", " \n", " ax[i].plot(dat.density.isel(time=i), dat.depth)\n", " \n", " xmin=dat.density.isel(time=i).min()-0.05\n", " xmax=dat.density.isel(time=i).max()+0.05\n", " \n", " ax[i].hlines(y=dat.mld_03.isel(time=i), xmin=xmin, xmax=xmax, colors='orange')\n", " ax[i].text(xmax-0.2, dat.mld_03.isel(time=i)-4, 'MLD', color='orange', fontweight='bold', fontsize=18)\n", " \n", " ax[i].set_ylim(120, 0)\n", " \n", " ax[i].set_xlabel('Density')\n", "ax[0].set_ylabel('Depth (m)')\n", "\n", "plt.savefig('../figs_submission2/MLD_03_QI_bad_profiles.png')" ] }, { "cell_type": "markdown", "id": "d1de8d7a-ef42-456c-9281-2a4f3ccf12c1", "metadata": {}, "source": [ "##### def rmse(data, data_mean):\n", " \n", " return np.mean(np.sqrt((data-data_mean)**2))\n", " \n", "rho_mld =\n", "\n", "QI = 1 - ( rmse() )" ] }, { "cell_type": "code", "execution_count": 3, "id": "8b8a8a43-5bab-4b78-8168-d35e97061759", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:      (time: 369, depth: 1000)\n",
       "Coordinates:\n",
       "  * time         (time) datetime64[ns] 2018-12-10T03:00:00 ... 2019-03-12T03:...\n",
       "  * depth        (depth) int64 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999\n",
       "Data variables: (12/25)\n",
       "    temp         (time, depth) float64 ...\n",
       "    salt         (time, depth) float64 ...\n",
       "    lat          (time) float64 -43.01 -43.01 -43.0 ... -43.06 -43.05 -43.07\n",
       "    lon          (time) float64 7.793 7.768 7.791 7.796 ... 7.845 7.812 7.84\n",
       "    density      (time, depth) float64 ...\n",
       "    mld_01       (time) float64 18.5 16.33 13.5 56.0 ... 23.0 12.0 14.0 16.09\n",
       "    ...           ...\n",
       "    v10          (time) float64 -5.411 -7.296 -5.422 ... -5.715 -7.256 -8.585\n",
       "    emp          (time) float64 -1.539e-07 -5.5e-07 ... -5.452e-09 -1.584e-10\n",
       "    qnet         (time) float64 29.03 114.1 336.4 -24.47 ... 479.6 -19.01 -53.45\n",
       "    taux         (time) float64 0.01897 0.006262 0.154 ... 0.02407 0.0184\n",
       "    tauy         (time) float64 0.04108 0.07468 0.04125 ... 0.07386 0.1034\n",
       "    wind_dir     (time) float64 145.8 163.8 117.5 104.8 ... 139.6 150.3 157.1
" ], "text/plain": [ "\n", "Dimensions: (time: 369, depth: 1000)\n", "Coordinates:\n", " * time (time) datetime64[ns] 2018-12-10T03:00:00 ... 2019-03-12T03:...\n", " * depth (depth) int64 0 1 2 3 4 5 6 7 ... 993 994 995 996 997 998 999\n", "Data variables: (12/25)\n", " temp (time, depth) float64 ...\n", " salt (time, depth) float64 ...\n", " lat (time) float64 ...\n", " lon (time) float64 ...\n", " density (time, depth) float64 ...\n", " mld_01 (time) float64 ...\n", " ... ...\n", " v10 (time) float64 ...\n", " emp (time) float64 ...\n", " qnet (time) float64 ...\n", " taux (time) float64 ...\n", " tauy (time) float64 ...\n", " wind_dir (time) float64 ..." ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dat_saz" ] }, { "cell_type": "code", "execution_count": null, "id": "13c38687-b051-445d-a779-7895f16641ee", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "duplessis2021_JGR", "language": "python", "name": "duplessis2021_jgr" }, "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.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }