global_precipitation_energy_constraint.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# L δP = δRtoa - δRsfc - δH (O'Gorman version => fluxes positive upward)\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import esmlab\n",
    "from pathlib import Path\n",
    "import matplotlib.pyplot as plt\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "dloc = Path(\"/Volumes/Jedha/CMIP6\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rlus_Amon_IPSL-CM6A-LR_amip_r1i1p1f1_gr_201501-201712.nc\n",
    "def cmip_loader(loc, v, t, m, c, r):\n",
    "    g = \"_\".join([v, t, m, c, r, \"*.nc\"])\n",
    "    fils = sorted(list(loc.glob(g)))\n",
    "    if len(fils) == 0:\n",
    "        print(f\"No files found: {g}\")\n",
    "        return None\n",
    "    elif len(fils) == 1:\n",
    "        return xr.open_dataset(fils[0])[v]\n",
    "    else:\n",
    "        return xr.open_mfdataset(fils)[v]\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get one variable, global average, and put into dictionary\n",
    "def get_glbavg(VAR):\n",
    "    X = {}\n",
    "    shortnames = {\"CESM2\":\"cesm\", \"CNRM-CM6-1\":\"cnrm\", \"HadGEM3-GC31-LL\":\"hadgem\", \"IPSL-CM6A-LR\":\"ipsl\",  \"MRI-ESM2-0\":\"mri\",\n",
    "                 \"amip\":\"amip\", \"amip-lwoff\":\"amip-off\", \"aqua-control\":\"aqua\", \"aqua-control-lwoff\":\"aqua-off\"}\n",
    "    for M in [\"CESM2\", \"CNRM-CM6-1\", \"HadGEM3-GC31-LL\", \"IPSL-CM6A-LR\",  \"MRI-ESM2-0\" ]:\n",
    "        for EX in [\"amip\", \"amip-lwoff\", \"aqua-control\", \"aqua-control-lwoff\"]:\n",
    "                if MODEL == \"CNRM-CM6-1\":\n",
    "                    REALIZATION = \"r1i1p1f2\"\n",
    "                elif M == \"HadGEM3-GC31-LL\":\n",
    "                    if EX == \"amip\":\n",
    "                        REALIZATION = \"r1i1p1f3\"\n",
    "                    elif EX == \"amip-lwoff\":\n",
    "                        REALIZATION = \"r5i1p1f2\"\n",
    "                    elif 'aqua' in EX:\n",
    "                        REALIZATION = \"r1i1p1f3\"\n",
    "                else:\n",
    "                    REALIZATION = \"r1i1p1f1\"\n",
    "                fullfield = cmip_loader(dloc, VAR, \"Amon\", M, EX, REALIZATION)\n",
    "                if fullfield is not None:\n",
    "                    w = np.cos(np.radians(fullfield['lat']))\n",
    "                    globalmean = esmlab.weighted_mean(fullfield.mean(dim=('lon','time')), dim=['lat'], weights=w).values.item()\n",
    "                    k = shortnames[M]+\"-\"+shortnames[EX]\n",
    "                    X[k] = globalmean\n",
    "    return X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlut' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No files found: rlut_Amon_MRI-ESM2-0_aqua-control_r1i1p1f1_*.nc\n",
      "No files found: rlut_Amon_MRI-ESM2-0_aqua-control-lwoff_r1i1p1f1_*.nc\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlutcs' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No files found: rlutcs_Amon_MRI-ESM2-0_aqua-control_r1i1p1f1_*.nc\n",
      "No files found: rlutcs_Amon_MRI-ESM2-0_aqua-control-lwoff_r1i1p1f1_*.nc\n"
     ]
    }
   ],
   "source": [
    "rlut = get_glbavg(\"rlut\")\n",
    "rlutcs = get_glbavg(\"rlutcs\")\n",
    "\n",
    "lwcre = {}\n",
    "for i in rlut:\n",
    "    if i in rlutcs:\n",
    "        lwcre[i] = rlutcs[i] - rlut[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlds' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No files found: rlds_Amon_MRI-ESM2-0_aqua-control_r1i1p1f1_*.nc\n",
      "No files found: rlds_Amon_MRI-ESM2-0_aqua-control-lwoff_r1i1p1f1_*.nc\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rldscs' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No files found: rldscs_Amon_MRI-ESM2-0_aqua-control_r1i1p1f1_*.nc\n",
      "No files found: rldscs_Amon_MRI-ESM2-0_aqua-control-lwoff_r1i1p1f1_*.nc\n"
     ]
    }
   ],
   "source": [
    "rlds = get_glbavg(\"rlds\")\n",
    "rldscs = get_glbavg(\"rldscs\")\n",
    "\n",
    "lwcre_sfc = {}\n",
    "for i in rlds:\n",
    "    if i in rldscs:\n",
    "        lwcre_sfc[i] = rldscs[i] - rlds[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cesm-amip': -27.466924652403804,\n",
       " 'cesm-amip-off': -22.797060945353394,\n",
       " 'cesm-aqua': -40.24013889325323,\n",
       " 'cesm-aqua-off': -35.21459054792166,\n",
       " 'cnrm-amip': -25.148659947460658,\n",
       " 'cnrm-amip-off': 0.0,\n",
       " 'cnrm-aqua': -31.480549057314988,\n",
       " 'cnrm-aqua-off': 0.0,\n",
       " 'hadgem-amip': -25.996770555103183,\n",
       " 'hadgem-amip-off': 0.0,\n",
       " 'hadgem-aqua': -33.61595103938225,\n",
       " 'hadgem-aqua-off': 0.0,\n",
       " 'ipsl-amip': -29.469329833984375,\n",
       " 'ipsl-amip-off': 0.0,\n",
       " 'ipsl-aqua': -43.71173095703125,\n",
       " 'ipsl-aqua-off': 0.0,\n",
       " 'mri-amip': -27.455639596185165,\n",
       " 'mri-amip-off': -23.948713445837143}"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lwcre_sfc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAFVCAYAAAAUpK/DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAyBklEQVR4nO3deZwlVXn/8c93WBNBFhkXYGCAgArI5ogoibsBI4pL1FHEDQMajGBwA39xzSTighETNShBRQxBAQUBFQlLkM0ZHBiGER1llAmgg2EzIGF5fn+cukxNz+3uWs503zl8369Xv7pv3apzn66+/VTdU6fOo4jAzMzKMmO6AzAzs/yc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7jYtJL1J0qUN171I0ls7vs5XJP19l21bvMZBkn6wJl/DrC0n98JVSXSRpHsk3SrpC5I2rT3/YUlfH2fbZZLulfT7atuvSNpoyoKfApK2k3SxpLur3/cNk6w/W1JIWnewLCJOiYg/X/PRmjXn5F4wSUcBxwLvATYB9gG2Bc6XtH7DZl4SERsBewB7AkevgVCn0z8Ay4DNSfvn+mmNZgTUD1y29nJyL5SkRwMfAf4mIr4XEfdHxDLg1aQE//o27UXErcD3SUm+aQx/LukGSXdK+nx1hjy0e0XSMyX9uFr3x5KeOWaVHSRdVT3/HUmb17b9ZvXJ4k5Jl0japcWv9gCwvNo/t0bE/EnWv6T6fkf1ieYZY7uYqjP7v5b08+oTwcck7SDpckl3STqtfnCVdICkhZLukHSZpN3Ge3FJn5V0U9XOAkl/Vi3fsvqUVd8ve0q6TdJ61eO3SFoi6XZJ35e07ZiYD5f0c+DnE71W9dwfSfpq1dYSSe+VtLz2/JaSTpe0QtKNkt45yX61zJzcy/VMYEPgjPrCiPg9cB7wwjaNSdoaeBGwtOH6WwDfIp3pPwa4oYpp2LqbA+cAx1frHgecI+kxtdXeALwF2JKUkI+vPXcesCPwWOBq4JSmvxdwFfBuSfs3XP9Z1fdNI2KjiLh8nPX2B55K+jTwXuAE4CBgFrAr8FoASXsB/wYcRvrd/xU4S9IG47T7Y9IBdnPgG8A3JW0YETcDlwOvrK37OuBbEXG/pJcBxwCvAGYC/wX8+5i2XwY8Hdh5oteqnvsQMBvYnvReevhkQdIM4GzgGmAr4PnAkZL2G+d3sjUhIvxV4Bfpn+3WcZ77OHB+9fOHga+Ps94y4PfA3UAAF5CSWpPXfwNwee2xgJuAt1aP3wRcWv18MHDVmO0vB95U/XwR8PHaczsD/wesM+R1N61i3aR6/BXg78eJcV/gRuDZwHJgv2r5jsBtgIZsM7tqf93asod/l+pxAPvWHi8A3ld7/Gngn6qfvwB8bMxr3AA8u+F+vh3Yvfr5rcB/jtnfz6oenwccUttuBnAPsG0t5ue1eK1fDvZX7bWXVz8/Hfj1mG2PBk6a7v+LR9KXz9zLdRuwxTj9p0+onm/iZRGxMfAc4EnAFg2325KUXACI9B++fIJ1fzVm2a9IZ30DN415bj3S77eOpI9L+oWku0gHJBrG+Q7g5Ii4GHg5cHJ1dvlM4IIq5q5+U/v53iGPBxemtwWOqrpk7pB0B+nsfsthjUo6quoGubNadxNW/q7fAp4haUvSJ4wgnaEPXueztdf4H9IBYLx9PNlrbTlm/frP2wJbjvmdjgEeN+x3sjXDF07KdTlwH+lj+GmDhZIeRepeOaZNYxFxsaSvAJ8ifXyfzC3A1rXXVf3xGDeTEkLdNsD3ao9njXnuftIB6nXAgcALSIl9E9IZphrEuC6pi4eI+LGkuaRurNuBl46zTe5pVG8C5kXEvMlWrPq830fq5lgcEQ9Jevh3jYg7lIZkvhp4MvDvtQPU4HUm6rJ6+Heb7LVY+fcdXICu/31uAm6MiB0n+51szfGZe6Ei4k7SBdXPSdpf0nqSZgPfJJ1Bn1xbfYakDWtf4/X3/hPwQkl7wMPDLJeNs+45wFMkvaz69HA48Phx1j0X2EnS6yStK+k1pK6X79bWeb2knSX9MfBRUl/yg8DGpIPY74A/Jo1+aeqbwDslPavqJ76FdIB4HOmTwTArgIdIfc05fAl4m6SnK3mUpBdL2njIuhuTDkYrgHUlfRB49Jh1vkHqEntl9fPAF4GjVV1slrSJpFdNENdkr3Va1d5mkrYifQoauAq4S9L7qguv60jaVdLTJt4VlpOTe8Ei4hOkM/RPAXcBV5LOqp4fEffVVn0tqatg8PWLcdpbAXwN+Ltq0SzgR+OsexvwKuATpMS7MzCflIjHrvs74ADgqGrd9wIHVG0MnEzqP7+VdKF4MPria6Rumv8mnUVeMSyecWI8DXg/6WLnHaQLjJ8hDR39rqRthmxzDzAP+FHV5bBP09cbJ4b5wF8B/0z6xLCU1Ic/zPdJfec/I/3Of2BMVwpwFumawW8i4pra65xJGhZ7atV9dR3pE9x4Jnutj5JOEm4EfkjqErqveq0HgZeQLsbeSPqE9WXSpyqbIurXrWiPZFUXwBERsaTBujNIyeCgiLhwjQdnU0rS24G5EfHs6Y7FEp+5W2cR8ecTJXZJ+0natOrmOYbUX9v4zNpGl6QnSNpX0gxJTyR96jpzuuOylXxB1dakZ5D6fdcndZm8LCLund6QLJP1SWPytyN1aZ0KfH46A7JVuVvGzKxA7pYxMyuQk7uZWYFGos99iy22iNmzZ093GGZma5UFCxbcFhEzhz03Esl99uzZzJ8/2WR8ZmZWJ2nstB0Pc7eMmVmBnNzNzArk5G5mViAndzOzAjm5m5kVyMndzKxAkyZ3SbMkXVhVZFks6Yjac3+jVAB5saRP1JYfLWlp9ZzrJpqZTbEm49wfAI6KiKurAgILJJ1PKmhwILBbRNwn6bEAknYG5gK7kEpx/VDSTtUcz2ZmNgUmTe4RcQupQg0RcbekJaS6i39FKlo8mKD/t9UmBwKnVstvlLQU2JtU9s1srTT7/ed02m7Zx1+cORKzZlrdoVqVaduTVNHnk8CfSZpHqtLy7oj4MSnx1+fsXs6qRXjNzLIZlQPvqMQx0Di5S9oIOB04MiLuqupibgbsAzwNOE3S9gwvTLzavMKSDgUOBdhmm9WqmZnZiBu1ZGarajRaRtJ6pMR+SkScUS1eDpwRyVWkosFbVMvrldC3JlW3X0VEnBARcyJizsyZQ+e9MTOzjiY9c5ck4ERgSUQcV3vq28DzgIsk7USqzHIbqUDvNyQdR7qguiOpGrqZZdDljNlny488Tbpl9gUOBhZJWlgtOwb4N+DfJF0H/B/wxkhlnRZLOo1UVu0B4HCPlDEzm1pNRstcyvB+dIDXj7PNPGBej7jMzKwH36FqZlagkSjWMQrcj2lmJfGZu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcijZWwozxtitnZzcreR5oOMWTdO7iPGyczMcnCfu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcjJ3cysQE7uZmYFcnI3MyuQk7uZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBZo0uUuaJelCSUskLZZ0RLX8w5L+W9LC6usvatscLWmppBsk7bcmfwEzM1tdk2IdDwBHRcTVkjYGFkg6v3ruMxHxqfrKknYG5gK7AFsCP5S0U0Q8mDNwMzMb36Rn7hFxS0RcXf18N7AE2GqCTQ4ETo2I+yLiRmApsHeOYM3MrJlWfe6SZgN7AldWi94h6VpJ/yZps2rZVsBNtc2WM/HBwMzMMmuc3CVtBJwOHBkRdwFfAHYA9gBuAT49WHXI5jGkvUMlzZc0f8WKFW3jNjOzCTRK7pLWIyX2UyLiDICI+E1EPBgRDwFfYmXXy3JgVm3zrYGbx7YZESdExJyImDNz5sw+v4OZmY3RZLSMgBOBJRFxXG35E2qrvRy4rvr5LGCupA0kbQfsCFyVL2QzM5tMk9Ey+wIHA4skLayWHQO8VtIepC6XZcBhABGxWNJpwPWkkTaHe6SMmdnUmjS5R8SlDO9HP3eCbeYB83rEZWZmPfgOVTOzAjm5m5kVyMndzKxATu5mZgVycjczK5CTu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcjJ3cysQE7uZmYFcnI3MyuQk7uZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBXJyNzMrkJO7mVmBnNzNzArk5G5mViAndzOzAjm5m5kVaN3JVpA0C/ga8HjgIeCEiPhs7fl3A58EZkbEbdWyo4FDgAeBd0bE99dA7GaNzH7/OZ22W/bxF2eOxGzqTJrcgQeAoyLiakkbAwsknR8R11eJ/4XArwcrS9oZmAvsAmwJ/FDSThHx4BqI38zMhpi0WyYibomIq6uf7waWAFtVT38GeC8QtU0OBE6NiPsi4kZgKbB31qjNzGxCrfrcJc0G9gSulPRS4L8j4poxq20F3FR7vJyVB4N6W4dKmi9p/ooVK9pFbWZmE2qc3CVtBJwOHEnqqvkA8MFhqw5ZFqstiDghIuZExJyZM2c2DcPMzBpolNwlrUdK7KdExBnADsB2wDWSlgFbA1dLejzpTH1WbfOtgZtzBm1mZhNrMlpGwInAkog4DiAiFgGPra2zDJgTEbdJOgv4hqTjSBdUdwSuWgOxm5llUeKIqiajZfYFDgYWSVpYLTsmIs4dtnJELJZ0GnA9qfvmcI+UMTObWpMm94i4lOH96PV1Zo95PA+Y1ysyMzPrrMmZu61lSvyIaWbtFJHcuyQzJzJ7JPP/TPk8t4yZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBXJyNzMrkJO7mVmBnNzNzArk5G5mViAndzOzAjm5m5kVyMndzKxATu5mZgVycjczK1ARxTpsNLkilNn08Zm7mVmBnNzNzArk5G5mViD3uZtNIRemtqniM3czswI5uZuZFWjS5C5plqQLJS2RtFjSEdXyj0m6VtJCST+QtGVtm6MlLZV0g6T91uQvYGZmq2ty5v4AcFREPBnYBzhc0s7AJyNit4jYA/gu8EGA6rm5wC7A/sDnJa2zJoI3M7PhJk3uEXFLRFxd/Xw3sATYKiLuqq32KCCqnw8ETo2I+yLiRmApsHfesM3MbCKtRstImg3sCVxZPZ4HvAG4E3hutdpWwBW1zZZXy8a2dShwKMA222zTMmwzM5tI4wuqkjYCTgeOHJy1R8QHImIWcArwjsGqQzaP1RZEnBARcyJizsyZM9tHbmZm42qU3CWtR0rsp0TEGUNW+Qbwyurn5cCs2nNbAzf3CdLMzNppMlpGwInAkog4rrZ8x9pqLwV+Wv18FjBX0gaStgN2BK7KF7KZmU2mSZ/7vsDBwCJJC6tlxwCHSHoi8BDwK+BtABGxWNJpwPWkkTaHR8SDuQM3M7PxTZrcI+JShvejnzvBNvOAeT3iMjOzHnyHqplZgZzczcwK5ORuZlYgJ3czswJ5PveMPFe3mY0KJ3czmzYuor7muFvGzKxATu5mZgVycjczK5CTu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcjJ3cysQE7uZmYFcnI3MyuQk7uZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBXJyNzMr0KTJXdIsSRdKWiJpsaQjquWflPRTSddKOlPSprVtjpa0VNINkvZbg/GbmdkQTc7cHwCOiognA/sAh0vaGTgf2DUidgN+BhwNUD03F9gF2B/4vKR11kTwZmY23KTJPSJuiYirq5/vBpYAW0XEDyLigWq1K4Ctq58PBE6NiPsi4kZgKbB3/tDNzGw8rfrcJc0G9gSuHPPUW4Dzqp+3Am6qPbe8WmZmZlOkcXKXtBFwOnBkRNxVW/4BUtfNKYNFQzaPIe0dKmm+pPkrVqxoF7WZmU2oUXKXtB4psZ8SEWfUlr8ROAA4KCIGCXw5MKu2+dbAzWPbjIgTImJORMyZOXNm1/jNzGyIJqNlBJwILImI42rL9wfeB7w0Iu6pbXIWMFfSBpK2A3YErsobtpmZTWTdBuvsCxwMLJK0sFp2DHA8sAFwfsr/XBERb4uIxZJOA64nddccHhEPZo/czMzGNWlyj4hLGd6Pfu4E28wD5vWIy8zMevAdqmZmBXJyNzMrkJO7mVmBnNzNzArk5G5mViAndzOzAjm5m5kVyMndzKxATu5mZgVycjczK5CTu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcjJ3cysQE7uZmYFcnI3MyuQk7uZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBZo0uUuaJelCSUskLZZ0RLX8VdXjhyTNGbPN0ZKWSrpB0n5rKngzMxtu3QbrPAAcFRFXS9oYWCDpfOA64BXAv9ZXlrQzMBfYBdgS+KGknSLiwbyhm5nZeCY9c4+IWyLi6urnu4ElwFYRsSQibhiyyYHAqRFxX0TcCCwF9s4ZtJmZTaxVn7uk2cCewJUTrLYVcFPt8fJqmZmZTZHGyV3SRsDpwJERcddEqw5ZFkPaO1TSfEnzV6xY0TQMMzNroFFyl7QeKbGfEhFnTLL6cmBW7fHWwM1jV4qIEyJiTkTMmTlzZtN4zcysgSajZQScCCyJiOMatHkWMFfSBpK2A3YEruoXppmZtdFktMy+wMHAIkkLq2XHABsAnwNmAudIWhgR+0XEYkmnAdeTRtoc7pEyZmZTa9LkHhGXMrwfHeDMcbaZB8zrEZeZmfXgO1TNzArk5G5mViAndzOzAjm5m5kVyMndzKxATu5mZgVycjczK5CTu5lZgZzczcwK5ORuZlYgJ3czswI5uZuZFcjJ3cysQE7uZmYFcnI3MyuQk7uZWYGc3M3MCuTkbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBXJyNzMrkJO7mVmBnNzNzArk5G5mVqBJk7ukWZIulLRE0mJJR1TLN5d0vqSfV983q21ztKSlkm6QtN+a/AXMzGx1Tc7cHwCOiognA/sAh0vaGXg/cEFE7AhcUD2mem4usAuwP/B5SeusieDNzGy4SZN7RNwSEVdXP98NLAG2Ag4Evlqt9lXgZdXPBwKnRsR9EXEjsBTYO3PcZmY2AUVE85Wl2cAlwK7AryNi09pzt0fEZpL+GbgiIr5eLT8ROC8ivjWmrUOBQ6uHTwRu6PF7jGcL4LYRaGOUYnEbbmNticVtTG7biJg57Il1m7YgaSPgdODIiLhL0rirDlm22hEkIk4ATmj6+l1Imh8Rc6a7jVGKxW24jbUlFrfRT6PRMpLWIyX2UyLijGrxbyQ9oXr+CcBvq+XLgVm1zbcGbs4TrpmZNdFktIyAE4ElEXFc7amzgDdWP78R+E5t+VxJG0jaDtgRuCpfyGZmNpkm3TL7AgcDiyQtrJYdA3wcOE3SIcCvgVcBRMRiSacB15NG2hweEQ/mDryhHN0+ubqORiUWt+E2pqIdt5G/jVZaXVA1M7O1g+9QNTMrkJO7mVmBnNzNzArk5F5Tje6ZdpL2rb5vMM1xjMT+gJHaJ8dW3181zXGcXH0/ImOb60vaTdJTJK3fYrsLqu/H5oqli+l+bwyMzHu1xAuqkvYC/pR089SPBtMnNNhuQUQ8VdIFEfH8DHFsRhoKuuFgWURc0iKOqyNirwxxrNX7Y0wsvfeJpO2BzwLPAB4CLgfeFRG/bLDtImAv4MoMccwE3gfszKr75HkNtr0eeBFp6PFzGHPzYET8T8tYXgx8EfhF1dZ2wGERcV7DWN5ebf+6IbE0er/1NXhvSDo5Ig7O0N4rWPl/c2lEnNlwu6z/v101vkN1bSHpg6RhmYObrU6S9M2I+PsGm8+Q9CFgJ0l/O/bJMeP8J4vjrcARpJu4FpImXbscmPQfF7hf0knA1pKOHxLHO1vEUcL+gIz7BPgG8C/Ay6vHc4F/B57eYNvvkW4jf5Sku2rLlcKIR7eI4xTgP4AXA28j3S+youG2X6xi2R5YwKoJNarlbXwaeG5ELAWQtANwDjBpcgc+SJo4cOuqnbGxNP0bI2kf4HPAk4H1gXWA/224X9eX9EbgmVViXkXtBswmcXwe+BPS+wLgMEkviIjDG2ye873aWXHJHXgtsGdE/AFA0seBq4EmyWwuaQK0dYGNe8ZxBPA00jw7z5X0JOAjDbc9AHgB6Z9iQc84StgfkHefKCJOrj3+uqR3NNz2/0XEeyR9JyIO7BnHYyLiRElHRMTFwMWSLm647dkRcbykL0TE23vGAfDbQWKv/JKVd51P5paIeJGkD0bER3vG8c+k9903gTnAG0hJtom3AQcBmwIvGfNcsPIEp4lnA7tG1bUh6avAoobb5nyvdlZicl9G+oj7h+rxBqSPmk3sHxHHStogw5v0DxHxB0lU7f1U0hMbbvueiHifpG0i4quTrz6hZaz9+wPy7pMLJb0fOJX0T/8a4BxJm8OkXRqXk7pl7ppgnabur77fUnWL3Ew6+23iW8BTgZ0yxAGwWNK5wGmkffIq4MeDM+BJznqPr2J5GdD3fUJELJW0TnXz40mSLmu46RMi4u2SflLNXdXHDcA2wK+qx7OAaxtum/O92llxfe6Svk06Qzyf9CZ9IXAp1VnIRB+JJC2MiD0y9eueCbwZOJJ0BL8dWC8i/qLBtjn7db/NWr4/qu1z7pMbJ3g6ImLcLg1J1wGfJHVFvGfIxm0++h8A/BcpcXwOeDTwkYg4q8G2PwG+DbwV+MyQOBp3mVXtnTTB0xERb5lg2ytIU4H/BambaezGbboRLyGd9X4ZuBW4BXhTROzeYNtBn3uO9+vFpP+bwdQpTyMd2O8BiIiXTrBttvdqHyWeuZ9ZfQ1c1GLbJZKWATMl1Y/Sg/7U3Zo2FBGD/twPS7oQ2ITUR9rE2H5dkRJzl37dEvYHZNwnEdFnFFC2j/4R8d3qxzuB57aMI2eXGRHx5h6b5+yGOJjUz/4O4F2kA98rG277u+q9tZ2k1Q6QEyXkIT7YYt2xcv7/dlbcmXtfkh4PfB9Y7Y0QEb9afYtx29lm2PKI+HWLNnL06/YySvujaqfzPpH0vIj4z2EX26pY2px1HxIRJ3aJo9bGSQyfDnvcs+QhbbyoyYiWCbZ/b0R8QtLnxomlzVn37hFxTddY+lIavrkXcDLpE80qqusaUxnPtP7/FpPcJZ0WEa+uPhINe5M2Psus2luflf2ZN0TE/ROtP2T7QRwi9XlvV7WzS8t2Hkf6SAjpY16j0RSl7o+qra775CMR8aFxuiAm7HoY0tb6pLP4Z1WLLga+2Ga/SKqfkW5IGr1zc8uEugnwoTFxfDQi7my4/Usi4uxqlMlq2vQZS9qa1L20L9XwQeCIiFjeoo0bGf5+bTz6R9LMiFghaeO0afy+xbaXRsSfSrp7TBydzrq7vldzKCm5PyEibpG07bDnW55lPhv4GulipEgfDd8YDcdkj9PmXqRxw4e12OZVwKdIXSkC/ox0seZbE21XbVvc/qi267xPcpL0ZWA9VpaaPBh4MCJWO2Ns0eYM4IfRYJx7bZvTgevGxLF7RAz9dLImSTqfNMx0MBLp9cBBEfHCFm08pvZwQ9KF3c0jonE3iaRdqxg2J71HVpDer9c1bSOHaX+vRkSRX6SLU5sPvlpuuwB4Yu3xTsCCDDFd3XL9a4DH1h7PBK55pO6PXPuE1F/+TuA40kiP44Hj28bRZFnLNp8ILG25zcImyxq0M4d0beZq0qiQa4FrM+yT1rEMaePSlutfRhqzP3j8HOCyDq+7GbAbqatnL2Cvtvsj1/9vl6/iLqhKOow0HOteVn6santTx3oR8XBN14j4mVI1qjZx1G/6mUF6c7T9SDYjIupjjX9HyykjCtsfkGGfAOcCV5DGLT/UIQaAByXtEBG/gIfvem1Vt6D20X9wwe1W0h2rbdwr6U8j4tKqzX1Jf+u2TiGN/umzT1ZIej0rb/x5Lenv01j1iW5gBumg0/aC8aMi4sLBg4i4SNKjWsbxMeBNpPH+g/3R6oYs8rxXOysuuQPvBnaJiD7FaOcrFfYefLw8iPajAOpvyAdId/ud3rKN70n6Piv/WV5DSkxtlLQ/IM8+2TAiVrvjtqX3kMbL/5KUnLclDfVsLCJ6j3Ih9ft/rep7F/A/rKyQ1saKaDAEcxJvId2E9BlSIrysWtbGp2s/P0DqCnx1yzZ+KenvWLV7aKLhr8O8GtghIv6v5XZ1Od6rnRXT5z4g6XvAKyLinh5tbAAcTppXQsAlwOcj4r48UbaKZTC/hYBLouH8FrXti9ofVTx998m7gN8D3wUe/h2i/XwsG5C6UgT8tO3+GHOWuppoMSeLpEdX23S6uUrS80ln2hew6j5pc1fnSFCaw+gjrHyPXEy6f+D2Fm2cDrx9zJl3l1h6vVd7vXaByX1P4CTgSlZ9k3aaz2FwYbLDdhOeBUW7MbdI2qvNP3ttuyL3R9Vm131yODAPuINaV1W0GJExpr0PR8SHO2x3Bal76lrSP/9upL/T/VU8bboAkHRCRBzaNo5q268DTwIWU+uGiBYjiMa0992IOKDDdhN+ooqWN2d1JWkOqS70daz6f9P6fVq11+m92keJ3TL/Cvwn/foO684h/QO2dSPweODr1ePXkj5ifr9jHF/uGEep+wO675O/Bf6kZ1dV3UuBD3fYbhnwVxGxCB4e5fHuiHhTxzjmdNwO0gibp/TYfqytOm43hzR0cHAy8BLSJ8WbujTW427VrwLHku//put7tbMSk/sDGfpT6zT5KkPtGRHPqj0+W9IlEXHMFMdR6v7oE8tiqtvIM+kax5MGiR0gIq6TtEePOPp0IVwhaeeIuL5HG3U/6bjdFqRRKXdD+lQEfDO6DzHt+re5LSJWm9Gxh65xdFZicr9Q0qHA2fToT635UsftZkraPqo5wpUKX8zs2Ba0m0GxrtT9Ad33yYPAwupW9d5dVXQ/I1tSjZf/Oql76PWkOVo6iYj9u25L6hd+Y3UT0X2svGmny81uTwI+I2n9DhcktwHq2/wfMLtlG3XndNxugaR/JH2CqL9HunatdH2vdlZin/uwq+Kt+1PVscBFbfv9gRNIQ6kgvUEPi4hW3RDqUeCi2r6o/VG11Xef5Lgbc3s6FvyotbEhqcjF4BPNJcAXopqeuWEbnQt+jGknx81unQt+1Nr4AGmkypmk99rLgdMi4h+atpFDdeAfq8t1kMEF1aBFwY8cikvuOWj1AhcvI300bDIHer2dDUhnMdBtNMXQAhdt32B9jcr+qNoYlX1yBangx2CY21zgbyKiScGPnHH8gDQT47upFfyIiLbj5QftPZZVDxJt5kL6KXBAjCn4ERFPmnjL1dp5KikhQhph0qqLR/0KfmSj1Qt+vAb4RTQr+NH/9UtM7tWFqbFnMl9rsf0SVi1w8UekuymfPMVxLGJlgYs9VBW4iIjXTHEcI7E/qjZ67xNJOwL/OCSWNvOXXDk2kUu6IiL2meI4BiXdrh10oUi6OCKe3bSNapuXksaYb0nqu98WWBIt5v6prqE8q/ZYwMVjrrU0bavPQWY+Qwp+RMQHWsbwYmCXMXE0nq9e0mJWLfgxA1jUZp/2UVyfu1JZuOeQ/mHOJdWZvJQ0N0pTy+he4CJnHH0LXJS2PyDDPiENDf0Q6Wab55JuPmp7watPwY+ccfQp+FH3MdKnoB9GxJ6Snksa0dRGn4IfwNCDzDbAT0lJtrHoXvBjEMcXgT8m/V2+DPwlK+d2b6pPwY/epuxW2Cn0l8DzgVsjzVG9OykZtXEf6Y36FaUZBK8Dfi/peA2pibgG41guaVNSUYbzJX2H9M/bRkn7A/Lskz+KiAtIn1x/VY1Rb9ut8xrgMOBC0sRQbyfdjbkAmD+Fcfy90t2pR5G6Zr5Mmge9rfsj4nekurkzIt2+v0fLNjYEfkMqUfcc0vQSm5OGMzYd8z44yPws0rz7LwB+1DKOe6oLuwslfULpprVW0w8Az4yINwC3R8RHSNdWZrVs4zGki+YXSboIuJ40sOAsTXLfRw7FnbkD90bEQ5IeULpr77e0LxZ8Jt0LXGSLI/oXuMgSByOyPyDbPvlD9RH550q1U/8beGzLOPoU/MgZR5+CH3V3SNqIdFH3FEm/Jd3+3yaWPgU/Bu6PiN9JevggI+nYlm30KfgxMJif5x5JW5LmhWn7N+9T8KO3EpP7/OrM7kuks6jf0/LjVJtRE2syDq1a4GIw6uXxQJsCF8XsD8i2T44kfeR+J+lM8bk0nI9FGQt+9ImjFs9J9Cz4UTmQlNDeRZo7aBMa1kNVxoIf5DnIDLpB7qX7EMTvVu/XT5JmygxaDgOOKS4OMlaRF1QHJM0GHh0Rjfq5lLnARdc4attlK3DRJY5R2x/Vtrn3yeMj4tYW62cr+NEnjtp2vQt+DGnzgNongibr5yz48ShSUp7ByoPMKVWXUdM2bqRnwY8x7W1AmmyuaQGUrAU/uio9ubea80MZC1z0iWOCdjoVuOgax6jvj6qtvvukdzHlHHLFoQ4FP9ZULH21PcjUtutd8GNMe53n7JlOpSf3Xm/Sqm/44a6rhiMgsseRq60S90ff9iT9JCL27LDdpqQhdrNZdZ90nZCtUxxD2nkiaWz5n/Roo+s+mQN8gDSMsr5Pun7Cy/l/c2lE/Onka+aNQ+mGu1msuj+mZAKxEvvc6zrN56A8BS5yxJGrwEXfOEZif1Sx5N4nXadTyFHwo3ccylPwY6xOn4LIU/Cjruv7NUfBj7pOc/YoT8GPzko/c58REa3fZJJ+DjwjMs0a2COOD9UeDgoXnB4tbk/PFMdI7I9q2yz7pO8ZVcZulGk7s6tef+iF4VosjS8Q9zk7Hqe9vSOiy0X3+tQBg/fIp6JWTWwqSLoBeEr0K/jR/fVLS+5KE1L9Dat/XG48D7PyFLjoHUcO3h9DYxmcUf2CVedzb3xGpQwFPzLF0avgxzgXhmubN79ArB4FP3IeZPqQ9E8RcaSksxl+UbbN/02Wgh9dlZjcrwFOZMxHwzbDkpShwEWmOHoXuChpf1Tt5Ngnvc+olKHgR6Y4shb86EM9Cn5kPsh0Lvgh6akRsUDS0OkbWv7fZC340VaJfe5/iP7zMOcocJEjjhwFLkraH5Bnn1wHbEq/+c9zFPzIEccyMhT8kHQE6QB+N6n/fy/g/RHxgxbN7B4dC35EnhugBjoX/KgS+zqkffr6nnHkLvjRSoln7q8jTQf7AzrOwyzpsoh45gjEscpETOMtm4I4RmJ/VO3k2Ce9z6iqTxBze3ZV5YhjYUTsMdmyBu1cExG7S9qPVC/374CT2lxXkPQl4DPRo+BHjoOM0kyZr4yVBT82Js1i2ni+e6XC1i/p+amq9QRuOZV45v4U0u3Hz6P7FeocBS5yxJGjwEVJ+wPy7JMcZ1Q5Cn7kiCNXwY/ByJS/ICX1ayS1Ha2So+DHWyLis9VB5rGkydROIp0UNJWj4Mcy4EfVQfx/Bwsn6tIZInfBj1ZKTO4vB7bveYX6ddX3o2vL2g79yxHHu4CLJK1S4KJlGyXtD8izT3KUUPt29dVHjjjeTJq07Ijq8SXAFzq0s6A6490OOLo62217wOlTCWogx0HmZOAqSfWCH21nH725+hKwUcttBwb3CtSngfZQyK4k/QepaMK0XKHOHYf6F/woan9UbfXdJ8eRzqSm5Yxq1OKoYplBmgXylxFxh9JdnltFtyki+szFfhKpuPZ2pJlD1wEuiointoyhb8GPpwHHsOrorrafQqZVicn9ItKIgR/T4wq1+he4cBxrII5MsfQuoaY8hTZGIo5aW6+gR0k45Sn4MSoHmRtIUyhfx6qju1pNuaGeBT/6KLFb5kOTrzIx5Sks4Tgyx5ErlojoMzXuQO9CG6MSB4BWLwl3mKQXRLuScL0LfkSaFno28HpJuQ4yXQp+rIiIs9u87pA4chT86C4iivoiTco/o/p5J+ClwHot21hEum35murx44Cz17Y4qm2vy7BPRyKOjH+bx5HG3J9XPd4ZOKRlGwsG8dSW/dfaGEe1zWKqT/K1v9nilm3Mr75fU3vvX9Wyjc+TLp6+ufr6HvAvLdu4hlQo4yfV4+cCJ7Rs4/mkhPxa4BWDr5ZtXDvm+0bAD9r+bbp+zaA8lwAbStqKdKfcm4GvtGzj3ki3x/cpLDHtcVTbXqNV50DvYlTi6B1L5SukcfFbVo9/RppbvY1VCm1IejktC22MUBywsiTcQJeScHdo1bnYP0vLudhJVZz2i4iTIuIk0oXV57RsI0dVqTdX2+xPGiffpprUwNiCH/fTvuBHZyV2yygi7pF0CPC5SEUEFrZsI0dhiVGJ4wmkEnlXseqQrjZ93aMSR65YtoiI0yQdXcXwgKQHW7ZxJD0LbYxQHLCyJNxgXz4NuLwaCtj073QgHQt+1OSoOzr2INO64Ac9bsiq6V3wo48SL6j+BPhrUh/kIRGxWNKirn8odS+0MSpx9L6NehTj6BnLRaSya+dHxF6S9gGOjY43nKh7oY2RiKPadsLXbPt3Uve52C8mHVhWOcgA91RxNJleIkfBj943ZI1pr1XBjyymqv9nqr5IH+vOAt5XPd4eOL5Hex9em+NYA/t3JOLo+bfZi1R0+c7q+8+A3XrEcfXaHMca+tt03SfPnuirQ3sHdIxjCenmpxtInxwWUfWdd2yvVZ9/jq/iztwHJD0qIv538jUnbadvgYtpjaMa3nYsqR9W0K/U16jE0SeWatt1gSdWcdwQEff3iOMn0bHQxnTHoTVUEq7PPsmpx/t122HLo3v1sSmvblVcn7ukZ5BGIGwEbCNpd1IZtr/u2uTaHAfwCdIcGV1uSR/lODrHotWnl91J0p2kESddbrLqWmhj2uOIav71iOhTzGKYVncNr6mDDB3fI12T+ASm/CbC4s7cJV1JGk961uDMQdJ1EbFrx/a6FrgYlTh+FBH7dnnNUY6jZyznAM8ABjcRPYdUVWkn4KMRcXLDdvoW/BiJOPoacpBaRUzRXOzDqGPBjxIUd+YOEBE3adXpKFqNQNCYwhKDtqLlyI4RiWO+0q3/36ZlAYVRiyNjLA8BT46I31RtPo40H8vTSSMsJk2qGqfQBu3mDRmVOPp6yQTPBTBVhTaGHmQkbQ1rZ8GPPkpM7jdJeiYQktYnDRFr2xXwbVKXytl0n61vVOJ4NGmkwZ/XlrX9hxuVOHLFMnuQUCu/BXaKiP+R1LTP+9XADtFvIrRRiaOXyDsXex8jcZBh5UH5U1P0ekOVmNzfBnyWNPnQctLdbm37uXMUlhiVOGYAR0TEHfDwR/hPr6Vx5IrlvyR9F/hm9fiVwCXVELo7GraRo9DGqMSRhfIU/OhsVA4ykbfgR2cl9rl/FTgyIm6vHm8GfDralenKUeBiVOJYbdRC25EMoxJHxlhESqT7ki64XUoqst34n0F5Cm2MRBy5KEPBj0xxTOtBphZH74IffZR45r7bIKECRMTtSjVA28hRWGJU4pghabPaQWZz2v/dRyWOLLFUyfNb1VdXvQttjEocGeWYiz2HHAU/clhG/4IfnZWY3HMkkRyFJUYljk8Dl0n6FikJvppU2HltjKNXLEOG2a2i5XC7zoU2RiWONSBHwY8cRuUgk6PgR2clJvccSeQa+vdjjkQcEfE1SfNJZ7YizWzX9pbqUYmjVyyDsdySPgrcSrrwJdIt6m3HeXcuoTYqcawBh7ByLvZ7lOZin45+8FE5yJzLkIIftJ9vp5Pi+twBJO3MyiRyQdskonwFLkYijr5GJY5csUi6MiKePtmySdq4cMjiiHaFNkYijpzUs+BHphiyFfzoGUeWgh9dlXjmTpVE+0z486GS4shgVOKAPLE8KOkg4FRSEnotLe9BiDyFNkYljiyUp+BHb5Gh4EcmvQt+9FHkmXtf1VC0e6s3yU6kep3nRY95P9bWOKqzoGuj4521oxhL9Y//WdIolSBN2nVkRCxr0cbjgH8AtoyIF1Wf0p4RESeubXHkImkxsOtgtE/191oULcrsZYpj7EHmNcAvpvogI+n5pAP2BfS4ca/z6zu5r07SAuDPgM1It4PPB+6JiIMeoXGcAhwdLWpQlh6LpPNIIzA+UA3/W5dU+afvHOBrZRxVLGcA7xp0OyhNvvXxiGhVai9DHKNykPk66YRsMbWRXW2GQ/dRZLdMBjkKbZQUR65CGyMRi6QNSRf/xhYubvNP17vQxqjEkVGOgh855Cj4kcPu03GQHXByH05KszoeRPrnA1jnERzHR6bhNceTI5aTSQWT9yONXDiI9lND/G91oW5wdrgPaV72tTGOXD44Ta871qgcZK6QtHPHUWG9uVtmCKWqNEcBP4qIYyVtT+oLfecjMY7SDO6MlXRtROwmaT3g+y1HuuwFfA7YlTQaYibwl21GZIxKHKVR5qpSPeJYAuwA3Ejqcx9MX7zblLy+k/v4lKnQxtoeh9ZAoY3pjEXSVRGxt6RLSPP93ApcFRGtCm2rZ6GNUYmjL625udjXaspc8KP16zu5r061QhsRkaPQxtoex1LyF9qYtlgkvRU4nTSVwVdIdw/+XUT8a4s2hk0v26rQxqjEUQofZFbl5D6EMhfaKCCO7IU2usoRi1Kx4leS7hxcr1ocEdH4zkFlKLQxKnFYmXxBdRzRs9BGYXFkKbQxQrF8h3R2u6DeRku9C22MUBxWICf34XIU2igpjlyFNkYllq0jYv+eceQotDEqcViBnNyHy1Foo6Q4chXaGJVYLpP0lIhY1COOHIU2RiUOK5D73IdQhkIbhcXxk8hQaGO6Y5G0iHSWvy6p4Mcv6ThETepeaGNU4rCy+cx9uByFNkqKI1ehjemO5YBcQVTJs2uhjVGJwwrm5D7cqCSzUYkjV6GNaY0lx/jiIcPsxr7GpMPtRiUOK5uT+3CjksxGIo7IV2hjrY8l8hbaWOvjsNHlPvdxqGehjdLisFUpQ6GNkuKw0eMz93FE/0IbRcVhq+ldaKOwOGzEzJjuAMzWUq8jdZP9pvp6VbXskRqHjRh3y5iZFcjdMmYdKE+hjWLisNHjbhmzbk4GHk8qtHExsDVw9yM4Dhsx7pYx6yBHoY2S4rDR4zN3s24Gk3LdIWlXYBPS1L2P1DhsxLjP3aybE6q5fv4fcBZVoY1HcBw2YtwtY9ZBjkIbJcVho8dn7mbd5Ci0UVIcNmJ85m7WwXSUOxzlOGz0+IKqWTeXSXrKdAfB6MRhI8Zn7mYt5Cy0UUIcNrqc3M1akLTtRM/nmKt9bYrDRpeTu5lZgdznbmZWICd3M7MCObmbmRXIyd3MrEBO7mZmBfr/nVIrRYf7WwoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the δOLR\n",
    "# use rlut for control, but rlutcs for lwoff\n",
    "olr_correct = {}\n",
    "for i in rlut:\n",
    "    if 'off' in i:\n",
    "        olr_correct[i] = rlutcs[i]\n",
    "    else:\n",
    "        olr_correct[i] = rlut[i]\n",
    "fig, ax = plt.subplots()\n",
    "bot=200\n",
    "ax.bar(np.arange(len(olr_correct)), np.array(list(olr_correct.values()))-bot, bottom=bot)\n",
    "ax.set_xticks(np.arange(len(olr_correct)))\n",
    "ax.set_xticklabels(olr_correct.keys(), rotation=90)\n",
    "ax.set_title(\"OLR, global & time average\")\n",
    "fig.savefig(\"/Users/brianpm/Dropbox/Projects/extremes/cmip_lwoff_OLRcorrect_bars.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note:\n",
    "# Some models (CESM2 for sure, likely MRI) have diagnostic LW fluxes that **do** include clouds,\n",
    "# so you should use the clear-sky fluxes for analysis/budgets.\n",
    "# BUT you need to use the clear-sky fluxes only for the lwoff experiments!\n",
    "VARIABLES = [\"rsdt\", \"rsut\", \"rlut\", \"rlutcs\", \"rsds\", \"rsus\", \"rlds\", \"rldscs\", \"rlus\", \"hfss\", 'hfls', 'pr', 'prw']\n",
    "\n",
    "# MODEL = \"MRI-ESM2-0\" # \"IPSL-CM6A-LR\" # \"CNRM-CM6-1\" # \"CESM2\" # \"CESM2\"  #  # \n",
    "CASE0 = \"amip\" # \"aqua-control\" #  \n",
    "CASE1 = \"amip-lwoff\" # \"aqua-control-lwoff\" #  \n",
    "# REALIZATION = \"r1i1p1f1\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CESM2\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rsdt' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rsut' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlut' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlutcs' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rsds' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rsus' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlds' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rldscs' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'rlus' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'hfss' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'hfls' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n",
      "/Users/brianpm/miniconda3/lib/python3.7/site-packages/xarray/conventions.py:498: SerializationWarning: variable 'prw' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  decode_timedelta=decode_timedelta,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Surface imbalance: (control) -3.962, (perturbed)  8.788, difference = 12.750\n",
      "-12.668118990253674 + 16.995954466865214 + -0.5517143161004014 = 3.7761211605111384\n",
      "δP = -3.6368985116724284\n",
      "δE = -3.6938362920707704\n",
      "δPRW = -1.1652383888711597\n",
      "CESM2 & -3.6 & -3.7 & -12.7 & 17.0 & -0.6 & -3.8\n",
      "CNRM-CM6-1\n",
      "Surface imbalance: (control)  0.253, (perturbed) 15.752, difference = 15.498\n",
      "-15.479893088902116 + 19.84039397308817 + 0.07609687212636018 = 4.436597756312413\n",
      "δP = -4.362056636814254\n",
      "δE = -4.418114290675007\n",
      "δPRW = -1.3950838162766672\n",
      "CNRM-CM6-1 & -4.4 & -4.4 & -15.5 & 19.8 & 0.1 & -4.4\n",
      "HadGEM3-GC31-LL\n",
      "Surface imbalance: (control) -0.905, (perturbed) 12.215, difference = 13.120\n",
      "-13.072275774673017 + 17.18274474505199 + -0.754578423765917 = 3.3558905466130575\n",
      "δP = -3.238369789511023\n",
      "δE = -3.3081098119939725\n",
      "δPRW = -1.4990674513702231\n",
      "HadGEM3-GC31-LL & -3.2 & -3.3 & -13.1 & 17.2 & -0.8 & -3.4\n",
      "IPSL-CM6A-LR\n",
      "Surface imbalance: (control) -2.931, (perturbed)  9.568, difference = 12.499\n",
      "-12.512710571289062 + 18.1861572265625 + -1.408651351928711 = 4.264795303344727\n",
      "δP = -4.225937843322754\n",
      "δE = -4.278312683105469\n",
      "δPRW = -1.5664501190185547\n",
      "IPSL-CM6A-LR & -4.2 & -4.3 & -12.5 & 18.2 & -1.4 & -4.3\n",
      "MRI-ESM2-0\n",
      "Surface imbalance: (control) -2.082, (perturbed) 13.852, difference = 15.933\n",
      "-15.908681658484 + 21.150947889806275 + -1.7477667650959319 = 3.494499466226344\n",
      "δP = -3.3736416212680544\n",
      "δE = -3.469716677272501\n",
      "δPRW = -2.028473406871921\n",
      "MRI-ESM2-0 & -3.4 & -3.5 & -15.9 & 21.2 & -1.7 & -3.5\n"
     ]
    }
   ],
   "source": [
    "L = 2.501e6\n",
    "\n",
    "# hold on to global avgs for plots\n",
    "toa_sw = {}\n",
    "toa_lw = {}\n",
    "sfc_sw = {}\n",
    "sfc_lw = {}\n",
    "sfc_sh = {}\n",
    "sfc_lh = {}\n",
    "\n",
    "for MODEL in [\"CESM2\", \"CNRM-CM6-1\", \"HadGEM3-GC31-LL\", \"IPSL-CM6A-LR\",  \"MRI-ESM2-0\" ]:\n",
    "    print(MODEL)\n",
    "    if MODEL == \"CNRM-CM6-1\":\n",
    "        REALIZATION0 = REALIZATION1 = \"r1i1p1f2\"\n",
    "    elif MODEL == \"HadGEM3-GC31-LL\":\n",
    "        if CASE0 == \"amip\":\n",
    "            REALIZATION0 = \"r1i1p1f3\"\n",
    "            REALIZATION1 = \"r5i1p1f2\"\n",
    "        elif CASE0 == \"aqua-control\":\n",
    "            REALIZATION0 = REALIZATION1 = \"r1i1p1f3\"            \n",
    "    else:\n",
    "        REALIZATION0 = REALIZATION1 = \"r1i1p1f1\"\n",
    "    data0 = {}\n",
    "    data1 = {}\n",
    "    has_variables = True\n",
    "    for v in VARIABLES:\n",
    "        data0[v] = cmip_loader(dloc, v, \"Amon\", MODEL, CASE0, REALIZATION0)\n",
    "        data1[v] = cmip_loader(dloc, v, \"Amon\", MODEL, CASE1, REALIZATION1)\n",
    "        if data0[v] is None:\n",
    "            data0.pop(v)\n",
    "            has_variables = False\n",
    "        if data1[v] is None:\n",
    "            data1.pop(v)\n",
    "            has_variables = False\n",
    "    if not has_variables:\n",
    "        continue\n",
    "\n",
    "    \n",
    "    glb0 = {}\n",
    "    glb1 = {}\n",
    "    for i in data0:\n",
    "        w = np.cos(np.radians(data0[i]['lat']))\n",
    "        glb0[i] = esmlab.weighted_mean(data0[i].mean(dim=('time','lon')), dim=[\"lat\"], weights=w)\n",
    "    for i in data1:\n",
    "        w = np.cos(np.radians(data1[i]['lat']))\n",
    "        glb1[i] = esmlab.weighted_mean(data1[i].mean(dim=('time','lon')), dim=[\"lat\"], weights=w)\n",
    "\n",
    "    # L δP = δRtoa - δRsfc - δH\n",
    "\n",
    "    # Atmospheric energy => positive fluxes are \"into\" the atmosphere\n",
    "    # L δP = -δRtoa - δRsfc - δH \n",
    "    \n",
    "    # use clear-sky for lwoff experiments ... (*1)\n",
    "    drtoa = (glb1['rsdt'] - glb1['rsut'] - glb1['rlutcs']) - \\\n",
    "            (glb0['rsdt'] - glb0['rsut'] - glb0['rlut'])\n",
    "    dh = glb1['hfss'] - glb0['hfss']\n",
    "    dpr = L*(glb1['pr'] - glb0['pr'])\n",
    "    dlhflx = (glb1['hfls'] - glb0['hfls'])\n",
    "    # I think that for AMIP/AQUA, the fixed SST/ice means surface imbalance is probably not as constrained, so look:\n",
    "    rsfc0 = (-glb0['rsds'] - glb0['rlds'] + glb0['rsus'] + glb0['rlus'])\n",
    "    rsfc1 = (-glb1['rsds'] - glb1['rldscs'] + glb1['rsus'] + glb1['rlus'])\n",
    "    drsfc = rsfc1 - rsfc0\n",
    "    ressurf0 = rsfc0 + glb0['hfls'] + glb0['hfss']\n",
    "    ressurf1 = rsfc1 + glb1['hfls'] + glb1['hfss']\n",
    "    print(f\"Surface imbalance: (control) {ressurf0.values:6.3f}, (perturbed) {ressurf1.values:6.3f}, difference = {(ressurf1-ressurf0).values:6.3f}\")\n",
    "    print(f\"{drtoa.values} + {drsfc.values} + {dh.values} = {(drtoa + drsfc + dh).values}\")\n",
    "    print(f\"δP = {dpr.values}\")\n",
    "    print(f\"δE = {dlhflx.values}\")\n",
    "    print(f\"δPRW = {(glb1['prw'] - glb0['prw']).values}\")\n",
    "    print(f\"{MODEL} & {dpr.values:3.1f} & {dlhflx.values:3.1f} & {drtoa.values:3.1f} & {drsfc.values:3.1f} & {dh.values:3.1f} & {(-drtoa - drsfc - dh).values:3.1f}\")\n",
    "    \n",
    "    toa_sw[\"_\".join([MODEL,CASE0])] = (glb0['rsdt'] - glb0['rsut']).values.item()\n",
    "    toa_sw[\"_\".join([MODEL,CASE1])] = (glb1['rsdt'] - glb1['rsut']).values.item() # net upward sw flux\n",
    "    toa_lw[\"_\".join([MODEL,CASE0])] = (glb0['rlut']).values.item() # upward lw flux\n",
    "    toa_lw[\"_\".join([MODEL,CASE1])] = (glb1['rlutcs']).values.item() # upward lw flux **use clearsky for LWOFF\n",
    "    sfc_sw[\"_\".join([MODEL,CASE0])] = (glb0['rsus'] - glb0['rsds']).values.item() # net upward sw flux\n",
    "    sfc_sw[\"_\".join([MODEL,CASE1])] = (glb1['rsus'] - glb1['rsds']).values.item()\n",
    "    sfc_lw[\"_\".join([MODEL,CASE0])] = (glb0['rlus'] - glb0['rlds']).values.item() # net upward lw flux\n",
    "    sfc_lw[\"_\".join([MODEL,CASE1])] = (glb1['rlus'] - glb1['rldscs']).values.item()\n",
    "    sfc_sh[\"_\".join([MODEL,CASE0])] = glb0['hfss'].values.item()\n",
    "    sfc_sh[\"_\".join([MODEL,CASE1])] = glb1['hfss'].values.item()\n",
    "    sfc_lh[\"_\".join([MODEL,CASE0])] = glb0['hfls'].values.item()\n",
    "    sfc_lh[\"_\".join([MODEL,CASE1])] = glb1['hfls'].values.item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cesm-amip': 238.0526155592835,\n",
       " 'cesm-amip-off': 235.43461851563595,\n",
       " 'cesm-aqua': 227.29511520788338,\n",
       " 'cesm-aqua-off': 225.07610839946426,\n",
       " 'cnrm-amip': 238.1980624145999,\n",
       " 'cnrm-amip-off': 256.3543160735536,\n",
       " 'cnrm-aqua': 233.02272962854997,\n",
       " 'cnrm-aqua-off': 256.19433943648625,\n",
       " 'hadgem-amip': 240.78337220838014,\n",
       " 'hadgem-amip-off': 260.44294916634226,\n",
       " 'hadgem-aqua': 231.98320319818112,\n",
       " 'hadgem-aqua-off': 259.30868889414154,\n",
       " 'ipsl-amip': 237.94119262695312,\n",
       " 'ipsl-amip-off': 260.0086975097656,\n",
       " 'ipsl-aqua': 227.0421905517578,\n",
       " 'ipsl-aqua-off': 260.163818359375,\n",
       " 'mri-amip': 239.09466166534497,\n",
       " 'mri-amip-off': 240.00644632818089}"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rlut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cesm-amip': 262.14375951216164,\n",
       " 'cesm-amip-off': 258.24031363881045,\n",
       " 'cesm-aqua': 259.2399412574777,\n",
       " 'cesm-aqua-off': 255.86578116259952,\n",
       " 'cnrm-amip': 260.38499190240776,\n",
       " 'cnrm-amip-off': 256.3543160735536,\n",
       " 'cnrm-aqua': 257.8907395343314,\n",
       " 'cnrm-aqua-off': 256.19433943648625,\n",
       " 'hadgem-amip': 263.4922761850152,\n",
       " 'hadgem-amip-off': 260.44294589139076,\n",
       " 'hadgem-aqua': 260.9797557150099,\n",
       " 'hadgem-aqua-off': 259.3086892021664,\n",
       " 'ipsl-amip': 262.8774719238281,\n",
       " 'ipsl-amip-off': 260.0086975097656,\n",
       " 'ipsl-aqua': 260.50408935546875,\n",
       " 'ipsl-aqua-off': 260.163818359375,\n",
       " 'mri-amip': 263.11159822657555,\n",
       " 'mri-amip-off': 260.09405285754366}"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rlutcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'SFC Balance')"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAGXCAYAAABbU1CIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABFG0lEQVR4nO2debhcVZW33x9BCPMYBhMgzJMySJgEFEE60IxONIMKTtiKMkg3iK0iKNjOINL6pVsUFVFQkUEkzChKAgkgGAZBmQcJYRYCBn7fH3tXbqVS9+aG3H1O5Z71Pk89t84+p2qtfc+pdfZZe+21ZJsgCIKgWSxStwJBEARB9YTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IOgIJK+IOkndesRBJ2E8Q8WGiQ93/Z6VdKLbdsH52M2kXShpGckPSfpaklv7vJdO0uypGPnIXPnLKsl52FJJ5bqYxBURRj/YKHB9tKtF/AAsHdb29mS1gX+ANwGrA28HjgfuEzS9h1fdwjwZP47Lx5pk7sj8CFJ+w1Rt4KgFsL4B8OJLwDX2/4v20/afs72t4EfA19pHSRpSeDdwOHA+pLGDVaA7XuBPwKbtH3faZIelPSspKmSdurv85LOk/RYfjL5naRN2/b9UNIZkn6Tn1om5xtaa/+mki6X9KSkv0v6TG5fRNKnJf1V0gxJ50pacbB9CppJGP9gOLEbcF6X9nOBHbLRB3gX8Hw+diLw/sEKkLQ+sAMwqa35RmALYEXgp8B5kkb28xW/BdYHVgFuAs7u2H8gcCKwAnAPcHKWuwxwBXAp6YlmPeDK/JkjgP2At+Z9TwFnDLZPQTMJ4x8MJ1YGHu3S/ijpWl8hbx8C/Nz2KyRjfaCk1w3wva+X9LSkZ4G/AJOB61o7bf/E9gzbs2x/A1gc2LDbF9k+Mz+RvER6Utlc0nJth/zK9g22Z5FuDFvk9r2Ax2x/w/bM/B2T876PAv9l+6G27323pEUH6FPQcML4B8OJJ4DVu7SvDrwKPCVpDeBt9I24LwBGAnsO8L2P2F7e9rLA8sCLwFmtnZKOkXRHduU8DSxHuhHNgaQRkv47u2eeBe7Lu9qPfazt/QvA0vn9GsBf+9FvLeD8fIN6GrgDeAVYdYA+BQ0njH8wnLgCeE+X9v1JcwEvAO8jXfcXSXoM+BvJ+A/K9WP7GdLTwt4A2b9/XJaxgu3lgWcAdfn4QcC+wNtJN4ixub3bsZ08CKw7wL498g2q9Rpp++HB9CloJmH8g+HEicCbJZ0saUVJy0j6JMmwH5ePeX8+bou217uAPSWtNC8BkpYGDgCm5aZlgFnAdGBRSZ8Hlu3n48sALwEzgCWBU+ajbxcDq0k6StLiuW/b5n3fA06WtFbWcZSkfefju4MGEsY/GDbYvpsUirk5yaXyKMmwj7f9B0nbkUbbZ9h+rO11IWly9cB+vvr1rTh/4H7SxO7Bed9E0iTuX/K+maSReDd+lI95GLidOSeN59W350gT2nuTXEN3k9xXAKcBF5JCWp/L37ttt+8JghaKYi5BEATNI0b+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDWShWf698sore+zYsXWrEQRBsFAxderUJ2yP6mxfaIz/2LFjmTJlSt1qBEEQLFRIur9be7h9giAIGkgY/yAIggay0Lh9Fla+dflfiss4ercNissIgoGI63zhI4x/EAwhYQSbxcJ8vhth/BfmE7Qg1Nnvpsquk+h3OXqx3wtK+PyDIAgaSBj/IAiCBrLAxl/SGpKuzmXspkk6MrevKOlySXfnvyu0feZ4SfdIukvS+AXVIQiCIJg/hmLkPws4xvbGwHbA4ZI2AT4NXGl7feDKvE3edwCwKbA78D+SRgyBHkEQBMEgWWDjb/tR2zfl98+RikePJtUqbRW5PgvYL7/fF/iZ7Zds30uqoLTNguoRBEEQDJ4h9flLGgtsCUwGVrX9KKQbBLBKPmw0c5a5eyi3BUEQBBUxZMY/F7b+JXCU7WcHOrRLW9dakpIOkzRF0pTp06cPhZpBEAQBQ2T8Jb2OZPjPtv2r3Px3Savn/asDj+f2h4A12j4+Bnik2/fanmB7nO1xo0bNlZQuCIIgeI0MRbSPgO8Dd9j+ZtuuC4FD8vtDgAva2g+QtLiktYH1gRsWVI8gCIJg8AzFCt8dgPcBt0m6Jbd9Bvhv4FxJHwIeAN4DYHuapHOB20mRQofbfmUI9AiCIAgGyQIbf9vX0d2PD7BrP585GTh5QWUHQRAEr41Y4RsEQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvjPgeexHQRBMDwI45+55ppruHTiRPoMvrl04kSuueaaGrWqiibf9Jra9+h39+3mUJvxl7S7pLsk3SPp03XpkTAzX5rJ5EmTZ98ALp04kcmTJjPzpZkM5wukyTe9pvY9+t2sfvdHLcZf0gjgDGAPYBPgQEmb1KFL1ojdx49n2+22ZfKkyZx44klMnjSZbbfblt3HjwdUgQ51jEiae9Prjb439Zw3td+9xaI1yd0GuMf23wAk/QzYF7i9Jn1o3QAmT5o8u6Uqw3/NNdcw86WZbfLShTly8ZHsvPPOBSUry4TJkybP7nu1Nz1IPzwNsF2Cevve1HPe1H4n6rjO+6cu4z8aeLBt+yFg286DJB0GHAaw5pprvmZhRy/6y3keY8zEiRM5atE+47/xlQ8xfvx4NKgTdPxrkm3Mxq9MZNKUyWy86C8ZP348EydO5Jkpk9luu20Zv+gTg5D/2mS3OGpPc+KUvmNP2HNTxK8G9dkFlX3Ntdcwc+bM2f/n1nkYOXIkO79156KyYUH6/tplL/g5j37Pr+wWdfQbylznC0pdxr/bGZ7rucv2BGACwLhx44o9l7VOxKRJ+SLMF+WkPDoY/A1g/hFifB6RTJo0ebbMlh6l5LZo9b2diRMnViZ75syZc/yf28+DcVEd6up7U895U/td93XeH3UZ/4eANdq2xwCP1KQLQowcOXKOi7B1kY4cObL4iWnJm9Tmcqryx1DHTQ/qNQa90vemnvMm9bvum15/1GX8bwTWl7Q28DBwAHBQTboAsPNbd57jDtw6YVWcmDpHYnXe9Fo61GEM6u57U895U/td13U+ELUYf9uzJH0CmAiMAM60Pa0OXdrpPBFVGv66RmJ13vSgXrdTXX1v6jlvar+h3uu8P+oa+WP7EuCSuuT3CnWPSFo6DLRdirqNAdTT96ae86b2uxeu827UZvyDPuoefddFLxiDumjqOW9iv3v1Og/j3yPUNfqumyYagxZNPedN7HcvXueR2yeonSYag6B59Np1HsY/CIKggYTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IMgCBpIGP8gCIIGEsY/CIKggYTxD4IgaCBh/IMgCBrIAhl/SV+TdKekWyWdL2n5tn3HS7pH0l2Sxre1byXptrzv25K0IDoEQRAE88+CjvwvB95gezPgL8DxAJI2AQ4ANgV2B/5H0oj8me8ChwHr59fuC6hDEARBMJ8skPG3fZntWXlzEjAmv98X+Jntl2zfC9wDbCNpdWBZ29fbNvAjYL8F0SEIgiCYf4bS5/9B4Lf5/WjgwbZ9D+W20fl9Z3tXJB0maYqkKdOnTx9CVYMgCJrNovM6QNIVwGpddv2X7QvyMf8FzALObn2sy/EeoL0rticAEwDGjRvX73FBEATB/DFP42/77QPtl3QIsBewa3blQBrRr9F22Bjgkdw+pkt7EARBUCELGu2zO3AcsI/tF9p2XQgcIGlxSWuTJnZvsP0o8Jyk7XKUz/uBCxZEhyAIgmD+mefIfx58B1gcuDxHbE6y/e+2p0k6F7id5A463PYr+TMfA34ILEGaI/jtXN8aBEEQFGWBjL/t9QbYdzJwcpf2KcAbFkRuEARBsGDECt8gCIIGsqBun4WDtx1ftwb1EP1uFtHvYD6IkX8QBEEDCeMfBEHQQJrh9gmqJx7Fg6CnUd+6rN5m3LhxnjJlSt1qBEEQLFRImmp7XGd7uH2CIAgaSBj/IAiCBhLGPwiCoIEsND5/SdOB+ysStzLwREWyQnZvyA/ZIXu4yl7L9qjOxoXG+FeJpCndJkhC9vCVH7JDdhNktxNunyAIggYSxj8IgqCBhPHvzoSQ3Tj5ITtkN0H2bMLnHwSFkPRD4CHbn61blyDoJEb+wUKHpB0l/VHSM5KelPQHSVvnfYdKekXS822v77R9dhtJl0h6On/2Bkkf6EdO53f9TdLHqupnEJQkjH+wUCFpWeBi4HRgRWA0cCLwUtth19teuu31ifzZ7YGrgGuB9YCVSJXl9hhA5OzvAt4NfFXSlkPdryComjD+wcLGBgC2z7H9iu0XbV9m+9ZBfPZrwFm2v2L7CSem2t5/MIJt3wTcAWzcapN0nqTH8lPI7yRt2u2zklaQdLGk6ZKeyu/HtO2/RtIX81PMc5Iuk7Ry2/7W087Tkh6UdGhuX1zS1yU9IOnvkr4naYnB9CdoNmH8g4WNvwCvSDpL0h6SVhjMhyQtCWwP/OK1Cs6upQ2A9gyDvwXWB1YBbgLO7ufjiwA/ANYC1gReJNXAbucg4AP5uxYD/iPLXTPLOR0YBWwB3JI/85Ws0xakp5nRwOdfS/+CZhHGP1iosP0ssCNg4H+B6ZIulLRq22Hb5RFy67UdsALpen90PkW2vut54Abgx8Ddbfqcafs52y8BXwA2l7RcF71n2P6l7RdsP0eqb/3WjsN+YPsvtl8EziUZdICDgSvy084/83fdIknAR4CjbT+Zv/cU4ID57GPQQML4Bwsdtu+wfajtMcAbgNcDp7YdMsn28m2vScBTwKvA6vMprvVdSwOrAZuSDCySRkj6b0l/lfQscF/+zMqdXyJpSUn/T9L9+djfActLGtF22GNt718Als7v1wD+2kW3UcCSwNTWjQ64NLcHwYCE8Q8WamzfCfyQdBMY6LgXgOuBdy2ArL8DvwT2zk0HAfsCbweWA8bmdnX5+DHAhsC2tpcF3jLAsZ08CKzbpf0Jkvto07Yb3XL5RhUEAxLGP1iokLSRpGNak6WS1gAOBCYN4uPHAodK+k9JK+XPby7pZ4OUvRLwDmBablqGFGU0gzQCP2WAjy9DMtRPS1oROGEwMjNnA2+XtL+kRSWtJGkL26+SXF/fkrRK1nG0pPHz8d1BQwnjHyxsPAdsC0yW9A+S0f8zaWQ9ILb/COySX3+T9CRpteUlA3xs+1acPynSZzrwybzvR6RMsw8DtzPwDehUYAnSaH0SyT0zKGw/APwrqY9PkiZ7N8+7jwPuASZld9IVpCeMIBiQWOEbBEHQQGLkHwRB0EDC+AdBEDSQBTb+ktaQdLWkOyRNk3Rkbl9R0uWS7s5/V2j7zPGS7pF0V0xOBUEQVM8C+/wlrQ6sbvsmScsAU4H9gEOBJ23/t6RPAyvYPk7SJsA5wDak+OwrgA1sv7JAigRBEASDZtEF/QLbj5JXTdp+TtIdpCXm+wI758POAq4hRSbsC/wsr4i8V9I9pBvB9QPJWXnllT127NgFVTcIgqBRTJ069YluNXwX2Pi3I2kssCUwGVg13xiw/WgrDpl0Y2gPiXsotw3I2LFjmTJlyrwOC4IgCNqQdH+39iGb8JW0NGn141E5/0q/h3Zp6+p7knSYpCmSpkyfPn0o1AyCIAgYIuMv6XUkw3+27V/l5r/n+YDWvMDjuf0hUq6SFmOAR7p9r+0JtsfZHjdqVKQrCYIgGCqGItpHwPeBO2x/s23XhcAh+f0hwAVt7QfkPORrk9Lh3rCgegRBEASDZyh8/jsA7wNuk3RLbvsM8N/AuZI+BDwAvAfA9jRJ55KWw88CDo9In+HHHRttPO+DFpCN77yjuIz5Jfpdjl7s98LMUET7XEf/mQl37eczJ5PymQ974kcRBMOXhfn3HSt8gyAIGkgY/yAIggYypHH+QW+xMD+SBkFQljD+wbAjbnpBMG8aYfzDGARNoKnXeVP7vaCEzz8IgqCBhPEPgiBoIGH8gyAIGkgY/yAIggYSxj8IgqCBhPEPgiBoIGH8gyAIGkgY/yAIggYSxj8IgqCBhPEPgiBoIGH8gyAIGkgY/yAIggYSxj8IgqCBNCKrZ53sf3z5f/FtxSXMP03tdxAsLITxD4IhpKk3veh3OUr1O9w+QRAEDSSMfxAEQQMJ4x8EQdBAwucfBMECc9uaB9WtQi0szP0O41+YOi+OOiejFuYfxYLQ1H4HCx+NMP4L84x8MP/E+Q6CedMI4x8ETaCpN72m9ntBiQnfIAiCBhLGPwiCoIGE8Q+CIGggjfD5NzUCI/odNIE436+NRhj/IGgCYQSD+SHcPkEQBA0kjH8QBEEDCeMfBEHQQML4B0EQNJAw/kEQBA0kon2C2pnx8PM8eNeTvPziKyy2xAjW2HBFVhq9dN1qBcGQ0mvXeRj/HqHXLoyqmPHw89x72xO8+ooBePnFV7j3ticAGtH/oBn04nUexr8H6MULoyoevOvJ2f1u8eor5sG7nhz2fW/yDb9p/e7F6zx8/j3AQBfGcOflF1+Zr/bhQuuG3+pn64Y/4+Hna9asLE3tdy9e57UZf0m7S7pL0j2SPl2XHr1AL14YVbHYEiPmq3240NQbflP73YvXeS3GX9II4AxgD2AT4EBJm9ShSy/QixdGVayx4YosMkJztC0yQqyx4Yo1aVQNTb3hN7XfvXid1zXy3wa4x/bfbL8M/AzYtyZdaqcXL4yqWGn00qz9xpVn3+gWW2IEa79x5WHvA27qDb+p/e7F67yuCd/RwINt2w8B23YeJOkw4DCANddc8zUL+9asd73mzw6WoxdE9qowY9bzPHTnk7w88xUWGzmCMRutyEqrLg2zCsteQIZE9qr51c4w7/eMdZ/nvrZJfkg3/LHrrszVswY2CNHv1yZ7QejF63xBqcv4q0ub52qwJwATAMaNGzfX/uHESqOXHvaj3aCP1rme64Y/zK+Bpva7F6nL+D8ErNG2PQZ4pCZdgqAWmnrDb2q/e426fP43AutLWlvSYsABwIU16RIEQdA4ahn5254l6RPARGAEcKbtaXXoEqTY63gMD4JmUdsKX9uXAJfUJb8bTTSCMx6ecwLu5ZmvcF9DVhcHzaOJv/H+iPQOmaYawYfu7L7o5qE7h396BWiuMWhiv5v6G++PSO+QGcgIDmdentnPopt+2ocTLWPQ6mvLGAz3VANN7XdTf+P9EcY/01QjuNjIfhbd9NM+nGiqMWhqv5v6G++PMP6ZphrBMRt1X108ZqPhv7q4qcagqf1u6m+8P8L4Z5pqBFcavTRj37jy7B/AYiNHMLYB6RWgucagqf1u6m+8P2LCN9PklYdNXXQzZqMVu6YaGO7GoKn9bvJvvBth/NtoqhFsKk01Bk3tN8RvvJ0w/kGjaaoxaGq/gz7C5x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDSSMfxAEQQMJ4x8EQdBAwvgHQRA0kDD+QRAEDWTRuhWogqN326BuFYIKaer5jn4H80MjjH9TiR9FEJRlYf6NhfEvzMJ8cSwI0e8g6G3C+AdBsMDETW/hIyZ8gyAIGkgY/yAIggYSxj8IgqCByHbdOgwKSdOB+ysStzLwREWyQnZvyA/ZIXu4yl7L9qjOxoXG+FeJpCm2x4Xs5sgP2SG7CbLbCbdPEARBAwnjHwRB0EDC+HdnQshunPyQHbKbIHs24fMPgiBoIDHyD4IgaCBh/IMgCBpIGP8gCIIGEsa/RiTtkP8uXrcuVdLUfgNIujL//UrdulSJpCPz3x3q1qVKJL0n/127bl06CeOfkTRS0qck/UrSLyUdLWlkYbHfzn+vLyynX5ra7xaSVpO0j6S9Ja1WgcjVJb0V2EfSlpLe1P6qQD6S1pF0kaQnJD0u6QJJ6xQW+4H89/TCcvpF0kqSTpd0k6Spkk6TtFJhscfnv78sLGe+iWifjKRzgeeAn+SmA4EVbL+noMxJwB3AnsDPOvfbPqKU7DYdGtnvrMeHgc8DVwEC3gqcZPvMgjLfDXwI2BG4McttYdu7lJLdpsMk4AzgnNx0APBJ29sWlHkOsD0wCvhr+y5SvzcrJbtNh8uB39F3rR8M7Gz77QVlXgGMALYAft+53/Y+pWTPizD+GUl/sr35vNqGWObKwNuBr5CM0BzYPquU7DYdGtnvrMddwJttz8jbKwF/tL1hQZk72P6DpM/bPqmUnHnoMLnT0EuaZHu7wnJXAyYCcxk828XzdkmaanurjraiqRYkLQa8Cfgx8OHO/bavLSV7XkQxlz5ulrSd7UkAkrYF/lBY5n/aPk7SmlUZvC40td8AD5Geelo8BzxYWOa3ga2A/YBajD9wtaRPk566DPwb8BtJKwLYfnKoBUq60vaukiZWYej74WpJBwDn5u13A78pLPP7tt8n6X/rNPTdiJF/RtIdwIbAA7lpTZJr4lUKPZZKuo00KphsuxJ/bxcdGtnvrMePgDcCF5CM4L7ADcBfAGx/s4DMlsvrX4Gfd+6vyNV37wC7bXvI/f+Sbgc+BnwPOIg53V3YvmmoZXbR4TlgKeCVLH8R4B99KnjZAjJvB/YALgR2Zu5+D/mNdrDEyL+P3WuQeSkptetSkp4l+z9bf0tcjF1oar8h+Z7b/c8X5L/LFJS5F8nltQswtaCcfrFdR+TJ54FPA2OAzpuqSf+PotgueV7743uk630d0vmeY44nt9dC40f+kpa1/WzrkbeTKu7Mki6wvW9pOR0yG9nvXkHS5rb/VLHMXWxfJemd3fbb/lUFOnzO9hdLy+mQuZHtO/uLpqroqeO7tj9WWs78EMZfutj2XvlRuDX6bFHkEbgfPVYFts6bk21PLyyvkf3OMk+1fZSki0h9n4MqIjAkjSGFPe6QdbgOONL2QwVlnmj7BEk/6LLbtj9YSnaHHvsAb8mb19i+uLC8CbYPk3R1l92VRFhlPTYHdsqbv7N9axVy+9Wn6ca/F8gLQb4OXEMywjuRJkV/Uadepamr35K2sj01x9vPRRUTczns8KekKBCA9wIH296ttOw6kfRlYBvg7Nx0IDDF9vH9f2rhR9IRwGFA6+nqHcAE2/Wtewjj34ekzYCxtM2FVPQo/CdgN9uP5+1RwBUlwy075Dey33XST4jtLba3qED28sD7mfucVzHZfCuwhe1X8/YI4OaK4vxHkNaWjGXOfg/5xH4X2bcC29v+R95eCri+in73R0z4ZiSdCWwGTCNFukB6HC9uBIFFWgYwM4OKVl83td8AkvYCvgisRfotVDnhPF3Se+lbaHUgqf9VcAkwCbiNvnNeJcsDrTml5SqUexEwk3r6LVKUUYtWxFFthPHvYzvbm9Qk+1JJE+kzBP9G+oFWQVP7DXAq8E7gNlf/CPxB4DvAt0g32z/mtioYaftTFcnq5MuktSVXk4zfW+hLgVCaMTWOtH8ATJZ0ft7eD/h+TboA4faZjaTvA9+wfXtN8t9JWvIv0mTQ+fP4yFDJbWS/s+yrgV1bLogqkTTS9syq5WbZRwPPAxcDL7Xaq4o5l7Q6aZJfpEn+xyqS+xXgStuXVSGvi/w3Mee1fnMdeszWJ4x/QtJbSI+Fj5F+EFXmHPkg8Hvbd5eW1UV2I/ud5W9Ncvtcy5xGsAof8D3A30n5Xn4H/MH2M6XlZtmHAycDT9MX7VRJhJekH5P6+3vbd5aW1yH7HaS8PosA/6RCN5+kk0jn+o8tv3/dhPHP5B/jp+jwB1axFD1fGDuSfM9TSRfJ723fUoHsRvY7y7+MNALu7PuJFclfkxThtANpxe/TFU34/hXY1vYTpWV1kb0L6ZzvRFrgdAtpFHxaBbL/RnK3VO7mywOdHUnJ7Z4j3/RtXzDgB0vqFMY/IemqquJ9B9BhCeAjwH8Ao22PqEBmI/ud5RZN6jUP2WNIBvCtwOakCdDrbH+5AtkXAgfYfqG0rH7kjyC5fd4G/Dvwou2NKpA7EdijDjdfmw6rAfuTrvUValp1DMSEbzt3SvopyQXS7gKoIuTxs6TR39LAzaQLY670r4Voar8BrpD0LzX5gB8gpXQ+xfa/Vyz7FeCWPOfRfs6rCPW8kpRf53rSud66I+KrJI8C10j6LdW7+f4P2IQ+V9+7geIriwcijH8fS5AuiH9pa6sq5PGdwCxShsFrgUkVTgY2td8AhwPHSnqJin3AwJYkN8BBOcPm3cC1tquIAPl1ftXBraSspm8AngGelnS97RcrkH1vfi2WX1WyEimv/9Okp7wnbM+qWIc5CLdPjyBpGZIx2JH0WPh32zvWq1V5mtpvAElL0+f/fi/pxjO2VqUqIvf9A6SnvdVsN6Kkp6SNgfHA0cAI22Pq0iVG/hml0oUfAjYFZpcxrCLfiaQ30Of/HUfKKV+J+6Op/W7TYQVgfebs++8qkDsFWJwU338d8JYqJtmz7PVJ8fabMGe/q4j2+QTpnG8F3A+cSXXX+ijgWOa+1quonrYXqd9vAVYgVY+r9FrvJIx/Hz8G7iTdlU8ilXi7oyLZXyGFv30buNH2PyuSC83tN0plHI8kpRm+BdiO5IuuYgJ8jyqS2PXDD4ATSAvM3kYagVe12nQJUkrnqTW4Pc4m1VDYizTRfAhQ1TnYg3Stn2b7kYpkDki4fTKSbra9paRbbW8m6XXAxLojYUrT1H4DraIyW5PmGraQtBFwou1/KyhzwJW1FU0+TrW9laTbbL8xt/3e9k7z+uwCyOyaOrxFFQvM2vp9a2sdi6RrbXdN8DfciZF/H61R59PZHfEYKQFUMbLx6Xb3rWyhFc3tN8BM2zMlIWlxp5zvxer3ZmoL7WtjpqRFgLuzG+ZhYJXCMqfSd847nzKqKmrSutYflbQn8Ajpqa8YStXDOq/1Vv+rCi7oShj/PiZk/+9nSSXXlgY+V1jmXoW/fzA0td8AD+UMl78GLpf0FMkgFKOqBWTz4ChgSeAI0grnt5FcIMVwPdXDOvmSpOWAY0i1FJYlTbwWo844/nkRbp9BIukQV1BsXNJeLlzcYn5oSr+VcvsvB1xq++XctoLtpyqQfZNrrGXciaTTbX+yAjlfsP2F0nIGi6TjK1pkd5jtCaXlzIvK0ucOA46sSM5JFckZLI3ot+1rbV/YMvyZKysSX2tq3y7sUJGc4hXT5pP3VCSn6kV9XQnjP3iq+oH2miFoar+hOp1+U5GcXqPXznmjrvUw/oOnKv/YRyuSM1ia2m8o3HdJqynVs52cc740jZ5xdWWqutZ7Ys4rjP/gKXa3ljRS0qck/Qo4TtLRefFVL9DUfhclrzG4gZTi4t3ApJz5sRcoec7XkXSRpCeAv0u6QFIVkT6DoWS/V5J0uqSbgAslnSZppVLyBkMY/8Hzh4Lf/SPSqsPTSdWdNqavsHfdNLXfUPbx/D+BLW0favsQ0orX4wrKGxBJa7Vtlkyv/FPgXGA14PXAefRVcqub8wp+98+Ax4F3kW7200kLzmojon2AvLhnNKmq0PNt7bvbvrQC+d2Kec/VVhWSPmD7BxXI6al+Z/kP2F4zv1+x1OKjnN1yj7bIosWAS2y/vYS8Nrnbk67139l+XNJmwKeBnWyvUVJ2lj/Z9rYdbZNsb1dY7nhSLv/RJPfOI8AFVfy+s/yptrfqaKstpTjEyB9JRwAXAJ8E/ixp37bdp1Skxs2SZl/8kral7Ih7XlQVi95r/Ya20X7hVacPk3z9X5B0Aqmg+j3ZDVakvq6kr5Fy6bwL+E2WezkwmZTfqAqulvRpSWMlrSXp2KzLivNaBfxakXQqKWrtWuCrwNfy+yMkFS8ik7la0gGSFsmv/al5or/xI/+82nR7289LGgv8Avix7dNaqQ8q0OEOYENSjneANUn5dV6l0IpXSbf2twvYoIosi3X0exA6zR75F5ZzwkD7SywGk3Q78Ka8qnkF0uh3M1dYRlPSvQPsdonkcpL+YnuDLu0C/mK7+I0vr/RdilRLQaSBd6ucYy0rfWOFb0qr+jyA7fsk7Qz8IvtAqwrJ2r0iOe2sSkrm1rmISaRMk1VQR78Hyq8j0grn4tS00vfFVr0E209JuqtKw5/l1rHSd6akbWzf0NG+NVBJ/YheXOkbxh8ek7SFc93Y/ASwF+nx+I0lBUta1vazpJqec1HY7XAxsLS71MuVdE1BuXX3GwbOr1PUDSDpVNtHSbqILqGFtksufFpXqYRji7Ht2yVlS9rF9lWS3tltv8tWjjsU+K5S7YiHctsawLN5XzEkbZRzRnUNa7VdWzWvcPukWqqzbD/WZd8Otov5oCVdbHuv/Chs5nzSKPII3Av0cr8lHWX71ILfv5XtqTmdxFzYvrag7AGzVxaWfaLtEyR1CySwq6kfsRppwlfAQ91+8wVkTrB9mFLJzE5cZ/bcMP49kGq2V+iVnCN1UpXPv9eQ9POSqax7DfVYXqE6COMvvUp6FGwVlqhlFJpD7sbS5oor/CjcTYfKE4z1Qr879HmwopDHvUgZNdci9b3K+sHd9Klqont54P3Mfc6LF4/v0KPSa13SCGBP5u538foN/RE+/7TAaGdSiOE5wHWu+I4o6UxgM2AaKdIFqiuiPocqlQrrnX63U9W5P5W0uve2qq+3mrmEFNZ6G33nvA6qzq9zEWlyue5+z6bxI3+YHfK1M3AgsA1wGfBd2wOFpQ2l/Nttb1KFrHnoMcb2Q/M+csjk1dJvdS+wAckgLGG7+KAo+4B3tV2ZIehv0pHU74ttr16BDj2RvlrSIhX/72dXD+sVYuRPetYmLcK4GTiA9Dh+N/C/FalwvaRNbN9ekbzWgqo7bD8raQnSKs835VjwU2w/U4Ealfcbeibs7ljgEknXAi+1Ggu7Ab4xwL47C8pt58eSPkKKNmvvd9G5tbzCdwxwpe37WoZf0gdtn1lSdua3kv7F9mUVyBoUjR/5S1oK2Bf4N2AUyeXwc9sPVqjDW0iPhY+RfhDFyxlKmgZsbnuWpAnAC6QFbrvm9q4heUOsQ+X9HkCXSie7JV0GPE+HG6Cm+P/KkHQ4cDLwNH1PX0Xn1iSdAuwI3ATsDZxq+/S8r5InEUnvAH5CWtz1T2qe44Ew/kj6B2mUfw5wDx3ugComHyXdA3yKuQ3B/QVl3mF74/x+jh+ApFtsb1FKdpucyvs9gC5VTwDWmtelTY8Jtg+rUN5fgW1tP1GhzNtISfRm5QnnnwJ32T66wlX8fyPlFuqZOZ5w+6RMfgY2yq92qpp8fMD2hfM+bEj5s/oSuP1J0jjbUyRtQF+h69LU0e/+qHoC8IoecQNUfQOaRnrKrJJFbc8CsP20pL1JtavPAxarSIe7gT/3iuGHGPn3BJL+B1ie5AJp94MWu/EoFbI+DdgJeIJUWOPB/DrC9p9KyW7TofJ+D6BL1ZPdrVwvL1GjG0DSpbYrS7Mh6XxSGu+rmfOcFwv1lHQx8LXORWySvgR8xnbxBJeSfgisA/yW6uZ4Btap6cY/jwJubbkaJH2elPXwfuDIKiJ+al71uAzpolyUtOrx76Vltsmurd/9oYrSWfcSknYADrJ9eAWyDunWbvusgjKXyDJe7LJvtO2HS8luk9M1kV+dczxh/FN2y+1sv5AX3nyTFPK5JfAe2+NrVbAQsbK5O1Wu8FXKrLk+MLt6me3fVSR7C+AgYH/gXuBXrUnQ4cYAIa5Avfl16iR8/mmk2fJBvhP4vu2pwFRJH69CAaXShR8iPQ63G4KSI+AnGGBlM+lpoCg19Xte6axXLSm7TYcPk3LMjwFuAbYDrgeK5XrJ8zkHkAY3M0iVpGT7baVkdtFhfeDLwCbMec5LXm9TSHMN01tqtO0zBf/nLSSNIoX3dl7rteX2aXwxF9Iar6UlLUIKc7yybV9V9WR/TCprN55UZGIM/WS8HEJOJ6VzvhQ4BFjH9tr5VVVitTr6DcnAv58U9tf5mlGBfEiGf2vg/mx8t6TPOJXiTtI1vrftHfNI/5XCMjv5AfBd0qDjbaRSnqVLdx4DPAO8mOXvbftt+VWV8T2b9P9fm1Qs6T7gxopkd8d2o1/AB0khnjcBl7a1b0laEFKFDjfnv7fmv68DrqpArkg/wAmk0edXgbUr/N/X1e/vAzv2s++nFfX9xvz3FmDx1vvCMt9BGu0/SFrAuCtwb1XnO+swNf+9ra3t9xXJXhv4DKly2bnAFjX0+9a2tmur/N93vhrv9rF9pqSJwCpAe4TLY8AHKlKjFVr5tKQ3ZNljSwt1ugLrXNlcV78/NMC+g0rLzzyUY85/DVwu6SlSZa1i2D4fOD8vbNwPOBpYVdJ3gfNdTdjpzPyUfbekT5DKWa5SgVxs3yvpAmAJ4H3ABqSbbxW0rvVHJe1JOtdjKpLdnTrvPL3wAt7b9n6Hjn2fqEiHDwMrAG8B/gY8Dny0sMylSBN+F5Aqd/0HsEbF//vK+z2ALivWITfLfiuwD7BYW9sKVfUbOIwKnriyvK1J1dLGkFwwvyQFXJSUuQ59I/5fAO8GRlZ8jvcClgPeQApznQrsU6UOna+I9mlb2dllpWuvJKE6xEMcCtcLK5vnRYl+5+/dAfg/0qriDwJfAtYluZ32t339UMucX6q49iStB2xOyvFUaX6l/pB0uu1PDvF3vgrcShroPMvc13ptsfYtJB1v+8tVyowJ3zln/jtXeVa96rM/jizwnecBN5NWNe/FnJOeexWQ91oo0W+Ab5FCHD8M/AY40WmSe1/g64Vkzi9Dfu1JulrSyvn9+0jplfcAfi5pSA3uArBDge88CTifdLNfmlTGs/3VC7ynaoGN9/kz5yig8zGoVx6LhtwQ2D50qL+zAKVuvq+zfRuApOm2r4MU791aENQDlLj2Rrkvp84RwPa2Z0hakpRjf1jG+XvhqNhV+UAzRv6wkaRbc/Kn1vvW9oZ1K5cZckMg6VOS5pr4lPRJSUcNtbzXSKmbb/t1f3zHvqpyvdTBPyWNzu+fB/6R378EjKhHpfJI+kheX4ASZ0p6Jv/Oiyd1GySVDzRj5J98nquSwt/aWYvC0RfzQYlRwQdJ+Xw6mUCKPz61gMz5pdRo6HOSlrT9gu1fzxYmrUuKO+8FSvT9aOAySb8kLXq6StKlpPxOvZLSokS/jwR+mN8fSPrNr0MK5/42qf91EyP/GvgW8Kzt+9tfpMyD36pZtxZ/KPCdtv1yl8ZWXv1eoES/IVVqW6pL+3MkY1ALkh5o29x1qL/f9jXAm4FHSaGHU0mj/k/arm2uQ9JabZunFRAxy3Yr1HIv4Ee2Z9i+gu7XQR2cV7XAiPaR/mz7Df3su832GwvL3wgYDUy2/Xxb++62Ly0o9zbg7e5I5CZpVeCK0v0eQK/iidWUitdc2hnRJOlg0uKvj5WUP4BelRSPrwtJ25Ou9d/ZflzSZqQKcjuV7Lekm0jF058iJWzcxfa0vG92XYuC8seT1lWMJrl3HgEuKPn7HpReYfx1j+315nffEMk+AjgcuAPYgpRF9IK8r2ion6T3kyb9jiGtbgbYirTK94wSIZaD1Kt4YjUNUDtY0jTbm5aU3x+l+y5pwNoJtvcpKPtrpFH3LcB6pDKOHwdOAf6f7ZkFZe8F/D/SvMZFtj+S298KHGt7z4KyTyUtJvsRKZcWpDUO7wfutl0qom2ehM8fbpT0EdtzrGrNk6FTC8v+CLCV7ecljQV+IWms7dMo7Hqx/SNJ00lhcG8gjUimASfY/m1J2T2QWG2g/21RV6ikT/W3ixSGWJLtSXNb55AWPFXp3tuTVE1rZs5m+giwme27Swu2fXF2LS1j+6m2XVNI5VtL8q+2N+hslPRz4C+UC2eeJ2H84SjSkveD6TP240hRH+8oLHtEy9Vj+z5JO5NuAGtRwQ8zG/mihr4fViUlc3uqo12k1caleVzSNrZvmEO4tDXlk6sNFFdewt/dzmrAbqRJz4NIaxzOablACvNia3Rv+ylJd1Vh+AEkHWv7q8BTkt5j+7ysxz+U6vt+pqD4md2uNdJK52JPO4Oh8W6fFpLeRhoBA0yzfVUFMq8CPmX7lra2RYEzgYNtVxp+V9WKZknfB37Qiq/v2PdTF86vI2kbUmKvHzLnDf/9wAG2J5eUP4BeR9k+tSJZi5NuAl8DTnLhXP6SngbaaxW8pX27sMuptlX8SrUEvku66bfcPmuQVhp/3Cl9fC2E8a8RSWNIkQiPddm3g+1S0S796XOzKyhm3Qvkie2P03bDB75j+/EadapivmNxkgvmQFISvQuBM124mlX2r/eLO0osDrHs2dd15zVe1TUvaTXShK9IFfPm+s1XTbh96uUF6Leq1h0V6wLJDVALkg6zPaEqeTnKqWtpvRop6uqTdBbpZvdbUkqLP5eU185Axj37v4sZf3pgFX829o8BSPoC8IUq5A5EjPxrJCec6realqsrqlI7VSbRk7QvMMb2GXl7MjAq7z6u5ROumgqifV6lb1Vv+w+/luLxs4WX7/crpH6LlM65VblPpOyeryslux99eiJhZIz86+V0YGfSYqZzgOtc0d1Y0nN0MQDUZwiqjDw5llS/oMXipAm4pUgrXYsZ/y7/99m7SIapGLYbuaiz6rmzQdATiyjD+NeI7SMliXQDeB9wuqTLgO/avrew+CtJ0R+/An5m+4F5HF+avSuUtZjt9nQe19meAcxQKnRSDNu1ZZHM0Uwrd4byStobeKTk5KP6L6IuUirtSqnazdjBVjXJnYNw+/QISlWdWtW0PtO57qCQzOVIResPINUr/jnpRvBkYbnbknLIP6uURfPTpDxDtwOn2H6msPyBFvb91fa6JeV3kVmJIZJ0DXCo7fs62tcDJrhgPVtJVw+03xUWkYfK3YzjSQu7rmz/30v6oO0zq9Chq15h/OsjjzL3JS00GUUahf+8Y1RahR6LZB1OJxnfosUtJE0DNrc9K6daeIFUYWnX3P7OwvLPBq7psrDvo8DOtg8sKb+LPlWF2PabrkTSn2xvXlqHXqHCKJ9TgB1Jq+j3Bk5thdXW7fsPt0+9PM7c1bS2zo/nxatpSXozKeRvJ+A64B22f19SZmYR261J7nFtP4DrJN1SgfyjgV9LOog5U1ssTsrBUjVV+YAHmlOoPMGZpAm2D6tabqYqN+PepJXNs3KUz08lrWP7aGr2/Yfxr5fzSAZ/o/xqx6QngSJIug94GvgZqYbrrNz+JkiFTUrJBv7clsDtT5LG2Z4iaQP6Cl0XI8fyv1nSLkArj89vqljY1w9VGaIrJJ0MfLY9sEDSiUAdfR9XhZD+Umqk6bbiZRwXbQ10bD+d51cmSDqPmmtHhNunoWT/b+vkt6J8Wriw/3c5UiqDnYAnSP7+B/PrCNt/KiV7AJ1Wq3vhjQpnNM1uxv8DtiElWIOUUPBG4MNuyypbBZIutb17BXIGXM9h+8SCsi8Gvta5zkHSl0hze7VFYIXxr5E8CrjVqX4Akj4PvIuUdvbICiJ+akXSMqSiGouSVj3+fR4fKalL7bHXVazwzXLWoe+JZ5rtv5WWOYAuOwAH2T68JvlL2f7HvI98zd+/BIDtF7vsG116ZfVAhNunXk4GtgNaaWffS/LBbwl8j5T8rFIk7UZKc7tbQRntK5of7GwvHW3UD5X4X1V/RlOysf9b1ucLVLzaVNIWpMRy+wP3UtC92SZzNLA6abD1sqRVSEkdDwVeX1D0xll+f/vD+DcU226tNnwn8P0caz1V0sdLCs7+7u+RLvxfk/Kq/4hkhE4uKZvk6ul3ZTPpaaBqiofWZurOaNrJPlRg/PN8zgGkwc0MUlixqgjxVKpJ/V+koIrFJZ0GfJN0vZeOuZ9CyhvVyhbbea0Xc6/OizD+9SJJS5NCHXcF/qdt38jCsr9Bmui9HtgDmAR8zqmWQGlqW9k8AD+rSM7FwNJuy+TaIs/DVE1VESd3Ar8H9rZ9D4CkoyuSfRiwoe0nJa1Jugm8xfakCmQfQ3Llvki6xs6vem6lPxq53LuHOJU08TaFtOhpCoCkLUl1Vkti29fYfsmpiPn0igw/TtWLtiBFO70PuFnSVyWtXYV8STtIukPSNEnbSrocmCLpQaVSg8Ww/SF3SWWd9xVNZd0PWwFIKl3D912kxGZXS/pfSbtS3Y1nZsuVmFey/6Uiw4/tb9neEfgEKZXzlZLOza6vWokJ35rJvshVgD/ZfjW3rQ68rmTKBUl/A/6jrenr7dul1xi06bE81a9svgH4EKly1kXAfravy2Gup9veobQOHfqsWNM8R7sOVU02L0VaS3EgyeVxFmk0fFlBmY8z55PdAe3bto8oJbtDj02z7PeR5tXOrUJuv/qE8a8PSe+1/ZP8fo78/ZI+Yfs7BWUPFFJo2x8sKLvWlc2aM7/7HAW8S0f95OiW/wNeBT4IfAlYl5TfZn/b15eSPQ+9Ki8enyf4300qoFMytPiQgfa7YL3qHFl1AOl6f5B007nYBWsWD5Yw/jWiGisM1YmkfzD3yubZVLCyeXYqA0n7ZbdXa9+fbb+h3w8vuOzanjrUvW4EJPfLn2yPKSW7iy7rAZuT3J23VyW3apTSaN8KXECq3tV5rRdNpTIQMeFbL+rnfbftoRUsvX+A3bb944Lia1vZnPmcpCVtv9Bh+NclRYCU5HW2b8vyprf8/7ZvasWEF2Qqcy/oa/FyScFKid3eY/sJSe8DPkcq4/iFnOahWBlJSRcOtN8FS0gCJ9Fn8JcuKGe+CeNfL3VWGNq6S5tIqQZGA8WMv+1DS333IOV3NQa2/wp8tbD49iCL4zv2FV3ub7uSCfV+GGX7ifz+CGB72zMkLUmKNCtZQ3h7ksvlHGAyFebUsf2FqmTNL2H862WjvOhHwLptC4BE4Vh3259svVdagXIwcBzph1g0zj/nWnnG9vc72j8JjHDhIuY1jwRre+qQtBbwtHPKbElvI02+3gecYbvk6P+fbStan6evothLQOliK6sBu5EmmQ8ilSs9x/a0wnKR9BFSBtm78+/s+/St4j/E9s2ldehXt/D514ek9UmLfjonOtciFde4p7D8RUkrHI8hjYi+bPuukjKz3D8Db+o0NkrFxW+0vVlh+dMZYCTYmYdliGWPBJaxPb2jfRXg2ZITgUrlKt9h+5EcangF8GVgM+Cftj9cUPbOwBnAL4EVSfmcLiXld5pou3SoaUuPxUk3ga8BJ5V0N2V5fyZl9fynUhbZY4B/Ia3iP8H2TiXlD0TE+dfLt0g/+PvbX6RFX98qKVjS4aTiKVsBu9s+tArDn3G3Uabtl6jmkXw14DOkYuankUaFT9i+tqThz3ybZPA62Y3C5xxYwvYj+f17gTNtfwP4ACnZWzFsXwO8mbR+5Z+k+YeXgE9WYfglLS7pncBPgMNJ56GKcOZZtluZavcCfmR7hu0rqCGNdjsx8q+RgSJLNEDhjSGS/SqpnsB0uswvlBx9S7oNeLs7ErlJWhW4omS/u+hS9Ujwdtub9LNvmu1Nu+0bItmzrylJNwHH256Yt28t/cRVF5LOIt3of0uqVPfnCmXfBOxJSudxP7BLy93UGWZcNeHzr5eBUjiUjvzYkAFcToVlfw34jaRjmLOYyldJi82Kk43+niTDP5bqRoIDPdmUfhK/StK5pNH3CuQc/nlRYelonzrnWd5HmmPYADiiLcmakmgvW1D250kr+EcAF7YZ/reSk+vVRRj/erlR0kc6V7VK+hDpsbgk3yKtqL2/Q/aovK9YgRHbP8p+95NIIzKTkl+d4I7i4iXoGAmeWOVIEHhc0ja2b+jQaWv6kn+V4ijSwrrVgR3b3BGrkRKflaTOiJva3Nu2L84T7cvYbk/mN4V0Lmoj3D41kt0c55NGXS1jP44U8vcOFywuUqfLqW6yy6sVbdL+Ayg+EpS0DXAu8EPmPOfvJ610nVxKdpa/H7AecFvL5VMFkkbQF3GzGRVG3NSJpGNtfzW/f4/t89r2nWL7M7XpFsa/fnLIXcsQT3MF5QQl3WN7vfndV1CfYbuiuZN80/84becc+I5TecmScr8LbEJKHb0rcJHtL5aU2Y8elc6z1Ekvr+IPt08PYPtq4OqKxdbpcupGZW4ASbu0brCS1nZbxTRJ73Th9BJ5onvA0oKF2AnY3PYreXHV70kJ9SqhxnmWOqltFf+8COPfXI4Czpd0MF1cTjXo85sKZX2dFGcOKe68ffT1WQoaJEn7AmNsn5G3J5OS2wEc1+4WKMDLtl8BsP2C1H95qaGm5nmWOqlzFf+AhNun4dThcqqbjqyes9932y4g+w8k3/6DefsWkgtmKeAHtnctKPsFUiI9yKvK83ZrrqNkeG9t8yx1IukVUr9FiuBrVe4TMNL26+rSLUb+DacOl5Ok5+hiAKjOENQ5GlvMc6auvs72DGCGUqrrktQWU15nxE2d2C6duuI1E8Y/qIMrSeGFvyItuilWtKYf1slx52p7T94unfxshfYN259o2xxFQTrDeqskh7Ku3BnKK2lvUiqTOuaZakHSYbYn1K5HuH2COpC0HKlo/QGkxW4/J90Iile0ygts+qVwbp+zSYm+OifaPwrsbPvAgrI7n7hm76J8iOs1wKG27+toXw+Y4ILFXHqNuqN8ZusRxj+oE0mLkBa7nA6c4hqLW1RBTuD2a1Jem/bVzYuTCrv8vZ+PLtQMtHZEbcV1mkDpeaXBEm6foBYkvZkU8rcTcB1pUdvvK5I9UMTNsbZ/UUp2juV/s6RdgFYen980YKJ9oHQltSY4q4Fiq+fnhxj5B5Uj6T7gaVI906uAWe37bd8096eGVH5tETf96LNaydXcvYCk7wEzgM+6zehIOhFY3fZhtSlXkFy7ol/qfNKNkX9QB/eRfM/jSbnN2+PNDZT2/9YZcdONS5hzrcFw5BhS4fp78s0WYAvgRqBYHYEeYJm6FeiPGPkHjWMeqS3+anvdivXpCR9wFUhahz531zTbtWa2rBNJS9n+x7yPLEMjY2+D3kTSbpIur0DUZKXyep3yPwrc0OX40vzvvA8ZHtj+m+2LbF9ESmY37JE0WtI4SYvl7VUknQLcXateMfIPqiZPdn4PeD0p8uUUUv1aASeXzq3TaxE3klasIsS11+iVkMeSSDqKlC77HtL1dRrwTdL1/lXbj9amWxj/oGok3QwcDVwP7EH6IXzO9mkV69EecVNVNtUdSL7vV4EPAl8ipVl4HbC/7etL69ArNMHdJel2Uu2EJyWtSboJvMX2pJpVC+MfVE+X1LaV+tklrTjQ/pKjcEk3AB8ClgYuIj1pXCfpTcDptncoJbvXkLSI7Vclfd32f9StTwm6XOv91tGomoj2CepgeaVi2i3Uvl3a7UPKYtrKJbQmqb6qgOWBByib4uF1tm8DkDTd9nWQwlsllS7d2VPYfjW/3R8YlsYfGCPp223bq7Rv2z6iBp2AMP5BPVzLnAtd2rdN4RzvtteG2bHnF9q+JG/vAby9pGzmDLI4vmPfYoVl9yq15rUvzH92bPdMDqNw+wSNRdJU21t1tE2xPa6gzH2AK2y/0NG+LvAu55J/w40BXG0C/mR7TJX6BGH8gxqQNFCIn23/uCI9JpKqWf2E9MTxXtJk3Pgq5DcJSffS52rrxLbXqVilSmjLGNsV2/tUpUsnYfyDypHUrWarSK6f0bYrcUfm0egJwFty0+9IVaZKTvj2rDEIhh5J04EHgXOAyXTc/EpmkJ0XYfyDWsmlBA8GjgNuJ8X531qvVuXoZWNQEklrAU/bfiZvvw3Yj5Tq4wzbL9enXTkkjQB2IyUx3IxUrvQc29NqVYww/kFNSFoUOJSU82Uy8GXbd1WswyjgWFKs/8hWe8nc8r1sDEqSM6e+w/YjkrYArgC+TPof/NP2cM7vA8wuYH8g8DXgJNvdnoArI6J9gsqRdDhwJKmi1+41Vpg6m1REZi/g34FDgOklBToVUL8UuLTNGFwjqXZjUJglbD+S378XONP2N3I9h1vqU6s8+TzvSTrXY4FvUziibTDEyD+onFzM+3GSoZ3rAnTBQuIdeky1vZWkW1syJV1re8BKX0Mgt9MYXEgyhg+XlFsn7cVcJN0EHG97Yt6+tapzXjWSzgLeAPyWVKnuzzWrNJsY+Qd1sCGwKsn33c5awCNzH16Mf+a/j0raM8suGnLYYQxO7CVjUJirJJ0LPEqqY3wVgKTVgWHp78+8D/gHsAFwRJriAioonTkvYuQfVI6ki4HPdE7sShoHnGC7kkpHkvYihXquQSojuSzJIA8YkbOAMl8lGQOY86mndmNQkjyx/2/A6sC5raccSVsCq7SeAoLqCOMfVM5A+U0GqvUaLNxI2g9YD7gtjH39hNsnqIORA+wrnt8mrzPod9RTMt+KpF1a2UMlrW373rZ976wgr1EtSPousAnwR+CLkrax/cWa1Wo0UcwlqIMb+ymm8iGqyX0yJcuZCuzT9r71KsnX297/smPfZwvLrpOdgF1sHw/sTIrxD2okRv5BHRwFnC/pYPqM7ThSYrN3lBZu+6zWe0lHtW9XgPp53217OPFyDnPF9gtqm/kM6iGMf1A5uVLWm/Mqz5bv/zdVFFPppk6N8jplD+cJuI0ktSb4Baybt1sT3cMy1LOXCeMf1Ibtq4Gr69ajYtbJ+X3U9p68XbKOQN1sXLcCwZxEtE/QOCQ9R98oe0mglV65eLilpAEXkA3X3D5B7xHGPwiC4nTccOfYxTBe39DLhPEPggqRtC8wxvYZeXsyMCrvPtb2L2pTLmgUEeoZBNVyLCmXT4vFga1J4Y8fq0OhoJnEhG8QVMtitttzGl1newYwQ9JSdSkVNI8Y+QdBtazQvmH7E22bowiCigjjHwTVMrmf1c0fBW6oQZ+gocSEbxBUiKRVgF8DLwE35eatSL7//fICuCAoThj/IKgBSbuQykcCTKtpdXPQYML4B0GFSFpxoP22n6xKl6DZhPEPggqRdC9psZOANYGn8vvlgQdsD+cUD0EPERO+QVAhtte2vQ4wEdjb9sq2VyIVkR+WufyD3iRG/kFQA63i8R1tU2yPq0unoFnEIq8gqIcnJH0W+AnJDfReYEa9KgVNItw+QVAPB5IWdZ1PCv1cJbcFQSWE2ycIgqCBhNsnCGpA0ihSkrdNaStob3uX2pQKGkW4fYKgHs4G7iRV7zoRuA+4sU6FgmYRbp8gqIFWtI+kW1v1ayVda3vASl9BMFSE2ycI6uGf+e+jkvYEHgHG1KhP0DDC+AdBPXxJ0nLAMcDpwLLA0fWqFDSJcPsEQRA0kBj5B0GFSDqd7oXMAbB9RIXqBA0mjH8QVMuUtvcnAifUpUjQbMLtEwQ1Ielm21vWrUfQTCLOPwjqI0ZeQW2E8Q+CIGgg4fYJggqR9Bx9I/4lgRdauwDbXrYWxYLGEcY/CIKggYTbJwiCoIGE8Q+CIGggYfyDIAgaSBj/IAiCBhLGPwiCoIH8f23WirmKEtvxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "fig, ax = plt.subplots(nrows=2, sharex=True)\n",
    "ax[0].bar(np.arange(1, len(toa_sw)+1), toa_sw.values(), alpha=0.5)\n",
    "ax[0].bar(np.arange(1, len(toa_lw)+1), -np.array(list(toa_lw.values())), alpha=0.5)\n",
    "ax[0].scatter(np.arange(1, len(toa_lw)+1), np.array(list(toa_sw.values()))-np.array(list(toa_lw.values())), marker='x', color='black')\n",
    "\n",
    "ax[1].bar(np.arange(1, len(sfc_sw)+1), sfc_sw.values(), alpha=0.5)\n",
    "ax[1].bar(np.arange(1, len(sfc_lw)+1), sfc_lw.values(), alpha=0.5)\n",
    "ax[1].bar(np.arange(1, len(sfc_lw)+1), sfc_sh.values(), bottom=list(sfc_lw.values()))\n",
    "ax[1].bar(np.arange(1, len(sfc_lw)+1), sfc_lh.values(), bottom=np.array(list(sfc_lw.values()))+np.array(list(sfc_sh.values())))\n",
    "ax[1].scatter(np.arange(1, len(sfc_lw)+1), np.array(list(sfc_sw.values())) + np.array(list(sfc_lw.values()))+np.array(list(sfc_lw.values()))+np.array(list(sfc_sh.values())))\n",
    "\n",
    "ax[1].set_xticks(np.arange(1, len(sfc_sw)+1))\n",
    "ax[1].set_xticklabels(sfc_sw.keys(), rotation=90)\n",
    "\n",
    "ax[0].set_title(\"TOA Balance\")\n",
    "ax[1].set_title(\"SFC Balance\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'LW_sfc, global & time average')"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAGXCAYAAACqWAsFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABCpUlEQVR4nO3deZxcVZn/8c+XhB0XkLBIgLAJorIGBFlEFoVhZ0RBwKAojhuLOoj+VMBRcVxBxsGJgkZEBBcIi4IQdoVAwh5AQVklkBBBNkGW7++PcyqpVKq7K50693Z1P+/Xq19VdW91P+d2VT117rlnkW1CCCH0nsXqLkAIIYTBiQQeQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj4oEHkIIPSoSeGhL0r6SHpL0jKRNC8UYJ8mSRnfw3OMl/WyQcQ6VdO1gfnchYqyR/1ejSsYJoVkk8CFO0v2Sdm56PDonii2bth2UE2HrtrsXIfS3gE/YXs72zYvwd4Y8SUtK+qmkJyXNknRyB78z3+ti+8H8v3q5bGlDmCcSeI+x/RJwHfD2ps3bA3e32Xb1IoRaE5ixCL/fSw4FNgPWyj/n1VmYoUBJ5IchLl6g3nQ1KUE3bAf8d5tt/SZwSVtKmibpKUmPSfpOro0+A4wCbpX0l/zc1SX9RtJsSXMk/U8nBZW0lqSrJT0t6TJJ3++rKUTS6yWdL+nvku6V9OGWpywl6ez8t26StHHT7x4r6S95352S9u2kfNlLwD9sP2H7WdtXDHBMZwBrABfks6FjWpuDJF0p6SuS/pifc4Gk10k6M/+/b5Q0rulvbiDp0nzsf5L0nn7if0DSXflY/yrpI0377pK0R9Pj0ZIel7RZfrxVLtOTkm6VtEPTc6+U9FVJfwCeA9buL1b+nWMkzZT0iKQP5f/BunnfkpK+JenB/P76gaSlB345Qsdsx88Q/gHuB3Zu2fZ24O+kL+AVgQeAZYDHmra9AqwxwN++Djgk318O2Kppn4F18/1RwK3Ad4FlgaWAbTss/3Wk5pglgG2Bp4Cf5X3jcpzR+fFVwP/mv78JMBvYKe87HngReDewOPAZ4D5g8bx/f+D1+fjfCzwLrJr3HQpc208ZNwZeBr482NelzbFcCdwLrAO8BrgT+DOwMzAa+Cnw4/zcZYGHgA/kfZsBjwNv6iP27vnvKr8XngM2y/u+BJzZ8ty78/3VgDnAv+X/0y758ZimMj8IvCmXY/EBYu0KPJqfvwxwRsv75iTgfGAF4FXABcCJdX+mhtNP7QWInwFeoPYJfCng+Zx49m18YIHrm7bd18Hfvho4AVixzb7mD+LWpGQ6eiHLvgapdrtM07af0SaBA6vnJPqqpueeCPwk3z8euL5p32LATGC7PmLfAuyd7x9KHwk8J5cHczKaChzftO9vwFs6eV1on8D/X9P+bwO/a3q8J3BLvv9e4JqWv/9/wHEd/p/PA47M99cFnm78z4EzgS/l+58Fzmj53UuACU1l7vdLrCXW6c0JOcd2vhXpS3Sdpv1bd/K+jJ/Of6IJpQfZfh64gdRksj1wTd51bdO2Ttq/DwPeANydT+n36ON5qwMPOLW/L4zXA3+3/VzTtocGeO7TTdseINUaF/hd268AD+ffQ9L7Jd2SmwaeBN5MOhMZyP7A/bYvBnYD3q3U42Uc6Uvijg7+Rl8ea7r/zzaPl8v31wTe2ih7Lv9BwCrt/qik3SRdn5tbniTVqFcEsH0vcBewp6RlgL2AnzfF2b8lzrbAqk1/fr7Xp79YpP/9Q3387hhSrXx6U6yL8/bQJQN23wpDVqMdfC3gR3nbNcDBedupA/0B2/cAB+aLVfsBv5L0OtvPtjz1IWANSaMXMonPBFaQtExTEl+9j+c+kp/7qqYkvgapFtww93dzmccCj0haE/ghsBNwne2XJd1CqgUOZDTpLAHbf5f0TuAPwIHAV5yrjm10cxrPh4CrbO8y0BMlLQn8Gng/MNn2i5LOY/5jPYtU/sWAO3NSb8Q5w3brtYVmc4+rg1gzSa9BQ/Nr+zjpS+pNtptfw9BFUQPvDYtLWqrpZzQpgb+D9KG5Mz/vWmAHUvvxgDVwSQdLGpNrs0/mze26wd1A+rB+XdKyuQzb5L/RuHg3rvWXbD8ATAOOl7SEpK1JTQcLsP0Q8EfgxPz3NyKdIZzZ9LTNJe2Xj/8o4AVSs9GypMQzO5fpA6QaeCd+C2wh6SOSFie1s/+RdGbySj+/9xiwdocxBnIh8AZJh0haPP9sIemNbZ67BLAk6VhfkrQb8M6W5/wib/so82rfkJqv9pT0Lkmj8v95B0ljaW+gWOcAH5D0xlzb/1JjR35P/RD4rqSVACStJuldnfxDQmcigfeG35JqM42f40lJ5jXA1EYt0fYc0odtVq5dD2RXYIZSr5OTgQNy88x8nPo270lq23yQ1HTx3rx7dVJTR1+1rINIbZ9zgK8AZ5MSbzsHktqSHwHOJbUBX9q0f3KO+wRwCLCf7Rdt30lqY76OlFjfQqpFD8j2faSmk/eTao3XkS6Obg98Q9KuffzqicAXcvPAZzqJ1U8ZniYlxgNIx/4oqVfRkn089whS8nwCeB/pQmHzc2bm43gb6f/d2P4QsDfwedL75CHgP+kjDwwUy/bvgO8BV5Au2F6XdzVe38/m7ddLegq4DFh/wH9I6Jj6PkMMYWCSvgDMtv1/HT7/bFKviOPKlixULZ8x3AEsOYjrJWEQIoGHoiRtQeryeB+plnkesLWH+ejOkUKpv/1FpGasScArtveptVAjSDShDHOSfqc0kKT15/MVFWEVUve0Z0in2x+N5D2sfITUHPMX0vWTj9ZbnJElauAhhNCjogYeQgg9KhJ4CCH0qEoH8qy44ooeN25clSFDCKHnTZ8+/XHbC4xirTSBjxs3jmnTplUZMoQQep6kB9ptjyaUEELoUZHAQwihR0UCDyGEHhUJPIQQelQk8BBC6FGRwEMIoUdFAg8hhB4VCTyEEHpULKkWQhjxxh17UfEY9399967/zaiBhxBCj4oaeAhhrtI10RK10JEsauAhhNCjogbegV5tHwuDE6936BWRwEOfIpGFMLRFE0oIIfSojhK4pKMlzZB0h6SzJC0laQVJl0q6J98uX7qwIYQQ5hkwgUtaDTgCGG/7zcAo4ADgWGCK7fWAKflxCCGEinTahDIaWFrSaGAZ4BFgb2BS3j8J2KfrpQshhNCnARO47b8B3wIeBGYC/7D9e2Bl2zPzc2YCK7X7fUmHS5omadrs2bO7V/IQQhjhOmlCWZ5U214LeD2wrKSDOw1ge6Lt8bbHjxmzwJqcIYQQBqmTJpSdgftsz7b9IvAb4G3AY5JWBci3s8oVM4QQQqtOEviDwFaSlpEkYCfgLuB8YEJ+zgRgcpkihhBCaGfAgTy2p0r6FXAT8BJwMzARWA44R9JhpCS/f8mChhBCmF9HIzFtHwcc17L5BVJtPIQQQg1iJGYIIfSoSOAhhNCjIoGHEEKPigQeQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj4oEHkIIPSoSeAgh9KhI4CGE0KMigYcQQo+KBB5CCD0qEngIIfSoSOAhhNCjIoGHEEKP6mRNzPUl3dL085SkoyStIOlSSffk2+WrKHAIIYSkk1Xp/2R7E9ubAJsDzwHnAscCU2yvB0zJj0MIIVRkYZtQdgL+YvsB0kr1k/L2ScA+XSxXCCGEASxsAj8AOCvfX9n2TIB8u1I3CxZCCKF/HSdwSUsAewG/XJgAkg6XNE3StNmzZy9s+UIIIfRhYWrguwE32X4sP35M0qoA+XZWu1+yPdH2eNvjx4wZs2ilDSGEMNfCJPADmdd8AnA+MCHfnwBM7lahQgghDKyjBC5pGWAX4DdNm78O7CLpnrzv690vXgghhL6M7uRJtp8DXteybQ6pV0oIIYQaxEjMEELoUZHAQwihR0UCDyGEHhUJPIQQelQk8BBC6FGRwEMIoUdFAg8hhB4VCTyEEHpUJPAQQuhRkcBDCKFHRQIPIYQeFQk8hBB6VCTwEELoUZHAQwihR0UCDyGEHhUJPIQQelSnK/K8VtKvJN0t6S5JW0taQdKlku7Jt8uXLmwIIYR5Oq2BnwxcbHsDYGPgLuBYYIrt9YAp+XEIIYSKDJjAJb0a2B44DcD2v2w/CewNTMpPmwTsU6aIIYQQ2umkBr42MBv4saSbJf1I0rLAyrZnAuTbldr9sqTDJU2TNG327NldK3gIIYx0nSTw0cBmwKm2NwWeZSGaS2xPtD3e9vgxY8YMspghhBBadZLAHwYetj01P/4VKaE/JmlVgHw7q0wRQwghtDNgArf9KPCQpPXzpp2AO4HzgQl52wRgcpEShhBCaGt0h8/7JHCmpCWAvwIfICX/cyQdBjwI7F+miCPbuGMvKh7j/q/vXjxGCKH7Okrgtm8BxrfZtVNXSxNCCKFjMRIzhBB6VKdNKCFUKpqOQhhY1MBDCKFHRQ08hDAklD7rGo5nXFEDDyGEHtUzNfBoEw0hhPlFDTyEEHpUJPAQQuhRkcBDCKFH9UwbeAgjRfTGCJ2KGngIIfSoSOAhhNCjIoGHEEKPigQeQgg9KhJ4CCH0qEjgIYTQozrqRijpfuBp4GXgJdvjJa0AnA2MA+4H3mP7iTLFDCGE0GphauDvsL2J7cbKPMcCU2yvB0xhIVaqDyGEsOgWpQllb2BSvj8J2GeRSxNCCKFjnSZwA7+XNF3S4XnbyrZnAuTbldr9oqTDJU2TNG327NmLXuIQQghA50Ppt7H9iKSVgEsl3d1pANsTgYkA48eP9yDKGEIIoY2OauC2H8m3s4BzgS2BxyStCpBvZ5UqZAghhAUNmMAlLSvpVY37wDuBO4DzgQn5aROAyaUKGUIIYUGdNKGsDJwrqfH8n9u+WNKNwDmSDgMeBPYvV8wQQgitBkzgtv8KbNxm+xxgpxKFCiGEMLAYiRlCCD0qEngIIfSoSOAhhNCjIoGHEEKPigQeQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj4oEHkIIPSoSeAgh9KhI4CGE0KMigYcQQo+KBB5CCD0qEngIIfSoSOAhhNCjOk7gkkZJulnShfnxCpIulXRPvl2+XDFDCCG0Wpga+JHAXU2PjwWm2F4PmJIfhxBCqEhHCVzSWGB34EdNm/cGJuX7k4B9ulqyEEII/eq0Bn4ScAzwStO2lW3PBMi3K3W3aCGEEPrTyar0ewCzbE8fTABJh0uaJmna7NmzB/MnQgghtNFJDXwbYC9J9wO/AHaU9DPgMUmrAuTbWe1+2fZE2+Ntjx8zZkyXih1CCGHABG77c7bH2h4HHABcbvtg4HxgQn7aBGBysVKGEEJYwKL0A/86sIuke4Bd8uMQQggVGb0wT7Z9JXBlvj8H2Kn7RQohhNCJGIkZQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj4oEHkIIPSoSeAgh9KhI4CGE0KMigYcQQo+KBB5CCD0qEngIIfSoSOAhhNCjIoGHEEKPigQeQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj+pkUeOlJN0g6VZJMySdkLevIOlSSffk2+XLFzeEEEJDJzXwF4AdbW8MbALsKmkr4Fhgiu31gCn5cQghhIp0sqixbT+THy6efwzsDUzK2ycB+5QoYAghhPY6agOXNErSLcAs4FLbU4GVbc8EyLcr9fG7h0uaJmna7Nmzu1TsEEIIHSVw2y/b3gQYC2wp6c2dBrA90fZ42+PHjBkzyGKGEEJotVC9UGw/SVqVflfgMUmrAuTbWd0uXAghhL510gtljKTX5vtLAzsDdwPnAxPy0yYAkwuVMYQQQhujO3jOqsAkSaNICf8c2xdKug44R9JhwIPA/gXLGUIIocWACdz2bcCmbbbPAXYqUagQQggDi5GYIYTQoyKBhxBCj4oEHkIIPSoSeAgh9KhI4CGE0KMigYcQQo+KBB5CCD0qEngIIfSoSOAhhNCjIoGHEEKPigQeQgg9KhJ4CCH0qEjgIYTQoyKBhxBCj4oEHkIIPaqTFXlWl3SFpLskzZB0ZN6+gqRLJd2Tb5cvX9wQQggNndTAXwI+bfuNwFbAxyVtCBwLTLG9HjAlPw4hhFCRARO47Zm2b8r3nwbuAlYD9gYm5adNAvYpVMYQQghtLFQbuKRxpOXVpgIr254JKckDK3W9dCGEEPrUcQKXtBzwa+Ao208txO8dLmmapGmzZ88eTBlDCCG00VECl7Q4KXmfafs3efNjklbN+1cFZrX7XdsTbY+3PX7MmDHdKHMIIQQ664Ui4DTgLtvfadp1PjAh358ATO5+8UIIIfRldAfP2QY4BLhd0i152+eBrwPnSDoMeBDYv0gJQwghtDVgArd9LaA+du/U3eKEEELoVIzEDCGEHhUJPIQQelQk8BBC6FGRwEMIoUdFAg8hhB4VCTyEEHpUJPAQQuhRkcBDCKFHRQIPIYQeFQk8hBB6VCTwEELoUZHAQwihR0UCDyGEHhUJPIQQelQk8BBC6FGRwEMIoUd1sqTa6ZJmSbqjadsKki6VdE++Xb5sMUMIIbTqpAb+E2DXlm3HAlNsrwdMyY9DCCFUaMAEbvtq4O8tm/cGJuX7k4B9ulusEEIIAxlsG/jKtmcC5NuV+nqipMMlTZM0bfbs2YMMF0IIoVXxi5i2J9oeb3v8mDFjSocLIYQRY7AJ/DFJqwLk21ndK1IIIYRODDaBnw9MyPcnAJO7U5wQQgid6qQb4VnAdcD6kh6WdBjwdWAXSfcAu+THIYQQKjR6oCfYPrCPXTt1uSwhhBAWQozEDCGEHhUJPIQQelQk8BBC6FGRwEMIoUdFAg8hhB4VCTyEEHpUJPAQQuhRkcBDCKFHRQIPIYQeFQk8hBB6VCTwEELoUZHAQwihR0UCDyGEHhUJPIQQelQk8BBC6FGRwEMIoUctUgKXtKukP0m6V9Kx3SpUCCGEgQ06gUsaBXwf2A3YEDhQ0obdKlgIIYT+LUoNfEvgXtt/tf0v4BfA3t0pVgghhIHI9uB+UXo3sKvtD+XHhwBvtf2JlucdDhyeH64P/GnwxV1oKwKPVxgvYkfsiB2xS1jT9pjWjQMuatwPtdm2wLeB7YnAxEWIM2iSptkeH7EjdsSO2MMldrNFaUJ5GFi96fFY4JFFK04IIYROLUoCvxFYT9JakpYADgDO706xQgghDGTQTSi2X5L0CeASYBRwuu0ZXStZd9TSdBOxI3bEjthVGPRFzBBCCPWKkZghhNCjIoGHEEKPigQeQgg9KhJ4F0jaJt8uWXdZqjSCj3tKvv3vustSJUlH5ttt6i5L1STtn2/XqrsszYZdApe0lKRPSfqNpF9LOlrSUoXDfi/fXlc4Tp9G6nEDSFpF0l6S9pS0SgUhV5X0dmAvSZtK2qz5p4L4SFpb0gWSHpc0S9JkSWsXDvuBfHtK4Tj9kvQ6SadIuknSdEknS3pd4bCfy7e/LhxnoQy7XiiSzgGeBn6WNx0ILG97/4IxrwfuAnYnzQkzH9tHlIrdVIaRetwfAr4EXE4aHfx24Mu2Ty8Y893AYcC2pPEQzaOSbXvHUrGbynA9aTK5s/KmA4BP2n5rwZhnAVsDY4C/NO8iHfdGpWK3lONS4GrmvdcPAnawvXPBmJeRuktvAlzTut/2XqVi92c4JvBbbW880LYux1wR2Bn4b1IymY/tSaViN5VhpB73n4C32Z6TH78O+KPt9QvG3Mb2HyR9yfaXS8UZoAxTW5O1pOttb1U47iqksR8LJCzbD5SM3VSG6bY3b9lWdGh7Hqy4GXAG8KHW/bavKhW7P4syF8pQdbOkrWxfDyDprcAfCsf8T9uflbRGFUmrDyP1uB8mnXk0PA08VDjm94DNgX2AWhI4cEWeg/8XpDmI3gtcJGkFANt/73ZASVNs7yTpkqqSdR+ukHQAcE5+/G7gosIxT7N9iKQf1pWs2xmONfC7SLMePpg3rUE6zX+FQqd5km4nfTtPtV1JG2ibMozU4/4p8BZgMimR7Q3cAPwZwPZ3CsRsNB39G3B26/6Kmo7u62e3bXe9PVzSncBHgR8A76NlQjvbN3U7Zh/leBpYFng5l2Ex4Nl5xfCrC8S8k7T2wfnADix47F3/wuzEcKyB71pDzItJU0suK+kpcptg47bEG6qNkXrcf2H+9tjJ+fZVBWPuQWo62hGYXjBOn2zX0RviS8CxpInrWr8YTfp/FGe75Gvblx+Q3u9rk17z+a575O2VGzY1cEmvtv1U4xSyVRXfkJIm2650UYuRetxDgaSNbd9accwdbV8uab92+23/poIyfNH2f5WO0ybuBrbv7qunTxVnAJJOtf3R0nE6NZwS+IW298inlo1aYEORU8o+yrEysEV+ONX27MLxRupxn2T7KEkX0H4e+uK9AiSNJXWp2yaX4VrgSNsPF4x5gu3jJP24zW7b/mCp2C3l2AvYPj+80vaFFcScaPtwSVe02V1J759cjo2B7fLDq23fVkXctmUZLgl8KMid/b8FXElKpNuRLvT9qs5ylVbHcUva3Pb03B97AVVcaMrd2X5O6pkAcDBwkO1dSseuk6QTSUsqnpk3HQhMs/25vn9reJB0BGmFscaZzr7ARNu19I0flglc0kbAOJra+Cs6tbwV2MX2rPx4DHBZya58LfFH5HHXpY+um7fY3qSC2K8F3s+Cr3cVF1BvAzax/Up+PAq4ucJ+4KNIYw/GMf+xd/2CdZvYtwFb2342P14WuK6qY2817C5iSjod2AiYQeqBAen0tngiAxZrJLFsDhWNdh3Bx70H8F/AmqT3c5UXUGdLOph5g2kOJB17FX4LXA/czrzXu0qvBRrXV15TcewLgOep59hF6v3S0OgJU4thl8CBrWxvWFPsiyVdwrwP9HtJH7QqjNTjPgnYD7jd1Z9OfhD4H+C7pC/LP+ZtVVjK9qcqitXqRNK4gytIyWt75g01r8LYumq8wI+BqZLOzY/3AU6rqSzDrwlF0mnAt23fWVP8/UhDrEW6wHHuAL/Srbgj9bivAHZqnM5XSdJStp+vOm6OfTTwDHAh8EJje1X9kSWtSrpoLdJF60eriJtj/zcwxfbvq4rZEn8z5n+v31xHOWB4JvDtSadYj5Le2JXN0yDpg8A1tu8pHatN7JF63FuQmlCuYv5EVkV76L3AY6S5Ma4G/mD7H6Xj5tgfB74KPMm8XjiV9DqSdAbpeK+xfXfpeG3i70uaB2Ux4EUqbDaT9GXS6/3HRjt4nYZjAr8X+BQt7WNVDP3NL+62pPbY6aQX+hrbt1QQe6Qe9+9JNdHW4z6hdOwcfw1Sr5ttSCMzn6zoIuZfgLfafrx0rDaxdyS93tuRBrDcQqqJnlxR/L+Smi4qbzbLlZVtSZN6PU3+8rY9ud9fLFWeYZjAL6+qP2g/ZVga+DDwGWA126MqiDlSj7voJEYDxB5LSmJvBzYmXdS71vaJFcQ+HzjA9nOlY/URfxSpCeUdwH8A/7S9QUWxLwF2q6PZrKkMqwDvIb3Xl69pdOiwvIh5t6Sfk5oTmk+pq+hO9wVSTWw54GbSi7vA1JOFjNTjvkzSO2tqD32QNJ3s12z/R8WxXwZuydcAml/vKroRTiHNRXId6XXeoqUXUmkzgSsl/Y7qm81+BGzIvKazdwOVzAHTznBM4EuTXtR3Nm2rqjvdfsBLpJnRrgKur/Ai10g97o8Dx0h6gYrbQ4FNSafT78szA94DXGW7il4J5+WfOtxGmo3xzcA/gCclXWf7nxXFvy//LJF/qvQ60rzgT5LOuB63/VLFZZhr2DWh1E3Sq0gf6m1Jp1iP2d623lKVN4KPeznmtQcfTPryGFdroSqSj/0DpDOuVWyPmKX1JL0ReBdwNDDK9tg6yjHsauBKy4gdBrwJmLukWBVzREh6M/PaRMeT5qWupClhpB53jr88sB7zH/fVFcSdBixJ6v99LbB9FReNc+z1SP2xN2T+466iF8onSK/35sADwOlU+3qPAY5hwfd6FSsh7UE69u2B5UkrQVV27K2GXQInzUtxN+nb8cuk5Zbuqij2f5O6V30PuNH2ixXFhRF63EpLqh1JmuL0FmArUttsFRd0dys9aVc/fgwcRxpE9A5STbiqEYFLk6aTnV5T88GZpHnY9yBdQJ0AVPU67EZ6r59s+5GKYvZp2DWhSLrZ9qaSbrO9kaTFgUvq7qFR2gg+7ttJvSGut72JpA2AE2y/t2DMfkdAVnQxbbrtzSXdbvsteds1trcb6HcXIWbbKYsbKhxE1Dj22xrjHCRdZbvtxGbD2XCsgTdqf0/mU/tHSZPeFJOTSLtvwioXex2px/287eclIWlJp/mii62HmdXSZazF85IWA+7JTRp/A1YqHHM6817v1tp+lYsaNN7rMyXtDjxCOgMrRmkVoNb3euN/UNVF8wUMxwQ+MbeJfoG0/NFywBcLx9yj8N/vxEg97ofzzHznAZdKeoL0gS6mqkFCAzgKWAY4gjQS9R2kpoRiXM8qQO18RdJrgE+T5mN/NeliYjF19fMeyLBrQhmIpAmuZrX0PVzBJPedGgnHrTQ3+GuAi23/K29b3vYTFcS+yTWtC9qOpFNsf7KCOMfbPr50nIUh6XMVDaY63PbE0nH6U8mUn0PMkRXFqWu18r4M++O2fZXt8xvJO5tSUfjaphTtwzYVxSm+8tEg7F9RnKoHby1gJCbwqj5oQ+0DHcdd1kUVxRlqhtrrDSPovT4SE3hVbUYfqShOp+K4C5C0itL6kFPz/BgjzZBpNmpS1Xu99mtAIzGBF/vWlLSUpE9J+g3wWUlH5wE2Q8FIPe5ich/0G0hTCbwbuD7PVjcUlHy915Z0gaTHgcckTZZUVQ+UTpQ89tdJOkXSTcD5kk6W9LpS8QYyEhP4Hwr+7Z+SRoedQlqp5Y3MW/C2biP1uEue5v4nsKntQ21PII1M/GzBeP2StGbTw5JTu/4cOAdYBXg98EvmrcY0FPyy4N/+BTAL+HfSl/Zs0qCiWgyrXih5EMdqpBVCnmnavqvtiyuI326R2wW2VUXSB2z/uII4Q+24H7S9Rr6/QqkBJnlWvt2aerwsAfzW9s4l4jXF3Zr0Pr/a9iylxayPBbazvXrJ2Dn+VNtvbdl2ve2tKoj9LtJc4KuRmkoeASZX8fnO8afb3rxlW21TGg+bGrikI4DJwCeBOyTt3bT7axUV42ZJc9/Ekt5K2ZrvQKrqrzzUjnturbvw6MC/kdq+j5d0HGmR4Xtzc1KR9SolfZM098i/AxfluJcCU0nzwVThCknHShonaU1Jx+SyrDDQaM1FIekkUm+qq4BvAN/M94+QVMliEqRjP0DSYvnnPdR4AXvY1MDzqMCtbT8jaRzwK+AM2yc3hplXUIa7gPVJ80QDrEGaj+QVCo1MlHRbX7uAN1QxQ1wdxz1AeebWwAvHOa6//SUG/Ei6E9gsjz5dnlQD3cgVLmcn6b5+drvUhFqS/mz7DW22C/iz7eJfYHlE5rLMW41+MaCxtFrlIzKH00jMUY1mE9v3S9oB+FVuF6yqu8+uFcVptjJpAqvWwSoizZJXhcqPu58arkijUIuraUTmPxtzrdt+QtKfqkzeOW5dIzKfl7Sl7Rtatm8BVDL//FAbkTmcEvijkjZxXocx18T3IJ1uvqVkYEmvtv0UaY28BRQ+jb8QWM5t1p+UdGXBuHUfd38fpKKn05JOsn2UpAto02XNdsnBLesoLafWMK75ccnYkna0fbmk/drtd/nVnw4FTlWae/7hvG114Km8rxhJG+R5dtp2m7Rdy6o8w6kJZSzwku1H2+zbxnaxNllJF9reI59amvlr/MVOKes2VI9b0lG2Tyr49ze3PT0P3V+A7asKxu53xr3CsU+wfZykdhfG7Qrmns/lWIV0EVPAw+0+8wViTrR9uNISdq1c16yfwymBD4mpLoeCoTBHQ52qagMfaiSdXXIa3aFIQ3AulioNpwT+Cum0qjHBfC21wdylaxxNzVMVnFq2lqHyiZWGwnE3leWhirrT7UGaCXBN0nFXuR5nu/JUdfH2tcD7WfD1Lr6gcpuyVPpelzQK2J0Fj734HPDtDKc28FOAHUjd184CrnXF306STgc2AmaQemBAdQsLz1eUSoMNneNuqOp1P4k0CvP2qt9rNfstqcvk7cx7vetS9XwkF5AumA6FYx8+NXCY251oB+BAYEvg98Cptvvr9tTN+Hfa3rCKWAOUY6zthwd+ZtfiVX7caj/BPqQP9NK2i1dOcnvoTrYr+yD3dRGNdNwX2l61gjIMmalzJS1W8f9/7ipAQ8FwqoGTa0FXSLoZOIB0ensP8MOKinCdpA1t31lRvMagmbtsPyVpadKIvM1yf+Gv2f5HBcWo/LiHSHeuY4DfSroKeKGxsfDp9Lf72Xd3wbjNzpD0YVIPqObjLn6dKY/EHAtMsX1/I3lL+qDt00vHB34n6Z22f19BrAENmxq4pGWBvYH3AmNIp+9n236owjJsTzrFepT0xi6+tJikGcDGtl+SNBF4jjSIaae8vW2Xry6XofLj7qMclV68lfR74BlaTqdr6h9eGUkfB74KPMm8s6Di15kkfQ3YFrgJ2BM4yfYpeV8lZwWS9gV+RhrA8yJ1X/cYRgn8WVJt+yzgXlpOr6u4oCbpXuBTLPiBfqBgzLtsvzHfn+9NLOkW25uUit0Up/Lj7qMcVV/Qqm0OjJZyTLR9eIXx/gK81fbjVcXMcW8nTR72Ur6Q+nPgT7aPrnC09V9Jc7EMiesew6kJ5ZekpL1B/mlW1QW1B22fP/DTuuoOzZu06lZJ421Pk/QG5i3+Wlodx91O1Re0Lhsip9NVf4nMIJ3pVW207ZcAbD8paU/SWrC/BJaoqAz3AHcMheQNw6gGPhRI+l/gtaTmhOa2wWJfHkqLu54MbAc8Tppg/6H8c4TtW0vFbipD5cfdRzmqvnjbmBfjBWo8nZZ0se3KpjOQdC5p+uArmP/1LtqNUNKFwDdbBytJ+grwedvFJ+eT9BNgbeB3VHfdo+/yDJcEnr+Nb2uctkv6EmnGtgeAI6voiVLnCLU8vHht0lnVw7YfKx2zKXatI/NaqaJpdIcSSdsA77P98QpiTWi33YUXzc4X6bH9zzb7VrP9t5Lxc5y2E5jVdd1jOCXw24CtbD+XB1h8h9SdcFNgf9vvqrWAhcQI1AVVORJTaUbA9YC5KxDZvrqi2JsA7wPeA9wH/KZxUW846qcLJVDffCR1Gk5t4LbdaJfbDzjN9nRguqSPVVEApWXEDiOdXjZ/oEvWRB+nnxGopFp5UXUct/qfRnflUnFbyvAh0vzUY4FbgK2A64Bi82LkaxsHkConc0irwcj2O0rFbFOG9YATgQ2Z//Uu/V6bRmp/n90oStM+U/D/3iBpDKn7aOt7vZa5UIbNgg6kcTzLSVqM1IVuStO+qtZnPIO0zNS7SBPNj6WPmfq66BTSVLIXAxOAtW2vlX+qmkyqjuNemTSce882P3MKx244kjSV6QM5gW7KvORSyt2k9/eetrfNNe6XC8ds9WPgVFKl4R2kJfWqWELv08A/gH/mMuxp+x35p6oEeibpNViLtGDK/cCNFcVekO1h8QN8kNR98Cbg4qbtm5I6/VdRhpvz7W35dnHg8griivRBmkiqCX4DWKvC/33lxw2cBmzbx76fV3TcN+bbW4AlG/cLx9yXVOt+iDRAbSfgvqpe61yG6fn29qZt11QYfy3g86RViM4BNqnh2G9r2nZVlf//5p9h04Ri+3RJlwArAc09Lx4FPlBRMRrd9p6U9OYce1zpoE7vojpHoFZ+3LYP62ff+0rGbvJw7o98HnCppCdIK+QUY/tc4Nw8cG0f4GhgZUmnAue6mi6Nz+cz3XskfYK0tNxKFcQFwPZ9kiYDSwOHAG8gfYlWofFenylpd9LrPbai2Auq65ujwDfjwU33t2nZ94mKyvAhYHlge+CvpNWrP1I45rKkC1mTSSvwfAZYveL/feXH3Uc5Vqg6ZlPstwN7AUs0bVu+quMGDqeCs70cbwvSqkdjSU0ZvyZ1ICgdd23m1bx/RVoVfqmKX+c9gNcAbyZ1o5wO7FVlGZp/hlMvlLmj8NqMSBwSk+9ImuAud7UaCiNQB1LouLcBfkQa+flB4CvAOqTmm/fYvq6b8QajivedpHWBjUnz4VQ2F01/JJ1i+5MF/u4rwG2kyspTLPher6UvdjNJn7N9YlXxhtVFzD7ut3tclyML/M1fAjeTRp/uwfwX8/YoEG8wShz3d0nd5z5EWhX8BKeLtnsD3yoQbzC6/r6TdIWkFfP9Q0hTu+4GnC2p60lzkLYp9He/DJxL+tJejrSsXvPPULB/lcGGTRs4838bt55WDJXTjK5/oG0f2u2/WUCJL9DFbd8OIGm27Wsh9QVuDPgYAkq878Z43hwkRwBb254jaRnSHN3Dth+4e2PlnUori8OpBr6BpNvyhDeN+43H69dduKzrH2hJn5K0wAU9SZ+UdFS34w1SiUTW/N79XMu+qubFqMOLklbL958Bns33XwBG1VOkakj6cO6DjpLTJf0jf86LT2TVoUori8OpBr4xqW9w6/Sxa1K4Z8BCKPHt/EHS/CetJpL6p55UIObCKnHcX5S0jO3nbJ83N5C0Dqlf8lBQ4riPBn4v6dekQS2XS7qYNBfOUJk+oFQt9EjgJ/n+gaTP/NqkrsLfI/0P6hY18EH6LvCU7Qeaf0izpn235rI1/KHA37Ttf7XZ2JiXeygocdy/J/XAafU06cNcC0kPNj3cqdt/3/aVwNuAmaQubdNJte9P2q6t7V/Smk0PTy4U5iXbjW58ewA/tT3H9mW0fy/U4ZdVBhtOvVDusP3mPvbdbvstheNvAKwGTLX9TNP2XW1fXDDu7cDObpm8StLKwGWlj7ufchWdUEpp8YqLW3vZSDqINMDno6ViD1CuShZUroukrUnv86ttz1JazPpYYLvSxy3pJtKCwk+QJqnb0faMvG/uvPgF47+L1Pd+NVJTySPA5JKf7wHLNIwS+L22113YfV2KfQTwceAuYBPS7IeT876iXckkvZ90MevTpFGoAJuTRmN+v9vd9xaiXEUnlFI/63BKmmH7TaVi96eC4+533nXbexWM/U1SzfcWYF3SkmofA74G/J/t50vFzvH3AP6P1NZ/ge0P5+1vB46xvXvB2CeRBgz9lDT3EKR+8O8H7rFdoqfVgIZTG/iNkj5se77Rh/kC3/TCsT8MbG77GUnjgF9JGmf7ZAo3Y9j+qaTZpC5WbybVDGYAx9n+XcnYNU8o1d//tWjToKRP9bWL1L2tpK1J13nOIg1oqbKZbHfSijjP51kYHwE2sn1PFcFtX5ibal5l+4mmXdNISymW9G+239C6UdLZwJ8p01V2QMMpgR9FGmJ8EPMS9nhSj4R9C8ce1Wg2sX2/pB1ISXxNKviA5URdNFn3YWXSBFZPtGwXaVRoSbMkbWn7hvkCS1tQfkKp/vocl2r/bVgF2IV0Ee99pD7wZzWaEgr7Z6OWbfsJSX+qKnkDSDrG9jeAJyTtb/uXuSzPKq2X+fmC4Z9v934jjUoteubRn2HThNIg6R2kmijADNuXVxDzcuBTtm9p2jYaOB04yHal3buqGnkq6TTgx40+2C37fu6Cc5JI2pI0kdFPmP8L+/3AAbanloo9QLmOsn1SRbGWJCXybwJfduG5wCU9CTTPdb598+OSzTc5fm2jrZXmIj+V9OXdaEJZnTQi9GNOU1dXbtgl8DpIGku6Qv5om33b2C7RC6O/8tzsChZ4rVu+UPsxmr6wgf+xPavGMhVfTCIn7t1JyXsccD5wuguvSJPbmvvklqXOCsSf+75ufY9X9Z6XtArpIqZIK18t8Jmv0nBqQqnTc9Dn6jh3VVwWSKfVtZB0uO2JVcTKPW/aLnFVo6JNZpImkb6wfkeaPuCOkvGa9Zegc1tw0QTOEBhtnRP2owCSjgeOryJuX6IG3gV5kp0+V8VxdQsr1K7C5pu9gbG2v58fTwXG5N2fbbSPVq2CXiivMG/0ZfOHt5YFlecGr+bM42XSsYs0lWxjBS6RZiVcvGT8NuWpfZK8qIF3xynADqQBK2cB17qib0alldEX+CBT3we6ql4Rx5DmPm9YknRBaVnSiMRiCbzN/3zuLlJiKcYVrLw+VFV9LakDtQ+UiwTeBbaPlCRSEj8EOEXS74FTbd9XOPwUUs+E3wC/sP3gAM8vbc+K4ixhu3nahGttzwHmKC12UIzt2ma+y71sVmztIippT+CRkhfT1PeiwiJN41u5Kpvs2ti8prhzRRNKlymt0NJYFefzrf3SC8V8DWkh5wNI63+eTUrmRVekl/RW0jzUTynNAHgsaV6WO4Gv2f5Hwdj9Ddz6i+11SsXuI2YliUTSlcChtu9v2b4uMNEF14aUdEV/+13hwsoNVTZj5JGYY0lLNN7ftP2Dtk+vogwLlCkS+KLLNb69SYMJxpBqw2e31BCrKMdiuQynkBJo0QnuJc0ANrb9Uh7a/hxppZSd8vb9CsY+E7iyzcCtjwA72D6wVOw+ylNV23+f00JIutX2xqXLMJRU2Pvka8C2pNHOewInNbpt1tkWHk0o3TGLBVfF2SKf7hZfFUfS20hdyrYDrgX2tX1NyZjZYrYbF27HN72Jr5V0S+HYRwPnSXof808hsCRpvoqqVdUe2l8be+UTOkmaaPvwquM2qarJbk/SKNSXcu+Tn0ta2/bR1NgWHgm8O35JStob5J9mJtXIi5B0P/Ak8AvSuogv5e2bQVrgoFRs4I6mSatulTTe9jRJb2De4q9F5L7eb5O0I9CY9+SiKgZu9aGqRHKZpK8CX2i+UC7pBKCOYx9fVaC+pjBIl5+KL6k2ulFZsf1kvuYwUdIvqXH++WhC6XG5TbTxIjZ6nzS4cJvoa0hDx7cDHie1fz+Uf46wfWup2H2UZ5W6B1ao/CyMy5LWAt2SeSuxb0Ka+/1DbpoJswqSLra9a0Wx+u3zb/uEgrEvBL7Z2hde0ldI17pq6R0UCbwL8rfxbU7zjyPpS8C/k6a8PLKCnii1kvQq0sT6o0mj0x4b4FdKlaP2frlV9IfOcdZm3pnHDNt/LR2zn7JsA7zP9sdrLMOytp8d+JmD/vtLA9j+Z5t9q5UeBduXaELpjq8CWwGNKS8PJrVJbwr8gDThU6Uk7UKaYnOXgjGaR54+1Lq9dC+YdkWqJEi9szACkBP2X3N5jqfiEYGSNiFNpvUe4D4KNhO2xF0NWJVUYfqXpJVIE9kdCry+YOg35vh97Y8E3sNsuzEqbD/gtNwfd7qkj5UMnNuAf0B6855Hmpv5p6Rk8tWSsUnNJn2OQCXVyqtUvMtmVucsjO3sRQUJPF/bOIBUOZlD6q6qqroPKq3x+v9IHQWWlHQy8B3S+710n+xppLl2GjNdtr7XizVV9icSeHdI0nKkbnQ7Af/btG+pwrG/Tbp4eR2wG2ll8i86zUVeWm0jUPvwi4riXAgs56bZJxvyNYmqVdUL4m7gGmBP2/cCSDq6otiQ3ufr2/67pDVIiXx729dXEPvTpGbRf5LeZ+dWfb2hnRE7LLfLTiJdUJpGGtgyDUBppeyZhWPb9pW2X3Ba3Hd2Rckbp1VINiH1wjkEuFnSNyStVTq2pG0k3SVphqS3SroUmCbpIaVlv4qxfZjbTKGb9xWbQrcfmwNIKr0m5r+TJnK6QtIPJe1EtV3onm80y+URx3+uKHlj+7u2twU+QZpGdoqkc3JTUm3iImaX5La5lYBbbb+St60KLF5yeLukvwKfadr0rebHpfugN5XjtVQ4AlXSDcBhpBVwLgD2sX1t7j55iu1tSsZvU54Vamjzby1DVRdQlyX1tT+Q1HQwiVQj/X3huLOY/yzrgObHto8oGb+pHG/KsQ8hXWc6p4q4bcsSCXzRSTrY9s/y/fnm/5b0Cdv/UzB2f13WbPuDBWPXNgJV888NPd+CtqV7o+ReFz8CXgE+CHwFWIc0H8h7bF9XKvYA5ap8QeV8wfrdpEU0irYDS5rQ334XXP819/o5gPR+f4j0xXGhC68DOmC5IoEvOtW4UkidJD3LgiNQ5ypZ+28eNi5pn9x81Nh3h+039/nLix67ttq/2s85D6kp41bbY0vFblOWdYGNSc2Gd1YVtw5K0/jeBkwmrcLT+l4vOm1FX+IiZneoj/vtHnc3cFqVvi+2fUbB8LWNQAW+KGkZ28+1JO91SL0SSlrc9u053uxGe7jtmxr9hQuazoIDthr+VTKw0mRW+9t+XNIhwBdJS6odn4fUl17S7fz+9rvskm5fZl7SLr1wdccigXdHnSuFbNFmm0hDu1cDiiVw24eW+tsdxG77Ybb9F+AbhcM3X/z/XMu+osOqbRe/QNyPMbYfz/ePALa2PUfSMqTeT0UTOLA1qfniLGAqFV5AtX18VbEWRiTw7tggD+4QsE7TQA9RuC+07U827iuNMjgI+CzpA1W0H3iem+Iftk9r2f5JYJQLLu5bc22sttq/pDWBJ52n6lVaxHsf4H7g+7ZL1sJfbBp1+AzzVgZ6AahisYVVgF1IF0/fR1o68CzbM0oHlvRh0uyX9+TP2WnMG209wfbNpcvQtlzRBr7oJK1HGtzRevFuTdIk+/cWjj+aNBLt06SayYm2/1QyZo57B7BZa9JQWnT3RtsbFYw9m35qY61zVnQ59lLAq2zPbtm+EvBUyQtbSkvH7Wv7kdyF7TLgRGAj4EXbHyoYewfg+8CvgRVIc99cTJoL5xLbpbsxNpdlSVIi/ybw5Qqab+4gzUb4otIMmJ8G3kkabX2c7e1Kxu9L9APvju+SPrgPNP+QBvZ8t2RgSR8nLaCwObCr7UOrSN6Z29X4bL9A+dPbVYDPkxb4PZlUM3vc9lUlk3f2PVLSarULhV9vYGnbj+T7B5NWo/828AHSBFfF2L4SeBtpbMOLpPb4F4BPVpW8JS0paT/gZ8DHSa9FFV1lX7LdmGFzD+CntufYvowapvFtiBp4F/TX60H9TMDfpdivkOYjn02b9vbCteDbgZ3dMnmVpJWBy0oed0u8qmtjd9resI99M2y/qd2+LsWe+36SdBPwOduX5Me3lXy96yZpEukL+3ekFafuqDD2TcDupOkTHgB2bDTdtHZjrVK0gXdHf8PlS/dKWJ9+mm8Kx/4mcJGkTzP/ogrfIA0oKion7t1JyXsc1dXG+ju7KH1We7mkc0i14OXJc4DnQWOle6HUed0B0sCZZ4E3AEc0TSxVxQLeXyKNtB4FnN+UvN9OnlSsDpHAu+NGSR9uHX0o6TDSaWZJ3yWNfHygJfaYvK/YQgO2f5rbor9MqhmZNOHPcW5ZdLfbWmpjJ1RZGwNmSdrS9g0tZdqCeZMdlXIUaeDUqsC2Taf1q5Ameiqptl4gAK5pzu0c+8J8AflVtpsnMZtGej1qEU0oXZCbDM4l1YAaCXs8qUvZvi64yECdzTd1yk1HjV4QzW/i4rUxSVsC5wA/Yf7X+/2kEYlTS8XO8fcB1gVubzSfVEHSKOb1AtmICnuB1E3SMba/ke/vb/uXTfu+ZvvztZQrEnj35C5djWQ6wxUs76X+V2fvc1/B8gzbkafN8pf2x2h6vYH/cVrqrWTcU4ENSdPW7gRcYPu/SsbsoxyVXneo21AdbR1NKF1k+wrgiorD1tl8005Viyrs2PiClLSWm1Y9krSfC0/ilS/c9rvEVyHbARvbfjkPoLmGNIFYJWq87lC32kZb9ycSeO87CjhX0kG0ab6poTwXVRTnW6R+yJD6JTfXgL5A2YWk9wbG2v5+fjyVNJkXwGebT68L+JftlwFsPyf1vURMt9V83aFudY627lM0oQwTdTTf1KllNsK599s9LhD7D6S27ofy41tIzRnLAj+2vVPB2M+RJg6DPPI3P260/ZfsNlrbdYe6SXqZdOwi9SxrrMAlYCnbi9dRrqiBDxN1NN9Iepo2H2Sq+UDXWSNawvNPmXut7TnAHKUpdkuqpb8x1NsLpG62q5gqYKFFAg+LYgqp+9pvSAMrii1c0cbauV+ymu6TH5ee8Gn55ge2P9H0cAwFtXYXrVLuJrliaxdRSXuSpoyo45pLbSQdbntirWWIJpSwKCS9hrSQ8wGkAU1nk5J50dVp8gCKPhWeC+VM0sRGrReOPwLsYPvAgrFbz3rm7qJ898krgUNt39+yfV1gogsv6DDUDIUeV5HAQ1dIWow0oOEU4GuuaYL7KuRJq84jzQPSPAJ1SdLiDo/18as9rb9xBWpaYGOkKH2tpRPRhBIWiaS3kbqUbQdcSxq4dE0FcfvrCXKM7V+Vip37er9N0o5AY96Ti4b7hWP6nxaitgmdalRslHOnogYeBk3S/cCTpPUBLwdeat5v+6YFf6trsWvrCdJHeVYpOeJ2KJD0A2AO8AU3JQ5JJwCr2j68tsIVlue+71NdZ5xRAw+L4n5Se+y7SHMjN/dJNmnF8lLq7AnSzm+Zvy/6cPRp0mLO9+YvTIBNgBuBYvOQDxGvqrsA7UQNPPSkAaYQ+IvtdSouT+3toVVRWqG90XQ0w3Zts/ENBZKWtf3swM/svhHbrzOUI2kXSZcWDjNVaZmr1tgfAW5o8/zSfjjwU4YH23+1fYHtC0gTeI0IklaTNF7SEvnxSpK+BtxTW5miBh4GK1/E+wHwelKvjK+R1oQU8NWS85EMtZ4gklYo3XVyKBoKXemqIOko0nS995LeYycD3yG9379he2Yt5YoEHgZL0s3A0cB1wG6kN/MXbZ9cYRmae4JUNQPkNqS24FeADwJfIQ1pXxx4j+3rSpdhqBgpTUeS7iTNv/53SWuQEvn2tq+vtVyRwMNgtZlWs7K2Z0kr9Le/ZG1Y0g3AYcBywAWkGv+1kjYDTrG9TanYQ42kxWy/Iulbtj9Td3lKafNe73Me/ipFL5SwKF6rtMBsg5ofF57SdTrz5l1Zg7RWoYDXAg9Sdjj94rZvB5A02/a1kLpNSiq9hN6QYvuVfPc9wLBN4MBYSd9rerxS82PbR9RQpkjgYZFcxfyDGZofm4JTutpeC+b2TT7f9m/z492AnUvFzZov/n+uZd8ShWMPVbXNiV2R/2x5PCTmfYkmlNDTJE23vXnLtmm2xxeMuRdwme3nWravA/y789Jbw00/zVYCbrU9tsryhEjgYRFI6q8LmW2fUUEZLiGtSvMzUq3/YNLFpXeVjj3SSLqPec1WrWx77YqLVJmm2S7bsr1XVWVpFgk8DJqkdusgitSMsprt4k10uVZ4HLB93nQ1abWYkhcxh+SHOZQjaTbwEHAWMJWWL7GSs1/2W65I4KEb8tJeBwGfBe4k9QO/rd5SlTFUP8ylSVoTeNL2P/LjdwD7kKZU+L7tf9VXurIkjQJ2IU3cthFp6cCzbM+otVyRwMOikDQaOJQ0T8ZU4ETbf6ow/hjgGFJf8KUa20vOTT1UP8yl5Rkf97X9iKRNgMuAE0n/gxdtD/f5UIC5CzsfCHwT+LLtdmeilYheKGHQJH0cOJK0Ms+uNa0WcyZpEYk9gP8AJgCzSwZ0WlT4YuDipg/zlZJq/TBXYGnbj+T7BwOn2/52ngv+lvqKVY38Wu9Oer3HAd+jYE+rjsoUNfAwWHmR21mkhLnAG8kFF9htKsN025tLuq0RT9JVtvtdsacLcVs/zOeTEtrfSsatU/OCDpJuAj5n+5L8+LYqXu+6SJpEWjT8d6QVp+6ouUhA1MDDolkfWJnUHtxsTeCRBZ9exIv5dqak3XPcot3ZWj7MJwyVD3MFLpd0DjCTtC7o5QCSVgWGbft3dghpVfo3AEekSz5ANQt49ylq4GHQJF0IfL71YqWk8cBxtouvWCJpD1I3wtVJy7m9mpRU++0psogxXyF9mGH+M49aP8yl5QvV7wVWBc5pnG1I2hRYqVEbD9WJBB4Grb/5IPpbPzH0Lkn7AOsCt0fCrl80oYRFsVQ/+4rOCZL7oPdZ+yg5N4WkHRuzHkpay/Z9Tfv2KzwHTG0knQpsCPwR+C9JW9r+r5qLNaLFgg5hUdzYx6IKh1F+rohpOcZ0YK+m+42fkr7VdP/XLfu+UDh2nbYDdrT9OWAHUh/wUKOogYdFcRRwrqSDmJc0x5MmdNq3ZGDbkxr3JR3V/LgC6uN+u8fDyb9yF0psP6emK3mhHpHAw6DlVW/elkfkNdrCL6piUYXWotQYrzX2cL6otIGkxgVrAevkx42Lt8O2G+FQFQk8LDLbVwBX1F2OCq2d50NR033y45LzkNftjXUXIMwveqGEniTpaebVdpcBGlO7Fu/KJ6nfQULDdS6UMPREAg8hdKTlS3O+XQzj/u9DWSTwEBaSpL2Bsba/nx9PBcbk3cfY/lVthQsjSnQjDGHhHUOa+6RhSWALUte6j9ZRoDAyxUXMEBbeErab53+51vYcYI6kZesqVBh5ogYewsJbvvmB7U80PRxDCBWJBB7CwpvaxwjUjwA31FCeMELFRcwQFpKklYDzgBeAm/LmzUlt4fvkAU4hFBcJPIRBkrQjaSk3gBk1jEANI1wk8BAWkqQV+ttv++9VlSWMbJHAQ1hIku4jDWgRsAbwRL7/WuBB28N5OH0YQuIiZggLyfZattcGLgH2tL2i7deRFlYelnOBh6EpauAhDFJjQeWWbdNsj6+rTGFkiYE8IQze45K+APyM1KRyMDCn3iKFkSSaUEIYvANJA3fOJXUrXClvC6ES0YQSQgg9KppQQhgkSWNIE1u9iaYFnm3vWFuhwogSTSghDN6ZwN2kVXhOAO4HbqyzQGFkiSaUEAap0QtF0m2N9SAlXWW73xV7QuiWaEIJYfBezLczJe0OPAKMrbE8YYSJBB7C4H1F0muATwOnAK8Gjq63SGEkiSaUEELoUVEDD2EhSTqF9ov7AmD7iAqLE0awSOAhLLxpTfdPAI6rqyBhZIsmlBAWgaSbbW9adznCyBT9wENYNFEDCrWJBB5CCD0qmlBCWEiSnmZezXsZ4LnGLsC2X11LwcKIEwk8hBB6VDShhBBCj4oEHkIIPSoSeAgh9KhI4CGE0KMigYcQQo/6/+XdhUB+8hdeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# visualize the δRsfc\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(np.arange(len(sfc_lw)), np.array(list(sfc_lw.values())))\n",
    "ax.set_xticks(np.arange(len(sfc_lw)))\n",
    "ax.set_xticklabels(sfc_lw.keys(), rotation=90)\n",
    "ax.set_title(\"LW_sfc, global & time average\")\n",
    "# fig.savefig(\"/Users/brianpm/Dropbox/Projects/extremes/cmip_lwoff_OLRcorrect_bars.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('rsdt', 340.35496689055265), ('rsut', 91.19573888473442), ('rlut', 235.43461851563595), ('rsds', 200.56782329133765), ('rsus', 26.108805379914813), ('rlds', 334.21612020749205), ('rlus', 393.8249469077779), ('hfss', 21.097049177031543), ('hfls', 79.74416502699901), ('pr', 3.1894594603906825e-05)]\n"
     ]
    }
   ],
   "source": [
    "print([(i, glb[i].values.item()) for i in glb])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "79.76838110437097"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "glb['pr'].values.item()*2.501e6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_values([241.63964891654496, 249.15922800581822, 238.68703648594055, 241.36339705599215, 240.16539001464844, 249.72018432617188, 241.13717530852193, 246.22788484223662])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "toa_sw.values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/brianpm/miniconda3/envs/py38/lib/python3.8/site-packages/xarray/conventions.py:492: SerializationWarning: variable 'pr' has multiple fill values {1e+20, 1e+20}, decoding all values to NaN.\n",
      "  new_vars[k] = decode_cf_variable(\n"
     ]
    }
   ],
   "source": [
    "cc = xr.open_dataset(dloc / 'pr_Amon_CESM2_amip_r1i1p1f1_gn_195001-201412.nc')"
   ]
  },
  {
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}