ccarouge/mrsos.ipynb

## mrsos.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model:  CanESM2\n",
      "Read data...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'mrsos' (time: 1740, lat: 64, lon: 128)>\n",
       "[14254080 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * time     (time) object 1861-01-16 12:00:00 ... 2005-12-16 12:00:00\n",
       "  * lon      (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2\n",
       "  * lat      (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86\n",
       "    depth    float64 ...\n",
       "Attributes:\n",
       "    standard_name:     moisture_content_of_soil_layer\n",
       "    long_name:         Moisture in Upper Portion of Soil Column\n",
       "    units:             kg m-2\n",
       "    comment:           the mass of water in all phases in a thin surface soil...\n",
       "    original_name:     WGFL\n",
       "    cell_methods:      time: mean (interval: 15 minutes) area: mean where land\n",
       "    cell_measures:     area: areacella\n",
       "    history:           2011-04-13T22:40:08Z altered by CMOR: Treated scalar d...\n",
       "    associated_files:  baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation..."
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import time\n",
    "import os\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "from glob import glob\n",
    "from scipy.stats import rankdata\n",
    "#===============================================================================\n",
    "start_time   = time.time()\n",
    "\n",
    "## directories and list of CMIP5 models\n",
    "c5model_list= ['CanESM2','CSIRO-Mk3-6-0', 'GFDL-CM3','GFDL-ESM2G','GFDL-ESM2M','MIROC5'] #'HadGEM2-ES','MRI-CGCM3','CNRM-CM5','BNU-ESM'\n",
    "#idir        = '/short/w35/dh4185/mrsos_merge/'\n",
    "idir = './'\n",
    "odir        = idir\n",
    "ndays       = 30\n",
    "\n",
    "# normalisation parameters\n",
    "climstart = 1861\n",
    "climend   = 2005\n",
    "C0        = 2.515517\n",
    "C1        = 0.802853\n",
    "C2        = 0.010328\n",
    "d1        = 1.432788\n",
    "d2        = 0.189269\n",
    "d3        = 0.001308\n",
    "\n",
    "################################################################################\n",
    "#for model in c5model_list:\n",
    "model = c5model_list[0]\n",
    "print('Model: ',model)\n",
    "\n",
    "infile  = idir+'mrsos_Amon_CanESM2_historical_r1i1p1_18610101-20051231.nc'\n",
    "outfile = odir+'mrsos_Est'+model+'_historical_r1i1p1_'\n",
    "\n",
    "# read in aggregated PET file\n",
    "print('Read data...')\n",
    "ds  = xr.open_dataset(infile).sel(time=slice('1861-01-01', '2005-12-31'))\n",
    "mrsos = ds[\"mrsos\"]\n",
    "mrsos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why add data by working on day of year when the initial data is monthly? Just do the calculation by month then add the day of year to the end result if you want daily data. But since you don't interpolate the values between months, all days within 1 month have the same values...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original version with the loop for reference. Simply changed day of year for month per comment above.\n",
    "# I would not try it on the full array but might be useful to compare on subsets if you're curious.\n",
    "## loop through month and create ranking for each grid cell\n",
    "#start_time   = time.time()\n",
    "#\n",
    "## make climatology for each day from 1867 to 2005\n",
    "#mrsos['month']   = ds['time'].dt.dayofyear\n",
    "#rank           = xr.full_like(mrsos,fill_value=-9999.)\n",
    "#rank['month']    = ('time',mrsos['month'])\n",
    "#months            = np.arange(1,13,1)\n",
    "#for m in months:\n",
    "#    for lt in range(len(mrsos_t['lat'])):\n",
    "#        for ln in range(len(mrsos['lon'])):\n",
    "#            t = mrsos['month']==m\n",
    "#            rank[dict(time=t,lat=lt, lon=ln)]    = rankdata(mrsos[dict(time=t,lat=lt,lon=ln)],method='dense')\n",
    "#print(\"--- %s seconds ---\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- 8.823814392089844 seconds ---\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fb4c0838a58>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEXCAYAAAC+mHPKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsXXd4VMXXfie9NyAJSQghkEDooFKkqigqSFFERGkiKoqKiCKKoiiIoqiAUhRERVAEpSNNUUCRIiDSS0JJg/Te5/tjl8+8M0s2m6KE332fJw+cuXfnzr2795w5XUgpYcCAAQMGDFiD3X+9AAMGDBgwUDNgCAwDBgwYMFAuGALDgAEDBgyUC4bAMGDAgAED5YIhMAwYMGDAQLlgCAwDBgwYMFAuGALjGoUQIkYI0eO/Xsf/MoQQrwshllTT3MOFEDurY24DBqoLhsAwUC4IIcYIIfYJIfKFEIstHH9UCHFaCJElhPhRCBFkZb5BQohjQohsIcQZIUSXUscGmo9lCiGOCiH6lTq20XyNK38FQojDpY7HCCFySx3fXEWPoMIQQoQJIaQQwuE/uvbPQogcIcTxsjYhwoR3hBDJ5r93hRCi1PFr7tka+HdhCAwD5UUcgLcALFIPCCG6AZgGoC8APwDRAJZdbSIhxO0A3gEwAoAngK4AzpqPBQNYAmAcAC8ALwBYKoTwBwAp5V1SSo8rfwB+A/Cdcol7Sp1zR8Vv+brAMgAHANQC8AqAFUKIOlc59zEA/QC0AtASQG8AjyvnGM/2fxiGwKgBEEI4CyE+FELEmf8+FEI4m491F0JcFEI8L4S4JISIF0KMqOo1SCm/l1KuApBs4fA9AL6TUh6RUhYAeBNAVyFEw6tM9waAKVLK3VLKEillrJQy1nwsBECalHKjNGE9gGwA2lxCiDAAXQB8Vamb+2e+BkKIX8yazRYAtZXjHYQQvwkh0oQQh4QQ3Usd2y6EeFsIsUcIkS6EWC2E8DMf/tX8b5p5Z96x1OfeE0KkCiGihRB3VcV9lJo7EkBbAJOllLlSypUADgO47yofGQbgfSnlRfP38T6A4VW5JgM1G4bAqBl4BUAHAK1h2v21AzCp1PFAAN4AggGMBPCxEMLX0kRCiE/MDM/S318VXJ8w/5WmAaC5hevbA7gRQB2zCeuiEGKOEMLVfMo+AMeEEH2EEPZmc1Q+AEtrGwpgh5QyWhn/WghxWQixWQjRyob7WApgP0yC4k2YGOiVdQcDWA+TluUHYDyAlcpufSiARwAEASgCMMs83tX8r495Z/67mW4P4IT5eu8CWFjaBFQaQoh1ZXxv665yP80AnJVSZpYaO2Qev9r5h6ycW9Fna+B6gJTS+LsG/wDEAOhh/v8ZAHeXOtYTQIz5/90B5AJwKHX8EoAO1bSutwAsVsZuA5AEkxnDFcB8ACUAHrTw+SAAEibBUBcmZrkLwNRS54wEkAUT080B0OsqazkNYLgy1sm8BjcAEwEkwMSord1XqPl67qXGlgJYYv7/BABfKZ/ZBGCY+f/bAUwvdawpgAIA9gDCzPdc+jsaDuB0KdrNfE5gFX5XQwDsVsamqt9fqWPFAJqUoiPMaxKVebbG3/XzZ2gYNQNBAM6Vos+Zx64gWUpZVIrOAeDxbywMAKSU2wBMBrDSvLYYAJkALlo4Pdf872wpZbyUMgnATAB3A4DZKfsuTILQCUA3AJ8JIVqXnkQI0RkmzWqFspZd0mR+yZFSvg0gDSazlTUEAUiVUmaXGiv9zOsDuL/0zh5AZ5iE3hVcUD7rCMWspSCh1LpzzP+tyu8tCyY/UGl4wfTdlOd8LwBZUpqkRSWerYHrBIbAqBmIg4lhXUGoecxmCCHmKVFGpf+OVHSBUsqPpZQRUkp/mASHA4C/LZyXCpMguVqZ5NYAfpVS7pMm/8ZeAH8AUKN7hgH4XkqZZW1pYHPZ1RAPwFcI4V5qLLTU/y/ApGH4lPpzl1JOL3VOPeWzhTBpXpUuCW0hOqz038arfOwIgHAhhGepsVbm8aud36qc5wLlf7YGrhMYAqNmYBmASUKIOkKI2gBegymSyGZIKZ+QpaKMlL+r2bYhhHAQQrjAZGKxF0K4XAkTNf+/uTksMxTAAgAfmYWDJXwO4GkhhL/Z1zIWwBU7/F4AXa5oFEKINjDtYv/fh2H2d9wPYLGyxlAhRCchhJN5TS/gH5NX6fDWMAvP5RxMZrI3zJ/vDJMz/wqWALhHCNHT7FtxMQcchJQ652EhRFMhhBuAKQBWSCmLAVyGyUQXfrXnaw1SiQ5T/iw6y6WUJwEcBDDZvN7+MJkNV17lMl8CGCeECBamsOjnYX7G1p6tgf8R/Nc2MePP8h/Yh+ECkwM13vw3C4CL+Vh3ABev9tkqXM/rMO0oS/+9bj7mAxNDz4bJzPI2APtSn30ZwMZStCOAT2AyaSSUvh/z8TEw+ScyYQq3fV5Zy4MwmXyEMt6s1DqSAWwDcGOp413Mz8bxKvcYDmAHTKaZLQDmwOzDMB9vD+AXACkwCYH1AELNx7ab73sPgAwAawHULvXZKebPpMEUwDAcwE7l+hJAoyr+3sLMa8uFycHeo9SxLjCZnK7QAiZzYIr57138478o89kaf/8bf1d+DAYMXPcQQkwCcFlKOb8a5t4Ok3D5rKrnNmDgWsG/nnlqwMB/BSnlW//1GgwYqMkwfBgGDBgwYKBcMExSBgwYMGCgXDA0DAMGDBgwUC5ctz4MJ+EsXeBu/UQDVuHTrIjonGInorMLmS4psr4PCfe6TLSzYE03tpDzzQIcMrQ5TmVxDb0oT56zUFGei5T90YUEpQafBWXbsVYBD8QoaQd2PGexqz3RJUyaLqOM2RUqdDHTIi0btsKhCV+kOJpp6cSvvsjl+5TOjtqcBT58r06xtq3Lv3meNpZYwN+zOFmgnfNvIBOpSVLKqxVltIqet7jL5JRi6ycC2P9X/iYp5Z0VvdZ/ietWYLjAHe3Fbf/1Mq4LnJ/SgmgnRSB45CrMxU7nvLVqc3JxcpIn0R4/s9BplMtz+PyuJ40HJyYRfW5IW6LdO/DxVyI3EN3COZHop5v21K6hpkcWt+UaiE0/5ry23j4HiX5+9mPalPb5TOf5MV1vEz8rue8wbEWn75jxPuW3j+iD+fz8Z13kvMiigSX6pJ7M3IvEadsWdVTfSEQE+itzJmjn/BvYKlecs37W1ZGcUow9m0KtnwjAvu6psrL/r2lctwLDQNVhUgtmtK9uG0B0o2XMnE4/rO9OM7NdiJ5587dEv+jMBVQjJjHTLElK0ea08+IqGvXfZ2YtXJyJXlB8A09QzDvC4uwcaJDMOO3+OkX0xjNNiT4+hF+pwvH6lOnNWKXwqsM79dN1mDE3ZF5fLjztt5/om5Y/T7RdEN+r18+sjbt30HfLeb7M8H1O2CYw2h0o1Mb+TFOESHebprxmIAGUwIKQvc5gCAwDVjFvIguIqF1cHLYkjc1FEYWNtTmyQ9yInnRgKNHuinUjqxkLGDd3pgHgfE8for3O8QvrGc2Tnnyc51jbYxbRB/NCoGLxY32JLlBMTn4/MB0/mrWcsMUx2pyp3eoT7R7H64p72lq1E+t4qMntRNfrxgIgqRkLiOzbWUDf13iPNmcfLxbI4xZ31M4pC3uevEEbK3JnFuRQsYo3/zkkJApl+UxSNRmGwDBgFcktmCl6HWTmX9SEGe2Z+3UNw47dIGj8YSwPZDNzz2sdRnSBP18TAAYP2UZ0UiFrHL/PuInoqDfjiX5q1bO8xgLdlOY4mT/T3pdNYw1dLhEdX8hCbFcfvRqIQxGbXYpmsoWi8JRaL7ACcOLv4NwAFqYDW/9GdHevY0S3dtI1OkfLldfLjU0rvtDGGi3n/kwNt1TqEv8pDA3DgAEAQwYwY14Q0pXoyIVskgpbq3t6W007QPRqbyo+C7cTwUTXX847zbRP9DnbuZ0hetzhgbyup/j4/u4NiHZI5Z+/u4XN7bNBe4leOupuove68LrSGrEvxv+3NG1Ox/dZg7h5+m6iv/7Dtp27JST1a0J05Fy+pvfCXKJfP9GHaJe5ejsVVaA6Yq92Tlm4M/RGbaxxqCI8bZrx2oGERHEVpSgIIRbB1O3wkpSyuXJsPIAZAOpIU6VnCCEmwtQSoBjAM1LKTVWyEAuoNoEhhKgHUzGzQJgKry2QUn5k7kL2LUw1bmIADJTmInVXu3EhxA0wFUFzBbABwLPSSCD51zCx1kmin76Lexm9cxMzgtQiXRtYs5PPiWihOLHZSoMlY9SuqzpWZLIDuqCImXfGiyyEmhxn/4NwdSUaBXqEzpcKI73UiXfurl3YsZ6exqae1BbsXAYAXOTnc+4vXufa+2YSPQ62CxAHJSAprQlrXzv7RhHtFsGaUW4dXUB7n2ShY+sLuCJmpzaWWMIi4snQTjbOeu2gpPJFia9gMUx1zL4sPWjmqbcDOF9qrCmAQTDV+goCsFUIESll9djHqlPDKIKpaNyf5vLK+4Wp7eVwANuklNOFEC8BeAnABCs3PhemfsO7YRIYdwK4WklnA1WMXrewQ7rYmxmtXS47M4+P9tbm8Ipm5+blExxR4hHPv+/BZ0cSfXqIPmexJ38m5Ec2mRR4sYnAzZmd4IXhAUT/+K3WrhwzUiKJvtWdTTcPrH+K6LqNWICknOZrAMCAvmwOWreIW0o835A1OFOVdNvw84zZRL+c2J7ov4fXJfr8XywgHEN0P4r7XNaeLEQMl4nbJo3Vxt57bZ6Ns1ybkACKq0hgSCl/tVRRGcAHAF4EsLrUWF8A30gp8wFECyFOw9SR83cLn680qk1gSCmvVFaFlDJTCHEMphaiffFPLMQXMFXSnICr3LgQIgaAlzS3tRRCfAlTo3pDYPxLSLiVQx8LlE2zRxy/KA4Zuq27SPFZ57djk0haGjOjsAj2aRSf1fsKeQYwU+v1JptIVs7kUFGH3CCiiydxe/KF6bzTB4DHfFibevim/kR7PcBs02MMm8HEAD20/+uQDkQPf3Q70YcH8joyOnN+SXnQvymHlBenpxN96RkW2BHz/yS6xIK2ld+btURbBcbDL27QxsbOGE10bfymnVNTYIOGUVsIUTr2bYGUckFZHxBC9AEQK6U8pHTxDYZpI30FF81j1YJ/xYdhlpZtYGqEE2AWJpBSxgshrnCjq914Ibhz21UfiBDiMZg0EbhAN4sYqBgc72GG5bCcnbRFzoqAsPDeFN/IUThd63Gk1eS6PzIdx76CmOJAbc6sc+wcnp/Wjeha97FAcBrNoaRxq5lprpijh+J/NpyjpFxvZK0meA07PkracksRNUkPAKJe5uvsLmQtJGFQBNF1YLvAkA05EEH8xffuP5v9JnZKPkT8A7qzvlZvxYy41rY19fXQezFFjGcfxgdzo7RzagIkgMLyW8mTpJS6Q+cqMPdXeQXAHZYOX2U51YJqFxhCCA+YGraMlVJmXKXHPXD1Gy/3AzFL6QUA4CX8DB9HFaFZLY4U+vlGzjR7rccqotVIIQBYF0u+O/y8nxnr74ktiXZvx6Yd/zA9aifClxnppVzWQhLWsUDovJ6d9wO8OBdkgNML2jVCVpwnOqEXzxnXibWW9Q/OIPpgvr63WZPMDv8/45i5PxG1nuh1n+gOaGs48Tj7UvreyK/D74ns/0lO42cn7PQsbv83aikj57VzysJjA0ZrYydH8MYuAn/YNOe1AglZZSYpC2gIoAGAK9pFCIA/hRDtYNpAl+70GIIKduMsD6pVYAghHGESFl9LKb83DycKIeqatYu6AK7EJV7txi+a/6+OG/iXkDCUd/eR3rxb/XQHm2lcRulfT/Je3kXf0ZNj+n0cec7vN91MdK2/9ZfxbCEzsJSm7Cept5fNXl96s5nm0zrdiW48j9cEAHBnhuZ3lL3Jqe35mgPfZ6FT9wt9V13UkqO1QjM49fvHVM6st9wavWy81m010cdz2WfRP5Tv9WQt/n5OpeumtPazOVFvT2vbjFL2KboQEoXXSfkeCRRXk7yQUh4G8P8qoNlMf6OUMkkIsQbAUiHETJh8vxEwNfGqFlRnlJQAsBDAMSll6bCPNTD1Y55u/nd1qXHtxqWUxUKITCFEB5hMWkMBsEfPQLWi8wouTdHZnaOmPO30GkEqHCPYAe0imM4s4Z/iG0Ospzfvy2eG1cmFmfcdW4cTfWTUx0RfLGIGdncg5wQAQP3RrMUktFecMTlsogpaw/klsp5uSnNI4uuKFA69lXlK7ZAK4NsHWDjaJfE1Chrxuuzy+T7iB+gm3eX5LHTCbPQ3TN26TBt7cN+jNs1xrcKU6V01EEIsg8nPW1sIcRHAZCnlQovXlfKIEGI5gKMwBRo9VV0RUkD1ahidAAwBcFgIcWU78zJMgmK5EGIkTDrt/YDVGx+Nf8JqN8JweP+rGODFDtGGjsxMVmbx7nR3Fps7AOBgCptmoi+yzTx4DTP/1MZM5zbVhdLtjY8T7RPAJqenP1/O5x9lf0RrP2buwe/pCYdrD3BIe4uPGxH92q1sjnPvwcz+LjeuVwUAA28fQnTmzaxxeO5WfCmKw7o8qD2PtbwZIWzmSi5h4Tr4AEelRY5jMyQA5ERWuDbfVfFBG/6OPkDN9GEAAsUWree2Q0r5oJXjYQo9FcDUKrm4FVy3/TC8hJ80ig9WDRyC2U5/bkgY0S6pfH6xzndRoETF+u/nUFGnDKZju7FQ6txfNxftuKA7ZkujsJCFTgN/doKLO5iZyyI9fFU48M2UdGRfTE4Ah+qmRSi1kSy8XnmBZe9FJ/RcQ/SKKP+rnHl1pA/l3A2/PxUfkGI/yQ1nv5PrzhPanCVKrS1Lz6tMdGilDXX/lJ3v21u4auf8G9gqV+y3xRGtonlLJ7lyfflqCjYJja/Utf5LGJneBqwirn8Y0TlNOOTSfxnvrJwv60X87OKZWctMJc6/ATt+ndKV5LYneBcOAP7+zFzsC5kRn32Y1xX3IzusvdYxs5df6ow5JYrnqL+R7y21sVLyuz0zZvv1SilaAK6X+TPJnfl5/pSi7rKTYSu8ohWNLEYpxWLPwtTpOPsniqUu1FThaSu++O4Tbczfnr/n7WhTqWv8V5AACv4H2gsZAsOAVdz6CEeu9PNhE9UjSU8QHfG6Hp567mmOgvI7zmZWt1hmxD5nePcq7fSXscSZxxyzOWvY1YvnCL2bzTTHjtYjOjhPZ5JFniwwcl9n89CvTZkJvhjLkY87Q3WBYdeYQ4wdT3BiS8ZUNUnRdoERrQQkifP8/O3CWWCPaMp5Xsez2F8BACGKKmmr0/uRO0ZoYyJDTRCsufEsJbJqTFLXMgyBYcAqnq39C9GhDszgTgxmptmo1ihtjjHtOM9ifTybdp6o/xPR317izOSUfN1UsT7yM6IdBTOwNdlsIohwYhPUy08PI7rYnTUOAPAMZ5+E81QOcR08metXFbzP5rt6FhzYOX9zZFAyPwpc+ohfS79e2hRWEf7QIaJLunIor9NRZsw/XeY1RU9XFgWg0RdJyshJ7ZyysHzLV9rY6mx+Xl82rqedUxNgyvQ2BIYBA3gs/Bai7etwOKuszfbvCE+9hNyn57nBmD1HvGLijUoG9TrOC6i1X3GUALjPk6Oams7+m+g/3mMzsfdqztpOGcjM3zFLdzgE3s+2/PR7mfG291Ya/kxhen+yzgB9n1CK+OVw0ED4bTxHRfbc4gbOc3E8yBnoUukFItpxKK//Pl3bmrie63u9Fa77JMrCH/l6CO2qy22VEduTFK8FSAgUGyYpAwaAN0+xuaKeA++at+Wwb+Bwjs4kY5dzXkWxojDY72etpXgAm2GKh+ox/JdWsbkneks7oh0Vn7jd3WyWKbHnHWGer75D9FKaMHmfYHPSoZfZ5p76FB93/UZPukvhRwG7h7lE+i/72IdRkWS2xA78PB1bcqOnzLvZFPTlDZ8TPXr6M9qco75mO1d9G8NqzxboPqJjl5RouRoqMADDJGXAAADg9duVsuHfcobvH+9z3wlLfSVCVnGdJ/tazEjTuzB39/7Sesc9/3xOaLP3Ztt/XgcuHNh9xi6idz/AZhdRou+q3dZzjaucB1jTuTiQr2mXx+e7Wcjmqv8UV83dv4MbTnlXgRk/R3FBuCUyM4usw4zZUbDG8dnED7Q5x4x/VhuzBfszw7Qxx5/0opI1ERICBZbqwFxnMASGAasoUbrd7fyEd/KZSlBPcUM9Z2L821zldeJS7rg3oO8Oor/eydtwez+96qtM4HV168QJhjd5cbrOR8u5VPmyjcwUv0jWS2s/qfhvtm1l5j73BJvjAt5njST6Hn3XGWnPzvgX+nIux8zl/bTP2ApHxZe84DnuLjjuea6y+8oOdpTIAj1k1j2jcgnE/Wrt18YeGsfa61tzbDNzXSswJe5d/yYpIw/DgFWc+ZCrq/oe5hfDNUXpje2kM0mHPP6duVxms1ZmfWb+BR58De8YvXqqKOE5PSezxpH/JDPz/LpspnE5ozhx8/VrfL+HcyISi1kYvhHPvpkCJWM9ZoberjbPl3eiyUpBw4ivFEf5bxZKlliBmjtz/qEwou0UeZDVgLUrUaR/h25xPFb3PdtMUoVb6mtjwR5KFd2OesOpfwOVzcNo3NJFzl2j358l3NbgpJGHYeD6Rcg2pa/Ez+xcTruHTTu+v+tZwkXnLvCAEufvvZ8FxskZ7FyuN0T5PIDD+9iM5bmcM8yDCtk3kB7O5qL3P+Xig9klfBwAsiQLkSeasoA49To7i+1zmak66PICec3Z47+/G0eZrb4ljOhlTfQQV2uIu4/nCJrBCXKnFiv5Dor93f1v/VnUOla5fngF8/X7SD6iRqb9NwKjspBSoFhe/xqGITAMWIXzWjZFxI5jc9Gh8XOJHnlBN+389NcNREfNYh9F8WEu8xG5mI+f7aaXG/FUdsne55ihHR3PfpLIx5hpvnCCnbiOBziSCABgz0xAeDKDC9rBgu/iAF6U/XGl9hQA56Ps8R8Q8ADRuYu4jIpnBXrh+J7gdZx/nbVED7beISuCn53orDPuwP4clpyyXjulTHz+/vvamJrsNq5+5dvT/lcoMcJqDRgAYr/nEM36L3LYZ49DjxBd6KH/rJzasRnm9GCOcPI5xYzCTtnMhiyP0RfmoRTIK2bm3TiewzhFG3a2qC7uk5/opUbszjNzb/gNM9Ki0WzWct7Fvpa8Jro/p18zzpFYdYTt9k6NmfFYaPJqFfMXfMjXXMBVdHPasJYzvAVHYp3L1RMOozM4nNrRxoTCwVP08vG1l/ypjFgvZHktwpSHYWgYBgwgZBq/CHEz2Vwxpzkn0A3ZwJnfADC863aib7PQTIfm2MFVTC/frld9LcljIdS5OUcfOQhm5pcUIZUTway4WyM2tQFAv/YHiJ71I2sDCZdZg5CN2P/gv0VPBrz9Zr7O3TezAGl7C2tXgybbvuuOcOQ8lukjviD6Dlf2HezM4/s4l8saCQCkf8+aT23oGf1lIbWZ7i/Ne5bzMILerZkd9yQECuX1z06v/zs0UGlcmMgvets6rGH8mRdGdIMmug9jf6qSq5HBzOdChlL87qTCaNtmaHM+feN2olOLWKNYtOlWogNbs06RNpBDiSz1gHj+Z84GL+nPc9gp4aoz+nI284Wb1KZDwJjVrJF5n+I5fpvEEU0VQZcxjxEdfx/7YgY0ZUHYwo19RCUW7PHu/ZXvda52SplwTdBNNkG/6vk1NRXFRh6GAQNA6GAuTHfoad4VnrqVGW3mLp3x5tblSKDA35Q6UEodp/onlcS9aaw9AMAaL/ZrXByltEcN4Tm//pBt6H/ms9aihvoCgHRnYdloKTO4i7exlvLxKM5ZcYrRzTaNldztU4+z8IxawyGvFUnci+V25nA7zKa1X9exBrF3JAv0xK1cDBIA0I61ElsbRz81Yo02NseeQ52Dd2un1AgYmd4GDJjR5Dd2KKR9wEw08QI7l4c/sF2b4/M/2BF+uTW/XG4JSv8L3oRjTCOl2iqA7/p2IbreWo6KSmvNtaSGbRlHtIsS7RXYVc89ONeH11nsyhVbQ2ZwQqKdkhme314Pk0pqwef4nJIKrYf32grPIDZr1V3Ez/diDxZ0vo+wcA2qr9RuAXC6ceVKj3/3zJ3aWPCWmmmCsgRLWtn1BkNgGLCK4zfzz8RNYawNVvD5u1/UNYzG2ZxLYKe0PlW7zJXMZOfnmx9xrSkAkBNZkDV5h81WtZ+MIfrwUY6TDytijcR5q57vEHVYuZcCZuYFnTms9szdLAwCd+vZ45k38L3JYpXRsD/Bdwdshp2dklehyEKlIy7S27FG4ZqkCy375MqVN4/prbMb32dZoNbqrffhqAkwnN4GDJhx4Tk2QQkrvSgDC/UQWKcTbIYprseM2P4SM/vsNmwuss/RX8bgjYrNWKkNVTiG/SJNirisx+nJSgTUduulHWQA+yRCprO5rrFSZ2tfczb1AEBJLK9LOPEDTWrP5jvfRVaXpUHtwyFKuLRKnQOsQZwezqxgVDsuowIAhUrpi13QczXKQpNW57WxewK4IOQPKF8TomsNEsLwYRgwAADBv/J2NDuId8BxvXmnX/SXznhTe3PIag4nIqP+WiU5cK3SFvZshDZnejO1DpEXUSlPsL/B53M+v1lQDNEuP+nVVPf9rmgla1kgJI5QorccWLC5hPGaAODGCdFEHzzPu/uqYDy5d7AAfmz8OqLr2LPJakrMPUSveYcDBgAgLYLXFWpj8cGLq8O0sbngsUAb57xWICWMKCkDBgDg1Eg2RQQpHdUbP3nU6hwuEWFElxziz2T35/4XJeGc6Oe9h8t+AIDvHjbtnFaYt98KFhBOaczsC4azhlGwSDfDhK3jzyS1YGE5ZPR2oufsZG/zpz109cDLjte91ZtNY42cOUHuM4WplgdNA3iOts4c2XZBiSh7qwHXsxpxL0eHAcDdYZxc+fdk29ZU0ClTGysquF5YkDAS9wwYAACRxRrDsDc52mWRS1+iCzz1F6dQ2bzXy2aNw66IHb/2uaxx5EbpZSUywliQFYUyI94w7GOiu+zh3A77nazmOHClEABA2r1KGG1tNuWsnMId9jwf5J39qQI9f2T20e5EO+5gLSRgj9ri1vZaUqmT2RR26z3jibbLY02oxImff+Bves5E8BuLjg3JAAAgAElEQVRs0vsbuvZUFqa0XKuNpRWzL2sFbO9ffi1AAkZpEAMGAOD2jmxnXtWLq9X6lbB/ot9GjhwCgKFebL8+/TR7YV+Mvo/odxusJLrP6ue0Oe0VvnpXFCcDLsloQnT3UPY3RK/nvh0FgToDDFzPkVcFDZihxTzOGkjglxx9tHYHrwEA6rZkx3haIz5+ajAfj6iAlSbmHhamAUp1kUt9WfCNb7OF6MI7dbPiT0nqvSRq55SFd98ZrI01ekR1ctvejvZaQVU5vYUQiwD0BnBJStncPDYDwD0ACgCcATBCSplmPjYRwEgAxQCekVJuqpKFWIAhMAxYxbbtXLqieDzvuqNmcm+FcCdmsoDePjXKken1kYqdC7zz3NVfr0P0ViLb2bee5YibDTkcwQRl02yvmNo8zuuakWsQ+xfUKq8lSsFb8QQz0exCvdxIXCd+7QLbckLcZ5HfED1hDJvrygOHTGZe9gXsSA/9nJ//l+vYh/HzjNnanE/5xBB9F9po55SF3F7p2tjNvizE10FvOFUTICGqsoHSYgBzAHxZamwLgIlSyiIhxDsAJgKYIIRoCmAQgGYAggBsFUJESimLUQ0wBIYBq5h336dEH8pjc8cn3t2InnmBzTQA8G531kKSH+cChvmK/zo7lH/vfmF6MTzfGWzn8mzE/gW17Lrr05zL4daQuf/fdjpzd8piJpA+gO3wTZ7jdcksdrSff0w3pUW8z72wk+9ih/4bj3Iym607eQBwUK1ao1mISyXUzXs8P8ub39I77hUryfcBNjqoo/z1jcRgL+6Tsg43a+fUBEhUndNbSvmrECJMGdtcitwNYID5/30BfCOlzAcQLYQ4DaAdUIGKleWAITAMWIXqpL3JlaN85ndgc9OcWL0PyanpzAgKazGzdr3Iu/0mL7F5Ka+T0qUJwDMLlxI98S9uPNQymGsdxQ9nn8URJax2xwPvadd45PWeRNfexDWazg1nIWPfngVI6MNszgMA6cT36nmOzVontnBYcr0KCIyiNlz2JHMNC648Jb3EuxWrX+4J+gY1qWXlOspldNHNTUP91WZRulCpGRAoLr/Tu7YQYl8peoGUcoENF3sEwBWPWzBMAuQKLsL2JPxywxAYBqxi/FguVeH2IxfLky14h3z5Jr2+ar3THIHkHM879Zj+nN+QeQdHDnntYiEFAB8PupfogYu4o1thCTO4weu47sTrL4zkNXXWbdCnP2Pm7XCMd+L1NrFGkZDHORYn5+vFB+3seXfv6sICucEoFsAV6UIxLIrLiTx9MwuuY4X8bGYlcHTX7h26gHbTS4TZBIcgPQAg9v4wogM+rJkCQ8KmTO+kijZQEkK8AtNP4usrQ1dZTrXAEBgGrCLxIXaQDn6TI4HSizhENqVAz2f45QQLlahXlRyJ02y7VgVE5s0NtDk997OJ6fcOSkKcUqZjW0+Oknps2g+8xlw9Q93+BN9LQSQ/i+Q4Ph68hXfRYrkeSpp8G+d25PjzO39iItvnGo5THCXlQBu3GKJHxPQmelo9DqPtU4sjsXbVUTzxAELm8TpsFWQFEbrAyGjD2pXeiLfmwAYNo0IQQgyDyRl+m/ynVepFAKWjN0IAVEFXeMswBIYBqyg5w2aYrT7sXPZ1ZiZ6crvuC4hcrzDOQmY3Xt+xdpDej/MwcuvouzcPL2bWJUGspSR2UJi9ovgsHXU30Wq2MwBEdOHdvl1vzpg+Na0l0ZlhvIbw5RbMOIOY8RbuZUHllF55xvPkxuFERzXnarTL0m8i+mk/fv7+AbqDOrYPm/QCPoyxaU2p4/XKtN6b9Wq+NRFSimqtJSWEuBPABADdpJSlPVRrACwVQsyEyekdAaByzdfLgCEwDFiFaxTb5YPdWcMYHMCmnuV36E1wEreyhlDQkkNanVKZaSYNZK+twwEWWgBw5kFmNur72vB1zhY/+xqXOJmwmEuRv9uImT8AyM4cIVbQiXfJjul80fB5Z4m+8JAuPPO5sjjcOMgMPtGVa4UKAH6HWVDFnQwj+paxHJXWdiM7uYWDbtUYOIyd3Ac+1E4pEyUbdOEQsIc3EtVmS/kXUFV5GEKIZQC6w+TruAhgMkxRUc4AtgghAGC3lPIJKeURIcRyAEdhUvqeqq4IKcAQGAbKgUlNNxC96H7emU/ox1nBeaF61Vf/BvxTc8xiO35OAAuE4E95jkIP/R14+b3FRN/pyn6SB7ty2O17AdxE6LWXlUS+fjq7Wj5rJtEju3AugezK/sXcL9kMtjJyhjZnr2XceS5oC9vtE7tzrkdF9uBdHuNN5roTHGL83DvcnrbpJjbvFcXozZEON1KF31ntnLJQ97uT2tjle9hU6btPO6VGwNRAqXJBAf8/l5QPWhheWMb5UwFMrZKLW4EhMAxYhYuSfHB6CNvYw9uy2ebBID1xr/Bmfpnc7JQifVmsgaw9zgzui44c2gsAT8xjZ/yTUSww2keyH2RRXGei7Udy9NHFM3qW8bAH+Bol9cveRbqO4QioPoP1tqQOSvHGMw/zdR11t4fN2LiBkyvD17M5KElJoZBJ7Huxc1Pa3wIoiUvQxmzB+RGR2lihYiasmVkYV5zeRmkQAwbw7A7e8EzstZro9w9zhM1hH735jtq3+uZGvDs9P539Io13M7Of5smZ4ADQZDE3VTqwjyOakl5h57LddjZRuXbgNTXw1LWYjHAOvW369GGic1LZTXu2hOkGM/Sw2hPTmhNd66BQaDb5WSkObBHhS9lPUuzN95GntOzO78hZ3NH36rvlht/yxkF9ntZweKzeou9UIYf/jnm1k3ZOTYFR3tyAAQBRbzLzWenHeRZ1QpkZ/RSkZyY7KQFIJ7dx2Ka9G5uD8lpycqDfZN1EMjSQberRjZgL5vlx1FSikguiGsyL/HVTmt9uZuZTg7jqwoAZzxMdfpB34fFDdb9Ik3diiE7vyPfqM5vnSKkADx2xejPR2SV6eG9pRIzkaz455yntnNhuLETqbbdtTT2GPKKN2Wer/ppD2jk1AVWc6X3NwhAYBqwi8XaOjunxBDNqfyfeEX+8Tc/0rvsr08ktlBatmfyySXs27ZzfrYd55sxje0ZAJhfHk4qQSm3M3oCPHuVcqcWX2GQFABlK9viQh58mOv9FjpqS4/mVqnVUDwBI7caaz4hXWWPbncG+Ar5C+bBwMJf6sL/EzybxLhZSSsoKCi04Tmp3qlwihuPOI9rY6alsGwuvoS1aAaDE0DAMGACWvcyO2/7z2C7vlshbdXGbziQTlIZ5Id+wQOg57Rei+3hxXsDsRD17fHMd1lLuiOJd8q4VLGRKlH4/77VgO7+dv9pfA4ADc1KnAnZQe87iDOqMZny+9z6dyTodYX/CimUceZXfuynRzhWIkoyfxOa1F6O4UkQzZy7u6GnHO/3duXrjp6lLHiC6HvRkyrJQslFvjuSclmHhzJoHKaumj8m1DkNgGLCKpxt2J7p+QAzRWTdyiGzj8TqTPDaN/RrOyeyg3jGQ/QmfD2UB4XfMQsBlB7bu/5HATM7/Du6h4e+q1HnqyMze50ndh1HizhpG+k3soPY6w3NeepzXFHeX7ki3T+XXrs5+vrc7J7I6tmutbZ3tACD7LAu/yFbs4A+0ZwHx3HmuX5X8ki4wHCuUm/wPziboiZEz2y0neg70Rlk1ARICRaqadh3CEBgGrCLxSU7y8j3Btv6cx9jcMf59vQn1pL+5Z0bseN6NhUxjG3ujr5WM6US9DpHnOWZq6Q05xqYkjk1WJ8JZq/n15Q+IdtupM+YXEphL/vUcCzb7WF6X/c4wokNP6TkVxY5KaZDLLDz/uF0NGrC9XMaTPdnX4ihYGL4SxzWyDv3IQQd2z+mJewMbcbmRXR/YJsgcTrlqY6da1eTcbkZ1Z3pfCzAEhgGrCFrHlQaarOCsYZUZ+dlx5AsAFO1lB3RuKDPSU88yfeiWz4jOh7777zqXmwKNfZjLXbyzkYXU6UEcpTM7jcM8w5z0EhxrtrLZ6rYZ7JTddpoZLcBZ75A6k7S7hb0ShQ58b/kb2JTmP8d2gbHtMkc9LVx6J1+jjhJ71YTX/XlrLuwIAM98yLkb/jZWq22wQvfGzAnlXJkI1MxEDCOs1oABM5Ln8M581U7WOFSH9lZfvUR12BZODEu8nRPeXFLZYXjPN+xcdonXy0r4TGMzy8zlXPl0xTDWII4UMpOc/1UvoiPuOqNdY9sg9t903ziO6GO9PyH6pYQORK+L42cFAI77OZorYCMnXuT76dFatuLCujCic8NZKEU8ydqCuJHzXqYsGKHN6Vy/cnnYl27208bGtPuR6E02dvG7dlC9pUGuFRgCw4BVJJ3kkJnQLUr71Fpsu+3wONclAoA3J28nemceO0AP5CgF+YrZ3PGGv77zbLGDq83aKe/rgFXPEr26PwsQj1hmgEkfhmnXGH1kONENA1kTapHCJTXCX+WkxYatdW2r4G129NrfwgLirdD1RE8Jt61REQAE/cpCqOgAP0+H8DD+QIqSLWinVxx2j68cQ0xppwvCCKV/ec0VGDB6ehswAABP3s4x/Y538G511kE2KxyYoTO4PoWtiU7ooDCfEI6sui+KCy5FbX5CmzNwC2s+Pus5bLMkgn0cY1dxbkHuWLbTb73xc+0aLdeyQPA+xsIxYC8/izNfclKe/RndJNXdmzWZXzbysxkewfkK4VCKT5UDl29ghq82P3IM41Dp2gf4WSS10gWG3/FcbcwWuJx31MYmfcxlZQJtNHNdK5BSL6d/PcIQGAasYlUsO3qTdnF0kTNX+UAGKwsAALdLvJv3iubdWJ0VvHP/uWlHoj199d2bKFa6xm3gn/NNPruIHub9N9E9D7LZRRUOAPBid97tzz/C0USOozmU9/kQ1q42hXBfDwD46ScWEEUhvPN2OaELGVtR+34u19LOjxMf7/TiDPSpfTmbX1hIL7+g+JlC9diGMtH+rr+1sQPfNLdwZs2DkbhnwIAZvYL4RV93WAkt/Z2ZUeAq3Qzzyy62kTd+mwvRpfZgR2+KUvG21o8c3goACYqrxCOfGe2yGI5w+moTO359TjKjdreQDL181V1EB2xjk9OpSC7DvnYoM8Dox/RqtRtHsl+kRHENzEvuSvTfb+jrsoa8j9hHtKQ3C/l9z7IGIQtY66kTzSHJAOC3UP9ebcG8ej9pY2ueZMH1mQWzYE2BYZIyYADAmjc5J8LrjJJs5cWVZhOG6qXIG6cxQ8pvyWqI7y5mUNKeczsy79UTvEQu2+UvbmMTlHs8c2Lvy7xDdnqR80XsB+iO9eJUJbzUns0OjV8/QXRSH45OCtqpqF8AHt33HNFpDfk1fHz0GqL/hp7wZg0pjXmd3Vtw7+w9X7Egc3Vm4RnsqYfVOgvuL5LT1bbWsX2C9UQOvcih2oy8ZsCIkjJgwIysIPY3eB9V7BVKtu6psXoZDzRgZlN/LjPzmCHM7JUCuciJ04WQcyAzl6EPbiH6eBbvqrcf5hDYqGd4TTJMz/Q++xILgOI6vLDwJXwffss57FaonngAUHp62z3Hu/3PZrLZq1YF7PqWCv0RlLy8nsc5Yiz1Q92u6PrDH9qYLXCoq3fcgzOrdSUWyqrXFBhRUgYMAChSNoEiiRsqxd/HAmJcP94hA0A/D96Jr2rFzHv2EmaSOS3YJNVkmgVzSCJ3Hpr/AZtyOkWyVvNEx+1Ef+7NfpIRTbl8BgCI+9mBL7JZSOUtYiaRsIr9E84ZeijqpY4scO0P8hzueodbm/HERb637TH8HbUJ5jDnSWHsq3lu5EBtzhx/tgHWmm+bIItcp+e57E9mTdJVL0NWMyANH4YBAwCA+h+xnfnka+yPKHFhBvjRcmb+ALD0IEcTJTVnk0leI852rluHTSIJt+gZwdte5oZINy/kyrEp09iU82sK/9wberMQWtmVy7QDQJ25vOO9P5Cd2lN+5yJ/TX7h5LSB3/2szVnLga+7NJFzN/bG8Pa/IrnQn4Swwz86kLWvAQc5JPl7F/bF9AhmAQ8AKxuyELK1sdOjtXUv+bBaLCwnQK90XBMgARQZGoYBA0DLHczgzkezCSr0FTbTnH1Dd1DHRyhF+X5hB7V0YH9EfAGzo1q99N3pt5lcd6jbXVywsE5fzi34dmMXorvcyoJQjtTDPr+csoLodt9z4l7DlSzoig8fJ/qLsbrwvHQDX6fpXdzXQ8Tqz89WXCxif4yfHT//Wc2/5TU48fljL3BXRQCQTpVL3Hvm1APaWPRxNhtGoHJmr/8Khg/DgAEzfv6Id5ZRj54mOmQJ15KK/ZrLaQBA8ELOkZCFLGRqtWTmn9CR7TL29XRmNWNLb6IjF7Ig26z4KEJuYCf37+u4V4Xve7oTd0cuR4SFKwIiO4gFnetW1g7yFuuv2KxHuKz6lgyOrDobXfkCfHfNf5HosEXcsOrY21yvyv4SC7HQtmyyAoAuHY4SHaedUTZcHtKrGItXrh8mawgMAwYALHqD+1o3clB6VRRxJNDmTuwoBoBFz3Ff8FOF7MTOlmzqeWbvIKId1uqVThuc5OsW+fDOXE0UK/iFna4OSiHZnA268ecNDOV1juad+Nft5xD9+jnWKDxG6SXAR+3kZDUnVxae08d9SfTcuRaCCKxg6nCew/MRZtadXZhu9t0YosM89bpPQc5sJoyDbYlqx1/RQ4wbL+DNhl4xrGbAyMMwYMCMvjufJLpZCO/UO/gxUwz+UDftdL2be2hIR0VjUAKvRAm/fIE/6QX4vBYxs/kglJ3tu/NYAHwRz63rxoVwBvuINY9r14j8jK+RnsgVcV9czAX5vljA5UeWZehZ74ePsxYS/h4LoXGjHyK6ImaaucMGEG13kPNeXu3HyZg9xnJ019SgrdqcdkqewSB01M4pC2v6fqCNvdeOq+bGddBOqTH4X8jDEFJWzi55rcJL+Mn2Qm+6Y8B2lHRrS/QzC9n+fYsL70Y97HQb/A7FGnGmgLf3y5QifZ804mvcP50FDgB4xPN+1DmFd+pnBinhq15sTvJfz+tM6GwhvdmNr9Hwcz7nzGClt8XvvOv2Oa7nFcR2Z1OZdzTP6fMTm/yKLuv+G2tI38BmrcvJHLrrdpjvvUiJzLLjRwUAyA0uu4ChNWQ+qEsDac9M1muJHqn2b2CrXLFfSlnhjh/ejQNkhwUPWj8RwObuH5V5LSHEIgC9AVySUjY3j/kB+BZAGIAYAAOllKnmYxMBjIRJQXtGSrnJwrRVgmrTMK5y060AzAPgAdNNPySlzDAfawlgPgAvmPabN0kp84QQNwBYDMAVwAYAz8rrVcrVEOxSnM3rUzju1tLu9ExBA6LfX3Qf0VkRzOx7/sURTw0HcKkLAMgtYoEQ/xM7UIUHS6lQfxZsTcdxWY/+vnrRxAnTHyM6cwKH8ka+zozXcybPqfYdB4AmjqwtnS/i/I9RmziCKeJJ2wWG5zuKgHBhQeaylxP5ctux2cspTZcYhd62N3IqjeCn9GrABcVKtNySSl3iP4MEUFRSZVFSiwHMAVDarvgSgG1SyulCiJfM9AQhRFMAgwA0AxAEYKsQIlJKWS3Wveo0SS2GftOfARgvpfxFCPEIgBcAvCqEcACwBMAQKeUhIUQtAFc4yFwAjwHYDZPAuBPAxmpctwEFTqeZCR5oyy+GQwPWFtqN5UgiAHBJ5M/4/8UCIuhdLrmR8TCbO+RK3YfhkcbF8OyVHhsuG1lgpLXggnunzrJ5aUxPDi0FgOLb2FxUx461gdjuzJizj3D0V1tvXdA9t4b9IuE/8DqjLnFTJr0Fk3V8tWQ20ckKM/smlQMT3vDnsh2FFvhNh/0PE+1v4z4267ZMbaykwIIqUwNRlT4MKeWvQogwZbgvgO7m/38BYDuACebxb6SU+QCihRCnAbQDUC2qWrUJjKvcdGMAV7onbAGwCcCrAO4A8JeU8pD5s8kAIISoC8BLSvm7mf4SQD8YAuNfxf3beOf9bduGRJck8I7Z7SIzZgAoac8RTLEteGcZ6Mrx987pzLBc39LbvvYP4Cqu01ax1tJwOb/AvgdZw6j1KUdFRTjq1Vi3nuUEQ7e3lWzwiTxHwAxuFLV5TTdtznrFLAIyGijmoQlsxqrFwWDlwtDB7MR+9YvFRMfkcthyo9VcDdjtnO7QtquI5Cr9+boWMkr8OPih5M8j+jk1BLL8AqO2EKJ0vf4FUsoFVz3bhAApZbzpOjJeCHFllxYM02b6Ci6ax6oF/7bT+28AfQCsBnA/gCtpnpEApBBiE4A6MEnMd2G68dJFhsp8GEKIx2DSRuACtUaNgYriUDZn4wb9zLvVDl4csplZwhoJAHw7g52bl9vzTt3zWe7i17U25yZcKtD7JHw2qT/RHkrU05SVnNi3Jp0d0F//xpnLUVN5DQAQ7suRWIW1+ZXxcOLjWS9wFrz3czrjlee4blbK45wd7vaN3mjIVtjv58S7Ca+wQz+5OX+HTZaw2Uuk6bW7ioNZy7PVLlwU6KONnRnI72nDP22c9BqCDU7vpMr4SxRYumi1mez/bYHxCIBZQojXAKwBcEUfdQDQGcBNMFUf2yaE2A9A/9WW8TDMUnoBYHJ6V+G6/6ex+nc21fgeYSaYNIh3iVlvqT2pgbzm/LsOW8UCo+RTNg+taHM70X7HdOdx3BD+im9owUJmYybnWXx9gLWYO248TPT+Bfq6U08w83aL4/so2cwRT3V/Z/NSibveNOjk+5wp76DU+St2VhzB2gzWceY1FkIje3Omd34J+38mDOckxtNF+rpzlM9MamAbz7t0o17zxDPGpimuWUhZ7XkYiUKIumbtoi7+afR+Ef9svAEgBLanyJQb/6rAkFIeh8n8BCFEJIArFc8uAvhFSplkPrYBQFuY/Bql3+JqfRgGLKPxp2x7jh7AO8X4L9ihndJLl9VOHJ0Kx3S2XefVZWbyykucRzBuK4eaAoBLAguug9nsuP07m53z/e/hqJ6jI9jc5NREZ811lAhhz3Nstjr9KO/UPe7mhLes13SFOOILFir2x2L4BHuesyq8l709WSA0c+SdfdP53AtEttD9DXs7fKaN2YIFz8/Sxt6N5fLxGfopNQQCxVXn9LaENQCGAZhu/nd1qfGlQoiZMDm9IwDsqa5F/KsCQwjhL6W8JISwAzAJpogpwOTLeFEI4QaT1tENwAdmaZophOgA4A8AQwHMtjS3gepDxjtsdvGfzyxMNSflpLC2AADF6Wz7z53CyuPFOLapv/QNO4YdLbyLRUqfIRHCWki3cNY4dsxhR2/aIN4ROukVvRHyHjvjhVJptr47N0g62oX9O7Xq6qG6aiVeryxOKBSFiohIUaRtORDxAYfmvqiEfB6bwPY7XyXNpWiX3nHvJjxKdCgOa+eUhaFf6g2qwhfqQQE1FTb4MMqEEGIZTA7u2kKIiwAmwyQolgshRgI4D5NJH1LKI0KI5QCOwhQf8VR1RUgB1RtWa+mmPYQQV/pkfg/gcwCQUqaaJeRemExOG6SUV8pnjsY/YbUbYTi8/3X4jOBIoWOvclE/1y1hRNdfzVE+AFDkyzvzM77sAA35mbWSRMXaYcnhOuX+pUT/mMolNuJzWWNIbsXXcArmGllBH+mvw4UJnB/if5AXMmnWIqKfWsphuIkWEtFc6/Hu3XEhC9N8L5aOnieY+ZcHeS3ZVOZyMIboxvNZo1M1PIds/YHHu+hCxBbUPqTzsaw2rIG5nNcbN9UEVGUtKSnl1RI6LCaWSSmnAphaJRe3guqMkrraTX90lfOXwGSCUsf3Abg++jjWUAz/ZTfR93mwzI4u4p39qO6DtTkyv2ENwkUp2+Scymaa8JXMzC/01HtVTNhxP9ENv2KGFN+Ro4/6DGBNfeu3rHFIe73WkdqEyeVn7j44ceooXsOv7PAvOs0BAQDgEKB455WSGyKK80kqgh4zuTLsxvimRLu8qRR7HMb3Pr3N99qc3Vw5B2XQ27Zlertf1P1Q9hd4zkoGYv13kCY/xvUOozSIAatwFPwa37CP9wL+Hko12yM6w1vxKpeFeOK1sUSXTOIonVkRy3jOIl1g5EhuvhPTkTWfjzazffzQyxwl5RzG8xV66q+D73HWjOKeYGdyZriSCT6Smf+hH5U+sgDqr+FIqhJXNnNdast0cAXydr9dyJvRoaN+JHppU25X67iPSCx4hO8TAN56iIMIbG3sNOSrDdqYlyKk5zSqfOHF/wr/C6VBDIFhwCrSitlcMb4xR9ykFHOU1Je5ekjzy82ZgXm0Z6e344ts7hjQjUuB5FnoUtpgMvsX7BtylzjvW9m0kxPAL3Se4moJWKPb02Pm8IWbBrBfZP9pvuakkHVEp43Qs97Hxj1FtBoBpkZaVQR1f2ahtOHwLUTnKHJs/gju0Jf9uJ7VPT9WadE637Y1vXmolzZW72NmQXaomXG1ElXnw7iWYQgMA1Yx53R3or9tyXZ7RyXSOaauzt0fObKT6Ht+5oq2r3VgxvppTGeie/jHaHOe78Yhr9lPlF09NSuIX+jb72cTlccgvf+2k9IxL30ih96K+/kVGj2WNSenDN3I4uLNWkl2EJvOil0qH1YbOJ99Ab/siyLahyuVY9jP7NAOC9WLPQa4sSZpa/dtcUxvsxv9GGtwDX+xcdJrBgLFJYbAMGAAnnMVc5DSLvrZGM6wPnRM7we9NpVzIDyb8A54xhHOu1AdiJt/0juxhfdi/0D44hii5/kvJLrn59wjYssK9mGEbNNDSZPasOaTx3IML9zOFXJju7LaclyJgAKAA/s4kqrHzRzyeuZFvTy8rYjrwFFoEYK1MTtXFlJ1V3OotPTTxVSanZpHoQc3lIUH+2/Xxlq6slY3F7aXcr9WYGgYBgwAmDL7U6JVjaKND4fVxu/S+x6ksM8V+YeZQfkd4zm10iDndOaUN4N3wSe6st199Dlm1vXqKTkUg9nsEjxb711xk+KQ3vkShz1t68U79+Nr2AafHaXXSop6mwXd+Tz2WTT/iQXIUb3EleZ5RykAACAASURBVFXYe/BuvvDGSKKjn+Dn274+t6L1d9abSR3P4Mi2ku62rWnH03rI2OqmPEltG/0i1wqkNASGAQMAgBde5zpD9vnM3L23cBkKnwzFgwoAgzlO1iOWGWlGfWbeqY35pxl+TDcXOdTnkiXH7+U53GI5ZLPwRjapeO9kE9av2ZyBDQDSge+1Qa6y7i4syLIXsXblkKD7Ao5N55BXFLCv5fhBvmYEWDsoD2Qkr8PtDc53DfiIj5+TnMT4ezed+W26932in0Qn7ZyycGa4nkwzscMqolfM9dfOqSkwGigZMAAgYCTvvLvV5mY8Pd/lgnG17PRktWy5negPLt1K9E9reRu9c+R7fI0bRmhzun/KWkqTTzjS6tjLbErb0oFtaU9PeYTorAi91lHsLcwEln/9MdHvXObQ0uRYNmv1aM3CFACW/8m5HY3nc56LmrhXkSys7PqsYcgXmFlntWNhOflZzqyftJATJwHg9VjVaZ2mnVMWXE/rwvPr77myonP1JSlXO4ywWgMGAMxssJLohCKOgtqRw2aYd3ZyOCsANPyGhYhTIu/2w+K4P8NQJdM7ZxDncQCA1xguw5EpmSm2cmHH7+pM1iB6fMs79+8n36Fd4/4unIMyMY6jvbbHsM19vSKU3orTn0VUQ153Tl12pLvu1IWMrbio3IrwYLNXveXsjH/hhyFER36j9/T+O4vtiv42mo+8o3WO2uZ1joo6utamKa8ZSAiUVG9pkGsChsAwYBVPhXcn2j6MTUGZnzAjCKqn94NOeErpwreXBcDkEZxsMGHrA0Tf0JLDWQFgdv1V2lhpdP5uPNEFauvYYhZiXumclAcAy+9hzSfqJWakHr15J//MU32Jzu6g+3Ncz7OwhFLzMOtWdnq7/mB7i1a1XtWDizkPo97NbEobO58z1I+9oUe6OVSyipvPX/rvoqcPlxc5iijtnJqC/wEFwxAYBqzD9xeOmEl/nBmv5wg2w8TdF6bNUXAzm10cFOvE4m7sEG3sz3PmZOqVTgc1e47o8/2U0h8hfM2jr3JCoUM6//wjP+WsYwCIekNxttspjaCWszmuuDH7J9wP6Vz25DMscEWhEkYbw+crJbPKhfw6nNS4fEB3ok9PYgHue4mfnSjhzwOV74dx5kFdS/ygYc0VEATD6W3AgAl//MVml7APuPxFySzeIhdaaEUiz/NgwD7mPsWXmTEn380789oH9JfR/U82Odl3URzOvkojIkXzSYnjqJ+8j/WIJoceHPZp58yM9uxr3O+8wF/pJLhVL5nuv1cNGjhO9OkXmYnqbNY67Ar4Grn12Z/T6A32P1zspdSWspAo2ak7a2BxH9q2JgcLiRv23ryu4nQLFSBrCv4HVAxDYBiwDjdm7q7DlXKrRTFEFnjoZpjiBGb4RUrPh8yR7Ai+b8w2ohtZCPM8nc8Mv73kdSzZ0J3oS3asKTm3YOYUt0tn7gMOssN/97Nsxy8M4egtewc2c/kc1CvNxgzgSKCE3lwe3tW98kxTLS9SfxEXMDwxkb8jKXjd0ktXJ3b+ypV5w23sAhr6wUFt7PwzHAodPL1mhtUChoZhwAAA4JPOXBPyqencm8LrD951B33LUVQAEPcg5wG4JSgMyo6jdr4+xWG47YM5TwAA6ruyxtDGjc/Z1JqFTHwCJ9W5rufdbfrtbMICgO1TuIZGQQOlkuwhZhIvPv4t0UuL9RatIVOZKebcy0mJrpcrv1UtUHItL85nlcHuON9HVEfODUnN1w1hN7TkfJuj2hll4/RnjbWxiNGKSc/GOa8lGFFSBgwAmHSsH9HOJ1lAtB/KO8fofdYLyKU24R2wSwq/bYGPM7M/00HJ/APwRzP++S4OZXbz1Z3ziC6JZCY5LIsdve+0XQ0VU3dz5V3vaL5GRgMWdK/v70N0vbm6ozc1l53aAfdypFBxJ47mqsi+tdEXSsvVQtYYchqxlphwnLWc8FG60G/lrgoM26rqRryt90xHqDLH4ZppkpISkEaUlAEDgP9D7LOIXsge69aebOd3teQLyGQzzG112G7/8Z/difbby3WiLvTQX0avM0w7X2bmPepzLvKXF8DM3t6P1/n5vdx3HACCU2OILglgLcW+KW/lndwUE9VU3QOR/iC/dvkT2RxX6M3CM7wC9ZXqLmZn+/bd3CGg0TJm3oUe/OwSpnH5EgCYchf7siJgW/TWyZF6nkvkQttyOa5lGBqGAQMAEMS+Ahzh+kqr3+xKdIG/7vV2+pUdppuLmZE2bqKYgy4ww/M6o3thszow09vT9ROib3/jeaLr7maBcXYg//xnrF+sXWOd0hf870zOHs/P4hDZ3kEcJjq3ty6EWjRhv8ijt3Hvij1Z7F/Yg7KLKlrCH7EcAOAazwLXIYmj0Fzc+VnEdtWjpBz0Uls2oVtHPWw5sgdrkj8116PhagwMgWHAAHBsHO+i7XLY/2AXxzWdcpvpJqm0RziaqM5e5j4inqOkilrybtZ/j+5fKLFn5uLdnU1lXgM5Z6LwE64t1asNm9IGHxipXcPDhTWGseHsjP/LlUNkb3XnBMSPnfRkwMI+LOje6s9Jin4HVLOMrd4CoN7DLJTsarPGBgd+9Z2PxxNdP1PXjPp98TPRP0y0EEpVBj6tt0Mba/LrcKIb4JBNc147EIbT24ABALDP5J9J2Ho25RQ05dyDPB/9xak/iO1Hdz/DO/Hz+cygsopZCA3244xrAHjyKPsXGq3imlfDOzOD2vQ438emn1iIiWJ93XVXsmBbiHuIzg5lbWpN/S5E1zutRxvFfsFayujGnN489wRrbIHsQioflKRE6c7rLPZl2iGazY72sXqxx9Wt1Sgy2/p29Dii30j9ebZrT9csDA3DgAFgWf/ZRGf05Z38rmyOgFq8i5kmALwRxP0wDuSwySRVKTfirGSJPTr7GW3OkPUsVBzbMsO/7U6OwNlqx1E6oVvY8SvtdIFxchRHC0VEsKksxJnpjCFsrmuwnHfuAPCgIvweO8BlORzs9VpctuL0VO4uWOLE3Ex6873X+56/j77TuEkWAEQooc22dsc7l6BrLXIABz9EbLdpymsHRuKeAQMmfJHMVUmDndlRufQ4h8BGLtQztOZ+wjvzpJvYeVx7L+crFPkwow7K0J2jd/7AtaA2PMzNKob8wv22/Xazs77LdC5019RNz8pu6sxmrWlduABfzC0sLPN6M9NI/DRIm/NXTy43EtCTr5G8gTWQiiByFkc0pXZi7cA9jgVIgeKP3pGsC4M5Ss8SW53eg1voVXdXL9E3FzUWhsAwYADYepZ35s67uX5SnRh2Jl++STczuPfnnXbeZn65UtqwAPE9xHZ8uzw98mrJu3cT3fQT1iju9eZktTOtOVLL15H9IqsS9T7W30zlazjnsmnN8xybZewe5sZFTh8qvgMAiTfxrjrmYh2iHW1zDViEWn7EVUmcbD+ee26owlIVlADQJJx9L4NH6/3Ky8KuCXo/jHoxrCVWsvrIf4sqNEkJIZ4D8Kh51sMARgBwA/AtgDAAMQAGSin1zNBqhCEwDFjF5g4cfQTlvb/vMJcJD3RXiusB6FqbiwdeGsyhuOdzmLEm5rBQirush2Q6HWcm+Md6Dh09dYxzN7x3xfAEzhwJFH+3vrPPVTrsObbish0Dhm4neuVZFjpOT+t5Bf7z+F7OhfJ9BO6pfPra7P6LyjzuqKTItXZmQfdWIvtRAOBouto98KJ2TlnIrquzm9w6LCy9TpzWzqkxqCKBIYQIBvAMgKZSylwhxHIAgwA0BbBNSjldCPESgJcATKiaq5YPhsAwYBX3/MkJbpmp7G9o8j47hhM7hmlzLGzIY1oORYZSOfYUz1nPX/+pFnrwZ/J8OXQ0r5bSA+IBDlfNZV6FLUPe1a7RfQsXODx4J5cvH3RG6RHxCwuDzI668FSWjUZfM/NObqokNWozWMeYtdw/xCWUn2fxYY58C93I2lZcVz28tXlfLruupySWDWHBNZPYhe/da4l+To2ARFWbpBwAuAohCmHSLOIATATQ3Xz8CwDbUQGBIYS4H8CPUspMIcQkAG0BvCWl/NPKRw2BYcA6ZrZYTvSoraxRJHZhZ6ZTpr7Vernf90QfzmGb+oFkpqN3827/48ELtDlr2TOTG/HXMKI9P2em6JLEZi3HC8zy+lzmnt8A8Pij7Py9txP3L89uwTkqHk+xYzj7uCKVAIz76Cuie7mxUCmUzET7fMI+ovKg8SyOepq1/WuieySPJbroLbZsfBPBWfIA4ATm+GNs7Lh3qZOuOTmkXD8sqKoS96SUsUKI9wCcB5ALYLOUcrMQIkBKGW8+J14IUdH2hK9KKb8TQnQG0BPAewDmAmhf9scMgWGgHJgWzXb84R054ml7ODtI50cu1eeI50ZCO85wJnFJOjuk7dz57Xv9BT1HwmsHqymBGTE8Zz77FxKeY5u7Q3MWShmtdD9JY2f2vcx/g001Qonm6urDJdIvB7JpDQCWX+bM7lXKHD//yUX+bHUuA0DMDL7u/dNfIPqJJ1kQfrWEHdrPv8RaJQAU+rIJzwEWWvGWgVZRej2wQ9F6wccai5Jyaxi1hRClH94CKeX/74iEEL4A+gJoAFNbw++EEA9X2Tr/KdnVC8BcKeVqIcTr5fmgITAMWIWD0nJ12Q/dic4P4RDNAdkcnQQArk58jiogHNLZUe4axVFRJU31nt6BE1ggvBfMiWUlilG5x7ssMFJbM6OOmqobWT5NZ4dNg9bMFKbM/4zomEL2WHfz0bvnTdnF9aZq/cYmKGel5XdF4L6ew3s9Yvn5fzOPBUROU979nxmoC7pixZYWoUfelonPwn/Qxm76+zkLZ9ZMiPJrGElSyrLUxh4AoqWUlwFACPE9gJsBJAoh6pq1i7oALpUxR1mIFULMN1/nHSGEM4ByFcKyKjCEEJEwqSsBUsrmQoiWAPpIKd+q4GIN1DDEbmEO1mAtJ3WpvbA9ftFzD2Q2R9j4ufBvXfjyHDlNWNuO76Q7vRM+Zyb3cG5/oovqs7koKJWvme/L15CX9WS141O5UGCTeSzIhvzIyYIOGSz4vPVGgYh4iB3+n9zOFW535/Lz/nIyRzyVB3V2saYjE5iuG8dF/5a88DnR96xlkxUA2OVVfXG9b+/iHumTYLv57ZqARFVGSZ0H0EEI4QaTSeo2APsAZAMYBmC6+V+9Wmb5MBDAnQDek1KmmYXPC1Y+A6B8Gsan5snmA4CU8i8hxFIAhsD4H4G9srlXu7VNbrOC6DNKnwoA2HCRzSy5W9m2rybhFXgz47XTrUVosTKG6E0XmLlnH2XHbd3fed15YTzpiSl697e/+n1EtGN/XldqCT+cM0r3qLgiDhe2hB057IwvlJXPfl7/M/dhP1XIfpIxZ7gFbp/VvNO3tFuWDpXjiPeN1oWQW0yGMnJcO6dmQFSZ01tK+YcQYgWAP2GKND4AYAEADwDLhRAjYRIq91dw/hwhxBkAPYUQPQHskFJuLs9nyyMw3KSUe4Sgh1Gjw6UN2Ia2A7loXN9aHEzxaSzb9VMXcNYwAPgsV+zdg9l0k9OAGatLCv/EpNrTFcDqNWxiCvydzS4Zd7EJJVeJmrohgntAHDyv29PveXQM0W5HWXsqrM+C7/RIZvYzO32jzfnBeO4n4pDDmpKace4IPeHNGsI3s8+nRQMOgfVyYkHXrzNfw9tBL0X+2yOccGir+HA/lqQPZupRZDUWVZiHIaWcDGCyMpwPk7ZRKQghngUwCsCVSJQlQogFUsrZZXwMQPkERpIQoiHMj0MIMQCAbnMwcN3i8ELOb/i1Bec3eJ1kRlz3gM4YhC9HLBW6M1O8pOyV3uzM9u6vO3LVWACQYZxFbZfEOQ9iCDNzj1j+uXs4MNOc156jlwDgje+Y8RY0YDOWw17eEYe8ytrVp3fpxQcvDVH8NYlMZykWqAabtCmsQ2FeqR8q7WuVApLpStRaVj09mDe7i5IvYpvPG8mzdM0pwpe1vEsdbZvzmkLNqSU1EkB7KWU2AAgh3gHwO4AqERhPwaQONRFCxAKIBlCVHnsD1zhqPcg2d7/JSmbycS5Dkd1Bb9GaNJUdu8UlSq3sbNYgXt3O4av2L+vMxvsUM7BCdxZKDR/mmk1+O1mLSbyb55xwnx4ZNOEdjvh66cdBRIeDzWBZq5jRunvpiXsNvlN8JUq13+ImuoZmKz7qsozoAzfwnD++zVph1BSusvtusG6huHXfozzwgW1rqvWczlEve6v+mRraH0PCliip/xoC3NywGOXs02VVYEgpzwLoIYRwB2AnpaxkVXwDNQ2pS/ilLmzFvy275lxPSY0+AoCA75V+2ulshkmvzz/F9Bt45+l5Vv+pesWwCerCML7uqdkcvuoczzWuCmfznJ4WWlTPeolt/Xe9fIDoY2u5O17aDbymEWM4BBkA5vzAyX7hUzjc1D6To78qkvddz4Ejvlr5cH7IH4+HEX3yHfYxDf2Dvy8AcO/mrY3Zgou99LSBXeNY6twXYjUV4JqFDVFS/zU+B/CHEOIHmARFXwALy/PB8kRJ+QAYClP9EocrvgwppV4+1MB1iVTeRMOLWy3A9TKbNwIb6yUjUrdz5E9CO6XrXKBSPXUN7/49d7O/AQD8f2AB0E4pbxHmwqaxtYls1rr0VRjPt02vn5TfkH0tJ55nxpryPNvgo8Yys9/optdPqt2Mn1fJRr7GoyE/Ev1ZJK+zPJh6sTfRJ7/nXBk167puDGtCFwfq1wz5hr8DWx2Zn4z5WBvr/B472wPwm3ZOjUENERhSyplCiO0AOsMkMEZIKQ+U/SkTymOS2gBgN0wFsCpfd9lAjYNjJmsUeUo9vdQbeQ9cz1HPmTjdnP0c7z3wBdE/pXOEUoLS+vStemu0OafFcze7nVOZOW9343V7RTMzr5Ok+FpydUfv2SE8h5uXUmzwV34YeQtYWF7arJf0rnOIhaN4kU1lE55mc1wE9mtzWMPFedyAylGp9KEEc6Hv19wH9pcULjgJABcusNBxX5GgnVMWpnXUuw8G+fF3UPkqWgasweyTPiKl/FMI0R1AFyFEtJTSqj2wPALDRUo5rrKLNFBz8fdo3hnaKTk+q7KZG50t0E0PKR25XPnY39nU43Cebf+FisbxwGI9TDz4Yd7xdn6FfRabFnIUVYkTr/vcvbyzL3LTy8T+1YNNJq1WPstrOMMsruQDLtAXfFm34Mbcw0l19k1ZUHUK4gz2imRn9XtpK9FLF3KinuNt7Ed59487ifbZo0elBR6pHHO/3EvvE97/uZ+I3t7CVTunpqAGmaRWArhRCNEIwGcA1gJYCuDuMj+F8gmMr4QQowCsgymsCwAgpbS19piBGoqR57sRfSyFI4FGNthF9NyN+k6y2I2V0029ZhKdUMyZxSP3cNvSfF92mgPA4bMc2TO+G4cTvfASmze+y2Rfy9evsC9h4FtsCgKA2alsxgr/nn0rF3qwoGvwNrd9zbuFI8wAoN5W1lIc34khOsVb9R/Y7gj+v/auOzqq6uvuk0oaKUACoYUSiqh0QaVKE0FBxC6CICgIKogde0OwUUQEqSIKShEQVECRoqiAdEIPEAIECARCenK/P2by+7LvDUmGmtG713LJeXnz5s6beefc0/ZZm8iFByH7NAr6kuz5RP/EvFxeJ8zHO/42NobhO3a5tKaT9UyN+sWvrUmuDnOyotvAfeZh5CilskSkG4BRSqkxInLJQlIZAEYCeBn/H6VTAMxSGIt/JfZr3c4hx3hHPKJ7F5KvaaIlOQAMrcjKPMyTH66N6RyCykpjA9HkTnOudYAnK+/HpvUnOXI1h8aO12UupJBsVqLjv9KYZwEsf3wkyfUnx5I85hCXxaeu5hxHto+pRA71Y8UZUYoNsKdGxeJnVuYWik17uVDhtlfYkK39nHsqEl7ie5V03Ayl1XqWK6lc9TA+7GhS0Q7964J6z4ofFNwpYJ8pIvfDkZvOnWxm7sjyQVEMxhAA1ZVS+XTdWPwXkFiTfyY1n+fYdc5+Vk4ZT5pDgx55rhfJladyUltyWIlWKMl/P7lHn8UAHPfl33h4OU7DpmleSYgWPsrUchyp0WY7ebfnh5L8xGvM3Ht4fhTJJf35PXwTObQGADWGs8GVc/yarH2xxmtchYcPX3P5D2wgmvbhZsy2oTx8qmk9LqUGgEGRzFKMRNdm97z/Rg/jWE4b8/64K9woJPUIgMcBvKOU2i8iVQAUiVi+KAZjGwBz5qbFfwYp9VnBbYrjhrmbqnMuYfMt3NgHAGWW8NOUorGHhC3luL3nMa03QUweo3EHVpHcb9f9JCes4pBVud/ZIBxuzT//htVNNtVjkRwqm9qPvamSYayYT1zLhi7ib3PbuXeoRlkSwUrTYwKXlvrPdZ2t9vtmPPRqe2PmjjqphQBfWduVZL9As3BB3cmhsopbjVMKxJBXvzaOPb/qX+JhAG5RJSUingBeUkr9r5dOKbUfDn6qQlEUg5ENYKOI/ArOYdiy2v8Iar7M8Wx1hstX/+nHBiKjCcfDAcBzJSvJpBs5jv/AS8zq+tVwzr8lR5qhnSe6c/jHN1VTvFqEyfsUK8GaH7GnFH+jmZQ9Xo8NlZc2iC5AGwMe+Tu/x9EbOAwGALXfZWOojmrJ5DNmSM9VdP+bmxAHXLOS5BLC98pvF68zM9Bcd1bVfAi9XMCUpvWNYx6vXjxvVrGBGxgMpVS2iJQRER+llMtfaFEMxnznfxb/UQxazsngrzT+hjjmuUP5qWaFTdwt/DTVepuNzuJY3gEHp3MXXcgN3CAHAJ77maEmqSWXklb8mDmvPEryrvr47VwmGrbZ5DUKiOdHxGsze1N7XmCjlVSdFa3yNLWI3sDmoVWVlV2lUbtvMvM3hWF6Q2af9Rb2hKad5OFH7bv9RXJMc/M7TG9exzjmCqr+bAYq+gVzJOSzp6sb57gDRLlVSCoWwBoRWQAHAy4AR39GYS8sSqf3tMLOsfh3Y8DPvUge134qyS8MWkLyQ5t4PCgAPFmdDcC8+jz7ekFt7gNI11KqzbQSWQDo3oUTuWuer0KyqssGIbEGG4xsX/Zaun7J8zQA4NvDDfiab7EXUvsmNiBbtjMFR61PzQqnmMeZqn1cx6l8gjYi4uNqJotuYXipF3sYPge4jDb2QR4eFbZDy73caLblybPcLQ4XOa5WfNfQOBZ6378o2u0+1CDxzv8EDgbcIuO8BkNEZiul7hGRLTCdLaWUquvyMi3cErVf5fzC639w8jPgYe6QDiphxr/nHWYD0b8yG4gdmRyaeOAPnjMx8+HRxjUb+vJrmjzJivXUaQ6DRU1mD3zEm5NJnnrcHDl69FfOg/hqm+yUcVEkV9YpTz4wFW/4TK2JcR5Ts7UbxeGjC0GlEVzy2rhkLMkjfuUigtId+O9b1plFkFHv6eNmzcR4QRjY02y+TMjUS4jdN0TlRh7GYgAvwcne4TymALxZ2AsL8jByO5R2gIdrCIARLi/Rwm0x45+FBf69+afPFHqNih9xeOi9Pg+QXGq7xp/Ui5Xq6wd4Sh0APFCOk8G3lueyz7fqc1b2/boc7hg4hstw/U6aT3yZU6zws715F6n3bmw4wx5GFX+zuHBqU24QDHuO8zdffsVNduUvgC4j5iPu/9idwZYuIoDvb8YH3DBXM4sJJQGg4necy9p/g2trauNvTh8cdPfj2pEtrl20OMF9DMYMAEMBbIWLxcDnNRi5w8bhKKml8hERqZXPSyz+pbj3/gEkH6/PyiXgDD8pIbtMig3R+JDC/+IOaI803v377eY+gKxeWnYZwHTRWF0V//ZbdtHCMkmc6A0LYGMQ1zI/Rlx+RCIWcQhq2j6mI/H9kmk+1lU2q7u63Mv5grgf+TXqmD5UyHX4nmZPJ7k8f47QKWxsM1qwB+h5zix33d/E/A5cwd4ss9x6V3/O+UT/ZZziHnCvHMZxpVTBu8DzoKCQVH8AAwBUFZHNef4UBGBN/q+y+Ddi791aIjeAlUnZ5fwz8jyVT1z6GI8IPdeKu67T+nAYa1BVDl88/oRJDKhPkbtnJNOHlDzESvNUTe7KXvMqh7nqzBpkvMd1Wr/C1q4cyhlW4weS177AXsyvo8wBDz8u4q25t8Ye4n0JxpN5J7EB9izD31FOSzYQ3pu5Muv0reae8ImvuOx4ek3XRsfOP2nmMLzjzeS628J9DMZrIvIFgOXgyte553+JAwWFpGYCWALgPQAv5Dl+1tKC/LfwZodvSX59PYeHgjSacTljVhtlXM8J6YCDbFSCXuCd+PtPdiT5k2AzL5J5hBn0/G5hxlWvz9gTSi2t5TxGcmX4sMe1ci8A0w+xB5GYyDnC4bt5nREBrP2zu5uPSqW32XCFfciEhS3DOHQz7xOT46ow7OvG92Zqd+7LiM3ka2YqVgXRPszxBAAVvbhcejpcMxhbPjSHYHn+m2IV7mMwHgFQC47u7ly3XOH/J/CdFwWFpJIAJAG4/3znWPw3oA8zer81dzs/37c7yWFro4xrZGvEAz7J7LWkhGtjSbXQv+92k7mgmka3fawzJ2ozgtjDyAjlJzq1AoewRo3lzwkAZ6L5nBovbib53HzuQIw9xWGXCn3zoQ4M4GT8qZYcgpoPLjEGXO+GLnGC7+fArfwYNynLCeukTDau90eZHp23BBjHXMGxO0yj7xPjvmSDOtwoJFVXKWXWqRcBRenDsPiPo72WPH7zC05Y11rI2r3/Qg7TAMDudFas33zIBIVeWhQrciXvZk9eZyqrHS9zDkOyWbkP7MDtQ97CsZ7XlnCX8XfPMW8UANz2JYe5do7jJkX/ZM7XlB/Bj1TseJPS5MN67LGdzGKv5btjjUhOaaGVsxYB52qykdFHHy3bzfTlNcvze7R80QzPhf2tTQqEa+SD5eaY4aeSGzkvcgmicVcP7mMw1orINUoplxt8rMGwKBR/zuB4d+ldrIxO3MgJ6qeWmhN8Q7ZwOMhDyy97aXnyw8+x6vD0MPsZPGK5nsc4pgAAIABJREFUJLP8r/zEfluXY+YHF0eRXE0zSh3TTAr1Civ5syY8xgutPIj5lFQSewtlppo0KcPK8pzwxGa881YZHJ6rAdcNRsh69sjKbOQwWEon3tnvSOby4fYDNxnX/CM+iuSyXY1TCoT/UdPD8JrC1XFZLY1T3APulfRuBqCniOyHI4chcLRKmDFDDdZgWBQKfaLe2YGsFDtW4I1KqyAubwWANl04PLQxgxXFQxu4t8N/ERuDEy3MsMzCO3lWxduNecpczCzeRac1ZWWf1IrdmtvLmuRIKw4wr1PlpzlPcvDBKJL9TrDWSAszm7lSw/kc/x0cnvO+6eJThMH7+H4N/pJ5nN4cxs2VJXqzUfrjG94kAEDg4YujY8160/xcsd9zGLEsXBvKVKzgPgbj1sJPyR+XzWCIyGQAnQEkKKWudR6rC2A8HN2FsQAeVEqdEZF2cJBf+cBBp/6sUuoX52saApgKwA+OhpOnlFLu89X8C/Dl+x+S/N0ZVibTvuW+gcyuZnnq4/O4WqjGOKb1qJzOO/Vdg9lg1I4ySzo7/8xJ66qz+WcR/iLH6SN787r29+Bu58VeemMaELmHq42yItmbqriEleChjloOYyx3owPAnte55zUojtd99BgPWDLHURWOQ+350R40jw3yqhEfkLwtg4NW395nNln8vIV7OaJnubamemFmXmRZC63RuFByimKMS6iVnKOxvwBwrfPKvQHsBDALjoa7WAD3KKVcowwGoLdJuILL6WFMBTAWwPQ8x74AMFQp9ZuI9IajIfAVACcA3K6UiheRa+EgHch9mj8D0A+OMbGL4bCOzEVhcVlx6+qBJJderA0NWrqH5L3tTMUbXEtTrN2Y8VbTV8jx5t3szo08ExwAHmrF1d0Lymh5vDNsdFR3puSovIDDXDsHmCwJbT9iRtxy3vyajydwovytvtNJbjDQ3DF3GssGN2Q7ey0BR/j+Xghqfs5lzDklOH8wsyPfq2hf9jCOpOod2ID4XNwA1cULmxjHarVmFgF3JQoRXPKQ1CgAPyqluouIDwB/OLqzlyulhovIC3BUrz5/Sd+1EFw2g6GUWikiUdrhmgByeQ+WwmEYXtEGkG8DUEJEfAGEASiplPoDAERkOoCusAbjiiLwd453N3yaZ0wfe5yVy+lhZrmlbwQrLB+NlC94HMfMj/VipZrOvW0AgKXxXJO56QYOu1RdyrtqqcRGaP5iVu6b8+HuHPr0EyT7nOZQT7ZGcTWxI087ygkwlX9AXf7se+7j+5cdxOuMXm6uqzBsH8o3TDI5LzI6kDuqq3lzGW7rat8b13wlkKlTXM2YSj4RrZSWx82D7giV/+e7EIhISQAtAPQCACerbIaIdAHQynnaNAAr8G8xGOfBVgB3APgewN1AvoXcdwH4RymVLiLlAeQtUo/D/3seBkSkHxzeCErA/3ynWbiIyB95l7wpnpV51LPcN3Aq2lSSSa05f9CmOr9mVRW+ZiZHZRCyy9y+pfzCnkz0dqaZEB9+jfJieUkKexxHs1gGgDZvsYfRK4Rbkb8+w5TdN/Vlb+v1fbdDx51l/ya5e0kOW8VrVVNvw3XatpBNnPROaswJZ73rusOvPNzosUbM9QUAi3ayV1IVRZrq+T/06G5avkUxPKI1cJY7j2gt8pmlRWRdHnmCUmpCHrkqgOMApjjD+OvhoGqKyGXgUEodEZELiVZeFK60wegNYLSIvApgARz5iv9BROoAeB9A7jYtP/rH834tzps+AQBKSpjNc1wiZJbjeFFCQ84F5CRpcf0eHGYAgNBX2M6vvokNRHAsb89O1+Adse9pc/tWej6XdYo3K8kd73MYK3gdezkfLX2Q5NEfmwSHd2rhuGGdeDLdEyGcKK875ymSq8/i5D4AfPcyG5kp5zi/E7KIS4iDwUy/RYFPMv/8a47idXz82X0kB7TgezNnZVvjmqve4LxHT5hkjQWhko/Jq5V4DX/PLlGnFjcUXeOcUEo1KuDvXgAaABiklPpTREaBm6evGq6owVBKxcBpDESkBoD/jbgRkQoA5gF4WCmVq3HiAOSt96sABy2vxRVE/M0ckqr+OSeTk+uzMRjysdkwOvodVkBxy5l6vPbTrHhzFCuSKY/xTh8Axifx+46ZwR3o99Rl0r5qTbiJ7vOPuC50aB/mzAKAytoY1w4zuCTWZzvfi+jKXKob20VzlQBkHOCqqJvrs7f1eyMOtQVz5KxIONOZu+19zrAqLjuEmx53b+RqpbTSJgdWz3uf0I6YCf2CMO1Rk0Ay4273GYRdGC5hDiMOQJxSKpfw6zs4DMYxESnn9C7KAcinK/Ty4ooaDBEJV0oliIgHgGFwVEzlVgT8AOBFpdT/MpnOG3NWRJoC+BOOoeVjruSaLYDKUzjMkjqDDci5dE4Eb0wzE9Qn3mNqkIwOrCj2vs39Cl4pnGCt2t2cCeF/gD2d8I2cX/hzI1f6/O7Dyv90c37Chz1nlv209OMY+4ADXLqbpRm2LfH8HulJZjnwvY04JOXvycmT9SfN3g1XEbiEDdXhtnw/D+/hpseBrZeS/MMRZrsFgMz6Wjd+O+OUArGnt1k9Fx6hNwO6MS6RwVBKHRWRQyJSUym1E0AbOFJG2wH0hKOitCccof0ristZVvs1HAma0iISB+A1AIEikrtNmQsgdyzYQADVAbwiIq84j7VXSiUA6I//L6tdApvwvuKI+YCbuip/oD34g7gvIyHDrLBJ7Ms73poP8c5cymsd0ZncuFfz6Xzq82vxrvhIa85BnKnJSlIyWeE90GI1yaNj2xhvsbAkK7Ttc7i3Q8+1VGsdS3KDGiZN+G9vcKb8VDTfz9ADF7/r9k5l7VX9azZKOT5s6Ja/x4OivEqZwaGDHTgvWAmuVWe2vibGOPZOJE9hcjXMVWxwCZPeTgwC8JWzQmofHPxPHgBmi0gfOIaRXPGB6JezSup8HFSj8jn3bQBvn+c66+CoRba4SvDfwknsM5oDUepd9jh+rtHMuMY/7zL5Xdx2Dt1EePJ76BVLLz76mHHNwy04tFNlJPM8lbqRQzsJ9TlOvySOd/IlppqlWPGL2cMody2TC+55gPMNBxK1azxqFgAc7cMGomxD7kmZUovHlg6Y5boSDV3NhkoFswFQu2JZ1ujnPdabyr3qUTbqrtJ4rPrFbCTuM0nP2+4zznEbXMKsqVJqI4D88hzmruYKwnZ6WxSKMu244eqtaszRpLOYVvAyeZ9SctgChHlwgvqvdJZLCIdyBo43w0Uj9zIfVcQy3jUffIVEaJEfhL3PO2bvw6YXc+BJTs5/P4D5psI8+D3P5mi9Cibpq4EN6ayIu6zjOR7lwYn2okAlsWFTJ7gPRmXzOnO2sIHwKmMy5OYcvbiQuW4YAaBkC/ZOE93UwQDcihrkgmENhkWhOHaG4y4DxnLyU5/n0KivmQzdfprJB+N28c6y3CoOF3lkFv70haziaqzVz2udyHEcToo4xwYiI5Q9jtgndIo+IHoYG5H+a7hqSlZz/4hnLZ6HkVrVLNVND2YPI+EG/qw5fhcf24h5h72rutfHkhzmy0b+wEu8+991n5n0Ht16Jsljq0cb5xSE76+ZaRxr/xpPawy7gOmCxQbWYFhYAD7aRJ87e60geXEcK+pD95nl4UcG8C5a75E4cqv2HnW5xv/vE2Yi/ZQvK+cf7+Gyz+9uZe9g5p7GJPsu4FBayM9m786OZznsVfNzVrSnenD3cpYvG75sfrkD7Xm336XCbpK3ntbpzV3H0juYYyPCgx/1BjMGkxyg1RSU+cPUfu/9+jDJQS6W+75/3BwmdUYbHW7O5HMTKFiDYWEBAOHvstab27QVycH7ObyxY4h5jYkdJ5I89jCHYnMUK9oBpVeSPNWDBxkBwOoz3Lj3ZLN7Sd7fk43MzD5MVjgnikPEfz9qEu6lB7N3dfoaTuiX+om9nITObMTyQ1IcezLzY9mQ1XovluQLofy+fRIz71aZxmHF6Ex+j/29o0jOb755yBamLXKVKCTU+5xxzMOt+cz/H4L8m8b+bbAGw6Jw/Mk0EgGvc4ns4RhW3BWWm8rmUV9mR62whEMe2drOvE8y74D9DppT/JIb8c/3UBsmDvCqwK+p7MWhnidKcVfxTYNMwr1ag5l5N+ZtHi3b8TmelrfmcebIumaMyYC75lM2EOkh/NlPteb7GzTTdQbXJh35fX8ryy6EZzAndIK1arBEH27GBICAI5ybMoNWBWNjkuklpkW6PhyquOISV0kVS1iDYVEodk/mksvr/FhJ+vzIgYQSv7KBAYAKwvQWYU9xSWZJH+5ETuzG1UUqRRuYAaDWeO7hXL2VY+rqKIeYOrzKrs/ZylreJNh84lNv5GtGT+d1zolkr8SzHr/nslmmERo37DOSf09hr2TiJq4yCzJD/4UixJtp/MIqcq9MUgx/Z6c1G6/yuRc1P+Dk++6CepXzQVIHcx5GVMt/0b7chqQsLICba3Hj3u5T7FEs+YIpNRou58QwAJT+jfejB+Zx8Lr8Alb+OUlaRU62qcASbmTlXdOLk+2ekZxoL/stK83eZbh7XK/2AoA7dz5HcmJX1gqfXfcNyWOeZ+6ohGZmPufpkf1J9tJ6JnwqXrwS7RbKBJG+WuxnhQcbwtQMrlLLzGd06tqjUSSXwk7jnIIQ/JOZ0GlekjmrVvzgxiNbrcGwsACuCeRyyK3HOYF9w0wtaRFuhhlCt/OO12svx9Qz6mjhioq8A/Y6a+5OUytzb8Fjw78j+fXveQuc+ionk99Zy9VKUtakZf9kCXsDjyzgfpBhO5heJKAWV0Vl3G5OCkw+wHmQHH/OBrzUjEfcfveW6xxzo+OZiqVBMDdK/l6P57InZPP3E9jIVA3ewvfrjnzbBM6P7ccjjGMJWnjTC+uMc9wC7jVx74JhDYZFoZjyAyeoo5vGkrwjiJO44eE83wEAEl5iZZN0iqkpwkpxvqF2KZ7PsGY9d1gDgG9ZVnJvb7mN5M5tmIIjvCPX/06fy58rcrXJb/7aXuY/qjJPO0c4rr+/L/+9+rsmHcZu5jxE1Sr8WdMU7/YvBOf6cw5i1WH2MFZ78/3PrsjK3CPdNPrKW/8srvWHpOwyy5bf/JwbOt+sWt84x21gDYaFBfBsF27Ua+y3n+QlpbiGf0iYqUgePcg73jUnOW7v78OKNjmTwxfBMabixS6uYCr7G+/mn1/E4Q5f4bDY0pu4V+FEXbPhcGy1BSTP+pBzEjX9OSEdpTGyzvzYrO7yW8jG74DiRPnOEH2Gt+ldFYYDd7CHluXPskSzgQ4K4PDeyTizf8QngdVF5Q2uranGyN3GsYFx3NMT7sZ9GDbpbWEBYNR0Drs81oNDJqtvZ8U7t43JXhC6i5PWQY24aS5lLYeLTmscTaX2m/WXJY5wzmHvfazkWk7n0tLIVbxr9t/E9BllvzILRZ+N6U5yUgon4w+VYSqQ2PnajOo/zLxIRheWJYNzFqtmNuRrXIASrbCc39czmQ3CuWp8r1q8wX0vD19rzqWI9uYQYIdXXJvTcaxrDfNgW23Ot8kw7zawISkLCwAVP+Kt5A8jOaa+dzjvkL2itNZvAJ53scIK9GBKiLTvOSQSuZKvcaKuSRNe9VWutBqmJbF3Z3Cu5a3AO3kNj/Oue1iFxcZ7/H6WPaGtQ9ibOh3CoZ0MLihDwg2m1xL5Oxs/zzTemmaXcLVgNR+8y2Wyi2otJLnJP9yzsuTT5iR/VV8bJQhgRecPjWOuILOTmc8Zey1TvlzIsKhiAdu4Z2HhwIm5nJDuU5V3vJ9/xjvksl+aO3U5xKEalcZhlrMDzIRoXpypZh77/RAnTPe/yr0G17zDZIQ/dmWFd8dU9kDeXvmA8R6hu1m5x9/Kj0zFxlxiHJjOHkh+oZ3MVrz7D57FxvCIRrsezbq+SDiymL+zlmN5GiH6cBPeXU/yNLzqvnpYDHgt/jbtyBnjnIIwuJY5ca++jzlgym1hDYaFBdCvKtOAv/dXR5KDbuGdY3Kc6Q0cfpyPRazhXbTeBLy3O4c/aow352btepzDWH47+Zx9bVh5D6rJSvPhSazApiwyQ2krhjC5cu99HJ7L6M0eRFJrJu2rNdvM5xzvzlQqOd7adLzJ7F1dSGg8R3uyfQZwpVv6TB4+NaMTNxP6/WhS1Jf64uLGp76xyhygtKWhnghxtX+8eEBgQ1IWFgCAbw5z+WTNfqwE09pqYZpqZoLaM5hzGGtGTCa59hROftb4kCk3do8y+ZV2txhPctVSj5LsfZSVYmYF9mq6an0XYXWZyhwA7nlImzKnKQWP0qzgzrTha5b+x+xuLrWJDYLnAVbmWRqz7IUgoi17Pr0rriH5lcacm/m7yQSSM28wtd/uF7VZ41VdCx/p1WAAsHgfU8xXgtn06TawBsPCAsj5QCu5LM/egf82VgT+68yu7KzR3IjXojs3r2U15320/xx++gISzUqh5gOZBrxGHJfZHrmZy1NDdnOifWxQS5KbVeDqLwA48gbvtE9/xDmLwN+5vyEnXmNw9TApTRIasbdV2o8fQ48/uCxZZbnuY/So8CfJcZmcrylXhUOEMRncMOcj5k6/77oeJFfGZuOcgnBicQXzoFlp655QgOT8+y2GNRgWhcJ/DydQj9zKSW7vFH5QfB8wuY8OH+NscE4aK6TqM1jedoIVb7mWvGMGgPveYYLCWwK4bHPiSR6uMGcps6WWncE75rWVTPJBL832JXZjw7XnMx5tuj79R5LTupo9FR8c4jkeW+qyFxJ6LYeHSn/mepXU7O6tSM6J4cFEwSXZi3kzkfsfPOqaY2Iz+l5cF3bQQdPwlYy5OELD4gQbkrKwAIAsfoxLaMRD1z/NMyGSMk3FkrK6YMpupREDltnM7/n6I2bmV9/xvneK521HT+PejujNHO7IaML9ENXvY6JBAHiqLBuElzozxXf1nL4kl1rNXkxAglkOHLCNva3a6WwMT7Uww1iu4t1F00kesvsekgdXWUby6Ee5asrjjOnRVZ13cdSyJ+qa1V++SexteZlcje4DazAsLIDs0hyWOd6ZK1uW7WbFq06YnEEV41nZBKzlHW+O1mkc34xjFYPf5xAWAFSbzxxXCGDiP90TSunIuZZPHphEcpCYFTvb0vkaTWey0fkyjKksGnsMIrnFdeYwKX22x+za35L8bFwnkhO+Ni5RKLqt4NxLpTmsrJ+9oSfJ2bextoviVhsAgOfy9eZBF9C47XbjWEwDLtEOW2qc4jawHoaFBYAnvplL8rCt3HkWOolDO8camTvJE3X4WEYQh5z8jnNTXWhM4XH7XR9zUjskmBPOo2ozD1SkF+cTHnppKMlth3I1GADM+4b7Eyp/x+G2aU+2IDkolj/nwiM8YAkAcipynOuBD58mucRBnVrF7F8oDD7+Wuf8Y+wx+Kxg6hDPNC6NPvKU6WFkdOTejCovuBYqW7umtnHsrjY8hOkf4ww3gjUYFhbAsHG8Gw1sz0nu5Me0yPN6c5aCPgs7QCNk/fYsx8xTcji0MzXGpNjwX8PhjHOBHArrndSL5JGNeSefUoYX8cMENg4AUPmeWD4wj72YWloyP7Mce0aef5lhrns3cXI9sR6X5iZksEf3zwXQK4Uu4Gue7cbfUeX5bPhy4ln2WGQSMeIcf++uBqie6/S9cezrOM7XeOOAcY5bwJIPWlg4MLjfHJIb+8WSXNub+x0mVeOdPwCEeLBy1qtyGpTga55VfM1yE1sZ1zyq2ZDve7NRumdTH5IndmA+q2wO6+c7ZW5EFf7stX/lda1JY49i/mlO7v8SZ869nt1EC9llawa3IofBgF3GNQpDw6e5v+GPSUw3cqIZGyX/49w/4nPKJGL0PKvlpo6fMM4pCGMmdTWONezOSQuz28Y9ILj0XFIi4glgHYDDSqnOIhIGYBaAKACxAO5RSp06/xUuPazBsCgUE2N5oM+bh7kBq/YI/s3GjzR/VqNXs7Lw4gpY+CWyss5mBwNhx82uYq9zvJvv9il3bqdcxzmJ5AHM++RZiUNUr/dlmgoAeG7/XSTv2MpT/YJ3cc9J5AKmbQ+taTYxIpuT3DFjmYur/AK+pr/ppBSK3Y04pFRa46MSL67e2vUJG7oVd4w1rpmujdEdWPlm45yCcK6iqVHHV/yFZFcp04sV1CV3MZ4CsANArnV/AcBypdRwEXnBKT9/qd+0IFiDYVEogl7lMIx3BCuO/fdw4jInzew9GN6Xq3aGLuQKp8gfOdyh/DSLkc/DqO/oyv/CsX+P2byOnDDOtejzuT8v28p4j4Rz/Jpq33KuJZ3bG5ARZYbjdOybzASF0AYoVXiGk/mJnEIqEvSy2BJjuDR6SIWfSX54OZcUdxrLg6MAoNJXsdoR1/yBJ9r+bBx7/+T1+ZzpnriUISkRqQCgE4B3AOQOnOkCoJXz39MArIA1GBbFDbseYYNRbRYrzaBDvCNOS2YlCwBDTjFP08QuE0mudBcr+2re/J7vn+RdOAD8dSqK5JTOvGt+qDwnVDMV/9wnH+Qk7pFkkw5jYp0vSQ6YzpH7jekcfltwkhVvdIA2ORDAHSU5tRsgBWcDBsK1nTwAjF7A97fDjzwjvfdapkkZ0ZVLsUZEtjeueeoQV3cFfe2awVjW3fQe5Iy+uXDToJRr5IOlRSRved0EpdQE7ZxPADwHIK+LGqGUOgIASqkjIuL6ZK2LhDUYFoUidDMbhCNNWfbUCmpK3cZhGQAo350rfUae40yu1OFYv8cJPj/lOjMv8vY4fsZq+XD1UZPfuLT0vcbzSP6sJg/Lvv0HrlYCgKM1tQR0ShTJU7ZxIuXDRjz1z1/MaqNHt3Ivhz4sas+nXE1UEmz4ioKn27AHVyOM741s4ibHKe8wzcfpV0xPyac7U5oEuVjum1rZbOs+eBu7aNWfdFODAZdyGCeUUueNvYlIZwAJSqn1ItLqEiztksEaDItCERjPSdmkZpwQ3daKd7N7ssxpbdv/5sa9HWmc2H0whK8R4cE/zTGnzNDFlyd45/3POFZ63lU4dDYslr0cPY9SeaO50x8Uwq+J+pyv6fEEv+adXczompFPtVHl+5gn69Yw7u1YNpA1T/wM4xKFIq0KK+Kjj3M+x2s15yyaP8g9FqnDzHsR1zaffIwLONDFLLeuOZ5zU+48g+gSJr1vBnCHiNwGoASAkiIyA8AxESnn9C7KATDd18sMazAsCkXgYB40VO0FVoINtz5JcoWl5jwMtV5jbb2RDcDvGVzF47GbOZqSbjVr+H1PsyHz1DbFtVqzYt58gL2U0os4T+K/gd8TAKJ/45BJ9jm2MtX2s+HLjtcJ9mKNayZvv47k6RlsTHO8WbF6wMXRdgBK7OQy2cp9uUdFAtmD2/s5e40+qWaTXo0tPF/E1bLayovMmM3QuTxbfES1a128ajGBwiVLeiulXgTwIgA4PYyhSqmHRGQkgJ4Ahjv/b9YpX2ZYg2FRKPqWZ86m9eN5DsWa5zks43nWJB/US0Vz1seQ7KH9fdcrnLSNWmSGdk4/y8p83DUcYvL3YE/nI992JD/wLtN1Lz1jKqtFsznPUeJmTh5HDCx4nkPOjeY1fQ4yj1Pyddzlrof8qvCk2SIhowqHt30OaQy43vzon7mdDXjiNaY3cGvnv0je3tA4pUAMGz3ZOPbyG0ytEnwB4bfigivQhzEcwGwR6QPgIIC7L/s7arAGw6JQ6GGWp6vzHIkhk3iuZrAH9yoAwIA4Dh/tfIMVVOkXuZnNaznvZk/VNK+5ph7nMK6fyp5OmX84RpBahpXgnnieS3HkflP5r9QaDns8xNQfMYO5VFd587rzi7EsvJ2D/3W0BP/RbDaEPV90Pel9ujrfr9J/c/TCI5hzM8Hb2EglVTVzGCcy9GIG1zrQ+67qZRyrEcOf1a173y7D4pVSK+CohoJS6iQAc2jLFYQ1GBaFokU5Du3U8NGStJlcnTTisFlhk5jGTV/dRv5E8qks7kzeU58VVqmfTINxWy+u9Mlpx09sWik2ED+//AHJPerwICjfmmxAAGBbY07U9prAUYCSnmxkPt7PzYFnZpvJ+jkteGs+tAdXVu19ke9nFTC5Y1GQrPEX+szlsNe5H9mr8daiiJW/MsNzW7Vu/HAXZ43XftNs9Iu/jT3L8L+MU9wCdoCShYUTMWdYudy3jedQ+OxjZV7lI3PKXOYt3HvwRRTPRqiwmHfAEXvZSHlGccMcACCcd8nZQbxLjljB4aPucVwFdeg1Pj8o1nyLsl6sSd/rxzQp+7jgCdEV2Zgm1jcJu39vwMZx1yStgzpV4025AJS5kYcyZeew8UzRJuJ6sqOEJnPNYVLHjmrzyc3evgKRFW5WSUX+yLmWi+PDvYpQ6nI07hU7WINhUSh2reOhQRHr+cEokci5guwknTwP8NXOOduZYzWJxzmRHhTOFTnH6puU6bf04iFBh79loj91lJWev1bzHxjNRkxXogBw5x/sxVRfx011/o01r0Sza0F7zOmDCf2YPylYyy9H/M5G6kLU0E1lOMQ390fOxZRfzarZZzGz7u4fq3UkAvDuYB5zBSevCzCOlTzIhQfee/YZ57gLLjU1SHGENRgWhUKV47DLc28yid8bO3gORfjz5mS1vQd4hxv8D4ddTt2hxbIDOXFeK9gMkawaewPJ527UurDn8I729BxeV7feK0heeMhMUJ89x96T70KWU7U8yb4NbDFqfKVRsAPIqMXrOFWLr7nvLjaWVViXFwnlfZmuxU8rwMwI1GaqN+OS5JQAUzVkBF2c5xOxyDQG2cfZC3TnPboNSVlYAKjwDSv3ceO5OCMilXer2YFmA2qtY1xvn71DI9QbwwrssFZ5lTTBjH+XCeRy3zr9+D2Ss5jkz3c1Vwr9uYSbBVM+1OhIAHzakCuvorw17ymKxZ5DniF50jpuFgSAt4/dQvKSNdzEWKKCSa3iKiaP55kaohU9haxg5Z3UnL2tHHNQIMKn8WwPVzfU+x6vahz76RGEu7FBAAAgAElEQVRudOxbqZlxjltAAbAjWi0sgGM9ebfvuZ53wB3u5tDQnrNms9oHUcz62mHhEJIraoNzvLWmuvTZZtNY1hhO5Ca/xwnmDz75lF+wgMW1qdVI/uxrVrIA8MJizteEL9AMXZDGT/UBh5PueZI/JwD4H+YPVzGCVa/3WdNwuYrg/VrI6TR7X7Hj2KhnxbBFqdXM9AYyt2qZ9C0xxjkFoU0ns7ej9XyeSVIda41z3Ab/fnthDYZF4ciM49izaPqsawg3ltUPN8tTl6SwMv/p9o9IPtiRw0d9VzxCctCn5ojXgKXcIX2qO4dVXm/H/OXZYfw5lJfGgdXTTLlO6z2O5Ieqculuta85pFJmEr/HAdMGISCWjV+n+wuuNrqQeRgPjljE6/Ji7+vz7sw4vLMvJ+fvjDBHGd2shZRcZavdNszs1vfslU/PjpvChqQsLADUmMT19pUmcT6hz2wen1p5iWkwfGKY0junLM9fyAjXBhPFcegnrhOfDwBHX+Vy1NF388jVF8bwPIzyi7gi52A37pko84f5xL8yiQ1X9QOck4i/uzrJkQv43gQ/ZhIx1mvAfEmrhnP4LbsnG6FgMO9TUTB2AtPJ63xfwRXYOEasZQ9jZDlucgSASV+xMfTHn8Y5BeFwKzPOVSYkMZ8z3RS2SsrCAig1gRXtoAhu3GvalUtg69xrkg8ezWYP4sUvuDw1pTZrtAfr8hCItt5MbQEAPYN5+I5ONvj5IG7sSxzAyjtDsYfR3M+c9nb/dl5nwDvsKUUu5Wzywfs4bFPpMTO0c6w09x6kN2Vl7TNDN46uG4yz13IIKnIJf9bEvpwnSd8SQvKN5c0igy292cvzd5F2vW4z83Ns2MBhwerYa5zjFlC2SsrCAgBwcwjvqm9fyP0M4zpOJflkjlk+OT6uJckplTgE4h3HCeoZGbzr7trADJE8fDsTA1bbx+d8FMHJ5ZT6XB584nre8fo8ygluAGgewQqs4aRYksM8WfHqfRsf3WXu1OPO8WOXuJONac3xfI0L0UMr2n1McscQrckxiw1IdjUODWUpkxpkTf2vSHZ12NGm1eb0wVe6cm7ra5ihR3eAo3HPehgWFhi5lMtmo6/l8NKgP1lx+68zeyYytZx1mMZiXfIg74j9dnDj2cbrmF0VAM624Z+vb2MOMZXczYo3sTYbCI+bONS2L8NM1v/24Y0kr8xhOdOPS02DY9lTimuljWMFELpTUyztmP03xy+fEiUX0fInnn/hkcIGovwKXoOeiI8P5lAbADSrxJ3eYS52envlk654dx5PNKzi4jWLFayHYWEBXFs3luTO4ZtJHr6FE6gN7uFQEQAMiOBRnLW92cPw9+BMevWFj5EctNP8qVb4jsMmZxtyuEgy+T3KT+MwV0wt3vH+FMAKEQBaPsNkeFs6czgp/s4okjMDeZ3hG8xO79CnOfR1Rxg32XWfxUUEriaXASBkExudDwZr4bmOHJ7r6M8d6s1HsMEBAN9bte7vScYpBaJjFzPnMX9143zOdE9YD8PCAoCPJyu94T/fcZ4zHTj2qBlWGBbICWjPXazsJYy5KWpnsAsS847Z29G3/wqS55zgEMmBETVJTm7OncqP3cy1vN2CuM8AAB7VymJ9qmu5gTmco0i7lpvySmyMNa6Z9Q+H7Fb7Mt356nR9uqCZWykMZSey0floDnNcwZc9n/c7sLENjTVnmiQEsAfmajJ+7nrTSyy39uJpUIoFXJu457awBsOiUMSP4cSkz91mAjov8pus5pXMVTm73+PdfdXZbJQkm5++ksFaYwaAt0cwkVPEfO6R8GjK13j3Gd4SN/PlUtMZZ9nAAMDkMVz+e+ts7hvo0Io7qodFcEy+3ZhnjWvq4bnKP/D93DuI8wfVHnTdYJy6l2txU7tx1dkvjb4g2QOsuA9kmTmM633YC+z4nmv1vjX7mwb5zD2u5UGKLxTENu5ZWJi5gOA+nMOIeYd3xId7mGW10S+xwqq4mImbfGM4Z3H0jiiSlzUwme7O1eOgcefg50jeOIQb96rP58TvPTdyiGTWOjM8MnYfe1NVh3Nj2X5tjscDtbkgoOJhLpEFgDO12JtKiWRqkBJbTP4pV/HMMKZQ/yKOO6h73MrlwpOX8KyKyvlohk4dH9SObHdpTaceMo3DiQasZIPMugP3gQ1JWVgA6WW4R6JEKit7vyOs4ILXmEnvs9ooztRZ/JrMFlyOeqopJ48bLXnKuKZHqjYlTtu5P3OEDUBYJU5yH9fmO7zWXGsFB7CufhTJVR/ia0yYy9VcVb7jv8cMZkZdAPD048/msZ/vV8jui1c8M48yEePu7Vo5cC32KO5+ikNvXunmGkp4FuxZFobEfIbpVfvWHIzllrBltRYWDpyswwlU70GsbPyFeZ6qdTZHDcd8wQll71RWSJ6aggpfzuGPUkvMeHlWNVaCrSby7n/ex1xWW2Y9ezmHt7GymlXHnE1zuB33JyTdxbvqqrO58Ux5sxEru8z0FrJ9+bE7eT1/9tTSfH95BUVDz3JrSB6yn3MrPd9aQvKnn3OjX8nb2OMDAB9fLnNKa2mcUiBqjD5kHMsJZYPq1jrXehgWFkDYTs4/HAvi4UY9unMj31dfmYo3hyM3kHoc5lrSeDzJncZyeKnREnNO+PZTvJufvIQNRNTegnevaR05Bu+RYT7w2lwnHH2BCfTiunM4qdQWVnmhf3LTIwAktNbmgAdxgvlcxYsPSU2J5xCUxyk2+nMfbE1y9Gg2yO1LmTNN3vudeU6iYX62grB7QCXjmN6BXsn1WVHFB/9+e2ENhkXh8DnFfQJRH/Iu+5uzbCBSqpulpHo46Nx6NjqtD3MyWUXxNQ6nmPvs+FOcXA+I5515x3ErSP72vQ4kn7iNcy1Vx5pPfOWRrDgPDmY6Ek+ttyC+I6+75E5zUmDpGVzBFP4jV2+pVL6oeTcLx9ZNUST7VWKD++0Czlkk5rDRCvLIxzO6ib2SeTDpWgqCVz4RraBD/x4ta8tqLSwAeKSxMhGNobXc6L9JjvQxG8+Sb+XSUT9NDWb7srJPLcUKK2W6NisbQEZfDltl1GZPaPb7bCCOd+DtbElt5sahNuaAoIYa5fmNAStInrqeBxM1jOaKpkrTTK6kfcmsaOPO8v1MTOTPWr0HV2IVBVHXcFnyma85JHXvY/mwIuZd00M1zINaBWxZF5vsPPIZp5dc/uLDb8UCCkC2NRgWFsgqyTX7cV04LFPlGy2xe0hr4wYQtIOVXnYM041ktuc51xkluaxzznfceAYAHyVySGnOlFYkl1oeS3LwdF6Xhy/v/g99be6qW4Rwqe68m7n0ttp1HILa8CgPUNqfDx3G889wKVCJSDbIPyQy6y639RUNXm25z+XkJC5UUB5sEK7rw82W+3jkuvM1F7CQPHivzzTj2NiHul/cRYsJBMp6GBcDEZkMoDOABKXUtc5j9QCMB1ACjvG9A5RSf4mIN4AvADRwrmm6Uuo952saApgKwA/AYgBPKfUf+GaKEY41YsVaqzXzK8WA+zSyAjjcBAD+h7WdZA0uLQ1cy2oxy78Kyd27PWpc80AnLosK0Fo1ZvzJPRHLUrihcNxgpj+veD93sAPA/Cj2IBLu5gbC0zX5p1jzfS0J7mmWGE9pwOVCMaNqk/x3h1EkPwBeQ1GQ3YYNcPQkNkon6rIXmPAQez3pT5uBsLmdR5P8/CtNjHMKwhsxnY1j/uX4t2XW17kRLpFaEpGKAKYDKAtHHcAEpdQoEQkDMAuOsV2xAO5RSrnufl4ELqeHMRWOMfHT8xwbAeANpdQSEbnNKbcCcDcAX6XUdSLiD2C7iHytlIoF8BmAfgDWwmEwbgXAwVSLy4oBfbjcNCmbH+sje9njCHjYZKtN2skhEb8jmnb342sevZcVbVaqyckEcG4lnbkFMess7+5HrLyNZOnAW2bV3az7HHgDU5p0COScxnP7mQvpQEIUyd75DM8rHcpK0qMEx2qeO9xee4XrE/i8fmPjd+AVLjH24FuHrHD2EqvOM+NHd6dwj0lV/GGcUxBaRprjageMWsnyPNdpUIoNLt0+NgvAM0qpDSISBGC9iCwF0AvAcqXUcBF5AcALAJ6/VG9aFFw2g6GUWikiUfphALm/zGAA8XmOB4iIFxybjAwAZ0SkHICSSqk/AEBEpgPoCmswrihmDe5Iss+PnLMIFm6A81pu5hu8DnO8O7sV00ScasWVQzlaVKtX29+May76hOs6y/zE+YP55bgSCAO1ok1/TSkmmbmX5bfz0J+lh9nD2P0+V/6EaPu9iPmmkkxuyt6T3yYOhe37kj0OH/D9Lgr2TGPjVzqUeaASDrGHd7Qp99ro3egAkFPO9JZcQUwnk9yx3Svc/xHt4oyNYgOFS1YTrJQ6AuCI899nRWQHgPIAusCxwQaAaQBW4N9iMM6DpwH8JCIfAPAA/udrfwfHzTgCwB/AYKVUoog0ApC3rTgOjhuXL0SkHxzeCErA/3ynWbiI+OasSEuFMWNrcAw35R1palKDVL6PQxzNS60gefJOvmbNMNa8lXzMjuk6/Xi3f+/LrGze3csehd9Z3gGmxXGy2be8Wcaz4xUO1Xj6sqKt+hkn0hNraUEVDzPw75vI2/tsX37N3C849HNfBb43RUHICr7m488sJvlsFHs5d9/GlW+BYqqGTE0j3gfX1pVwmznTu/IP+WTC3RQu5DBKi8i6PPIEpZSZpAPg3HTXB/AngAinMYFS6oiImARrlxlX2mD0h8MYzBGRe+Dgu2wL4AY4qgcjAYQCWCUiy2DUZQAooNrZedMnAEBJCbN5jkuEzFBW9oldubooeDifX/KQqQSOfca76qU7WRGXC+KQ09Frokge5c0yAGS24ka8+Kf5PfxzWMFdM5ZzL/tDuSoq+x1ODAPAvu68+/+r3RiST9zM79FhBY9wHTzU9A72pvP7/NmZFWm7Jr1ILoWdxjUKg94YOe1ppjjx38Jhw897cNVU53vNCqilEziXUtrFKqnQXSa/ecUPuP8j7geXLlmMoICcIrsYJ5RShZJoiUgggDkAnlZKnRG5+kSNV9pg9ASQy/HwLRyJbgB4AMCPSqlMAAkisgZAIwCrAOQNflfA/4exLK4S6pdnZdPtm/UkvzT//kKvkfEgG5WWkTEkz/udnyevs2YFU8MIbhzb/yobgLNrOATybKm/SK4XyZ5R71e7Ge8ROSaK5Pun9CP52A0cuwk/w4p65ltmO/SeRzhkVy2YE+UhH5sDqFzFa29wn8XAtTyzpOw8DqWtH8RezdvHTWLBXoNYmy/6LNQ4pyAkVzR7Ulbu1UgoYQ7KcgsoXNJOb2ch0BwAXymlcmcbHhORck7vohwAk1LhMuNKG4x4AC3hiL3dgv+fPXkQwC0iMgOOkFRTAJ84b8xZEWkKh0v2MIAxxlUtLi88+UE4/moUyVN3c/jD+xFzJ/TyA7NITslhj+KL97uQ7F+WrxF2i0lVMaAcJ6T7H2OlOLo3e/n9fuBKqzLVOcyVEGv2YUTH8a74RF0tuK9FnJIjed1pnc18zvUteVedMpaNoW8CG5ALCdoMnsF08qX283d4pAuH0mr/zMSMNceY3kDi9ZwYD3Ex6X3/sMXGscl7XA+3FVtcohyGOFyJSQB2KKXy0iUvgGPTPdz5/+8vzTsWHZezrPZrOBI0pUUkDsBrAPoCGOVMbqfBmW8A8CmAKQC2whGGmqKUyi3z6I//L6tdApvwvuLwTNIU2l4etrPjZU5Y15ho0njMGMn9Czn1uA/AszIrtIqjmSPCYwrnGwDgkY968Tme/MRW8uKQ1d3NmWtqwVwOsVRbaSZ141vwbj+9Aec5qo1gdZ4dwI1+ex4yE+mJb3A5l08S5w/iBnI3eeQI1zeSz9/LJcWlvLjSSu/1WH+cq9iSapml0Trfl6tYcmOUcSynj2mk3RWXsA/jZgA9AGwRkVxO+JfgMBSzRaQPHJvsuy/VGxYVl7NK6nxxiYb6AaVUMs7z4ZVS6wDkw3NpcaVQdi0r4sNdue7AS2tmTiuXD1ttY1aCwfu4LyC1DG/V5+7k4UYZ+TyMq1JZ8b6+nuP0008xk+yxdN4hL+w7kuQKj5ulu20HPUFyiT/ZeO58ShsTG8vXqDXW7PTe9QiHct77lA3GzKPsxaSMMC5RKN766U6S/Y7w/S0/nPMP4VHscYxe+YFxzdhMLmYY8bVrj+WozaaHMSmRjfo/HxmnuA8ukcFQSq1G/vlbADCJ2q4gbKe3RaE4VYOVZAmtdNRP2wB/OpYbzwAgJYd32vcs609ytW9YYfVuz/Macnab/c4qm5Pxvi+woerVjD2KNK1VeeTRdiSX8TE9o6M38WuaNeMy2QdDOCH91WROHt8xe7Vxzeb+HJL6KbkOyUnpHOu/kAnfP9/5Iclnc/hRv/5J9oSScji81HABl7sCQOl1fC9CXQxJvXuko3Fs12nOMwVir3GOW0ApINutuXaLBGswLApFWjjvnKY8xmmke3/j+PdTvQca1/BMZeVe+xTvLNU+prLQ92rpHcwk7Nn+fI3eVX4k+eFnnyG59IBYkrdolN+eCaZqXnwv77TDPHnj93ALzpskdmBlP+2d241rfj+Xw22pbXin7r+fP9eFkA/e8bc2E/37fBor8qDUUjbItc6abLUqk8NvrqrHliFmtdfeM2boy23xHyCgsAbDolAEVGUF9vLDfUmudZrj4wc7m3HpVndxJdXJdM5JbD/OpaVnjnPuIGCv+VOt9Cw/oEv3cVl6YDp7GBkLuTendmmustr5NPNAAcCtKweRHD2KQ2m7XtFCUB9yIj2+rakQwxYzxd47lT4nebnmcay5nr2BosDPl9eZcLNmdrLZ8D3+Elcn9SpZeN6kQ2TdQs/Ji1Hj7jKOlfnHTK67LazBsLAAKvTX+JFS+CHfP4QVHOqYoZ3Ff3MOo2X9HSS3rcS7z9LV2AhN9OH5DgAg07gsVkqyEZJU/nnvfpUZc5s25zVUSjYrtjtHbiF5eSJ3ft92HZMTLn6G38OXHScAgBrIu/1301iRZpXWvQHXh0SEv8ByxFgeq7t9N+eh3lzD+Z/3g80CgEBtUmAYdhnnFIQ3nzTJB996l+eyh65y6ZLFBwqAneltYQGkTudd9Le155L82lFWLuH55AJWBHO9/W87tCqpkxwOyg7ggIfnWbNjusYinvTXsiT3ctwRwOv49DRPfJsfzzvktCnalCcAn2neUrOZHGP/8xgn3kuFsxHL2GLOjEiqw0nv09X5s6VGcyd4NDtKRULH2fyiEsIex95VUSTrHIkhzc1kffwhvheu1jdtTzPv77MvaLPHJ0cZ57gHFKBsDsPCAmNqfEOyLzgJnprNyr6aL5fdAkC1SnysaU3eerfTkqwvtlpI8vYUU9ksXMXNfeFt2EBMOxJFcuxsZtUN38AEiKEHmIsKAEJ/ZSVwLIT7KrLasUEou5rDd1klzZ16eijfr0ytYji/XIqr0A3EytNsoKW2ZtT/Ya8meKipGjyjL47ffMWDjY1jE/pw5Vp1XIB1LC6wISkLC+B0NidyN6axkszRqo8+3t3WuMaoa9joBGhFg5u7cqfxlDOs4BYtN5VN+D/8gP424waSgz7hHMW5ZtxD4bWU5ZwIs3P58C1cShq5ghVtxDiN+kMbHrV3HM8yB4DgdWxwPaP5mg/X5I70Fc+5Tvo9vwW/b05FpiMJrcVWKlVzhHY+b3ab+wVwmLD8PNfWlFna/Bx33LyO5O3GGW4CBVslZWEBAO804lnZmdcxZ9PhFmxQMmubicwnRnM/g2jty5FLOcmqDrOyr1bLpPjeNZBDZSUSeR0zqnCP5/slOf/waxWm0k4ub+7sQ/ZysrjOBFZpOx6uTnLsXZzkvq7KPuOaWW9z8v1MPJeWrvxBN1yus8Re9zOH6+bvZg/NawNb7Io/s2eUVMc0GN+35Oq4AXCNijxpiBmqXDWBjXwpF/mpihWsh2FhAeweG0WybwlOfs5s8BnJTzz/FHQcaaWV1dbmfEJ8ew5N/NKAG/cazzb7Aq6rwqWg/cf9SvI133F5b4QW7QiJ44aSgx1ND8OnPIettjzFRid4PCeT07axR5LyqhlKS2zHhm36C9wz0edRTgQH38Z9G0XBpps1tuYX2aNo1Z2r1n66nj2SWTeOM67ZfcSzJIe7qNyT85lpUuMBNqhpnxunuAmUNRgWFgCgjrKCS/XiB7/bUWZofevNb41rJGazwuoXzMljb+EwzZo0VnjfdfvEuObJbN4FN/HlENMnnb4keXpDpgLpVJqHDG1J4b4MAPhlMk+VS9B4Cu4uxf0K6305CZ49zKRlDx3O8Z8HP+Z+kdAYzj9cCPa9rHXW72VllpXD91v/jnt/xt8pAATcxjM1MNa1NZULPmMcS8v+l6ggBVfYat0W/5Jvy+Jy4vaWHKf/ObYWyQ/X5DkU72271bhGye/Yg1jyA4d24h/h0ty0mzgElXHCjH97hXGoJvo5Vs5Z2mzx5/awB/LOPu7KPrCdR7gCwIdPTSf5kyHcqPddH56OJ/fxI3Vykzm+pWQwe1vlF3H4LaX6xTez5fiygfA/zu95JJW/D59THKIKOGrulhOO82vMqScFIyXT7CfJSuf8l1u38VkPw8ICCPTkEJTfz6w4Jh5iepv7W60xrjGzKce7JZunypX9nePbmdvYi/HINHfdyRXYC4nrzh5HhTt4x9d/NnNL/fAAc0n9VYkpvwFgWRIbMp/TvI74ZryGwEqneQ3vmeWpOnKSeOddIp4NyAXtWyPZmHo+xcb0+mCmqN/lxY2TjZ7cYFxyw4dmt70rOJtPSCothY2INRjFG9ZgWBSK35/lsEzpX9njeGATK6M5h0zFErqNd5InOBWAbF8OWfkf46y48jS52NLu4RxE6nbuDIifG0Vy9WkcPuq/jMMuHhmmao4fzD0Rs2ZwbH94PPMjxWqsvDteNSk56taJJfm1Suz5vBzLxIE5rbT550XA6w15Dvv9QZwEv3M3rzu9DHsgg8OZOh4A2t7CX1rQ18YpBWJlE3OonD6/PM44w02glMFt9m+ENRgWhaLrKE5AN/fnDt9fzrG3cPS4GawoqW0uy6/knbrfDp53cbYhh3IO8bRVAEBVf1akGWfZYPhow4wSvuIeirdrzyD52c95hgQAtKrEobMh9zNpot7dm34NG8bab5qkiRl+HF67v+dgkktqsytCLmBm2LvT7yW5rDYbxN+L7/+Kzpx4b/flc8Y1a83hSipXPZ+Rx28yjg0py7+tIS6OfS1WsJ3eFhbAx38xq+uo07w7VXo/V2gGdGhkqTjUjg/U2M2yV4o222KR2TTm9RSf0+cBJh9s4seJ9bOKE7tpGoNueHtzf7u3H5fN5gTxOrICed3pd2iNezHmACXPzcx4W+m1WP57IHtbF7JvDWnOjZJrznFfy9jK3Bh5yzqeJPjiXTxPAwBGpjCFSQUzalUgNrcxq9CG1NIMMDYa57gNbEjKwgKA8INQsx53aZ/NYPehVAkzhLIpm4n9PBJZWR+4h8tPK8/mXbU6qXGqA4ivwvmFaa04t7LsRW1sshbWOv4+f66F108x3iNjPp/zzEEOFx06w0SCFR7RDEY1s6z2wLQokqNKcZ4je6gWyV/HfFZFwZw6nKy/aRmXOv/6HXNzVdjDa5jt3cq4pnR2eRmEnOpmFVrmG/y9erczTnEPKJdmerstrMGwKBTRn/Met+zHnKQ9m8GNZzt+YwoOAAjS9L1ehRNwmJO0J1pwxVLQITMdGvkTx+W97+R1ndBmY4cM4M8R3psN2/31nzbew0ujZT/ahL2Ub57gUE7vjtwvcrxRPkrkDD92O86yUblnInc//3MBueb7BnCYK9KXPSO/g5ycP9aKmX5Lrzeb7F7rPZPkL4ZHubSm9DJmpVvKV+xNhcCkZ3EbWA/DwgI42oQrgToE8u7/1xhO9M7t8bFxjeraL23MKU6g9gtleu2H93QnecdeszwV6RziGFyGy3uXHefcyu4eXAlUQmsrCDxiBn8O9eDke9Xx3MV+eyU2Ml412Yup9q1Z3RU3iI9VHs5/35zOZctADFzFwe78Wfx28RfQ7GW+5s1B3By4O52pRADg1a+5pLiSi417X3xm/i7uf32oS9covrBJbwsLAEC5NbzbnBjCfRZdOjP3Uc+NPC0PALx/4NBNxAKO48/pyqW5wfs4D1KZ+8wAAH6ruepp4uPcV/H3YOanuiWVE8Gf1JxFcmNfkxpkbxb3g3Tax93OtT5lL+dIG27K2/eocUkEruHKKTWCQ3z7E9kQVmLbWSS0rMkGYPUxDt/tPcdeoY4eYeY0ve6PcKNj39dMyvmCMPZ4S+PYiRZsPEMnu3TJ4gNLb25h4UDUWC15/CF7BytOMh9Qha/NOQkHHmWDkVmDPYbwlax4lQ//NE9fZ1Ze3f4nezp1SzD3x3fJvEueWpvj+tHeHA5JV6Y38EjMQyTf3J6V5qQ+3HNyKofDXDnG7ECgx7u9SI6pZ+7mLxYJafzZvM+x53OuExvkf85wyOqfJjxFEQDgoZc2u5ag/n5dA+OYT6jrPFnFFpbe3MICeCaCSx+fe4J3yFvXcqin4g9m0nsf8wBiyFSOhy86xVQWP67mwL1HmtmH8fVH7On4D+Xegzq+3Jy2MJmHG02Yy683KL8B/N30C5LvfJCreqI78TVzSrCBqPSTqUS6zfmJ5JhfeB0Lb2calAspNT06K4rk/oN+IDn6YW4O/CyuNcnb1mtcVACmdeUelDerupZcqf2OWR68c4iZCHdHKADKehgWFkCnPwaQHP0aczZV3cW5A7PzAKhag9lox3zJ4aMcrbsZE7m7/OZq7OUAwBsPLCa5184HSfZ9nkMg3y7geEe/R3jiXt8DZrPHjaO1HomybACCYtmQlVvEpIpZB81S3UWLOYEf1Z6veagje2MXgqRmnGuZvp+bL/tU5fzDvGi26Jsrm6XReh7KVZy+0TQO0dPZSMlR3JMAAArBSURBVLvtHl3ZAUoWFgCAeTeOJ3ny1xy73nqaq3w8OnJ4CQBiu3MVTtTUWH5NaW66C/qTK2riJnA/BAC0fIiVee2hHAp75O/fSH4tgRvHXgrneaBeHmbS8ubuHHYpo00TXDCtOcnn6vK9OPiimawP28AJmWSNkeSpb7iBMOoCKL/rVOBGyN3L2Av8KYhzGqNm84hWVdukky83jSvEfPCXcU5BCF660ziWcBcXJpRysbejOOG/4GGI+peWgonIceCS1+iVBmBqw6uP4riu4rgmoHiuqziuCSie67pca6qslCq4EqAAiMiPcKytKDihlDIZOt0A/1qDcTkgIuuUUo0KP/PKojiuqziuCSie6yqOawKK57qK45r+S7i4Ib0WFhYWFv8ZWINhYWFhYVEkWIPhGkx+5uKB4riu4rgmoHiuqziuCSie6yqOa/rPwOYwLCwsLCyKBOthWFhYWFgUCdZgWFhYWFgUCdZgFAEiMktENjr/ixWRjc7jUSKSmudv4wu71iVc0+sicjjPe9+W528visgeEdkpIh2u1Jqc7z1SRGJEZLOIzBOREOfxq3avnO9/q/N+7BGRF67ke2vrqCgiv4rIDhHZJiJPOY+f9/u8QuuKFZEtzvde5zwWJiJLRWS38//mBKTLu6aaee7HRhE5IyJPX+179V+GzWG4CBH5EECSUupNEYkCsEgpde1VWMfrAJKVUh9ox68B8DWAGwBEAlgGoIZS6opwL4tIewC/KKWyROR9AFBKPX+V75UngF0A2sExNvpvAPcrpbYX+MLLs5ZyAMoppTaISBCA9QC6ArgH+XyfV3BdsQAaKaVO5Dk2AkCiUmq408iGKqWev0rr8wRwGEATAI/gKt6r/zKsh+ECRETgeLC/vtprKQBdAHyjlEpXSu0HsAcO43FFoJT6WSmVO0RiLYDiwC53A4A9Sql9SqkMAN/AcZ+uOJRSR5RSG5z/PgtgB4B8hn0UC3QBMM3572lwGLarhTYA9iql3HjCkvvDGgzX0BzAMaVU3mEDVUTkHxH5TUSan++FlwkDnaGfyXnCBeUB5GXAi8PVU0i9AeRltbta96o43ZP/wel11QeQy96Y3/d5paAA/Cwi60Ukd8B3hFLqCOAwdADCz/vqy4/7wBu1q3mv/rOwBsMJEVkmIlvz+S/vTvR+8I/2CIBKSqn6AIYAmCkiPFj68q3pMwDVANRzriN3VqjJA458hjJcvnXlnvMygCwAXzkPXdZ7VdiS8zl2VWOxIhIIYA6Ap5VSZ3D+7/NK4WalVAMAHQE8ISItrvD7nxci4gPgDgDfOg9d7Xv1n4Vlq3VCKdW2oL+LiBeAbgAa5nlNOoB057/Xi8heADUArMv3Ipd4TXnWNhHAIqcYB6Binj9XAGAOIriM6xKRngA6A2ijnEmyy32vCsFlvyeuQES84TAWXyml5gKAUupYnr/n/T6vCJRS8c7/J4jIPDjCeMdEpJxS6ogz95JQ4EUuHzoC2JB7j672vfovw3oYRUdbADFKqf8NOBCRMs5kHESkKoBoAPuuxGKcD3Au7gSw1fnvBQDuExFfEaniXJNrPNQXt65bATwP4A6lVEqe41ftXsGR5I4WkSrO3ep9cNynKw5nHmwSgB1KqY/yHD/f93kl1hTgTMBDRAIAtHe+/wIAPZ2n9QTw/ZVakwby7K/mvfqvw3oYRYceQwWAFgDeFJEsANkAHldKJV6h9YwQkXpwhFZiATwGAEqpbSIyG8B2OEJCT1ypCiknxgLwBbDUoRuxVin1OK7ivXJWbA0E8BMATwCTlVLbCnnZ5cLNAHoA2CLO8mwALwG4P7/v8wohAsA85/flBWCmUupHEfkbwGwR6QPgIIC7r+CaAAAi4g9HdVve+5Hvb9/i8sOW1VpYWFhYFAk2JGVhYWFhUSRYg2FhYWFhUSRYg2FhYWFhUSRYg2FhYWFhUSRYg2FhYWFhUSRYg2FhYWFhUSRYg/Efh4gkF/L3EBEZkEeOFJHvnP+udyHU0k566qEunn9YRN48z99jRaS0q+twN4hIKxG5KY88VUS653NeNXHQfhf43VpYuAprMCwKQwiA/xkMpVS8UipXSdUDcKVmEXyslHr1cr6Bk/6lWMK5tlYAbirkVCil9iql6l32RVn852ANhgUABxmeiCwXkQ3iGKSTSyQ4HEDujnWkOAYhbXVSbLwJ4F7n3+7VPQfneVHOf78sjgFGywDUzHNONRH50cmSukpEahVhraVE5GdxMN9+jjzkgiLykIj85VzT53noSPqIyC4RWSEiE0VkrPP4VBH5SER+BfC+kyZjsoj87bx+F+d5ns7P/7c4WFIfcx4vJyIrne+3VQpg4RWRZBF5R0Q2ichaEYlwHq/svPebnf+vlM/aZgF4HMBg53vlvk8LEfldRPbl521YWFxKWINhkYs0AHc6GUtbA/jQyXv0AhxzCOoppZ7NPdk5V+JVALOcf5t1vguLSEM4qFXqw0Hg2DjPnycAGKSUaghgKIBxRVjrawBWO5lvFwDIVbC1AdwLB/NqPTgoSB4UkUgArwBoCgfNhG6UagBoq5R6BsDLcAyAauy8DyOd/Ep94Bic1di5/r7i4Op6AMBPzverC2Ajzo8AOKhS6gJYCaCv8/hYANOVUtfDwe47Op+13QVgPByeVj2l1Crn38sBaAYH2ePwItw7C4sLRrF1wS2uOATAu+Kgtc6BY15ExCW6dnMA83LJCEVkgfP/gXCEWL518hgBDh6qwtACDsMDpdQPInLKebwNHGzCfzuv5wcHw+oNAH7L5a4SkW/hUMS5+DYP31Z7AHfk8ZRKwGGQ2gO4Ps8uPhgOAsW/AUwWBwPtfKVUQQYjA//PrLoeDuMFADfmfh4AXwIYcZ615Yf5SqkcANtzPRYLi8sFazAscvEggDIAGiqlMsUxsrOEi9fIAnuteV+fH2mZB4DTFxhvz+96AmCaUupFOihyZyHXOqdd4y6l1E7tGgKHJ/ST8aYOI9sJwJciMlIpNf0875OZS/cOh/dzvucv72c7d55zcpGurd3C4rLBhqQschEMIMFpLFoDqOw8fhZA0Hleo/8tFkADABCRBgCqOI+vBHCniPiJg0b7dgBwDg7aLyJ3O18jIlK3CGtdCYeBg4h0BJA7cW05gO4iEu78W5iIVIaD3r2liIQ6k8d3FXDtnwAMchoIiEj9PMf7Oz0JiEgNZ76jMhz3bSIctOUNirB+Hb/DEbKD83OtPs95BX0XFhaXHdZgWOTiKwCNRGQdHEorBgCUUicBrHEmdEdqr/kVwDW5SW84hgKFiYO2uz+AXc5rbIAjabvRec6qPNd4EEAfEdkEYBuKNmv7DTiSvRvgCBUddL7PdgDD4Bg1uhnAUgDllFKHAbwLxyjUZXBQvyed59pvAfAGsFlEtjplAPjC+boNzuOfw+EhtAKwUUT+gcMQjSrC+nU8CeAR55p7AHjqPOcthMPw5k16W1hcMVh6c4tiDxF5HUCyUuqDi7hGoFIq2elhzINjJsa8S7XG4ggRSVZKBV7tdVj8e2A9DAt3QDKAfnKexr0i4nWn57MVwH4A8y/JyoohnKXKGwEcK/RkCwsXYD0MC4tLDBH5E2a1Vw+l1JarsR4Li0sFazAsLCwsLIoEG5KysLCwsCgSrMGwsLCwsCgSrMGwsLCwsCgSrMGwsLCwsCgS/g92CReAnOvWNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# group by month and create ranking for each grid cell\n",
    "start_time   = time.time()\n",
    "\n",
    "# Call the rankdata function on 1D slices\n",
    "def ranking(array, axis):\n",
    "    return np.apply_along_axis(rankdata, axis, array, method='dense')\n",
    "\n",
    "rank2 = mrsos.groupby('time.month').reduce(ranking, axis=0)\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "rank2.sel(lon=20., method=\"nearest\").plot() # Roughly the longest path through Africa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating P...\n",
      "--- 107.12161040306091 seconds ---\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([11239942.,   369868.,   380397.,   341512.,   336447.,   352657.,\n",
       "          323388.,   330942.,   293787.,   285140.]),\n",
       " array([0.0046102 , 0.10369504, 0.20277988, 0.30186472, 0.40094956,\n",
       "        0.5000344 , 0.59911925, 0.69820409, 0.79728893, 0.89637377,\n",
       "        0.99545861]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEuVJREFUeJzt3X+QXWddx/H3x8TOCIWWIQsDaSEBUyA6FGmslQEtgpIUMTKDYwPYoVPtVGhHcRxbGcUfjEoHf4FtibGT6aDYqKVCgUB1dLBoW2069kcCFGJa2phqU2qpLdWS5usf90Svt5vs2eTu3t1n36+ZO7nnnOee832ym88+ee49z6aqkCS15VsmXYAkafwMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuWrSS7Epy5qTrkBYiw10LVpJ7krx+ZN87kvw9QFV9R1V9boZzrEpSSZbPYanSgmO4S8fAHxpaqAx3LVrDI/skpyfZkeSRJP+e5He7Zjd0fz6c5NEk35vkW5L8UpKvJnkgyUeSnDB03nO6Y19L8ssj1/nVJNck+ZMkjwDv6K59U5KHk9yf5LIkxw2dr5K8M8lXkvxnkvcleXH3mkeS/Plwe2kcDHe14oPAB6vqmcCLgT/v9n9f9+eJVXV8Vd0EvKN7vBZ4EXA8cBlAkrXAFcDbgOcBJwArR661EbgGOBH4KPAk8G5gBfC9wOuAd468Zj1wGnAG8AvAlu4aJwPfCWw6hr5LTzHRcE+ytRs57ezR9veS3NY9vpzk4fmoURP38W5E/HD3Nb/iMO2+CXx7khVV9WhV3XyEc74N+N2q2lNVjwK/CJzdTbG8BfhkVf19VT0BvBcYXYDppqr6eFUdrKrHq+rWqrq5qg5U1T3AHwLfP/KaS6vqkaraBewE/qq7/teBzwDf1f+vRJrZpEfuVzEY0cyoqt5dVa+oqlcAfwBcO5eFacH40ao68dCDp46IDzkPOAX4UpJbkvzwEc75fOCrQ9tfBZYDz+2O3XfoQFV9A/jayOvvG95IckqSTyX5t26q5jcZjOKH/fvQ88en2T7+CPVKszbRcK+qG4CHhvd1c5GfTXJrks8neek0L90EXD0vRWpRqKqvVNUm4DnApcA1SZ7OU0fdAPuAFw5tvwA4wCBw7wdOOnQgybcBzx693Mj2h4EvAWu6aaH3ADn63kjHbtIj9+lsAS6qqtOAn2fkv+FJXgisBv52ArVpgUry9iRTVXUQODRl9ySwHzjIYG79kKuBdydZneR4BiPtP6uqAwzm0t+U5FXdm5y/xsxB/QzgEeDRbjDy02PrmHSUFlS4d//QXgX8RZLbGMxdPm+k2dnANVX15HzXpwVtPbAryaMM3lw9u6r+q5tW+Q3gH7p5+zOArcAfM/gkzd3AfwEXAXRz4hcB2xiM4v8TeAD47yNc++eBt3Zt/wj4s/F3T5qdTPqXdSRZBXyqqr4zyTOBu6pqNNCH2/8z8K6qunGeStQS1g04HmYw5XL3pOuR+lpQI/eqegS4O8mPAWTg1EPHk7wEeBZw04RK1BKQ5E1JntbN2f82cCdwz2SrkmZn0h+FvJpBUL8kyd4k5zH4mNp5SW4HdjH4TPEhm4BtNen/bqh1Gxm86boPWMNgisfvOS0qE5+WkSSN34KalpEkjcfEFj1asWJFrVq1alKXl6RF6dZbb32wqqZmajexcF+1ahU7duyY1OUlaVFK8tWZWzktI0lNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDZrYHarHYtUln57Yte95/xsndm1J6suRuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBs0Y7km2Jnkgyc7DHE+SDyXZneSOJK8cf5mSpNnoM3K/Clh/hOMbgDXd43zgw8deliTpWMwY7lV1A/DQEZpsBD5SAzcDJyZ53rgKlCTN3jjm3FcC9w1t7+32SZImZBzhnmn21bQNk/OT7EiyY//+/WO4tCRpOuMI973AyUPbJwH7pmtYVVuqal1VrZuamhrDpSVJ0xlHuF8HnNN9auYM4OtVdf8YzitJOkrLZ2qQ5GrgTGBFkr3ArwDfClBVm4HtwFnAbuAbwLlzVawkqZ8Zw72qNs1wvIB3ja0iSdIx8w5VSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNahXuCdZn+SuJLuTXDLN8ROSfDLJ7Ul2JTl3/KVKkvqaMdyTLAMuBzYAa4FNSdaONHsX8IWqOhU4E/idJMeNuVZJUk99Ru6nA7urak9VPQFsAzaOtCngGUkCHA88BBwYa6WSpN76hPtK4L6h7b3dvmGXAS8D9gF3Aj9TVQdHT5Tk/CQ7kuzYv3//UZYsSZpJn3DPNPtqZPsNwG3A84FXAJcleeZTXlS1parWVdW6qampWRcrSeqnT7jvBU4e2j6JwQh92LnAtTWwG7gbeOl4SpQkzVafcL8FWJNkdfcm6dnAdSNt7gVeB5DkucBLgD3jLFSS1N/ymRpU1YEkFwLXA8uArVW1K8kF3fHNwPuAq5LcyWAa5+KqenAO65YkHcGM4Q5QVduB7SP7Ng893wf80HhLkyQdLe9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoN6hXuS9UnuSrI7ySWHaXNmktuS7Eryd+MtU5I0G8tnapBkGXA58IPAXuCWJNdV1ReG2pwIXAGsr6p7kzxnrgqWJM2sz8j9dGB3Ve2pqieAbcDGkTZvBa6tqnsBquqB8ZYpSZqNPuG+ErhvaHtvt2/YKcCzknwuya1JzpnuREnOT7IjyY79+/cfXcWSpBn1CfdMs69GtpcDpwFvBN4A/HKSU57yoqotVbWuqtZNTU3NulhJUj8zzrkzGKmfPLR9ErBvmjYPVtVjwGNJbgBOBb48liolSbPSZ+R+C7AmyeokxwFnA9eNtPkE8Joky5M8Dfge4IvjLVWS1NeMI/eqOpDkQuB6YBmwtap2JbmgO765qr6Y5LPAHcBB4Mqq2jmXhUuSDq/PtAxVtR3YPrJv88j2B4APjK80SdLR8g5VSWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNahXuCdZn+SuJLuTXHKEdt+d5MkkbxlfiZKk2Zox3JMsAy4HNgBrgU1J1h6m3aXA9eMuUpI0O31G7qcDu6tqT1U9AWwDNk7T7iLgY8ADY6xPknQU+oT7SuC+oe293b7/lWQl8GZg85FOlOT8JDuS7Ni/f/9sa5Uk9dQn3DPNvhrZ/n3g4qp68kgnqqotVbWuqtZNTU31rVGSNEvLe7TZC5w8tH0SsG+kzTpgWxKAFcBZSQ5U1cfHUqUkaVb6hPstwJokq4F/Bc4G3jrcoKpWH3qe5CrgUwa7JE3OjOFeVQeSXMjgUzDLgK1VtSvJBd3xI86zS5LmX5+RO1W1Hdg+sm/aUK+qdxx7WZKkY+EdqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWpQr3BPsj7JXUl2J7lkmuNvS3JH97gxyanjL1WS1NeM4Z5kGXA5sAFYC2xKsnak2d3A91fVy4H3AVvGXagkqb8+I/fTgd1VtaeqngC2ARuHG1TVjVX1H93mzcBJ4y1TkjQbfcJ9JXDf0Pbebt/hnAd8ZroDSc5PsiPJjv379/evUpI0K33CPdPsq2kbJq9lEO4XT3e8qrZU1bqqWjc1NdW/SknSrCzv0WYvcPLQ9knAvtFGSV4OXAlsqKqvjac8SdLR6DNyvwVYk2R1kuOAs4HrhhskeQFwLfATVfXl8ZcpSZqNGUfuVXUgyYXA9cAyYGtV7UpyQXd8M/Be4NnAFUkADlTVurkrW5J0JH2mZaiq7cD2kX2bh57/JPCT4y1NknS0vENVkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoOWTLmCxWXXJpyddwry75/1vnHQJ825SX+el+HetuWG4a0ZL8QfapPh3PX9a/0FquEtakib5g3Q+frA45y5JDeoV7knWJ7krye4kl0xzPEk+1B2/I8krx1+qJKmvGcM9yTLgcmADsBbYlGTtSLMNwJrucT7w4THXKUmahT4j99OB3VW1p6qeALYBG0fabAQ+UgM3Aycmed6Ya5Uk9dTnDdWVwH1D23uB7+nRZiVw/3CjJOczGNkDPJrkrllVO7ACePAoXreYLcU+w9Ls91LsMyyxfudS4Oj7/MI+jfqEe6bZV0fRhqraAmzpcc3DF5PsqKp1x3KOxWYp9hmWZr+XYp9hafZ7rvvcZ1pmL3Dy0PZJwL6jaCNJmid9wv0WYE2S1UmOA84Grhtpcx1wTvepmTOAr1fV/aMnkiTNjxmnZarqQJILgeuBZcDWqtqV5ILu+GZgO3AWsBv4BnDu3JV8bNM6i9RS7DMszX4vxT7D0uz3nPY5VU+ZGpckLXLeoSpJDTLcJalBCzbcl+KSBz36/Laur3ckuTHJqZOoc5xm6vNQu+9O8mSSt8xnfXOlT7+TnJnktiS7kvzdfNc4bj2+v09I8skkt3d9nsv37uZFkq1JHkiy8zDH5y7HqmrBPRi8cfsvwIuA44DbgbUjbc4CPsPgM/ZnAP846brnoc+vAp7VPd+wFPo81O5vGbxx/5ZJ1z1PX+sTgS8AL+i2nzPpuuehz+8BLu2eTwEPAcdNuvZj7Pf3Aa8Edh7m+Jzl2EIduS/FJQ9m7HNV3VhV/9Ft3szgfoLFrM/XGeAi4GPAA/NZ3Bzq0++3AtdW1b0AVbXY+96nzwU8I0mA4xmE+4H5LXO8quoGBv04nDnLsYUa7odbzmC2bRaT2fbnPAY/8RezGfucZCXwZmDzPNY11/p8rU8BnpXkc0luTXLOvFU3N/r0+TLgZQxugLwT+JmqOjg/5U3MnOXYQv1lHWNb8mAR6d2fJK9lEO6vntOK5l6fPv8+cHFVPTkY0DWhT7+XA6cBrwO+Dbgpyc1V9eW5Lm6O9OnzG4DbgB8AXgz8dZLPV9Ujc13cBM1Zji3UcF+KSx706k+SlwNXAhuq6mvzVNtc6dPndcC2LthXAGclOVBVH5+fEudE3+/vB6vqMeCxJDcApwKLNdz79Plc4P01mIzeneRu4KXAP81PiRMxZzm2UKdlluKSBzP2OckLgGuBn1jEI7hhM/a5qlZX1aqqWgVcA7xzkQc79Pv+/gTwmiTLkzyNwUqsX5znOsepT5/vZfA/FZI8F3gJsGdeq5x/c5ZjC3LkXgtvyYM517PP7wWeDVzRjWQP1CJeSa9nn5vTp99V9cUknwXuAA4CV1bVtB+nWwx6fq3fB1yV5E4G0xUXV9WiXgY4ydXAmcCKJHuBXwG+FeY+x1x+QJIatFCnZSRJx8Bwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEudZIsyPs+pKNhuKtZSVYl+VKSK5PsTPLRJK9P8g9JvpLk9CS/mmRLkr8CPpLkO5L8U7eO+h1J1nTn+rnuHDuT/Gy37+lJPt2tP74zyY9PtMPSEG9iUrOSrGJw5993AbsY3AJ/O4NF136Ewd2AtwFvAl5dVY8n+QPg5qr6aHeb/DJgLXAVg/W2A/wj8HYGa5Ovr6qf6q53QlV9fb76Jx2JI3e17u6qurNbOnYX8DfdwlR3Aqu6NtdV1ePd85uA9yS5GHhht//VwF9W1WNV9SiD9X1e053j9UkuTfIag10LieGu1v330PODQ9sH+b+1lR471KCq/pTBqP5x4PokP8D0y7LSLd52GoOQ/60k7x1v6dLRM9ylIUleBOypqg8xWLHv5cANwI8meVqSpzP45SGfT/J84BtV9SfAbzP4dWrSguCnA6T/78eBtyf5JvBvwK9X1UNJruL/1hW/sqr+OckbgA8kOQh8E/jpiVQsTcM3VCWpQU7LSFKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoP8BWtVVw9Sb8wwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# For the rest, you definitely don't need any loops. That's simple array comprehension will do.\n",
    "# I've replaced numpy functions with the equivalent in xarray so you keep xarray DataArray throughout.\n",
    "\n",
    "# Define sample size wrt length of climatology + year in question\n",
    "N = climend - climstart + 1\n",
    "den = N+0.33  # denominator is constant, so calculate it once only.\n",
    "\n",
    "# Empirical Tukey plotting position (Wilks 2011)\n",
    "print('Calculating P...')\n",
    "P = (rank2 - 0.33)/den\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "P.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating W...\n",
      "--- 284.2396261692047 seconds ---\n",
      "<xarray.DataArray 'mrsos' (time: 1740, lat: 64, lon: 128)>\n",
      "array([[[-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        ...,\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621]],\n",
      "\n",
      "       [[-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        ...,\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621]],\n",
      "\n",
      "       ...,\n",
      "\n",
      "       [[-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        ...,\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621]],\n",
      "\n",
      "       [[-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        ...,\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621],\n",
      "        [-0.004621, -0.004621, ..., -0.004621, -0.004621]]])\n",
      "Coordinates:\n",
      "  * time     (time) object 1861-01-16 12:00:00 ... 2005-12-16 12:00:00\n",
      "  * lon      (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2\n",
      "  * lat      (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 79.53 82.31 85.1 87.86\n",
      "    depth    float64 0.05\n",
      "    month    (time) int64 1 2 3 4 5 6 7 8 9 10 11 ... 2 3 4 5 6 7 8 9 10 11 12\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEWCAYAAABliCz2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFJtJREFUeJzt3X+0ZWV93/H3JzNhrSAK6owGAR1MR9MxS4xOCVqNWG2cwZrRtdIVRo2BkrJIxJXYZSuJjbVltZWa2GoEyYTOIiYWkhpiCI7BNGlKIpAwGH6NCI78kMkQGUFAlAQHvv1j7zHHw517z5277z135nm/1jrrnr33s/f+nn32/dzn7HPOc1NVSJIOfd8z7QIkSUvDwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBr4NWkh1JTp52HdLBwsDXspXkriSvG5t3WpK/AKiqF1XVn82xjTVJKsnKRSxVOigY+NIC+IdEBxMDXwet0VcASU5Msj3Jw0m+muRDfbOr+p8PJnkkycuTfE+Sf5/k7iT3Jfl4kiNHtvv2ftn9SX55bD/vT/LJJL+d5GHgtH7f1yR5MMm9ST6a5LCR7VWSn0vypSTfSHJukh/o13k4ye+OtpcWi4GvQ8WHgQ9X1dOAHwB+t5//o/3Po6rqiKq6Bjitv70GeD5wBPBRgCTrgAuAtwJHA0cCx4ztaxPwSeAo4BPA48C7gFXAy4HXAj83ts4G4GXAScC/A7b0+zgO+CFg8wIeuzSRqQZ+kq19D+uWCdr+9yQ39Lfbkzy4FDVq6j7V95wf7J/zC/bT7tvAP0qyqqoeqaprZ9nmW4EPVdUdVfUI8IvAqf3lmZ8A/rCq/qKqHgPeB4wPOHVNVX2qqp6oqker6vqquraq9lbVXcCvA68eW+e8qnq4qnYAtwCf7ff/EPAZ4IcnPyTSgZl2D/9iup7PnKrqXVX1kqp6CfBrwGWLWZiWjTdV1VH7bjy557zPGcALgC8muS7Jv5hlm88B7h6ZvhtYCTy7X3bPvgVV9S3g/rH17xmdSPKCJFck+dv+Ms9/oevtj/rqyP1HZ5g+YpZ6pUFMNfCr6irggdF5/bXNP0pyfZI/T/KDM6y6GbhkSYrUQaGqvlRVm4FnAecBn0zyFJ7cOwfYDTxvZPq5wF66EL4XOHbfgiTfBzxzfHdj0x8Dvgis7S8p/RKQA3800uKYdg9/JluAd1bVy4B3M/YSPsnzgOOBP51CbVqmkrwtyeqqegLYd7nvcWAP8ATdtfp9LgHeleT4JEfQ9ch/p6r20l2bf2OSV/RvpP5H5g7vpwIPA4/0HZSfHeyBSQNaVoHf//K9AvjfSW6guxZ69FizU4FPVtXjS12flrUNwI4kj9C9gXtqVf1df0nmPwOf698HOAnYCvwW3Sd47gT+DngnQH+N/Z3ApXS9/W8A9wF/P8u+3w28pW/7G8DvDP/wpIXLtP8BSpI1wBVV9UNJngbcVlXjIT/a/q+Bd1TV1UtUohrWd0IepLtcc+e065EWYln18KvqYeDOJP8SIJ0T9i1P8kLg6cA1UypRDUjyxiSH9+8B/ApwM3DXdKuSFm7aH8u8hC68X5hkV5Iz6D4yd0aSG4EddJ953mczcGlN+2WJDnWb6N7Y3Q2spbs85Dmng97UL+lIkpbGsrqkI0laPFMb+GnVqlW1Zs2aae1ekg5K119//deqavWBrDu1wF+zZg3bt2+f1u4l6aCU5O65W83MSzqS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktSIqX3TVgePNed8eir7vesDb5jKfqVDlT18SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqRFzBn6SrUnuS3LLfpYnyUeS7ExyU5KXDl+mJGmhJunhXwxsmGX5RmBtfzsT+NjCy5IkDW3OwK+qq4AHZmmyCfh4da4Fjkpy9FAFSpKGMcQ1/GOAe0amd/XzniTJmUm2J9m+Z8+eAXYtSZrUEIGfGebVTA2raktVra+q9atXrx5g15KkSQ0R+LuA40amjwV2D7BdSdKAhgj8y4G395/WOQl4qKruHWC7kqQBrZyrQZJLgJOBVUl2Af8B+F6AqroQ2AacAuwEvgWcvljFSpIO3JyBX1Wb51hewDsGq0iStCj8pq0kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaMVHgJ9mQ5LYkO5OcM8PyI5P8YZIbk+xIcvrwpUqSFmLOwE+yAjgf2AisAzYnWTfW7B3AF6rqBOBk4FeTHDZwrZKkBZikh38isLOq7qiqx4BLgU1jbQp4apIARwAPAHsHrVSStCCTBP4xwD0j07v6eaM+CvxjYDdwM/DzVfXE+IaSnJlke5Lte/bsOcCSJUkHYpLAzwzzamz69cANwHOAlwAfTfK0J61UtaWq1lfV+tWrV8+7WEnSgZsk8HcBx41MH0vXkx91OnBZdXYCdwI/OEyJkqQhTBL41wFrkxzfvxF7KnD5WJuvAK8FSPJs4IXAHUMWKklamJVzNaiqvUnOBq4EVgBbq2pHkrP65RcC5wIXJ7mZ7hLQe6rqa4tYtyRpnuYMfICq2gZsG5t34cj93cCPDVuaJGlIftNWkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjZgo8JNsSHJbkp1JztlPm5OT3JBkR5L/N2yZkqSFWjlXgyQrgPOBfw7sAq5LcnlVfWGkzVHABcCGqvpKkmctVsGSpAMzSQ//RGBnVd1RVY8BlwKbxtq8Bbisqr4CUFX3DVumJGmhJgn8Y4B7RqZ39fNGvQB4epI/S3J9krfPtKEkZybZnmT7nj17DqxiSdIBmSTwM8O8GpteCbwMeAPweuCXk7zgSStVbamq9VW1fvXq1fMuVpJ04Oa8hk/Xoz9uZPpYYPcMbb5WVd8EvpnkKuAE4PZBqpQkLdgkPfzrgLVJjk9yGHAqcPlYmz8AXpVkZZLDgR8Bbh22VEnSQszZw6+qvUnOBq4EVgBbq2pHkrP65RdW1a1J/gi4CXgCuKiqblnMwiVJ8zPJJR2qahuwbWzehWPTHwQ+OFxpkqQh+U1bSWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakREwV+kg1JbkuyM8k5s7T7J0keT/ITw5UoSRrCnIGfZAVwPrARWAdsTrJuP+3OA64cukhJ0sJN0sM/EdhZVXdU1WPApcCmGdq9E/g94L4B65MkDWSSwD8GuGdkelc/7zuSHAO8Gbhwtg0lOTPJ9iTb9+zZM99aJUkLMEngZ4Z5NTb9P4D3VNXjs22oqrZU1fqqWr969epJa5QkDWDlBG12AceNTB8L7B5rsx64NAnAKuCUJHur6lODVClJWrBJAv86YG2S44G/AU4F3jLaoKqO33c/ycXAFYa9JC0vcwZ+Ve1Ncjbdp29WAFurakeSs/rls163lyQtD5P08KmqbcC2sXkzBn1VnbbwsiRJQ/ObtpLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjDHxJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWrERIGfZEOS25LsTHLODMvfmuSm/nZ1khOGL1WStBBzBn6SFcD5wEZgHbA5ybqxZncCr66qFwPnAluGLlSStDCT9PBPBHZW1R1V9RhwKbBptEFVXV1VX+8nrwWOHbZMSdJCTRL4xwD3jEzv6uftzxnAZ2ZakOTMJNuTbN+zZ8/kVUqSFmySwM8M82rGhslr6AL/PTMtr6otVbW+qtavXr168iolSQu2coI2u4DjRqaPBXaPN0ryYuAiYGNV3T9MeZKkoUzSw78OWJvk+CSHAacCl482SPJc4DLgp6rq9uHLlCQt1Jw9/Kram+Rs4EpgBbC1qnYkOatffiHwPuCZwAVJAPZW1frFK1uSNF+TXNKhqrYB28bmXThy/2eAnxm2NEnSkPymrSQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9JjTDwJakRBr4kNcLAl6RGGPiS1AgDX5IaYeBLUiMMfElqhIEvSY0w8CWpEQa+JDXCwJekRhj4ktQIA1+SGmHgS1IjVk67AE1mzTmfnnYJkg5y9vAlqRH28OfJnrZ0aJjm7/JdH3jDVPZ7UAa+oSsdOvx9XjoHZeBLGpah2wYDX1pGDF4tJt+0laRGTNTDT7IB+DCwArioqj4wtjz98lOAbwGnVdXnB65VWjL2tHUomrOHn2QFcD6wEVgHbE6ybqzZRmBtfzsT+NjAdUqSFmiSHv6JwM6qugMgyaXAJuALI202AR+vqgKuTXJUkqOr6t7BK1Yz7GVLw5ok8I8B7hmZ3gX8yARtjgG+K/CTnEn3CgDgkSS3jSxeBXxtgnqmyRqHYY3DsMZhLHmNOW/eq4zW+LwD3e8kgZ8Z5tUBtKGqtgBbZtxJsr2q1k9Qz9RY4zCscRjWOIyWapzkUzq7gONGpo8Fdh9AG0nSFE0S+NcBa5Mcn+Qw4FTg8rE2lwNvT+ck4CGv30vS8jLnJZ2q2pvkbOBKuo9lbq2qHUnO6pdfCGyj+0jmTrqPZZ5+ALXMeKlnmbHGYVjjMKxxGM3UmO6DNZKkQ53ftJWkRhj4ktSIJQ38JM9I8sdJvtT/fPoMbV6Y5IaR28NJfqFf9v4kfzOy7JRp1Ni3uyvJzX0d2+e7/mLXmOS4JP83ya1JdiT5+ZFli3Ick2xIcluSnUnOmWF5knykX35TkpdOuu5QJqjxrX1tNyW5OskJI8tmfM6nUOPJSR4aef7eN+m6S1jjvx2p75Ykjyd5Rr9sqY7j1iT3JbllP8uXw/k4V43Dno9VtWQ34L8B5/T3zwHOm6P9CuBvgef10+8H3r0cagTuAlYt9DEuVo3A0cBL+/tPBW4H1i3Wceyfqy8DzwcOA27ct7+RNqcAn6H73sZJwF9Ouu4S1vgK4On9/Y37apztOZ9CjScDVxzIuktV41j7NwJ/upTHsd/PjwIvBW7Zz/Kpno8T1jjo+bjUl3Q2Ab/Z3/9N4E1ztH8t8OWquntRq/pu861x6PUH2UdV3Vv9AHZV9Q3gVrpvPy+W7wzBUVWPAfuG4Bj1nSE4qupa4KgkR0+47pLUWFVXV9XX+8lr6b5TspQWciyWzXEcsxm4ZBHqmFVVXQU8MEuTaZ+Pc9Y49Pm41IH/7Oo/n9//fNYc7U/lySfK2f3Lm62LcblkHjUW8Nkk16cbMmK+6y9FjQAkWQP8MPCXI7OHPo77G15jkjaTrDuE+e7nDLoe4D77e86HNGmNL09yY5LPJHnRPNddqhpJcjiwAfi9kdlLcRwnMe3zcb4WfD4O/g9Qkvwf4PtnWPTeeW7nMODHgV8cmf0x4Fy6B3ou8KvAv5pSjf+0qnYneRbwx0m+2P+1HsSAx/EIul+2X6iqh/vZgxzH8V3NMG/SITgmGppjABPvJ8lr6H7BXjkye1Gf83nU+Hm6y5yP9O+/fIpupNpldxzpLud8rqpGe7FLcRwnMe3zcWJDnY+DB35VvW5/y5J8Nf0omv1Lp/tm2dRG4PNV9dWRbX/nfpLfAK6YVo1Vtbv/eV+S36d7GXgVMJ/HuKg1JvleurD/RFVdNrLtQY7jmIUMwXHYBOsOYaIhQJK8GLgI2FhV9++bP8tzvqQ1jvzhpqq2JbkgyapJ1l2qGkc86VX6Eh3HSUz7fJzIkOfjUl/SuRz46f7+TwN/MEvbJ13368NtnzcDM76zvUBz1pjkKUmeuu8+8GMjtcznMS5mjQH+J3BrVX1obNliHMeFDMExybpDmHM/SZ4LXAb8VFXdPjJ/tud8qWv8/v75JcmJdL/H90+y7lLV2Nd2JPBqRs7PJTyOk5j2+Tinwc/HxXjneZZ3pJ8J/Anwpf7nM/r5zwG2jbQ7nO4EPnJs/d8CbgZuonsCjp5GjXTv3t/Y33YA751r/SnU+Eq6l6E3ATf0t1MW8zjSferhdrpPOLy3n3cWcFZ/P3T/TOfL/f7Xz7buIp2Dc9V4EfD1kWO2fa7nfAo1nt3XcCPdG3mvWG7HsZ8+Dbh0bL2lPI6X0A3R/m263vwZy/B8nKvGQc9Hh1aQpEb4TVtJaoSBL0mNMPAlqREGviQ1wsCXpEYY+JLUCANf6iUZ/Jvn0nJi4OuQlWRNki8muSjdmOyfSPK6JJ9L978ETkz3vwG2JPks8PEkL0ryV+nGGL8pydp+W/+m38Yt+Yf/z/CUJJ/uBzG7JclPTvUBS3Pwi1c6ZKUbJXQn3UihO+i+Mn8j3bcZfxw4ne7bi28EXllVjyb5NeDaqvpE/7X6FcA64GK6MdNDN+ro2+i+7bihqv51v78jq+qhpXp80nzZw9eh7s6qurmqnqAL/T+prpdzM7Cmb3N5VT3a378G+KUk76EbkfJRumEqfr+qvllVj9CNbfKqfhuvS3JeklcZ9lruDHwd6v5+5P4TI9NP8A+jxX5zX4Oq+l90vf9HgSuT/DNmHi6X6gazehld8P/XjPyrQWk5MvClEUmeD9xRVR+hG1juxXRDzr4pyeH9yIRvBv48yXOAb1XVbwO/Qvev6qRly08lSN/tJ4G3Jfk23f9T/k9V9UCSi4G/6ttcVFV/neT1wAeTPEE32uHPTqViaUK+aStJjfCSjiQ1wsCXpEYY+JLUCANfkhph4EtSIwx8SWqEgS9Jjfj/9kbcwyCsAeMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "print('Calculating W...')\n",
    "W1 = P.where(P>0.5,xr.ufuncs.log(1-P)) # Contains the value we want where P<=0.5, else contains P\n",
    "W  = W1.where(P<=0.5,xr.ufuncs.sqrt(-2. * xr.ufuncs.log(P))) # We take W1 values where P<=0.5, else the alternate calcultion.\n",
    "W.plot()\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "print(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating EDDI...\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(array([   23829.,    48203.,    47981.,    47575.,   120386.,   169270.,\n",
       "          472297., 11738625.,   211420.,  1374494.]),\n",
       " array([-18.1018103 , -16.07346968, -14.04512906, -12.01678843,\n",
       "         -9.98844781,  -7.96010719,  -5.93176656,  -3.90342594,\n",
       "         -1.87508532,   0.15325531,   2.18159593]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFgFJREFUeJzt3X20ZXV93/H3xxmxVVSMMxgenTEFE7BoZIpotWJ9GjBmtI0VxCAGO4tEXIldrooxMbastFITW43gdEJZxISCqUEywTGYpjUkCgmDPAwjT8ODMg6BUYsIEnGcb//Ye8zJ5dx7z7333If5+X6tddY9e+/fOft79rnzmd/9nb1/J1WFJKktT1jsAiRJ42e4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHDXPivJtiQnLHYd0lJkuGvJSnJPkldNWHd6kr8CqKqjq+oL0zzHqiSVZPk8liotOYa7NAf+p6GlynDXPmuwZ5/kuCRbkjyU5P4kH+mbXdX/fDDJw0lenOQJSX4tyVeTPJDkk0mePvC8p/Xbvpnk1yfs54NJPp3kD5I8BJze7/vqJA8muS/Jx5PsN/B8leSXktyR5DtJzknyE/1jHkryh4PtpXEw3NWKjwIfraqnAT8B/GG//l/0Pw+oqv2r6mrg9P72CuA5wP7AxwGSHAWcD5wKHAQ8HThkwr7WAZ8GDgAuBn4AvBtYAbwYeCXwSxMesxY4Fjge+PfAxn4fhwHPA06Zw2uXHmdRwz3JhX3P6eYR2v7XJDf0t9uTPLgQNWrRXd73iB/s3/PzJ2n3feCfJFlRVQ9X1TVTPOepwEeq6q6qehh4H3ByP8Tyc8CfVNVfVdVjwAeAiRMwXV1Vl1fVnqp6tKquq6prqmp3Vd0D/Hfg5RMec25VPVRV24Cbgc/3+/828Dngp0c/JNL0FrvnfhFdj2ZaVfXuqnpBVb0A+B3gsvksTEvGG6rqgL03Ht8j3usM4Ejg1iTXJvmZKZ7zYOCrA8tfBZYDz+q33bt3Q1V9F/jmhMffO7iQ5MgkVyT5236o5j/R9eIH3T9w/9Ehy/tPUa80Y4sa7lV1FfCtwXX9WOSfJrkuyV8m+ckhDz0FuGRBitQ+oaruqKpTgAOBc4FPJ3kKj+91A+wEnj2wfDiwmy5w7wMO3bshyT8GnjlxdxOWPwHcChzRDwv9KpDZvxpp7ha75z7MRuBdVXUs8B4m/Bme5NnAauD/LEJtWqKSvDXJyqraA+wdsvsBsAvYQze2vtclwLuTrE6yP11P+1NVtZtuLP31SV7Sf8j5H5g+qJ8KPAQ83HdGfnFsL0yapSUV7v0/tJcA/yvJDXRjlwdNaHYy8Omq+sFC16clbS2wLcnDdB+unlxVf9cPq/wm8MV+3P544ELg9+nOpLkb+DvgXQD9mPi7gEvpevHfAR4AvjfFvt8DvKVv+7vAp8b/8qSZyWJ/WUeSVcAVVfW8JE8DbquqiYE+2P564J1V9aUFKlE/wvoOx4N0Qy53L3Y90qiWVM+9qh4C7k7yJoB0nr93e5LnAs8Arl6kEvUjIMnrkzy5H7P/LWArcM/iViXNzGKfCnkJXVA/N8mOJGfQnaZ2RpIbgW105xTvdQpwaS32nxtq3Tq6D113AkfQDfH4O6d9yqIPy0iSxm9JDctIksZj0SY9WrFiRa1atWqxdi9J+6TrrrvuG1W1crp2ixbuq1atYsuWLYu1e0naJyX56vStHJaRpCZNG+7TTe6V5NQkN/W3Lw2euihJWhyj9NwvYurJve4GXl5VxwDn0E0fIElaRNOOuVfVVf1VpJNtH7xS9BoGJl2SJC2OcY+5n0E3N/VQSdb335azZdeuXWPetSRpr7GFe5JX0IX7eydrU1Ubq2pNVa1ZuXLaM3kkSbM0llMhkxwDXACcWFUTv9hAkrTA5txzT3I43bci/XxV3T73kiRJczVtz72f3OsEYEWSHcBvAE8EqKoNdN8x+Uzg/CQAu6tqzXwVLEma3ihny0z5rexV9Q7gHWOrSNKCW3X2Zxdt3/d86HWLtu+WeYWqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBk0b7kkuTPJAkpsn2Z4kH0uyPclNSV44/jIlSTMxSs/9ImDtFNtPBI7ob+uBT8y9LEnSXEwb7lV1FfCtKZqsAz5ZnWuAA5IcNK4CJUkzN44x90OAeweWd/TrHifJ+iRbkmzZtWvXGHYtSRpmHOGeIetqWMOq2lhVa6pqzcqVK8ewa0nSMOMI9x3AYQPLhwI7x/C8kqRZGke4bwJO68+aOR74dlXdN4bnlSTN0vLpGiS5BDgBWJFkB/AbwBMBqmoDsBk4CdgOfBd4+3wVK0kazbThXlWnTLO9gHeOrSJJ0px5haokNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KCRwj3J2iS3Jdme5Owh25+e5E+S3JhkW5K3j79USdKopg33JMuA84ATgaOAU5IcNaHZO4GvVNXzgROA306y35hrlSSNaJSe+3HA9qq6q6oeAy4F1k1oU8BTkwTYH/gWsHuslUqSRjZKuB8C3DuwvKNfN+jjwE8BO4GtwC9X1Z6JT5RkfZItSbbs2rVrliVLkqYzSrhnyLqasPxa4AbgYOAFwMeTPO1xD6raWFVrqmrNypUrZ1ysJGk0o4T7DuCwgeVD6Xrog94OXFad7cDdwE+Op0RJ0kyNEu7XAkckWd1/SHoysGlCm68BrwRI8izgucBd4yxUkjS65dM1qKrdSc4CrgSWARdW1bYkZ/bbNwDnABcl2Uo3jPPeqvrGPNYtSZrCtOEOUFWbgc0T1m0YuL8TeM14S5MkzZZXqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWrQSOGeZG2S25JsT3L2JG1OSHJDkm1J/mK8ZUqSZmL5dA2SLAPOA14N7ACuTbKpqr4y0OYA4HxgbVV9LcmB81WwJGl6o/TcjwO2V9VdVfUYcCmwbkKbtwCXVdXXAKrqgfGWKUmaiVHC/RDg3oHlHf26QUcCz0jyhSTXJTlt2BMlWZ9kS5Itu3btml3FkqRpjRLuGbKuJiwvB44FXge8Fvj1JEc+7kFVG6tqTVWtWbly5YyLlSSNZtoxd7qe+mEDy4cCO4e0+UZVPQI8kuQq4PnA7WOpUpI0I6P03K8FjkiyOsl+wMnApglt/hh4WZLlSZ4MvAi4ZbylSpJGNW3Pvap2JzkLuBJYBlxYVduSnNlv31BVtyT5U+AmYA9wQVXdPJ+FS5ImN8qwDFW1Gdg8Yd2GCcsfBj48vtIkSbPlFaqS1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaNFK4J1mb5LYk25OcPUW7f5bkB0l+bnwlSpJmatpwT7IMOA84ETgKOCXJUZO0Oxe4ctxFSpJmZpSe+3HA9qq6q6oeAy4F1g1p9y7gj4AHxlifJGkWRgn3Q4B7B5Z39Ot+KMkhwBuBDVM9UZL1SbYk2bJr166Z1ipJGtEo4Z4h62rC8n8D3ltVP5jqiapqY1Wtqao1K1euHLVGSdIMLR+hzQ7gsIHlQ4GdE9qsAS5NArACOCnJ7qq6fCxVSpJmZJRwvxY4Islq4OvAycBbBhtU1eq995NcBFxhsEvS4pk23Ktqd5Kz6M6CWQZcWFXbkpzZb59ynF2StPBG6blTVZuBzRPWDQ31qjp97mVJkubCK1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1aKRwT7I2yW1Jtic5e8j2U5Pc1N++lOT54y9VkjSqacM9yTLgPOBE4CjglCRHTWh2N/DyqjoGOAfYOO5CJUmjG6XnfhywvaruqqrHgEuBdYMNqupLVfX/+sVrgEPHW6YkaSZGCfdDgHsHlnf06yZzBvC5uRQlSZqb5SO0yZB1NbRh8gq6cH/pJNvXA+sBDj/88BFLlCTN1Cg99x3AYQPLhwI7JzZKcgxwAbCuqr457ImqamNVramqNStXrpxNvZKkEYwS7tcCRyRZnWQ/4GRg02CDJIcDlwE/X1W3j79MSdJMTDssU1W7k5wFXAksAy6sqm1Jzuy3bwA+ADwTOD8JwO6qWjN/ZUuSpjLKmDtVtRnYPGHdhoH77wDeMd7SJEmz5RWqktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUEjnQopSa1ZdfZnF23f93zodfO+D3vuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBTvkrLSGLOQ2t2mLPXZIaZLhLUoMMd0lqkOEuSQ3yA1VJi8oPkefHSOGeZC3wUWAZcEFVfWjC9vTbTwK+C5xeVV8ec63SgjFwtK+bNtyTLAPOA14N7ACuTbKpqr4y0OxE4Ij+9iLgE/1PadYMWGn2Rum5Hwdsr6q7AJJcCqwDBsN9HfDJqirgmiQHJDmoqu4be8X4j16SpjNKuB8C3DuwvIPH98qHtTkE+AfhnmQ9sL5ffDjJbTOqdmZWAN+Yx+efraVY11KsCaxrJpZiTbA061r0mnLu0NWj1vXsUfYxSrhnyLqaRRuqaiOwcYR9zlmSLVW1ZiH2NRNLsa6lWBNY10wsxZpgada1FGuC8dc1yqmQO4DDBpYPBXbOoo0kaYGMEu7XAkckWZ1kP+BkYNOENpuA09I5Hvj2fI23S5KmN+2wTFXtTnIWcCXdqZAXVtW2JGf22zcAm+lOg9xOdyrk2+ev5JEtyPDPLCzFupZiTWBdM7EUa4KlWddSrAnGXFe6E1wkSS1x+gFJapDhLkkN2qfDPcmbkmxLsifJmoH1pya5YeC2J8kLhjz+g0m+PtDupHmua1WSRwf2t2GSx/9Ykj9Lckf/8xnzWNOrk1yXZGv/819O8vgFPVb9tvcl2Z7ktiSvneTxYz9WE57/UwOv+Z4kN0zS7p7+GN6QZMs4a5hkfyO9H0nW9sdve5Kz57mmDye5NclNST6T5IBJ2i3IsZrutfcngHys335TkhfOVy39/g5L8n+T3NL/zv/ykDYnJPn2wPv6gVnvsKr22RvwU8BzgS8AayZp80+BuybZ9kHgPQtVF7AKuHmEx/8X4Oz+/tnAufNY008DB/f3nwd8fYkcq6OAG4EnAauBO4FlC3Gspqj1t4EPTLLtHmDFfO17Nu8H3QkQdwLPAfbrj+dR81jTa4Dl/f1zJ3svFuJYjfLa6U4C+RzddTrHA389zzUdBLywv/9U4PYhNZ0AXDGO/e3TPfequqWqprvK9RTgkoWoZ68R65rKOuD3+vu/B7xhvmqqquurau81CduAf5TkSXPd31zrojsGl1bV96rqbrozsY6bpN1Yj9UwSQL8Gxb4d2mOfjh1SFU9BuydOmReVNXnq2p3v3gN3fUui2WU1/7DaVOq6hrggCQHzVdBVXVf9RMqVtV3gFvoruSfF/t0uI/ozUz9D/Ks/k+yC8f9J/0kVie5PslfJHnZJG2eVf11Av3PAxegLoB/DVxfVd+bZPtCHqvJprSYaKGO1cuA+6vqjkm2F/D5fmhr/SRtxm2692PUYzgffoGuVzzMQhyrUV77oh2fJKvo/mr+6yGbX5zkxiSfS3L0bPex5OdzT/K/gR8fsun9VfXH0zz2RcB3q+rmSZp8AjiH7pftHLo/u39hHuu6Dzi8qr6Z5Fjg8iRHV9VDo+xznmra+9ij6f6Ufs0kTRb6WI00pcU4jFjfdH8B/vOq2pnkQODPktxaVVfNV12M9n6M/RiOcqySvB/YDVw8ydOM/VgNK3XIullNmzJuSfYH/gj4lSH/9r8MPLuqHu4/R7mcbrbdGVvy4V5Vr5rDw09min+QVXX/3vtJfhe4Yj7r6nvE3+vvX5fkTuBIYOKHSvenn1Wz/zPxgfmqCSDJocBngNOq6s5JnntBjxWjT2kxq2M1k/qSLAf+FXDsFM+xs//5QJLP0A0LzCmwRj1uU7wfY58WZIRj9TbgZ4BXVj+IPOQ5xn6shliS06YkeSJdsF9cVZdN3D4Y9lW1Ocn5SVZU1YwnOmt2WCbJE4A30Y21TdZmcHztjcBkPfxx1bQy3fz4JHkO3f/Idw1pugl4W3//bcCUve451nQA8FngfVX1xSnaLeixojsGJyd5UpLVdMfqbyZpN9/H6lXArVW1Y9jGJE9J8tS99+n++pnv36VR3o9Rpg4ZZ01rgfcCP1tV352kzUIdqyU3bUr/uc3/AG6pqo9M0ubH+3YkOY4uo785qx3O56fD832j+6XeQdcbvh+4cmDbCcA1Qx5zAf1ZGcDvA1uBm+je6IPmsy66Me1tdJ/cfxl4/SR1PRP4c+CO/uePzWNNvwY8AtwwcDtwsY9Vv+39dGc83AacuFDHakiNFwFnTlh3MLC5v/+c/j29sX9/378Av/tD34/Buvrlk+jOyrhzvuui+9D73oHfow2LeayGvXbgzL3vJd2wzHn99q1McsbdGOt5Kd2wz00Dx+ikCTWdNZAR1wAvme3+nH5AkhrU7LCMJP0oM9wlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEu9/kpUqQmGu5qVbv78W5NckOTmJBcneVWSL6ab//24dPOib0zyeeCTSY5O8jf9XNo3JTmif65/1z/HzUl+pV/3lCSf7Sd5ujnJmxf1BUsDvIhJzepn3ttON/veNrpL0m8EzgB+lu6L3G8AXg+8tKoeTfI7dFc2X9xftr6Mbm75i+jm/A7dTH5vpbvacm1V/dt+f0+vqm8v1OuTpmLPXa27u6q2VtUeuoD/8+p6NFvpvjwFYFNVPdrfvxr41STvpZud71G6y8Y/U1WPVNXDwGV0UwBvBV6V5NwkLzPYtZQY7mrd4Nz0ewaW9/D3s6I+srdBVf1Pul79o8CV6b52cNjUsFTV7XQzRW4F/nPm8pVo0pgZ7tKAfrbOu6rqY3QTch1DNx3tG5I8uZ/J8I3AXyY5mO77Av4A+C1gXr+DU5oJzw6Q/qE3A29N8n3gb4H/WFXfSnIRfz/l8AVVdX26L+3+cJI9wPeBX1yUiqUh/EBVkhrksIwkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ36/9zoLPZ0CdfBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Inverse normal approximation (Vincente-Serrano et al 2010)\n",
    "print('Calculating EDDI...')\n",
    "mrsos_est1 = P.where(P>0.5,W - (C0 + C1 * W + C2 * W**2.) / (1. + d1 * W + d2 * W**2. + d3 * W**3.))\n",
    "mrsos_est = mrsos_est1.where(P<=0.5, -1. * (W - (C0 + C1 * W + C2 * W**2.) / (1. + d1 * W + d2 * W**2. + d3 * W**3.)))\n",
    "mrsos_est.plot()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-19.04]",
   "language": "python",
   "name": "conda-env-analysis3-19.04-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}