gp1_arm_site_data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:43:27.559694Z",
     "start_time": "2018-08-17T06:43:27.555107Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "\n",
    "import glob\n",
    "import xarray\n",
    "import matplotlib.pyplot as mp\n",
    "import matplotlib.cm as mcm\n",
    "import datetime\n",
    "\n",
    "import cmocean\n",
    "import netCDF4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:11:29.310677Z",
     "start_time": "2018-08-17T06:11:29.172210Z"
    }
   },
   "outputs": [],
   "source": [
    "file_list = numpy.array(sorted(glob.glob('/beegfs/DATA/pritchard/blangenb/ARM_DATA/*c1*.cdf')))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:15:19.332405Z",
     "start_time": "2018-08-17T06:15:19.211485Z"
    }
   },
   "outputs": [],
   "source": [
    "file_datetimes_strings = [i.split('/')[-1].split('.')[2:4] for i in file_list]\n",
    "\n",
    "file_datetimes = numpy.array([datetime.datetime.strptime(i[0]+i[1], '%Y%m%d%H%M%S') for i in file_datetimes_strings])\n",
    "file_datetimes_local = file_datetimes-datetime.timedelta(hours=4)\n",
    "\n",
    "fmaaso_indices = numpy.array([i.month in [2,3,4,8,9,10] for i in file_datetimes_local])\n",
    "file_list_fmaaso = file_list[fmaaso_indices]\n",
    "file_datetimes_local_fmaaso = file_datetimes_local[fmaaso_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:15:38.664382Z",
     "start_time": "2018-08-17T06:15:38.651804Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1,  2,  7,  8, 10, 11, 13, 14, 19, 20])"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unique_hours = numpy.unique([t.hour for t in file_datetimes_local_fmaaso])\n",
    "numpy.unique([t.hour for t in file_datetimes_local_fmaaso])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "group by unique hours above, take means, then plot PBL heights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:25:26.564692Z",
     "start_time": "2018-08-17T06:23:17.723903Z"
    }
   },
   "outputs": [],
   "source": [
    "pbl_height_h_array_list = []\n",
    "pbl_height_ll_array_list = []\n",
    "pbl_height_br50_array_list = []\n",
    "pbl_height_br25_array_list = []\n",
    "\n",
    "for h in range(unique_hours.__len__()):\n",
    "    specific_hour = unique_hours[h]\n",
    "    time_indices = numpy.array([t.hour==specific_hour for t in file_datetimes_local_fmaaso])\n",
    "    pblh_array_h = numpy.zeros(numpy.sum(time_indices))\n",
    "    pblh_array_ll = numpy.zeros(numpy.sum(time_indices))\n",
    "    pblh_array_br50 = numpy.zeros(numpy.sum(time_indices))\n",
    "    pblh_array_br25 = numpy.zeros(numpy.sum(time_indices))\n",
    "    \n",
    "    files_to_open = file_list_fmaaso[time_indices]\n",
    "    \n",
    "    for i in range(files_to_open.__len__()):\n",
    "        ds = xarray.open_dataset(files_to_open[i])\n",
    "        pblh_array_h[i] = ds['pbl_height_heffter'].values\n",
    "        pblh_array_ll[i] = ds['pbl_height_liu_liang'].values\n",
    "        pblh_array_br50[i] = ds['pbl_height_bulk_richardson_pt5'].values\n",
    "        pblh_array_br25[i] = ds['pbl_height_bulk_richardson_pt25'].values\n",
    "        ds.close()\n",
    "    pbl_height_h_array_list.append(pblh_array_h)\n",
    "    pbl_height_ll_array_list.append(pblh_array_ll)\n",
    "    pbl_height_br25_array_list.append(pblh_array_br25)\n",
    "    pbl_height_br50_array_list.append(pblh_array_br50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:41:26.048495Z",
     "start_time": "2018-08-17T06:41:26.031539Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/export/home/blangenb/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: Mean of empty slice\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/export/home/blangenb/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: Mean of empty slice\n",
      "  after removing the cwd from sys.path.\n"
     ]
    }
   ],
   "source": [
    "pblh_h_means = numpy.array([numpy.nanmean(height_vals) for height_vals in pbl_height_h_array_list])\n",
    "pblh_ll_means = numpy.array([numpy.nanmean(height_vals) for height_vals in pbl_height_ll_array_list])\n",
    "pblh_br25_means = numpy.array([numpy.nanmean(height_vals) for height_vals in pbl_height_br25_array_list])\n",
    "pblh_br50_means = numpy.array([numpy.nanmean(height_vals) for height_vals in pbl_height_br50_array_list])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-08-17T06:59:11.624796Z",
     "start_time": "2018-08-17T06:59:11.423370Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,1.02,'PBLH at Brazil ARM site (2014-2016)')"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAEbCAYAAAC/TMN/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXV4k1cfhu8TqQGFAqXFrVgpDhsMGO4M9+IMd9twGWy4bTiDYcWdYUOGw3At7lpKkdJSS3K+PxL4SqmkpWna8t7XlYv2fY88KXmfHP0dIaVEQUFBIbmhsrYABQUFBUugmJuCgkKyRDE3BQWFZIlibgoKCskSxdwUFBSSJYq5KSgoJEuSvLkJIaQQQmd66YUQj4QQ3Uz3Kka4HyyEOCGEKBju/vUoyl0qhBgTRX05YqnxvhCidCTXo9UXHwghFgshfjL9LIUQrtGkHWdKkz/C9YOmv+0HnUFCiP+EECVM99ub8g2OkE8thHghhLgfC73ZhRABpp+j/P8xo5z6QoifTT/XEkJcE0K8F0KcDf9/IYT4SQjx0vQaG0k5nYQQ86Ooo6AQIiSqv6kQQiWEmC6E8DO9VgghUpnupRVC7BRCBAohrgohykWS/0jEz40QorXp8/RWCLFZCJFGCOEkhNgcu79Q8ifJm5uJLFJKDaABWgBTwxnEDSmlxnQ/NfAPsMZKOiPDLH1CCHVcCpdSdpJSTjYzeSvgEsa/YUQ8w+lMZ9K5Itx9f6BxhDzfA7ax1PtASpkyNnkiYvpbjQQWCCEcgbXAJIy6ZwCbhBB2QogqQB+gPFAC8BRCNDCV4SaEGAD8Fk0diwCbaKS0B+oD3wG5gFTAL6Z7s4AXQGZgPLBBCGFjKruGyVA/MTwhRFFgOtASyGKq+ycp5WvgvhCinll/oK+E5GJuAEgjx4DrQIFI7ocAXkDu+K5bCJFHCHHU1Dp4KITobbq+A8gOHBVCfBuD/k/0mVpEfwshtgFbTdeGCyGeCyEChBD7hBBZTNdvR2jBStP1pUKIIWboLwsEYzSFljHofI/RgHOFu3wKyPlBj4nGwI5o6uwjhHhmatX8aTKcHKYWrAuwH8gnhLhiSl9UCHHS9N73CiGyR1F0A+CilPINUAZ4IKVcKqUMklKuBAyAO9AamC+lvCalfADM4//GnseU5kkUdQwATkb13kxUA+ZKKW9IKd8Cs4EqQghboAkwQkr5Rkq5GvAFKpvylcP4RR0QobyOJr0npJTvgA7ActO9PzH+3ymYSFbmZuoGlAPyAhcjuW8PdML40MQ3vwEHgDRAQ2C6ECK1lLIO8AAoJ6X8L7oCotBXC1gC1DO1RnsBpYEMGB+I/gBSSjdTq8oO+A/4PZb622B8UHYDGYQQxWPQ2cb0fj+gx2jAjUxpBEaTibS7JITICYwGimM087wYH15M78cHqIKxZeshhEiB0Sh/A1yB08DSKCQ2DqftFEYj+VBvNiA98BwoDJwPl+8q4Gaqf5eU8kdgWyTa82I0lpjMZGQEjSWBp6b3GiClfBxF3SNNdftFKK8EoBJCnBNCvMTYGn1kynMVyCqEyBSDpq8GjbUFxBP3jc8SAuOHp7+U8pYQIjOQVwgRbEr3oQsRbcskHCPMafWYGAk8BFIA9oAEnIC3MeSLSd9pKeUWACHEI4xdqMdANiAUY1crPGMAB+AnM3Vj6g41BopJKUOFEFtMGs6FS7ZcCLE0nE49UDRCURuBYRiNtQzGVs+jqKo16axtylfTdD1DFOnrAheklNtMmscAb4UQTqZuWXhKA5MBTPdem/JUwGg2s6WUT01d1vD/P+8Axyjqx1SGClgM9JVSBpo+d5EipbxtyqMFhgCDgOqmOiJ+LmKsG6MpN8H4t3gFrMfYzR1oun8ZYxd4QwzlfBUkl5ZbDimlnZTSVkqZU0q5ONy9m6Z7doAWY0voT1PLISbGf8gbroyo+Abjt+8BoB3Gh9ccYtL3JlxaG2AmcAeYj9HgPmJ6ePsALU1dXHOpg9EkL5haBE2B5uLTJ7dtOJ12wB8mDeE5ALibupSNgU1RVSilvGuqpylGA9yOsfseFdmBGqYuazDGMT4N4BJJWhdMhgYghHAQQiwy6RkvpRxkuvUG4xfRBxzC54uCnsAdKeXeiDeEEPvDDQ1UMF0rClzA+DcuY2q9R6zX3LoBFkgpb0spX2Ecf6sc7t4rQGm5mUgu5mYWUkq9lHIPcAPwiK9yhRAOGAeXa0spiwF9icPf1gx9/TF+w+eQUlbH2P38oMEJ4wD/ACnltVhW3Rpji6uo6ZUf48P32QyeSWcoxtZWlgjXwzB2HRti7J5GaW5CCDfgpZSyJkYzugGMi0bjC2BjOIP9oO9mJGk1mL5cTAP/u4CMQP4IX3zX+fTv7E4kwxkRqAi0DGeyYOw5/CClrPJh0kVKeUgIUQg4hPFLoIyU0vtDesBJCJE2lnXfx/gF+AEVEBTud4lxPFGBr8zchBBaIURdjGMeH7pcGiFElggvbTTFRIba9HIQQjjz/xm2DzOFOoxdirjoC48G44fbzjS26AloTS2sRcB/Uso/YyNcCJEGY9dwhZTy8YcXRpOKrvuuJ/LPz0ZgKMYxpVvR5M8BrBZC5MH4UNpiHEMMjw5wNP1/7AK+F0KUFkLYAT9jHKyP7GG+z/+761UwGlsjKWXE8r2AHkKIzEKIfEA/YFU0mpFSNjb1EMK35HNIKbdHknwoMENK+YcMF35HShmIcSxvjBDCXgjRCeNMarRjshgncXoIIfKbjHEAxpngDzgBPjGU8dXwNZhbvg9dBSAQmIix2/ZhFiw3xm5R+FeR2FRgmrkaBOzDOEB9A+MHcaUpyd8Ylx+UiIO+8MzCuHTADxgO9ADqAW0xdgMbhesW6aKZTQxPM4wzixHr2w40FUJENS4bAGQ0DdCH5x+MD1mUrTYAKeU+jGNGRzE+kBkwTjCE5wrGlslp0wRDa4yzgq8wdt8jW7ICcBwoZPq5MMaB+vcR/jYVpJS7MU6iXAaOAH9IKQ9GpzuWFAZGRqj3tuleX4wt5FcYJ4maSCn1MZS3FOOM7n7AG+Okyuxw9z2IeQb3q0Eo8dwUkhtCiCZAfSllG2trSShMLc+VUspS1taSWPgaWm4KXx9bME5sOFlbSALSGWOrX8GEYm4KyQ4ppY7/d9uTPUKI1EBBKeVGa2tJTCSIuQkheplWovsL44r7jJGkqSGEuCWMe+3WCdMePAWFuCCl3C2l/NXaOhICKeVbKWUta+tIbFjc3IQQHsBYjKvVM2FcyjA5QhonYDUwGOPyAkH0ywIUFBQUoiUhWm7fA39LKf+TUgZgnH6PGPWiHnBOSrnFtKJ8IlHPhCkoKCjESEKY2yKMg50IIdJhXJt1PEKaiHv8rgAuQogvig6hoKDw9WLxvaWmVesIIToDCzFGnqgWIZkjxs3lH/KEmNZ9ORIhMoIQogvQBSBFihQl8uf/JPSYgoJCMuTs2bMvpZTOscmToOvcTNEkOgOjgIzhjG8KECqlHG763QYIARyklEFRlVeyZEl55swZywtXUFCwKkKIs1LKkrHJkxATCsOFEK0BTEY1D0jLpxEQIu7xKwDcis7YFBQUFKIjIcbcfIEhQohcpj2Bg4DLUsrwsaq2AOWFEOVNM6e/EcMePwUFBYXoSAhzWwLsxTiJ4INx9rSxEKLth312JqNrh3Gf3yOMESAmJIA2BQWFZEqS3luqjLkpKHwdJMoxNwUFBQVroJibgoJCskQxNwUFhWSJYm4KCgrJEsXcFBQUkiWKuSkoKCRLFHNTUFBIlijmpqCgkCxRzE1BQSFZopibgoJCskQxNwUFhWSJYm4KClbCYDCwfflJXjx5Y20pyRLF3BQUrMROr9NM6rOW4W3+QqeL6bB5hdiimJuCghV4/TKAuaO2kSFLGq6de8iKafusLSnZoZibgoIVmDNiK+8DQpi2oSvVmhZn6ZR/uH7+obVlJSsUc1NQSGDOHrrF7jVnaNW3Mjnzu9J/SmPSZkjFuK5ehASFWlteskExNwWFBCQkOIypA9eTOWd62g2sCoBjGgeGzWnJg5svWPDLDisrTD4o5qagkICsnLGfR7d9GTS9Cbb2Nh+vl6qUj8ady7Fu3mHOHb5lRYXJB8XcFBQSiIe3XrByxj6qNS1OqUr5PrvffewPZHVzZnyPVQS8VQ5++1IUc1NQSACklEzpvx5bext6/9og0jR2DjaMXOCJ3zN/Zg3ZnMAKkx+KuSkoJAC7V5/m/NHbdB9bl7QZUkWZzr1EdtoMqMqu1ac5tP1SAipMfijmpqBgYd74BTB7xDYKfZuDH9qWjjF9+5+qk7dIFqb0W8erF+8SQGHyRDE3BQULM3fUdgL8gxg8oxkqVcyPnEarZuQCT94HhDCpz1qS8vGb1kQxNwUFC3L+6G12ep2iZe9K5HLPaHa+nPld6TqqDsd2X2XHylMWVJh8UcxNQcFChIbomDpgPRmzp6X94Oqxzt+0+/cUK+fGrKGbeXrfzwIKkzeKuSkoWIhVvx/gwc0XDJjaBDsHm5gzREClUjF8bktUKsGv3Veh1xssoDL5kiDmJoRoJ4S4L4QIEEIcEkJ8tsjHdD043GtJQmhTULAEj+74snzqXio3LEqZagXiXI5rtrT0m9iQiyfusm7uoXhUmPyxuLkJIfIAc4DugDOwG1gTSdJcQGoppZ3p1dHS2hQULIGUkqkD1qO11dBnQuRr2mJDzZalKF+nEAvH7eCu97N4UPh1kBAtt8rAv1LKXVLKIGAqUEQIkfZDAiGEHaCTUoYkgB4FBYuyd/05zh66RbfRdUnvmvqLyxNC8NPMpqRMbc+4rl6EheriQWXyJyHMbQvQO9zvJYD3gH+4azkAGyHEOSHEKyHEZiGEawJoU1CIV/xfB/L7sC24l8xO/Q5l4q1cJ+dU/DSzGbcuP2HJpD3xVm5yxuLmJqX0kVLeBxBCNAa2A0OklOG/fpwAb6A1RqN7C6ywtDYFhfhm3pi/eff6PYNnNDVrTVtsKF+nELU9v8Frxn6unLofr2UnRxJqQiG9EGIrMAvoJKWcHf6+lPKElLKalNJbSukPDAIqCyHsIymrixDijBDijK+vb0LIV1Awi4sn7rJ92Uma9ahAnkKZLVJH3wkNyZA5DeO7eREUqIziREdCTCikAI4CfkB+KeW2SNLUFkJUDXdJC+iAzyL3SSkXSilLSilLOjs7W0q2gkKsCAvVMaX/OlyyONFxSA2L1ZPC0Y7h81rx5J4fc0Ztt1g9yYGEaLl5As+klB2llAFRpEkLzBdC5BdCOAKTgA1SSuXUDIUkwZrZB7l/3YcBUxtjn8LWonUVK+dG854V2LL4GCf3XbNoXUmZhDC3wkAFIYQuwqudEOI2gJRyJbAaOAzcAdRArwTQpqDwxTy595K/Jv9DxXqFKVuzYILU2XlEbXIWcGVCrzX4vw5MkDqTGgkxodBLSqmSUmoivJZJKd3CpRsppcwgpXSWUnpKKV9bWpuCwpcipWTaoI1otCr6TmyYYPXa2mkZOd+TNy8DmDZoY4LVm5RQtl8pKHwB+zed59T+63QZURvnTGkStO68RbLQcUhN9m88z76N5xK07qSAYm4KCnHk3Zsgfh+6hfzFstLwx3JW0eDZrzIFS2Vn2sCN+D5VTq4Pj2JuCgpxZMEvf/PmZQCDZzRFrbbOo6TRqBkx35OwUB0Teiux38KjmJuCQhy4cuo+W5Ycp2m378lXNKtVtWTN7UyPX37g1P7rbFly3KpaEhOKuSkoxBJdmJ4p/daRIXMaOg2rZW05ADTsVJZvquRnzshtPLqjLG4HxdwUFGLN2rkHueP9jP6TG+GQ0rJr2sxFCMHQP5qjtdEwvpsXOp2yRFQxNwWFWPDswSuWTNxD+doelK9TyNpyPsE5UxoGTmvM1dMP8Jp5wNpyrI5ibgoKZiKlZPrgDahUgn6TG1lbTqRUbVycKo2KsWTibm5efGxtOVZFMTcFBTM5uO0iJ/65RufhtXHJ4mRtOVEycFpj0qRPybhuXoQEh1lbjtVQzE1BIQYC3gaxcuZ+JvddT97CmWnUxTpr2szF0SkFQ2e34N615ywav9PacqyGxtoCFBQSK34+/qyff5jNi48R6B/MN5Xz0W9yIzQatbWlxUjpqgVo0PE71s45RNmaBSlWzi3mTMkMkZQX/ZUsWVKeOXPG2jIUkhlP7r1k9R//stPrFGGheio1KIJn38pWX88WW4ICQ2hfbip6nZ5lx34ihaOdtSXFGSHEWSllydjkUVpuCgombl16wspZ+/l38wXUGhW1Wn1Dy96VyJo7acYNtE9hy8gFnvSo+Tu/D9vC0NktrC0pQVHMTeGrRkrJhWN3WDlzP//tu45DKlta9K5Es+7fx8vhLtbG45scePavwopp+yhXq2CiW75iSRRzU/gqMRgMHNt1lZUz93P19AOcnFPSdVQdGnQsS6o0n0W3T9J0/LkGJ/+5xqS+6/D4JgdOzqmsLSlBUGZLFb4qdGF6dq06RdsykxnquYRXL94xcFoTNlwaSZsBVZOdsQFobTSMXOBJoH8Qk/ut/2o21ystN4WvgqDAELYvP8maOQd58fgNud0zMmpRayo3LJokZj+/lFzuGek8sg5zR25j9+rT1Gr1jbUlWRzF3BSSNW9fBbJx4RE2LjzK21eBFCmTi8HTm1K6WgGEENaWl6A071GB47uuMnPIZoqVc8M1W9qYMyVhlKUgCskSn8evWTvnINuWnST4fShlaxbEs19lCpfOZW1pVuXpfT/alZtCgWJZmbm1e7yfrWoplKUgCl8992/4sGrWAfasOwMSqjYpjmffyuRyz2htaYmCTDnS0fe3Bkzss5YN84/QrEcFa0uyGIq5KSQLrp55wMoZ+zmy4zK29loadipLi54Vk33XKy7UafMtR3ZeYf7YvylVOR8587taW5JFULqlCkkWKSWnDtxg5Yz9nD96m1RpHGjSpRyNupTHKX1Ka8tL1Lx68Y62ZSaRIYsTC/f1Q6NN3JMqSrdU4atArzdwcOtFVs7Yz63LT3DOlJrev9bnh3ZlEk3wyMRO2gypGDyzGcPb/MXSyf/w4/DEEVE4PlHMTSHJEBIcxq7Vp1n9+788ufeSbHkyMOSPFtRoXgKtjfJRji0VfihMzRYlWTF9H2VquFOwZHZrS4pXlG6pQqz48HlJyGUUAW+D2PLXcdbNPcSrF+8oUDwbrftVoVwdD6udOpVcCHgbRNuyk7G10/LX4UHYOdhYW1KkKN1SBYsSFqpjQKMFXL/wiFwFXMnlnpHc7hnJ5Z6J3AUzkjptinitL2LIoVKV8jJ6URuKf+/21a1RsxQpU9szfE4r+tafy7zR2+k/pbG1JcUbirkpmM3cUds5f/Q21ZoWx/fpWw5tu8T2ZSc/3k/n6khu94zkLGAyvYIZyZnPBVv72LUGIoYcqli/MJ79qpA/iYUcSiqUqJCHZt2/Z928w5Sr7UGpSvmsLSleSJBuqRCiHTAWSA+cBbpIKW9ESFMDmA1kAnYAnaSU76IrV+mWJhz7N51ndMflNOv+PX0mNASMXVS/5/7c8X7GHe9n3DP9e/+GD6Gm8NYqlSBzrvSmFl5GcrtnIlfBjGTKke6zLmXEkEM1W5aiVZ/KSTbkUFIiJCiUjhWm8T4ghGXHf8IxjYO1JX1CXLqlFjc3IUQe4DzQFDgI9AOaSSmLhUvjBNwBOgKHgIXAEyllv+jKVswtYbh/w4fOlafjVigzf2zvGeOyAb3ewJO7L7nj/ZS7V42Gd9f7GU/u+X0cs7O115Izvyu5CmQkR34Xzh6+xX/7rmOf0paGHb+jWfcKpM+Y9EMOJSWun39I12qzqNKoGKMWtra2nE9IrObWFagrpfzB9LsWCAHSSylfma61A9pIKauafi8B7JBSRru6UDE3y/P+XTCdq8zA//V7lhwaiHOmNHEuKygwhPs3fLhz9Rl3vZ8aW3vXnvPqxTvSpE9Js+7f06BT2UTXavia+GvSHhZP2M0vS9tRuUFRa8v5SGKdUNgC7An3ewngPeAf7lphjK27D1wBXIQQKaWUAZaXqBAZUkom9V3Ho9u+zNjc/YuMDYyRYQsUz0aB4tk+uf7GLwCHlHbY2CpDwNamzYCqHN/jzdT+6ylcOmeSDthp8Xl0KaWPlPI+gBCiMbAdGCKl1IVL5gi8DZcnBNCZrn+CEKKLEOKMEOKMr6+vRbV/7WxceIT9m87TeWRtSlTIY7F60qRLqRhbIkGjVTNygSfBQWFM7L02Scd+S5BFQkKI9EKIrcAsjBMFsyMkeQPYh0tvg7FV+TpiWVLKhVLKklLKks7OykCzpbhy6j5/DN9KuVoeePatbG05CglItjwZ6DG2Lif3XmPbshPWlhNnLG5uQogUwFHAD8gvpdwWSbLrgEe43wsAt6SUQZbWp/A5r33fMbL9UlyyODF8XqskExZHIf5o1LkcJSvmZfbwrTy599LacuJEQnxqPYFnUsqO0YyfbQHKCyHKm2ZOfwNWJYA2hQjo9QbGdFrB21fv+XV5h2QZdlshZlQqFUNnt0CtVjG+mxd6vcHakmJNQphbYaCCEEIX4dVOCHEbQErpB7QDlgOPgBfAhATQphCBP3/bxdnDtxg0rQl5Cme2thwFK+KSxYn+Uxpz+b/7rP7jX2vLiTUJMaHQS0qpklJqIryWSSndwqXbLqXMKaVMKaXsYJpUUEhAju66wopp+/ihXWlqeyb/GPsKMVO9WQkq1i/Cn7/u4tblJ9aWEyuUwRQFAJ7cf8n4rl7kLZKFfpMaWVuOQiJBCMGg6U1wdHJgfFcvQkN0MWdKJCjmpkBIUCgj2ixFqFSMX9YeWzuttSUpJCLSpEvJkD+ac8f7GX/+tsvacsxGMTcFpg/exK3LTxi10JNMOdJZW45CIuS7GgX5oV1pVv/+LxdP3LW2HLNQzO0r5+/lJ9mx8j/aDapGmeru1pajkIjp/WsDMmZPy/huXrx/F2xtOTGimNtXzI0Lj5g+eCMlK+al49Ca1pajkMhxSGnL8HmteP7wNX+M2GptOTGimNtXiv+b94xou5Q06VMy+s82SkRbBbMoUiYXrfpWYvuykxzbfdXacqJF+UR/hRgMBsZ39cL32Vt+WdpOOSlKIVZ0GlqL3O4ZmdRnLW/8Em9cC8XcvkJWztjP8T3e9P6tAR6lclhbjkISw8ZWw4gFnvi/fs/U/usT7eZ6xdy+Mk4fvMGfv+6iWtPiNPqxrLXlKCRR8hTKzI/Da3Fw2yX+WXfW2nIiRTG3r4gXT94wttMKsudz4aeZzZRDVhS+iJa9K1GodE5mDN6Iz+PPAvhYHcXcvhLCQnWMbL+UkBAdvy7vgH0K5fBihS9DrVYxYl4r9HoDv/VcjcGQuDbXK+b2lTBnxDaunn7AsDktyJYng7XlKCQTMudMT69f63P20C02LTpqbTmfoJjbV8DeDefYsPAIzXtWpFL9xBMXXyF5UK9dGcpUL8Dc0X/z4KaPteV8RDG3ZM7da8+Y1GcthcvkpPuYutaWo5AMEULw8+/NsXewYVxXL3RhemtLAhRzS9a8fxfMiLZLcUhpyy9L2sV4JJ9C1Oi9LxAyfiC6/w5ZW0qiJL1ragbNaML1849YPm2vteUAirklW6SUTOi9hsd3fBmzuK1yBmgckaEhhG1YStiiaUi/F+hPHrS2pERLpfpFqd6sBMum7OX6+YfWlqOYW3Jl/bzD/LvlIl1H16F4ebeYMyh8huHRfUKnjUJ/bD/qirVQV62HfHAH+e5tzJm/UvpNbkRal1SM6+pFSFCoVbUo5pYMuXjiLnNGbaN8nUK06qOcXBVbpMGAbt92QmeOQYYEoe0+BG39VqiLfgNSove+YG2JiRbHNA4Mm9OSBzdfMH/sDqtqUcwtmfHqxTtGdVhGxmxpGT63pbJQN5YYXvkSOuc3dDvWoSpUAtvBv6HOWxAAkSkbpEmH4er5GEr5uilVKR+NO5dj/fzDnD10y2o6FHNLRuh0esZ0Wk7A2yDGL+9AytTKyVXmIqVEf+YYoVOGI588QNuqK9p2vRAp/h9UQAiB2qMYhhuXkaHW7XIldrqP/YGsbs782nMV795Y54ROxdySEX/+uotzR24zaHpT3DwyWVtOkkG+DyRsxVzCvOYjMmbFZvCvqEuVi7TVqypYDEJDMdzytoLSpIOdgw0jF3ji98yfWUM2WUWDYm7JhCM7LrNyxn7qdyhDrZalrC0nyaC/5U3IlGEYLp5GU7spNr2Go0oX9Q4OlVsBsLXDcPVcAqpMmriXyE7bgVXZveYMh7ZdSvD6FXNLBjy+68uv3VeRv1hW+kxoaNG6dIf/QXdoT6INc2MuUhdG2LbVhM2biNDaYNN3FJpq9RCq6B8JodGiyl8IvfcFZCLbS5kYaTe4OvmKZmFyv3X4+fgnaN2KuSVxgt+HMrzNUlRqFeOWWvbkKv2Vc+g2r0C3ZSW63dbpasQHhmePCZ0xBv2/O1GXqYTNwHGosuUyO7+6YHF4+xr5+L7lRCYTNFo1Ixe0JigwhEl91ibol6LG3IRCCHcgExEMUUr5T3yLUjAPKSXTBm7grvczJq/rTMbsaS1X1+uXhK1eiMiSE1XmrOj/2YJQqdDUsGxLMT6RBgP6o3vRbV8LdnZofxyAumCxWJejci8CQqC/ej5Wpvi1kiOfC91G1+X3YVvYseI/6rYtnSD1mmVuQohpQHsg4pleElDMzUpsX36SXatP0+Gn6pSpVsBi9Ui9jtDlc0CvR9uuJyKtM0hpbL2pVGiq1bdY3fGFfPuasNWLMNy4jMq9KNoWPyJSxW3XhkiRCpEzr3FJSK3G8aw0edKkW3mO7rrCrGFbKP59ngQ5QtLcbmk7wF1KWSrC6xtzKxJCHBFCRGrZQohDQojgcK8l5pb7tXL9/ENmDN7IN5Xz0f7nGhatS7dzA/L+bbTNf0SV3sXYYmv+I6qSZdHt3IBu33aL1v+l6C+dNk4a3L2Bpkl7tD8OiJOx+fsGfNwUrvYohnzyAPn6ZXzLTZaoVCqGzW2FSiX4tfsq9HrLj1eaa263AZe4VCCEqCGEmA+UiyZZLiC1lNLO9OoYl7q+FvxfBzKi7VKcMqRi1KLWFj25Sn/tIvoDO1B/Vxl1sW8/XhcqFdqWXVDj7vtOAAAgAElEQVSV+A7djnXoDlh3NXpkyOAgwtYsIuyv3xFOztgMGo+mbJU4LWy+ceI+P387k6lNlhH4OghVweIA6JUFvWbjmtWJfhMbcvHEXdbOtXwAAnPH3H4Bjggh9gKfTHmYYUTlTPVEekyOEMIO0EkpQ8zU8lVjMBj4pYsXL5/7M3d3b9Kks9zJVfLNK8K8FiAyZUVT3/Oz+x8MLsxgQLd9jbGLWrGWxfTEBsP9W4StnI985Yu6Wj00NRoi1GYPMX/C42s+LOy0iiLOcM/7MZMaLqGvV2tSOrtiuHIeylWLZ/XJl5otS3Fk5xUWjdvBt5Xzkbug5dZjmvuVPwXYC5wFbkR4RYuUcqSU8kfAL4okOQAbIcQ5IcQrIcRmIYSrmbq+OpZP3cfJvdfoO7Eh7iWyW6weqdcTumIuhIWibdcbYWMTaTqhVqP17Iaq6Lfotq5Cd2i3xTSZg9TrCNu1kdDfx4E0YNNrONraTeNsbC8fv2Ftm8VUSBVIDvGeOhVdeePzjon1FhOYMR+G297IYOuswE+KCCEYPLMZKVPbM66rF2GhOovVZa65ZQbaSCknSCknhX/FgwYnwBtojdHo3gIrokoshOgihDgjhDjj6+sbD9UnHU4duM7iCbup0awEDTp+Z9G6dHs2I+/eQNu0A6oMGaNNK9RqtK27oSpcCt0WL3SHrTPHZPB9Tugf49H/swVVibLYDPoVVa58cS7P9+pj9jadRUFNAKndXMhQvgDBNx8zcJWxFbts0QPQ6zFcvxxfb+GrwCl9Sn6e1ZzbV56yZOIei9UTm5bbSCFEqvgWIKU8IaWsJqX0llL6A4OAykKISDdGSikXSilLSilLOjs7x7ecRMvzR68Z++NKchZwZdCMphbdEK+/cQX9vm2ov/0edclPj/+TUka6VkmoNWjb9kBVqKRxLdzRfRbTFxEpJbqTBwmdOgLp+xxt217YeHZF2DvEqTyDTs/NxQc43WkeKfUhpGtclu+X9yJ7i7IYQnTYvnrDkG2deGPjQmCoCt89/8bzO0r+lKvtQZ3W3+I1cz+X/7tnkTrMNbchptcbIUSQEOK96d/ALxUghKgthKga7pIW0AHKzmQToSHGk6vCwnSMX97eoidXSf83xj2WGTKhadjmk3v6kDBO9VjEqR5/ooskVpfR4Hqi8iiObuMydMcPWEznR70B/oQtmYlu7WJU2d2MUTzCTXzEljdXHnK83WzuLtrHs/cCl/6NKDW4DkKtIm3RHNg4peD5v1dIlzk1P23uxIMQZzQPvDmw5GQ8vquvgz6/NcAlixPju63ifUD8D7mba27DgZwYZzXzAwWAksCJeNCQFpgvhMgvhHAEJgEbpJSJIxB7ImD2iK1cO/uQ4XNakc3NcidXSYOBsJXzIDjIGBHD1u6T+9embef1+fu8vnifC0O9MIR9Pl4iNBq07Xqjci+Kbv1f6E5YrlWjv3aRkMnDMFy7hKZ+K7TdfkKkidtC5rCAYLynbOVk5wW8ffSa474acg2ozzctS35MI9QqMlRwx/fYDfTBYaRwsqdgl8aktDVwatYWNk3Yn+S3pSUkKRztGDa3JU/v+zF31LZ4L99cc+sKdJJSPsA4WzoIOI0ZEwqRIYRoK4S4DSClXAmsBg4DdwA10Csu5SZH/ll3lk2LjtKydyUq1Cts0br0+7ZhuOWNpnE7VBmzfHLv0ZZTPN52hlztKuIxtCEvT97i0i8bIt1fKTQatB36oCpQBN26JfF+7oAMDSFs4zLCFk5FpEyFzYCxaCrWinFfaKRlScnzA1c42mIGDzefQu+Wlb/vQrEuFanU/vMABK6VPNAHhfLylDFOmU3hoqBWU6tqSnbNPspf/bcmmgNSkgLFyrnRvGcFtiw5zom91+K1bHM/DeWBqkKIHcAtjK24ElLKnuZWJKXMIaU8afp5uZTSLdy9kVLKDFJKZymlp5Qy8R1fbQXuej9jcr91FCmTi66j61i0LsPta+h2b0JVoizqb8p/cu/NlYd4T91O+tJ5yNOlKll+KEnenjV5vvcS16b/HfkYnEZrNLj8hdGtXYz+1JH40fn4PqHTR6E/ug91hZrY9B+LKlO2OJUV5POGc4NXcGHYKmzSpsS2cUW2HPDh2+bFaPBTpUjzpC2RC62jPT7/XgVA2NmjcnMnf3p/6g2qyIn1F5nTcQ3Bgcqoirl0HlGbnAVcmdh7DW9fffFI10fMMjcp5SugCsbtVjeB+lLK+LVZhU8I9A9meNu/cEhpy9i/2qLRWO7kKhngT+iKuYj0Lmibtv9ksiLE7x3nh67CLoMjhcc2R5gWDOdq8z05PcvzcMNJbv+5P9JyhdYGbce+qPIUJGzNIvSn435orzQY0O3/2xj6OzgIbbef0TbwRGgjX6ISbVl6A/fXHONoi5m8OnOHfL1rYd+8CqtnnqBQlTy0mfRDlBM2Ko2aDOUL8OLItY/dcpVHMeSLZ9TxzEvbKT9w9eAdpjVbxju/+HtQkzO2dlpGzvfkrV8g0wdtjLdyozU3IcQ1IYS3EMIb4xq3fEBp4G646wrxjJSSCb1W8/SeH7/81Y70rpY7uUoaDIR5zYf3gZ+Nsxl0ei4MX02YfxDFJrXGJvWns495e9Ukc90S3Fl8gAfrIx9+FVobtJ36o3JzJ2z1QvRnj8de46uXhM79Dd3fa1F5FDdOGuTziHU5AP43nnLix3lcn7kDp6I5KLuqH7o82fmzz2ZyFMlE1/lNUGui/853qeyBLiAYv9N3AFC7GzffG66eo3yr4vRY3Jwn118wsd4SfB8onRBzyFskCx2H1GD/pvPs3RA/sfJiWtnYLV5qUYgVa+cc4uC2S/QYV4+iZXNbtC79gR0Yrl9G07QDqsyfLgq+8ftOXl+4T+ExzXDM8/laNyEEBYc0IMz/PdembUeb2oFM1Yt8ns7GBu2P/QlbNM1opEKgLl7GPH1njxO2YRlIA9pWXVGVLBunZTC6oFBuL9rHg7XH0aZ2oMi4FrhWLcTTm77Mbr+atJlS03t5K2wdYm4Jpi/lhiaFLc//vYLzd/kQadMjMmVDf+U8mkp1KFo9HwPWtGV2+9VMrL+YPis8yV4o+rWCCtCqb2WO7b7KtIEbKPpdLpwzpfmi8qL9ipJSHorp9UW1K3zGxeN3mDd6OxV+KEzLXhUtWpfh7k10uzagKvot6jKfjjE93XWeB+tOkL35d2SqWTTKMlQaNUXGtcCpWE4uj12P74mbkaYTNrbGDeu58hHmNR/9hf+i1SaDAgldMZewlfMQGbNEG/o7JnyP3+Boy5ncX3WUzD+UoPya/mSsVpjXT/2Z5emF1lZDv1WtSZXWvHVxKhsNzmXz8+LwNQw64+SByqM48t5NZOA7ANxKZeXnrR3R2GiY0ngp3ofvxFr314ZGo2bEfE90YXom9FrzxTPPSrDKRISfj7/x5Koc6Rg6u4VFF+rKwHeErpiDSJsebfNOn9Tlf+MpVyZuxqlYTvL1jnmvqNpWS/EpbUiV25XzQ7x4fTnyA3mFrR02nQcicuQhbMVc9BdPR5pOf/saIZOHY7hwCk3tJtj0HBZt6O+oCPF7x4Xhqzk7YBlqOy3fzu+Cx5CGaB3tCXwTxMzWXgQHhNB3pSfps8auleBa2YOwt+95fd64AFVdsBhIicH74sc0Gd3SM3RbJ9JnS8PvbVZxclPCh9pOamTN7UzPcfU4deAGmxcf+6KyFHNLJOh0ekZ3WE6AfzC/Lm9v0ZOrpJSErVoI7/zRtu2NsPt/XaFv33N+iBdaRweKjm+BysyJDG1KO0rMaIedsyPnBi7j3Z3nkab7aHDZ3QhbPgf9pTP/1/Uh9PfcCQit1hT6uz5CHbvJFGkw8GjLKY40n8GLI9dw61KVsst741Q0h/E9BoUxu8MafO+/osfi5mQtGPutzOlL50Ftp+X5gSvG95UlBzimQR/hbIU0rqn4aWMH3EplY3HvzfwzP/Zjjl8bDTp+xzdV8jNn5DYe3n4R53IUc0sk/DVxDxeO3+Gnmc0sGikBQH9oNwbvC2jqtUSVNcfH61Jv4OKotQS/9KfYBE9s08Vut51tulSU/L0jKlstZ/ot5f3TyAfThZ09Nl0GIbLlJGzZbPRXzmF4/oTQmWPjHPr7AwH3fPiv2yKuTtyCY75MlF3ZB7eOlVHZGIeXDXoDi3pt4s7ph3Sc1ZD8ZXPGug4AtZ0NzmXz4XPIG6k3IFQq1AWLYbh+GakL+yStQ2o7+np5UqKuO+vH7WXtmD0YDMpi36gQQjD0j+bY2GoZ320VOl3c1g0q5pYIeHDTB69ZB6jRvCQ1mpeMOcMXYHhwG932tagKlURd/tNQPbcW7sXvv1u4D6xHGo+scSrfIZMTJWe2xxASxpk+SwjxexdpOqPBDUZkyU7Y0t8JnT4S+fYV2k790Tbt8NnuiJjQh4Rxa8FejrWZTeD9F3iMaEyp2Z1IkS39xzRSSryG7eTC7us0/6UmpeoVjNN7/IBLRQ9CXwXw+tIDAGOMt5BgDLevf5ZWa6uhy7wmVO74DfsWnWRx702EhVguIkZSxzlTGgZOa4z3mQd4zYzbNj7F3KyMlJIZP23Czl5Lz3H1LFvX+0BCl81BpHEyhtkON872/N8r3F12iCz1S5G1wZcdDZgqtyvFp7Uj5KU/Z/svIywgONJ0wt4Bm24/ocqZF1X+wtj+NAG1R/FY1+d35g7HWv/Onb/+JWO1wpRf258sdUt8Nma5Y+ZhDq88S61eZanSMe77Tz/g/F0+VDYafP41dk1VedxBaxPlsX8qlaDFLzVpNKwKp7Zc4fc2qwh6p4QxjIqqjYtTpXExlkyMWxgtxdyszL9bL3Lm4E06j6hN2gzxHnTlI1JKwtb8CW9fo23bE+GQ4uO9gHs+XB63gdTuWXAf+EO05RgMkUcFiYhToWwUneDJuzvPOTd4BfqQsEjTCfsU2PQchk3HfrEO/R36JpBLv2zgdK/FICUlf+9I4dFNsXH6PIDnkVXn2Dr1IGWaFKHhkCqxqicqNClsSV86Dz4HvZEGA8LGBlU+D/RXzkf5NxJCUKtnOTrMbMCt/x4wpfFSXj1N2CPvkhIDpzbGyTluz4ViblbkfUAIfwzbQt7CmWnQqWzMGb4A/dF9GC6fQVO3GarsH3e+ERYQzLmfvVDbaik20fPj2FRk6ML0TGu6jGnNlpu1f9K5TF4KjWrC6/P3uDhyzcdlE1+KlJInu85zpMUMnu25QK52FSm7si/pv3GLNP2Ff26w4ue/8ajkRtupUe8+iAsulTwIfvGWt96PAeOSEN74IZ9GPmP8ge+aFqHX0pa8uOfH2GrzOPu3sh4+MhydUjB0dos45VXMzYosnfIPvk/fMmBqE4ueg2B4dB/d1lWo3IuiDhcGXBoMXP5lPUFPXlH015bYZYi+5bRt6kFunnzAjeP32firefHaMtUoSoEBdXlx+BpXJ2754rVLgY/8ONN7CZfHridF1vR8t7wXebtXRx3Fea23Tz9iYfcNZC+cka4LmqLRxu82tgzl8iM0ap6b9pqq3YuCEMaTsWLAo5IbI3Z3JUOOtMzvup4l/bYo3dRI+LZK/jjlU8zNSty/4cPaOQep0/pbPL7JYbF6ZHAQYcv/gFSOaFt1+aTVcnfpQV4cvka+PrVIWzz6mcnrx+6xe85Ryrcq9nFQ/OwO81ob2Zt9R+5OlXny91luzo7b+IkhTMedpQc51noWb689xn1wPb5d0IVUuaNexvHslnH3gZOrI32Wt8IuRez3ocaENpU96UrlxufAFaSUiFSpEdlyo79i3hYi19zp+HlLR+r0Lc/JjZf4pfp8bp+OvtWnYB6KuVkBKSUzBm/EPoUt3cbUtWg9YesWI1+9xKZNT0SK/49d+B6/wa1F+8lYoyjZm0UfsjzwdRCL+2wmQ850NB9bk6Yjq5OzWGaWDdzGi3uvzNLi9mMVsjX+lnteR7i74nCs3sfrSw843n4Ot+b/g3PZ/JRb059sjUtHG+Lo9TN/ZrZaiVqrMu4+SJciyrRfimslD4Kevcb/xlPAdOzfo3vIt+btK9Vo1TT4qTI/bWoPwORGS9ky+YASOukLUczNCuzfdJ6zh2/RdVQdnNJb7vQq/Yl/MZz/D02tJqhy5f14/f1jPy6OXksqN1c8hjaIdgxKSsnyn7bzzi+QznMaYetgg8ZGTdcFTVGpVczrso7QoMgnC8IjhKDAwB9wrVqIm3N283jbmRjzhL0L4uqkLfzXdSG6gGCKT2lDsd9aYefsGG2+92+DmdXai/f+wfRd6YlzdqcY6/oSMnxfAKFWfQyDFNdj/9xKZWPUP90o06QIO2YdYVKDJTy/E9W5SgoxoZhbAvP+XTB/DN9K3iJZqNfevM3jccHw9CG6zStR5S+EuvL/Y8HpgkI597MXQgiKTfREbRd9V+3o6vOc23mNhj9XJnvh/y8uTpc5NZ3+aMhjbx9Wj9xlliahUlF4dFPSfZuHKxM343PwaqTppJQ833+Zoy1m8mjrabI3/45yq/uRoXyBGOsIC9Yxp+Mant95SY8/m5PNw/Ib1m3SpCBtsZw8/9fUNXXNjEjrbNa4W0TsU9nSYUZ9ui1oyov7rxhXYwGHV55VIvzGAcXcEpglk/bwyucdA6dZbhJBhgQTtmw2OKRA26rrx+6blJIrv20i4K4PhX9pjkPm6ENyP7/9kjWjdlOgXE6qdf2861qoch5q9ynP0dXnOb7uglnaVFoNxSZ6krpAFi6OWovf2U83lAc9f8O5QSu4MHw1tulSUWZJDwr0q4PGIeZzIwx6A3/23sTNkw/oOKshBcrHfodDXHGp7MH7hy8JuOuDEAKVR3EMN68iQyJf4xcTJeq6M2Zfd3KXzMqKn/9mToc1+L9U4sPFBsXcEpC73s9YP+8wddt+S8GSljlzVEpJ2PqlxlOg2vT4ZO3YgzXHeL73Enm6VsO5dN5oSgFdqJ5FvTahtdXQYWYDVKrIu671B1Uk33c58Bq6g8fXfMzSqLG3ocT0djhkTsu5wSt5e/0JBp2e+6uPcrTlTF6dvUO+PrUpvaQ7qfNnNvt9rxm1m3M7r9F8TA2+qR+3eG9xxaWCOwjxca+pyqMY6MIw3Iy8dWoOThkd6beqNc1GV+fqoTuMqTKPS/sij7qi8DmKuSUQUkqmD95ICkc7uo6yXMhw/akjGM4eQ1OjIWq3/3fj/M7e4cbs3bhUcCdXuwoxlrNl8gEeXn5Gu2n1cMoY9RiXSq2i85zG2DvaMb/LerOXMtikdqDkrA5oHe05238pJzvN4/qsnaQtnpNyq/uRs1U5szftA+z84yj/Lj1Nje7fUbVzabPzxRe26VLhVCT7x662Klc+sHOIcreCuahUgmpdyjB8Z2ccnVPwR7vVeA3dQYgZ45xfO4q5JRB7N5zjwrE7dB1dhzTpLDOJYHj2GN3GZajyuKOuVv/j9SCfN1wcsQaHLOkoNKpJjItYvQ/fZc+843zfugTFasa8xih1hpR0mdcE3wevWD54u9njQ3YZUlPq944AhLx8R5HxLSg+tS32GWM3AXBs7Xm2TDpA6UaFaTSsaswZLIRLJQ8C7vgQ8MAXodagKlAY/dULkR6iE1uyFHBh+N+dqdalNAeXn2F8jQU8uPQ0HlQnXxRzSwAC/YOZM2IrBYpno24by7QqZGiIcZzN1g6tZ/eP42z6kDAuDF2FPkRHsUmeaFJEvyH93av3LOm3GVe39DQbU8Ps+vOWzk6DnytzZvtV/l0aeZy2yEiRLT3l1vSn/IaBZKxaONa7By7tu8nywdtxr5CbdtPqRdl9TghcKxk34n+YNVV7FIcAf+TDu/FSvtZOQ7PRNRiwpg3BgaFM+GExO/84gkH/5eaZHFHMLQFYPGEXr14EMGBqY4tNIug2LUe+eIq2dTdEamPgRSkl16Zu5633YwqPbkrKHNEHfJRSsnzQNgJfBxmXfdhHvuo/Kmp0L0vhqnlZN3YPd889NjufTWoHNPaxX2B799xjFnRdT9aCrnRf2BSNjeUO0TEHuwypSe2R9f8b6fMXBpX6sxhvX0qB8rkYs687RWvkZ/PEA0xtuoyXj97Eax3JAcXcLMztK0/ZuPAo9dqXpkDxuB1BFxP6M8fQ/3cYddUfUOcr9PH6oy2neLz9DLnaVzQOeMfA4ZVnubDnBg2HVInTEgqVStBxZgPSuDqyoNsGAl69j3UZ5vL89kv+aLuK1K6p6LPCE7uUMc+mJgSulTzwv/GU909eIRxSoMqVF4OZuxViQwone7ouaEKHmQ14dPU5v1Sbz8mNl5QlI+FQzM2CGCcRNpAytT1dRlpmEsHw4hlh6/9C5MqHpkajj9ffXHnItWl/k750XvJ0jnkc6tktX9aN2YN7hdxfNCCfwsmebgua4u8bwJK+WywSlPHN83fM9FyJUKvo59Uax/SW230QW1w+dE0/TCx4FEc+f4LhZdwjykaFEILvmhZh1N5uZM6fgcV9NrOox0YC3wTFe11JEcXcLMietWe4dOIe3cfUJXXa+H8AZWgoYcv+AK0WmzY9Pobj/njWqEtqCo9t9vGs0agIC9GxsMdGbBxs6DCj/hePW+Uokolmo2tw+cAtds2O+1mlkfHe37j7IOB1EH1XtCJDjujX6iU0DpnS4pgv0/+XhBT8cOxf7Bf0motzNicGbWhPg58qcW7nNcZWnc/1Y/csVl9SQTE3C/HuTRBzRm7DvWR2arf+xiJ16LZ6IZ8+Mi7UTWN8yD85a3Si52dnjUbG5on7eeztQ/vp9UnjEj8x5Sq2K0mp+h5snfJvvD1oYSE65nZcy7NbvvT4s9knOyYSEy6VPHh79RHBL96iSu+CcMn8xUtCYkKtUVGn7/f8vLUjNvZapjdfzvpx/3zV0X4TzNyEEEeEEJH2d4QQNYQQt4QQgUKIdUIIy0VtTCAWT9jFm5eBDJzaGFU0G7zjiv78f+iPH0BdqbYxzI6J67OMZ416DGsY6VmjEbl66A57F56kUvtSFKkW/cLe2CCEoO3kurjkSseinht54xN5uHFzMRgki/ts5saJ+7SfUR/37y17nuuX4FrZuID4QxgklUcxDHduIIMsv8MgZ9HMjNzThe9bl+Cf+Sf4re6fPLkR/13ipIDFzc1kXPOBclHcdwJWA4OBLIAAxllalyW5dekJmxYdpUGn78hXNG5nEUSH4aUPYWv/RGR3Q1On6cfrT3ad5+H6E2RvUZZMNaI+a/QD7/wC+avfFjLlc6bJiGoxpo8tdilt6bawKcEBoSzqsRG9Lm5LFqSUrB29m7N/e9N0ZDVKNyocz0rjlxTZ0pMyt8vHWVO1R3Ew6DFcS5ij/WwdbGg9sS69/mrBW593jK+1kP2L//vqDqVJiJZbOYwn2wdEcb8ecE5KuUVK+RqYCMQt9GYiwGAwMH3wBhydHOg8ona8ly91Ycb1bGo1Nm17ItTGyLn+N55y9cNZo71qxlyOlCwdsI3At0F0nt0Ym1gu+zCXzPky0HpiHW6efMDWKXE76GP33GMcWHKKal1KU71b9OGZEgsuFQvy+uIDQvzeIbLlhpSOsY4S8qUUqZ6PMfu7U6BcLtaM2s3vbbx48/zLWtBJCYubm5RypJTyRyCq2C2FgfD/61cAFyGE5WIBWZDdq89w+b/7dB/7A45pzDvBPDbotq1GPr6PtkVnRFrjyU4fzhq1SZ2Cor+2NGvb0sFlp7m07yZNhlcji7tLvOsMT5kmRfjeszi7Zh/j4t7Y7Y08vu4Cm37bzzcNPGgysrqFFMY/rpU9QEp8Dl41HvvnXhSD90WkPmHHwBydU9J7WUs8J9Th1skHjK06j3M7ryWoBmuRGCYUHIG3H36RUoYAOtP1zxBCdBFCnBFCnPH19U0giebh/+Y9c0dvx+ObHNRq9WUnSEWG/tIZ9Ef2ov6+BupCJQDTWaMj1xD80p+iE1phmzbm74QnN16wftxePCq7UbmjZSY7ItLil1pk83BlSb/NZi84vXzgFssGbaNA+Vx0mBH15v3ESMpcLjhkS///GG8exSD4PYa7NxJcixCCim1LMnJPV9JlTcO8zutYOnArwQHJO6R5YjC3N8DHI8+FEDYYu7GRhjGVUi6UUpaUUpZ0dnZOIInm8eevu/B/FcjAqU3ifRLB8MqXsDWLEFlzovnh/732mwv24nfqNu6D6pGmYMzje2HBOhb13IhdSls6TK8fr4elRIfWTkPXBU2RBsn8rutjnMW7d/4J87usJ0sBF7ovamb13QexRQiBayUPXp2/R+ibQFR5PUCjxXAlYbum4XF1S8+QrZ2o3bscx9de4JfqC7hz5pHV9FiaxGBu14Hw8WkKALeklElqJeKNC4/YsvgYjX4sR57C5oXpMRep0xnH2aRE27YXQmMcZ3t+4Ar3lh8iS4NSZK1vXktx4297eXLtBR2m18fROWF7/hlypKXDjAY8uPiUdWP3RJnO564fv7ddReoMKem70hP7VIlj90FscalUEKk38OLwNYStHao87hiunLPqLgKNjZqGQ6owaGN7DHoDkxv9xbZpB+M82ZOYSQzmtgUoL4Qob5o5/Q1YZWVNscI4ibCR1OlS0GlYrZgzxBLdjnXIh3fRNv8RVXrj/tCAez5cHr+B1AWz4j4g+rNGP3B5/y32Lz5FlU7fUKhKnnjXaQ7FauanetcyHFx2hlNbr3x2/+2LAOPuAwF9V3omuAHHJ475MmGf0YnnH/aaehRHvvJFPn9iZWWQ99vsjNrbjW8bFmb79ENMarDE7PMwkgpWMTchRFshxG0AKaUf0A5YDjwCXgATrKErruxceYqrpx/Q45d6pEpjH3OGWKC/eh79wV2oy1ZBXdQ4PvbJWaMTWkV71ugH/H0D+GvAVjIXyEDjYfG/7CM2NBxaBbdSWVk+aBvPbv1/3DToXQizWnvx7mUgvZe3wiVXOiuq/HKEELhU9sDv9B3C3gV9XI9oyd0KscHB0Y6OsxrQZV4TfO768Uv1+RxZZd2WZdV1uCEAACAASURBVHySYOYmpcwhpTxp+nm5lNIt3L3tUsqcUsqUUsoOpkmFJIH/60DmjfmbwmVyUrNFyXgtW772I2zVQkTm7GjqtzJeMxi4PNZ01uhvrWI8axSMyz7+GmAcQO48uzFau5jN0JJotGq6zGuCjb2WeV3WE/I+lLAQHfN+XMvTGy/otrAZOYvGb9feWrhWKojU6Xlx9DoiTVpE1pzxHiXkSylVryCj93UnZ7EsLB+8nbk/ruOdBYMeJBSJoVuapFk4bicBb4MYMCXmIJCxQep1hK6YC3od2rY9EVpjSKA7Sw/y4sg18vWpTdpiOc0q68CSU1w5cJumI6uROX/0YY8SCqeMjvw4uzHPb/my8ucd/NV/K9eO3qPdtHp4VIr85PikSGr3LNhlSP3/Bb0FiyEf3EG+extDzoQlbSZH+q9pQ9OR1bhy4BZjq8zjyr+3rS3ri1DM7Qu4fv4hW/86QaPO5XDziN99jrrdm5D3bqJt2gFVBuM2Kt/jN7i9aD+ZahYlezPzTs567O3DhvF7KVw1LxXbxf/ylC/B/ftc1B1QgZObLnF66xUaD69KmSZFrC0rXhEqFS4V3Xl58ha6wBBUHsVBSvTe5h2ok5CoVILq3b5j2N8/ksLJnlmtvVg1YqdZRzcmRhRziyMGg4FpAzeSNkNKOg2NeUdAbNBfv4R+33bUpSuiLmFckR/4yHTWaB5XCg6J/qzRD4QGhbGw5//aO/OwKKs2Dt/PDCDuG+57KoviBpJLGmqaVi6V7ZlrmaV+pbZ+babti/VppVnutmqi2WJmmpVaimgCbrivKIIIKjDMzPn+eAcFBBmUYQY893V5ybzzzjs/5Z3fPOc553nOd5SvUpahU/oX27KPwtD3iRvpcn87+k0Ip/djJaP6oLDU6h6M3WIlYf0upG5DqFLdY/JuedGgZW1e/GkkPR/uwJo5m3jtlpkcijnublmFRpvbFfLDgn/YEXWI0ZMHUKFy0U0iqDOnyfxiBlKnPl53DAKMvUa3POf8XqNZLH7tV47vTmDYh7e7dMf1q8FkNjHkvf70H9/NI823KKjauhE+1SoQvyYGEcHcsi32XdEoi8Xd0vLF29eLe1/tw7ivBnE+JZ03+n7Oik/WlaiW5trcroAzSeeYMfEH2nZuSq+7Q4rsuspuJ3PhdLBk4D1kDOJTJsdeo20m30e5us71L/t35S7WzN1Er5EdaRnuuR00rgXEbKJWeAtOrd+FLd1iDE0tFuxx290trUBa3NiUiaseo02vAL57fRXv3zufxKOelS/MD21uhcBut7Ny0WYe7fkh51LSGf/ewCKNNqy/RGDfswPvgUMx1TJmCw985dhrdFQv/Do4tzYt+UQqcyd8T4OWtbnjuZuKTJ/myqndIxhbeian/o7D1CwIyvh69NA0OxWqlWPUzLsZOmUAB7cd59We0/lnabS7ZRWINjcnUErxxw/RDO3yHpMeWUiZsj68u+gRrmtR+H0G8sMWF4vt12WY2nfBfH1XwNhrdPfHK6jVrSXXDS54r1Ew+p7NGbcMy3kLj3x8J95l3LvsQ2NQtV0TvCuXI351DOLljSmgFbbtW4pk27/iQES44d62vPLrKOo2r8Hno5fw2ZglnD+T7m5p+aLv/MuglGLj6l189tpP7NxymAbNavDq7MF0v71NkdaO2hNPkjnvI6RGHbzvGgIYe41ufeEryjWoTquXnF9msuqzv9m+di+D3rqNOs09q/b2WsbkZabmjUHE/xaD3WLFHByCfdsm1JEDSMPr3C3PaWo0qsrTS4bx87Q/Wf7BWvZsPMTw/91OQKfG7pZ2CTpyy4d/N+xj7G0fM2HgpySfOsvzH93Hgr+f5aY72xWpsamMdDJnfQg2O94jnkTK+GLLyGTLc19gt9ho9/YgvMo7V1t5KCaeiLd+o23vAG4cFFpkGjVFQ+3uwdjOZ3Bq4x5MQW1ApNh7vBUFZi8TfceF89yyEXh5m3n/7nl898YqrBabu6XlQJtbLnZuPcyEuz5l9C3TOLz3JOPeHciXkf/ltkEd8HKiT1phUHY7mV98ioo/gveQ0Zhq1kEpxfb3vidlx1FaT7ybCo2ci74y0jL5bPR3lK9alsHveeayj2ud6mFN8argy4k1MUiFikgT/xKTd8uLJu3q8dLKR+nyQAgrPl7Hm/0+z1FO5260uTnYt+M4Lzw0h4e7TWHH5kM8Pqkf32x5kYGPdMHHRXkr68ql2KMj8ep/P+ZAo3X24YiNHF2+mabDulPrxoL3Gs3i24m/cGLvKUb87w4qViv6Jpmaq8fk7UXNrkGc/HMHdqvNqFY4ehB1+pS7pV0xvuV9GPxOP0bPvpekYylM7jOTNXM3ekR96jWfczuyL4E5b/3CykVRlK3gw/DnenPv490oX8nXpe9r+3cTtl8iMIV1wRxuLAI+HX2IHVN+wK+TP80edn6Wc8uKnfyxcDO9H+tMUNeSk7+5FqnVvSXHft5C0uZ9VAsOgeVfY/nodUyt2mMKDsHUxP/CFo0liba9A2nSrj5zJyzjyxd+ZtuqOIZOGUDlmu7r6iKe4LBXSvv27VVkZOQVvfbEkdPMfXclPy3ciJePmYEju/LgEz1csr9obuxHD2KZOgmp0wCf0f9FvH3ISExl/dCPMfl40XnOaLwrObcw+PTxFF7tOYPqDarw/PcjSlxTx2sNW3omq295nTq92xD83B3Ytm7EtvEPY82bNRPKVcDUog3mliGYAlshvkXbZcbVKKX4fV4kiyavpEw5H4a814+2vQOv+roislkpVajOFNdc5HY6IZUFU1axdPZ67HbF7cM789CEnvjVLri7RlGgzqZgmfUhlC2Pz7AnEW+fC63CralpdPxslNPGZrcr5jy5lMwMK498dKc2thKA2debGjcEcHLtdtTTAzC3vR5z2+tRGenYd0Zji4nCvn0r9sh1YPbC5N8CU3Ao5pbtkMpV3S2/QESE7kPDCLyhMZ+PjeDj4d/Q9YF23DOxD77lnausKSquGXNLST7PV9PWsHjGH1jSrfS5P4xhz9xM7YbFt2O5slqxzJ0KqWfwGfsiUrkKAPvm/U5S1H6CXxxIRSf2Gs3i10/Xs+Ov/Qx+tx+1m/m5SramiKndPZj4VdEkbT1A9VAjjSBlfDG3CcPcJgxls6H27zaMLiYK66I5WBfNQRpehzk4BFNwCFK7vkdPGtVpXoPnvx/B9++vYcXH69i1/iAjpt3BdSH1i01DqR+Wnk9N59sZf/D1tDWcTUnnpoHtGPF8Hxo2K/7WP5mL5mBbvxrvQaMwh94AwOl/D7Lx8c+ofVMrWr96j9M37MFtx3iz3yza3BzAqJl3e/SNrsmJ9XwGq295g/r9QmnxVP/LnquUQsUfxR4ThS02CnVwLwBSrYaRo2sV6vF5ut1/H2TWfyJIjk+h77hwbh3bFbNX4eYyr2RYWmrNLSPNQsSs9Sz8YBXJiefockswD79wS5G3JnIW67pVWBfPw9zjNrwdG7xkpqSxbvA0xGzihvlj8Crv3CRGxnkLk3vPJOO8hYmrHqN81ZKVl9HAlmcXkhx7mG7fP4sUYt2kOpOMLdaI6C7m6cpjCmprRHUemqc7n5LOVy/8zN9LttE0tD7Dp95BzcbOj5p0zg3ItFj5ceE/zHvvVxKOnSGsuz+PvHgrLUIbuU2Tfc8OrEsWYgpqg9dt9wDGN3LMm0vISEihw8xHnTY2gG9eWcHJ/YlM+HaINrYSSq3uwZxYu53kmMNUbe38vSmVq+DVuQd07mHk6XZFY4t25Ok2O/J0zVtgCg7BHBziMXm6cpV8GTHtDlr1bM7C535g0s2fcv/kPnS+p63LRh2lxtxsNjsrv93M7LdWcPxgEq06NOalTwcR0tW9XV3tSQlY5k5D/Gri/dDjF76ljyzbxIk1sQSM6ePUlnxZbP5xO39+uYVbxnQhoHNjF6nWuJoaXQIRbzPxq2MKZW7ZkTK+mFuHYW7tyNMdiHPk6TZjXTwX6+K5Rp6uZTtMwaFIHffn6a4fEEzT9g2Y88RS5o7/nm2r4njo7b5UcMHazBI/LN24cSNrl29j1hsrOLDrBP6t6/HIi7fSsVeQ23+RKiMdy9RJqKREfMZNvNBRN3XfCTYM+5iqbZvQ/oMhTg9Lko6e4dWbZ1CzcTWeXTocL2/PzbNoCmbzhHmk7j1BeMTTRduiXinUiWPYYzZji8kjTxccgum6ALfm6ex2xa+frifi7dVUqFaOYR/cftnWXNdczi3IP1h1qj2U3duO0jigFiP+24fwfq2LfEPkK0HZ7WTOm4Y9ejPeI5+6UIFgS89kw4hPsCSd5YaF/6FM9YpOXc9uszPl3gUc+PcoL68cRc0mxTfLq3ENR37YTMxr39H2jQeoFd4CMbvmvlUpydhitxh5ut2xufJ07TAFtnZbnu5QTDyfj13C8d0J3DSiAwP/2zPPDYyuOXOr4FVL3Ro8luHP9abX3aGYXXRzXAnWFUuw/hKB14AH8Op2cS/T2HeWcXjJP4R+OJQaHf2dvt5P0/4k4q3VDJ0ygBvubesKyZpixnLmPH/d9wGW0+fwrlQWv47N8esUgF+H5pSp5pqV/UaeLgZbzGbs27fCubPZ8nTtMLcMQaoU7xenJS2TJW+u4rdZG6kbUIOHp91Jg5a1c5xzzZlbs8aBamdcrMcNz2z/biJz7lRMYV3wvn/khSFH/JoYtj7/JY0f7ErgWOc3b96/9ShvD5hNu1uCGDm9aBtkatxLZkoap/6JI2HDLk5t2I3l9DkQoVJgXWp08sevUwBVWtR3SVRn5On2GEYXE4U6dQIAadDk4nq6Og2K7X6L+X0Pc8ct41xyGnc824OeIzthMhnvfc2Z29WUX7kK+7FDWP73ao7SKoC0+GTWPTSVcvWr03Hmo5i8nZvLST+bweTeM7FarLz86yjKF/GmzxrPQdntpOw+zqn1u0jYsJvk2MNgV3hXKkv1Ds0Ns+vo75KoTimFOnkMe3SUkac7tBeUQqr5YQoOdeTp/BGza+cgU5POM//p5WxdsZOAzo0Z/uHtVKtXWZubu1FnU8iY8grYbZQZN+lCBYLdamPT6M9JiTtO53ljKd/A+Z3U54xbxobF//LU4iH4d3DfchZN8WM5c57EjXtI2LCbU3/vxpJ0FoBKQfVcH9WlJGPbvhV79GbscbGQmZWna+NYT+e6PJ1SinXfbOXrl37G7G1m0Ft9uX5AsDY3d6FsVizT30Yd3IvP2BcxZeuuGvfZKvbOWk3rifdQt4/z+bJN38cy87HF3PZEV25/pocrZGtKCE5FdR2aOz1BVaj3duTp7LFR2GK3wrlUR54uyFhP56I83ckDScwaG8G+qCN8fmyiZ5qbiPQGPgLqAj8CI5RSqbnOWQt0yHboS6XU8Mtd15PMLXPxXGzrfstRWgWQFLWPjWNmUbdPW1q/fLfT10s8ksyrvWZQu6kfz0QM87i8osa95BvVBWZFdf5UadmgyKM6Zbej9sddqJJQCfEASP0mmFsVfZ7OZrXz49Q/GDChu+eZm4hUBfYCw4G1wEzgqFLqyVznHQaaKaUynL22p5ibdd1vWBfPzVFaBcYNuG7QVMy+3nSeNwavcs61C7fb7Lx39zwOx8bz8spR1GjkGavMNZ7Jhahuw24jqos5dDGqu74ZNToF4Nex6KO6HHm6rLrXrDxdyxBMrbLW0119ns5Ty6/6A1FKqaUAIvIWRvR2wdxExBewFsbYPAWjtGpBjtIqcJRXvf4dltPn6Pj5KKeNDeCnaX8R94+x8YY2Nk1BiMlE5cB6VA6sR9Nh3Y0Z2I1xJKw3orr4VcY2fEUd1YkIUqseplr18OrZ72KeLiYK299rsP25EsqWc+TpQjEFFe96uuKI3N4H7Eqppx2PywDpQEWl1FnHsUDgN+AE0BgjwntMKRV/uWu3axuitmyNcqH6y2NPSsAy5RWkfEV8nnwFKXuxhOTg4g3seG85gU/cSuP7uzh9zb2bj/DOHbNp368lD390p172obkq3BbVZaRj3x2DPWaLsQnOuVQwmzE1y1b3Wog8nUfOlorIZ8BBpdRr2Y5lAo2UUsccjzsBk4AngCPAVKCeUqrX5a5do0w99cYj0+g+JIym7Yu3bi6/0iqA1LjjbBgxnertmxLy/mCndaWlZjDp5hkoBS+vfJRyLm51rrn2yIrqTm3Yzam/48hINFLflQLq4tfJnxqdAqjcsj6mItwMSdnt2epec+XpgtsZebq6DS/7OfFUc3sXsCilXnA89gEygHJKqbR8XuOHEcVVyH2OiIwERgLUrlIv9O5aY0hLzaB+i1p0HxJGhztbUaacazt+GqVVH2GPjsxRWgVgTbOwYdjHWFPTuWHhWHyqOr8madZ/IvgnIppnlgylWVhDV0jXaC6g7HZS4+KNSYkNuzgdXTxRnf3EMWPoGhOFOrgHlIKqfhcWDpuaXpqn81RzGwH0V0oNcDxuAyxSSvlnO+dWDANc5XhcBziAYYD5bobYvn179dfa9WyMiGbN3E0c2XGCshXL0PmeNoQPDqOOi7rTWn+JwLpiySWlVQAxbyzhyPLNhE0dRvUw5zqSZGZY+XXmBiLeWk2/CeH0H9/NBao1msvjlqgu9Uy2utcYYz3dhTxdCKagNohvWY81t+pAHDAAiAEWApuUUhOznTMImAj0BY5hLBsxK6UevNy1s8+WKqXYG3mYNfMi2fxDLLZMO0FdmtBtSBhtbg4odOfPy2Hftwvbtki8BjyQI5Q+vmob/774NdcNCcf/sd4FXicz3cpfX0fx07S/SI5PpdVNzRk9+74i1arRXAm5o7rkmMMom92I6sKaOczOv0ijOmXJMNbTxURh274Fzl7M05V57FnPMzcAEemHkUerASwCRgH3Ai8rpZo5zpkMPAoIsBIYo5Q6fbnr5rcUJCXhLH9+tYU/FkSSdCyFKrUrEj4olC4PhFClVtEvcgQ4f+w06wdPo3zjGnSYMfKy326Z6Vb+/CqKnz8yTK1ZWAP6je9GUNcmegJB45FkpqblWFeXccq1UZ2Rp7tY9+r7wnueaW6uoqB1bjarnejfdrNmXiTb1+7F7GUi5NYgug0Jo3mHyycwC4PdauOfUTM5t/8knReMpVzdvGeBLGmZ/PllFCs+WUdyfCrNOzSk37hwArtoU9OUHJRSpMYdN4xu/cWozquiL37XN3dJVOeRw1JXUphFvPF7E1m7IJL1327l/Jl06gXWpNvg9nQc2BrfCs6vQcuL3dNXsm/e77SZfB91erW+5HlLWiZ/fLGZFZ+s48yJs/h3bES/8eEEdG6sTU1T4sk3qvOve8HoKgc3uKqoTpubE2SkZbJxaTS/z93EoZh4fCv40OmuNoQPbk+9gMLviJUYuZdNY2dTv28owS/cecl7/bEgkl+mr+fMybP4d2pE//HddHtwTaklR1S3YTfJ0YeyRXXNjH51HZvj61epUNfV5lYIlFLsizrK7/M2Ebk8FqvFRkCnxnQb0p62fQKdquW0nD7Luoem4VXel05zR+NV1liCkpGWydr5kfwyfR0pCecI6NzYiNQ6Nb4irRpNSSUzNY3ETXsumN2VRnXa3K6Q1MRz/PX1FtbOjyTxyBkq16rAjQ+GcuODoVSpnXfeQClF1FPzSdy0l46zHqNS8zpknLewdn4kK6avJ/XUOQJvaEK/8eH4d9StijQapRSpe+JJWL+r0FGdNrerxG6zE716D7/P30Tsmj2ISWjXJ4huQ9pfkh878M06dn7wI0ET+lH7tlB+n7eJX6avJzXxPEFdr6Pf+HCaX68X4mo0+ZFfVFfRv87FGtjghpi8zNrcipKTB5JYuyCSdV9v5VxyGnWa+9FtSBid7mpD5rFENjw8nWphzUhs2pSVn24gNfE8LW40TE1XF2g0hSMrqjNqYHeRvM0R1VXwxa9DM9q98aA2t6LGkpbJpu9jWTNvEwf/PUaZct70aFOJsomJrDnly+mkdFqEN6XfuHCahTm//6hGo8mfzLPpjhnYXZz6O44ePzyvzc2V7N96lN/nRbJx2Ta8Mm00C29O3ydvpGl7bWoajatQSmEymbS5FQdnk85z9nQatZs6vxeCRqO5cjy1WWWpo0K1clSoVq7gEzUajdvQFdoajaZUos1No9GUSrS5aTSaUok2N41GUyrR5qbRaEol2tw0Gk2pRJubRqMplWhz02g0pRJtbhqNplSizU2j0ZRKtLlpNJpSiTY3jUZTKtHmptFoSiXa3DQaTalEm5tGoymVaHPTaDSlkmIxNxHpLSJxInJORL4VkUv2y3PmHI1Go3EWl5ubiFQFvgKeBuoDAkwu7DkajUZTGIojcusPRCmlliqlTgNvAfddwTkajUbjNMVhbq2BLdkexwC1RKRCIc/RaDQapykOc6sEnMl6oJTKAKyO44U5R6PRaJymOHa/SgbKZj0QER/H+54u5DlZz40ERjoeZohIjAs0FzV+wCl3i3CCkqCzJGgErbOoCSjsC4rD3HZi5NSyCALilFJphTwHAKXUTGAmgIhEFnYvQ3egdRYdJUEjaJ1FjYgUeoPi4hiWLgW6ikhXx6zoG8CXV3CORqPROI3LzU0plQgMAeYDh4GTwJsiMlhE9lzuHFdr02g0pZdi2XFeKbUcWJ7r8HzHn8udUxAzr1JacaF1Fh0lQSNonUVNoXWKUsoVQjQajcat6NpSjUZTKimR5iYi14nI74461O0i0s3dmvJCRJ4WkSMictahN9DdmrIjIn+KSEfHzyIiU0TkjIgcE5FR7taXRS6d5URkroikiEi8iPxPRLzdrRFy6sx1/HURWeEOTXmRW6eIDBKRA47ffYSIVHGnPoem7L/zGiKy1KEvwfE7Nxd0jRJpbsC3wGqgJvA28IV75VyKiHQAxgNdgerAP8AnbhXlwNGkYAbQJdvh4UBPoAXQG3hNRNq5Q18W+eh8DgjE0NkWuB543A3yLpCPzqzn2gETil/VpeSlU0TaAlOA+zHqun2AZ9yjMN//y0kYa2FrY/zObwIeLOhaJc7cRCQEY+HhZKXUOaXUPOAeZ5y8mLEBiouTNmYgyX1yctAFQ9fZbMcGAe8opY4qpaIxGhnc6w5x2chLZy/gXaXUEaVUPDAL42Z3J3npxBFRfobnJO3z0jkcmKGU2qCUSgWGkW2izw3kpdGK4VUmjKYaTn2WimW2tIgJBeKAeSJyK7AfeFwpZXOvrJwopSJF5Htgt+NQGsa3jttRSr0EICI9sx3OXd8bixHJuY18dD4MHMr2uD1wrDh15SYfnQDPA78BOwD/4taVm3x0hgJrRCQKaIixYmGMG+QB+WqcBERx0fB+A34s6FolLnLDiNpuAn4BGgBzgKUi4utWVbkQkZuAO4AQoCIwFw8cPmcjR30vkIoH1vYqpWKVUqkiUl5EpgF3YqQmPAoRCcaIfCe6WUpB+AF3AfdgGHBDDDPxJKZgpHWqYKQkGuFEKqIkmhtAtFLqC6VUmlLqE8CbK6g9czG3AwuUUluUUmcxckXtRaSam3XlR476XqAcedT2egIi0gMjGroOCFFK7XezpByIiAmYDYzJq4TQA/lUKbVHKZWEYSQ93C0oF7cDE5VSZ5RSu4BpGHnhy1ISh6UHMMwMMGb5MEza026i8xi5gSysGHk4T9OZxU4gGGPID0bC/l/3yckbEekNfA2MUkp94249+VARaAf8bNyemAGTiMQrpWq7VdmlHCDb5wnP/SxlJxM4V9CLSmLktgKoKyIPi0hZjPzAQS5+KD2FH4BBIhLi0DkRWOHB3+RfAE+JSDXHTO+DwCI3a8qLV4BxHmxsOCIMb6WUr1LKF3gE+NUDjQ2ML4rHRSTQMaoYD3ja/+1y4BURqSIijYDHgO8KelGJMzdHp96bgUeBExgh613Kw0otlFJ/YiSUvwXiMXIFI9wq6vJ8BmzF+KJYDDyplPK0LwwwJj4+FxFrtj+/uVtUCWYuMB0jSb8d2AR85E5BeTAeIwe8F/gL+EoptbigF+nyK41GUyopcZGbRqPROIM2N41GUyrR5qbRaEol2tw0Gk2pRJubRqMplWhz02g0pRJtbporQkS6ichOd+sAEJGnHP3dhhdwXmMRSS8uXRr3os1N43FcQfuq24CnlVKzXaFHUzLR5qa5Gkwi8omjI/J+EekMICJeIvKuiBwXkaMi8rKjBhhHx9fsXWAPiEhHR1QVLyKvAQm530hEKojIbBE55XivRx3HZwHdgE/yitxE5BEROSEiiRhlO1nHTSIyVUQSHR1evxORSo4ypHQRKZ/t3COOYn1NCUKbm+ZqaI5ReO2Hsc9sVqucp4EOGP3rbgAGYvRhK4iaGM0c6ubx3PtAeaAZxgbeL4pIL6XUCGAt8GDuyE1EmgHvOt6/KUYzgCxuwTDFrBY6NYAhSqmdGD0Cezqu0RbwBf5wQr/Gg9DmprkaUjC64qYBPwGNHceHYLSoOaGUOgC8g2EwzjBJKZUjL+aI+h4CnldKJTs6Bc9w4pr3AouVUn8ppZLJ2adsI9DX8W+oBaRjtIMHiMAY6uL4e7lSyuqkfo2HUBJbHmk8h/hsDQtsXLyfGgF7sp2XgBGV5UX2L1iLUip3exscry2DEVFlv2ZoAfrqk7Nrb/aOvZUwGp3Ww+gNVzXbcxHAEsfPtwJvFfA+Gg9ER26aqyG/rgsJGMaSRdbwFcCO475zTBxUp2CSHO9VJ59r5scJjM6yWWT/eRKwXinVVCnVF6PjRBaRDn3dMXrc/eqERo2Hoc1N4woWAS+LSHURaQL8h4ubjhzB2PcA4D6Mjr+XRSmVCSwDXnck/dsCg4EFBbz0O2CgiHRybFf3QrbnvIAyIuIrIn0xOrt6i4g4otFlwAcYPfj08pESiDY3jSt4GWMIGIfR+/4rpdRSx3NPASNFJBojoR/j5DVHYQwdj2I0An1FKbXlci9w5OYmYPSn24Wx70YWkzF20koEHsBoejoOo18cGEPTNhj9+DQlEN3PTaPJyYBQIgAAAGFJREFUAxGpg9G8sY6O3EomOnLTaHLhWOM2AiPi1MZWQtGzpRrNpURgTD4UuMOSxnPRw1KNRlMq0cNSjUZTKtHmptFoSiXa3DQaTalEm5tGoymVaHPTaDSlEm1uGo2mVPJ/8DnIhutVMsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 324x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fontsize = 12\n",
    "\n",
    "fig = mp.figure(figsize=(4.5,4))\n",
    "\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(unique_hours+0.5, pblh_h_means/1000., label='h', c=mcm.magma(0.2))\n",
    "ax.plot(unique_hours+0.5, pblh_ll_means/1000., label='Liu and Liang', c=mcm.magma(0.3))\n",
    "ax.plot(unique_hours+0.5, pblh_br25_means/1000., label='Bulk Ri 0.25', c=mcm.magma(0.5))\n",
    "ax.plot(unique_hours+0.5, pblh_br50_means/1000., label='Bulk Ri 0.5', c=mcm.magma(0.7))\n",
    "\n",
    "ax.tick_params(labelsize=fontsize)\n",
    "ax.set_xlim(6,18)\n",
    "ax.set_xticks([i for i in range(6,19,2)])\n",
    "ax.set_ylim(0,3)\n",
    "ax.set_yticks([i for i in numpy.arange(0,3.1,0.5)])\n",
    "ax.set_xlabel('hour of day', fontsize=fontsize)\n",
    "ax.set_ylabel('km', fontsize=fontsize)\n",
    "\n",
    "ax.text(s='PBLH at Brazil ARM site (2014-2016)', x=0, y=1.02,\\\n",
    "        ha='left', va='bottom', transform=ax.transAxes, fontsize=fontsize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}