{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Reformat-for-Florida-USGS\" data-toc-modified-id=\"Reformat-for-Florida-USGS-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Reformat for Florida USGS</a></span><ul class=\"toc-item\"><li><span><a href=\"#NoData-mask\" data-toc-modified-id=\"NoData-mask-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>NoData mask</a></span></li><li><span><a href=\"#Solar-Zenith-Angles\" data-toc-modified-id=\"Solar-Zenith-Angles-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Solar Zenith Angles</a></span></li><li><span><a href=\"#Pixel-index\" data-toc-modified-id=\"Pixel-index-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Pixel index</a></span></li><li><span><a href=\"#(Potentially-useful-code)\" data-toc-modified-id=\"(Potentially-useful-code)-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>(Potentially useful code)</a></span></li><li><span><a href=\"#Write-text-outputs\" data-toc-modified-id=\"Write-text-outputs-1.5\"><span class=\"toc-item-num\">1.5&nbsp;&nbsp;</span>Write text outputs</a></span></li></ul></li><li><span><a href=\"#Albedo\" data-toc-modified-id=\"Albedo-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Albedo</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results\n",
    "## Reformat for Florida USGS\n",
    "*Write to text file structure required by the Florida USGS evapotranspiration model.*\n",
    "\n",
    "Glob for the result dataset(s) and print the header for one of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Black sky albedo outputs for 2018:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['result/black_albedo_MCD43A1.2018_1.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_2.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_3.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_4.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_5.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_6.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_7.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_8.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_9.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_10.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_11.nc',\n",
       " 'result/black_albedo_MCD43A1.2018_12.nc']"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import warnings; warnings.filterwarnings('ignore')\n",
    "sort2 = lambda ls: sorted(ls, key=lambda x: int(x.split(\"_\")[-1][:-3]))\n",
    "\n",
    "bsa_result = sort2(glob.glob(\"result/black_albedo*.nc\"))\n",
    "wsa_result = sort2(glob.glob(\"result/white_albedo*.nc\"))\n",
    "alb_result = sort2(glob.glob(\"result/blue_albedo*.nc\"))\n",
    "\n",
    "print(\"Black sky albedo outputs for 2018:\"); bsa_result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the twelve black sky albedo netCDFs as a multidataset and print the header:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "np.set_printoptions(suppress=True)\n",
    "\n",
    "ds = xr.open_dataset(bsa_result[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NoData mask\n",
    "Make an array of booleans that represent masked pixels by reading the first timestep of the input dataset and checking for `np.nan`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False, False, ..., False, False, False],\n",
       "       [False, False, False, ..., False, False, False],\n",
       "       [False, False, False, ..., False, False, False],\n",
       "       ...,\n",
       "       [False, False, False, ..., False, False, False],\n",
       "       [False, False, False, ..., False, False, False],\n",
       "       [False, False, False, ..., False, False, False]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# temporarily open the input dataset that we processed previously\n",
    "with xr.open_dataset(\"data/MCD43A1.2018.nc\") as tmp:\n",
    "    \n",
    "    # select timestep 1 of band 1 of the input dataset\n",
    "    tmp_array = tmp[\"BRDF_Albedo_Parameters_Band1\"][0].sel(Num_Parameters=0)\n",
    "    \n",
    "    # get the inverse of the boolean array of invalid pixels\n",
    "    mask_array = np.logical_not(np.isnan(tmp_array.data))\n",
    "\n",
    "mask_array"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make an `xarray.DataArray` for the mask and add it to the results dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fade0f08f60>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApsAAAJRCAYAAAAK4Xy2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdebxkZX3v+8+PQcChmTQJDXglgt6rHFsjB4h6E4WEbnNyxESNRI/ikMPRYNSY5KgnA3HKCTHROAUuURS85nYI0SMxYAcjxBiRSW0mBzo40GkTLzYgRkW6+3f+qKeguqhh1d61qmrV+rxfr/XaVavWs/aza9eu+u7fep61IjORJEmS6rDHvDsgSZKk5WXYlCRJUm0Mm5IkSaqNYVOSJEm1MWxKkiSpNoZNSZIk1cawKUmStEQi4tyI+FZE3DDk8YiId0bEloi4LiJ+os7+GDYlSZKWyweADSMefzpwVFlOA86qszOGTUmSpCWSmZ8Cto/Y5GTg/Oz4LHBARBxSV3/2qmvHy+YBsU/uy4Pm3Y2xHvW47827CwvjK9c9cN5daKxleh35OpDa4y5uvy0zHzbr77v+aQ/Kb2/fObPvd+11d98I/KBn1TmZec4EuzgUuLXn/tay7ptT6N79GDYr2pcHcVycOO9ujLVp0+Z5d0FaGOvXruO4mHcvJM3KJ/LCr8/j+357+06u2vTwmX2/PQ+5+QeZecwqdjHonbG265d7GH2JbNpm0JQkSWNtBQ7vuX8YsK2ub2Zlc0kYNCVJmo8EdrFr3t2YxEXAKyJiI3AccGdm1nIIHQybS8GgKd3f+rXr5t0FSZqLiPj/gKcCD42IrcAZwN4AmXk2cDHwc8AW4HvAi+vsj2GzwQyZ0mAGTUmzlezMxalsZuYvj3k8gdNn1B3HbDaVQVMazr8PSVocVjYbyA9SSZIWR2fMZm2TuRvPymbDGDSl8TyMLkmLw8pmgxg0JUlaTA2bjT5Ths0GMGRKkqSmMmwuOIOmNBkPoUuatSTZmY7ZHMYxm5KWhkFTkhaPYVOSJEm1MWxKWgpWNaX52rRt825Dv9o2DGwXObOlaRyzKanxDJrSfLU5ZGo8K5sLzg9RSdIiGxYu2/T5lcBOcmZL0xg2JTVamz7QpHnrD5ajqpj9h9XVXh5GlyRJY3WD46Ztm1m/dl3lINm/3bL+g9jEsZSzYmVTUmMt64eWtOhWU7G04tk+VjYXnH+QkqR5muXnULdq2jQJntR9BMOmpMZq6geT1AR1h0yLKe1h2JTUWAZNqR4GwcntmncHFphhU5IkLYxBE5C6/1i2ZbLRsjFsLrhhf2CSPL+fVIdF+LwZd4qlbhgdFEjnIRt6/stZMWxKWjrT/LA0uKpNFiFojjPsb3LTts3seciMO6NKDJsNYYVTmo86/uYMsFpETfl8GX1E4+bZdqYrYaeFzaEMm5I0Y9P6UDe0alqaEjQH8e9g8Rk2JamhDK2ahiYHzUWROBt9FMOmJLXcJGHDYKpF4uuxGQybDeJ/n5LmbRrvQwaExdHkzxVfR81h2GyIJr8hSFKv1byfGTAEi/g6CHYS8+7EwjJsNoBBU5I6DKry99g8tYbNiNgX+BSwT/leF2bmGT2P/ybwVuBhmXlbWfd64KXATuCVmbmprH8i8AFgP+Bi4FWZmRGxD3A+8ETg28BzM/Nrpc2pwO+Ub/fmzDyvrD8C2AgcBHwOeEFm/rCmp2HFDJmSND0G1d018TNmUX8PCezy1EdD1V3ZvBs4ITO/GxF7A5+OiEsy87MRcTjws8A3uhtHxGOAU4DHAmuBT0TEozJzJ3AWcBrwWTphcwNwCZ1gentmHhkRpwBnAs+NiIOAM4Bj6LwOro2IizLz9rLN2zNzY0ScXfZxVs3PhSSpoVYbzBY1JDWJz2Fz1Ro2MzOB75a7e5elm/3fDvx34KM9TU4GNmbm3cBXI2ILcGxEfA1Yk5lXAETE+cAz6YTNk4HfL+0vBN4dEQGsBy7NzO2lzaXAhojYCJwAPK+0Oa+0X5iw2cT/NiVJw630fd2A1RyO2Ryu9jGbEbEncC1wJPCezLwyIp4B/Etmbu7kwnsdSqdy2bW1rLun3O5f321zK0Bm7oiIO4GDe9f3tTkYuCMzdwzYV3/fT6NTTWVfHjjBT71yBk1JUteoz4TVBNGmfdYYuput9rBZDoE/PiIOAD4SEY8Dfhs4acDmg/4tyBHrV9Jm1L52X5l5DnAOwJo4qPbRGE3745ckzd5qg1fTPmuaEDQTK5ujzGw2embeERGX0znsfQTQrWoeBnwuIo6lU2U8vKfZYcC2sv6wAevpabM1IvYC9ge2l/VP7WtzOXAbcEBE7FWqm737moum/eFLkubDoKkm2qPOnUfEw0pFk4jYD/gZ4POZ+SOZ+YjMfASdUPgTmfmvwEXAKRGxT5kxfhRwVWZ+E7grIo4v4zFfyH1jPS8CTi23nw18sowV3QScFBEHRsSBdCqpm8pjl5VtKW17x43OVNP+8CVJ89G2oNk0uzJmtjRN3ZXNQ4DzyrjNPYALMvNjwzbOzBsj4gLgJmAHcHo5DA/wcu479dElZQF4H/DBMploO53Z7GTm9oh4E3B12e6N3clCwGuBjRHxZuDzZR8z5x++JGkWmvh5Y1VzedQ9G/064AljtnlE3/23AG8ZsN01wNED1v8AeM6QfZ8LnDtg/S3AsaP6Vbcm/uFLkjQLTQuajtkcrdbD6JIkab6aVtxoWtDUeF6ucg6a9ocvSWqmJn7edPvcpNCZBDut3w3lMyNJ0oJbSfBqYtDsalLQ1HiGTUnSUlm2oGLQbAZnow9n2JyD9WvXNfaPSZIW2bK9t7YtaGo5OWZTkrQ0liloTRo0l+FnX7Z/FtRh2JyDZXhDkCTVp0mha/3adVP5XGvSz9zPUx+NZticAcOlJKlO8/qcmVbQ1HIzbNZo2f4AR/3X2XuqikG3JUnVNOnwue/xXcHOdBrMMIbNmjT1D3ClhzF62w273dXU50aSFs2yvJ82+RC6xjNs1qApf/zz+uN2dqUkLZ55HY1ahqCZwC5P8DOUYXPKFiUUdd80luGPGO57M1qU51eS6jTova73/byO90KDpupi2JyieQ7QnmR9kxk6JS27Ye9vy/a+t2yfUc5GH86wOSWzeBNYtj/M1eh/LpbtTViSpGVh2FyBuoONoXJydR9ekiRNz7J9zmU6G30Uw+aEJg0yy/YH1QQeapckaXEYNmtiyJw/z/MpSYtnWT8fdzlmcyjD5oQGVc2W9Q9nGRg4JWlx+HnZTobNih71uO/B9ffd9w+mORzPKanp+v9x9h/pxdK5NrpjNofxmanoK9c9cN5dkCS1VO9lgHvvN4lFmvaysqlWafIbtSRpUTkbfRSfGbXS+rXr/C9bkmbE99t2s7IpSVJDeFRmMXlt9NF8ZtRqVjglqV6+x8qwKeGboSRJdfEwuiRJqkWb/pHfmZ7UfRgrm1LRpjdFSZJmxbAp9TBwStJ0tOn9NAl2ssfMlqZpXo+lmrXpDVKSpLoZNqUBDJyStHJtfA/dlXvMbGma5vVYmgHPZSdJ0nQYNqUB2vhfuSRNQxvfPxMcszlC83osSZIWUhuDpsYzbEpD+Ka5PPxdSqpTEuzM2S1N40ndpRG6IcUxnM2wkt/XuCDq716qxn/qNIxhU9LS6f/Q6w+MfihKmrZdHiweyrApVdAbTqx0LaZRAXI14XL92nX+zqUx/AdOoxg2JTWaH3KS5i0Tdjbw/JezYtiU1DizDphWN6Xh/IdP4xg2JTXGPD/UDJyShgt20bxZ4rNi2JS08KycSIvJv01VYdiUtLAW7YPMU2FJ0uQczSpNaNECkGbP14Dk30GvpDNBaFZL0zSvx5KWih9YkrTcPIwuae7GnYR9ETWxz9K0+E/i/e20fjeUYVNaAWcmT0/v89j9APODTJKWh2FTWiErW9PXfQ6bGDb9B0Rt0cS/z7olwa701EfDWPOVJElSbQyb0pT43/50NPl5XL92XaP7L43j63u4newxs6VpPIwuaap6z0U56VADP8gkafkYNiVNRX9QHBQcR50UfZmCpuM3pXZJYFcDz385K4ZNSSu20oC4TMFSkjSaYVNapUGHi5dd235eSf7djxbsxNnow1jzlVao/zQ9bThs6gQYSdKkrGxKK9CmamZbfs5pasM/HmoX3wdGc8zmaIZNaQKjTjq+bJNC/HCRBL4XaPUMm9IERr3pNj1ojvtAaVM1d7WW7R8PSeM5ZnM4w6ZU0bKGrSohU1I7LeN7nmbPsClVtGxVzUlCph84k2ni60GS6mLYlFqmSnAcNTZVUjv4919dZjhBaATDprRK86pi9X4QTPsykH7IrIwVTS0L3wM0TYZNaRXmHTS91rgkLYadVjaHMmxKI4yaFDSPoFmlmmnAlCQtEsOmNEITgmbVcLmss+klTZfvE5NLYJenPhrKmq80xCKMv+teHnJQ0Kx66Ugn+0hSu0TEhoj4ckRsiYjXDXh8/4j4m4jYHBE3RsSL6+yPlU1piHlWNUcFw0kDpiFzdjyZu5rO94uVioUZsxkRewLvAX4W2ApcHREXZeZNPZudDtyUmf85Ih4GfDkiPpSZP6yjT4ZNaYBhh5xXGiSqnqNzGm/0VjLnq+rkLUmqybHAlsy8BSAiNgInA71hM4GHREQADwa2Azvq6lCtYTMi9gU+BexTvteFmXlGRLyJzg++C/gW8KLM3FbavB54KbATeGVmbirrnwh8ANgPuBh4VWZmROwDnA88Efg28NzM/FppcyrwO6U7b87M88r6I4CNwEHA54AX1JXmtTwmCQ+TBL1pBczufgyZkjRbCezKmY7ZfGhEXNNz/5zMPKfcPhS4teexrcBxfe3fDVwEbAMeQic77aqrs3VXNu8GTsjM70bE3sCnI+IS4K2Z+bsAEfFK4PeAl0XEY4BTgMcCa4FPRMSjMnMncBZwGvBZOmFzA3AJnWB6e2YeGRGnAGcCz42Ig4AzgGPovA6uLWXk28s2b8/MjRFxdtnHWTU/F2qI3lBZNWBWqYLWEQKtYkpSK92WmccMeWxQ6s2+++uBLwAnAI8ELo2If8zM70yxj/eqNWxmZgLfLXf3Lkv2/TAP4r4n4WRgY2beDXw1IrYAx0bE14A1mXkFQEScDzyTTtg8Gfj90v5C4N2lLLweuDQzt5c2lwIbSjn5BOB5pc15pb1hU8B94+6qzt7ubltlu9WEwkGH2w2Zi2uSk+5Li8D3k9XZuThzrrcCh/fcP4xOBbPXi4E/LDltS0R8Ffg/gavq6FDtYzbLQNVrgSOB92TmlWX9W4AXAncCTyubH0qnctm1tay7p9zuX99tcytAZu6IiDuBgxlcRj60PHZHZu4YsK/+vp9Gp5rKvjxwkh9bDbeSw+CjAsW0Q6YkTZPvLUvlauCoMmTwX+gcMX5e3zbfAE4E/jEifhR4NHBLXR2qPWyWQ+CPj4gDgI9ExNGZeUNm/jbw22WM5ivoHPIeVvodVRKetE2V8nK37+cA5wCsiYMGbiONM+mb+CwOv0uSpieJWY/ZHKoU3l4BbAL2BM7NzBsj4mXl8bOBNwEfiIjr6eSi12bmbXX1aWaz0TPzjoi4nM5Yyxt6HvoL4G/phM1hpd+t5Xb/enrabI2IvYD96cyq2go8ta/N5cBtwAERsVepbg4qL0sTG3Tam5UGTQPm8vAQuqRZy8yL6cxv6V13ds/tbcBJs+pP3bPRHwbcU4LmfsDPAGdGxFGZeXPZ7BnAl8rti4C/iIi30ZkgdBRwVWbujIi7IuJ44Eo6h9/f1dPmVOAK4NnAJ8ss9U3AH0TEgWW7k4DXl8cuK9tuLG0/WtuToFZZ7eFyQ6YkNdOuxRmzuXDqfmYOAS6LiOvojCG4NDM/BvxhRNxQ1p8EvAogM28ELqBzLqiPA6eXw/AALwfeC2wB/pnO5CCA9wEHl8lErwFeV/a1nU6Z+OqyvLE7WQh4LfCa0ubgsg9pqqpOGuoyaC4fq5pqAt97VLe6Z6NfBzxhwPpnjWjzFuAtA9ZfAxw9YP0PgOcM2de5wLkD1t9C56SnUq36D4v3Vy99k19eBk1J6rDmK9VkUJA0XLaDQVNql0zYmTGzpWm8XKVUI8Nl+xg01SS+R2kWDJuSJEmrtCinPlpEHkaXpCmxqilJ92dlU5KmwKCppvEQ+vR0Tupu/W4YnxlJkiTVxsqmJK2SVU1JOwdeDVtgZVOSVs3DkWoaX7OaJcOmJE2BH95SeyWd2eizWprGsClJU2LglKT7M2xK0hQZOLXofI3WoTMbfVZL0zSvx5IkacWc0KZZcza6JEktYmWzHrucjT6UlU1JmiKrRpK0OyubkjRF/VUjw6cWiVXNemTCzgbOEp8VK5uSVCM/3CW1nWFTkmpkZVOLwn98NC8eRpckSVqlJp6SaFZ8ZiSpJlY1JcnKpiRJ0qokzbyM5KxY2ZQkack5XlPzZGVTkiRplTyp+3BWNiWpBo7XlKQOw6YkTZlBc7l5SFr9EtiVMbOlaQybkiRNoGn/TBiONW+O2ZQkSVolz7M5nM+MJE3Z+rXrrCZJUmHYlKSaGDg1b74GZ2SG4zUdsylJ2o0f9pLazjGbklSjpk0mkTS5xPNsjmJlU5JqYtDUPFlV16KwsilJNTBoSu3SxLGUs2JlU5IkSbUxbEqSJKk2HkaXpCnzELrULt3LVWowK5uSJC0ZJwdpkVjZlKQpGlbVHPfhbzVUajYrm8MZNiVpSgYFxqoVpkHbGUAlLQPDpiTVwMOYUnskzbyM5Kw4ZlOSVqG3+jjNgGlVUyvlPzpaNFY2JWkV+j/Yp/FBb9CUmsfLVQ5nZVOSVmHawdCgKWnZWNmUpFXw0Lkk0tnoo1jZlKQFYNCUtKwMm5K0AJzUITVX9wpCs1qaxrApSQti/dp1hk5JS8cxm5K0QDycLjVTEyuOs2JlU5KkJWFlXIvIyqYkLQirmlIzeQWh0axsStICMGhKWlaGTUmaM4OmpGXmYXRJkqRVSg+jD2VlU5IkSbUxbErSnDmDWGq+XcTMlqYxbEqSJKk2jtmUJElahUxP6j6KlU1JkpaAwzG0qKxsStKceeojqfmcjT6clU1JkiTVxsqmJM2RVU1pGXi5ylEMm5I0YwZMSW1i2JSkGTJoSsvJMZvDOWZTkmbEoCmpjaxsStIMGDSl5ZV4ns1RrGxKUs0MmpLazMqmJEkN5wnd5yw7VxHSYLVWNiNi34i4KiI2R8SNEfGGsv6tEfGliLguIj4SEQf0tHl9RGyJiC9HxPqe9U+MiOvLY++MiCjr94mIvyzrr4yIR/S0OTUibi7LqT3rjyjb3lzaPqDO50FS+2zattmKpiRR/2H0u4ETMnMd8HhgQ0QcD1wKHJ2ZjwO+ArweICIeA5wCPBbYAPxZROxZ9nUWcBpwVFk2lPUvBW7PzCOBtwNnln0dBJwBHAccC5wREQeWNmcCb8/Mo4Dbyz4kqRZWnSS1Wa1hMzu+W+7uXZbMzL/LzB1l/WeBw8rtk4GNmXl3Zn4V2AIcGxGHAGsy84rMTOB84Jk9bc4rty8ETixVz/XApZm5PTNvpxNwN5THTijbUtp29yVJU9MNmVY4peW3i5jZ0jS1TxCKiD0j4gvAt+iEvyv7NnkJcEm5fShwa89jW8u6Q8vt/vW7tSkB9k7g4BH7Ohi4oyfs9u6rv++nRcQ1EXHNPdxd7QeWJKxmSlJX7ROEMnMn8PgyLvMjEXF0Zt4AEBG/DewAPlQ2HxTXc8T6lbQZta/+vp8DnAOwJg5y6K+kiVjRlNoh8aTuo8zs1EeZeQdwOWWsZZmw8/PA88uhcehUGQ/vaXYYsK2sP2zA+t3aRMRewP7A9hH7ug04oGzbvy9JkiRNUd2z0R/WnWkeEfsBPwN8KSI2AK8FnpGZ3+tpchFwSplhfgSdiUBXZeY3gbsi4vgy5vKFwEd72nRnmj8b+GQJr5uAkyLiwDIx6CRgU3nssrItpW13X5K0as5El9om2JWzW5qm7sPohwDnlRnlewAXZObHImILsA9waTmD0Wcz82WZeWNEXADcROfw+unlMDzAy4EPAPvRGePZHef5PuCDZZ/b6cxmJzO3R8SbgKvLdm/MzO3l9muBjRHxZuDzZR+StCoGTEm6v1rDZmZeBzxhwPojR7R5C/CWAeuvAY4esP4HwHOG7Otc4NwB62+hczokSZoaZ59L7eVJ3YfzcpWSNGXORJek+3i5SkmSpFVyNvpwVjYlSZJUGyubkjQljtWU2inTyuYoVjYlSZJUGyubkjQFVjWldmvi+S9nxcqmJEmSamNlU5JWyaqmJM+zOZyVTUmSJNXGsClJkqTaeBhdklbBQ+iSwFMfjTIybEbEQRX2sSsz75hSfySpMQyakjTeuMrmtrKMiut7Ag+fWo8kqQEMmloU69eum3cXWi8JK5sjjAubX8zMJ4zaICI+P8X+SJKkigyaaoJxYfMnK+yjyjaStBSsaEoaxDMfDTcybGbmD7q3I+JA4PDeNpn5ud5tJGmZGTQlaXKVZqNHxJuAFwH/zH3hPYET6umWJC0Wg+Z8rV+7zt9BHw+hL5B0NvooVU999EvAIzPzh3V2RpKkQQyaUnNVDZs3AAcA36qxL5K0kAw6ksZy0OZQVcPm/wQ+HxE3AHd3V2bmM2rplSQtAEOmJK1e1bB5HnAmcD2wq77uSJK0u+7YRMO/FpljNoerGjZvy8x31toTSVogBpvF4e/iPk4KUhNVDZvXRsT/BC5i98Pon6ulV5I0R4YbLRpD5uLLBRqzGREbgHfQucrjezPzDwds81TgT4G96RQVf7qu/lQNm92rCB3fs85TH0laKobMxeMpjwyamkxE7Am8B/hZYCtwdURclJk39WxzAPBnwIbM/EZE/EidfaoUNjPzaXV2QpKkQdocNA2ZzZEs1JjNY4EtmXkLQERsBE4GburZ5nnAhzPzGwCZWevZhvYY9WBE/Py4HVTZRpIW2aZtm1sdarR4DJoa46ERcU3PclrPY4cCt/bc31rW9XoUcGBEXB4R10bEC+vs7LjK5lsj4l+AUXH9D4CPTa9LkjQ7hkwtGoOmKrgtM48Z8tigzNY/onQv4InAicB+wBUR8dnM/MoU+7jbNxvl34C3jdnm5in1RZKkVjNoNlQCi3MYfStweM/9w4BtA7a5LTP/Hfj3iPgUsA6YfdjMzKfW8U0lad6saGrRGDQ1JVcDR0XEEcC/AKfQGaPZ66PAuyNiL+ABwHHA2+vqUNXZ6JIkqSYGzeZblFMfZeaOiHgFsInOqY/OzcwbI+Jl5fGzM/OLEfFx4Do6F+t5b2beUFefDJuSWsWKphaNQVPTlpkXAxf3rTu77/5bgbfOoj+GTUmtYdCUVJsFqWwuopGnPuqKiAdGxO9GxJ+X+0d5yiNJTWLQ1CKyqqk2qBQ2gffTuUzlT5b7W4E319IjSZoyg6YWkUFzmQSZs1uapmrYfGRm/hFwD0Bmfp/R596UpIVg0NQiMmiqTaqO2fxhROxHGZEQEY+kU+mUpIVl0JQ0M47ZHKpq2Px94OPA4RHxIeDJwIvr6pQkrZZBU4vKqqbaplLYzMy/i4hrgePpHD5/VWbeVmvPJGmFDJqSZipp5FjKWakUNiPi7zPzROBvB6yTpIVgyNSis6qpNhoZNiNiX+CBwEMj4kDumxS0Blhbc98kSVoaBs0l55jNocZVNv8b8Go6wfJa7gub3wHeU2O/JGkiVjUlaTGNDJuZ+Q7gHRHxa5n5rhn1SZLupxsme6tDBkw1hVXNNnDM5jBVJwi9KyKOBh4D7Nuz/vy6OiZJ/QyaktQ8VScInQE8lU7YvBh4OvBpwLApqXabtm02aKqxrGq2hGM2h6p6BaFnAycC/5qZLwbWAfvU1itJGsKgqSYxaErVw+b3M3MXsCMi1gDfAn68vm5J0n26H9gGTTWNr1mpeti8JiIOAP6czqz0zwFX1dYrSRrAKpGaxtdsi+QMl4apOkHoV8vNsyPi48CazLyuvm5J0v1ZJVKTGDSljqqVTSLi0Ih4EvBw4ICI+Kn6uiWp7QyWkhojgYzZLQ1TdTb6mcBzgZuAnWV1Ap+qqV+StBvDpyQ1U6WwCTwTeHRm3l1nZySpy1Mdqck8hN4+2cCxlLNS9TD6LcDedXZEksBgqeYxWEqjVa1sfg/4QkT8PXBvdTMzX1lLrySpMHxq0fkaFdDIWeKzUjVsXlQWSaqNVwrSMrDSKe2u6qmPzqu7I5K0fu06A6YazaDZYg2cJT4rI8dsRsQF5ev1EXFd/zKbLkpqi/7KZpcf4JI0fxFx0IB1R4xrN66y+ary9edX0ilJqqJbzRxW2exfN0kF1GqppFmIdozZ/JuIeHpmfgcgIh4DXAAcParRyLCZmd8sX78+rV5K0iCThELDoyTNxR/QCZz/CXg0cD7w/HGNRobNiLiLEfOrMnPNhJ2UpN1MGhy7oXRcOO0eep+kArqS/khSU69ZPqnM/NuI2Bv4O+AhwDMz8+Zx7cZVNh8CEBFvBP4V+CAQdFLsQ1bbaUntNu7web9uIKwyhrMbSKvu25ApSYNFxLvYPU6voXMO9l+LiLGnwqx66qP1mXlcz/2zIuJK4I8m6q0k9ekPg8PC4bCJQ6O2NUBqlpzI1mbNvGb5BK7pu3/tJI2rhs2dEfF8YCOdZPvL3HeNdEmaWDcI9gfCqkGz/7Eq+5EkTW7QKTAj4kDg8Mwce3aiqperfB7wS8C/leU5ZZ0kTWyScZRVq0W9h9gHzV6vElglScNFxOURsaacAmkz8P6IeNu4dmMrmxGxJ/ALmXnyFPopqeUmGZs5rf2P+55WQSWtWgsmCAH7Z+Z3IuJXgPdn5hlVzrs+trKZmTsBg6akVTPUSVKj7RURh9A52v2xqo2qHkb/p4h4d0T83xHxE91lRd2U1EqTBM2VhNJJD51L0+RrTfee/mgWy/y8EdgEbMnMqyPix4Gxpz6qGjafBDy2fJM/Kcsfj2sUEftGxFURsTkiboyIN5T1zyn3d0XEMX1tXh8RWyLiyxGxvmf9E8tlM7dExDsjIsr6fSLiL8v6KyPiET1tTo2Im8tyas/6I8q2N5e2D6j4PEhagXHhcVAw3LRt873LpJyNLknTl5l/lZmPy8xfLfdvyQ72yBIAACAASURBVMxnjWtXaTZ6Zj5thf26GzghM79bTgL66Yi4BLgB+EXg/+nduFz26BQ6wXYt8ImIeFQ5lH8WcBrwWeBiYANwCfBS4PbMPDIiTgHOBJ5bBq+eARxD5/+AayPiosy8vWzz9szcGBFnl32ctcKfUdIIwwJfb7ic9JyYk3wPA6ekmWjBmM2I2JdOZnossG93fWa+ZFS7SpXNiNg/It4WEdeU5U8iYv9x7bLju+Xu3mXJzPxiZn55QJOTgY2ZeXdmfhXYAhxbxgesycwrMjPpXB7pmT1tulPyLwROLFXP9cClmbm9BMxLgQ3lsRPKtpS23X1JmpFBpyoaFgxHra97wpEk6V4fBH6MTsb6B+Aw4K5xjaqeZ/NcOtXIXyr3XwC8n051cqQym/1a4EjgPZl55YjND6VTuezaWtbdU273r++2uRUgM3dExJ3Awb3r+9ocDNyRmTsG7Ku/76fRqaayLw8c+XNK2t2kV+2ZRhVy0KF4q5uSapcs+0ndu47MzOdExMmZeV5E/AWdMZwjVQ2bj+w7Jv+GiPhClYblEPjjI+IA4CMRcXRm3jBk80G/qRyxfiVtRu1r95WZ5wDnAKyJg1pQIJfmp44JRFY0JWmq7ilf74iIo+lcyvwR4xpVDZvfj4inZOanASLiycD3J+ldZt4REZfTGWs5LGxuBQ7vuX8YsK2sP2zA+t42WyNiL2B/YHtZ/9S+NpcDtwEHRMRepbrZuy9JUzDtSmL/+M559EGSRol2lKTOKVcO+l3gIuDBwO+Na1R1NvrLgfdExNci4uvAu4H/Nq5RRDysVDSJiP2AnwG+NKLJRcApZYb5EcBRwFWZ+U3grog4voy5fCHw0Z423ZnmzwY+WcZ1bgJOiogDyxNzErCpPHZZ2ZbStrsvSQtuXLWyGzI99ZFmxdeZ2iIz35uZt2fmP2Tmj2fmj2Tm2ePaVZ2N/gVgXUSsKfe/U7FfhwDnlXGbewAXZObHIuIXgHcBDwP+NiK+kJnrM/PGiLgAuAnYAZxeDsNDJ/B+ANiPziz0S8r69wEfjIgtdCqap5Q+bo+INwFXl+3emJnby+3XAhsj4s3A58s+JDXEsEtSen10SXPTgspmKSC+kM6h83szZGa+cmS7TqFv7M73p3MaoZ8qq/6BTni7c4X9bZw1cVAeFyfOuxtSY0wr6K20amTQ1KxY2Vwcn8gLr83MY8ZvOV37PPzwXPtbr57Z9/vaK39zLj9nRHyGzkTu64Fd3fWZed7QRsxgNrokrYQhU01g0FTL7JuZr5m0UdUxm4/MzDPKmeJvycw3AD8+6TeT1B5VPoSHbWPQlKSF9MGI+K8RcUhEHNRdxjWa2Wx0Se0z7hyXg8ZdroQhU9K8tWQ2+g+BtwK/zX2jVJMxBciqYfPldCb6dK8adDvwosn7KKktRgXAQacyMmhK0sJ7DZ0Tu982SaO6Z6NL0v2stqJpwJSkubgR+N6kjSqFzYj4A+CPMvOOcv9A4Dcy83cm/YaStFKGTEkLqx2Xq9wJfCEiLgPu7q4cd+qjqhOEnt4NmmWntwM/t5JeSmqHqtVKZ/OqqXztqoX+F/AW4DPAtT3LSFXHbO4ZEftk5t1w79WA9llhRyUJqP5hbUVT0kJLWnFS93Hn04yIv87MZ/Wvrxo2/1/g7yPi/XSezpcAI7+hpPaqEg6rBE1DpiQ1ysBZ6VUnCP1RRFxH59rmAbwpMzdNsXOSlkz/jHMnAWmZeAhd99OCymYFA5+FqpVNMvPjwMcHPRYRV2TmT66wY5KWTP8H8SQfzIZMSVouVScIjbPvlPYjqaG6IbE3LE4aHA2akpoqcnbLAhs4Jb9yZXOMxf7RJdWuv3o5yaFzQ6aaxEPo0lCvHbRyWmFTknYLmAZNSa2yxGW3iLieET9hZj6ufP27QY9PK2y24kymkgbrP3TuTPPRuteMH3fteC0eq5pqqZ8vX08vXz9Yvj6fClcUqnoFoVcAHyoncx/kBVX2I2n5+WE8XPe5GTS+VVLDLXFlMzO/DhART87MJ/c89LqI+CfgjaPaV61s/hhwdUR8DjgX2JSZ9z6tmXnDZN2WtCw8rVF1bf25JS2NB0XEUzLz0wAR8STgQeMaVT3P5u9ExO8CJwEvBt4dERcA78vMf15FpyU1nKc1ktR2DZglPi0vBc6NiP3L/TvoXOhnpEnOs5kR8a/AvwI7gAOBCyPi0sz87yvosKQWMGBqWThERG2XmdcC6yJiDRCZeWeVdlXHbL4SOBW4DXgv8FuZeU9E7AHcDBg2pZZYyWFzSVp6ufxzpSNiH+BZwCOAvSI6P3NmTmXM5kOBX+wOEO3KzF0R8fMRceCIyUOSlsgintLIWd2qm/9gSQB8FLgTuBa4u2qjqmM2f2/EY18sE4d+ouo3lbS85hH6DJqSNBOHZeaGSRtN63KVy187ljQ21Bn6tIysaqqSnOEyP5+JiP8waSMvVympslEfugZNSVp6TwFeFBFfpXMYPejMIX/cqEZerlLSRLzyjSTdX0tOffT0lTTyMLqkiSxS0PTwpurma0y6T2Z+vUwW/z4THNivXNmMiD2BH+1tk5nfKDdPnKi3khppUUJm16L1R1KLtaCyGRHPAP4EWAt8C/g/gC8Cjx3Vrup5Nn8NOAP4N2BXWZ3A4wAyc/uKei1JkqSmeBNwPPCJzHxCRDwN+OVxjapWNl8FPDozv72KDkpqsCpVxO4hRyuOklqlPZervCczvx0Re0TEHpl5WUScOa5R1TGbt9I5iackDWXI1DJxvKZ0P3dExIOBTwEfioh30LmE+UhVK5u3AJdHxN/Sc8b4zHzbSnoqaTn54SyptdpR2TwZ+AHw68Dzgf2BkZeqhOph8xtleUBZJOlehkxJWn6Z+e89d8+r2q7q5SrfABARD+ncze9O1j1Jy8zD51o2/gOliS1xZTMi7mLwT9g9qfuaUe0rjdmMiKMj4vPADcCNEXFtRIyc5i5puSzSh+8i9UWSll1mPiQz1wxYHjIuaEL1w+jnAK/JzMsAIuKpwJ8DT1pxzyVphaykqk7+M6OVaMls9BWpOhv9Qd2gCZCZlwMPqqVHkiRJWhqVZ6NHxO8CHyz3/wvw1Xq6JGkRWU2UJK1E1crmS4CHAR8GPlJuv7iuTkmSJGk5VJ2Nfjvwypr7ImmBrV+7zuqmJGliI8NmRPxpZr46Iv6GAVPeM/MZtfVMkqQZc3KQVswJQkONq2x2x2j+cd0dkbT4rG5qmRk0pXqMHLOZmdeWm4/PzH/oXYDH1989SYvEoKll5utbK5adUx/NammaqhOETh2w7kVT7IckSZKW0Lgxm78MPA84IiIu6nnoIcC36+yYJEmz5GF0rUoDK46zMm7M5meAbwIPBf6kZ/1dwHV1dUrSYvHwoiRppUaGzcz8OvD1iHg+sC0zfwAQEfsBhwFfq72HkuauW/HpD51OGJKkwsrmUFXHbF4A7Oq5vxP4q+l3R1KTGDS1LDyELtWnatjcKzN/2L1Tbj+gni5JWkQGSy0zX99ajcDZ6KNUDZv/f0TcewL3iDgZuK2eLkmSNFtWNqX6VLpcJfAy4EMR8R46oxK2Ai+srVeSJM2IQVNT0cCK46xUvTb6PwPHR8SDgcjMu+rtliRJs7Fp22YDp1SjSofRI+JHI+J9wF9l5l0R8ZiIeGnNfZMkSVp8XkFopKpjNj8AbALWlvtfAV5dR4ckLR4nT2iZWdWU6lU1bD40M+89/VFm7qBz+iNJkiTlDJeGqRo2/z0iDqb8iBFxPHBnbb2SJGkGrGpK9asaNl8DXAQ8MiL+CTgf+LXaeiVJkqQViYgNEfHliNgSEa8bsd1/jIidEfHsOvtTdTb65yLip4FH0zl36Zcz8546OyZpMTheU5IqWJDD2xGxJ/Ae4GfpnKry6oi4KDNvGrDdmXTm5NRqZNiMiBMy85MR8Yt9Dz0qIhLYDnw6Mx2/KUmSNH/HAlsy8xaAiNgInAzc1LfdrwF/DfzHujs0rrL508Angf885PGDgd+hk54lSWoMx2tqmmZ8SqKHRsQ1PffPycxzyu1DgVt7HtsKHNfbOCIOBX4BOIF5h83MPKN8ffGwbcr5NyUtoe4h9PVr13k4XZIWx22ZecyQx2LAuv4o/KfAazNzZ8Sgzaer0pjNiNgfOAP4qbLqH4A3ZuadmenJ3aUlZvVHkipYkDGbdCqZh/fcPwzY1rfNMcDGEjQfCvxcROzIzP9VR4eqzkY/F7gL+KWyfAd4fx0dkrRYNm3bbFVTc1PXPzv+E6UldjVwVEQcEREPAE6hc0ahe2XmEZn5iMx8BHAh8Kt1BU2oWNkEHpmZz+q5/4aI+EIdHZIkqTt0w3901AgLdLL1zNwREa+gM8t8T+DczLwxIl5WHj971n2qGja/HxFPycxPA0TEk4Hv19ctSVJbdauOjhWWViYzLwYu7ls3MGRm5ovq7k/VsPky4PwydhPgduDUerokaRH4Ia958bWnJprxbPRGGRs2I2IP4NGZuS4i1gBk5ndq75kkSZIab+wEoczcBbyi3P6OQVOS1GRODlItcoZLw1SdjX5pRPxmRBweEQd1l3GNImLfiLgqIjZHxI0R8Yay/qCIuDQibi5fD+xp8/pyLc8vR8T6nvVPjIjry2PvjDJfPyL2iYi/LOuvjIhH9LQ5tXyPmyPi1J71R5Rtby5tH1DxeZAkSdIEqobNlwCnA58Cri3LNSNbdNwNnJCZ64DHAxsi4njgdcDfZ+ZRwN+X+0TEY+hM0X8ssAH4s3LtToCzgNOAo8qyoax/KXB7Zh4JvJ3OdT4pYfgMOmfNPxY4oyfUngm8vXz/28s+JElLzqqm6hI5u6VpKoXNcj6m/uXHK7TLzPxuubt3WZLONTrPK+vPA55Zbp8MbMzMuzPzq8AW4NiIOARYk5lXZGYC5/e16e7rQuDEUvVcD1yamdsz83bgUjphN+hcnunCAd9fkiRJU1T1CkL7Ar8KPIVOWPxH4OzM/EGFtnvSqYQeCbwnM6+MiB/NzG8CZOY3I+JHyuaHAp/tab61rLun3O5f321za9nXjoi4k8412wddG/TQ8tgdmbljwL76+34anWoq+/LAcT+qJElqqwZWHGel6mH08+kc2n4X8G7gMcAHqzTMzJ2Z+Xg6l0s6NiKOHrH5sOt5jrrO56RtqlwztLMy85zMPCYzj9mbfQZtIklqCA+hS/NR9Tybjy7jLrsui4iJToSWmXdExOV0xlr+W0QcUqqahwDfKpsNu57n1nK7f31vm60RsRewP7C9rH9qX5vLgduAAyJir1LdHHTNUElTNOxD3vMpStLyq1rZ/HyZ2ANARBwH/NO4RhHxsIg4oNzeD/gZ4Et0rtHZnR1+KvDRcvsi4JQyw/wIOhOBriqH3O+KiOPLmMsX9rXp7uvZwCfLuM5NwEkRcWCZGHQSsKk8dlnZtv/7S2J1IbA/WA4Kml6GULNmVVO1muVpjxp4uL5qZfM44IUR8Y1y/+HAFyPiejrzgB43pN0hwHll3OYewAWZ+bGIuAK4ICJeCnwDeA6dHd0YERcANwE7gNMzc2fZ18uBDwD7AZeUBeB9wAcjYgudiuYpZV/bI+JNdC5ID/DGzNxebr8W2BgRbwY+X/YhaUqGBUxJUvtUDZsbxm9yf5l5HfCEAeu/DZw4pM1bgLcMWH8NcL/xnmWS0nOG7Otc4NwB62+hczokSTNg0JS0zILBE0LUUSlsZubX6+6IpOazoqlF5CF0ab6qjtmU1CLT/HD2g15SKzhmcyjDpqSpGBUqDZyS1F5Vx2xKapFpH/r2ULrmxX90NCtNvIzkrFjZlCRJUm2sbEpatVHVI6uamhermpopK5tDWdmUVBuDpiTJsCnpfiapCFk90iLydamZczb6UIZNSbWwqqlehj+pvRyzKWnFhgUIg6Z6rV+7zteElls6G30Uw6ak3XRDQTdIDgoJTgjSJOp6TQx6jVpBlRaPYVPSUP3BU1oUg6qlvk41V1Y2hzJsSrrXpm2bd/vAnvTD26qmZsXXmtQcThCSBKz+w9sPf82TVU1pcVnZlHSvlX5gGzQltZ0ThIazsinpfofPJ2ln0NQsWcGUmsfKpqShnOWrlRp1NoO6vpc0V1Y2h7KyKbXcsDDQv37cfamrrqA5KFQaNKXFZ2VTarGVntrIoKlR6gyZVtu1qByzOZxhU9JEDJqaNV9zUrMZNqWWszqkJvJ1q4WSOGZzBMdsSi3m4XM1kUFTahbDpiRpYfQHyXH3pYWRM1waxrAptVjVmejSrHjWA2n5GDalFrNKpEUz6jXp61WLKujMRp/V0jSGTamFvPKPFpWvS2n5OBtd0m78sNcisqqphdfAiuOsWNmUWsoPby267mvU16rUbFY2Jd3LqqYWyaZtmw2aaoxIS5vDWNmUWsgPcEnSrBg2pZaZpHppKFXd1q9dd+8y6DFJzWfYlATcP4T2ftD7oa+6DPvnx9ecGmWWJ3Rv4NF6w6bUIisZk+mHvmbN15y0XJwgJLXMoA/yUVXNQY9L09T7+jJoqqmaeLL1WbGyKbVE1UtTGjQ1LwZNaTlZ2ZRapr+KNO4Dvvu4oXMxrV+7zt+NtAisbA5lZVNqsW5I6V6+0tDSPMvyO7OqKS0vw6bUEsM+zPvDyrKEFzWHQVPLIHJ2S9MYNqUWWc2HuoFA09I7fMPXlbT8HLMp6V5Vxm9a+dRq+RrSUmpgxXFWDJtSy6y2kmTg1DRY0ZTaw7ApaSIGTUnq09CxlLPimE1JE7EipdXyNSS1i5VNSdLMGDS1tKxsDmVlU5JUG8OlJCubkqTaeN1ztUHgmM1RrGxKkmpn0JTay7ApSaqVQVNqNw+jS5qIpz6SpAHS4+jDWNmUJNXGqqYkK5uSpFoYNNUmThAazsqmpMo8hK6qDJqSugybkqQV6w2VBky1Vs54aRgPo0uqxKqm+nXDZfdr9zVi6JTUy7ApSZrYsIBp0FRbxa5592BxeRhd0lhWNdWv/zWxadtmg6akgQybkqSJ9Fc1e9dJreWYzaEMm5JGMlCoa/3adbsdLvfQuaQqDJuSKvNwenutX7tu4KFzSR2Rs1uaxrApaSjDhAYdMt+0bbMzzyVVZtiUJA3VP/Gn975BUyqSzrXRZ7U0jGFT0kBWNdXlzHNJq2HYlDSWwUKSRnPM5nCe1F3S/TgRRMP4j4ekSRk2Je3GYKlhDJrSCA2sOM6Kh9ElSZJUm1rDZkQcHhGXRcQXI+LGiHhVWb8uIq6IiOsj4m8iYk1Pm9dHxJaI+HJErO9Z/8Sy/ZaIeGdERFm/T0T8ZVl/ZUQ8oqfNqRFxc1lO7Vl/RNn25tL2AXU+D1JTWNWUJE1b3ZXNHcBvZOb/BRwPnB4RjwHeC7wuM/8D8BHgtwDKY6cAjwU2AH8WEXuWfZ0FnAYcVZYNZf1Lgdsz80jg7cCZZV8HAWcAxwHHAmdExIGlzZnA2zPzKOD2sg9J0hAeQpeGC5wgNEqtYTMzv5mZnyu37wK+CBwKPBr4VNnsUuBZ5fbJwMbMvDszvwpsAY6NiEOANZl5RWYmcD7wzJ4255XbFwInlqrneuDSzNyembeX77OhPHZC2ZbStrsvqbWsaqpXb7g0aEpajZmN2SyHt58AXAncADyjPPQc4PBy+1Dg1p5mW8u6Q8vt/vW7tcnMHcCdwMEj9nUwcEfZtn9f/X0+LSKuiYhr7uHu6j+sJDXIoOude4UgaQKzPKG7J3UfLCIeDPw18OrM/A7wEjqH1K8FHgL8sLvpgOY5Yv1K2oza1+4rM8/JzGMy85i92WfQJpLUeL0nae+tcBs0JU1D7ac+ioi96QTND2XmhwEy80vASeXxRwH/qWy+lfuqnACHAdvK+sMGrO9tszUi9gL2B7aX9U/ta3M5cBtwQETsVaqbvfuSpNYZFDQlTaaJYylnpe7Z6AG8D/hiZr6tZ/2PlK97AL8DnF0eugg4pcwwP4LORKCrMvObwF0RcXzZ5wuBj/a06c40fzbwyTKucxNwUkQcWCYGnQRsKo9dVraltO3uS2olQ0a7bdq2+X6vAauakqal7srmk4EXANdHxBfKuv8BHBURp5f7HwbeD5CZN0bEBcBNdGayn56ZO8t2Lwc+AOwHXFIW6ITZD0bEFjoVzVPKvrZHxJuAq8t2b8zM7eX2a4GNEfFm4PNlH5IkDJrSiljZHKrWsJmZn2bwGEmAdwxp8xbgLQPWXwMcPWD9D+hMMhq0r3OBcwesv4XO6ZAkqXXWr103tJpt0JQ0bV6uUpIkaZUcszmcl6uUpJaxqilplqxsSpIMmtJqJLDL0uYwVjalFhs0C1n1WsRQt4h9krQ8rGxK0gwZ7qUlZWFzKCubktRiVjUl1c2wKbWYQaPd/P1L0xM5u6VpDJtSi3lIt1kMh5KayLAptZRBs3mm+TszuEqaFcOm1FKGjcVXx+9o/dp1/u6lOmTObmkYZ6NLAoYHGyug89X9vXR/D4MuNTnq8pOSNG+GTUlWuhZQ7++kN0gOCpUGTWn+mjhxZ1YMm1JLGVAW27R/P/5DIWleHLMptViV8Xvdxw0r89P9PQ36Xfh7kRZAzngZIyI2RMSXI2JLRLxuwOPPj4jryvKZiKj1jcTKptRSk4QUA81i6FY7N23bfL+xnJIEEBF7Au8BfhbYClwdERdl5k09m30V+OnMvD0ing6cAxxXV58Mm5IqM+DM3qDn3N+DtFgCiMWZJX4ssCUzbwGIiI3AycC9YTMzP9Oz/WeBw+rskGFT0sSc/Vy/3kDZP9yht7JZdT+SlspDI+KanvvnZOY55fahwK09j21ldNXypcAlU+7fbgybkiZm0Jyu/mC5advmgSGzf/thM9b7t5M0A7tm+t1uy8xjhjwWA9YNLLtGxNPohM2nTKtjgxg2JU1sXMhRdf1V4nHn0+xuMypIGjKlVtsKHN5z/zBgW/9GEfE44L3A0zPz23V2yNnoklbFK9KszrCwXmX9oHGckuYjMme2jHE1cFREHBERDwBOAS7ara8RDwc+DLwgM79SyxPSw8qmJDXEoKrnsJO8Gz6ldsrMHRHxCmATsCdwbmbeGBEvK4+fDfwecDDwZxEBsGPEYflVM2xKWjUPpc/WqOfbkCnNQcXzX85KZl4MXNy37uye278C/Mqs+uNhdEmrZsCZvf4TvXfXSdKisbIpSQ0yLFAaNKV5Slic82wuHMOmJDXAoDDp8AVJTeBhdEmrMu6ckFq9cdVMn3dJi8zKpiQ1mEFTWgzhUfShDJuSVqU/7Hgpy+kxSEpaBh5GlzRVBs3pGBc0fZ6lBZM5u6VhDJuSpsYANB1VKppWPSU1hYfRJWmBGCKlBkqIXfPuxOKysilJC8KgKWkZWdmUNBUeQl85Q6a0BBo4lnJWrGxKmgoD08r4vEladoZNSVNjcJqMz5e0RHKGS8MYNiVNjYfSqzNoSmoLx2xKkiStUjhmcygrm5I0Y1Y1JbWJYVPS1BiixvM5kpaUVxAayrApSTNi0JTURo7ZlDQ1ThAazJApLbkEvILQUFY2JU2FQXOwSYOmz6OkZWPYlKSarCRoWgWVtGw8jC5JC8CgKTVXkJ76aAQrm5JUg0mCo4fOJS0zK5uSNGWTViitaEpLwMrmUFY2JUmSVBsrm5IkSatlZXMoK5uSJEmqjZVNSZoix19KLeRJ3Ueysilp1ZxN3bF+7TqfC0nqY2VT0qoYrnbXX9nsfX6sekrLy/NsDmfYlKQp6QbLbqjsvy9JbWTYlLQqHjruGFXRlNQCVjaHMmxK0ioNC5pWNCXJsClpldpewesNlG1/LqT2SiubIxg2JU3MUOVEIEmqylMfSVJF3RBp0JS0m6RT2ZzV0jBWNiVV1uaK5qCgaciUpPGsbEqqpM1BcxCDpiRVY2VTkioYFrYNmpIAL1c5gmFT0lhtrGpWOW+mQVOSxjNsShqpjUGzCoOmpF5ernI4w6akodoaNEcFSUOmJE3GsClpIIPmfbwikKSxrGwO5Wx0SSoMmpI0fVY2JS299WvXjazUDguSm7ZtNmRKGi+BXVY2h6m1shkRh0fEZRHxxYi4MSJeVdY/PiI+GxFfiIhrIuLYnjavj4gtEfHliFjfs/6JEXF9eeydERFl/T4R8Zdl/ZUR8YieNqdGxM1lObVn/RFl25tL2wfU+TxImr1BJ2Ffv3bdyPC4advm3RaDpiStXt2VzR3Ab2Tm5yLiIcC1EXEp8EfAGzLzkoj4uXL/qRHxGOAU4LHAWuATEfGozNwJnAWcBnwWuBjYAFwCvBS4PTOPjIhTgDOB50bEQcAZwDF0/ue4NiIuyszbyzZvz8yNEXF22cdZNT8XkqasGwYHBcP+df23eyudntZI0uo08zKSs1JrZTMzv5mZnyu37wK+CBxKJ/ytKZvtD2wrt08GNmbm3Zn5VWALcGxEHAKsycwrMjOB84Fn9rQ5r9y+EDixVD3XA5dm5vYSMC8FNpTHTijbUtp29yVpxiYNdd3q5LAgCeMPf4+b/GTQlKTpmdmYzXJ4+wnAlcCrgU0R8cd0Au+TymaH0qlcdm0t6+4pt/vXd9vcCpCZOyLiTuDg3vV9bQ4G7sjMHQP21d/n0+hUU9mXB07y40qNVsdM9GEnSZ/0e01axRzUfpJ+SlIlVjaHmsls9Ih4MPDXwKsz8zvAy4Ffz8zDgV8H3tfddEDzHLF+JW1G7Wv3lZnnZOYxmXnM3uwzaBNp6Uw7aA4bJ9k/pnLQGMtB7asGyf6fozsOU5I0W7VXNiNibzpB80OZ+eGy+lTgVeX2XwHvLbe3Aof3ND+MziH2reV2//reNlsjYi86h+W3l/VP7WtzOXAbcEBE7FWqm737klqvdxzkLA2qdPb3ZVjQHDYZqCqrmZJWzcrmULWGzTI+8n3AFzPzbT0PbQN+mk74OwG4uay/CPiLiHgbnQlCRwFXZebOiLgrIo6ncxj+hcC7etqcClwBPBv4ZGZmRGwC/iAiDizbnQS8vjx2Wdl2Y2n70en/wgFCsAAAEU9JREFU9FL7rDS0jTq/ZffxScJv/2F1K5qSND91VzafDLwAuD4ivlDW/Q/gvwLvKJXIH1DGRWbmjRFxAXATnZnsp5eZ6NA59P4BYD86s9AvKevfB3wwIrbQqWieUva1PSLeBFxdtntjZm4vt18LbIyINwOf577D+JL6rKbqNyj0TVJ9rHL6oUHb9M9SHxZWrWhKmgrPszlSpGXfStbEQXlcnDjvbkgLqxvspnF+yt59jDqEPsnVfQYdnpe0XD6RF16bmcfM+vvuv8+P5ZMO/S8z+34f/+qfzOXnXCkvVylpKoZN8Kmqd8xm1cPmwyYf9U8GMmhK0vx4uUpJC6Hq7PRhBh2iN2RKmo2E3DXvTiwsK5uS5mpUBXPcqYx6DQua4y5RKUmql5VNSQurf6JPb/VyXIA0YP7v9u4/WK6yvuP4+9MEggH5rZZfklgpFRCphIAKguG30ym0RQQdCPgX1Y5Ih9ZSqAK1DjK1KIUBMzYkaAVKkDZVIQIC8sMgISRA+NHENJaUtAxNDAQmCeF++8d5NpzsPWfv7t49e+/ufl4zO/fsc55zzrNflpvvfZ5znsfMusrPwJRyz6aZjamyeT2LhsBHmtC93VWJzMysOk42zayjRkr0yvbnh75bfaI9P8VR7VxmZl1Tm/qoW68e42TTzDqqmQnaRzq2lYeC3ItpZja++Z5NM6tcJ3saGz0YZGY2ZnzPZin3bJpZVzXTGzlSHQ+Xm5n1DvdsmllHjXS/ZTMrAZUd38rT6GZmXeWezVJONs2so9pJAps9xgmmmVnvcbJpZuNG2VKVTjLNbHwL92w24Hs2zWzccKJpZtZ/3LNpZh3X6jyZ+WOcXJpZzwlgyGujl3HPppl1XDMJo+fHNDMbDO7ZNLOuK0o03aNpZj3N92yWcs+mmXVF0b2Y7t00M+t/TjbNrCvyPZeeK9PMbHA42TSzyrkH08z6XkT3Xj3G92yaWWXqVwYyM7PB42TTzDrKS0qa2eAJGOq9Hsdu8TC6mXWUE0wzM8tzsmlmbfF9mGZmSUDEUNdevcbD6GbWEieZZmbWCiebZl3QDw/K5D9DO8tRmpn1Nd+zWcrJpllFinoAezlJK5on08zMbCS+Z9Osg2oJZqOh5gUvLfVQtJlZv/E8m6WcbJp1QD6BHG0i2e1E1ImvmZlVycPoZm3odILW7bkp8+33kLiZ2ShFwFDvPSXeLU42zVrUiURzrBI8J5lmZtZtTjbNWtDJHs0qE7/6czvJNDOrWA/eS9ktTjbNRtAocRuNTp6r6Nw1vfwEvJmZ9T4nm2YF6pPAWmLY6eRwtAln0bH1iaUTTTOz6oXv2SzlZNMsGWm6oiqudfLeH2o7GSxKiM3MzMYbT31kA6+dRLIssWsl4etkcuhE08zMxiv3bNrAyt/LWL/dzLGtlFfBCaaZ2XjRm5Otd4t7Nm2g1O67rCWXzSSHVcx/Wf8Aj5mZWb9ysmkDoejhnqL39U9ut5qQttKe2k/3UJqZ9bgAhqJ7rx7jYXTre/mHcTrx5HfZ0Hszw/B+UtzMzAaNeza7yMOl3ZeP+UjxHynxa5QoNjq29sS5E0szsz4WQ9179Rj3bHaRk43ua2UYvGjy9namF/J/ZzMzs7c52bSBVTbsXT9cbmZm1kgA0YP3UnaLk02rTLfX4y562KbZidq9driZmVk1nGxaXyiayqjZIW/3ZJqZ2ahE9OS9lN3iB4SsMt1O3oruuWwkP3+mE00zM7NqONm0jipK8DqRyI20Yk99wlg2f6WfDDczsyrEUHTt1Ws8jG4dVf/QTbtJXdHk6kV16ut7HkszMxt0kk4Bvg1MAL4bEVfV7Vfa/0ngDeC8iFhcVXucbFolRpvktTrFkJNKMzMbU+Pknk1JE4DrgROB1cDjkuZHxLO5aqcCB6TXkcAN6WclPIxuYya/ZGN9uSfANzMza8t0YEVErIyIzcCtwGl1dU4Dbo7MQmBXSXtV1SD3bDbpNda9cm/M+3V6uyfwyli2p4dtjd2EvQCWb/1ZM2Hr1305Noy/e+1z7Nrn2I2O49e+VmO3f1UNaeQ11i24N+bt2cVL7iBpUe79rIiYlbb3AV7M7VvN8F7Lojr7AGs63VBwstm0iHhXbVvSooiYNpbt6VWO3eg4fu1z7Nrn2I2O49e+XoldRJwy1m3IUUFZ/VNFzdTpGA+jm5mZmfWP1cB+uff7Ai+1UadjnGyamZmZ9Y/HgQMkTZW0PXAWML+uznzgXGWOAtZHRCVD6OBh9HbNGrmKlXDsRsfxa59j1z7HbnQcv/Y5di2KiC2S/gxYQDb10eyIWCbpgrT/RuAnZNMerSCb+uj8KtukiN6bHNTMzMzMeoOH0c3MzMysMk42zczMzKwyA5VsSjpM0kJJSyQtkjS9bv97JW2QdHGu7HBJT0taIenatMQTkiZJui2VPyZpSu6YmZKWp9fMXPnUVHd5Onb7VK507hWSnpL04apj0aqy2EmansqWSFoq6Y9yxzh2NIzdiZKeSDF6QtKM3DGOXdIgfntIuj/9P3td3TGOH41/50m6JLX9BUkn58odOyC1t/a7bZWkJal8e0k3pRgtlXRc7hjHLmkQv+0kzU1xek7SJbljHL9+FRED8wJ+Cpyatj8JPFC3/w7gduDiXNkvgY+QzUl1V+74zwM3pu2zgNvS9u7AyvRzt7S9W9r3L8BZaftG4E9zbbkrXeMo4LGxjlWzsQMmAxPT9l7Ay7n3jl3j2P0+sHfaPgT4b3/vWorfjsDRwAXAdXXHOH6NY3cQsBSYBEwFfgVMcOxK4/hN4Ctp+wvATWn73cATwG85dk3H7zPArWl7MrAKmOL49fdroHo2ySYs3Tlt70JuTilJp5N9UZflyvYCdo6IX0T2Lb0ZOD3tPg2Ym7bnAcenv8JOBu6JiLURsQ64Bzgl7ZuR6pKOzZ+ra8tGtakwdhHxRkRsSeU7pHqO3bbKYvdkRNS+g8vIVoSY5NgNUxa/1yPiYWBjvrLjt42y33mnkf2Dvyki/pPsidTpjt1w6XOcCdySig4C7gOIiJeB3wDTHLtiBfELYEdJE4F3AJuBVx2//jZoUx99CVgg6e/JbiH4KICkHYEvky1af3Gu/j5kE5/W1JZzqu17EbZOM7Ae2IPyJaD2AH6TS8wKz1W3r7I5r9pQGDsASUcCs8mWCTsnxcOxe1tp7HL+BHgyIjY5dsM0E788x+9tZbHbB1iYq1dr+5s4dvWOAf43Imrr5y4FTpN0K9mk2Ienn0M4dkXq4zePLOFbQ9azeVFErJU0Dcevb/VdsinpXuC3C3ZdChxP9sW+Q9KZwD8BJwBXANdExIZ0i8jW0xWcJ0bY12r5SNfpmjZjR0Q8Bhws6QPAXEm1IYp6jl1d7NKxBwPfAE6qFRWcp29jB6OLX9HpCsr6Nn5txq6Tn7cvYxcR/5a2z+btXjnI/rD+ALAI+DXwKLCFAYsdtB2/6cBbwN5kQ98PpfMMXPwGSd8lmxFR+o+QpJuBC9Pb24Hvpu0jgTMkXQ3sCgxJ2kh2D+e+uVPkl3OqLfW0Og0H7AKsTeXH1R3zAPAKWXf9xPTXVtG5iq7TNW3GLn/8c5JeJ7v/cDWOHdA4dpL2Be4Ezo2IX6XigYodjP67V2eg4tdm7Mra7tjlpM/5x2S9l7VjtgAX5eo8CiwH1jFAsYP24kd2z+bdEfEm8LKkR4BpwEMMWPwGyaDds/kScGzankH2C4KIOCYipkTEFOBbwNcj4rrIlm56TdJR6R6Qc4HaX2vzgdqTb2cAP0v3mSwATpK0m6TdyHqrFqR996e6pGPz5+raslFtKoydsif+Jqbt/YEDgVWO3TbKYrcr8GPgkoh4pFbZsRumMH5lHL9tlMVuPnCWsnuEpwIHAL907IY5AXg+IrYO70qanG69QtKJwJaIeNaxKzQsfsB/ATNS+3cke0jnecevz8U4eEqpWy+yJ1efILvn5jHg8II6l7Pt0+jTgGfInta8DrauurQDWU/BCrIn6N6XO+ZzqXwFcH6u/H2p7op07KRULuD6dI2ngWljHatmYwecQ/ZwyxJgMXC6Y9d07C4DXk+xq73e7dg1//8t2ZOsa4ENZD0WBzl+Tcfu0tT2F0hP/Tp2w+I3B7igrmxKitlzwL3A/o5dS/HbKX2eZcCzwF84fv3/8nKVZmZmZlaZQRtGNzMzM7MucrJpZmZmZpVxsmlmZmZmlXGyaWZmZmaVcbJpZmZmZpVxsmlmZmZmlXGyaWYjknSlpIarhZQcN0XSZ3Lvp0m6trOtG18kbUg/95Y0bxTn+ZKkyS0e84CkFyT9YQvHvEPSEkmbJe3ZekvNzBrzPJtmBoCkCRHxVofPeRzZIgl/0Mnzjhe55fDyZRsiYqcOnHsV2YTTr7RwzANk8V7UjeuZmTXDPZtmfS71Lj4vaa6kpyTNq/WYSVol6SuSHgY+JekwSQtTvTvTEnBImiPpjLR9uKQHJT0haYGkvVL5+yXdK2mppMWSfge4Cjgm9ZxdJOk4ST9K9XeX9K/pWgslHZrKL5c0O/XSrZT0xQaf7Yh0/A6SdpS0TNIhDer/paSnUxuvSmVln7ms/AFJX5f0IHChsiVbfyHpcUl/Wxf3Z9L2eZJ+KOluScslXZ2rd4OkRantV6SyLwJ7A/dLuj+VnZSus1jS7ZJGTGhTW6+R9HNJz6V4/TC14WsjHW9m1glONs0Gw4HArIg4FHgV+Hxu38aIODoibgVuBr6c6j0NfDV/EknbAf8InBERhwOzgb9Lu/8ZuD4iPgR8FFgD/BXwUEQcFhHX1LXpCuDJdK2/Tteu+T3gZGA68NV03WEi4nGytY6/BlwNfD8inimqK+lU4HTgyNTGWsJX9pkbxWLXiDg2Ir4JfBu4ISKOAP6n6NrJYcCngQ8Cn5a0Xyq/NCKmAYcCx0o6NCKuJVvX/BMR8Yk0vH0ZcEJEfBhYBPx5g2vlbY6IjwM3kq0P/QXgEOA8SXs0eQ4zs7Y52TQbDC9GxCNp+/tka2bX3AYgaReyJOrBVD4X+HjdeQ4kS1TukbSELAHaV9I7gX0i4k6AiNgYEW+M0Kajge+l+j8D9khtAPhxRGxKQ7ovA+9pcJ4rgRPJ1lW+ukG9E4Cbau2KiLVln7mJWNyW2/4YcEva/l6D698XEesjYiPZmtD7p/IzJS0GngQOBg4qOPaoVP5IivvM3PEjmZ9+Pg0si4g1EbEJWAnsV36YmVlnTBzrBphZV9TfnJ1//3oL5xFZwvKRbQqlndtokwrKau3alCt7i8a/q3YHdgK2A3ag/POI4XFoV/01mjnvsM8kaSpwMXBERKyTNIfsM9QTcE9EnN1GW2vXHaprwxD+N8DMusA9m2aD4b2Sagni2cDD9RUiYj2wTtIxqegc4MG6ai8A76qdS9J2kg6OiFeB1ZJOT+WT0n2hrwHvLGnTz4HPpvrHAa+k87RqFvA3ZMP432hQ76fA53L3q+5e9pmbjEXNI8BZafuzLbZ9Z7LEdb2k9wCn5vblY7cQ+Jik96e2T5b0uy1ey8xsTPivWrPB8BwwU9J3gOXADSX1ZgI3poRsJXB+bl9ExOb0oNC1aah5IvAtYBlZQvYdSVcCbwKfAp4CtkhaCswhGyquuRy4SdJTwBvp2i2RdC6wJSJ+IGkC8KikGWlYfhsRcbekw4BFkjYDPyG7V7TsMzeKRd6FwA8kXQjc0Ur7I2KppCfJ4reSLHGtmQXcJWlNum/zPOAWSZPS/suA/2jlemZmY8FTH5n1OUlTgB9FROlT2k2c49+Bf4iI+zvVLus8eeojMxuHPIxuZg1Jmg1MpmDo3cadtcActTGpO9k9r0OVtczMBpZ7Ns1s3EtT9NxXsOv4iPi/urofZPhT4Zsi4siq2mdmZuWcbJqZmZlZZTyMbmZmZmaVcbJpZmZmZpVxsmlmZmZmlXGyaWZmZmaV+X+78zYOYif9DQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 792x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "ds.coords[\"land_mask\"] = xr.DataArray(\n",
    "    data=mask_array,\n",
    "    coords=[ds.y, ds.x],\n",
    "    dims=[\"y\", \"x\"],\n",
    "    attrs=dict(\n",
    "        grid_mapping=\"crs\", \n",
    "        flag_values=\"0 1\", \n",
    "        flag_meanings=\"nodata data\"))\n",
    "\n",
    "ds.coords[\"land_mask\"].plot(x=\"x\", y=\"y\", figsize=(11,10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solar Zenith Angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<cartopy.mpl.feature_artist.FeatureArtist at 0x7fae3871b0f0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGaCAYAAABJ15eSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd5wkVbnw8d8zMz3Tk3fSzuacM7AEBQQkrhJFlCiCXhFBReTeC+YXwzXe1+t9laBiABFEWcnIisCCihI3787mNLMTd3Lq8Lx/dM/aO8zO9O52d3V1Pd/Ppz8z3V1VfapOdZ0+zzl1jqgqxhhjTKplOZ0AY4wx3mQFkDHGGEdYAWSMMcYRVgAZY4xxhBVAxhhjHGEFkDHGGEdYAXSURORrIvKA0+k4HCIySUQ6RSTb6bQYY7zLCqDDICKni8gehz77XhHZJCJhEfnoYa67Q0TOGniuqrtUtUhVQwlPaAKISJ6I/FxEdopIh4i8JSLLBi1zpohsFJFuEXlBRCbHvHdG9LU2EdkxxPaXiMjL0ff3iMhXRkjPlOj2uqOfeVbMe2NF5HERqRURFZEpcezfldF96xKRP4pIecx7HxKRv0U/68U4tjXccRAR+Y6INEcf3xUROZL9HCndQ2wrT0TuE5F2EdknIrcOen+JiLwR/aw3RGTJSPtqMo8VQO6xCvgU8KbTCUmBHGA3cBpQCnwZ+N3AxV1EKoFHo6+XA68DD8es3wXcB/z7Ibb/ILAyuu5pwI0icuEw6fkt8BZQAXwR+L2IVEXfCwPPApfGs2MiMh+4B7gGqAa6gZ/ELNIC/BD4dhzbGuk4fAK4GFgMLALOB24YZpOH3M840j3Y14CZwGTgDOA/ROS86LZygceAB4Ay4FfAY9HXjZeoqusfwA4iF5vVRC4+PyfyJXkG6AD+DJTFLH8hsA5oBV4E5g7a1m3RbbUR+UL7gUKgh8gFpzP6GEfki/Y74NfRz1oHLE3ivr4CfPQwlr8/muaeaJr/A5gCKJATXeZF4BvA36LLPEHkIvQboB14DZgSs805wAoiF8tNwIdSkMergUuj/38C+FvMewN5M2fQOmcBO4bYVjcwL+b5I8Adh/jcWUAfUBzz2svAJwctlxM9plNG2I9vAQ/GPJ8O9MduP/r6x4EXR9jWsMchmp+fiHn/Y8CrR7Kf8aY75v29wDkxz78OPBT9/5zo+xLz/i7gvGSfR/ZIr0cm1YAuBc4m8kW6gEjh8wWgkkhN7zMAIjKLyC+9W4Aq4GngiUG/vj4EnAdMJfLL8aOq2gUsA2o1Er4qUtXa6PIXAg8Bo4DHgf93qESKyGoRaT3EY7hflEdEVa8h8uW+IJrm7x5i0cuJ/LodT+Ti8nfgF0R+WW8AvhpNfyGRwudBYDRwBfCT6C/kdxCRnwyzv6vj2QcRqSaSr+uiL80nUiMc2McuYGv09Xj8EPiIiPhEZDbwLiI/UoYyH9imqh0xr606jM8aanuxad9K5EI+KwHbGnwcDnqfQekWkSdF5PaYZYfbz2HTLSK3i8iT0f/LiPw4O9RnzwdWq2rsOGCrOfJjalwqx+kEJND/qmo9gIi8DDSo6lvR58uBM6PLfRh4SlVXRN/7PvBZ4N1EagIAPxooXETkCWCk+PQrqvp0dPn7iRRuQ1LVRYe/aynxi+hFBRF5hkgN4c/R548Q+QULkTDODlX9RfT5myLyB+CD/KuAOEBVP0UkdHhERMRHpCb2K1XdGH25CGgctGgbUBznZp8kUmO9DcgG7lTV1w6xbFF024M/a3ycnxXv9uJN++BtDXccBn9WG1AkIqIR58eRrvEjvF8MoKqxIcOimPfjSdfg941HZFINqD7m/54hng98KcYBOwfeUNUwkfaG2AvKvpj/u2PWPZTBy/tFxG2Fe7zHbzJwYmxNBrgKGJPoBIlIFpEQYj9wc8xbnUDJoMVLiIRAR9pmOZE2mzuJhFYnAueKyKei76+TSA/BThE59Sg/69SYbQ0Uzke8vSGMtK3B75cAnYNqHke6rcHvD97WwPtHuy2TwTKpAIpXLZGLKBDpKUTkIrQ3jnWPeujwQRe4wY+7j3b7h5DIIc93Ay+p6qiYR5Gq3jjUwiJy9zD7+44aU8x6wr/a8i5V1UDM2+uINKwPLFtIJGx4yO3FmAaEVPXXqhpU1T1EwqfvA1DV+TEh1pej25wmIrG/zhfH81mq+nLMtgbCS4PTPg3IA2riSPtgIx2Hg94fId0j7Wfc6VbV/UDdMJ+9Dlg0qEfeomHSZjKUFwug3wHvl0j3VR/weSKNr3+LY916oEJESo/0wwdd4AY/Pnmo9UQkV0T8gAA+EfFHawgD3cOHK2TqiVx4E+FJYJaIXBNtQ/GJyPEiMneohVX1k8Ps73Ax/7uAuUTarnoGvbccWCAil0aPyVeItClshEjNKfq6L/JU/DFtfDXR166MLjeGSFh2FUNQ1RrgbeCr0e1cQuRi+YeBZaKflRd9mhd9fii/AS6I1o4KidTEHh1oexGR7Oj6OUBW9DN9h9jWsMeBSJjxVhEZLyLjiJzrvzzC/Rw23UP4NfAlESkTkTnAv8V89otACPiMRLprD9Ru/3KIbZlM5XQviEQ8iPRcOyvm+QPA12Kefxz4c8zzS4D1ROLOLwHzh9nW14AHYp7fBzQT6UE3boj3pxDTwyyB+/hidLuxj9Oj711DTG+oIda9iEhHhFYi7R4HpTG67Y/HLP8N4Jcxz88CtsQ8nw08RaT9oZnIhWNJAvd1cjR9vfyrx2EncNWgNG0kEh58kYN76Z0+xLF6Meb99xLp2ddGJHz6U6BgmPRMiX5GD5Fef2cNen/wZ+kI+3dlND+6iHRHLo9576NDbO+Xw2xruOMgwHeJ9FZsif4f2/PsGeALh7Gfw6X7C8AzMc/ziHxX2on8ALp10LaOAd6IftabwDGJ/L7Ywx0PiZ4MxsVE5GfAI6r6J6fTYowx8bICyBhjjCO82AZkjDEmDVgBZIwxxhFWABljjEdJZKDiNSLytoi8Pui92yQywG5lsj5/2Jslp0yZojt37hxuEWOM8SRVPeTI4kfq3DMKtbklMYPUv7G670+qel4ci56hqk2xL4jIRCJDm+1KSGIOYdgCaOfOnVgnBWOMOZgcelaLo9LcEuKff5qUkG1lj918NDWX/0tk4OLHEpKYQ3DbcDHGGJOxFAgTTtTmKgeF1e5V1XuH+Mjnojey36Oq90pkapK9qroqWQXtACuAjDEmMzWp6tIRljlZVWtFZDSwQkQ2EpkL6pzkJ88KIGOMSSNKSBNWAxr506Kj/qtqQ3TWgNOITEMzUPuZQGTE+xNUdd+ht3RkrAAyxpg0EQnBpabdPTqmX5aqdkT/P4fI1CSjY5bZQWSCzaZDbOaoWAFkjDHeVA0sj9Z0cojMePtsKhNgBZAxxqSRBHZCGJaqbuPgKTOGWmZKMtNgBZAxxqQJRQl56NYXGwnBGGOMI6wGZIwxaSRVnRDSgRVAxhiTJhQIeagAshCcMcYcBhueLHGsBmSMMYfhW9/6VlK3byE4Y4wx73D33Xdz3333JW37CtYLzhhjzMEeeeQRvv71r/Pcc885nZSMYTUgY4wZwYoVK7j55pt57rnnmD59elI/K3UjwTkvaQVQeN+sZG3aGGMSImtMzYjL/OMf/+DKK6/k0UcfZfHiYQcOOGqKWi84Y4wxsGHDBi666CJ+8YtfcOqppzqdnIyTvBqQpyqSxhg3Gu4X+K5duzj33HP53ve+x/nnn5+aBCmEvFMBSl4BlMo5LYwx5kgc6gLY2NjIOeecw6233so111yTsvREpmPwjiTWgDxUjBtjMkZHRwfLli3jgx/8ILfccovTyclo1gvOGGOient7ufjii1m6dClf//rXHUiBEEIc+FxnWBuQMcYAoVCIq666ivLycn784x8TnagtpRQIeyh4lMQ2IA8dRWOMq6kqn/zkJ2lvb+fJJ58kOzvb6SR5grUBGWM874tf/CKrVq3i+eefJy8vz9G0WAjOGGM84gc/+AHLly/n5Zdfpri42NG0RKZjsALoqHnpbl5jjDv96le/4kc/+hGvvPIKlZWVTifHcywEZ4zxpGef6+X2O27nhRdeYOLEiU4n54CwWg3oqFknBGNMuvr7q33c8vlWnn32VebMmeN0cg7wWgjOxoIzxnjK2nUB/u2GVu768SiOP/54p5PjaUkMwRljTHrZti3I1R9p4VvfLOGUU53t7TYURQh5qF5gnRCMMZ6wb1+IK69u4ZbPFbHsfH/aXqOsDSgBvDSiqzEmvbW1hrnmqv1cfkU+l19VkLbXJ2sDMsaYDNLTo3zsuv2ccmouN95U6HRyTAxrAzLGZKxAQLnxhv1MmpLN7V8uRkXSNPA2QAipd+oFSWwD8k410hiTfsJh5d9vbSMrW/jGd0ehWULI6USNIDIfkBVAxhjjWqrKN7/WTl1tmPseKMfnsx/E6Sh5Ibj0rucaYzLYT/6nk9f+0c+vH64g1y+uuh55KXpkIThjTEZZ/kg3y3/fw/1/qKSwNDvtw26xVK0NKCGsADLGpFpPT5gffrudu+6vpHx0jqsKHy+yNiBjTMZ45IEulizNZfY8n9NJOWJhD/14T2IbkHcOojHGeT09YX55dyc/ub/StdefyI2oFoI7ahaCM8ak0iO/6WbRcXlMn5tnoTeXSGIB5J1S3BjjrN6eML+6u4P/+fVol197rBOCMca4yvIHO1l4bC6z5uU6nZSjYjeiJohbY7DGGHfp7Q1z/93t/PcvRtt1x2WsDcgY42rLH+xk3pI8ps/3Z0TbT8hDhWgSp2PwTjXSGOOMvt4wD9zdxnfuG5sR1xyvTUjnnT01xmScJ37bztzFecyan36zm5qRJXE6BivbjDHJ09cX5sF7WvnWz8Zl1PUmnAE1uXhZG5AxxpWe+G07sxb6mb4gPyPafsBuRE2YTIjHGmPSU39fmN/evZ//c+8Eu9a4mN0HZIxxnWcebmPmAj+zFvqdTkpCKWK94BLBSwPqGWNSp78vzMN3N/PVeyZk5HUmk9qzRmJD8RhjXOWZh/czbZ6faQsLM6btx6usDegQQiElOzvzfl0Z42aBvjC/u7uJL9w12fXXmKGouv/aeTg82QYUDiutTUEaawPRRz9NtQEa6wKRv7X9dLSGOPeKcq67fSz+Au+cEMaksxWP7Gfq3HxmLipwOilJIhkZVjyUjLwPqKcrdKAgaawN0FTbT1Nd9HldgOZ9AQqLs6ka56NyXC5V43xUTchj7vFFB577coWf3VnLLRds5pYfTGLWkkLH9scYE6n9PHJXI7ffNcVT7STJJCI7gA4gBARVdamIfB24CAgDDcBHVbU2GZ+fxBBcckrxUFDZ3xCgqe7gwiXyN1KT6e8LRwqSsblURguZeScVUTk2l8pxuVSM8ZHnH/kE/vT3J/P3Z1r55id2cM4VFXzwpjHk+Lzz68SYdPLcI/uZPCc/0vajTqcmORRHQnBnqGpTzPPvqeqXAUTkM8BXgE8m44PTqhOCqtLdEYoUJgOFS20/zfv6D7y2vzFASXlOtDCJFC5jp/lZeHIJFdECprgsG5HhC4p4Gy9PWFbOzOOKuev2nXzhQ5u5+ftTGD89s7p+GpPuAn1hHr27ns//v2kZ38HJ6f1T1faYp4VEysWkSGkbULA/TEv9vwqWprpI7aV54HltPwpUjYvWVMblUjk2lyXvKaEy+n95tY+c3NRmUNloH3f8fDorftvEVy7fxGWfHss5V1eRlWW1IWNS4YU/NDNpVj4zFlso/DBUisjrMc/vVdV7By2jwHMiosA9A++LyDeBjwBtwBnJSuCIBdCmTZuYPn06OTnDL6qqNDc3s2vXLnbt2sVTbzTSXBcpZJqjhU17S5CyKt+BgqVyXC6TZhdw7BmjDtReCopHrr2EHap+n3VFNfNOKuXHt23ntefbufHbUygf4+4JsIxJd4G+MMvvqueWH03L+HHSFEnknEZNqrp0hGVOVtVaERkNrBCRjaq6UlW/CHxRRO4Abga+mqhExRqxAFq2bBl1dXVMnz6defPmMXfuXCZMmEBdXd2Bwmbg4ff7mTRpEpMmTaK3tI+KsXlMWVBExbhcKsbmMqoql+ycEQqXhO1aclRPLeCrD8/j8btruf2i9Xzky5N51/kVTifLmIz1wqNNjJ+Zz7QlJZ647yeVIbiBzgWq2iAiy4ETgJUxizwIPIVTBdC2bdvo7u6mpqaG9evXs2HDBl599VXGjRvHiSeeyGWXXcakSZOYOHEiRUVFB9a7f/NJQ24vIxoPs4ULb5rAgveM4p7btvL6n1u59mtTKCz1ZK92Y5Im2B/msbtruemHMzw1RE0qiEghkKWqHdH/zwHuFJGZqro5utiFwMZkpSGuK2ZBQQFLlixhyZIlyUqHK01bWMSdf1zA776/my9esIaP/9c0Fpxc6nSyjMkYLz/axLjp+cxYUux0UlJCSel0DNXA8miTRw7woKo+KyJ/EJHZRAJSO0lSD7iBD00Kr/TT9+VncdWXp7H4va387I6tHHt2OZfdNom8/Gynk2aMqwX7wzx+Ty03/vdMz1xPQFI2lY2qbgMWD/H6pSlJADYUT8LMfXc5X32smN/cuY2vXbKGj393JlMWeuNXmzHJsHJ5A2Om5jN1SWlmhO7NO9ho2AlUMCqXf/vvOfzjyUZ++IkNvPfqsSy7YeKIHS+MMQcL9od56u49fPz7sz11LUlxCM5x3tnTFDrx/Cq+tHwJNa+3850rV7Nve4/TSTLGVf7+WAPVU/KZcWyJ00lJuVA0DHe0DzewEFySlFbn8+mfLeSl39Ty7Q+v4kuPHUfZmDynk2VM2gsFwjx1126u+94cz19HMl1aDcWTcbLgPddMZO0r+9m+rouSMflOp8iYtPe3x+qpmlzA1GPLPHHfTyxV8VQILnltQNZn/4Ax0wup29LNwvfaMTFmOKFAmGfv3sk1357r2WuIl2p93tlTB1VPK2Tf1i6nk5H2ams6+fPPd7F/X6/TSTEO+efj9VRMyGf6caOcTopJAQvBpcDo6UW8/Nu9dkyG0LU/wJtP7+Off6yjvaGfGSeM4vn7dnHhv8/k+IvGjDguoMkcoUCYP929gyv/a55nvyuKt3oQJzEE580TaChV04qo39ZNMCQ2gjaRC83Gvzbz2h9rqflbC3NPq+R9n53JzJPKycoW9qxv57dfWMeq5xq57GtzKamyzhte8NrjdZSPL2DqcRWODTjsPPFUCM4GL0uB/GIf/qJs2vb1UjbOux0Rams6eG15LW8+uY+Kifkcf/E4PnznPPJLfActN2FeCZ/73Yk8d9c2vv+BV7nkjtksWVZttaEMFg4pK+7ZzuXfnO90UkwKJTEEZxeLWKOnF1G7tZuScZk6l/3Quvb389ZT+3j9j3vpbOnnuAvH8clfH8/oqf+a12Wonk6Sm825n53JvPeO5qE71rLquQYu+cpcispt+otMtPmfLeQV5jBlabnner7FityI6p1rp4XgUmT0tCLqt3Yz65TMPy6hQJhNLzfx+mN72fqPZuaeVsV5n5vNjBMryMqOfLniDbGMX1DGZ37/blb87xb+++K/c/GX5rLg7DFJTL1xwhtP7OOYC8bZdQNvtZ9bDShFKqcXUbehI6OPS93Gdt58bC9vP1VL5aRCjr14PJd+YyH+4kiITYl/KvRYWXk5nHvbHOacWc3vv7ia1SvqueAL8ygYZbWhTNDfE2LdX+o5+5ZZGf39MO9kbUApMnpaEauerHM6GQnX2dLHqqfqePOPe+huDXDsReO54f6TqJyc+KmTJx9Txqf/cAp/+uEmfnTJK1z81QXMOX10wj/HpNaGF+qZuGAUJVV+p5PiuATPiJr2LASXIpVTS2jY2kkoLK5vTA8GwtSsbOCtx/aw/bVmZp9WzXm3zWPqCRUHevklqxdTjj+L99++gHlnjuXRL69izYp6LvrqQnJ8dr651dtP1LHo/PF2zYjyztQTNhZcyvjL/SDQ3hSkqNKd3YrrNrbx9mO7WfP0XiqnFLHk4olc/M1j8BfFhNhS1H120tIqbvzD6fz0ypfZtaqNycfZtOhu1NXSx863Wrj0e8fZNcODbDqGVBGhcloxDds7Kah0T6ihq7mP1U/tYdVju+ntDLD4golcf/8plE/61/TrYYfS5ivwUTGliI6WfjvfXGrNs3XMfE81vgKfY+dROlHFU1OPWxtQClVNK6ZxWwdTjq90OinDCgXC1KysZ9Ufd7PjjSbmnDGGc/9zAVOWViBpdiNtYXkeXS19TifDHKHVT+7m9E/NcToZacXagBLAqtPvVDGtONIOlIbHRlXZt6GN1Y/vYt0ze6icXsyiCydx4X8dR15hJMQWhkicLY3kl/vpbA6k5TE1w2ve0UlrbQ+TTxxt+edRNhp2CpVPLWHzyvq0Ojadzb2se2o3qx/fRX9XkIUXTuLaB06nbMK/erGl87AoBWV5NG3rSKtjauKz5qk9zFs2AbKz0/ocS6VILzjvFMY2GGkKlU8rpWlbh+PHJtgfYuvKfax5bAd73mpm5hnjOOs/lzDxuMoDITa33I3uL8+n8/Vmx4+pOTyqytqndnPR9060vBvES/dCWRtQChWPyaevM0BfR4C8Yt/IKySQqrJvfStrH9vBhuf2UDW9lAUXTebC75xIboF7T4PCijy6rQ3IdfauaiHLl0X1XJt2IZYNxZMgXjqI8RMqppbQuK2DcYtS0224s7GH9U/vYt3jOwn0hlhw4WSueeBMSse7I8Q2kvwyP90tfXa+ucy6J3cx//2TULJQF59/5ugksRu2VauHUj4tUgCNWVSVtM8I9oXY+lIt6x7fQe2qJma8dzxnfuE4xh/zrxBbpnR59ZfnRwogO99cIxQIsXHFHq7+zdmWb+9gbUAJ4aW+7Idj1JQSmrd1JPz4qCr161pY/8QOap7bRdWsUcy7YCrv++678eVHsjkde7EdLV9xHv3dAfr7w2T7sp1OjonD1pfrqZhaStG4opTduOwmXrqnzb3Bf5eqmFbK2uVbE7a9zoYeNj69g/WP7yAUDDPv/Clc+ZtzKBmX+LHY0pFkCf5RefTs76NotLemunCrjU/vYM77JzudDJMGrA0oxUZNLaF5W/tRHZ9gX4htL+5lw5Pb2be6ielnTuS9XzqesUsqD4wz5+Z2ncNVUOans6mPgipvFLpu1tfRz85X6znjSyfYNWIINhJCgngpjnk4iseV0NXUS39PmBx//IdfValf28zGJ7ax+c+7GD2nnDnnT+W875x6IMSm4MkG3YLKfDqb+qi0cy7t1azYw4QTxpBb7PfUj6TD4aVrp4XgUiwrJ4vqBRU8fvMLLPrwbKaePoHsYUZy7mzoZtNT29n45HbCoTBzL5jG5b9ZRvFY+7U/oHB0Pl2NPU4nw8Rh0zM7WHz5bKeTYdKETUjngPN//F62v7CbVQ9v4uUfvMHcD8xg3iUzKazKByDYG2T7S3vY9MQ2GtY1M+3MSZz+lZOoXvSvEJtbbhRNhYLKAjoae+ycS3MddV20bGlj4injLa8OweYDShAvHcTDJTk5TDt7KtPOnkrz5v2se6SGhy57kgknjSWvyMe253dRNbeCWRdM4+zvnYbP7+0Q20jyqwpo3tRi51yaq3lmB1PPnIT4ciz8NgzrBZcAXopjHo2yGRWccse7OOHTx1Hz1FaCvSE+8JsLKBqTGTeKpkJ+ZQFdr9TaOZfGVJWap7dz6hdOsnwyB1gbUJrILcplwYfnOp0MVyqsKqCrsdvpZJhhNNe0EOwLUb3YplAfjg3FkyBeqkYaZ/mrCulu7LZzLo1tfno7M86bhkpWpt0LnXBeqiHaSAjG9XJH5dPb1kcgoGTleOfL6xbhYJitf9rOeXefZ9cFcxBrAzLul52FvzyfrqZeCquLRl7epFTta3UUVBdSPKnM2jNHotYLzhjXKagsoLux2wqgNLTt2S1MO2+608lwBcVbzRfWDdtkhPzKAroae6iw8y6tBHsC7H55N8fcfKJdE8w7WCcEkxHyqwrpbuqx8y7N7HxpF5ULR5NXUZAxU4Akm5cKaqsBmYzgryygq6Hbzrs0s+PZLUxZNtPyJU5e64ZtPQVMRsivLKSnqcvpZJgYPc3dNK1rYPx7pjidFJOmrBecyQj+yiJ6GnvsvEsjO1dsY9wpU8jKy7Xeb4fBSzUgC8GZjJBXWUhPY5edd2lk57NbWHCjdT44HDYYaYJYY7BJpdzKInqau+y8SxPtO/bT09RF5bHjLU/MIdl9QCYj5JbkEeoLEewNkOP3OZ0cz9v9pxomnjMTybaQ6OHyUoFtITiTIQR/RSHdjT0UTch1OjGepmFl13ObOem/ltl14HCpt66dVgCZjOGvLKSnsZuC8aOcToqnNa+uI6fAR9H0SrsOmGFZCM5kDH9VIb3WFdtxe/60iQnnzD4we6+Jn9fuA7IakMkYeRWFdFtPOEeF+kPUvrSV9/z8csuHI+Sl42YFkMkYeRVF9DZZAeSk+r/vpHhaJXmjS+zeHzMi64ZtMkZuZRGtNY127jmo9vkaxp452/LgCKX6PiAR2QF0ACEgqKpLReR7wAVAP7AVuE5VW5Px+dZH0mSMvMpC+po6nU6GZwW7+2l6bRfVp9nUC0dDVRLyOAxnqOoSVV0afb4CWKCqi4Aa4I5E7+MAC8GZjJFrIThH7Xt5O6MWjiOnpMDCby6mqs/FPH0V+GCyPssKIJMxfBVF9DV3EQpjPbAcUPf8JsacOce++0cpgeHLShF5Peb5vap676BlFHhORBS4Z4j3rwceTlSCBrMCyGSMLH8ekpNFf0c/vmK/08nxlP62HtrW1jH/K++37/5R0MTeiNoUE1Y7lJNVtVZERgMrRGSjqq4EEJEvAkHgN4lK0GDWBmQySl5FkbUDOaDl9Z2MWjKBnHwbhcJNVLU2+rcBWA6cACAi1wLnA1epatICqlYDMhklt6KQ3qZuCqbY+ZdKPXXtFEwqt+99AhxmB4IjJiKFQJaqdkT/Pwe4U0TOA/4TOE1Vu5OZhqQVQKk6iMbEyqsopq+p086/FOupa6N49hg77kctpd2wq4Hl0fbSHOBBVX1WRLYAeURCcgCvquonk5EAuw/IZBRfRRG9Ni1DyvXua6fytJBovsQAACAASURBVDl23F1EVbcBi4d4fcZI64rI43F8RIuqfnS4BWwsOJNRfCV++pqtDSjVeve1kTemxOlkZASX1CLnAh8f5n0BfjzSRqwNyGSUzm1NlCyaYOdfCmkoTF9TJ7lVpXbcj5KLBiP9oqq+NNwCIvJ/RtqItQGZjKGqtK3azYSr323nXwr1NXXhK81HfD6S11/KpBNV/V3scxEpVNWu4ZYZitWATMborW1Fw0ruWOuNlUo9dW3kVVvtJyEUVxXiIvJu4GdAETBJRBYDN6jqp+JZ3+4DMhmjfc1uShZNtFEQUqyvvo28MTYJYKKEkYQ8UuT/AucCzQCqugp4T7wrWwjOZIz2VbsoXjjJzr0U661rI3d0qR13j1LV3YN+9IXiXddCcCYjqCrta3Yz5vKT7dxLsd59bRQvnGTHPQEU1/143x0Nw6mI5AKfATbEu3ISa0DJ2rIx79S3r41wIETeuHI791Ksr76NirNG2XFPiNTOB5QAnwT+BxgP7AGeA26Kd2W7D8hkhI7VkfCbtf+kXn99pBOC8R5VbQKuOtL1bSQEkxHa1+6iaOFkO+9SLBwIEdjfRU5VqR37BHFTTVJEZgF3AdWqukBEFgEXquo34lnfOiGYjNC5ehfVl9r9P6nW39iBr7wIsrJddeFMZy47h38K/DtwD4CqrhaRBwFnCyCXxTGNi/XXtxLuD+KbUGnnXYr17mvFVz3Kjrt3FajqPweFvoPxrmxtQMb1OtfupHCBtf84ob++jdxquwcoUVRdVwNqEpHpRDrwISIfBOriXdl6wRnX61yzi8IFk+2cc0D/vtboPUBOpyRzuKw2eRNwLzBHRPYC2zmMTgnWBmRcr3PNTiovfpedcw7oa2ij+Njpduw9SESygKWqelbs5HaHsw0rgIyr9Te2Ee7tJ3dClZ1zDujf14qvapQd+wRyS21SVcMicjPwu8EDkcbL2oCMq3Wt3Unh/MnW/uOQQEMrudVlTicjo7isMF8hIrcBDwMHCiFVbYlnZesFZ1yta80OCuZPsfPNAeG+AKHOHrLKSuz4J4gibiuAro/+jR39QIFp8axsnRCMq3Wt20nZBe+y880B/Q1t5FSWgogdf49S1alHs761ARnXCjS1Ee7qtfYfh/TXt+Ibbe0/ieamslxEPjDEy23AGlVtGGl9awMyrtW9bgf586cgWTatlRMCDZECyCSQ++4D+hjwLuCF6PPTgVeBWSJyp6reP9zKVgMyrtW9bicF86fYueaQSA2ozI6/t4WBuapaDyAi1UTGhjsRWAk4VAAla8PGRHWv3cGoZSfaueaQQEMredPm2vFPNHcd0CkDhU9UAzBLVVtEJDDSylYDMq4UbG4n1NmDb2K1nWsOCTTsJ6fKakCJ5rLj+bKIPAk8En3+QWBl9MbU1pFWtjYg40rd63eQP8/af5wUbIyE4Iyn3QR8ADgFEOBXwB9UVYEzRlo5eQWQu6qRxmV61m4nf94UO88cEu7pI9wbILuk0PIgwdzUpV1VVUReB9pU9c8iUgAUAXENyWMhOONKPet3UHruiXaeOaS/oZWcqlFAlqsumOlOcde1U0T+DfgEUA5MJzI1993AmfGsb/EL4zrB/R2E2rvJnVTtdFI8K9hg4TcDREJwJwPtAKq6GRgd78o2EoJxne512/HPnQxiv76d0t+wn5yqUXb8E00BF9WAgD5V7R8Yi1FEcjiMoKyF4IzrdL+1Gf/cqXaOOSjY0EpOVbnlQRK4rFB/SUS+AOSLyNnAp4An4l05iZ0Q7MQ0idf29F/prdlFxVXL7BxzUKBhP3mzJlsemNuJjIawBrgBeBr4WbwrWzds4xrtz79G2zN/Y+xXPk52aZHTyfG0YON+fFXWBpQULqoBqWoY+Gn0cdisDci4Qucrb7P/0b8w9ksfJ6eyzM4vhwUbW8musDagxHPHdAwisoZhikpVXRTPduw+IJP2ul5bR8uDzzDmjuvxVVfYueWwcE8fGgyRVVRgeeFd50f/DswDNDDm21VAd7wbsU4IJq11r9pE032PUf0f1+EbP8Z+caeBtmf+in/hDOweoCRxwTFV1Z0AInKyqp4c89btIvJX4M54tmP3AZm01bN+G033/J7Rn7uGvCnjnE6OAfp376P9ub9T8ZELnU5KZopOx5CIR4oUisgpA09E5N1AYbwrWwjOpKXezbto/N/fUvXpK/DPmGTnUxrQUIimnz5K2WXnklNeanliINID7j4RKSVyRrTxr2m6R2QhOJN2+nbspeGHD1B5w2X450y3ME+aaHvmb2T58yg87Xj7fieTi853VX0DWCwiJYCoatvhrG81IJNW+vfU0/iDX1Fx7cXkL5xt51GaCOxrpP2plxjztZsQxPIlqdK/cBeR81X1yYHnqto+0jJDsfuATNoI7Gui4fv3UXbF+yhYOt/p5JgoDYdp/vmjlF70XnxV5U4nx6SH74nIXoYvLb8FOFUApX8pbtJHsGk/9d/9OaUXn0Xhu45xOjkmRudf/glhpfisd2Pf6xRwR+2yHvjvEZbZPNJGLARnHBfc3079d35GyXmnUnzaCXbupJFgYwutf/wzY+64AZEsy5tUcMExVtXTE7EdK4CMo0LtnTR872cUnbqUkrNOtvMmjagqzb9cTsk5p+AbO9ryxiSctQEZx4S6uqn/wc/JP3Y+peePOHuvSbGuV94g3NlNyXnvcTop3uG+6RiOio2GbRwR7umj4b9/gX/WNEZdcq6dL2kmuL+d/Y88Q/XnP45k51jtJ4W8dNuBDUZqUi7c10/jj36Jb8JYRl1+AZEBGJ1OlRmgqrTcv5yi00/CN3Gc5Y05JBEpAD4PTFLVfxORmcDskbpfD7CheExKaSBI00/uJ7uslPJrLmFgJkWTPrr/uYpgQxOl73+v00nxJk3QIzV+AfQB74o+3wN8I96VrROCSRkNhmi650EkN5eK6y6zXlVpKNTRyf6HnqDq5muRHAu9OcJd4ejpqvphEbkCQFV75DB+VVobkEkJDYdpvu93aCBI1U3XIll2cUtH+3/7BIUnHUve1MmWPyYe/SKST/RsEZHpRGpEcUlaASR28pooVWX//Y8Sautg9KevJ8satdNS99vr6N+xm7Ff/px9fx2UymMvIjuADiAEBFV1qYhcBnwNmAucoKqvD7OJrwLPAhNF5DfAycBH4/1864ZtkkpV2f/w4wRq6xl9y8fJyvU5nSQzhHB3Dy0P/pHKj11OVl6u08nxrtS23ww4Q1WbYp6vBT4A3DPSiqq6QkTeBE4iMkzGZwdta1jWBmSSqnX5s/Rt2UH15z5BVl6enRdpquWRJ8lfPA//rOmWRx6nqhuAYTsIicixg16qi/6dJCKTVPXNeD7L2oBM0rQ9/Tw9q9dTfeuNZOXb9M3pqmd9Db0bNjPuy5+3763jJNV5oMBzIqLAPap6b5zr/WCEbcbVhdJqQCYp2p9/mc5XX6f61hvJLiy08yFNhXv7aP7N76m48lKy/H7Lp3SQuDyoFJHY9pt7hyhgTlbVWhEZDawQkY2qunLEJKomZOgSawMyCdfx8qu0/+Vlxnz+RnJKS5xOjhlG62PP4J85jfz5c5xOikm8JlVdOtwCqlob/dsgIsuBE4ARC6ABIvKBIV5uA9aoasNI61sNyCRU5z/eoO3pP1P9uU+SU1Zm50Ea692yje631jD2S5+3fEonKcoLESkEslS1I/r/OcCdh7mZjxG5CfWF6PPTgVeBWSJyp6reP9zKVgCZhOl6azWty59i9GduwFdZaedAGtNAkOYHHqH8QxeTXWDtc2kldXlRDSyPdjbIAR5U1WdF5BLgf4Eq4CkReVtVzz3ENsLAXFWtBxCRauAu4EQiNSmnCiBrzPSSnnUbaHl4OdU3/Ru5Y8bYBS3N9WzYTHZxMQWLF1leeZSqbgMWD/H6cmB5nJuZMlD4RDUAs1S1RUQCI61sbUDmqPVs2kzT/Q8x+obryZ0w3unkmDj0rN9I/sJ5TifDDOa+6RheFpEngUeizy8FVkZDeq0jrWwjIZij0rttB02/eIDR138E/xQbvsUNVJXedRsYfcP19j1NQy7Lk5uIFDonE7kR9dfAH1RVgRF7ylkbkDlifbv20PDTX1B59RX4Z9gNjG4R2NeAquKzUKk5StGC5vfRx2Gz6RjMEemvraP+np9TcfllFMyzLrxu0rN+A/nz5thUGOnKRdMxiMgHRGSziLSJSLuIdIhIe7zrWwFkDlugoZH6u35K+QcupHDRAqeTYw5Tz/pNFMyb63QyTGb4LnChqpaqaomqFqtq3Df/WRuQOSyBlhb2/fgeypadS/Gxx1gIx2XCvb307dxF/swZ9h01iVA/MHbckbBu2CZuwbY29v34HkrPOJ3ik06ywseFejZtIW/yJLJybdiddOWyHwavi8jDwB+JmQdIVR+NZ2Xrhm3i1vrcnylYuIDS95zqdFLMEerZsNHCb+nOXT/eS4BuIqMoDFDA4QLIXaW4GYGGw3SvXsPYT99seetSqkr3+g2MOf00y0OTEKp63dGsbwWQiUvv1u1klxTbEDsuFqjdh+Tk4KussjxMV85MSHfERMRPZDy4+YB/4HVVvT6e9a0TgolL96pVFC5abPnqYj0bNlAwZw5ZiKsucp7jrry5H9gInEtkINOrgLg7JVg3bDMiDYfpWrOGwsXvGDbKuEj3hg0UzLX2n3QnmphHisxQ1S8DXar6K+D9wMJ4V7YQnBlR7/YdZBcWkmuhG9cK9fTQt3cv/mk2YoVJqIEBR1tFZAGwD5gS78pWAJkRda1eTeHCxZanLtazqQb/1Klk+XItH9Odu/LnXhEpA74EPA4UAV+Od2VrAzLD0nCYrtWrGPeJT1qeuljPxg0Uzp5reegGLsojVf1Z9N+VwLTDXd/agMyw+nbtJCu/gNzR1U4nxRwhDYfp3rSRgjnW/mPSi42EYIbVuXoNRQsXWX66WH9tHVl5fnzl1oU+3aW4A4HjrA3IHFKk99sqxl7/CctPF+vesIGC2XMtD93CQz/2rA3IHFLf7t1Ibi65VdWWny7WvWkD5Weda3lokkJE3k2k59uB8kRVfx3PujYWnDmkzrWrKFq42OaNcbFQVxf99fvInzrd6aSYeLnoh4KI3A9MB94GQtGXlcjMqCOyEJwZkqrSuWYVY6/5mOWli3Vv3kT+tOlIdo7lo0u4rKa6FJgXnRn1sFkIzgypb89uJCeHvOqxlpcu1r1pAwUzrfu1SZq1wBig7khWthqQGVLn2lUUzV+M2LhhrqXhMN01G6k4c5nloZu4IK9E5AkiKS0G1ovIPzl4PqAL49mOtQGZd1BVOtauYtyVRzXSunFYX+1usouK8ZWVO50UEy/3dMP+fiI2YjUg8w59tXsREXKrx1k+uljXpg0UzrLu1ybxVPUlABH5jqr+Z+x7IvId4KV4tmNtQOYdOte+TdH8xTZsv8t11Wyg6pwL7LvoNu7Kr7OB/xz02rIhXhuSDcVjDqKqdK5bTfF8m3rBzYKdHQSaG8mfNNXppJjDpQl6JJGI3Cgia4DZIrI65rEdWB3vdqwNyBykv74WNEze2AlOJ8Uche4tGymYNhPJznY6KSYzPQg8A/wXcHvM6x2q2hLvRqwNyBykZ/cu8qfMsN5vLtdVs5HCGdb+40YuCZmqqu4QkZsGvyEi5fEWQtYGZA4SbN2Pb1S55Z+LaShE19ZNjD7nQstHkywPAucDbxD5mRM7XIoS59QMVgMyBwm0tlA4fY7ln4v17NmFr7SMnOJSy0eTFKp6fvTvUTUyWhuQOUikBlTmdDLMUejavCESfjPu5LIfDSIyHpjMwYORroxnXasBmYMEWlvwlZRZ/rlY15YNjD7vEstDN3LPjajAgXt+Pgys5+DBSJ0tgNx0EE2EhkIEuzrxFZda/rlUoKONQGsLBeMnWx6aVLgYmK2qfSMuOQSrAZkDgh3t5BQWIZJt+edSXVs2UjhttuWhm7kr37YBPmLGgTsc1gZkDgi0t5JTUup0MsxR6NqygaJZ851Ohjka7iqAuoG3ReR5Dh6M9DPxrGwhOHNAqL0NX/EoyzuX0lCIru01jDnvUstDkyqPRx9HxEJw5oBAeyu+4lGWdy7VvXs7ueVV5BQUWx66lOCuH++q+isRyQcmqeqmw13fCiBzQLC9lRwrgFyra9smCqfaPVyu56L8E5ELiEzNkAtMFZElwJ3xzgdkg5GaAwIdbfiKrQ3Irbq211A0dZbTyTDe8jXgBKAVQFXfBuK+OdXagMwBwWgIzvLOfYLdnfTvb6RgnHW/djWX3QcEBFW1TSR2JJ7463AWgjMHRGpAFoJzo64dmymYMA3JyrH8czt35d9aEbkSyBaRmcBngL/Fu7IVQAYADYcIdnWQU1BieedCXdtrKJoy2/LOpNqngS8S6YL9IPAn4Bvxrmz3ARmASOGTX2jzx7iQqtK5YxMVJ5zudFJMIrjrR8RxwFdU9YsDL4jIscCb8axsbUAGGGj/sSF43Ki/uQEB8spGW/5lAJfl4Z+A10TkQ6paH33tZ8Cx8axsITgDREdBsPYfV+rcUUPhlNk2iaBxwibge8CLIvIxVf0bB88NNCwLwRkgpgOCcZ3OHZsYNe84p5NhEsVdPyJUVZ8UkU3AwyJyH+nQC85l1UjPC3a04iuyEJzbhENBuvdsY8K5l1veZQLFbQWQAKjqZhE5BfglsCjelS0EZ4BIDSi/epLlm8v01O4kt6ySnPwiyzuTcqp6TMz/3cCHRGRSvOtbAWQACERrQJZv7tK5YxNFk6z7dSZxe01WVXfFu6wNxWMACHa24iuyNiC36dxVQ9FkG34no2iCHnEQkR0iskZE3haR16OvlYvIChHZHP1blsC9O0jy2oCStWGTcBoOE+zqwFdUYvnmIsHeLvpa6ikcO9XyzRyNM1S1Keb57cDzqvptEbk9+vw/B68kIlnAB1X1d0f6wRaCMwS7Osj2F5Blw7i4SteuLRSOm0ZWtuVbJkmDENxFwOnR/38FvMgQBZCqhkXkZiD9CqA0OIgmTpEecDYIqdt07txE0aRZlm+ZJnH5WTkQVou6V1XvHeLTnhMRBe6Jvl+tqnUAqlonIqOH+YwVInIb8DDQdWCjqi3xJNDuAzKRDgg2DYOrqCodOzdRueQ9TifFpK8mVV06wjInq2pttJBZISIbD/Mzro/+vSnmNQWmxbOyheAMgc42fIU2CoKb9Lc2oeEQeWXVlm+ZJMX3AalqbfRvg4gsJzK3T72IjI3WfsYCDcOsH/fcP0OxAsjQHw3BWZ65R8fOTRRPsuF3Mo2Qug5cIlIIZKlqR/T/c4A7gceBa4FvR/8+NsJ2FgDzAP/Aa6r663jSYG1AhmBnKwWV4y3PXKRz1yZGzVhieWaORjWwPDqZXA7woKo+KyKvAb8TkY8Bu4DLDrUBEfkqkQ4L84CngWXAK4CzBZBxj/7ONrsHyEU0FKJz71YmnPEhp5NikiFFPypUdRuweIjXm4Ez49zMB6PbeEtVrxORaiKjYcfFQnCGQGertQG5SNe+neSWVODLL7Y8y0Auq9X2RLtjB0WkhEh7UVwdEMBCcJ4XuQm1ndxCG4jULdq3r6VkonW/Hk6ov5f+9hayfD6yfH6yc/OQbB/RcJNJnNdFZBTwU+ANoBP4Z7wrWw3I44LdnWTn5dtNqC7R01xLy8bXmP3BWy2/gFCgj779DfS01NHbsi/y2L+PUG83ucVlhENBQv29hAN9aDhMti+PrFz/gb/nbnmW4uJiSkpKhvx7qPeSykX5qqqfiv57t4g8C5So6up417c2II8LdEXDbybtaSjErr88xLgT30ducdKG50pL4WCAvtYGegYKmWhBE+hqI2/UaPLLx+AvH0Pl/HfjrxhLbnEZkZFiYrYRChIO9BHq7yMc6CXU38fnbngf7e3tdHR00NHRQXt7O42NjWzbtu3A64Pf7+joSO7OuqAAik67fcj3VNWm5DYjC3S0kmvzALlC/Vt/IcdfSMWckzI2vzQUoret8UAB09NSR+/+evo7WsgrqcBfPhZ/WTUVs4/HXz6GvJIKJCv7EBs7+Gl2Vg7ZeTn48goPvHbeeecdUTotlMcPhnlPgffGsxELwXlcv3VAcIWe5loa16xk9gduzYh7fzQcpq+9id79kdpMT/RvX3sTuUVl+MsiNZqyaUsiBU1pVWTMuyE3ltq0J5W648e7qp6RiO1YCM7jAl1tFoJLcxoKsfPFhxh7QmaE3nqaa9n06P/FV1iKv2wM+eVjKJ00j+olZ+IfVUVWTq7TSXSWCwqgASLiA24EBsaEepHImHKBeNa3GpDHBTpbyS8fa/mVxurfjobeZp+UEflU//ZfGHv8MqoXHyJKkwH7eDTcUAOKcRfgA34SfX5N9LWPx7OytQF5XH9nK7kFNhJ2uuppqaNx7UrmXHIrWRkQeuvv3E/77o1MOvlSO+cyw/GqGnsz619EZFW8K1sNyOMCXW34Cmwq7nSk4RA7X/ot45a+j9zCsozIo4a1r1Ax83iyffkZsT9J4a7jEhKR6aq6FUBEpgGheFe2NiAPU1UCPe34CpJ8X4M5IvtW/YXsvGjoLQOE+ntprvkncy7+nNNJSWsuqxn+O/CCiGwjMo7qZOC6eFdOYgjOXUfRi8L9vYhkkZOTC5ZfaaWnpY7GdSuZe9HnyIKMyJ/mTf+gZNwM/EVlGbE/BlT1eRGZCcwmUgBtVNW+eNfPGnmRI02ZPdL9EejuIGdgPDF7pM1DQyF2rHyIccfFhN5c/tBQiIZ1L1M9/zTH03LQI924bN9E5DIgNzr6wQXAb4e7SXUw64TgYcGeDnz+IsurNLNv9Qvk5BVQNfPEjMmb/TvXkltQQlHV5PS88KcTdx2fL6vqIyJyCnAu8H0iveBOjGfl5NWATNoL9nRGRlQ2aaNnfx3161Yy+eQPZdTd9vVrX4rUfkymGehw8H7gLlV9DIj7Ri7rBedhgZ4OfH4b0j9dhMMhtr/8EOOPex95GdLrDaCzYQfB3k5GTVyQMfuULILrokd7ReQe4CzgOyKSx2FUbCwE52HBnk4LwaWR+jXR0NuMzAm9AdSve4nquaeSJVlWAMXDXcfoQ8B5wPdVtVVExhLpGRcXqwF5WKCng/xRYyyv0kD3/jrq169k3vs/lxFjvQ3o62imY98Wpr77ctfuUygUYvfu3dTU1Bx4mAhV7QYejXleB9TFu77dB+Rhgd5OSqwNyHEaDrH9bw8x4Zj3kVfk/rHeYtVvfJnKGSeS7ctzOinDUlUaGxsPKmQGHlu3bqWyspJZs2Yxa9YsZs6cmdS0eOkWFgvBeViwp4PcvGLLK4fVrX2BnNzMC70F+3to2vo6C8+/LW32KxTopbe9id72Rno7Grn66jcPFDRZWVnMnj37QEFzxRVXMGvWLGbMmEFhYeFB2/n85z+fnASma/fwJLEQnIcFe7vIySu0vHJQd2sd+zasZP77Miv0BtBY8yqjxs8jtyC1032EQ0H6OpsPFDKxBU6ov4e84kr8JVX4S6o4++yzuemmm5g1axYVFRWpS6QBrAbkacH+Lny5hZZXDhkIvU1c8j78BZnT6w0iPfrqN73MrNOuT8r5pRqmv7stppBppCf6t7+7jdzCUeQXRwqZwrLxVExejL+4ityC0oNmSr322msTn7ij5KXvo7UBeZSGwwT7e8nJzXc6KZ5Vu/4FsnMLqJoe1z17rtKycxX+4ioKyycc8TZUlWBf14ECprejiZ5ogdPX0UR2bj7+kqoDBU3JmJn4S6rIKyw/9OR1bmAFUAJ4qCHNjYL93WTn5EV+DVpepVx3ax37Nq5kwXm3IJBReaCq1G14kQmLzo1rv0LBPno7mg4UMgdqMx2NoOAvqcRfHCloKiYtwl9chb+4cviODRl0PDOZheA8KtTbjS/Pwm9O0HCIba8+zMRFyzIu9AbQ0bCNcLCfsrFzDpxf4XCIvs6WSG2m41/hst6OJoL93fiLKg8UNCXV06mecRL+4ipy8goPPSJEhh23AV76TlonBI8K9nWRk1dg+eSA2vUvkuMrYPS0EzPy+NdteIm8wnJ2vvXkgQKnr6uV3ILSSLisuIqC0rFUTFg0ZLvMO2TgMRqWh/bXxYFSczSCfV3k5BY4nQzP6W6to65mJQvPuSWjxnqL5fNH7i3LzS+hZPR08oszoF3GJEXyQnDhZG3ZJEKwt5uc3ELLpxTScIit/3iYSQuW4c8vgww99tOXXjb0Gxm6vwmlFoJLDA8dRDcK9kW6YFs+pU7txhfJyS1g9NTMDL2ZBPHQuWGdEDwq2NeNL7fA8ilFutrqqK1ZyeKzbiErw244NeZIWVDWo4L9XeQVjHI6GZ6g4RBbXn+YyQuWkVeQWWO9mcRy4XQMR8XuA/KoQF8XObmFlk8psHfTi+T48hk95QQ73mZkHjpHLATnUcF+C8GlQlfbPmo3r2TJez9roTdjBrFOCB4V7OvC57P7gJJJwyG2vPEwk+ctIy8/8244NcnhpR+F1gbkUYH+SDdskzx7N79ETm4+1VNOcDopxi1sOobE8FIp7jaqSrC/m1yfheCSpat9H3u3rGTJ6RZ6M+ZQrBOCB4UCvUhWNllZOZZPSaDhEJvffJgpc8/Dnz/KjrE5LF66OdxqQB4UsnuAkmr35kivtzGTTrBjbA6fh86ZYUYANJkqEIiOgmASrn3/Lmq3/ZVZSy7L2LHejEkU6wXnQYG+bnKsB1zCBYO9bHrjt8xYdAl5/tROQ20yh5dqzRaC86DIOHAWgku0bWseZ1TFdKrGLLTCxxwZxVNthsmrAYW9cxDdJtjfHbkHyPIoYRrrVtG+fwfHnvxZO67GxMnuA/KgQH8XOT5rA0qU3p79bFn3GAuWXk92Tq7TyTEu56XIhLUBeVCgv5uCwirLowRQDbPx7YeZMOU9FJdOsGNqjp6HziFrA/KgYH8XvlGTLY8SYNfWF8mSLCZOfY8dT2MOk3XD9qBAoNu6UhLNRgAAEZ9JREFUYSdAe+su9u78K3MWfQgR+yqZozcwHUMiHm5gIyF4ULC/C19OvuXRUQgG+9iw6iFmzruYvLxSO5YmMVQ9dS5ZCM6DAoFILzjLoyO3ZcNjlJVPY/ToBZ6K2RuTSNYJwYMC/d34cgotj45Qw75VtLfuYulJn7FjaBLOSz8MLXDtMaFQANUw2dnWXfhI9Pa0snnj48xdcLkdQ5McmqCHCyQxBOeSI+Axwf7IRHRZ4KlYcyKohtmw9iEmTj6V0pLxdvxMRhCRbOB1YK+qni8ii4G7gSJgB3CVqrYn47OTOBJC0rZsjkKwb2AUBKdT4j47d7yEkMWkiafa8TNJ40AI7rPABqAk+vxnwG2q+pKIXA/8O/DlZHyw1YA8JtDfGe2AYPlzONradrNn9185funN0Qnm7PiZJFBSOpSTiEwA3g98E7g1+vJsYGX0/xXAn0hSAWRtQB4TCEZrQCZuwWAf6zc8zOzZF+H3lzqdHGPiVSkir8c8PjHEMj8E/oOD6/RrgQuj/18GTExWAq0XnMdEesDZVAyHo6bmCUaNmsboSutybVIgcedYk6ouPdSbInI+0KCqb4jI6TFvXQ/8SES+AjwO9CcsRYPYjaheo0pPTwsaDtnd+3Gob1hNe/sujj/uJjunTUqksA3oZOBCEXkf4AdKROQBVb0aOAdARGYRCdElhd2I6jHjxyyloXEtW7Y+w6zpSTuvMkJvbys1W55kyYKPkJOVa7Ufk1FU9Q7gDoBoDeg2Vb1aREaraoNEfqF+iUiPuKSwn8Aek5WVw6L5V9Gyfwu79vzV6eSkLdUw6zY+wqQJp1BSPMHp5BgvGRiO52gfR+4KEakBNgK1wC8Ssl9DsBCcB/my/SyZ/xFeX3Uv/rxSRlfOdzpJaUU1zMYtj5GVlc3k8SfbuWxSyonokaq+CLwY/f9/gP9JxecmLwRn90mktfzcUSyeezVvrf8VeTnFjCqZ5HSS0oKqUrPtKbq6Gjlm3rVkaZaF3oxJEqsBeVhJ4VgWzLiU1Rsf5LgFH6Mwv9LpJDlKVanZ8QztnXs4dt5HycnOtfPYpJaLhtFJBGsD8riKsplMn3QWb6//Nf39nU4nxzGqypZdK2ht38Ex864lJ8fvdJKMB0XmA9KEPNzA7gMyjB99HL29rby94QGOm3+9JwfZ3Lb7LzTvr+HYedfhy86389c4x0PNFzYUjwFg+oQz6O1rZc3m37Fk1hWeukdo296XaGhex9J515GbU2BhN2NSxDtXGTMsEWHetAsJhwNs3PE06pGL8I7aV6hrfJvj5n6UXF+R08kxxkJwCeGSA2D+JUuyWTTjw7y+/j521r7ClHGnOJ2kpNq171X21L/G0rnXkecrsnPWOM9jnRBsOgZzEF+Wn2NmXcVrG36O31fKmIqFTicpKXY3vMbOur+zdM51+H2ldr4a4wBrAzLvkO8r4ZgZV/L6pl9SWjiegrwyp5OUUHub3mR77UqOn/1RCnJLreZj0shRj2LgKsmrARlXK/RXENYQOdl5TicloWqbV7Fl7wssnX0tBf5yp5NjzDt4aRxNawMyQ9rfsZOi/Cpys/MzJi/3taylZvcKls76CIV5FRmzX8a4lRVAZkhNbVuoKJ6eMflY37qRjbuf5bgZV1Pkr8yY/TIZyEPnpnVCMENqbt/KvInnZ0Q+NrbVsGH3Uxw77UqK/dUZsU8mQ6m3xtG0+4DMO/T2t9MX6KC0YJzTSTlqTe1bWLv7cY6ZejklBWOdTo4xJob1gjPv0Ny+hYriqWQhrg4HNHdsZ83/b+9eY9s67zuOfx9eJJIiJZGSLFu+W7ZlJ3bsOkaUOG7TtE4btGmyoYO3IUOQN0uyBUE2IBmQIC+yGCiQYQi2vAicYvabDEjWtOuAOOlWD3WcDtlQ+W5l9urItlzfdTdFUhTJ8+yFXF/iu03ykDy/D0DAlojz/MnnOeev53ZO/7/xtfnraY50VPVnEQ/xUDvVHJBcYTDZR1tsYVXX4UjqGHv7/5WVc39IPDKrqj+LeIyHmqoSkFzGsQ5DySMsnfFI1dbhaOo4u/t/yoo5f0CiYU7Vfg6RWqd9QHKZsfQJwnVN1AdjbodyW8bSJ9nd/yHLZ/+Aluh8t8MRuWVemr5QD0guM5j8ktbogqqsv3OZM+w6+i/cPfN7tEVrZwm5eIyH2q2WYctlBpOH6Wr/dtXVX3LiLDv7P2DpjO8yLbq46uIX8SKtgpMLsvkU6ckR4uGZVVV/49khdh77gCXt32ZGbImn/oKUGmPx1B9PmgOSC4ZSR2hpmIvP+N0O5aalJofpOfY+i9seoqPpbrfDEbkjhup5lk8xaA5ILhgYP0xrQ/XM/6QnR+k59j4LWx9kZtPyqolbRKaUcA5IF4NqYq1lMHWYxS1fr4q6y+TO0XP8feYnupnduKIqYha5KR76Q0pDcALAuexp6vwRwsEmt0O5oYlckp7j7zO3+V7mNq9yOxyR4lICKgIPfYm1YGD8MG2R+RVfb9n8OD3HP2BW4z3Ma15d8fGKyLUpAQkAg+nDLEw8WNH1NllI03PiJ8yILWVBvLuiYxW5LVoFVyS6OFSNXGGCZHaQeP3Miq23yUKGnpM/YVqkk87m+ys2TpE75aVVcHocgzCU6ScRnonfV5lTgrnCBDtO/ZSW8FwWJdZijHE7JBEpAq2CEwZSR2gNzavIOss7WXac/hnx+g664t/AWNT7kdrmofZdwiE4Dw1kVjFrLYOZIyxoWl1xdZZ3Jtl55uc0BltZEn8Ig/XUySle5K02rjkgjxufHMCxDvW+cMXUmWPzDGSOcnhsB9FggrsS38JAxcQnF+WdLLsHPibnTBANthANJmgIJogGE4QDTfiMRvnl2ipz0F/KJuCrp7l+OtuObyIR6qA9sohp4QXU+cNljcNay2j2FCdTBzidPkQ02MLs6DJmRu/SnE+FcmyB3QMfU+9vYGHz/aRyI4znhhkZ3894bphsIUUk0HxZUmoIJmgIxCt2vtF1Hhti1hyQx4V9Me5tfZyck2Ugc4Qz6T4ODm+nqa6d9shCpoU7CfkbSlZ+KjfKyfRBTqUPYvDR0bCEB9r/lEigceoNHjshq4W1lt7hrfgJsDy+DmN8xIMzLntPwcmRyo8ynhsmlR/mdOoQqdww6cIY9b4oi5sfYEaky6VPUMEqayS8pDQEJwAETR0dkS46Il0UnBwD2X7OZPr47ejnRIMJpocX0h7uJPz7xHAHJgsZTmcOcSJ9kEx+jOmRxaxIPEpjcNrF3o7aT0U7NPY5qfwo97X+IQZz1frymwCNwVYag62X/dyxBfYN/5JcYUL17HFKQHIFvwkwPdTJ9FAnTnOBoezvOJ3poy/ZQ9jfSHt46ncNwfhNH9OxBc5OHOFk+iDD2RO0hubSGV1Na2jO5XffVrupeMfG93M68yX3t/4RfhO45Trz4SOTHyPWcI/q+yq8tA9IA7FyXT7jpy00j7bQPBz7MCOTJziT6eM3gz8n6AvRHl7A9PBCooGWK+ZqrLWMTp7iRPogZzJ9xIItdESWsDz+CEFfvUufSO7EmfN/iHS3/fC25wmtdRjPDxP7Ss9IzlMCKgIPfYle4cPQUjeLlrpZLG38BqO505zJ9LFzaAs+fLSHO2kPdRL01XEy/VtOZv4Pn/HTEe5iTdsfEw7ELh5M7aPqjEyeonf0V6xOPE7E33jbdZjKjVLnixAwQbUDjyvhIgQPzaR5kAHigXbisXa6og9wLj/AmYnD7B/ZSs5mmRFaxMrm79AYaLvYM1KbqFrj+RF2D3/CPU3raAq03lFdJifP0hhoUXu4GounFnCpByR3zABNgTaaom0sjnZf+Qa1haqWLaTZObKFxdH7aaubfcf1mcwNEgu0ql1clbc2omqXmIhcU96ZZOfox8wMLWFWeElRjnkuP0Qs0FKUY0l1Uw9IRK7KsQX2jP2SxkAbnZFVRTunk/lBYoGErhHX4qHvRRtRReQK1lp6x7dj8HFXw9qi3QR20slQsHnCRHWNuBYloDtnK+zGliJy8w6le0gVRljd+H0MxTufz+UHifkTgMV66EIrV6d9QCJymd9NHOB09jDdTY9PLZUuoqTmf65Pq+CKxENfokitODvZT196F/fFHqOOUNHP42R+iHhghq4P12Qr7rEopaRFCCICwGj+LF+kf82qhu8Q8cVKcg6fyw8zp+4uXR8E0EZUEQFShTF2p/+TZaG1NPlKs0nUsQXSzhhR06Trw/V4KDlrDkjE47JOhp3prSyqX0VbcHbJyhl3xgj7olM3MJWr0xxQkXgoi4tUq7zNsSu9lY5gJ7OCi0p63iYLQ8R82v8jF5VuGba62CIVzbEOe7LbiPniLAgsL/k5m8wPEzNxXRtupMwJ2hjjB3YAJ6y1jxljVgIbgRCQB/7SWvubUpStHpCIB1lr+d/J/8ZYw9K6bszUD0taZtIZZl5QCxBuqPzfz4vAAeD3T5v8O+BvrbW/MMZ87/z/v1mKgnUvOBEP6svtY9yOck/91/GZ0l8GrLUknRFivpt/iKGUnjFmFvB94J8u+bHlYjJqAk6WqnztAxLxmOP5Q5wqHOa+ukcJ2Ft/ountyNo0AHVOiKn7+sjVFfVu2K3GmB2X/P/H1toff+U9/wD8DXDJw7r4K+A/jDF/z1QnZU2xAvqqEg7BaZxXpNIMFE7wZX4v99U9Qj11ZTtPk4UhYiaO8djjBm6ZpZhL1Aettauv9UtjzGPAWWvtTmPMNy/51V8Af22t/ZkxZj2wCVhXrKAuVcJFCGpkIpVkzBmit/A/fM3/EGEbLeu92JLOCDHTrOtCZXkQePz8PE8IaDTG/DPwA6bmhQA+5PLhuaLSHJCIB6Rtkj2Fz7jb3z210bTMknaUmGkue7lVydrivG5YjH3FWjvLWjsP+BPgV9baP2Nqzueh82/7FnCoVB9VQ3AiNW7STrCr8CmdvmW0mRmunJtJO8J8luq6cDPcH6L8c+AfjTEBYAJ4plQFaQhOpIYVbJ7d9tdMZw4zWeDKeVmweSZIE3FiWC1AuAHrygIua+2nwKfn//1fwL3lKFc9IJEa5ViHfXxOlEYWuNj7GLejRIhNjffruiCX0E2ZRGqQtZaD7AJgCaswxrgWyzijxGhyrfyqYr31ME8lIJEadIQDJBnjXh4qy0bT60kyRlQJ6OZ5aPqiZAloq/NhqQ4tItexadMmfvSjXRz6/ADt7e1uh8PatWvZsGEDDz/8sNuhSIVRD0ikhnzyySe89tprfPbZZxWRfBzHYd++faxYscLtUKqH+6vgykYJSKRG9PT08PTTT/PRRx+xaNEit8MB4OjRozQ1NZFIJNwOpTpY66mH9WkjqkiNePfdd3nyySfp7u52O5QL9u7dq96PXJMSkEiNePnll3nvvfc4fvy426FcoAR0G8p0J4RKoAQkUiO6urp4/vnnefHFF2/85jLZs2ePEtAtso5TlFc1UAISqSGvvPIK+/fvZ8uWLW6HAkz1gFauXOl2GFKhlIBEakgoFOKdd97hhRdeIJVKuRrL2NgYAwMDdHZ2uhpHdSnS8JuG4ETEDevWrWPNmjVs2LDB1Tj27dvHsmXL8Pv9rsZRVSxTG1GL8aoCSkAiNeitt95i8+bN9Pb2uhaDFiDIjSgBidSg9vZ23njjDZ599lkclyaklYBuk3WK86oCSkAiNeqZZ57BcRw2b97sSvm9vb0sX77clbKrlWXqUTbFeFUDJSCRGuXz+di4cSOvvvoqZ8+eLWvZ1lq++OILli1bVtZypbooAYnUsBUrVvDUU0/x0ksvlbXc/v5+GhsbicfjZS236lmrITgRqR2vv/4627dvZ9u2bWUrs7e3V72f26QhOBGpGdFolLfffpvnnnuObDZbljL37NmjBCQ3pAQk4gFPPPEES5cu5c033yx5Wel0mo0bN7J+/fqSl1WTPDQEZ+x1dswaY44Cc8sWjYhIdei31s4r9kGNMf8OtBbpcIPW2keLdKySuG4CEhERKRUNwYmIiCuUgERExBVKQCIi4golIBERcYUSkIiIuOL/AXnOGBwUGLuEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cartopy.crs as ccrs\n",
    "\n",
    "sinu = ccrs.Sinusoidal(\n",
    "    central_longitude=ds.crs.longitude_of_central_meridian, \n",
    "    false_northing=ds.crs.false_northing, \n",
    "    false_easting=ds.crs.false_northing, \n",
    "    globe=None)\n",
    "\n",
    "fig = plt.figure(figsize=(8,7))\n",
    "ax = plt.axes(projection=sinu)\n",
    "ds[\"solar_zenith_angle\"][0].plot(ax=ax)\n",
    "ax.coastlines()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pixel index\n",
    "\n",
    "Make an array of index values for the permuted xy arrays. The evapotranspiration model takes inputs that sequence in column-major order:\n",
    "```\n",
    " 1  24.249  -79.398  0.000 0.000 0.000 ... \n",
    " 2  24.266  -79.398  0.000 0.000 0.000 ...\n",
    " 3  24.284  -79.398  0.000 0.000 0.000 ...\n",
    " 4  24.302  -79.398  0.000 0.000 0.000 ...  \n",
    "```\n",
    "\n",
    "But `xarray`'s mappings are more convenient. So, reshape the index array using the shape of longitude array in the results dataset:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Potentially useful code)\n",
    "\n",
    "Coordinate array for `xyindex`:\n",
    "```python\n",
    "xyindex = np.array(list(range(ds.x.size*ds.y.size)))    # seq 1 to npixels\n",
    "xyindex = xyindex.reshape(ds.lon.shape)                 # reshape with lon\n",
    "\n",
    "ds.coords[\"xyindex\"] = xr.DataArray(\n",
    "    data=xyindex,\n",
    "    coords=[ds.y, ds.x],\n",
    "    dims=[\"y\", \"x\"])\n",
    "\n",
    "ds[\"xyindex\"].plot(x=\"x\", y=\"y\", figsize=(6,5))\n",
    "```\n",
    "\n",
    "Dask piping:\n",
    "```python\n",
    "df = dd.from_pandas(pd.DataFrame(load_boston().data),npartitions=10)\n",
    "\n",
    "def operation(df):\n",
    "   df['new'] = df[0]\n",
    "   return(df[['new']])\n",
    "\n",
    "df.pipe(operation).to_csv('boston*.csv')\n",
    "```\n",
    "\n",
    "Simple progress widget:\n",
    "```python\n",
    "import time, sys\n",
    "from IPython.display import clear_output\n",
    "\n",
    "def update_progress(progress, nbar=20):\n",
    "    \"\"\"Simple ASCII progress bar for Jupyter environment.\"\"\"\n",
    "    \n",
    "    if isinstance(progress, int): \n",
    "        progress = float(progress)\n",
    "    if not isinstance(progress, float): \n",
    "        progress = 0\n",
    "    if progress < 0: \n",
    "        progress = 0\n",
    "    if progress >= 1: \n",
    "        progress = 1\n",
    "    block = int(round(nbar*progress))\n",
    "    \n",
    "    clear_output(wait = True)\n",
    "    prog = \"Progress: [{0}] {1:.1f}%\"\n",
    "    print(prog.format(\"#\"*block+\"-\"*(nbar-block), progress*100))\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write text outputs\n",
    "\n",
    "Now write each valid pixel's time series to output text files. We'll concatenate them column-wise in shell script. This could take a while so import the progress bar from dask. \n",
    "\n",
    "The for loop below is iterating over the output netCDF files and reformatting as text files where columns are timesteps and rows are pixels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask.diagnostics import ProgressBar\n",
    "\n",
    "\n",
    "for f in (bsa_result + wsa_result + alb_result):\n",
    "    print(f)\n",
    "    \n",
    "    # make a numpy array to serve as pixel index\n",
    "    xyi = np.array(list(range(ds.x.size*ds.y.size)))\n",
    "\n",
    "    # reshape band data to match USGS model requirement\n",
    "    b1 = ds[\"BRDF_Albedo_Parameters_Band1\"]\n",
    "\n",
    "    # reshape continue\n",
    "    b1stack = b1.stack(xy=(\"x\",\"y\"))\n",
    "    b1flat = b1stack.reset_index([\"xy\"])\n",
    "    b1flat.coords[\"xy\"] = xr.DataArray(data=xyi, dims=[\"xy\"])\n",
    "    b1tflat = b1flat.T\n",
    "\n",
    "    # convert to dask dataframe\n",
    "    b1srt_ddf = b1tflat.to_dataset().to_dask_dataframe()\n",
    "    \n",
    "    out = f[:-3]+\"-*.dat\"\n",
    "    with ProgressBar():\n",
    "        \n",
    "        b1srt_ddf = b1srt_ddf[b1srt_ddf.land_mask]\n",
    "        b1srt_ddf = b1srt_ddf.categorize(columns=['time'])\n",
    "        \n",
    "        b1srt_ddf.pivot_table(\n",
    "            values='BRDF_Albedo_Parameters_Band1', \n",
    "            index='xy', \n",
    "            columns='time'\n",
    "        ).to_csv(out, sep=\"\\t\")\n",
    "\n",
    "    #timeseries = None         # discard objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Albedo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t = ds[\"BRDF_Albedo_Parameters_Band1\"].groupby(\"time.month\").mean(\"time\")\n",
    "t"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}