{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radiative-Convective Equilibrium with CAM3 scheme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import climlab\n",
    "from climlab import constants as const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Here is how to set a simple RCE in `climlab`\n",
    "\n",
    "By initializing each component with the same state object, the components are already effectively coupled. They all act to modify the same state object.\n",
    "\n",
    "No extra coupling code is necessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial state (temperatures)\n",
    "state = climlab.column_state(num_lev=20, num_lat=1, water_depth=5.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Create individual physical process models:\n",
    "#  fixed relative humidity\n",
    "h2o = climlab.radiation.ManabeWaterVapor(name='H2O', state=state)\n",
    "#  Hard convective adjustment\n",
    "convadj = climlab.convection.ConvectiveAdjustment(name='Convective Adjustment',\n",
    "                                                  state=state, \n",
    "                                                  adj_lapse_rate=6.5)\n",
    "# CAM3 radiation with default parameters and interactive water vapor\n",
    "rad = climlab.radiation.CAM3(name='Radiation',\n",
    "                             state=state, \n",
    "                             specific_humidity=h2o.q)\n",
    "\n",
    "rce = climlab.couple([rad,convadj,h2o], name='RCM')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.process.time_dependent_process.TimeDependentProcess'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (1,) \n",
      "  Tatm: (20,) \n",
      "The subprocess tree: \n",
      "RCM: <class 'climlab.process.time_dependent_process.TimeDependentProcess'>\n",
      "   Radiation: <class 'climlab.radiation.cam3.cam3.CAM3'>\n",
      "   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(rce)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AttrDict({'Ts': Field([288.]), 'Tatm': Field([200.        , 204.10526316, 208.21052632, 212.31578947,\n",
       "       216.42105263, 220.52631579, 224.63157895, 228.73684211,\n",
       "       232.84210526, 236.94736842, 241.05263158, 245.15789474,\n",
       "       249.26315789, 253.36842105, 257.47368421, 261.57894737,\n",
       "       265.68421053, 269.78947368, 273.89473684, 278.        ])})"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Current state\n",
    "rce.state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1826 steps, 1826.2110000000002 days, or 5 years.\n",
      "Total elapsed time is 4.999422301147019 years.\n"
     ]
    }
   ],
   "source": [
    "#  Integrate the model forward\n",
    "rce.integrate_years(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AttrDict({'Ts': Field([276.77058287]), 'Tatm': Field([233.26150428, 215.95044945, 210.60589023, 211.45113437,\n",
       "       212.04321205, 216.46759988, 223.46187556, 229.6327028 ,\n",
       "       235.16956012, 240.2017945 , 244.8218937 , 249.09842749,\n",
       "       253.08373297, 256.81872305, 260.33602553, 263.66210602,\n",
       "       266.81874721, 269.82410613, 272.69348696, 275.43991673])})"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Current state\n",
    "rce.state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Field([1.87420483e-05, 9.64556905e-06, 5.58553481e-06, 6.57792987e-06,\n",
       "       7.30688198e-06, 1.30016777e-05, 3.02115363e-05, 6.04582768e-05,\n",
       "       1.08656640e-04, 1.80140927e-04, 2.80533313e-04, 4.15636693e-04,\n",
       "       5.91349145e-04, 8.13596531e-04, 1.08827995e-03, 1.42123518e-03,\n",
       "       1.81820176e-03, 2.28479967e-03, 2.82651219e-03, 3.44867359e-03])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Current specific humidity\n",
    "rce.q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'specific_humidity': Field([1.87420483e-05, 9.64556905e-06, 5.58553481e-06, 6.57792987e-06,\n",
       "        7.30688198e-06, 1.30016777e-05, 3.02115363e-05, 6.04582768e-05,\n",
       "        1.08656640e-04, 1.80140927e-04, 2.80533313e-04, 4.15636693e-04,\n",
       "        5.91349145e-04, 8.13596531e-04, 1.08827995e-03, 1.42123518e-03,\n",
       "        1.81820176e-03, 2.28479967e-03, 2.82651219e-03, 3.44867359e-03]),\n",
       " 'absorber_vmr': {'CO2': 0.000348,\n",
       "  'CH4': 1.65e-06,\n",
       "  'N2O': 3.06e-07,\n",
       "  'O2': 0.21,\n",
       "  'CFC11': 0.0,\n",
       "  'CFC12': 0.0,\n",
       "  'CFC22': 0.0,\n",
       "  'CCL4': 0.0,\n",
       "  'O3': array([5.38853507e-06, 9.86362297e-07, 3.46334801e-07, 1.90806332e-07,\n",
       "         1.19700066e-07, 7.69083554e-08, 5.97316411e-08, 5.27011190e-08,\n",
       "         4.80406196e-08, 4.44967931e-08, 4.18202246e-08, 3.99595858e-08,\n",
       "         3.83838549e-08, 3.66179869e-08, 3.42885526e-08, 3.18505117e-08,\n",
       "         2.93003951e-08, 2.69906527e-08, 2.49122466e-08, 2.28798533e-08])},\n",
       " 'cldfrac': 0.0,\n",
       " 'clwp': 0.0,\n",
       " 'ciwp': 0.0,\n",
       " 'r_liq': 0.0,\n",
       " 'r_ice': 0.0,\n",
       " 'emissivity': 1.0,\n",
       " 'S0': 1365.2,\n",
       " 'insolation': 341.3,\n",
       " 'coszen': 0.25,\n",
       " 'eccentricity_factor': 1.0,\n",
       " 'aldif': 0.3,\n",
       " 'aldir': 0.3,\n",
       " 'asdif': 0.3,\n",
       " 'asdir': 0.3}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Here is the dictionary of input fields for the CAM3 radiation module\n",
    "rce.subprocess.Radiation.input"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Latitudinally, seasonally varying RCE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial state (temperatures)\n",
    "state2 = climlab.column_state(num_lev=20, num_lat=30, water_depth=10.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.process.time_dependent_process.TimeDependentProcess'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (30, 1) \n",
      "  Tatm: (30, 20) \n",
      "The subprocess tree: \n",
      "Seasonal RCE: <class 'climlab.process.time_dependent_process.TimeDependentProcess'>\n",
      "   Insolation: <class 'climlab.radiation.insolation.DailyInsolation'>\n",
      "   Radiation: <class 'climlab.radiation.cam3.cam3.CAM3'>\n",
      "   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#  Create a parent process\n",
    "rcelat = climlab.TimeDependentProcess(state=state2)\n",
    "## Create individual physical process models:\n",
    "#  seasonal insolation\n",
    "insol = climlab.radiation.DailyInsolation(name='Insolation', \n",
    "                                          domains=state2['Ts'].domain)\n",
    "#  fixed relative humidity\n",
    "h2o = climlab.radiation.ManabeWaterVapor(name='H2O',\n",
    "                                         state=state2)\n",
    "#  Hard convective adjustment\n",
    "convadj = climlab.convection.ConvectiveAdjustment(name='Convective Adjustment',\n",
    "                                                  state=state2, \n",
    "                                                  adj_lapse_rate=6.5)\n",
    "# CAM3 radiation with interactive insolation and interactive water vapor\n",
    "rad = climlab.radiation.CAM3(name='Radiation',\n",
    "                             state=state2, \n",
    "                             specific_humidity=h2o.q,\n",
    "                             S0 = insol.S0,\n",
    "                             insolation=insol.insolation,\n",
    "                             coszen=insol.coszen)\n",
    "\n",
    "rcelat = climlab.couple([insol,rad,convadj,h2o], name='Seasonal RCE')\n",
    "\n",
    "print(rcelat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1826 steps, 1826.2110000000002 days, or 5 years.\n",
      "Total elapsed time is 4.999422301147019 years.\n"
     ]
    }
   ],
   "source": [
    "rcelat.integrate_years(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 365 steps, 365.2422 days, or 1 years.\n",
      "Total elapsed time is 5.9987591795252575 years.\n"
     ]
    }
   ],
   "source": [
    "rcelat.integrate_years(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_temp_section(model, timeave=True):\n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    if timeave:\n",
    "        field = model.timeave['Tatm'].transpose()\n",
    "    else:\n",
    "        field = model.Tatm.transpose()\n",
    "    cax = ax.contourf(model.lat, model.lev, field)\n",
    "    ax.invert_yaxis()\n",
    "    ax.set_xlim(-90,90)\n",
    "    ax.set_xticks([-90, -60, -30, 0, 30, 60, 90])\n",
    "    fig.colorbar(cax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAf9klEQVR4nO3df7AdZZ3n8feHAAEJapgENoSwiUzICFkHnJhipUYZUKGcKcPMllamVio1skONlVGwdFfQrbGs2WyhY7Hr1i4zdUtwUg7KZoSRrLuigZWxKIEYECQhZglmBi5ECAISZyQxyXf/6L7QXO75eU93P93n86q6dc/p2+f0c053f863n/P0bUUEZmaWrqPqboCZmXXnoDYzS5yD2swscQ5qM7PEOajNzBLnoDYzS5yD2sxsFiQdJ2mrpIck7ZD02Xz6SZK2SHo0/z2/8JhrJO2WtEvSxT2X4XHUZmbDkyTghIj4haRjgLuBK4E/AJ6LiGslXQ3Mj4hPSjoL+BqwGjgVuAM4MyIOd1qGK2ozs1mIzC/yu8fkPwGsATbm0zcCl+a31wA3R8SBiNgD7CYL7Y6OHnWjh7FgwYJYunTpy/d3PfrTUpd3+DiN9vnmjvTpGmvO3I4FQTIOH5hTdxOSMOdACc/5UjlH5yuW/4tX3b///vufjYiFs3nO8y84Ll547khf8z7y8K92AC8VJk1ExERxHklzgPuBXwf+R0TcJ+mUiNgLEBF7JZ2cz74YuLfw8Ml8WkdJBPXSpUvZtm3by/cvuORzpS7v+RXHjvw595/R30o3q9uJj43+QHr+roMjf84pd93+yVfdl/SPs33OF547wle/eUpf857zLydfiohV3ebJuy3OkfRG4O8krewy+0yVYtdPuSSCumpTG9UoA3tq43dgW6rKCGgoN6SbJiJekHQXcAnwtKRFeTW9CHgmn20SWFJ42GnAU92ed6z7qMvYwE587KjSdgizYZS5TTqkQdLCvJJG0vHAu4AfA5uBdfls64Db8tubgbWS5kpaBiwHtnZbxlhW1EXzdx0spSukuGO4yraqVVEsOKRftgjYmPdTHwVsiohvSroH2CTpcuBx4P0AEbFD0ibgEeAQsL7biA9wUAPlhfUUh7ZVpaqjOYf0KyLiR8C5M0z/GXBRh8dsADb0uwwHda7ssJ7i0LZRqrqbzQFdDwd1QVVhPWX6Tubgtl7q/P7DIV0fB/U0VYd1Uaed0AE+flL7QtohXS8H9QzqDOuZOMDbK7VAnolDun4O6g5SC+uZDLqTO9jL14TgHYRDOg1JBvXx258E4Jcru55VWbomhPUg2hYiVi6HdDqS3nOnArtO83cd9AZrYyeVbT6FDEhB0kFtZtVLJaTtFckHdSqfqN54bRyktJ2nsu+nIPmgTklKG7HZqHn7TlcjgtqfrGbjxfv8qzUiqFPiqsPayNt12hzUQ/BGbW2S2vbsavq1GhPUXnlmo5daSNvMGhPUqfEGbjZ6Lshm1qigTm0lOqytybz9NkejgtrMRiPFkE6tEEtJ44I6tZWZ4gZvZu3SuKBOkcPamiTF7TW1Aiw1jQxqr1Sz9vD+3Fsjgxq8cs3awPtxfxob1GY2uBS7Pay3Rgd1Sp/G3gHMBpPS/jsbkpZI+q6knZJ2SLoyn36OpHslPShpm6TVhcdcI2m3pF2SLu61jCSv8GJmo5dSMdGWkM4dAj4eEQ9IOhG4X9IW4PPAZyPiW5Lem9+/QNJZwFrgbOBU4A5JZ0bE4U4LaHRFDa1b4WbWMBGxNyIeyG/vB3YCi4EAXp/P9gbgqfz2GuDmiDgQEXuA3cBqunBFPUJtu8aiWRlSKK5+dngeX3n+7X3OvWmBpG2FCRMRMTHTnJKWAucC9wFXAd+W9AWyonhqgYuBewsPm8ynddSKoD5++5O1XwjXLGUpdXs00LMRsarXTJLmAbcAV0XEi5L+E/CxiLhF0geAG4B3AZrh4dHtuRvf9WFmzZFCNV0GSceQhfRNEXFrPnkdMHX7b3mle2MSWFJ4+Gm80i0yo9YEdSobgCsXs5mlso+OmiSRVcs7I+K6wp+eAt6Z374QeDS/vRlYK2mupGXAcmBrt2W0ouvDzNL3y5WL2xrW5wOXAQ9LejCf9ingj4EvSjoaeAm4AiAidkjaBDxCNmJkfbcRH9CioHYftZnVISLuZuZ+Z4Df6vCYDcCGfpfRs+ujy2DukyRtkfRo/nt+4TEDDeY2M7PO+umjnhrM/WbgPGB9PmD7auDOiFgO3JnfZ9pg7kuA6yXNKaPxKfLwPLPOfOQ7nJ5dHxGxF9ib394vaWow9xrggny2jcBdwCcpDOYG9kiaGsx9z6gbP8UrP037zzhS27JPfKw135ObDdZHPW0w9yl5iBMReyWdnM/W12BuSVeQd66ffvrpAzfc0lJnKM+kU3sc4NZEfW+10wdzd5t1hmmvGcwdERMRsSoiVi1cuLDfZryGq+n67D/jyMs/TdHENreN99nB9VVRdxjM/bSkRXk1vQh4Jp8+8GBua442BVzxtbjStpT1M+qj02DuzWRn3pD/vq0wfaDB3MPyJ3N12l6Ftv31pfYlt/fdwfRTUXcazH0tsEnS5cDjwPthuMHclqY2B1cnrrItRf2M+ug2mPuiDo8ZaDC3pWccQ3q6/WcccVhbEhq7FfrQqRxt7wIYlN+P8ngf7l9rTiG32XEYdTf1/rjCtjo0Mqj9STw6DujBOLCtDt7axpQP6WfH751VqXFB7Wp69hwyo+EPu9nz/tyfRgW1V+rsOVhGz++pla1RQW3Dc/VXLr+/w3MB1ltjgtorc3gOkOr4vbYyNCKoHdLDc3BUz+/54LyPd9fI4XnWm8OiXh7GZ6OU/FbkT9rBOaTT4XXRP+/rnSUd1E1bcSn8hzIHQ3q8TvrXtH2+Kkl2fXhlDceBkC7/g6f+Te3/x29/suaWpMNbTks4pNPndTSYX65c3IiiTdISSd+VtFPSDklXFv72EUm78umfL0y/RtLu/G8X91pGkhV1E9XZ7eEAaI46K+vnVxzL/F0Ha1l2yx0CPh4RD0g6Ebhf0hbgFLKLfb8lIg5MXVdW0lnAWuBs4FTgDklndvu//a6oG84h3TxeZ+0SEXsj4oH89n5gJ9kFvT8MXBsRB/K/TV2ucA1wc0QciIg9wG5gdbdluKI2q0FdlbWr6syLB4/jO4+v6Hf2BZK2Fe5PRMTETDNKWgqcC9wH/AXw25I2AC8Bn4iIH5CF+L2Fh03m0zpyUI9AXd0ersyazV8wNsazEbGq10yS5pFdBPyqiHhR0tHAfOA84G1kly58EzNfMSu6PbeDuqHaFtLzlv28499+secNFbbEbHCSjiEL6Zsi4tZ88iRwa0QEsFXSEWBBPn1J4eGnAU91e34HdQM1OaS7BfJsHtPUMHdV3XySBNwA7IyI6wp/+gZwIXCXpDOBY4Fngc3AVyVdR/Zl4nJga7dlOKhnqepuj6aF9DDBPIrlNCm4qw5r91OP3PnAZcDDkh7Mp30KuBG4UdJ24CCwLq+ud0jaBDxCNmJkfbcRH+CgbpSmhHRV4dxvG5oQ2q6smysi7mbmfmeAD3Z4zAZgQ7/LcFA3ROohnUI4d9KU0HZYWyfeKmYhhf/tUbd5y36edEhP17T2lsXbbrO4om6AFKvppofdVPtTq7BdVdtMvEUMqaqKJLWQbltFmuLrqWqdu6puDge19S21QBulFAPbbIqDegjjVk2PU4il8jpdVVuR+6gTlUJIpxJaVUul/9r91TbFW8GAxqUCGdeQLhqX92Bctukmc1AnqO5qelwCqh91vxd1bwuWBgf1AKqoPOreMesOphTV3UdfxTbhqjptDmoD6g+jJvD7Y3VxUPepzdW0A6h/db1XdR9pWb0c1GPOIT24toa1uz/S5aBORB0Vk0N6eH7vrEoO6j60sdJw0MxeHe+hq+rx1HdQS5oj6YeSvpnfP0nSFkmP5r/nF+a9RtJuSbskXVxGw9uk6mraIT06fi+tCoNU1FeSXQZ9ytXAnRGxHLgzv4+ks4C1wNnAJcD1kuaMprnVa1uF4WAZvarfU1fV46evoJZ0GvC7wJcKk9cAG/PbG4FLC9NvjogDEbEH2A2sHklrW6jKatohXR6/t1amfivq/wr8B6CYKqdExF6A/PfJ+fTFwBOF+Sbzaa8i6QpJ2yRt27dv36DtrkTZlYWHXLVLlWHtqnq89AxqSb8HPBMR9/f5nDNdOyxeMyFiIiJWRcSqhQsX9vnU1WnbhuqKrxp+n60M/VTU5wPvk/QPwM3AhZL+Bnha0iKA/Pcz+fyTwJLC408DnhpZi1vCXR42W66qx0fPoI6IayLitIhYSvYl4f+NiA8Cm4F1+WzrgNvy25uBtZLmSloGLAe2jrzlJWrTBuqQrp7fcxu12YyjvhZ4t6RHgXfn94mIHcAm4BHgdmB9RByebUPbpKpq2oFRn6ree1fV9ZO0RNJ3Je2UtEPSldP+/glJIWlBYdpAQ5gHunBARNwF3JXf/hlwUYf5NgAbBnnuVLRlw3RI12/esp/XfvEBq8Qh4OMR8YCkE4H7JW2JiEckLSErZB+fmnnaEOZTgTskndmtoPWZiRXzSI/xUsUHpqvqekXE3oh4IL+9n+x8k6mRbv+FbMRccUDFwEOYHdQFbdkgXU2b1UPSUuBc4D5J7wOejIiHps3W1xDmIl8zsWUc0umpoguk7OsrPr/iWObvOlja81ft8IE5g6yTBZK2Fe5PRMTE9JkkzQNuAa4i6w75NPCeGZ6vryHMRQ7qXJv/37TVz/3VjfZsRKzqNoOkY8hC+qaIuFXSvwKWAQ9JgmyY8gOSVjPEEGZ3fbSIq2krU1u6BkdNWRLfAOyMiOsAIuLhiDg5IpbmQ5sngbdGxE8ZYgizg5p2VNMO6fSVvY58xFab84HLyE4GfDD/eW+nmYcZwuyuDzPrW9v6qkchIu5m5n7n4jxLp90faAjz2FfUrqatSq6qbRhjH9RmVWv6B6v7qqs31l0frqbT8J7Tdw38mO88vqKElrRD2UP1rHpjG9RtqAqaGtLDBHOv52hacDd9uJ77qqs1tkFt1RlFMA+yjKaFdlM5rKszlkFdVTVdZrdHE6rpKgK623JTD+wyq2p3f7TL2K3JNnR5NEFdIT29DSm0o5smfOB24/2pGmMV1FVuVONaTacYjqm1pypVDdVzWJdvLLs+mizVkE49DFPuDmn6F4tWvrEJan/qlyP1gJ4u5cBuMn+xWK6x6PqoOqTH5eywpoV0UWptL+tIqcpt0cVQecYiqNsipW6P1IJuGG14DalxWJej9UHtDWf02hRwKb2WNlTVVo5WB3UdIV3WTpFKNZ1SsI1KiiNVmszF0ei1Nqi9sYzWOIRZCq8vlQ/k2fL+N1qtDOq6NpK2VtMpBFhV2vpa3f3RbK0Lan+Sj1Zbg6ubcXzNZfC+ODqtCuo2bhh1VtPjHFh1vva6j6BG6fkVx7Zyv6xaa4K67o2hbYeW4xzSU9r2HtS5jda9fzZdK4K6rRtBmyqrpqorrNu47tu6n1ah0UHtw6pytK2StHR4fx1OY4M6pRXepm4Ph/RrtamqTmFbTWnfbYrGBfW4VNF1HPo6pDvzezNabdqHJS2R9F1JOyXtkHRlPv0vJP1Y0o8k/Z2kNxYec42k3ZJ2Sbq41zIaEdRT4dymlZsaB1Fvfo9Gq0X78yHg4xHxZuA8YL2ks4AtwMqIeAvw/4BrAPK/rQXOBi4Brpc0p9sCkg3qpoRzGYeSbfwiqS2qDuu2dn9Macp+3k1E7I2IB/Lb+4GdwOKI+E5EHMpnuxc4Lb+9Brg5Ig5ExB5gN7C62zKS/H/UTV5pTeRK0VJQ3O+b+r+tJS0FzgXum/anDwH/M7+9mCy4p0zm0zpKMqjHWdXVtEN6cO85fVelFx4YxyvATIV2WYE95wCDXPx3gaRthfsTETExfSZJ84BbgKsi4sXC9E+TdY/cNDVphmVEtwY4qGchpUPIYTikh1d1WI9aU65SnsjR9bMRsarbDJKOIQvpmyLi1sL0dcDvARdFxFQYTwJLCg8/DXiq2/Onv6bMEuUPOgOQJOAGYGdEXFeYfgnwSeB9EfHPhYdsBtZKmitpGbAc2NptGa6oE1Jlt4dDplnGsfujQc4HLgMelvRgPu1TwH8D5gJbsizn3oj4k4jYIWkT8AhZl8j6iDjcbQF9BXU+/u9LwEqyvpQPAbvIOseXAv8AfCAins/nvwa4HDgMfDQivt3Xy22QJnd7OKRHp8ldIE3p/khdRNzNzP3O/6fLYzYAG/pdRr9r6YvA7RHxG8Bvkg0/uRq4MyKWA3fm94caI2hmZp31DGpJrwfeQdYHQ0QcjIgXyMYCbsxn2whcmt8eeIygVdft4Wp69Kp6Tz2+fnz1U1G/CdgHfFnSDyV9SdIJwCkRsReyAd/Ayfn8i4EnCo+fcYygpCskbZO0bd++fbN6EVVrcreHWZG35WboJ6iPBt4K/GVEnAv8E3k3Rwd9jRGMiImIWBURqxYuXNhXY212XE2Xx1W1lamfoJ4EJiNi6kybr5MF99OSFgHkv58pzD/QGMFx552vHfxBaGXpGdQR8VPgCUlTX21fRDasZDOwLp+2Drgtvz3wGMEmaeqhokPEOmnqNj1O+h1H/RHgJknHAj8B/ogs5DdJuhx4HHg/wDBjBK1cDunqVDFcz2Oqx09fQR0RDwIznUJ5UYf5BxojOM7c7WFmvXi0+wCaeIjoarp6TXzPm7htjxMHtVkD+UhsvDioa1T2ztbEyq4t/N7bKDmo++RDQ2s7b+PpclC3lCu6+nkd2Kg4qGviPkabLW9D48NB3UKu5NLhdWGj4KDug/vubFx4W0+Tg7oGZR6yuoJLT5nrxN0f48FBbWaWOAd1Dz4UtFFo0pGOt/n0OKgr5m4PGzV3f7Sfg9rMLHEO6i6adAjoajp9TVpHTdr2x0G//4/arDSXzf/+wI/5yvNvL6ElzeX/Ud1urqg7KKOicF/iKy6b//2Xf2bz+KZpUlVt/ZG0RNJ3Je2UtEPSlfn0kyRtkfRo/nt+4THXSNotaZeki3stw0HdAk3a+UcdsE0N7CZw90ffDgEfj4g3A+cB6yWdRXYR8DsjYjlwZ36f/G9rgbOBS4DrJc3ptgB3fVjpqgjS4jLcLWJVioi9wN789n5JO4HFwBrggny2jcBdwCfz6TdHxAFgj6TdwGrgnk7LcFDPoEndHqlX03VUu5fN/37SYV3WdRXdT92/OS8F83cd7Hf2BZK2Fe5PRMTETDNKWgqcC9wHnJKHOBGxV9LJ+WyLgXsLD5vMp3XkoLbS1NklMbXslAO7CfafcYQTHxv7HtJnI2Kma8a+iqR5wC3AVRHxoqSOs84wLbo999ivARu9lPqNU2nHdKkfCdlgJB1DFtI3RcSt+eSnJS3K/74IeCafPgksKTz8NOCpbs/voJ6mSV+gpLizpxiMKbapLB5ZVD1lpfMNwM6IuK7wp83Auvz2OuC2wvS1kuZKWgYsB7Z2W4aDugLjsvOkHIgpty11TSpeanI+cBlwoaQH85/3AtcC75b0KPDu/D4RsQPYBDwC3A6sj4jD3RbgPuqGSq2abkIQpvYlY1lfKlq1IuJuZu53Briow2M2ABv6XYYr6gJXDsNpQkhPSan/3KxfDmqblaaGXirtLuPIqIyuNhcx9XJQl6yMnSaVbo9Uws6s7RzUOVcMg2lDSLfhNdh4cFDbwNoUcCm8llSOkHpxMVMfB3WJ2tzt0SYphPWojcuQ0HHhoMaVwiDaGGpmqXNQW9/aHNJ1v7amHCm5qKmHg7pB6tyZ6w6yKozDa7RmclCXpE19hOMUYG16rW3aBsfd2Ae1D+UsFe7+sE7GPqiboq6duE0VZr/G8TVb2sY6qF0ZdDfOgTXOr93S01dQS/pYfnXd7ZK+Jum4UV5ht23cN2ip8LbYDj2DWtJi4KPAqohYCcwhu4LuyK6wW4cmVdN1dHu4oqznPXA/tc2k366Po4HjJR0NvI7ssjFryK6sS/770vz2y1fYjYg9wNQVds0axx9YloKeQR0RTwJfAB4nuyT6zyPiO0y7wi5QvMLuE4WnmPEKu5KukLRN0rZ9+/bN7lUMyNVAdw4n64f3o+r00/Uxn6xKXgacCpwg6YPdHjLDtNdcYTciJiJiVUSsWrhwYb/tNatc1R9co+7+cD918/XT9fEuYE9E7IuIXwG3Am9nhFfYrVLZVcCod4qq+yxdTdsgXFVXo5+gfhw4T9Lr8qvtXgTsZIRX2DVLnT/ArE49L24bEfdJ+jrwAHAI+CEwAcwDNkm6nCzM35/Pv0PS1BV2D9HHFXar4k//7hxGNoz9ZxzhxMfG+pSM0vX17kbEZyLiNyJiZURclo/o+FlEXBQRy/PfzxXm3xARZ0TEioj4VnnNb7cquz0c0r1V+R41ZZieVWNsPgarqKb9pY2lquxtc5yPViXdKOkZSdsL086RdK+kB/PRbasLfxv4hMCxCWrrzNV0//xe2Qz+muzkvqLPA5+NiHOAP8vvD31CoIPaLFFN6/4Y16o6Ir4HPDd9MvD6/PYbeGXk21AnBPb8MrENmrgBVbWTukIc3GXzv89Xnn973c2wWTjql7/i+O1P9jv7AknbCvcnImKix2OuAr4t6QtkBfHUBrMYuLcw34wnBE43FkFdBfdPW+rmLfs5v9jzhlKX0dIRIM9GxKoBH/Nh4GMRcYukDwA3kJ2T0tcJgdO17h2dronVdFVcTQ/P7531sI7s5ECAv+WV7o2hTghsfVA3UdP6Jq08TdwWXBwBWfi+M799IfBofnuoEwJb3fXhDaYzV4Sz575qA5D0NeACsr7sSeAzwB8DX8z/4+hLwBUw/AmBrQ1qh7SZVSEi/rDDn36rw/wbgA2DLMNdHyPgLxLHV9OOTKraVl0ojVYrg7rJG0kVfZJNC5dx18R+ahut1gV1k0ParE28L45Oq4LaG0ZvrqZHz+9pZ94nR6NVQV0H909b01S9zTqsZ681Qd2GjcF9kc1VdlXtbWO8tSKo2xDSVfAhutXF++jstCKozSx9DuvhNT6ovfItFT5i6c3763AaHdR1r/RRfilTdh+kQ6T5RrmN+EvwZmlsUNcd0mY2HO+7g2tkUHtFD8bVdHX8XvfH+/BgGvVPmbxyzWwcNaKi3n/GkeRC2n181kmZVXWb+qlT26dTlnxQj8PKLPOLRB+KW8pSLMJSlGxQewWajQ/v690lGdReadZ0PpIZnIuzzpIMahsNh4U1kQP7tRzUQ2jSiS7WPm36QrEbh/UrHNQt5Wq6fl4Hs+ewzjiozcwS56A2M0ucg7qFfMidjrLWhb/bSIekGyU9I2n7tOkfkbRL0g5Jny9Mv0bS7vxvF/ezDAf1gPxForVNyl8oNsRfA5cUJ0j6HWAN8JaIOBv4Qj79LGAtcHb+mOslzem1AAe1mdksRMT3gOemTf4wcG1EHMjneSafvga4OSIORMQeYDewutcyFBEjbPJwJO0HmlZeLgCerbsRA2pam5vWXnCbq7AiIk6czRNIup3sdffjOOClwv2JiJiY9nxLgW9GxMr8/oPAbWRV80vAJyLiB5L+O3BvRPxNPt8NwLci4uvdGpDKf8/bFRGr6m7EICRtc5vL1bT2gttcBUnbZvscEXFJ77lm5WhgPnAe8DZgk6Q3AZqpOb2ezF0fZmajNwncGpmtwBGyCn4SWFKY7zTgqV5P5qA2Mxu9bwAXAkg6EziWrHtpM7BW0lxJy4DlwNZeT5ZK18dE71mS4zaXr2ntBbe5Ckm1V9LXgAuABZImgc8ANwI35kP2DgLrIvtCcIekTcAjwCFgfUQc7rmMFL5MNDOzztz1YWaWOAe1mVniag1qSb8p6R5JD0v6X5JeX/jbwKdZVmWUp4aWTdKfS/qRpAclfUfSqYW/JdfeKZIuydu1W9LVdbdnOknHSdoq6aF8O/hsPv0kSVskPZr/nl93W4skvVHS1yX9WNJOSf865TZLulLS9vw9viqflmx7SxMRtf0APwDemd/+EPDn+e2zgIeAucAy4DFgTp1tLbT5d4A7gLn5/ZNTbjPw+sLtjwJ/lXJ787bNydvzJrJvyx8Czqq7XdPaKGBefvsY4D6yMbOfB67Op18NfK7utk5r90bg3+W3jwXemGqbgZXAduB1ZAMf7iAbJZFke8v8qbvrYwXwvfz2FuDf5LeHOs2yIiM9NbRsEfFi4e4JvDK4Psn25lYDuyPiJxFxELiZrL3JiMwv8rvH5D9B1s6N+fSNwKXVt25m+RHrO4AbACLiYES8QLptfjPZWXz/HBGHgL8Hfp9021uauoN6O/C+/Pb7eWUg+GLgicJ8k/m0FJwJ/Lak+yT9vaS35dOTbbOkDZKeAP4t8Gf55GTbS9pte5mkOfmpws8AWyLiPuCUiNgLkP8+ucYmTvcmYB/wZUk/lPQlSSeQbpu3A++Q9GuSXge8lywjUm1vaUoPakl35H1M03/WkHV3rJd0P3Ai2XhDGPI0y4raXDw19N+TnRqqOtvco71ExKcjYglwE/CnUw+rq719SLltL4uIwxFxDtnZZaslray5Sb0cDbwV+MuIOBf4J7KugyRFxE7gc2RH27eTdYEdqrVRNSn9hJeIeFePWd4DL5+987v5tKFOsxyVbm2W9GHyU0OBrZJmdWroKPTxHk/5KvC/yQbk1/oe95By214jIl6QdBfZP+B5WtKiiNgraRFZtZ2KSWAyr/wBvk4W1Mm2OSJuIO+qkfSfyV5Dsu0tS92jPk7Ofx8F/Efgr/I/DXWaZUW+wQhPDS2bpOWFu+8DfpzfTrK9uR8AyyUtk3Qs2f/v3Vxzm15F0kJJb8xvHw+8i+y93Qysy2dbR/Yf1JIQET8FnpC0Ip90EdkZcsm2uZARpwN/AHyNhNtblrpPIf9DSevz27cCXwaIiKFOs6zISE8NrcC1+Y55BPhH4E8g7fc4Ig5J+lPg22QjQG6MiB01N2u6RcBGZf/0/ShgU0R8U9I9ZN1hlwOPk333kpKPADflH4A/Af6IvP2JtvkWSb8G/IpsG31e0rWk295S+BRyM7PE1T3qw8zMenBQm5klzkFtZpY4B7WZWeIc1GZmiXNQm5klzkFtZpa4/w/fEy5lCd/iiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_temp_section(rcelat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Same thing, but also including meridional temperature diffusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Create and exact clone of the previous model\n",
    "diffmodel = climlab.process_like(rcelat)\n",
    "diffmodel.name = 'Seasonal RCE with heat transport'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3964424.9422310763\n"
     ]
    }
   ],
   "source": [
    "# thermal diffusivity in W/m**2/degC\n",
    "D = 0.05\n",
    "# meridional diffusivity in m**2/s\n",
    "K = D / diffmodel.Tatm.domain.heat_capacity[0] * const.a**2\n",
    "print(K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = climlab.dynamics.MeridionalDiffusion(K=K, state={'Tatm': diffmodel.Tatm}, **diffmodel.param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.process.time_dependent_process.TimeDependentProcess'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (30, 1) \n",
      "  Tatm: (30, 20) \n",
      "The subprocess tree: \n",
      "Seasonal RCE with heat transport: <class 'climlab.process.time_dependent_process.TimeDependentProcess'>\n",
      "   Insolation: <class 'climlab.radiation.insolation.DailyInsolation'>\n",
      "   Radiation: <class 'climlab.radiation.cam3.cam3.CAM3'>\n",
      "   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   Meridional Diffusion: <class 'climlab.dynamics.meridional_advection_diffusion.MeridionalDiffusion'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "diffmodel.add_subprocess('Meridional Diffusion', d)\n",
    "print(diffmodel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1826 steps, 1826.2110000000002 days, or 5 years.\n",
      "Total elapsed time is 10.998181480672276 years.\n"
     ]
    }
   ],
   "source": [
    "diffmodel.integrate_years(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 365 steps, 365.2422 days, or 1 years.\n",
      "Total elapsed time is 11.997518359050515 years.\n"
     ]
    }
   ],
   "source": [
    "diffmodel.integrate_years(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAf9klEQVR4nO3df7AdZZ3n8feHAAEJapgENoSwiUzICFkHnJhipUYZUKGcKcPMllamVio1skONlVGwdFfQrbGs2WyhY7Hr1i4zdUtwUg7KZoSRrLuigZWxKIEYECQhZglmBi5ECAISZyQxyXf/6L7QXO75eU93P93n86q6dc/p2+f0c053f863n/P0bUUEZmaWrqPqboCZmXXnoDYzS5yD2swscQ5qM7PEOajNzBLnoDYzS5yD2sxsFiQdJ2mrpIck7ZD02Xz6SZK2SHo0/z2/8JhrJO2WtEvSxT2X4XHUZmbDkyTghIj4haRjgLuBK4E/AJ6LiGslXQ3Mj4hPSjoL+BqwGjgVuAM4MyIOd1qGK2ozs1mIzC/yu8fkPwGsATbm0zcCl+a31wA3R8SBiNgD7CYL7Y6OHnWjh7FgwYJYunTpy/d3PfrTUpd3+DiN9vnmjvTpGmvO3I4FQTIOH5hTdxOSMOdACc/5UjlH5yuW/4tX3b///vufjYiFs3nO8y84Ll547khf8z7y8K92AC8VJk1ExERxHklzgPuBXwf+R0TcJ+mUiNgLEBF7JZ2cz74YuLfw8Ml8WkdJBPXSpUvZtm3by/cvuORzpS7v+RXHjvw595/R30o3q9uJj43+QHr+roMjf84pd93+yVfdl/SPs33OF547wle/eUpf857zLydfiohV3ebJuy3OkfRG4O8krewy+0yVYtdPuSSCumpTG9UoA3tq43dgW6rKCGgoN6SbJiJekHQXcAnwtKRFeTW9CHgmn20SWFJ42GnAU92ed6z7qMvYwE587KjSdgizYZS5TTqkQdLCvJJG0vHAu4AfA5uBdfls64Db8tubgbWS5kpaBiwHtnZbxlhW1EXzdx0spSukuGO4yraqVVEsOKRftgjYmPdTHwVsiohvSroH2CTpcuBx4P0AEbFD0ibgEeAQsL7biA9wUAPlhfUUh7ZVpaqjOYf0KyLiR8C5M0z/GXBRh8dsADb0uwwHda7ssJ7i0LZRqrqbzQFdDwd1QVVhPWX6Tubgtl7q/P7DIV0fB/U0VYd1Uaed0AE+flL7QtohXS8H9QzqDOuZOMDbK7VAnolDun4O6g5SC+uZDLqTO9jL14TgHYRDOg1JBvXx258E4Jcru55VWbomhPUg2hYiVi6HdDqS3nOnArtO83cd9AZrYyeVbT6FDEhB0kFtZtVLJaTtFckHdSqfqN54bRyktJ2nsu+nIPmgTklKG7HZqHn7TlcjgtqfrGbjxfv8qzUiqFPiqsPayNt12hzUQ/BGbW2S2vbsavq1GhPUXnlmo5daSNvMGhPUqfEGbjZ6Lshm1qigTm0lOqytybz9NkejgtrMRiPFkE6tEEtJ44I6tZWZ4gZvZu3SuKBOkcPamiTF7TW1Aiw1jQxqr1Sz9vD+3Fsjgxq8cs3awPtxfxob1GY2uBS7Pay3Rgd1Sp/G3gHMBpPS/jsbkpZI+q6knZJ2SLoyn36OpHslPShpm6TVhcdcI2m3pF2SLu61jCSv8GJmo5dSMdGWkM4dAj4eEQ9IOhG4X9IW4PPAZyPiW5Lem9+/QNJZwFrgbOBU4A5JZ0bE4U4LaHRFDa1b4WbWMBGxNyIeyG/vB3YCi4EAXp/P9gbgqfz2GuDmiDgQEXuA3cBqunBFPUJtu8aiWRlSKK5+dngeX3n+7X3OvWmBpG2FCRMRMTHTnJKWAucC9wFXAd+W9AWyonhqgYuBewsPm8ynddSKoD5++5O1XwjXLGUpdXs00LMRsarXTJLmAbcAV0XEi5L+E/CxiLhF0geAG4B3AZrh4dHtuRvf9WFmzZFCNV0GSceQhfRNEXFrPnkdMHX7b3mle2MSWFJ4+Gm80i0yo9YEdSobgCsXs5mlso+OmiSRVcs7I+K6wp+eAt6Z374QeDS/vRlYK2mupGXAcmBrt2W0ouvDzNL3y5WL2xrW5wOXAQ9LejCf9ingj4EvSjoaeAm4AiAidkjaBDxCNmJkfbcRH9CioHYftZnVISLuZuZ+Z4Df6vCYDcCGfpfRs+ujy2DukyRtkfRo/nt+4TEDDeY2M7PO+umjnhrM/WbgPGB9PmD7auDOiFgO3JnfZ9pg7kuA6yXNKaPxKfLwPLPOfOQ7nJ5dHxGxF9ib394vaWow9xrggny2jcBdwCcpDOYG9kiaGsx9z6gbP8UrP037zzhS27JPfKw135ObDdZHPW0w9yl5iBMReyWdnM/W12BuSVeQd66ffvrpAzfc0lJnKM+kU3sc4NZEfW+10wdzd5t1hmmvGcwdERMRsSoiVi1cuLDfZryGq+n67D/jyMs/TdHENreN99nB9VVRdxjM/bSkRXk1vQh4Jp8+8GBua442BVzxtbjStpT1M+qj02DuzWRn3pD/vq0wfaDB3MPyJ3N12l6Ftv31pfYlt/fdwfRTUXcazH0tsEnS5cDjwPthuMHclqY2B1cnrrItRf2M+ug2mPuiDo8ZaDC3pWccQ3q6/WcccVhbEhq7FfrQqRxt7wIYlN+P8ngf7l9rTiG32XEYdTf1/rjCtjo0Mqj9STw6DujBOLCtDt7axpQP6WfH751VqXFB7Wp69hwyo+EPu9nz/tyfRgW1V+rsOVhGz++pla1RQW3Dc/VXLr+/w3MB1ltjgtorc3gOkOr4vbYyNCKoHdLDc3BUz+/54LyPd9fI4XnWm8OiXh7GZ6OU/FbkT9rBOaTT4XXRP+/rnSUd1E1bcSn8hzIHQ3q8TvrXtH2+Kkl2fXhlDceBkC7/g6f+Te3/x29/suaWpMNbTks4pNPndTSYX65c3IiiTdISSd+VtFPSDklXFv72EUm78umfL0y/RtLu/G8X91pGkhV1E9XZ7eEAaI46K+vnVxzL/F0Ha1l2yx0CPh4RD0g6Ebhf0hbgFLKLfb8lIg5MXVdW0lnAWuBs4FTgDklndvu//a6oG84h3TxeZ+0SEXsj4oH89n5gJ9kFvT8MXBsRB/K/TV2ucA1wc0QciIg9wG5gdbdluKI2q0FdlbWr6syLB4/jO4+v6Hf2BZK2Fe5PRMTETDNKWgqcC9wH/AXw25I2AC8Bn4iIH5CF+L2Fh03m0zpyUI9AXd0ersyazV8wNsazEbGq10yS5pFdBPyqiHhR0tHAfOA84G1kly58EzNfMSu6PbeDuqHaFtLzlv28499+secNFbbEbHCSjiEL6Zsi4tZ88iRwa0QEsFXSEWBBPn1J4eGnAU91e34HdQM1OaS7BfJsHtPUMHdV3XySBNwA7IyI6wp/+gZwIXCXpDOBY4Fngc3AVyVdR/Zl4nJga7dlOKhnqepuj6aF9DDBPIrlNCm4qw5r91OP3PnAZcDDkh7Mp30KuBG4UdJ24CCwLq+ud0jaBDxCNmJkfbcRH+CgbpSmhHRV4dxvG5oQ2q6smysi7mbmfmeAD3Z4zAZgQ7/LcFA3ROohnUI4d9KU0HZYWyfeKmYhhf/tUbd5y36edEhP17T2lsXbbrO4om6AFKvppofdVPtTq7BdVdtMvEUMqaqKJLWQbltFmuLrqWqdu6puDge19S21QBulFAPbbIqDegjjVk2PU4il8jpdVVuR+6gTlUJIpxJaVUul/9r91TbFW8GAxqUCGdeQLhqX92Bctukmc1AnqO5qelwCqh91vxd1bwuWBgf1AKqoPOreMesOphTV3UdfxTbhqjptDmoD6g+jJvD7Y3VxUPepzdW0A6h/db1XdR9pWb0c1GPOIT24toa1uz/S5aBORB0Vk0N6eH7vrEoO6j60sdJw0MxeHe+hq+rx1HdQS5oj6YeSvpnfP0nSFkmP5r/nF+a9RtJuSbskXVxGw9uk6mraIT06fi+tCoNU1FeSXQZ9ytXAnRGxHLgzv4+ks4C1wNnAJcD1kuaMprnVa1uF4WAZvarfU1fV46evoJZ0GvC7wJcKk9cAG/PbG4FLC9NvjogDEbEH2A2sHklrW6jKatohXR6/t1amfivq/wr8B6CYKqdExF6A/PfJ+fTFwBOF+Sbzaa8i6QpJ2yRt27dv36DtrkTZlYWHXLVLlWHtqnq89AxqSb8HPBMR9/f5nDNdOyxeMyFiIiJWRcSqhQsX9vnU1WnbhuqKrxp+n60M/VTU5wPvk/QPwM3AhZL+Bnha0iKA/Pcz+fyTwJLC408DnhpZi1vCXR42W66qx0fPoI6IayLitIhYSvYl4f+NiA8Cm4F1+WzrgNvy25uBtZLmSloGLAe2jrzlJWrTBuqQrp7fcxu12YyjvhZ4t6RHgXfn94mIHcAm4BHgdmB9RByebUPbpKpq2oFRn6ree1fV9ZO0RNJ3Je2UtEPSldP+/glJIWlBYdpAQ5gHunBARNwF3JXf/hlwUYf5NgAbBnnuVLRlw3RI12/esp/XfvEBq8Qh4OMR8YCkE4H7JW2JiEckLSErZB+fmnnaEOZTgTskndmtoPWZiRXzSI/xUsUHpqvqekXE3oh4IL+9n+x8k6mRbv+FbMRccUDFwEOYHdQFbdkgXU2b1UPSUuBc4D5J7wOejIiHps3W1xDmIl8zsWUc0umpoguk7OsrPr/iWObvOlja81ft8IE5g6yTBZK2Fe5PRMTE9JkkzQNuAa4i6w75NPCeGZ6vryHMRQ7qXJv/37TVz/3VjfZsRKzqNoOkY8hC+qaIuFXSvwKWAQ9JgmyY8gOSVjPEEGZ3fbSIq2krU1u6BkdNWRLfAOyMiOsAIuLhiDg5IpbmQ5sngbdGxE8ZYgizg5p2VNMO6fSVvY58xFab84HLyE4GfDD/eW+nmYcZwuyuDzPrW9v6qkchIu5m5n7n4jxLp90faAjz2FfUrqatSq6qbRhjH9RmVWv6B6v7qqs31l0frqbT8J7Tdw38mO88vqKElrRD2UP1rHpjG9RtqAqaGtLDBHOv52hacDd9uJ77qqs1tkFt1RlFMA+yjKaFdlM5rKszlkFdVTVdZrdHE6rpKgK623JTD+wyq2p3f7TL2K3JNnR5NEFdIT29DSm0o5smfOB24/2pGmMV1FVuVONaTacYjqm1pypVDdVzWJdvLLs+mizVkE49DFPuDmn6F4tWvrEJan/qlyP1gJ4u5cBuMn+xWK6x6PqoOqTH5eywpoV0UWptL+tIqcpt0cVQecYiqNsipW6P1IJuGG14DalxWJej9UHtDWf02hRwKb2WNlTVVo5WB3UdIV3WTpFKNZ1SsI1KiiNVmszF0ei1Nqi9sYzWOIRZCq8vlQ/k2fL+N1qtDOq6NpK2VtMpBFhV2vpa3f3RbK0Lan+Sj1Zbg6ubcXzNZfC+ODqtCuo2bhh1VtPjHFh1vva6j6BG6fkVx7Zyv6xaa4K67o2hbYeW4xzSU9r2HtS5jda9fzZdK4K6rRtBmyqrpqorrNu47tu6n1ah0UHtw6pytK2StHR4fx1OY4M6pRXepm4Ph/RrtamqTmFbTWnfbYrGBfW4VNF1HPo6pDvzezNabdqHJS2R9F1JOyXtkHRlPv0vJP1Y0o8k/Z2kNxYec42k3ZJ2Sbq41zIaEdRT4dymlZsaB1Fvfo9Gq0X78yHg4xHxZuA8YL2ks4AtwMqIeAvw/4BrAPK/rQXOBi4Brpc0p9sCkg3qpoRzGYeSbfwiqS2qDuu2dn9Macp+3k1E7I2IB/Lb+4GdwOKI+E5EHMpnuxc4Lb+9Brg5Ig5ExB5gN7C62zKS/H/UTV5pTeRK0VJQ3O+b+r+tJS0FzgXum/anDwH/M7+9mCy4p0zm0zpKMqjHWdXVtEN6cO85fVelFx4YxyvATIV2WYE95wCDXPx3gaRthfsTETExfSZJ84BbgKsi4sXC9E+TdY/cNDVphmVEtwY4qGchpUPIYTikh1d1WI9aU65SnsjR9bMRsarbDJKOIQvpmyLi1sL0dcDvARdFxFQYTwJLCg8/DXiq2/Onv6bMEuUPOgOQJOAGYGdEXFeYfgnwSeB9EfHPhYdsBtZKmitpGbAc2NptGa6oE1Jlt4dDplnGsfujQc4HLgMelvRgPu1TwH8D5gJbsizn3oj4k4jYIWkT8AhZl8j6iDjcbQF9BXU+/u9LwEqyvpQPAbvIOseXAv8AfCAins/nvwa4HDgMfDQivt3Xy22QJnd7OKRHp8ldIE3p/khdRNzNzP3O/6fLYzYAG/pdRr9r6YvA7RHxG8Bvkg0/uRq4MyKWA3fm94caI2hmZp31DGpJrwfeQdYHQ0QcjIgXyMYCbsxn2whcmt8eeIygVdft4Wp69Kp6Tz2+fnz1U1G/CdgHfFnSDyV9SdIJwCkRsReyAd/Ayfn8i4EnCo+fcYygpCskbZO0bd++fbN6EVVrcreHWZG35WboJ6iPBt4K/GVEnAv8E3k3Rwd9jRGMiImIWBURqxYuXNhXY212XE2Xx1W1lamfoJ4EJiNi6kybr5MF99OSFgHkv58pzD/QGMFx552vHfxBaGXpGdQR8VPgCUlTX21fRDasZDOwLp+2Drgtvz3wGMEmaeqhokPEOmnqNj1O+h1H/RHgJknHAj8B/ogs5DdJuhx4HHg/wDBjBK1cDunqVDFcz2Oqx09fQR0RDwIznUJ5UYf5BxojOM7c7WFmvXi0+wCaeIjoarp6TXzPm7htjxMHtVkD+UhsvDioa1T2ztbEyq4t/N7bKDmo++RDQ2s7b+PpclC3lCu6+nkd2Kg4qGviPkabLW9D48NB3UKu5NLhdWGj4KDug/vubFx4W0+Tg7oGZR6yuoJLT5nrxN0f48FBbWaWOAd1Dz4UtFFo0pGOt/n0OKgr5m4PGzV3f7Sfg9rMLHEO6i6adAjoajp9TVpHTdr2x0G//4/arDSXzf/+wI/5yvNvL6ElzeX/Ud1urqg7KKOicF/iKy6b//2Xf2bz+KZpUlVt/ZG0RNJ3Je2UtEPSlfn0kyRtkfRo/nt+4THXSNotaZeki3stw0HdAk3a+UcdsE0N7CZw90ffDgEfj4g3A+cB6yWdRXYR8DsjYjlwZ36f/G9rgbOBS4DrJc3ptgB3fVjpqgjS4jLcLWJVioi9wN789n5JO4HFwBrggny2jcBdwCfz6TdHxAFgj6TdwGrgnk7LcFDPoEndHqlX03VUu5fN/37SYV3WdRXdT92/OS8F83cd7Hf2BZK2Fe5PRMTETDNKWgqcC9wHnJKHOBGxV9LJ+WyLgXsLD5vMp3XkoLbS1NklMbXslAO7CfafcYQTHxv7HtJnI2Kma8a+iqR5wC3AVRHxoqSOs84wLbo999ivARu9lPqNU2nHdKkfCdlgJB1DFtI3RcSt+eSnJS3K/74IeCafPgksKTz8NOCpbs/voJ6mSV+gpLizpxiMKbapLB5ZVD1lpfMNwM6IuK7wp83Auvz2OuC2wvS1kuZKWgYsB7Z2W4aDugLjsvOkHIgpty11TSpeanI+cBlwoaQH85/3AtcC75b0KPDu/D4RsQPYBDwC3A6sj4jD3RbgPuqGSq2abkIQpvYlY1lfKlq1IuJuZu53Briow2M2ABv6XYYr6gJXDsNpQkhPSan/3KxfDmqblaaGXirtLuPIqIyuNhcx9XJQl6yMnSaVbo9Uws6s7RzUOVcMg2lDSLfhNdh4cFDbwNoUcCm8llSOkHpxMVMfB3WJ2tzt0SYphPWojcuQ0HHhoMaVwiDaGGpmqXNQW9/aHNJ1v7amHCm5qKmHg7pB6tyZ6w6yKozDa7RmclCXpE19hOMUYG16rW3aBsfd2Ae1D+UsFe7+sE7GPqiboq6duE0VZr/G8TVb2sY6qF0ZdDfOgTXOr93S01dQS/pYfnXd7ZK+Jum4UV5ht23cN2ip8LbYDj2DWtJi4KPAqohYCcwhu4LuyK6wW4cmVdN1dHu4oqznPXA/tc2k366Po4HjJR0NvI7ssjFryK6sS/770vz2y1fYjYg9wNQVds0axx9YloKeQR0RTwJfAB4nuyT6zyPiO0y7wi5QvMLuE4WnmPEKu5KukLRN0rZ9+/bN7lUMyNVAdw4n64f3o+r00/Uxn6xKXgacCpwg6YPdHjLDtNdcYTciJiJiVUSsWrhwYb/tNatc1R9co+7+cD918/XT9fEuYE9E7IuIXwG3Am9nhFfYrVLZVcCod4qq+yxdTdsgXFVXo5+gfhw4T9Lr8qvtXgTsZIRX2DVLnT/ArE49L24bEfdJ+jrwAHAI+CEwAcwDNkm6nCzM35/Pv0PS1BV2D9HHFXar4k//7hxGNoz9ZxzhxMfG+pSM0vX17kbEZyLiNyJiZURclo/o+FlEXBQRy/PfzxXm3xARZ0TEioj4VnnNb7cquz0c0r1V+R41ZZieVWNsPgarqKb9pY2lquxtc5yPViXdKOkZSdsL086RdK+kB/PRbasLfxv4hMCxCWrrzNV0//xe2Qz+muzkvqLPA5+NiHOAP8vvD31CoIPaLFFN6/4Y16o6Ir4HPDd9MvD6/PYbeGXk21AnBPb8MrENmrgBVbWTukIc3GXzv89Xnn973c2wWTjql7/i+O1P9jv7AknbCvcnImKix2OuAr4t6QtkBfHUBrMYuLcw34wnBE43FkFdBfdPW+rmLfs5v9jzhlKX0dIRIM9GxKoBH/Nh4GMRcYukDwA3kJ2T0tcJgdO17h2dronVdFVcTQ/P7531sI7s5ECAv+WV7o2hTghsfVA3UdP6Jq08TdwWXBwBWfi+M799IfBofnuoEwJb3fXhDaYzV4Sz575qA5D0NeACsr7sSeAzwB8DX8z/4+hLwBUw/AmBrQ1qh7SZVSEi/rDDn36rw/wbgA2DLMNdHyPgLxLHV9OOTKraVl0ojVYrg7rJG0kVfZJNC5dx18R+ahut1gV1k0ParE28L45Oq4LaG0ZvrqZHz+9pZ94nR6NVQV0H909b01S9zTqsZ681Qd2GjcF9kc1VdlXtbWO8tSKo2xDSVfAhutXF++jstCKozSx9DuvhNT6ovfItFT5i6c3763AaHdR1r/RRfilTdh+kQ6T5RrmN+EvwZmlsUNcd0mY2HO+7g2tkUHtFD8bVdHX8XvfH+/BgGvVPmbxyzWwcNaKi3n/GkeRC2n181kmZVXWb+qlT26dTlnxQj8PKLPOLRB+KW8pSLMJSlGxQewWajQ/v690lGdReadZ0PpIZnIuzzpIMahsNh4U1kQP7tRzUQ2jSiS7WPm36QrEbh/UrHNQt5Wq6fl4Hs+ewzjiozcwS56A2M0ucg7qFfMidjrLWhb/bSIekGyU9I2n7tOkfkbRL0g5Jny9Mv0bS7vxvF/ezDAf1gPxForVNyl8oNsRfA5cUJ0j6HWAN8JaIOBv4Qj79LGAtcHb+mOslzem1AAe1mdksRMT3gOemTf4wcG1EHMjneSafvga4OSIORMQeYDewutcyFBEjbPJwJO0HmlZeLgCerbsRA2pam5vWXnCbq7AiIk6czRNIup3sdffjOOClwv2JiJiY9nxLgW9GxMr8/oPAbWRV80vAJyLiB5L+O3BvRPxNPt8NwLci4uvdGpDKf8/bFRGr6m7EICRtc5vL1bT2gttcBUnbZvscEXFJ77lm5WhgPnAe8DZgk6Q3AZqpOb2ezF0fZmajNwncGpmtwBGyCn4SWFKY7zTgqV5P5qA2Mxu9bwAXAkg6EziWrHtpM7BW0lxJy4DlwNZeT5ZK18dE71mS4zaXr2ntBbe5Ckm1V9LXgAuABZImgc8ANwI35kP2DgLrIvtCcIekTcAjwCFgfUQc7rmMFL5MNDOzztz1YWaWOAe1mVniag1qSb8p6R5JD0v6X5JeX/jbwKdZVmWUp4aWTdKfS/qRpAclfUfSqYW/JdfeKZIuydu1W9LVdbdnOknHSdoq6aF8O/hsPv0kSVskPZr/nl93W4skvVHS1yX9WNJOSf865TZLulLS9vw9viqflmx7SxMRtf0APwDemd/+EPDn+e2zgIeAucAy4DFgTp1tLbT5d4A7gLn5/ZNTbjPw+sLtjwJ/lXJ787bNydvzJrJvyx8Czqq7XdPaKGBefvsY4D6yMbOfB67Op18NfK7utk5r90bg3+W3jwXemGqbgZXAduB1ZAMf7iAbJZFke8v8qbvrYwXwvfz2FuDf5LeHOs2yIiM9NbRsEfFi4e4JvDK4Psn25lYDuyPiJxFxELiZrL3JiMwv8rvH5D9B1s6N+fSNwKXVt25m+RHrO4AbACLiYES8QLptfjPZWXz/HBGHgL8Hfp9021uauoN6O/C+/Pb7eWUg+GLgicJ8k/m0FJwJ/Lak+yT9vaS35dOTbbOkDZKeAP4t8Gf55GTbS9pte5mkOfmpws8AWyLiPuCUiNgLkP8+ucYmTvcmYB/wZUk/lPQlSSeQbpu3A++Q9GuSXge8lywjUm1vaUoPakl35H1M03/WkHV3rJd0P3Ai2XhDGPI0y4raXDw19N+TnRqqOtvco71ExKcjYglwE/CnUw+rq719SLltL4uIwxFxDtnZZaslray5Sb0cDbwV+MuIOBf4J7KugyRFxE7gc2RH27eTdYEdqrVRNSn9hJeIeFePWd4DL5+987v5tKFOsxyVbm2W9GHyU0OBrZJmdWroKPTxHk/5KvC/yQbk1/oe95By214jIl6QdBfZP+B5WtKiiNgraRFZtZ2KSWAyr/wBvk4W1Mm2OSJuIO+qkfSfyV5Dsu0tS92jPk7Ofx8F/Efgr/I/DXWaZUW+wQhPDS2bpOWFu+8DfpzfTrK9uR8AyyUtk3Qs2f/v3Vxzm15F0kJJb8xvHw+8i+y93Qysy2dbR/Yf1JIQET8FnpC0Ip90EdkZcsm2uZARpwN/AHyNhNtblrpPIf9DSevz27cCXwaIiKFOs6zISE8NrcC1+Y55BPhH4E8g7fc4Ig5J+lPg22QjQG6MiB01N2u6RcBGZf/0/ShgU0R8U9I9ZN1hlwOPk333kpKPADflH4A/Af6IvP2JtvkWSb8G/IpsG31e0rWk295S+BRyM7PE1T3qw8zMenBQm5klzkFtZpY4B7WZWeIc1GZmiXNQm5klzkFtZpa4/w/fEy5lCd/iiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAio0lEQVR4nO2df9AdVZnnP19IADFhxckAIQmVOBUUtFTcdyiU0nHAEQpd4o/CxVqZjLLDOhs1mcVVgq7ulJMtRp3MMLXlWNlBK9ZEMUOyQ0YZMSBqUQvBJBMnCS+M0SgGYiA7KGGRQMizf/S5Sefm/uh+33vffrrv86l66+0+9/Ttc7vP+bznPfecp2VmBEEQBH45oeoCBEEQBL0JUQdBEDgnRB0EQeCcEHUQBIFzQtRBEATOCVEHQRA4J0QdBEEwCSTNk3SPpHFJOyUtTemvkXSfpO2S/kHSabljlkvaJelhSZf1PUfMow6CIJg4kmYDs81sq6SZwBbgHcBq4KNm9j1JHwAWmNl/k3Q+8DXgQuBs4C7gXDN7ods5okcdBEEwCcxsr5ltTdsHgHFgDvBy4Psp20bg3Wl7EXCrmR00s93ALjJpd2XaMApellmzZtn8+fOP7P/LQ3s75jt8cjP/rrwwveoSBHXjxOerLsHU8YoFZx6zv2XLlv1m9puTec83vfkUe/JfDxfKu2P78zuBZ3NJq8xsVae8kuYDFwCbgB3AlcDtwFXAvJRtDnB/7rA9Ka0rLkQ9f/58Nm/efGT/997wp6WOP7DgRYMuUuU8PaeZf5SC8sx4tJhQ6srM3b/u+frG//PJY/Yl/Wyy53zyXw+z/o5ZhfKeO2/vs2Y21i+fpBnAOmCZmT2Vhjv+StKngA3Ac62sHQ7vOQbtQtSTpdeNrqvEizTOkHn9abqEob+Im4Ck6WSSXmNm6wHM7CHgren1c4G3pex7ONq7BpgLPNbr/Rsh6l50qyR1FXieiTTykPvwGAXp9mIUhNwJSQJuAcbNbGUu/Qwze1zSCcAngS+mlzYAX5W0kuzLxIXAA73O0XhRd6PJAu/FVMrEyx+FURfoIBlVGffhYuAaYLukbSntRmChpCVpfz3wZQAz2ylpLfAgcAhY0mvGB4ywqLsxqgIfBiHI+hJCLo6Z3UvncWeAm7scswJYUfQcIeqCNHEcPBhtQsb1IUQ9AIpU+JB5MJWEhJtFiHqKmGjDCcEHId0gRO2cKhtp/JE4lhBmUBUh6qArIaYg8IGP+VNBEARBV0LUQRAEzglRB0EQOCdEHQRB4JwQdRAEgXNczvqY9uNH++Y59Fs9w7cGQRA0BpeiLkIvmYfEg8A3RTpjwVFqK+pedKsEIfAgmDpCxoOjkaLuRgg8CAZLyHhqGClRdyOGUYKgOyHj6glR9yG+2AyaTEh48kiaB3wFOAs4TPbw25slvZbsqS6nkD0g4D+b2QPpmOXAtcALwEfM7M5e5whRD4BBVPapkn00zMHh5Q903NPKOQRcb2ZbJc0EtkjaCHwW+BMz+0dJV6T9N0s6H7gaeCXZo7juknRur6e8hKidEI2tfsQ9CwDMbC+wN20fkDQOzCF7svhpKdu/4egDbBcBt5rZQWC3pF3AhcB93c4Rog6CIBgQkuYDFwCbgGXAnZI+T7a48A0p2xzg/txhe1JaV0LUQRCMHL86fAp3PH1+wdx7Z0nanEtYZWar2nNJmgGsA5aZ2VOS/hT4YzNbJ+k9ZE8qfwudn69ovUoQog6CIOjNfjMb65VB0nQySa8xs/UpeTGwNG3/HfA3aXsPMC93+FyODot0JGJ9BEEQTAJJIustj5vZytxLjwG/k7YvAX6UtjcAV0s6WdICYCHwQK9zRI86OIKXWQxNIL5oHCkuBq4BtkvaltJuBP4QuFnSNOBZ4DoAM9spaS3wINmMkSW9ZnxAiLqRhHCrp+w9CLHXFzO7l87jzgD/tssxK4AVRc8Roq4pIeNmkb+fIe2gnRB1jQg5jwbt9znEHYSonRNyDlp1IIQ9uoSonRKCDtoJYY8uIWonhJiDosR49ugRop5iQsjBIOlWn0LgzSJEPSRCyEGVdKp/Ie/60lfUPWKtvhT4OjAf+CnwHjN7Mh1TKtZq3QkpB3UgZpPUlyI96m6xVv8AuNvMbpJ0A3AD8PGJxFqtGyHmoAmEuOtDX1H3iLW6CHhzyrYa+C7wcSYQa9UzIeVgVIjhEr+UGqNui7V6ZpI4ZrZX0hkpW6FYq5KuI619P+ecc0oXfBiElIPgWKLX7YPCou4Qa7Vr1g5px8VaTfFcVwGMjY31jMU6DELKQVCe6HVXQyFRd4m1uk/S7NSbng08ntJLx1odNiHlIBgeIe/hU2TWR7dYqxvIAmPflH7fnkv/qqSVZF8m9o21OghCxkHgh5jfPViK9Ki7xVq9CVgr6VrgEeAqmFis1W6EfIOgWRRt0yH0Yyky66NXrNVLuxxTKtZqOyHoIBhtwgHHEo/iCoIgmASS5km6R9K4pJ2Slqb0r0valn5+mhuRQNJySbskPSzpsn7niCXkQRAEk6PjokAz+/etDJL+HPhV2i69KDB61EEQBJPAzPaa2da0fQBoLQoEjkzIeA/wtZR0ZFGgme0GWosCuxI96iAIRo6nDp3Cd/a/omDuu2dJ2pxLWJXWgRxH26LAFm8E9plZ6ynkhRYF5glRB0EQ9Ga/mY31y9S+KDD30ns52puGgosC84SogyAIJkmXRYFImga8i2OfRl56UWCIOhg4Bxa8qOoiHGHm7l9XXYSg4fRYFAjwFuAhM9uTSyu9KDBEHZTGk4j70a+sIfJgAHRcFGhmd5DN7sgPe0xoUWCIOuhKnYQ8Ubp9xhB4UJReiwLN7A+6pJdaFBiiDo4wCmIuSqdrEfIOqiJEPaKElMvTfs1C3MFUEaIeEULMgyfEHUwVIeqGEmKeekLcwbAIUTeIkLMvQtzBoAhR15gQc73I36+QdlCGEHXNCDk3g5B2UIYQdQ0IOTebkHbQjxC1U0LOo0lIO+hEiNoRIecgT0g7aBGirpiQc1CEkPZoE6KuiBB0MFFC2qNHiHoKGUU5Pz2nmqe9zXj0cCXnnWpC2qNBiHrIjIKcq5JxL3qVqakSb9W1EHbzCFEPiSYL2qOYy9Be/qaJO3rZzSNEPUCaKue6i7kfTRZ39LKHj6R5wFeAs4DDZA+/vTm99mHgQ2QPCPimmX0spS8HrgVeAD5iZnf2OkeIegA0TdBNF3M/miju6GUPlUPA9Wa2VdJMYIukjcCZwCLg1WZ2UNIZAJLOJ3vyyyvJHsV1l6Rzez3lJUQ9QULOo0P+2jRJ2iHswWBme4G9afuApHFgDvCHwE1mdjC99ng6ZBFwa0rfLWkXcCFwX7dzhKhL0iRBh5zL0yRpj3Iv+9nnpzP+izOLZp8laXNuf5WZreqUUdJ84AJgE/A54I2SVgDPAh81sx+QSfz+3GF7UlpXQtQFCDkHnWiitEdN2AXZb2Zj/TJJmgGsA5aZ2VOSpgGnAxcBvw2slfQyOj9f0Xq9d4i6C02SM4Sgh03r+jZF2BDSLoOk6WSSXmNm61PyHmC9mRnwgKTDwKyUPi93+FzgsV7vH6Juo0mCDjlPPU3sZUNIuxeSBNwCjJvZytxLfw9cAnxX0rnAScB+YAPwVUkryb5MXAg80OscIWqaJWcIQXuhKb1siKGRPlwMXANsl7Qtpd0IfAn4kqQdwHPA4tS73ilpLfAg2YyRJb1mfMAIi7ppcoYQtFeaKGwIabcws3vpPO4M8L4ux6wAVhQ9x8iIuolihpBznWjSsAjEMyGnkkaLuqlyhhB03WlSL7tF9LaHR2NE3WQp5wlBN4smChs6t8eQ98SplahHRcadCEE3m6YKO0/Ie+IUFrWkE4HNwKNm9nZJLwW+DswHfgq8x8yeTHlLBRxpZ5SF3E4IerQYBWHn6dfWQ+QZZSywFBjP7d8A3G1mC4G70357wJHLgS8kyQcleHrOCSHpESbuf0Z02jIK9aglzQXeRjad5L+k5EXAm9P2auC7wMeZQMCR4ChNbpzPzB1uL/HUPc27dqPWww46U3To4y+BjwEzc2lnpqhRmNneVgg/CgYckXQdcB3AOeecU67UDaRJgh62kMuetwkCf3rOCSHrEaavqCW9HXjczLZIenOB9ywUcCRFn1oFMDY21jMgSdOps6SrknIZOpWxjvKO3vXoUqRHfTFwpaQrgFOA0yT9LbBP0uzUm54NtGKtlg44MqrUVdB1kHM/2j9DncQdwh49+tZOM1tuZnPNbD7Zl4TfMbP3kQUWWZyyLQZuT9sbgKslnSxpAQUCjowadfui6Jm5h4/5aSJ1/Hx1q0fBxJnMPOqbyOKrXgs8AlwFYGalA46MCnVqVHUS1qDJf/Y69LSjh918SonazL5LNrsDM/u/wKVd8pUKONJ06iLoUZZzN+ok7fjCsbnUamViHamDpEPQxWhdJ8/Cjt51MwlRD4kQdHOpQy87hN0sfNayGlOHL3jq9qWZZ7xfS+91sQlImifpHknjknZKWprS/7ukRyVtSz9X5I5ZLmmXpIclXdbvHNGjHiDeG4VnodQdz8Mi0bseOoeA681sq6SZwBZJG9Nrf2Fmn89nbguzcTZwl6Rze026CFEPgBB00OKZuYddyhpC2MMirdBurdI+IGmcDquxc5QOsxGingQh6KATnnvXELNDAOz5Ezj02KlFs8+StDm3vyqtrD4OSfOBC4BNZIsFPyTp98kij16fIowWCrORJ0Q9AULQQRE8Czt616XYb2Zj/TJJmgGsA5aZ2VOS/hr4DFkIjc8Afw58gIJhNvKEqEviWdIhaJ+EsJuPpOlkkl5jZusBzGxf7vX/BXwj7ZYOs+Gv5jjF82wO7zMPggzP98lz/faOJAG3AONmtjKXPjuX7Z3AjrRdOsxG9Kj74L3yem34QXe8f+EYvevSXAxcA2yXtC2l3Qi8V9JryYY1fgr8J5hYmI0QdRdC0MEwieGQ5mBm99J53PmOHseUCrMRom4jBB1MJSHsoAgh6oR3QUNIusl4Hw6BEHaVhKjxL+kQ9GjguXcNIewqGVlRe5dzi1GR9LSznymct8RChVriuXcNIewqGDlR10XQ0BxJl5HwIN6vCSL3LmsIYU8lIyPqEPTUMWgxT/b8dRW396GQFvm2FdIeDo0WdZ3k3KKOkq5azP2ou7jr0LtuEb3s4dA4UddRzlA/QXuXcy/yZa+LtOvSu24RvezBUntR11XMeeoi6TrLuRt1k3adetctQtqTp3aiboKY83iXdBPl3I26SLuOsm7R3n5D3MVwKeqmybgTIWjftD6/V2HXbSikG93aegj8WFyKuul4lvSoC7qdOgi77rLuxCh01soQV2OK8SrpaWc/E5Lugefr47VOBYMjetRThNfG5FU+XvHaw27KUEjQmbirU4BHSXvuIdYBr9fPY10LJk+Iesh4bDgeBVNXPArbY51rMpLmSbpH0riknZKWtr3+UUkmaVYubbmkXZIelnRZv3PE0McQ8dZgvAmlSUw7+xlXwyFN/ZLRKYfInjC+VdJMYIukjWb2oKR5wO8Bj7QySzofuBp4JXA2cJekc3s95SXu5BDw9mw8j72+JuLtOnurh03FzPaa2da0fQAYB+akl/8C+BjHPmV8EXCrmR00s93ALuDCXueIHvWA8dYwPIljVIjetX9OeK7UF6+zJG3O7a8ys1WdMkqaD1wAbJJ0JfComf0we/7tEeYA9+f293BU7B0JUQ8QT5IOQVeLt9khIetJsd/MxvplkjQDWAcsIxsO+QTw1k5ZO6RZh7QjxJ0bAN7+xQxJ+8HTcIinOto0JE0nk/QaM1sP/BawAPihpJ8Cc4Gtks4i60HPyx0+F3is1/uHqCeJp8rvSQrBsXi5L946FU1A2bjGLcC4ma0EMLPtZnaGmc03s/lkcn6dmf0C2ABcLelkSQuAhcADvc4Rop4Eniq8FxEE3fF0jzzV3QZwMXANcImkbennim6ZzWwnsBZ4EPgWsKTXjA+IMeoJ4amSe2r8QX88jV3HasbBYGb30nncOZ9nftv+CmBF0XPEHSpJSDoYBJ7unac6HXQmRF0CTxXaU0MPJoane+ipbgfHU0jUkl4i6TZJD6Vlkq+X9FJJGyX9KP0+PZe/1PJI73j7AsZTAw8mh6d76amOB8dStEd9M/AtM3sF8BqylTc3AHeb2ULg7rTfvjzycuALkk4cdMGnCk+VN2Z1NBNP99VbpyTI6CtqSacBbyKbfoKZPWdmvyRbBrk6ZVsNvCNtl14e6RFvFdZLQw6Gh6d77KnuB8VmfbwMeAL4sqTXAFuApcCZZrYXsrXuks5I+Qstj5R0HXAdwDnnnDPhDzAMvFVSTw14kJx31r5JHT/+izMHVBI/eFp+HrNC/FBE1NOA1wEfNrNNkm4mDXN0odDyyLRWfhXA2NhYz+WTU0lIevBMVshl37fuAvcka4jl5x4oIuo9wB4z25T2byMT9T5Js1NvejbweC5/qeWRHvAmaKivpIcl5omev47i9ihriN51VfS96mnJ488lvTwlXUq2omYDsDilLQZuT9ull0dWTUh68px31r4jP97wXLZeeKwDHtvKKFB0ZeKHgTWSTgJ+AryfTPJrJV1LFhT7KsiWR0pqLY88RIHlkVXhsdJ5bJzdqJv44Ngy16Gn7WklY4voXU89hURtZtuATmH+Lu2Sv9TyyKnGo6ChPpKuo6A70focdRG2J1lDjF1PJSMV68OroKEekm6KoNupi7C9yhqidz1sRkLUngUN/iXdVEG3U4dhEY+yhhD2sGn8VQ1JT5w6fgE3KDx/bs91xttCsabQyB51XSqK1wbnWVJTiechEa896xbRwx4sjbmKrb/kIemJM8o96F54vS4e61A7dWqTE0XSPEn3pIB1OyUtTemfkfTP6UEC35Z0du6YUoHrai/qOlYEjw3Mo4i84fEaeaxLnahbR6okh4Drzew84CJgSQpO9zkze7WZvRb4BvApmFjgutoNfdT9RntrWB7l4xmPwyHeh0HaybfhJgyNpJhHrbhHBySNA3PM7MFcthdzNJTGkcB1wG5JrcB193U7h2tR113K7XiSdAh6cngTdt1k3aK9jU+VuE98HmY8WtgvsyRtzu2vSrGKjkPSfOACYFPaXwH8PvAr4HdTtkKB6/K4FHXTBA0h6aZy3ln7QtYDxGnb329mnRb8HYOkGcA6YJmZPQVgZp8APiFpOfAh4NMUDFyXp/7/d9SAkHSz8XRNPdW1UULSdDJJrzGz9R2yfBV4d9ouHbguRD1kvDQcrzMXmoKn6+ulzo0KkkT2YJVxM1uZS1+Yy3Yl8FDaLh24zuXQR1Pw0mC8CGQU8DIU0oRhkBpxMXANsF3StpR2I3Btijp6GPgZ8EGYWOC6EPWQCEmPLiHr0cLM7qXzuPMdPY4pFbguRD0EPEg6BF0tXmaFhKybQYxRD5iQdJDHw73wUCeDyRGiHiAeGoQHMQTH4uGeeKibwcQJUQ8IDw3BgxCCzni4Nx7qaDAxQtQDwEMD8CCCoDce7pGHuhqUJ0Q9STxUfA8CCIrh4V55qLNBOULUk8BDhffQ8INyeLhnHupuUJwQ9QTxUNE9NPhgYni4dx7qcFCMEHVJpp39jIsK7qGhB5PDw7JzD3U56E+IugQeKrWHxh0Mlqrvp4d6HfQmRF2QqMzBMAlZB72IJeQF8FKJq27Mw+KSWQ/1z5TjO/tfMaSSjDax3NwvIeo+hKQHR1khl32fugvcQzCnkLVPQtQ9CElPnkHJeSLnqqO4Q9ZBJ0LUXQhJT5yplHMvWuWom7C9yBoIYTshRN2BkHR5vMi5E/my1UXaHmQN0bsugqR5wFeAs8geErDKzG6W9Dng3wHPAT8G3m9mv0zHLAeuBV4APmJmd/Y6R8z6aCMkXY5LZj3kWtLt1Km8XuqAlzbhmEPA9WZ2HnARsETS+cBG4FVm9mrgX4DlAOm1q4FXApcDX5B0Yq8ThKgTXhaygJ8G2os6Ca8TdSm/l7rgpW14xMz2mtnWtH0AGAfmmNm3zexQynY/2UNsARYBt5rZQTPbDewCLux1jhj6wFcl9NIwu1EHuZWhDuPYMQwyeE44eJiZu39dNPssSZtz+6vMbFWnjJLmAxcAm9pe+gDw9bQ9h0zcLfaktK6MvKhD0sVomqDb8S7skHWl7DezsX6ZJM0A1gHLzOypXPonyIZH1rSSOhxuvd57pEUdku5P0wXdjmdhe5I1xIyQPJKmk0l6jZmtz6UvBt4OXGpmLRnvAeblDp8LPNbr/UdyjNrTeDSEpD3i9bN7qiue2lCVSBJwCzBuZitz6ZcDHweuNLP8xdoAXC3pZEkLgIXAA73OMXI96qhc/fEqqanGc+/aCyM6FNLOxcA1wHZJ21LajcBfAScDGzOXc7+ZfdDMdkpaCzxINiSyxMxe6HWCkRG1V0F76iFBSLoTl8x6yJWsvQyBtBj1oRAzu5fO48539DhmBbCi6DkKDX1I+mNJOyXtkPQ1SadIeqmkjZJ+lH6fnsu/XNIuSQ9LuqxoYYZFSLoYIenueLs23uoO+G1nTaCvqCXNAT4CjJnZq4ATySZr3wDcbWYLgbvT/oQmcw8Tr5XHU0Ory5ziqvF2nTzVoRbevv9pCkW/TJwGvEjSNOBUsm8oFwGr0+urgXek7dKTuYeB5wrjqYF5Ek9d8HTNPNWlPF7bXl3pK2ozexT4PPAIsBf4lZl9GzjTzPamPHuBM9Ihc4Cf596i42RuSddJ2ixp8xNPPDG5T5HDs6DBV8PyJJy64enaeapTeby3xTrR98vENPa8CFgA/BL4O0nv63VIh7TjJnOnlT2rAMbGxnpO9i5CVIjieJJMnYlZIcUY9S8bB0GRoY+3ALvN7Akzex5YD7wB2CdpNkD6/XjKX3oy92So019tDz2fkPTg8XBNPdStftSprXqjiKgfAS6SdGqa2H0pWdCRDcDilGcxcHvaLj2ZuyytG16nm+6hIXkQSlPxcG091LEi1K3teqDv0IeZbZJ0G7CVbHL2P5ENWcwA1kq6lkzmV6X8pSdz9y1kzW+qhwbkQSRNx8N8a29zrHuRb9cxLNKbQgtezOzTwKfbkg+S9a475S81mfu4QtVczHlC0qNFyHpitLf5EPexjGSsj6kiJD2aeLjmHureZKjj8OYwCVE3GA/CGFXi2geDJEQ9JKru0YQoqqfqe1B1HQwGR4h6CFTdQKoWRHCUqu9F1XUxGAwh6gFTdcOoWgzB8VR9T6quk8HkCVEPkKobRNVCCLoT9yaYDCHqhhAi8E+V96jqTkSTkTRP0j2SxlM46KUp/aq0f1jSWNsxpUJBh6gHRJUNISRdH0LWjeQQcL2ZnQdcBCxJ4Z53AO8Cvp/PPJFQ0CHqARCSDsoQsm4WZrbXzLam7QNkITbmmNm4mT3c4ZDSoaBD1JMkKn4wEULWzUTSfOACYFOPbIVCQecZmWcmDoOqK3z0poNgYujg80z78aNFs8+StDm3vyqFaT72PaUZwDpgmZk91ev0HdJ6hnoOUdeUkHT9qTIuSB3jgVTIfjMb65VB0nQySa8xs/V93q90KOgY+pggMS4dDIIYAqk/KfzzLcC4ma0scEjpUNAh6gkQkg4GSci69lwMXANcImlb+rlC0jsl7QFeD3xT0p2QhYIGWqGgv0WBUNAx9FGSkPRguGLGgwN5nzuePn8g71M1MQxSX8zsXjqPOwP87y7HlAoFHaIuQfQ+Js6gxNzvfessbg+xrAOfhKhrQh1708OSc9Fz1lnaU030qn0TY9QFiSGP4lwx48FKJO21HGWI8eqgEyHqAoSki+FVjF7L1Y2QddBODH30ISTdmzoJsE7DIjFeHeSJHnUPonfRmzpJup06l33YRL33R4i6C1VXVs+96boNJXTD++eIIZCgRYjaId4l3TQ8f6aQdQAh6o7EuPTxeO99ThbPny9kHYSo24iKeTxeBTYMRumzFiXaRPXErI8cVVdIb73pUZVW63N7mhkSs0BGm+hRJ0LSxzKqks7j7RrEEMjoEqIOjsOboKokrsVRQtbVEaKm+groqTcdYjoeT9ek6rpSdVsZVUZe1FVXvKobXh5PQvKGp2vjqc4EU8NIizokfRRPIvJKXKOMqtvNKDKyoo7KdpQQUHG8zLeu+o98tJ+jSJon6R5J45J2Slqa0l8qaaOkH6Xfp+eOWS5pl6SHJV3W7xwjKWoPlazqhgZ+pFNHPFy3quuQh3bkhEPA9WZ2HnARsETS+cANwN1mthC4O+2TXrsaeCVwOfAFSSf2OsHIiToqVzAoPMi6aqI9gZntNbOtafsAMA7MARYBq1O21cA70vYi4FYzO2hmu4FdwIW9ziEzG0LRyyHpAPBw1eUoySxgf9WFKEndyly38kKUeSp4uZnNnMwbSPoW2ecuwinAs7n9VWa2qsv7zge+D7wKeMTMXpJ77UkzO13S/wTuN7O/Tem3AP9oZrd1K4CXlYkPm9lY1YUog6TNUebhUrfyQpR5KpC0ebLvYWaXD6IseSTNANYBy8zsKanb8247Pgi3Z4955IY+giAIBo2k6WSSXmNm61PyPkmz0+uzgcdT+h5gXu7wucBjvd4/RB0EQTAJlHWdbwHGzWxl7qUNwOK0vRi4PZd+taSTJS0AFgIP9DqHl6GPjuM9zokyD5+6lReizFOBt/JeDFwDbJe0LaXdCNwErJV0LfAIcBWAme2UtBZ4kGzGyBIze6HXCVx8mRgEQRB0J4Y+giAInBOiDoIgcE6lopb0Gkn3Sdou6R8knZZ7rdQSy6lE0odTuXZK+mwu3V2ZJX1G0j9L2ibp25LOzr3mrrwtJF2eyrVL0g1Vl6cdSadIekDSD1M9+JOU3nXZsAckvUTSbZIeSkueX++5zJKWStqRrvGylOa2vEPDzCr7AX4A/E7a/gDwmbR9PvBD4GRgAfBj4MQqy5or8+8CdwEnp/0zPJcZOC23/RHgi57Lm8p2YirPy4CTUjnPr7pcbWUUMCNtTwc2kS0f/ixwQ0q/AfizqsvaVu7VwH9M2ycBL/FaZrJFIzuAU8kmPtxFNkPCZXmH+VP10MfLyVbxAGwE3p22Sy+xnEL+CLjJzA4CmFlrbqTLMpvZU7ndF3N0Yr3L8iYuBHaZ2U/M7DngVrLyusEynk6709OP0X3ZcOWk/1jfRDaVDDN7zsx+id8yn0e2gu8ZMzsEfA94J37LOzSqFvUO4Mq0fRVHJ4HPAX6ey7cnpXngXOCNkjZJ+p6k307pbsssaYWknwP/AfhUSnZbXnyX7QiSTkzTsR4HNprZJuBMM9sLWQwI4IwKi9jOy4AngC9L+idJfyPpxfgt8w7gTZJ+Q9KpwBVkjvBa3qExdFFLuiuNMbX/LCIb7lgiaQswE3iudViHt5qyeYR9yjwNOJ3s39z/SjZPUlWWuU95MbNPmNk8YA3wodZhVZW3AJ7LdgQze8HMXku2suxCSa+quEj9mAa8DvhrM7sA+H+kiG4eMbNx4M/I/tv+FtkQ2KFKC1URQ1/wYmZv6ZPlrQCSzgXeltJKL7EcJL3KLOmPgPWWDZA9IOkwWXCXyspc4Bq3+CrwTeDTVHyN++C5bMdhZr+U9F2ykJX7JM02s71ty4Y9sAfYk3r+ALeRidptmc3sFtJQjaT/QfYZ3JZ3WFQ96+OM9PsE4JPAF9NLpZdYTiF/D1wCR/64nEQWecxlmSUtzO1eCbSCGLssb+IHwEJJCySdRBa7d0PFZToGSb8p6SVp+0XAW8iubbdlw5VjZr8Afi7p5SnpUrLVcW7LnHPEOcC7gK/huLzDouol5O+VtCRtrwe+DBNbYjmFfAn4kqQdZEM1i1Pv2muZb0oN8zDwM+CD4Psam9khSR8C7iSbAfIlM9tZcbHamQ2sVhbw/QRgrZl9Q9J9dFg27IgPA2vSH8CfAO8nld9pmddJ+g3gebI6+qSkjkuzm0wsIQ+CIHBO1bM+giAIgj6EqIMgCJwTog6CIHBOiDoIgsA5IeogCALnhKiDIAicE6IOgiBwzv8HJqhaogQLcaYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_temp_section(rcelat)\n",
    "plot_temp_section(diffmodel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inferred_heat_transport( energy_in, lat_deg ):\n",
    "    '''Returns the inferred heat transport (in PW) by integrating the net energy imbalance from pole to pole.'''\n",
    "    from scipy import integrate\n",
    "    from climlab import constants as const\n",
    "    lat_rad = np.deg2rad( lat_deg )\n",
    "    return ( 1E-15 * 2 * np.math.pi * const.a**2 * integrate.cumtrapz( np.cos(lat_rad)*energy_in,\n",
    "            x=lat_rad, initial=0. ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x16035fee0>]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAApFUlEQVR4nO3deVyVZf7/8deHHVlENhUQUXBDU1Tcy2W0sk1bpj2zZr7ZNq3TtE6TbdN32mZq6ltZZntpu2lTVtpiLomGC7ihgIgLKIIoO1y/PzjOj0wT5cB1ls/z8TiP4ED3eXd3ztvb+77u6xJjDEoppdyXj+0ASimlWkaLXCml3JwWuVJKuTktcqWUcnNa5Eop5eb8bLxodHS0SUpKsvHSSinltlauXLnHGBNz+PNWijwpKYmMjAwbL62UUm5LRPKP9LyeWlFKKTenRa6UUm5Oi1wppdycFrlSSrk5LXKllHJzWuRKKeXmtMiVUsrNWRlHrpRyTZU19ewoq2RHaeOj5GAtp6Z2JCU21HY09Ru0yJXyMsXl1SzP3eso6yoKHaW9s6yKkoM1v/r9x7/cwOmpnbhhXDL9EyLaPrA6Ji1ypbzIZ6t3cN/Ha9lfVQdAWJAf8RHBxEUEMzAxgs7tg//7fVxEEAG+Pry5LJ/XluTxRdYuTk6J5oaxyYxIjkJELP/XqEPExgpB6enpRm/RV6rtlFfV8sDcLD5aVcjAxAgeOKcv3WNCCA/yb/a///bybcxcnEtxeTUDukRww9hkTu3TER8fLfS2IiIrjTHpv3pei1wpz7Yyv4RbZ2dSuK+Sm37Xg5t+l4Kf74mNc6iqreeDldt56fstFJRU0iM2lOvGJDMpLQ7/E9ymaj4tcqW8TF19A88uzOG5hZuJ7xDMvy5OY3DXSKdte/7anbzw7RY27ConPiKYv5zei3MHxjtl++rIjlbkeo5cKQ+Ut+cgt87OJLOglPMHxfPgpL6ENfM0SnP4+fowOS2eSQPiWLihiGe/2cytszPZX1XLlSOSnPY6qnm0yJXyIMYY3l+5nelzs/DzEf596UDOGRDXaq8nIozv05FTesRww9ur+NunWfiIcMXwrq32murXtMiV8hD7DtZw78dr+c+6XQzrFsk/L04jLiK4TV47wM+H5y8fyPVvreKvn6zDR4TLhiW2yWsrLXKlPEJBSQUXvriUPQequWtib6aN7o5vG48mCfTz5YUrBnHtmyu59+O1+PrAxUO0zNuCXmZWys1V1tQz7c2VVNTU8dENI7l+bHKbl/ghgX6+vHjFYEb3jOHuj9byfkaBlRzeRotcKTdmjOHuj9awYdd+nrl0oEvceRnk78uMKYM5OSWaOz9cw4crt9uO5PFaXOQiEiQiP4nIahHJEpEHnRFMKXVss37M49PMHdw+oSfjesXajvNfQf6+vHxlOiOTo7jjg9V88nOh7UgezRlH5NXA74wxA4A0YKKIDHfCdpVSv2HZ1r08+vl6Tk3tyI3jUmzH+ZUgf19euXIIw7tFcfucTD7N1DJvLS0uctPogONbf8ej7e8yUsqL7Cyr5E/vrKJrVDuevmiAy94mHxzgy8yr0hmSFMltszOZt2aH7UgeySnnyEXEV0QygSLgK2PMcmdsVyn1a9V19Vz31ioqa+qZMWWwU2/0aQ3tAvx49aohpHeN5Jb3Mvl87U7bkTyOU4rcGFNvjEkDEoChItLv8N8RkWkikiEiGcXFxc54WaW80vS5WawuKOWpi9JIiQ2zHadZQgL9ePXqIQzsEsHN7/7Mkpw9tiN5FKeOWjHGlALfAhOP8LMZxph0Y0x6TEyMM19WKa/x7k/bePenAm4Ym8zEfp1sxzkuoYF+zLp6CIlR7bh9zmrKKmptR/IYzhi1EiMiEY6vg4EJwIaWblcp9Us/b9vHA59mMbpnDH8+rZftOCckLMiff12cxp4D1fz103W243gMZxyRdwYWicgaYAWN58jnOWG7SimH4vJqrn9rFR3bB/LsJWnWbvhxhv4JEdwyvgefrd6hI1mcpMW36Btj1gADnZBFKXUEtfUN3PjOKkora/jo+lFEtAuwHanFrh+bzKKNRfz1k3WkJ0US30ZzwngqvbNTKRf398/X81NuCf97fn9S48Jtx3EKP18f/nlxGg0NhjvmrKahQUcst4QWuVIu7NPMQmb9mMfVo5I8btGGrlEhPHBOX5Zu3cvMxbm247g1LXKlXFRZRS3T52YxKDGCe8/sYztOq7gwPYHTUjvyxJcbWb9zv+04bkuLXCkX9cw3mymtrOXhc/t57HqYIsJj559EeLA/t76XSVVtve1Ibskz3x1KubmcogO8sTSPS4Z0oW9ce9txWlVUaCBP/L4/G3eX8+SXG23HcUta5Eq5oEfmZxPs7+u248WP17jesUwZ3pVXFufqXZ8nQItcKRezaGMR324s5ubxPYgODbQdp83ce2YfuseE8Of39a7P46VFrpQLqa1v4OF52XSLDmHqyCTbcdpUcIAv/7o4jeLyau7Xuz6Pixa5Ui7kjaX5bC0+yH1n9iHAz/s+nofu+pyrd30eF+97pyjlokoO1vDM15s4pUc04/u4zmo/be36sckMSozgr5+so7C00nYct6BFrpSLePqrjRysqedvZ6ci4r5zqbRU07s+7/pgDcboXZ/HokWulAtYv3M/7yzfxhXDEunR0T3mGG9NXaNCuPuM3izO2cO8NboQxbFokStlmTGGh+dlEx7sz22n9rQdx2VcNqwr/eLDeWR+Ngeq62zHcWla5EpZtiB7N0u27OW2CT09YmZDZ/H1ER6a3I/d+6v59zebbcdxaVrkSllUXVfPo/PX0yM2lMuHJdqO43IGJXbg4vQuzFycy+bd5bbjuCwtcqUsenVxHttKKrj/7FT8PHQ+lZa6c2IvQgL9+NunWXrh8yj0naOUJUXlVTy3cDMT+sQyuqeuY3s0UaGB/OX0XizdupfP9MLnEWmRK2XJE19spKa+gfvOSrUdxeVdOjSRfvHhPKoXPo9Ii1wpC9ZuL+ODVdu5amQS3aJDbMdxeb4+wsOOC5/P6oXPX9EiV6qNGWN48LMsItsFcNP4HrbjuI2BiR24ZEgXXl2cyya98PkLLS5yEekiIotEZL2IZInILc4IppSnWrihiIz8ffz5tF6EB/nbjuNW7pzY23Hhc51e+GzCGUfkdcCfjTF9gOHAjSKiJ/2UOgJjDM8vyiE+IpgL0xNsx3E7kSEB/OX0XizbWsLc1Ttsx3EZLS5yY8xOY8wqx9flwHrAs1aJVcpJlueWsGpbKdeO6e6xy7e1tkuHJnJSfHsenb9eL3w6OPWdJCJJwEBg+RF+Nk1EMkQko7i42Jkvq5TbeH5RDtGhAVyU3sV2FLfl6yM8fG4/ig9U88zXm2zHcQlOK3IRCQU+BG41xvxqOWxjzAxjTLoxJj0mRsfMKu+zZnspP2zewx9P7k6Qv6/tOG4trUsEF6d34dUf89i4Sy98OqXIRcSfxhJ/2xjzkTO2qZSn+b9FWwgL8uOK4XorvjPcObE3oXrhE3DOqBUBZgLrjTFPtzySUp4np6icL7N3MXVEEmE6UsUpIkMCuHNiL5bn6oVPZxyRjwKmAL8TkUzH40wnbFcpj/HCt1sJ9PPh6lFJtqN4lEuGJNI/ofHCZ3mV9y7Y7IxRK4uNMWKM6W+MSXM8PndGOKU8QUFJBZ9kFnLp0ESiQgNtx/Eoh6a6bbzw6b13fOr4J6Va2cs/bMVH4JpTutuO4pEOXfictSTPa+/41CJXqhUVl1cze0UB5w2MJy4i2HYcj+XtFz61yJVqRa/+mEtNfQPXjUm2HcWjRYYEcIfjjk9vnOpWi1ypVlJWWcubS/M586TOdI8JtR3H413mxVPdapEr1UreXJrHgeo6bhirR+NtwddHeHCSd67xqUWuVCuorKnn1R/zGNcrhr5x7W3H8RqDu3bgwsEJzFycS06R91z41CJXqhW8t2IbJQdruGFciu0oXueuM3rTLsCXB+Z6zxqfWuRKOVlNXQMzvt/K0KRIhiRF2o7jdaJDA/nzab34MWcvn6/dZTtOm9AiV8rJPvm5kJ1lVdwwTs+N23L5sET6dA7nkfnZHPSCC59a5Eo5UX2D4YXvttA3LpwxPXWWT1v8fH14eHJfdpZV8dyiHNtxWp0WuVJO9J91O8ndc5Abx6XQOJ+csiU9KZLzB8Xzyg9b2VJ8wHacVqVFrpSTGGP4v0Vb6B4dwul9O9mOo4B7zuhDkJ8v0z38wqcWuVJO8u2mYrJ37ue6scn4+ujRuCuICQvk9tN68sPmPXyZ5bkXPrXIlXKSl77bQuf2QZybpkvWupIpw7vSu1MYD89bT2VNve04rUKLXCknWL9zP8u2ljB1ZBIBfvqxciV+vj48NLkfhaWVPO+hFz71HaeUE7yxNI8gfx8uGaKLKruiod0iOW9gPDO+30runoO24zidFrlSLbTvYA0f/1zIeQPjiWgXYDuOOop7zuhNgJ+PR051q0WuVAvNziigqraBqSOTbEdRvyE2PIi7Jvbih817eGtZvu04TqVFrlQL1NU38ObSfEZ0j6J3p3DbcdQxXDG8K2N6xvDo5+s9amy5U4pcRF4VkSIRWeeM7SnlLr5ev5vC0kqu0kWV3YKI8MTv+xPk78ttszOprW+wHckpnHVE/how0UnbUsptzPoxj/iIYCb06Wg7imqm2PAgHjvvJNZsL/OYecudUuTGmO+BEmdsSyl3kb1jP8tzS5g6sqveAORmzjipMxcMSuC5RTmszN9nO06Ltdk5chGZJiIZIpJRXFzcVi+rVKt5fUnjkMOL0nXIoTuaPimVuIhgbp+T6fYzJLZZkRtjZhhj0o0x6TExOiuccm8lB2v4JLOQ8wYm6JBDNxUW5M/TF6WxraSCR+Zn247TIjpqRakT8N6KbVTXNXCVDjl0a0O7RXLt6GTe/amAr7J3245zwrTIlTpOdfUNvLU0n5HJUfTqFGY7jmqh20/tSWrncO7+cA3F5dW245wQZw0/fBdYCvQSke0i8kdnbFcpV/RV9m52lFXp0biHCPDz4V+XpFFeXcfdH65xy7s+nTVq5VJjTGdjjL8xJsEYM9MZ21XKFc1akkdCh2DG65BDj9GzYxh3T+zNNxuKeG9Fge04x01PrSh1HLJ2lPFTbglTRyTpkEMPc9XIJEalRPHwvGzy3GxiLS1ypY7D60vyCPb31SGHHsjHR3jywgH4+Qi3zs6kzo3u+tQiV6qZGocc7uD8QfG0b+dvO45qBZ3bB/PoeSeRWVDK84u22I7TbFrkSjXTuz9to0aHHHq8cwbEcW5aHM8u3MyPOXtsx2kWLXKlmqG2voG3luUzKiWKHh11yKGne3ByP1JiQvmf1zNYvnWv7TjHpEWuVDMsyNrNzrIqrhrZzXYU1QbaB/vz9jXDiO8QzNWvrSAjz7WnktIiV6oZXluSS5fIYH7XO9Z2FNVGokMDeed/htEpPIirZq1g1TbXnVxLi1ypY1hXWMaKvH065NALxYYH8c41w4kKDWDqzJ9YXVBqO9IRaZErdQyvOYYcXqhDDr1Sp/ZBvHvNcCJC/JkycznrCstsR/oVLXKlfsPeA9XMXb2DCwbH0z5Yhxx6q7iIYN69ZjhhQf5cMXM52Tv22470C1rkSv2G2RkF1NQ1cOWIJNtRlGUJHdrx7jXDCfb35YqZy9m4q9x2pP/SIlfqKOobDG8v28bw7pH01CGHCkiMaixzf1/h8leWkVPkGmWuRa7UUXy7sYjC0ko9Gle/kBQdwjvXDEdEuPTl5WwpPmA7kha5UkfzxtJ8YsMCOTVVZzlUv5QcE8q71wzDGMNlLy+zPsmWFrlSR5C/9yDfbSrm0qGJ+Pvqx0T9WkpsGG//z3Bq6w1nPfsDTy/YyP6qWitZ9B2q1BG8tSwfXx/hsmGJtqMoF9arUxgf3zCSsb1ieXZhDqMfX8RL322hsqa+TXNokSt1mKraeuZkbOf0vh3pGB5kO45ycV2jQnj+8kHMu+lk0rpE8Nh/NjDmiUW8uSyfmrq2mQpXi1ypw8xdvYOyylqmDE+yHUW5kX7x7Xnt6qHMuXYEXaPacf8n65jw9Hd8/PN26htad/k4LXKlDvPWsnx6xIYyvHuk7SjKDQ3tFsmca0cw6+ohhAb6cdvs1ZzxzPd8mbWr1dYDddbiyxNFZKOI5IjI3c7YplI2rC4oZc32MqaM6IqIzquiToyIMK5XLPNuOpnnLhtIXb3h2jdXcu7/LWmVW/xbXOQi4gs8D5wBpAKXikhqS7erlA1vLM2nXYAv5w2Mtx1FeQAfH+Hs/nEsuG00j1/Qn70Hqgnwc/6JED8nbGMokGOM2QogIu8Bk4FsJ2z7FwpKKqiuqyclVu+yU86372ANn63ZwYWDEwgL0nlVlPP4+fpw0ZAuXDA4oVVm0HTGHw3xQEGT77c7nnO6pxZs5LR/fs+dH6xmR2lla7yE8mJzHPOqTBnR1XYU5aFaaxpkZxyRHynZr87oi8g0YBpAYuKJjc29/+xUIkMCeWtZPp9k7mDqiK7cMDaFDiEBJ7Q9pQ5paDC8tTyfoUmR9O4UbjuOUsfFGUfk24GmEzUnADsO/yVjzAxjTLoxJj0mJuaEXigqNJC/nZPKwjvGcE7/OF5ZnMvoxxfx3MLNVNTUnVh6pYDvNhVTUFKpR+PKLTmjyFcAPUSkm4gEAJcAc52w3aNK6NCOpy4awBe3jGZ4chRPLtjE6Me/5c2leW02AF95ljeX5RMdGsjpfTvZjqLUcWtxkRtj6oA/AV8C64E5xpislm63OXp1CuPlK9P58PoRdI8O4f5Ps5jw9Hd8mllIQysPwFeeo6CkgkUbi7hsaJdWGVGgVGtzyrvWGPO5MaanMSbZGPOoM7Z5PAZ3jWT2tcOZddUQ2gX4cst7mZz178Vs2OVaq3go1/TW8nx8RLhU51VRbspjDj9EhHG9Y/n85lN45pI09h6o5vKXl7vMxO/KNVXV1jNnRQET+sTSuX2w7ThKnRCPKfJDfHyEyWnxvDetceL3y15eTq7luYKV65q/Zif7Kmp18Qjl1jyuyA/pHhPKO9cMo66hceL3gpIK25GUC3pzWT7dY0IYmRxlO4pSJ8xjixygZ8cw3vrjMCpq6rnslWXsLNObiNT/t3Z7GZkFpUwZrvOqKPfm0UUOkBoXzht/GErpwVoue3k5RfurbEdSLuLNZXkE+/ty/qAE21GUahGPL3KAAV0ieO0PQ9i9v4rLX1nO3gPVtiMpy8oqavk0cwfnDoynfbDOq6Lcm1cUOTQOUXz1qiEU7Kvgipk/UVpRYzuSsuj9lQVU1zUwZbjeyancn9cUOcDw7lHMmJLOlqIDXPnqT9YWSlV21TcY3liaT3rXDqTG6bwqyv15VZEDjO4ZwwtXDCJ7x36unrWCg9U6R4u3WbShiG0lFVw1Ksl2FKWcwuuKHGB8n478+9KBZBaU8sfXV7T5itfKrteW5NEpPEjnVVEewyuLHOCMkzrz9EUDWJ5bwvVvr9S5WbzE5t3lLM7Zw5QRXfH39dq3v/IwXv1OnpwWz0OT+/HtxmJmLcmzHUe1gdeW5BHg58MlQ7oc+5eVchNeXeQAVwxLZEKfWP7xxQY279Z5WTxZWUUtH60q5Ny0OKJCA23HUcppvL7IRYTHzu9PaKAft89ZTW29zmfuqeZkFFBZW8/UkUm2oyjlVF5f5AAxYYH8/byTWFtYxr8X5tiOo1pBfYPh9aV5DO0WSd+49rbjKOVUWuQOE/t14vxB8Ty/KIfMglLbcZSTfbN+N9v3VXK1Ho0rD6RF3sT0SX3pGBbI7bMzdUiih3ltSR5x7YM4NbWj7ShKOZ0WeRPhQf48eeEAtu45yD++2GA7jnKSjbvKWbJlL1NGJOGnQw6VB9J39WFGpkRz9agkXluSxw+bi23HUU7w2pJcgvx1yKHyXFrkR3DXxN4kx4Twl/fXUFah87G4s30Ha/j450LOGxhPh5AA23GUahUtKnIRuVBEskSkQUTSnRXKtiB/X/55cRrFB6qZ/lmW7TiqBWZnFFBV26BDDpVHa+kR+TrgfOB7J2RxKf0TIrjpdyl8/HMhn6/daTuOOgF19Q28uTSfEd2j6N1JZzlUnqtFRW6MWW+M2eisMK7mxnEp9E9oz30fr9WVhdzQ1+t3U1haqbMcKo/XZufIRWSaiGSISEZxsXtcRPT39eHpi9KoqKnn7o/WYoxOrOVOZv2YR3xEMBP66JBD5dmOWeQi8rWIrDvCY/LxvJAxZoYxJt0Ykx4TE3PiidtYSmwod5/Rm4UbinhvRYHtOKqZsnfsZ3luCVNHdsXXRxdWVp7N71i/YIyZ0BZBXNnUEUl8lb2bh+dlMyo5msSodrYjqWN4fUnjwsoXpyfajqJUq9Phh83g4yM8ceEAfER0FIsbKDlYwyeZhZw/KJ727XRhZeX5Wjr88DwR2Q6MAOaLyJfOieV64iOCuXl8Cgs3FLFww27bcdRveG/FNqrrGrhKhxwqL9HSUSsfG2MSjDGBxpiOxpjTnRXMFV01shvdY0J46LNsqut0LhZXVOsYcnhySjQ9OobZjqNUm9BTK8chwM+H6ef0JW9vBa/8kGs7jjqCBVm72VlWpUfjyqtokR+n0T1jOC21I88tzGFnWaXtOOowry3JJTGyHeN6x9qOolSb0SI/AfefnUqDMfz9c50h0ZWsLihlRd4+rhyhQw6Vd9EiPwFdIttx7ZhkPlu9g2Vb99qOoxyeX5RD+2B/LhmqQw6Vd9EiP0HXj0kmPiKY6XOzqNN1Pq3btLucBdm7uWpkEqGBx7w9QimPokV+goIDfPnrWX3YsKuct5dvsx3H673w7RbaBfjqRU7llbTIW2Biv06MSoniqQUb2Xug2nYcr7VtbwVzV+/g8mGJOue48kpa5C0gIkw/py8VNfU8ucBjJ4F0eS9+vwVfEa45pbvtKEpZoUXeQj06hjF1ZBLvrShgzfZS23G8zq6yKj7I2M6F6QnEhgfZjqOUFVrkTnDLhB5EhQTywNwsGhp0qtu29MoPW6k3huvGJNuOopQ1WuROEB7kz10Te/HztlI++rnQdhyvse9gDW8v38bkAXF0idQZKZX30iJ3kgsGJTAwMYL//c8G9lfpgs1tYdaPuVTW1nP9WD0aV95Ni9xJfHyEByf1Ze/Bap79erPtOB6vvKqW15bkcXrfjjo5lvJ6WuRO1D8hgovTu/Dakjxyisptx/Foby/fxv6qOm4cl2I7ilLWaZE72V9O70W7AF+mz83WNT5bSVVtPa/8kMspPaLpnxBhO45S1mmRO1lUaCC3n9qTxTl7+DJrl+04HmlORgF7DlTr0bhSDlrkreCK4V3p3SmMh+etp7JGF6Bwptr6Bl76bivpXTswrFuk7ThKuQQt8lbg5+vD9El9KSyt5IXvttiO41E++bmQwtJKbhyXgohOVasUaJG3muHdo5g0II4Xv9vCtr0VtuN4hPoGwwvfbSG1czhje8XYjqOUy2jp4stPiMgGEVkjIh+LSISTcnmEe8/sg5+P8PD8bNtRPMIX63axtfigHo0rdZiWHpF/BfQzxvQHNgH3tDyS5+jUPoibfteDr7J38+3GIttx3JoxhucX5dA9OoSJ/TrZjqOUS2lRkRtjFhhj6hzfLgMSWh7Js/zh5CS6R4fw4GfZVNfphc8T9e3GYrJ37ue6scm6jJtSh3HmOfI/AP852g9FZJqIZIhIRnFxsRNf1rUF+vnyt3NSyd1zkFcX59mO45aMMTy3KIf4iGDOGxhvO45SLueYRS4iX4vIuiM8Jjf5nfuAOuDto23HGDPDGJNujEmPifGuC1Vje8VyampH/r1wM7vKqmzHcTs/5ZawMn8f00Z3x99Xr88rdbhjfiqMMROMMf2O8PgUQESmAmcDlxu9lfGo7j8rlboGw98/X287ilsxxvDUgk1EhwZw8ZAutuMo5ZJaOmplInAXMMkYo2PsfkNiVDuuG5PM3NU7WLZ1r+04bmP+2p38lFfC7af2Isjf13YcpVxSS/+e+hwQBnwlIpki8qITMnms68ckEx8RzPS5WdTVN9iO4/Iqa+r5+/z1pHYO16NxpX5DS0etpBhjuhhj0hyP65wVzBMFB/hy/9l92LCrnLeW5duO4/Je+n4LO8qqeOCcVB2potRv0CtHbez0vp04OSWap7/axJ4D1bbjuKzC0kpe/G4LZ/XvzLDuUbbjKOXStMjbmIgwfVIqFTX1PPHFRttxXNZjn6/HmMa7Y5VSv02L3IKU2DCuHpXEnJUFZBaU2o7jcpZv3cu8NTu5znFNQSn127TILbl5fA+iQwN54NN1NDToqM1D6hsM0z/LJq59ENeN0bU4lWoOLXJLwoL8uffM3qzeXsbbP22zHcdlzF5RwPqd+7nnzD4EB+hwQ6WaQ4vconPT4jmlRzR/n7+evD0HbcexrqyilicXbGRoUiRn9+9sO45SbkOL3CIR4fHf98ffV7h9Tib1Xn6K5ZlvNrOvooa/nZOq09QqdRy0yC3r3D6Yh8/tx6ptpbzoxasJ5RSV88bSPC4Zkki/+Pa24yjlVrTIXcCkAXGcdVJn/vX1JrJ2lNmO0+aMMTz4WTbBAb7ccVpP23GUcjta5C5ARHjk3H50aBfA7bNXU1XrXfOWf7O+iB827+HWCT2JCg20HUcpt6NF7iI6hATwj9/3Z+Pucp7+apPtOG2muq6eR+ZnkxwTwpUjutqOo5Rb0iJ3IeN6xXLZsERe/mGr18yQOOvHPPL2VvC3c/rqXONKnSD95LiY+87sQ2JkO+54fzXlVbW247SqovIq/v3NZib0iWVMT+9abEQpZ9IidzEhgX48fdEAdpRW8vC8bNtxWtXjX2ykpr6B+85KtR1FKbemRe6CBneN5LoxyczJ2M5X2bttx2kVizYW8cHK7fxhVDe6RYfYjqOUW9Mid1G3TuhJn87h3PPRGvZ62HS3uXsOcvO7P9Onczi3TOhhO45Sbk+L3EUF+Pnwr4vT2F9Zxz0frcVTlkM9UF3HtDcy8PMRZkwZTLsAP9uRlHJ7WuQurFenMO44vScLsnfz4apC23FarKHB8Oc5mWzdc5DnLhtEl8h2tiMp5RG0yF3cH0/uztBukUyfm8X2fe69vvXzi3L4Mms395zRm1Ep0bbjKOUxWlTkIvKwiKxxLLy8QETinBVMNfL1EZ66cADGGG6bnem2d31+nb2bp7/exHkD4/njyd1sx1HKo7T0iPwJY0x/Y0waMA/4W8sjqcN1iWzHYxf0JyN/H9e9tZLqOvcq8y3FB7htdiZ948J57PyTdGZDpZysRUVujNnf5NsQwDOuyLmgSQPi+Pt5J/HtxmL+9M7P1NY32I7ULPurarnmjQwC/Hx4aUo6Qf66WIRSztbic+Qi8qiIFACXo0fkrerSoYk8OKkvX2Xv5tb3Mqlz8TJvaDDcPjuT/L0VPH/5IF1/U6lWcswiF5GvRWTdER6TAYwx9xljugBvA3/6je1ME5EMEckoLi523n+Bl5k6Mon7zuzD/LU7+csHa1x6MYpnvtnM1+uLuP+sPgzvHmU7jlIe65iDeI0xE5q5rXeA+cADR9nODGAGQHp6uuu2jxu4ZnR3auobeOLLjQT4+vDY+Sfh4+Na552/zNrFM99s5veDE5g6Msl2HKU8WovuxhCRHsaYzY5vJwEbWh5JNceN41Korq3n2YU5BPj58NDkvi5zEXHz7nJun53JgIT2PHJuP5fJpZSnaultdf8rIr2ABiAfuK7lkVRz3XZqT6rrGnjp+60E+Pnw17P6WC/Nsspapr25kuAAX16cMlgvbirVBlpU5MaYC5wVRB0/EeHuM3pTXdfAzMW5BPj5cOfpvayV+abd5dz6XiYFJRW8O204ndvrxU2l2oJOdOHmRIQHzkmlpr6BF77dQpCfb5tPRNXQYJi5OJcnFmwkLNCPl6YMZkhSZJtmUMqbaZF7ABHhkcn9qKlr4J9fbyLAz4frxya3yWsXlFRwx/urWZ5bwqmpHXns/JOI1nU3lWpTWuQewsdH+McF/amtb+AfX2xgXWEZN41PoXen8FZ5PWMM76/czkOfNS5+8cTv+/P7wQnWz9Er5Y20yD3IoXlZEiPbMevHPOav3cnEvp24aXwKfePaO+119hyo5p6P1vJV9m6GdYvkyQsH6EyGSlkkNua5Tk9PNxkZGW3+ut6ktKKGV3/MY9aPuZRX1TGhT0duHp9C/4SIFm33y6xd3PvRWsqr67jz9F78YVQ3lxvDrpSnEpGVxpj0Xz2vRe7ZyipreX1JHjMX51JWWcu4XjHcNL4HgxI7HNd29lfV8uDcbD5ctZ2+ceH88+I0enYMa6XUSqkj0SL3cuVVtbyxNJ9XftjKvopaTukRzS3je5CeFIkxhrLKWorKqynaX01RedWvvt5SdIB9FTXcMDaFm8f3IMBPp7JXqq1pkSsADlbX8dayfGZ8v5W9B2uIDQuktKKWmiNMwBXs70tseCCxYYF0ah/MVSOTGNz1+I7klVLOc7Qi14udXiYk0I9rxyRz5Ygk3vlpG1mFZcSEBRITFkhseBCxYY3FHRseRGigvj2Ucgf6SfVSwQG+ulKPUh5CT3QqpZSb0yJXSik3p0WulFJuTotcKaXcnBa5Ukq5OS1ypZRyc1rkSinl5rTIlVLKzVm5RV9Eimlc49NVRAN7bIc4Dpq39bhTVtC8rc3V8nY1xsQc/qSVInc1IpJxpPkLXJXmbT3ulBU0b2tzl7x6akUppdycFrlSSrk5LfJGM2wHOE6at/W4U1bQvK3NLfLqOXKllHJzekSulFJuTotcKaXcnFcXuYjMFpFMxyNPRDIdzyeJSGWTn71oOSoiMl1ECptkOrPJz+4RkRwR2Sgip9vMeYiIPCEiG0RkjYh8LCIRjuddbt8eIiITHfswR0Tutp2nKRHpIiKLRGS9iGSJyC2O54/6vrDN8Zla68iV4XguUkS+EpHNjn+6xNqBItKryT7MFJH9InKrK+/fpvQcuYOIPAWUGWMeEpEkYJ4xpp/lWP8lItOBA8aYJw97PhV4FxgKxAFfAz2NMfVtHvKXuU4DFhpj6kTkHwDGmLtccd8CiIgvsAk4FdgOrAAuNcZkWw3mICKdgc7GmFUiEgasBM4FLuII7wtXICJ5QLoxZk+T5x4HSowx/+v4w7KDMeYuWxmPxPFeKASGAVfjovu3Ka8+Ij9ERITGD8S7trOcgMnAe8aYamNMLpBDY6lbZYxZYIypc3y7DEiwmacZhgI5xpitxpga4D0a961LMMbsNMascnxdDqwH4u2mOiGTgdcdX79O4x9GrmY8sMUY40p3n/8mLfJGpwC7jTGbmzzXTUR+FpHvROQUW8EO8yfHqYpXm/yVNB4oaPI723G9D/gfgP80+d4V96077Eeg8fQUMBBY7njqSO8LV2CABSKyUkSmOZ7raIzZCY1/OAGx1tId3SX88qDOVffvf3l8kYvI1yKy7giPpkdbl/LL/3E7gURjzEDgduAdEQm3nPUFIBlIc+R76tC/doRNtcn5subsWxG5D6gD3nY8ZWXfNoO1/Xg8RCQU+BC41Rizn6O/L1zBKGPMIOAM4EYRGW070LGISAAwCXjf8ZQr79//8rMdoLUZYyb81s9FxA84Hxjc5N+pBqodX68UkS1ATyCjFaMeM+shIvIyMM/x7XagS5MfJwA7nBztiJqxb6cCZwPjjeNijK192wzW9mNziYg/jSX+tjHmIwBjzO4mP2/6vrDOGLPD8c8iEfmYxtNXu0WkszFmp+O8f5HVkL92BrDq0H515f3blMcfkTfDBGCDMWb7oSdEJMZxwQMR6Q70ALZayncoU+cm354HrHN8PRe4REQCRaQbjVl/aut8hxORicBdwCRjTEWT511u3zqsAHqISDfHUdklNO5bl+C4jjMTWG+MebrJ80d7X1glIiGOi7KISAhwGo3Z5gJTHb82FfjUTsKj+sXfzl11/x7O44/Im+Hw82EAo4GHRKQOqAeuM8aUtHmyX3pcRNJo/Ot+HnAtgDEmS0TmANk0nsK40faIFYfngEDgq8YOYpkx5jpcc9/iGF3zJ+BLwBd41RiTZTlWU6OAKcBacQyTBe4FLj3S+8IFdAQ+dvy/9wPeMcZ8ISIrgDki8kdgG3ChxYy/ICLtaBy11HQfHvFz52p0+KFSSrk5PbWilFJuTotcKaXcnBa5Ukq5OS1ypZRyc1rkSinl5rTIlVLKzWmRK6WUm/t/awMJnuMkUosAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  Plot the northward heat transport in this model\n",
    "Rtoa = np.squeeze(diffmodel.timeave['ASR'] - diffmodel.timeave['OLR'])\n",
    "plt.plot(diffmodel.lat, inferred_heat_transport(Rtoa, diffmodel.lat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## If you want explicit surface fluxes...\n",
    "\n",
    "All the models above use a convective adjustment that simultaneously adjustments `Tatm` and `Ts` to the prescribed lapse rate.\n",
    "\n",
    "We can instead limit the convective adjustment to just the atmosphere. To do this, we just have to change the `state` variable dictionary in the convective adjustment process.\n",
    "\n",
    "Then we can invoke process models for **sensible and latent heat fluxes** that use simple bulk formulae. Tunable parameters for these include drag coefficient and surface wind speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.process.time_dependent_process.TimeDependentProcess'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (30, 1) \n",
      "  Tatm: (30, 20) \n",
      "The subprocess tree: \n",
      "Seasonal RCE with surface fluxes and heat transport: <class 'climlab.process.time_dependent_process.TimeDependentProcess'>\n",
      "   Insolation: <class 'climlab.radiation.insolation.DailyInsolation'>\n",
      "   Radiation: <class 'climlab.radiation.cam3.cam3.CAM3'>\n",
      "   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   Meridional Diffusion: <class 'climlab.dynamics.meridional_advection_diffusion.MeridionalDiffusion'>\n",
      "   ConvectiveAdjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "diffmodel2 = climlab.process_like(diffmodel)\n",
    "diffmodel2.name = \"Seasonal RCE with surface fluxes and heat transport\"\n",
    "\n",
    "#  Hard convective adjustment -- ATMOSPHERE ONLY\n",
    "convadj2 = climlab.convection.ConvectiveAdjustment(state={'Tatm':diffmodel2.Tatm}, adj_lapse_rate=6.5)\n",
    "diffmodel2.add_subprocess('ConvectiveAdjustment', convadj2)\n",
    "\n",
    "print(diffmodel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.process.time_dependent_process.TimeDependentProcess'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (30, 1) \n",
      "  Tatm: (30, 20) \n",
      "The subprocess tree: \n",
      "Seasonal RCE with surface fluxes and heat transport: <class 'climlab.process.time_dependent_process.TimeDependentProcess'>\n",
      "   Insolation: <class 'climlab.radiation.insolation.DailyInsolation'>\n",
      "   Radiation: <class 'climlab.radiation.cam3.cam3.CAM3'>\n",
      "   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   Meridional Diffusion: <class 'climlab.dynamics.meridional_advection_diffusion.MeridionalDiffusion'>\n",
      "   ConvectiveAdjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   SHF: <class 'climlab.surface.turbulent.SensibleHeatFlux'>\n",
      "   LHF: <class 'climlab.surface.turbulent.LatentHeatFlux'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#  Now add surface flux processes\n",
    "#  Add surface heat fluxes\n",
    "\n",
    "shf = climlab.surface.SensibleHeatFlux(state=diffmodel2.state, Cd=0.5E-3)\n",
    "lhf = climlab.surface.LatentHeatFlux(state=diffmodel2.state, Cd=0.5E-3)\n",
    "#  set the water vapor input field for LHF process\n",
    "lhf.q = diffmodel2.subprocess['H2O'].q\n",
    "diffmodel2.add_subprocess('SHF', shf)\n",
    "diffmodel2.add_subprocess('LHF', lhf)\n",
    "\n",
    "print(diffmodel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1826 steps, 1826.2110000000002 days, or 5 years.\n",
      "Total elapsed time is 16.996940660197534 years.\n"
     ]
    }
   ],
   "source": [
    "diffmodel2.integrate_years(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 365 steps, 365.2422 days, or 1 years.\n",
      "Total elapsed time is 17.99627753857577 years.\n"
     ]
    }
   ],
   "source": [
    "diffmodel2.integrate_years(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAf9klEQVR4nO3df7AdZZ3n8feHAAEJapgENoSwiUzICFkHnJhipUYZUKGcKcPMllamVio1skONlVGwdFfQrbGs2WyhY7Hr1i4zdUtwUg7KZoSRrLuigZWxKIEYECQhZglmBi5ECAISZyQxyXf/6L7QXO75eU93P93n86q6dc/p2+f0c053f863n/P0bUUEZmaWrqPqboCZmXXnoDYzS5yD2swscQ5qM7PEOajNzBLnoDYzS5yD2sxsFiQdJ2mrpIck7ZD02Xz6SZK2SHo0/z2/8JhrJO2WtEvSxT2X4XHUZmbDkyTghIj4haRjgLuBK4E/AJ6LiGslXQ3Mj4hPSjoL+BqwGjgVuAM4MyIOd1qGK2ozs1mIzC/yu8fkPwGsATbm0zcCl+a31wA3R8SBiNgD7CYL7Y6OHnWjh7FgwYJYunTpy/d3PfrTUpd3+DiN9vnmjvTpGmvO3I4FQTIOH5hTdxOSMOdACc/5UjlH5yuW/4tX3b///vufjYiFs3nO8y84Ll547khf8z7y8K92AC8VJk1ExERxHklzgPuBXwf+R0TcJ+mUiNgLEBF7JZ2cz74YuLfw8Ml8WkdJBPXSpUvZtm3by/cvuORzpS7v+RXHjvw595/R30o3q9uJj43+QHr+roMjf84pd93+yVfdl/SPs33OF547wle/eUpf857zLydfiohV3ebJuy3OkfRG4O8krewy+0yVYtdPuSSCumpTG9UoA3tq43dgW6rKCGgoN6SbJiJekHQXcAnwtKRFeTW9CHgmn20SWFJ42GnAU92ed6z7qMvYwE587KjSdgizYZS5TTqkQdLCvJJG0vHAu4AfA5uBdfls64Db8tubgbWS5kpaBiwHtnZbxlhW1EXzdx0spSukuGO4yraqVVEsOKRftgjYmPdTHwVsiohvSroH2CTpcuBx4P0AEbFD0ibgEeAQsL7biA9wUAPlhfUUh7ZVpaqjOYf0KyLiR8C5M0z/GXBRh8dsADb0uwwHda7ssJ7i0LZRqrqbzQFdDwd1QVVhPWX6Tubgtl7q/P7DIV0fB/U0VYd1Uaed0AE+flL7QtohXS8H9QzqDOuZOMDbK7VAnolDun4O6g5SC+uZDLqTO9jL14TgHYRDOg1JBvXx258E4Jcru55VWbomhPUg2hYiVi6HdDqS3nOnArtO83cd9AZrYyeVbT6FDEhB0kFtZtVLJaTtFckHdSqfqN54bRyktJ2nsu+nIPmgTklKG7HZqHn7TlcjgtqfrGbjxfv8qzUiqFPiqsPayNt12hzUQ/BGbW2S2vbsavq1GhPUXnlmo5daSNvMGhPUqfEGbjZ6Lshm1qigTm0lOqytybz9NkejgtrMRiPFkE6tEEtJ44I6tZWZ4gZvZu3SuKBOkcPamiTF7TW1Aiw1jQxqr1Sz9vD+3Fsjgxq8cs3awPtxfxob1GY2uBS7Pay3Rgd1Sp/G3gHMBpPS/jsbkpZI+q6knZJ2SLoyn36OpHslPShpm6TVhcdcI2m3pF2SLu61jCSv8GJmo5dSMdGWkM4dAj4eEQ9IOhG4X9IW4PPAZyPiW5Lem9+/QNJZwFrgbOBU4A5JZ0bE4U4LaHRFDa1b4WbWMBGxNyIeyG/vB3YCi4EAXp/P9gbgqfz2GuDmiDgQEXuA3cBqunBFPUJtu8aiWRlSKK5+dngeX3n+7X3OvWmBpG2FCRMRMTHTnJKWAucC9wFXAd+W9AWyonhqgYuBewsPm8ynddSKoD5++5O1XwjXLGUpdXs00LMRsarXTJLmAbcAV0XEi5L+E/CxiLhF0geAG4B3AZrh4dHtuRvf9WFmzZFCNV0GSceQhfRNEXFrPnkdMHX7b3mle2MSWFJ4+Gm80i0yo9YEdSobgCsXs5mlso+OmiSRVcs7I+K6wp+eAt6Z374QeDS/vRlYK2mupGXAcmBrt2W0ouvDzNL3y5WL2xrW5wOXAQ9LejCf9ingj4EvSjoaeAm4AiAidkjaBDxCNmJkfbcRH9CioHYftZnVISLuZuZ+Z4Df6vCYDcCGfpfRs+ujy2DukyRtkfRo/nt+4TEDDeY2M7PO+umjnhrM/WbgPGB9PmD7auDOiFgO3JnfZ9pg7kuA6yXNKaPxKfLwPLPOfOQ7nJ5dHxGxF9ib394vaWow9xrggny2jcBdwCcpDOYG9kiaGsx9z6gbP8UrP037zzhS27JPfKw135ObDdZHPW0w9yl5iBMReyWdnM/W12BuSVeQd66ffvrpAzfc0lJnKM+kU3sc4NZEfW+10wdzd5t1hmmvGcwdERMRsSoiVi1cuLDfZryGq+n67D/jyMs/TdHENreN99nB9VVRdxjM/bSkRXk1vQh4Jp8+8GBua442BVzxtbjStpT1M+qj02DuzWRn3pD/vq0wfaDB3MPyJ3N12l6Ftv31pfYlt/fdwfRTUXcazH0tsEnS5cDjwPthuMHclqY2B1cnrrItRf2M+ug2mPuiDo8ZaDC3pWccQ3q6/WcccVhbEhq7FfrQqRxt7wIYlN+P8ngf7l9rTiG32XEYdTf1/rjCtjo0Mqj9STw6DujBOLCtDt7axpQP6WfH751VqXFB7Wp69hwyo+EPu9nz/tyfRgW1V+rsOVhGz++pla1RQW3Dc/VXLr+/w3MB1ltjgtorc3gOkOr4vbYyNCKoHdLDc3BUz+/54LyPd9fI4XnWm8OiXh7GZ6OU/FbkT9rBOaTT4XXRP+/rnSUd1E1bcSn8hzIHQ3q8TvrXtH2+Kkl2fXhlDceBkC7/g6f+Te3/x29/suaWpMNbTks4pNPndTSYX65c3IiiTdISSd+VtFPSDklXFv72EUm78umfL0y/RtLu/G8X91pGkhV1E9XZ7eEAaI46K+vnVxzL/F0Ha1l2yx0CPh4RD0g6Ebhf0hbgFLKLfb8lIg5MXVdW0lnAWuBs4FTgDklndvu//a6oG84h3TxeZ+0SEXsj4oH89n5gJ9kFvT8MXBsRB/K/TV2ucA1wc0QciIg9wG5gdbdluKI2q0FdlbWr6syLB4/jO4+v6Hf2BZK2Fe5PRMTETDNKWgqcC9wH/AXw25I2AC8Bn4iIH5CF+L2Fh03m0zpyUI9AXd0ersyazV8wNsazEbGq10yS5pFdBPyqiHhR0tHAfOA84G1kly58EzNfMSu6PbeDuqHaFtLzlv28499+secNFbbEbHCSjiEL6Zsi4tZ88iRwa0QEsFXSEWBBPn1J4eGnAU91e34HdQM1OaS7BfJsHtPUMHdV3XySBNwA7IyI6wp/+gZwIXCXpDOBY4Fngc3AVyVdR/Zl4nJga7dlOKhnqepuj6aF9DDBPIrlNCm4qw5r91OP3PnAZcDDkh7Mp30KuBG4UdJ24CCwLq+ud0jaBDxCNmJkfbcRH+CgbpSmhHRV4dxvG5oQ2q6smysi7mbmfmeAD3Z4zAZgQ7/LcFA3ROohnUI4d9KU0HZYWyfeKmYhhf/tUbd5y36edEhP17T2lsXbbrO4om6AFKvppofdVPtTq7BdVdtMvEUMqaqKJLWQbltFmuLrqWqdu6puDge19S21QBulFAPbbIqDegjjVk2PU4il8jpdVVuR+6gTlUJIpxJaVUul/9r91TbFW8GAxqUCGdeQLhqX92Bctukmc1AnqO5qelwCqh91vxd1bwuWBgf1AKqoPOreMesOphTV3UdfxTbhqjptDmoD6g+jJvD7Y3VxUPepzdW0A6h/db1XdR9pWb0c1GPOIT24toa1uz/S5aBORB0Vk0N6eH7vrEoO6j60sdJw0MxeHe+hq+rx1HdQS5oj6YeSvpnfP0nSFkmP5r/nF+a9RtJuSbskXVxGw9uk6mraIT06fi+tCoNU1FeSXQZ9ytXAnRGxHLgzv4+ks4C1wNnAJcD1kuaMprnVa1uF4WAZvarfU1fV46evoJZ0GvC7wJcKk9cAG/PbG4FLC9NvjogDEbEH2A2sHklrW6jKatohXR6/t1amfivq/wr8B6CYKqdExF6A/PfJ+fTFwBOF+Sbzaa8i6QpJ2yRt27dv36DtrkTZlYWHXLVLlWHtqnq89AxqSb8HPBMR9/f5nDNdOyxeMyFiIiJWRcSqhQsX9vnU1WnbhuqKrxp+n60M/VTU5wPvk/QPwM3AhZL+Bnha0iKA/Pcz+fyTwJLC408DnhpZi1vCXR42W66qx0fPoI6IayLitIhYSvYl4f+NiA8Cm4F1+WzrgNvy25uBtZLmSloGLAe2jrzlJWrTBuqQrp7fcxu12YyjvhZ4t6RHgXfn94mIHcAm4BHgdmB9RByebUPbpKpq2oFRn6ree1fV9ZO0RNJ3Je2UtEPSldP+/glJIWlBYdpAQ5gHunBARNwF3JXf/hlwUYf5NgAbBnnuVLRlw3RI12/esp/XfvEBq8Qh4OMR8YCkE4H7JW2JiEckLSErZB+fmnnaEOZTgTskndmtoPWZiRXzSI/xUsUHpqvqekXE3oh4IL+9n+x8k6mRbv+FbMRccUDFwEOYHdQFbdkgXU2b1UPSUuBc4D5J7wOejIiHps3W1xDmIl8zsWUc0umpoguk7OsrPr/iWObvOlja81ft8IE5g6yTBZK2Fe5PRMTE9JkkzQNuAa4i6w75NPCeGZ6vryHMRQ7qXJv/37TVz/3VjfZsRKzqNoOkY8hC+qaIuFXSvwKWAQ9JgmyY8gOSVjPEEGZ3fbSIq2krU1u6BkdNWRLfAOyMiOsAIuLhiDg5IpbmQ5sngbdGxE8ZYgizg5p2VNMO6fSVvY58xFab84HLyE4GfDD/eW+nmYcZwuyuDzPrW9v6qkchIu5m5n7n4jxLp90faAjz2FfUrqatSq6qbRhjH9RmVWv6B6v7qqs31l0frqbT8J7Tdw38mO88vqKElrRD2UP1rHpjG9RtqAqaGtLDBHOv52hacDd9uJ77qqs1tkFt1RlFMA+yjKaFdlM5rKszlkFdVTVdZrdHE6rpKgK623JTD+wyq2p3f7TL2K3JNnR5NEFdIT29DSm0o5smfOB24/2pGmMV1FVuVONaTacYjqm1pypVDdVzWJdvLLs+mizVkE49DFPuDmn6F4tWvrEJan/qlyP1gJ4u5cBuMn+xWK6x6PqoOqTH5eywpoV0UWptL+tIqcpt0cVQecYiqNsipW6P1IJuGG14DalxWJej9UHtDWf02hRwKb2WNlTVVo5WB3UdIV3WTpFKNZ1SsI1KiiNVmszF0ei1Nqi9sYzWOIRZCq8vlQ/k2fL+N1qtDOq6NpK2VtMpBFhV2vpa3f3RbK0Lan+Sj1Zbg6ubcXzNZfC+ODqtCuo2bhh1VtPjHFh1vva6j6BG6fkVx7Zyv6xaa4K67o2hbYeW4xzSU9r2HtS5jda9fzZdK4K6rRtBmyqrpqorrNu47tu6n1ah0UHtw6pytK2StHR4fx1OY4M6pRXepm4Ph/RrtamqTmFbTWnfbYrGBfW4VNF1HPo6pDvzezNabdqHJS2R9F1JOyXtkHRlPv0vJP1Y0o8k/Z2kNxYec42k3ZJ2Sbq41zIaEdRT4dymlZsaB1Fvfo9Gq0X78yHg4xHxZuA8YL2ks4AtwMqIeAvw/4BrAPK/rQXOBi4Brpc0p9sCkg3qpoRzGYeSbfwiqS2qDuu2dn9Macp+3k1E7I2IB/Lb+4GdwOKI+E5EHMpnuxc4Lb+9Brg5Ig5ExB5gN7C62zKS/H/UTV5pTeRK0VJQ3O+b+r+tJS0FzgXum/anDwH/M7+9mCy4p0zm0zpKMqjHWdXVtEN6cO85fVelFx4YxyvATIV2WYE95wCDXPx3gaRthfsTETExfSZJ84BbgKsi4sXC9E+TdY/cNDVphmVEtwY4qGchpUPIYTikh1d1WI9aU65SnsjR9bMRsarbDJKOIQvpmyLi1sL0dcDvARdFxFQYTwJLCg8/DXiq2/Onv6bMEuUPOgOQJOAGYGdEXFeYfgnwSeB9EfHPhYdsBtZKmitpGbAc2NptGa6oE1Jlt4dDplnGsfujQc4HLgMelvRgPu1TwH8D5gJbsizn3oj4k4jYIWkT8AhZl8j6iDjcbQF9BXU+/u9LwEqyvpQPAbvIOseXAv8AfCAins/nvwa4HDgMfDQivt3Xy22QJnd7OKRHp8ldIE3p/khdRNzNzP3O/6fLYzYAG/pdRr9r6YvA7RHxG8Bvkg0/uRq4MyKWA3fm94caI2hmZp31DGpJrwfeQdYHQ0QcjIgXyMYCbsxn2whcmt8eeIygVdft4Wp69Kp6Tz2+fnz1U1G/CdgHfFnSDyV9SdIJwCkRsReyAd/Ayfn8i4EnCo+fcYygpCskbZO0bd++fbN6EVVrcreHWZG35WboJ6iPBt4K/GVEnAv8E3k3Rwd9jRGMiImIWBURqxYuXNhXY212XE2Xx1W1lamfoJ4EJiNi6kybr5MF99OSFgHkv58pzD/QGMFx552vHfxBaGXpGdQR8VPgCUlTX21fRDasZDOwLp+2Drgtvz3wGMEmaeqhokPEOmnqNj1O+h1H/RHgJknHAj8B/ogs5DdJuhx4HHg/wDBjBK1cDunqVDFcz2Oqx09fQR0RDwIznUJ5UYf5BxojOM7c7WFmvXi0+wCaeIjoarp6TXzPm7htjxMHtVkD+UhsvDioa1T2ztbEyq4t/N7bKDmo++RDQ2s7b+PpclC3lCu6+nkd2Kg4qGviPkabLW9D48NB3UKu5NLhdWGj4KDug/vubFx4W0+Tg7oGZR6yuoJLT5nrxN0f48FBbWaWOAd1Dz4UtFFo0pGOt/n0OKgr5m4PGzV3f7Sfg9rMLHEO6i6adAjoajp9TVpHTdr2x0G//4/arDSXzf/+wI/5yvNvL6ElzeX/Ud1urqg7KKOicF/iKy6b//2Xf2bz+KZpUlVt/ZG0RNJ3Je2UtEPSlfn0kyRtkfRo/nt+4THXSNotaZeki3stw0HdAk3a+UcdsE0N7CZw90ffDgEfj4g3A+cB6yWdRXYR8DsjYjlwZ36f/G9rgbOBS4DrJc3ptgB3fVjpqgjS4jLcLWJVioi9wN789n5JO4HFwBrggny2jcBdwCfz6TdHxAFgj6TdwGrgnk7LcFDPoEndHqlX03VUu5fN/37SYV3WdRXdT92/OS8F83cd7Hf2BZK2Fe5PRMTETDNKWgqcC9wHnJKHOBGxV9LJ+WyLgXsLD5vMp3XkoLbS1NklMbXslAO7CfafcYQTHxv7HtJnI2Kma8a+iqR5wC3AVRHxoqSOs84wLbo999ivARu9lPqNU2nHdKkfCdlgJB1DFtI3RcSt+eSnJS3K/74IeCafPgksKTz8NOCpbs/voJ6mSV+gpLizpxiMKbapLB5ZVD1lpfMNwM6IuK7wp83Auvz2OuC2wvS1kuZKWgYsB7Z2W4aDugLjsvOkHIgpty11TSpeanI+cBlwoaQH85/3AtcC75b0KPDu/D4RsQPYBDwC3A6sj4jD3RbgPuqGSq2abkIQpvYlY1lfKlq1IuJuZu53Briow2M2ABv6XYYr6gJXDsNpQkhPSan/3KxfDmqblaaGXirtLuPIqIyuNhcx9XJQl6yMnSaVbo9Uws6s7RzUOVcMg2lDSLfhNdh4cFDbwNoUcCm8llSOkHpxMVMfB3WJ2tzt0SYphPWojcuQ0HHhoMaVwiDaGGpmqXNQW9/aHNJ1v7amHCm5qKmHg7pB6tyZ6w6yKozDa7RmclCXpE19hOMUYG16rW3aBsfd2Ae1D+UsFe7+sE7GPqiboq6duE0VZr/G8TVb2sY6qF0ZdDfOgTXOr93S01dQS/pYfnXd7ZK+Jum4UV5ht23cN2ip8LbYDj2DWtJi4KPAqohYCcwhu4LuyK6wW4cmVdN1dHu4oqznPXA/tc2k366Po4HjJR0NvI7ssjFryK6sS/770vz2y1fYjYg9wNQVds0axx9YloKeQR0RTwJfAB4nuyT6zyPiO0y7wi5QvMLuE4WnmPEKu5KukLRN0rZ9+/bN7lUMyNVAdw4n64f3o+r00/Uxn6xKXgacCpwg6YPdHjLDtNdcYTciJiJiVUSsWrhwYb/tNatc1R9co+7+cD918/XT9fEuYE9E7IuIXwG3Am9nhFfYrVLZVcCod4qq+yxdTdsgXFVXo5+gfhw4T9Lr8qvtXgTsZIRX2DVLnT/ArE49L24bEfdJ+jrwAHAI+CEwAcwDNkm6nCzM35/Pv0PS1BV2D9HHFXar4k//7hxGNoz9ZxzhxMfG+pSM0vX17kbEZyLiNyJiZURclo/o+FlEXBQRy/PfzxXm3xARZ0TEioj4VnnNb7cquz0c0r1V+R41ZZieVWNsPgarqKb9pY2lquxtc5yPViXdKOkZSdsL086RdK+kB/PRbasLfxv4hMCxCWrrzNV0//xe2Qz+muzkvqLPA5+NiHOAP8vvD31CoIPaLFFN6/4Y16o6Ir4HPDd9MvD6/PYbeGXk21AnBPb8MrENmrgBVbWTukIc3GXzv89Xnn973c2wWTjql7/i+O1P9jv7AknbCvcnImKix2OuAr4t6QtkBfHUBrMYuLcw34wnBE43FkFdBfdPW+rmLfs5v9jzhlKX0dIRIM9GxKoBH/Nh4GMRcYukDwA3kJ2T0tcJgdO17h2dronVdFVcTQ/P7531sI7s5ECAv+WV7o2hTghsfVA3UdP6Jq08TdwWXBwBWfi+M799IfBofnuoEwJb3fXhDaYzV4Sz575qA5D0NeACsr7sSeAzwB8DX8z/4+hLwBUw/AmBrQ1qh7SZVSEi/rDDn36rw/wbgA2DLMNdHyPgLxLHV9OOTKraVl0ojVYrg7rJG0kVfZJNC5dx18R+ahut1gV1k0ParE28L45Oq4LaG0ZvrqZHz+9pZ94nR6NVQV0H909b01S9zTqsZ681Qd2GjcF9kc1VdlXtbWO8tSKo2xDSVfAhutXF++jstCKozSx9DuvhNT6ovfItFT5i6c3763AaHdR1r/RRfilTdh+kQ6T5RrmN+EvwZmlsUNcd0mY2HO+7g2tkUHtFD8bVdHX8XvfH+/BgGvVPmbxyzWwcNaKi3n/GkeRC2n181kmZVXWb+qlT26dTlnxQj8PKLPOLRB+KW8pSLMJSlGxQewWajQ/v690lGdReadZ0PpIZnIuzzpIMahsNh4U1kQP7tRzUQ2jSiS7WPm36QrEbh/UrHNQt5Wq6fl4Hs+ewzjiozcwS56A2M0ucg7qFfMidjrLWhb/bSIekGyU9I2n7tOkfkbRL0g5Jny9Mv0bS7vxvF/ezDAf1gPxForVNyl8oNsRfA5cUJ0j6HWAN8JaIOBv4Qj79LGAtcHb+mOslzem1AAe1mdksRMT3gOemTf4wcG1EHMjneSafvga4OSIORMQeYDewutcyFBEjbPJwJO0HmlZeLgCerbsRA2pam5vWXnCbq7AiIk6czRNIup3sdffjOOClwv2JiJiY9nxLgW9GxMr8/oPAbWRV80vAJyLiB5L+O3BvRPxNPt8NwLci4uvdGpDKf8/bFRGr6m7EICRtc5vL1bT2gttcBUnbZvscEXFJ77lm5WhgPnAe8DZgk6Q3AZqpOb2ezF0fZmajNwncGpmtwBGyCn4SWFKY7zTgqV5P5qA2Mxu9bwAXAkg6EziWrHtpM7BW0lxJy4DlwNZeT5ZK18dE71mS4zaXr2ntBbe5Ckm1V9LXgAuABZImgc8ANwI35kP2DgLrIvtCcIekTcAjwCFgfUQc7rmMFL5MNDOzztz1YWaWOAe1mVniag1qSb8p6R5JD0v6X5JeX/jbwKdZVmWUp4aWTdKfS/qRpAclfUfSqYW/JdfeKZIuydu1W9LVdbdnOknHSdoq6aF8O/hsPv0kSVskPZr/nl93W4skvVHS1yX9WNJOSf865TZLulLS9vw9viqflmx7SxMRtf0APwDemd/+EPDn+e2zgIeAucAy4DFgTp1tLbT5d4A7gLn5/ZNTbjPw+sLtjwJ/lXJ787bNydvzJrJvyx8Czqq7XdPaKGBefvsY4D6yMbOfB67Op18NfK7utk5r90bg3+W3jwXemGqbgZXAduB1ZAMf7iAbJZFke8v8qbvrYwXwvfz2FuDf5LeHOs2yIiM9NbRsEfFi4e4JvDK4Psn25lYDuyPiJxFxELiZrL3JiMwv8rvH5D9B1s6N+fSNwKXVt25m+RHrO4AbACLiYES8QLptfjPZWXz/HBGHgL8Hfp9021uauoN6O/C+/Pb7eWUg+GLgicJ8k/m0FJwJ/Lak+yT9vaS35dOTbbOkDZKeAP4t8Gf55GTbS9pte5mkOfmpws8AWyLiPuCUiNgLkP8+ucYmTvcmYB/wZUk/lPQlSSeQbpu3A++Q9GuSXge8lywjUm1vaUoPakl35H1M03/WkHV3rJd0P3Ai2XhDGPI0y4raXDw19N+TnRqqOtvco71ExKcjYglwE/CnUw+rq719SLltL4uIwxFxDtnZZaslray5Sb0cDbwV+MuIOBf4J7KugyRFxE7gc2RH27eTdYEdqrVRNSn9hJeIeFePWd4DL5+987v5tKFOsxyVbm2W9GHyU0OBrZJmdWroKPTxHk/5KvC/yQbk1/oe95By214jIl6QdBfZP+B5WtKiiNgraRFZtZ2KSWAyr/wBvk4W1Mm2OSJuIO+qkfSfyV5Dsu0tS92jPk7Ofx8F/Efgr/I/DXWaZUW+wQhPDS2bpOWFu+8DfpzfTrK9uR8AyyUtk3Qs2f/v3Vxzm15F0kJJb8xvHw+8i+y93Qysy2dbR/Yf1JIQET8FnpC0Ip90EdkZcsm2uZARpwN/AHyNhNtblrpPIf9DSevz27cCXwaIiKFOs6zISE8NrcC1+Y55BPhH4E8g7fc4Ig5J+lPg22QjQG6MiB01N2u6RcBGZf/0/ShgU0R8U9I9ZN1hlwOPk333kpKPADflH4A/Af6IvP2JtvkWSb8G/IpsG31e0rWk295S+BRyM7PE1T3qw8zMenBQm5klzkFtZpY4B7WZWeIc1GZmiXNQm5klzkFtZpa4/w/fEy5lCd/iiQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD8CAYAAABekO4JAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiuUlEQVR4nO3df/BddX3n8ecLEkBMWKERyC8m0UmQ6FTpfsuojNaCFQa7xB+Di7PSqGxZuyjQxRWCVrdjs0PVpmWn0zrfLTpxGsQsyZZUKRgQdZiFYBJjSfKFGg1iIALZooRFAiHv/eN8bnLyzf1x7q/veZ9z34+Z7+Tec8+59/O95/N5fT/53PP5XJkZIYQQ/Dqm7AKEEEJoL4I6hBCci6AOIQTnIqhDCMG5COoQQnAugjqEEJyLoA4hhD5Imi/pXkkTkrZLujptf6Ok+yU9JOkfJZ2UO2a5pJ2SHpF0QcfXiOuoQwihd5JmA7PNbIukmcBm4D3AKuCTZvY9SR8FFprZn0haAnwdOAeYA9wNLDazl1u9RvSoQwihD2a2x8y2pNv7gAlgLnAm8P202wbg/en2UuBWM9tvZruAnWSh3dK0YRS8W7NmzbIFCxYcuv8vD+9put/B4+v5d+Xl6WWXIFTNsS+VXYKp87qFpx1xf/PmzXvN7NX9POfb33GCPfOvBwvtu+2hl7YDL+Q2jZvZeLN9JS0AzgY2AtuAi4HbgUuA+Wm3ucADucN2p20tuQjqBQsWsGnTpkP3f++tf9bV8fsWvmLQRSrdc3Pr+UcpdG/G48UCpapm7vp128c3/J/PHHFf0s/6fc1n/vUg6+6YVWjfxfP3vGBmY532kzQDWAtcY2bPpuGO/yHps8B64MXGrk0ObzsG7SKo+9XuRFc1xIs0zgjz6qt7CEPnIK4DSdPJQnq1ma0DMLOHgXelxxcD70677+Zw7xpgHvBEu+evRVC306qSVDXA83pp5BHuwzMKodvOKARyM5IE3AxMmNnK3PZTzewpSccAnwG+nB5aD9wiaSXZh4mLgAfbvUbtg7qVOgd4O1MZJh7+KIx6eA7aqIZxB+cClwEPSdqatt0ALJJ0Zbq/DvgqgJltl7QG2AEcAK5sd8UHjHBQtzKqAT4MEZLVFYFcnJndR/NxZ4CbWhyzAlhR9DUiqAuq4zh4GG0RxtURQT0ARSp8hHmYShHC9RJBPUV6bTgR8CFCN0RQO1dmI40/EodFWIYyRVCHliKcQvCh/OunQgghtBVBHUIIzkVQhxCCcxHUIYTgXAR1CCE45/Kqj2k/ebzjPgde23b51hBCqA2XQV1EuzCPEA/BtyKdsXBYZYO6nVaVIAI8hKkTYTw4tQzqViLAQxisCOOpMVJB3UoMo4TQWoRx+SKoO4gPNkOdRQj3T9J84GvA6cBBsi+/vUnSm8i+1eUEsi8I+M9m9mA6ZjlwOfAycJWZ3dXuNSKoB6Dfyj7VQR+NczC8/IGO81m6A8C1ZrZF0kxgs6QNwBeAPzWzf5J0Ubr/DklLgEuB15N9Fdfdkha3+5aXCGoHoqFVU5y3AGBme4A96fY+SRPAXLJvFj8p7fZvOPwFtkuBW81sP7BL0k7gHOD+Vq8RQR1CCAMiaQFwNrARuAa4S9KXyCYXvjXtNhd4IHfY7rStpQjqEMLI+dXBE7jjuSUF994zS9Km3IZxMxufvJekGcBa4Boze1bSnwF/bGZrJX2A7JvK30nz71e0diWIoA4hhPb2mtlYux0kTScL6dVmti5tXgZcnW7/L+Dv0u3dwPzc4fM4PCzSVKz1EUIIfZAkst7yhJmtzD30BPA76fZ5wI/T7fXApZKOl7QQWAQ82O41okcd3Fy9UEfxgeNIOBe4DHhI0ta07QbgD4GbJE0DXgCuADCz7ZLWADvIrhi5st0VHxBBXTsRur4UPR8R6NVlZvfRfNwZ4N+2OGYFsKLoa0RQV1CEcf1MPqcR3CEvgroiIpxHSwR3yIugdi4COsDhehCBPZoiqB2KcA6tRGCPpghqJyKcQzfy9SVCu/4iqKdQhHEYhmb1KsK7XiKohyRCOZQpwrteOgZ1m7VWTwG+ASwAHgU+YGbPpGO6Wmu16iKUQxXElSTVVaRH3Wqt1Q8D95jZjZKuB64HrutlrdUqiVAOdRG97uroGNRt1lpdCrwj7bYK+C5wHT2stepVhHIYNdHr9qmrMepJa62elkIcM9sj6dS0W6G1ViVdQZr7fsYZZ3Rd8GGIYA7hSNHr9qFwUDdZa7Xlrk22HbXWalrPdRxgbGys7VqsgxaBHELvotc99QoFdYu1Vp+UNDv1pmcDT6XtXa+1OkwRyiEMV6s2FgE+OEWu+mi11up6soWxb0z/3p7bfouklWQfJnZca7VfEcYh+BMBPjhFetSt1lq9EVgj6XLgMeAS6G2t1ckieEOor27ad4R6pshVH+3WWj2/xTFdrbUaQgjNRKctE1/FFUIIfZA0X9K9kiYkbZd0ddr+DUlb08+juREJJC2XtFPSI5Iu6PQaMYU8hBD603RSoJn9+8YOkv4C+FW63fWkwOhRhxBCH8xsj5ltSbf3AY1JgcChCzI+AHw9bTo0KdDMdgGNSYEtRY86hDBynj1wAt/Z+7qCe98zS9Km3IbxNA/kKJMmBTa8DXjSzBrfQl5oUmBeBHUIIbS318zGOu00eVJg7qEPcrg3DQUnBeZFUIcQQp9aTApE0jTgfRz5beRdTwqMoA4Dt2/hK8ouwiEzd/267CKEmmszKRDgncDDZrY7t63rSYER1KFrnoK4kyJljTAPfWo6KdDM7iC7uiM/7NHTpMAI6tBSlQK5H81+zwjvUFS7SYFm9uEW27uaFBhBHYDRCeWiIryDJxHUIyqCuXsR3qEsEdQjIoJ5OCa/rxHcYRgiqGsqgrkcEdxhGCKoayKC2af8eYnQDr2KoK6wCOdqid526FUEdYVEMNdL9LZDURHUzkU4j4YI7dBOBLVDEc6jLUI7TBZB7USEc2gmQjtABHWpIpxDNyK0R1cEdQkioOG5uYP5cqEZjx8cyPNUTYT2aImgniKjFM6DCuFBvVbdw7xRtyKw6yuCeojqHs5TGcj9aFXOugV49LLrK4J6COoa0FUJ5qIm/z51Cu4I7XqJoB6QOoZz3YK5k7oGdwyNDJek+cDXgNOBg2RffntTeuwTwMfJviDgW2b2qbR9OXA58DJwlZnd1e41Iqj7EOFcb/n3og6hHb3soTkAXGtmWyTNBDZL2gCcBiwFftPM9ks6FUDSErJvfnk92Vdx3S1pcbtveYmg7kHdAjrCubO6hnYEdv/MbA+wJ93eJ2kCmAv8IXCjme1Pjz2VDlkK3Jq275K0EzgHuL/Va0RQFxThHBrqFNqj2st+4aXpTPzitKK7z5K0KXd/3MzGm+0oaQFwNrAR+CLwNkkrgBeAT5rZD8hC/IHcYbvTtpYiqNuIcA6dRGiPhL1mNtZpJ0kzgLXANWb2rKRpwMnAm4HfBtZIeg3Nv1/R2j13BPUkdQtniICeKo33ueqBDRHa3ZI0nSykV5vZurR5N7DOzAx4UNJBYFbaPj93+DzgiXbPH0FNhHMYrDr1siFCuxNJAm4GJsxsZe6hfwDOA74raTFwHLAXWA/cImkl2YeJi4AH273GyAZ1HcMZIqC9qVMvG+JDyBbOBS4DHpK0NW27AfgK8BVJ24AXgWWpd71d0hpgB9kVI1e2u+IDRiio6xrMEOFcBXUN7IZRDm4zu4/m484AH2pxzApgRdHXqHVQ1zmcIQK6iuo2LNIQwT1ctQjqugfyZBHQ9VC3XnZeBPdgVSqoRy2QJ4uArqc6B3ZDs7Yb4V2cy6Ae9UCeLAJ6NIxCYOd1aucR5IcVTgBJx0r6oaRvpvunSNog6cfp35Nz+y6XtFPSI5IuGEbBR8Fzc4+JkB5Bcd4z+xa+IjptSTc96quBCeCkdP964B4zu1HS9en+db0sOBKONCqN9Pl5/fUcT9xd7/dp1HrYobVCQS1pHvBusstJ/kvavBR4R7q9CvgucB09LDgSMnUL6H6DuN/nr0uQPzf3mAjrEVe0R/1XwKeAmbltp6VVozCzPY0l/Ci44IikK4ArAM4444zuSl0zdQjoYYdyL5qVqarhHb3r0dYxqCX9PvCUmW2W9I4Cz1lowZG0+tQ4wNjYWNsFSeqsqiHtMZiLmFzuqgV3BPZoKtKjPhe4WNJFwAnASZL+HnhS0uzUm54NNNZa7XrBkVFUtYCuajB3UtXgjsAeLR1rpZktN7N5ZraA7EPC75jZh8gWFlmWdlsG3J5urwculXS8pIUUWHBklFTpE/3n5x089DMqqvY7V6k+hd71cx31jWTrq14OPAZcAmBmXS84Mgqq0piqElBTIf9eeO9pRw+73roKajP7LtnVHZjZ/wXOb7FfVwuO1FkEdD003p8qBHaEdf24nJlYB1UI6Ajn7lWhlx296/qJoB4C7yEdAT0Y3nvZEdj14bOGVZT3D3aq9CFZlXh/Xz3XyTqQNF/SvZImJG2XdHXa/t8kPS5pa/q5KHdMV8tsRI96ALw3BM8hUieee9jRux6qA8C1ZrZF0kxgs6QN6bG/NLMv5XfuZZmNCOo+eQ7pCOhyRGCPljRDuzFLe5+kCZrMxs7pepmNCOoeRUCHTiKw/bKXjuHAEycW3X2WpE25++NpZvVRJC0AzgY2kk0W/LikPwA2kfW6n6HgMht5/mqQc1UYhw6+eD4n3uuzE3vNbCz30yqkZwBrgWvM7Fngb4HXAm8i63H/RWPXJoe3XUYjetQFea/MnsMg+O5dQ1x/3S9J08lCerWZrQMwsydzj/9P4JvpbtfLbPisNY5473F4v+IgHMnz+fJe172SJOBmYMLMVua2z87t9l5gW7rd9TIb0aNuwXuF9drYQzHPzzvouncNozt+3YNzgcuAhyRtTdtuAD4o6U1kwxqPAv8JeltmI4J6Eu8BDRHSdVGF4RCIwO7EzO6j+bjzHW2O6WqZjQjqJAI6lCUCO3Qy8kFdhYCGCOlR4Hk4BCKwyzSyQR0BHTzyHtYQgV2GkQrqqoRzQ11Detqc5wf6fF1MXKgE70MhDfn2FKE9XCMR1FULaKhHSA86kLt9naoHeBV61w3Ryx6u2gZ1FcMZqh3QUxXMRU0uTxWDuyq964boZQ9HbYK6qsGcV7WQ9hbMnVQ5uKvUu26I0B6cSgd1HcIZqhXQVQvndvK/SxVCu4ph3RCh3R/3QV2XMG6lKiFdp4BupvH7eQ/sqg2FNNOsTUd4t+cyqOsezg3eQ7ru4dxMVXrZVe5dNzMqbb5XLoO67iKgq8F7L7sOvetQTJzhKeY5pKfNeT5Cugnv74vnOhUGI4J6CnltUN6DyAvP75PXuhUGI4Y+poDXRuQ1dLzzOiRSt3HrcFgE9ZB5DOkI6MHwGNgxbl1PcTaHyFtIe/6ve5V5fF+91b06kzRf0r2SJiRtl3T1pMc/KckkzcptWy5pp6RHJF3Q6TWiRz0EHhuJtyCpo2lznnfXu46e9ZQ4QPYN41skzQQ2S9pgZjskzQd+D3issbOkJcClwOuBOcDdkha3+5aXOIsD5i2kPfb26szb++35Oxrrwsz2mNmWdHsfMAHMTQ//JfApjvyW8aXArWa238x2ATuBc9q9RvSoB8hbg/AUGKMmete+HfNiV+P4syRtyt0fN7PxZjtKWgCcDWyUdDHwuJn9KPv+20PmAg/k7u/mcLA3FUE9ABHQoRlvHzZGWPdsr5mNddpJ0gxgLXAN2XDIp4F3Ndu1yTZrsu2QOGt9ipAOnXg6JzEUMhySppOF9GozWwe8FlgI/EjSo8A8YIuk08l60PNzh88Dnmj3/NGj7pG3yu4pDMLRonddX8rGNW4GJsxsJYCZPQScmtvnUWDMzPZKWg/cImkl2YeJi4AH271GnKkeREiHXnk6V97qcYWdC1wGnCdpa/q5qNXOZrYdWAPsAO4Ermx3xQdEj7prniq3p0YfivP0QWNMkOmfmd1H83Hn/D4LJt1fAawo+hpxdgryNrYXIV1tHi/jC35FUBfgrRJ7auChP57Opbd6Hg4rFNSSXiXpNkkPp2mSb5F0iqQNkn6c/j05t39X0yM981Z5PTXsMBiezqm3/zmGTNEe9U3AnWb2OuCNZDNvrgfuMbNFwD3p/uTpkRcCfyPp2EEXfNi8VVhv/1UOg+Xt3Hqq+6FAUEs6CXg72eUnmNmLZvZLsmmQq9Juq4D3pNtdT4/0xFtAg79GHIbD2x9jj21hVBW56uM1wNPAVyW9EdgMXA2cZmZ7IJvrLqlxzWCh6ZGSrgCuADjjjDN6/gUGyWOl9NRwB+Gs058cyvNO/OK0oTxvGTxdFQJxzbUHRYJ6GvBbwCfMbKOkm0jDHC0Umh6Z5sqPA4yNjbWdPjlsHgMaqh/Swwrloq9V5fD2GNYQl/GVpUhQ7wZ2m9nGdP82sqB+UtLs1JueDTyV27+r6ZFlipAerKkM507yZaliaHsLa4jALkvHoDazX0j6uaQzzewR4HyyGTU7gGXAjenf29MhXU+PLIPXgIZqhbSnYG5ncjmrEtwewxpiOGSqFZ2Z+AlgtaTjgJ8CHyH7IHKNpMvJFsW+BLLpkZIa0yMPUGB65FSKgB6MqgR0K43yVyGwPYc1RO96KhQKajPbCjRb5u/8Fvt3NT1yKngOaKhOSFc9oCerSmB7W9QpLwJ7+Gq/1of3gAb/IV23cG6mKuPZXnvXEIE9TLUM6iqEc4PnkB6FgG7Gey/bc1hDBPYw1OqdrNoF+hHSvnl+DzzXnYZGe6xSm/Sq8j3qqlYCrw3NcziVwXPv2nvPOi962f2p3LuW/ysdIT04Z53+ZIR0G17fH491qZ2qt91mJM2XdG9asG67pKvT9s9L+uf0RQLfljQnd0xXC9e57lHX6WQ2eGxYHgPIq7NOf9Jd77pKPeu8ye27wr3tA8C1ZrZF0kxgs6QNwBfN7E8AJF0FfBb42KSF6+YAd0ta3O4yZpdBXceABn8hHQHdG4/DIVUN67xm7b4K4Z3WPGqse7RP0gQw18x25HZ7JYeX0ji0cB2wS1Jj4br7W72Gy6Cuowjp+vHWu65DWE82rE7bsS/BjMcLP/csSZty98fTWkVHkbQAOBvYmO6vAP4A+BXwu2m3QgvX5UVQTwFPIR0BPVjeetd1DGsH9ppZswl/R5A0A1gLXGNmzwKY2aeBT0taDnwc+BwFF67L8///ioqLkB4Nnt5bT3VuVEiaThbSq81sXZNdbgHen253vXBdBPUQeWkwXq9YqBtP77OXujcKJInsi1UmzGxlbvui3G4XAw+n2+uBSyUdL2khBRaui6GPIfHSULwExyjxMnYdwyBT5lzgMuAhSVvTthuAyyWdCRwEfgZ8DHpbuC6CeggipEOE9egws/toPu58R5tjulq4LoY+BixCOjR4OQde6mToXQT1AHlpEF4CIvg5F17qZuhNDH0MiIeG4CUUwpG8XMIXwyDVFT3qAYiQDkV4OEce6mroXgR1H6bNed5FxfcQAKEYD+fKQ50N3Ymg7pGXyu6h4YfueDhnXupvKCaCugdeKrmHBh964+HceanHobMI6i55qdweGnroj4dz6GX4LrQXQd0FLxXaQwMPg+HlXHqp26G5COqCvFRkLw07DI6Xc+qljoejRVAX4KUCe2nQYfC8nFsvdT0cKSa8dOCl4nppyMNw3qyHO+80yXf2vm4IJSlXrA8SWomgbsFLQEN9QrqXQO7mueoQ3hHWoZkI6iYipAdjkMHcy+tVNbgjrMNkMUY9SYR0/86b9fCUh7TncvTCy7n31B5GWfSoczxVSi8NtSjPgZgvW5V62Z561kD0rluQNB/4GnA62ZcEjJvZTZK+CPw74EXgJ8BHzOyX6ZjlwOXAy8BVZnZXu9eIHnUSId2bqvVaG+WtSpk91QVPbcSZA8C1ZnYW8GbgSklLgA3AG8zsN4F/AZYDpMcuBV4PXAj8jaRj273AyAe1t5lZnhpmO1UKu1aqUn5PdcJTW/HCzPaY2ZZ0ex8wAcw1s2+b2YG02wNkX2ILsBS41cz2m9kuYCdwTrvXGOmhD2+VzlODbKUq4VZU4/fxPiTiZRgE6vEh4zH7DzJz16+L7j5L0qbc/XEzG2+2o6QFwNnAxkkPfRT4Rro9lyy4G3anbS2NZFB7C2jwH9J1C+jJqhLYXozYuPVeMxvrtJOkGcBa4Bozeza3/dNkwyOrG5uaHG7tnnvkgtpjSHtW94CezHNge+pVN9Shdz0IkqaThfRqM1uX274M+H3gfDNrhPFuYH7u8HnAE+2ef2TGqL2NRed57U2PWkjneR2D91hXPLetqSBJwM3AhJmtzG2/ELgOuNjM8m/QeuBSScdLWggsAh5s9xojEdSeK5HHhuc1pMrg8X3wWGfAdzsbsnOBy4DzJG1NPxcBfw3MBDakbV8GMLPtwBpgB3AncKWZvdzuBWo99OG94nhscB6DqWznzXrY3VCIx2EQGLmxawDM7D6ajzvf0eaYFcCKoq9RqEct6Y8lbZe0TdLXJZ0g6RRJGyT9OP17cm7/5ZJ2SnpE0gVFCzMoVfivmLeQjl50ex7fH291KM97+6uajkEtaS5wFTBmZm8AjiW7WPt64B4zWwTck+73dDH3oFQhoMFfA/MWQJ7Fe1VcVdpjFRQdo54GvELSNOBEsk8olwKr0uOrgPek211fzN2vKlWICOnq8/SeeatPzVSpfXrVMajN7HHgS8BjwB7gV2b2beA0M9uT9tkDnJoOmQv8PPcUTS/mlnSFpE2SNj399NM9FT4qQO88/le+Sjy9f1UIa4j22o8iQx8nk/WSFwJzgFdK+lC7Q5psO+pibjMbN7MxMxt79atfXbS8h052FU+4lwblJWDqwMt76aVuFVHV9lumIkMf7wR2mdnTZvYSsA54K/CkpNkA6d+n0v5dX8zdTj6Yq3xyvTQkL8FSJ17eUy91rKg6tOupUiSoHwPeLOnEdGH3+WSLjqwHlqV9lgG3p9tdX8w9Wd1OoJcG5CVQ6sjLe+ulrnWrLh2yYel4HbWZbZR0G7CFbL76D4FxYAawRtLlZGF+Sdp/u6TGxdwHKHAxd515aThegqTOvFxv7fUa625EWB+p0IQXM/sc8LlJm/eT9a6b7d/Vxdx1FSE9eryEdaiXkZhCPsoipKeeh/fcSychDEYE9ZB4aCgeAmNUeXjvPdTBMBgR1EPgoYF4CIpR5+EceKiLoX8R1APmoWF4CIiQ8XAuPNTJ0J8I6gHy0CA8BEM4UpyT0K8I6hqJQPCr7HPjoRMRehdBPSBlN4SygyB0VvY5KruO1pWk+ZLulTSRloO+Om2/JN0/KGls0jFdLQUdQT0AZTeAsgMgFFf2uSq7rtbUAeBaMzsLeDNwZVrueRvwPuD7+Z17WQo6grpPZVf8sht+6F7Z56zsOls3ZrbHzLak2/vIltiYa2YTZvZIk0O6Xgo6groPZVf4sht86F3Z567sultXkhYAZwMb2+xWaCnovFp/Z+IwlV3Ry27oIVSZ9r/EtJ88XnT3WZI25e6Pm9n4Uc8pzQDWAteY2bPtXr7JtqOWgs6LoO5B2SEd6qHsdUHqsHjTFNlrZmPtdpA0nSykV5vZug7P1/VS0DH0UUHRm66Pss9ldDr6l5Z/vhmYMLOVBQ7peinoCOoulV2xy27YYfDKPqdl1+kaOBe4DDhP0tb0c5Gk90raDbwF+JakuyBbChpoLAV9JwWWgo6hjy6UXaHLbtDDdNGMHW0fv+O5JVNUknLEMEh1mdl9NB93BvjfLY7painoCOqCIqQHp1Mod3tMXUK87LAOfkVQFxAh3Z9egrmf569ycJcZ1tGr9iuCOgzFsMO56GtXObTLEGHtU3yY2EH0prtz0YwdpYb0ZN7KU0TZ57zsOh+OFkHdRtkVtuwG2w3vgei9fJOVfe7LrvvhSBHULZRdUctuqEVVLQCrVN6y60DZbSAcFmPUTZRdQctuoEVUJexaaZQ/xrBDFUSPepKyQ7oKqh7Sed5/l7L/aEd78CGCOnStTr3QOv0uob4iqENP6hBwdfgdpkL0qssXQZ0TFbI7VQ66Kpc9jJ74MDH0pRF43sd6GyKgQxVFUIeB8B7YEdD9iRmL5YqgTmLYYzC8BXYEdKiDGKMOQ3HHc0tKD8myXz+EQYkedRiqfFhORS87wnl4YvijPBHUxLDHVJkcooMI7gjmUDZJ84GvAacDB8m+/PYmSacA3wAWAI8CHzCzZ9Ixy4HLgZeBq8zsrnavEUEdShMhG2riAHCtmW2RNBPYLGkD8GHgHjO7UdL1wPXAdZKWAJcCrwfmAHdLWtzu67hijDqEEPpgZnvMbEu6vQ+YAOYCS4FVabdVwHvS7aXArWa238x2ATuBc9q9hsxsCEXvjqR9wCNll6NLs4C9ZReiS1Urc9XKC1HmqXCmmc3s5wkk3Un2exdxAvBC7v64mY23eN4FwPeBNwCPmdmrco89Y2YnS/pr4AEz+/u0/Wbgn8zstlYF8DL08YiZjZVdiG5I2hRlHq6qlReizFNB0qZ+n8PMLhxEWfIkzQDWAteY2bNSq++7bfpFuG17zDH0EUIIfZI0nSykV5vZurT5SUmz0+OzgafS9t3A/Nzh84An2j1/BHUIIfRBWdf5ZmDCzFbmHloPLEu3lwG357ZfKul4SQuBRcCD7V7Dy9BH0/Ee56LMw1e18kKUeSp4K++5wGXAQ5K2pm03ADcCayRdDjwGXAJgZtslrQF2kF0xcmW7Kz7AyYeJIYQQWouhjxBCcC6COoQQnCs1qCW9UdL9kh6S9I+STso9tlzSTkmPSLqgzHJOJukTqVzbJX0ht91dmSV9XtI/S9oq6duS5uQec1feBkkXpnLtTLO6XJF0gqQHJf0o1YM/TdtPkbRB0o/TvyeXXdY8Sa+SdJukhyVNSHqL5zJLulrStvQeX5O2uS3v0JhZaT/AD4DfSbc/Cnw+3V4C/Ag4HlgI/AQ4tsyy5sr8u8DdwPHp/qmeywyclLt9FfBlz+VNZTs2lec1wHGpnEvKLtekMgqYkW5PBzYCbwa+AFyftl8P/HnZZZ1U7lXAf0y3jwNe5bXMZJNGtgEnkl34cDfZFRIuyzvMn7KHPs4km8UDsAF4f7rd9RTLKfRHwI1mth/AzBrXRross5k9m7v7Sg5fWO+yvMk5wE4z+6mZvQjcSlZeNyzzXLo7Pf0YracNly79j/XtZJeSYWYvmtkv8Vvms8hm8D1vZgeA7wHvxW95h6bsoN4GXJxuX8Lhi8DnAj/P7bc7bfNgMfA2SRslfU/Sb6ftbsssaYWknwP/Afhs2uy2vPgu2yGSjk2XYz0FbDCzjcBpZrYHsjUggFNLLOJkrwGeBr4q6YeS/k7SK/Fb5m3A2yX9hqQTgYvIMsJreYdm6EEt6e40xjT5ZynZcMeVkjYDM4EXG4c1eaopu46wQ5mnASeT/Tf3v5JdJ6kyy9yhvJjZp81sPrAa+HjjsLLKW4Dnsh1iZi+b2ZvIZpadI+kNJRepk2nAbwF/a2ZnA/+PbOjAJTObAP6c7H/bd5INgR0otVAlGfqEFzN7Z4dd3gUgaTHw7rSt6ymWg9SuzJL+CFhn2QDZg5IOki3uUlqZC7zHDbcA3wI+R8nvcQeey3YUM/ulpO8CF5KmDZvZnknThj3YDexOPX+A28iC2m2Zzexm0lCNpP9O9ju4Le+wlH3Vx6np32OAzwBfTg91PcVyCv0DcB4c+uNyHNnKYy7LLGlR7u7FwMPptsvyJj8AFklaKOk4srV715dcpiNIerWkV6XbrwDeSfbetpo2XDoz+wXwc0lnpk3nk82Oc1vmXEacAbwP+DqOyzssZU8h/6CkK9PtdcBXobcpllPoK8BXJG0jG6pZlnrXXst8Y2qYB4GfAR8D3++xmR2Q9HHgLrIrQL5iZttLLtZks4FVko4l6/CsMbNvSrqfJtOGHfkEsDr9Afwp8BFS+Z2Wea2k3wBeIqujz0hqOjW7zmIKeQghOFf2VR8hhBA6iKAOIQTnIqhDCMG5COoQQnAugjqEEJyLoA4hBOciqEMIwbn/DwONWB1a7dQ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_temp_section(rcelat)\n",
    "plot_temp_section(diffmodel2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x16044e310>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAocklEQVR4nO3deVhXZf7/8ecbEJBFkdUVUUERdyWXTM20SSuzmjZbbdGabKqZaZrm17SM35qldZammbTSMpeWmca1NBvNzBVTUQQVFBWRzQUQZL9/f4BdZJooH7g/y/txXZ8r+EDnvDziy+M597lvMcaglFLKdXnZDqCUUqpxtMiVUsrFaZErpZSL0yJXSikXp0WulFIuzsfGTsPDw01MTIyNXSullMvasmVLgTEm4sz3rRR5TEwMSUlJNnatlFIuS0QOnO19vbSilFIuTotcKaVcnBa5Ukq5OC1ypZRycVrkSinl4rTIlVLKxWmRK6WUi7Myjlwp5XyMMRSVVZFTWMaRwlMcKSzjRGklVyZEERsZZDue+hFa5Ep5mEPHSlmbXsCRE7VlnVNURvaJU+QUllFSUf2D739peRrje7dl2uhYerVvbSGxOh8tcqU8hDGGeZsO8sKSVE5VViMCkcF+tGvdku5RwYzqHkm71v60C/Gv/W/rlnh7Ce+ty2TO+gMs25HD6B4RPHJFLIM6h9r+5ah6xMYKQYmJiUYf0Veq+RScLOc3nyTzZVoeI+LCef66XkSHBtDCu2G3yQpPVfL+ukze/WY/x0srGdo1lEdGxzE8NgwRaeL06jQR2WKMSfzB+1rkSrm3/6Xl8uQnyRSVVfHb8fHcMywGL6+LK9+S8irmbzrIjDX7yCsup1+nEB4ZHcuY+MiL3qZqOC1ypTzMqYpqXly2iw82HCS+bTB/vW0APdoGO2TbZZXVfLIli399lUHW8VPEtw3mN+PiGR0f6ZDtq7PTIlfKg+w8XMhjC7aSkV/ClBFdeOKqHvj5eDt8P5XVNSzens0bq9LJLCjhL7cN4Lp+7R2+H1XrXEWuNzuVciPVNYa31mTw2oo9hAf5MfeBIQyPDW+y/bXw9uLGgR0Z17stk2dt5hcfbsNbhGv6tmuyfaof0iJXyk1kHS/llx9uZ1PmMa7p244Xr+9NSIBvs+w7wNeHWZMv4Z53N/Hogq14e8G43lrmzUWf7FTKDezKLmL8X79m15EiXr25H29MGtBsJX5aoJ8Ps+8bTL+OrXlk3lZWpOQ06/49mRa5Ui7uRGkFD36QRICvN8seHcFPB3W0NiQwqK7Me3VozbR53/Jlaq6VHJ5Gi1wpF1ZdY3h0wTZyCsv4552DiA4LsB2JVv4teP++wfRs14qfffAtq3fn2Y7k9rTIlXJhr67YzZo9+Uyf2JuB0W1sx/lO65YtmHPfEOKigpg6Zwtr9uTbjuTWGl3kIuIvIptEZLuIpIjI7x0RTCn14z7feYQ3V2cwaXAnJg2Oth3nB1oHtOCD+4fQLSKIKe8n8U16ge1IbssRZ+TlwBXGmH5Af2CciAx1wHaVUuewN7eYX320nf6dQnj+ul6245xTm0Bf5j4whC7hgdz/3mbWZxy1HcktNbrITa2TdZ+2qHs1/1NGSnmIorJKHpyzhZa+3vzzzoFN8qCPI4UG+vLBA0Po1CaA+2ZvZtP+Y7YjuR2HXCMXEW8R2QbkAV8YYzae5XumikiSiCTl5+v1MqUuRk2N4ZcfbufgsVL+cftA2rVuaTtSg4QH+TFvylDah/gzedYmUrILbUdyKw4pcmNMtTGmP9ARGCwivc/yPTOMMYnGmMSIiAhH7FYpj/PGqnRWpuby9DU9GdI1zHacCxIR7Mf8KUMJ8vPhsQXbKKv84dzn6uI4dNSKMeYEsBoY58jtKqVgVVoer6/cww0DOjD50hjbcS5KZCt/Xrm5H+l5J/nTZ2m247gNR4xaiRCRkLqPWwJjAf0dUsqBMgtKeHTBVnq2bcUfbujj0nOAj+weweRLY5i9LlOHJTqII87I2wGrRCQZ2EztNfIlDtiuUoraOcAfnLMFby/hrbsG0dLXuW9uNsRT4+OJiwziiY+3c7ykwnYcl+eIUSvJxpgBxpi+xpjexpjpjgimlKpdnu3JfyezN6+Yv08aQKdQ+09uOoJ/C29ev7U/x0srePq/O7AxnbY70Sc7lXJi736TydLkI/z6qnhGxLnXIIHeHVrziyu7s2xHDv/59rDtOC5Ni1wpJ5VbVMarK3YzJj6Sh0Z1tR2nSTw4shuDY0J5blEKh46V2o7jsrTIlXJSf/48japqw3MTern0zc0f4+0lvHpLPwB+9dF2qmv0EsvF0CJXygltPXic/3x7mPtHdHGKGQ2bUqfQAH5/XS82ZR7jrTUZtuO4JC1ypZyMMYbpS3YREezHtNGxtuM0ixsHduDqPm15/Ys97DysT31eKC1ypZzMwm3ZbD14gl9f1YMgP89YjVFEePH6PoQG+vL4h/rU54XSIlfKiZRWVPGnz9Lo06E1Nw3saDtOs2oT6KtPfV4kLXKlnMi/VmeQU1TGcxMS8PJyzxucP2ZEnD71eTG0yJVyElnHS3lrzT4m9GtPYkyo7TjW6FOfF06LXCkn8afP0hCpLTJPVv+pz2cW7rQdxyVokSvlBDbtP8aS5CM8OLIbHUJcY47xptS7Q2t+fkUcS5KP6CWWBtAiV8qymhrD9CUptGvtz0OjutmO4zQeHNWVmLAAnl+UQnmVjmL5MVrkSln2yZYsdh4u4qnx8W4xs6Gj+Pl48/x1vdhXUMI7a/fbjuPUtMiVsqi4rJKXlqcxqHMbruvX3nYcp3N5j0iu6hXF379M5/CJU7bjOC0tcqUsemNVOgUnK3j22gS3nU+lsZ65NgGD4YUlu2xHcVpa5EpZkllQwqy1mfx0YEf6dQqxHcdpdWwTwM+viOOznTl64/MctMiVsuTFZan4eAtPjuthO4rTe2BEF7qEB+qNz3PQIlfKgm/SC/hiVy7TRscS1crfdhyn5+fjzXMTEthXUMLbX+uNzzNpkSvVzKqqa5i+eBedQlty/2VdbMdxGZf3iGRcr7b8/X979cbnGRpd5CLSSURWiUiqiKSIyGOOCKaUu1qSfITducU8Na4n/i10uOGFeGZCAgD/t1hvfNbniDPyKuBXxpiewFBgmogkOGC7SrmdmhrDP1dnEBcZxPjebW3HcTkdQlry8yvi+Dwlh6/0xud3Gl3kxpgjxphv6z4uBlKBDo3drlLu6Mu0PHbnFvPw6G4eObuhI+iNzx9y6DVyEYkBBgAbz/K1qSKSJCJJ+fn6N6nyPMYY/rEqnY5tWjKhrz78c7FOP/G5X298fsdhRS4iQcC/gceNMUVnft0YM8MYk2iMSYyIiHDUbpVyGev3HWXboRM8OKobPt46zqAxRnWP+O7GZ9bxUttxrHPIT5OItKC2xOcaY/7jiG0q5W7eXJVBeJAfNw/yrJV/msp3Nz71iU+HjFoR4B0g1RjzWuMjKeV+th86wdr0AqaM6KIjVRzk9I3P5Sm5rN6dZzuOVY44Ix8O3AVcISLb6l5XO2C7SrmNN1en08rfhzuGdrYdxa08MKILXfXGp0NGraw1xogxpq8xpn/da5kjwinlDvbmFrM8JZfJl8YQ5OdjO45bOX3jM/NoKTPX7LMdxxq946JUE/vn6gxatvBm8nB9irMpjKy78fnGqnSPvfGpRa5UEzp0rJSF27O5fUg0oYG+tuO4rdM3Pl9Ykmo5iR1a5Eo1oRlr9uEltddyVdOp/8SnJ051q0WuVBPJKy7jw6RD/HRgR9q11gWVm5onP/GpRa5UE3l3bSZV1TU8qAsqN4v6U9162hqfWuRKNYHC0ko+2HCAq/u0o0t4oO04HuPyHpH8JMHz1vjUIleqCby/PpOT5VU8fHms7Sge55lrE6gxhheXes4Tn1rkSjlYaUUVs9ZlckV8JAntW9mO43E6hQYwbXQsy3bksHZvge04zUKLXCkHW7DpEMdKKnj4cr02bsvUkV3pHBbAs4t2UlFVYztOk9MiV8qBKqpqmPn1PgZ3CSUxJtR2HI/l38Kb5yf0Yl9+Ce9+4/43PrXIlXKg/249zJHCMqaN1mvjto2Oj2Rszyj+9uVejhS6941PLXKlHKS6xvDPrzLo1b4VI+PCbcdRwHMTEqiuMby41L2f+NQiV8pBPtt5hP0FJUwbHUvt7M7Ktk6hAfzs8m4sST7Cugz3vfGpRa6UAxhjmLlmH13DA7mqly6q7EweGtWNTqEteXZhCpXV7nnjU4tcKQfYeugE27MKuXd4DN66qLJTOX3jMz3vJLPc9ManFrlSDvDeukyC/Xy4caAu4+aMxvSMYkx8JH9duZfcojLbcRxOi1ypRsotKmNp8hFuTuxEoC4c4bSem9CLyhrDdDdc41OLXKlGmrvxINXGcPcwXcbNmUWHBfDoFbEsTT7Cou3ZtuM4lBa5Uo1QXlXNvI0HuKJHJDE6OZbTe2hUNwZEh/C7T3eQ7UaTajmkyEXkXRHJE5GdjtieUq5i2Y4jFJys4J5LY2xHUQ3g4+3FX27tT1WN4YmPt1NTY2xHcghHnZHPBsY5aFtKuQRjDLO+yaRbRCAj9AEgl9E5LJBnr01gXcZRt3l83yFFboxZAxxzxLaUchVbD50gOauQyZfG6ANALubWSzoxtmcULy3fze6cYttxGq3ZrpGLyFQRSRKRpPx8z1tTT7mf2d/okENXJSL86ad9aOXvw2MLtrr80nDNVuTGmBnGmERjTGJERERz7VapJpFbVMayHTrk0JWFB/nx0k19Scsp5rUVe2zHaRQdtaLURdAhh+7hivgobh8SzYyv97E+46jtOBdNi1ypC6RDDt3L767pSUxYIL/6aBuFpyptx7kojhp+OB9YD/QQkSwRud8R21XKGS1Nrh1yOHl4jO0oygECfH14/db+5BaX89xC1xxB7ahRK5OMMe2MMS2MMR2NMe84YrtKORtjDLPX1Q45vCxWhxy6i/6dQvj5FbH8d1s2i13wqU+9tKLUBdAhh+7rkdGx9O8UwtOf7nC5FYW0yJW6ADrk0H35eHvx+q39qax2vac+tciVaqDTQw5vuUSHHLqrLuGBPHNtAt+kH2XWukzbcRpMi1ypBtIhh55h0uBOjImP5M+fp7Ht0AnbcRpEi1ypBqg/5LBzmA45dGciwp9v6ktUKz/ufmcjOw8X2o50XlrkSjWADjn0LOFBfsyfMpRg/xbc8fZGdmUX2Y70o7TIlTqP00MOYyODdMihB+nYJoD5U4YS4OvNHW9vIC3Hectci1yp8zg95PCeYZ11yKGHiQ6rLXNfHy/umLmRvbnOOVOiFrlS56FDDj1bTHgg86cMxctLmDRzI+l5J21H+gEtcqV+RJ4OOVRA14gg5k8ZChhun7mB/QUltiN9jxa5Uj9i/qZDVNUY7hqqQw49XWxkEPOmDKWqxjBpxgYOHHWeMtciV+ocKqtrmLfpACO7R+gshwqA7lHBzH1gCGVV1UyasYFDx0ptRwK0yJU6p5W7csktKuduPRtX9fRs14oP7h9CSUU1k2Zu4PAJ+/OyaJErdQ5zNhygQ0hLRsdH2o6inEzvDq354P4hFJ6qZNKMDdYn2dIiV+os0vOKWZdxlNuHROPtpUMO1Q/16diaOfcP4XhJBdf8bS1vf72Psko7a39qkSt1Fh9sOIivtxe3XtLJdhTlxPp3CuHjnw2jV/tWvLA0ldGvrGbBpoNUVdc0aw4tcqXOUFJexb+3ZHF1n7aEB/nZjqOcXHzbVsy5fwjzpgwhqpU/T/1nB1e+vobF27ObbSpcLXKlzrBwWzbF5VXcpbMcqgtwabdwPn34UmbenYivtxc/n7+Va/++llVpeRjTtIWuRa5UPcYY3l+fSc92rRgY3cZ2HOViRIQrE6JY9tgI/nJrf06WV3Hv7M3c8tZ6Nu0/1mT7ddTiy+NEZLeIpIvIU47YplI2bDlwnLScYu7WeVVUI3h7CdcP6MDKX47ihet7c+BoKbe8tZ7JszaRnuf4+VoaXeQi4g38AxgPJACTRCShsdtVyoY5Gw4Q7OfDxP7tbUdRbsDXx4s7h3bmq1+P5qnx8Ww7dIJTFY6/EeqIySMGA+nGmH0AIrIAmAjscsC2vyensIwaY2gf0tLRm1aKgpPlLNtxhDuGdCbAV+dVUY7T0tebh0Z1Y/KlMfi38Hb49h1xaaUDcKje51l17zncy8t3c/krq3lx6S6Ol1Q0xS6UB/tw8yEqqw136pOcqok0RYmDY87Iz3Yh8Qe3aEVkKjAVIDo6+qJ29Isr4xCBd9buZ8GmQzw4qiv3XdZFz55Uo1XXGOZuOMDw2DBiI4Nsx1HqgjjijDwLqP/UREcg+8xvMsbMMMYkGmMSIyIiLmpHHdsE8MrN/fj88ZEM7RbGKyv2MPKl1cxZn0lFVfMOwFfu5X9peWQXluksh8olOaLINwNxItJFRHyB24BFDtjuOXWPCmbm3Yn8+2eX0jUikGcWpjD2ta9YuO1wsw3AV+5lzoYDtG3lz9ieUbajKHXBGl3kxpgq4BFgOZAKfGSMSWnsdhtiUOc2fDh1KLPuvYRAPx8eW7CNa/++ltW7m34AvnIf+wtKWLMnn9uHROPjrY9WKNfjkIvLxphlwDJHbOtCiQije0QyKi6CxcnZvLJiN5NnbWZIl1BeuqkvncN0Hmn14+ZuOICPl3CbzquiXJTbnH54eQkT+3fgy19ezvSJvUjLKeY2J5r4XTmnUxXVfLwli6t6tyWylb/tOEpdFLcp8tN8fby4e1gM86cMpbSimtvf3kC2E0z8rpzT4uRsCk9V6k1O5dLcrshPS2jfivfvG8yJkkrueHsjeUVltiMpJ2OMYc76A3SPCmJIl1DbcZS6aG5b5AD9OoUw+75LyC0q4/a3N1Jwstx2JOVEtmcVsuNwIXcN1XlVlGtz6yIHGNQ5lHcnX0LW8VLufHujPhGqvjNn/QECfb25fkCTPIisVLNx+yIHGNo1jLfvvoR9BSXc/e4mCk9V2o6kLDteUsHi5GxuGNiBYP8WtuMo1SgeUeQAl8WF8687B5KWU8TkWZs4WV5lO5Ky6KOkQ1RU1XDX0BjbUZRqNI8pcoAr4qP4+6SBJGcVct+szZRWaJl7ouoaw5wNBxjcJZQebYNtx1Gq0TyqyAHG9W7LX27tT9KBY0x5P8naqtfKni9Tc8k6fop7L42xHUUph/C4IgeY0K89r9zcj3UZR/nZB1sor9Iy9ySz12XSvrU/VybovCrKPXhkkQPcOLAjf7ihD6t25/P4gm06N4uH2J1TzLqMo9w1LEbnVVFuw6N/kicNjub/XR3PZztzmLvxoO04qhnMXpeJn4+Xzqui3IpHFznAA5d1ZURcOC8uTWV/QYntOKoJnSit4NOtWdwwoANtAn1tx1HKYTy+yL28hJdv6kcLb+GXH22jqloXqHBXH24+RFllDffoTU7lZjy+yAHatvbn/67vzdaDJ/jXVxm246gmUFVdw/vrDzC0ayg927WyHUcph9IirzOxfweu7duOv6zcy87DhbbjKAdbmZrH4ROnmHxpF9tRlHI4LfJ6Xri+N6GBvvziw206vtzNzF63nw4hLRnbM9J2FKUcTou8npAAX166qS97807yyvLdtuMoB0k9UsSGfce4e1hnHXKo3JL+VJ/h8h6R3Dk0mne+2c/6jKO24ygHeG9dJv4tvLhVhxwqN6VFfhb/7+qexIQF8sTH2yku05kSXdnxkgo+3XqYGwZ0JCRAhxwq99SoIheRm0UkRURqRCTRUaFsC/D14dVb+nGk8BS/X7zLdhzVCAs2H6K8qobJOuRQubHGnpHvBG4E1jggi1MZGN2Ghy+P5ZMtWSxPybEdR12Equoa5qzP5NJuYTrLoXJrjSpyY0yqMcZt7wo+OiaOXu1b8dv/7CC/WJeJczVf7Molu7BMz8aV22u2a+QiMlVEkkQkKT8/v7l22yi+Pl785db+nCyv4rf/SdaJtVzMrHWZdGzTkjE9dZZD5d7OW+QislJEdp7lNfFCdmSMmWGMSTTGJEZERFx84mYWFxXMk1f1YGVqHh8nZdmOoxooJbuQTfuPcc+wGLy9dGFl5d58zvcNxpixzRHEmd03vAsrU3P5/eIUhnULo1NogO1I6jzeW5dJyxbe3JKoQw6V+9Phhw3g5SW8cnM/AH6/OMVyGnU+x0oq+O+2bG4c2IHWAbqwsnJ/jR1+eIOIZAHDgKUistwxsZxPxzYBPDomjpWpeazanWc7jvoR8zcdpEKHHCoP0thRK58aYzoaY/yMMVHGmKscFcwZ3Tu8C10jApm+eJcuD+ekKqtr+GDDAS6LDScuSoccKs+gl1YugK+PF89N6MX+ghLeWbvfdhx1FitScjmiQw6Vh9Eiv0CjukdwZUIUb/wvnZzCMttx1Blmr9tPdGgAo+N1lkPlObTIL8Iz1yRQVWP4w7JU21FUPdsPnWBz5nHuHtZZhxwqj6JFfhGiwwJ4aFQ3Fm3PZuM+nSHRWby5Op1W/j46y6HyOFrkF+lno7rRIaQlzy1K0XU+ncCe3GKWp+Qy+dIYgv11yKHyLFrkF6mlrze/u6YnaTnFzN140HYcj/fP1RkE+Hpz73Bdyk15Hi3yRhjXuy3DY8N4dcVujp7USbVsOXi0lEXbs7l9cDRtAnXOceV5tMgbQUR4fkIvSiuqeWWF204C6fTeWpOBtwgPjOhqO4pSVmiRN1JcVDCTL41hweZDJGedsB3H4+QVlfFxUhY/HdSRtq39bcdRygotcgd4bGwcYYF+PLswhZoaneq2Oc38eh9VNTU8NErPxpXn0iJ3gGD/Fjw1Pp5th07w7291qtvmcrykgrkbD3Jdv/Z0Dgu0HUcpa7TIHeTGAR0YEB3Cnz9Po0gXbG4Ws9dlUlpRzc8uj7UdRSmrtMgdxMtLmH5db46WVPDXlXttx3F7J8urmL0ukysTonQ9TuXxtMgdqE/H1tx2STSz12WyJ7fYdhy3Nm/jAQpPVTJttJ6NK6VF7mC/vqoHQX4+PL8oRdf4bCJlldXM/Ho/l8WG079TiO04SlmnRe5goYG+PPGT7qzLOMpnO3Nsx3FLH2/JIr+4nIdHd7MdRSmnoEXeBG4f0pmEdq14YckuSiuqbMdxK5XVNbz1VQYDokMY1jXMdhylnIIWeRPw9hKmT+xFdmEZb67KsB3HrSzenk3W8VNMuzwWEZ2qVinQIm8yiTGh3DCgAzPW7COzoMR2HLdQU2N4c3UG8W2DuUIXjlDqO41dfPllEUkTkWQR+VREQhyUyy38dnw8LbyF6Ut22Y7iFlbsyiE97yQPj47FSxeOUOo7jT0j/wLobYzpC+wBftv4SO4jspU/j4/tzv/S8vgyNdd2HJdmjOEfqzKICQvgmj7tbMdRyqk0qsiNMSuMMafv5m0AOjY+knuZPDyG2Mggfr94F2WV1bbjuKyv9xaw43AhD43qpsu4KXUGR14jvw/47FxfFJGpIpIkIkn5+fkO3K1za+HtxfMTenHwWCkz1+yzHcdl/WNVOm1b+XPDwA62oyjldM5b5CKyUkR2nuU1sd73PA1UAXPPtR1jzAxjTKIxJjEiIsIx6V3EZXHhXN2nLf9YnU7W8VLbcVxOUuYxNu4/xpSRXfHz8bYdRymnc94iN8aMNcb0PstrIYCI3ANcC9xh9FHGc3r6mgQA/rAs1XIS12KM4dUVewgN9GXSYF1UWamzaeyolXHAb4DrjDF6qvkjOoS05JHRsSzbkcPavQW247iM5Sk5rN93lMfHxhHg62M7jlJOqbHXyN8AgoEvRGSbiPzLAZnc1gMjutI5LIDnFu2koqrGdhynV1ZZzQtLU+kRFcztg6Ntx1HKaTV21EqsMaaTMaZ/3eshRwVzR/4tvHn22gQy8kt4b12m7ThO7521+8k6fopnJyTg463Pril1Lvqno5mN6RnFFfGR/GXlHvKKymzHcVo5hWX8Y1U6V/WKYnhsuO04Sjk1LXILnr02gcpqwx8/S7MdxWn9+fM0qmoMT1+dYDuKUk5Pi9yCmPBApo7syqdbD7Np/zHbcZzOlgPH+XTrYaaM6EJ0WIDtOEo5PS1ySx4e3Y32rf15duFOqqr1xudpNTWG3y9OIaqVHw/rWpxKNYgWuSUBvj787toE0nKKma03Pr/zybdZJGcV8tT4eAL9dLihUg2hRW7R+N5tGdszipeW79Y1PoHiskpe+nw3A6JDmNhPH8VXqqG0yC0SEf54Yx+C/Xz4xYfbPH5s+Rur0ik4Wc7zE3rpNLVKXQAtcssigv34w419SMku4m9f7rUdx5r9BSW8u3Y/Nw3qSD9dUFmpC6JF7gSu6tWWmwZ15M3V6Ww5cNx2HCteXLoLX28vnhzXw3YUpVyOFrmTeG5CAu1at+RXH23zuAWbv9qTz8rUPH4+Jo7IYH/bcZRyOVrkTiLYvwWv3NyPA8dKPWqGxMrqGqYvTiEmLIB7h8fYjqOUS9IidyLDuoVx//AufLDhIKt359mO0yzmrD9ARn4Jv7smQecaV+oiaZE7mSeu6kH3qCCe/CSZ4yUVtuM0qaMny3l95R5GxIUzpmek7ThKuSwtcifj38Kb127pz/HSCn63cCfuvFbHq1/sobSimucmJCCiww2Vulha5E6od4fWPD62O0uTj7Boe7btOE1iXXoBCzYd5O5hnYmNDLYdRymXpkXupB4c2ZWB0SE889+dHCk8ZTuOQx06Vsq0ed/SNSKIX/1Ehxsq1Vha5E7Kx9uL127pT2W14dcfJ1NT4x6XWE5VVPPgnC1U1Rhm3DWIIJ1PRalG0yJ3YjHhgTx9TU/Wphfw/vpM23EazRjDb/6dTGpOEX+7bQBdI4JsR1LKLWiRO7k7hkRzeY8I/vhZGul5J23HaZS3v97Pou3ZPPGTHoyO11EqSjlKo4pcRP5PRJLrFl5eISLtHRVM1RIRXvppX1r6evPo/K0Ul1XajnRRvt6bzx8/S+XqPm15+PJutuMo5VYae0b+sjGmrzGmP7AEeLbxkdSZIlv58/ot/dmTW8y9szZTUu5aj/AfPFrKz+dvJS4ymJdv6qdDDZVysEYVuTGmqN6ngYB73JFzQqPjI/nrbQP49uBxHngvibLKatuRGqS0ooqpc5KoqTHMuHuQLhahVBNo9DVyEXlRRA4Bd/AjZ+QiMlVEkkQkKT8/v7G79UjX9G3Hq7f0Y8P+o0yds4XyKucuc2MMv/4kmT25xfz99oF0Dgu0HUkpt3TeIheRlSKy8yyviQDGmKeNMZ2AucAj59qOMWaGMSbRGJMYERHhuF+Bh7lhQEf+dGMf1uzJZ9rcb516MYp/fbWPpclHeHJcPKO66++5Uk3lvP/ONcaMbeC25gFLgecalUid162XRFNRVcMzC1N4/MOt/O22Afh4O9cApNW783hpeRrX9m3HgyO72o6jlFtr1AVLEYkzxpxe1uY6IK3xkVRD3DUshvKqGl5Ymoqv93ZevaU/3k6yPFpmQQmPzt9Kj6hgXrqpr97cVKqJNfbO059EpAdQAxwAHmp8JNVQD4zoSnlVDS8v342fjzd/vLGP9bUuT5bX3tz08hJm3p1IgK/e3FSqqTXqT5kx5qeOCqIuzrTRsZRXVvO3/6Xj6+PF9Im9rJ0BHzhawi8+3EZ63knev28InUIDrORQytPo6ZIb+MWV3SmvquGtNfvw8/Hi6Wt6NmuZG2OYt+kgLy5NxdtL+NukAVwWF95s+1fK02mRuwER4anx8ZRX1fD22v34t/DmiauaZ1bB3KIynvwkma/25HNZbDgv39yXdq1bNsu+lVK1tMjdhIjw3IQEyquqeWNVOjuzC3lsTBwDots02T4Xb8/md//dSXlVNdMn9uLOIZ2tX6NXyhNpkbsREeHF6/vQOSyQt77K4IY31zGyewSPjYljUGfHFfqJ0gqeXZjCou3Z9O8Uwmu39NOZDJWySGwsJZaYmGiSkpKafb+e5GR5FXPWH2Dm1/s4VlLBiLhwHh8bx6DOoY3a7ld78nnyk+0cPVnB42PjeGhUN6cbw66UuxKRLcaYxB+8r0Xu3krKq5iz4QAz1tQW+mWx4Tw2No5LYi6s0EsrqvjDslQ+2HCQuMggXr+1P707tG6i1Eqps9Ei93ClFVV8UFfoBScrGB4bxmNjujO4SyjGGIrKqsgvLq99nSwnr6iM/JPl3723J7eYvOJyHrisC7/6SQ/8W3jb/iUp5XG0yBVQW+hzNxzkrTUZFJysIDLYjxOnKs86Z4uvtxcRwX6EB/sRFezHfZd1YWjXMAuplVJw7iLXm50eJsDXhykju3Ln0M7M23SQlOxCIoL8iAiue9X7uHXLFvp4vVIuQIvcQ7X09eb+y7rYjqGUcgAdbqCUUi5Oi1wppVycFrlSSrk4LXKllHJxWuRKKeXitMiVUsrFaZErpZSL0yJXSikXZ+URfRHJp3aNT2cRDhTYDnEBNG/TcaWsoHmbmrPl7WyMiTjzTStF7mxEJOls8xc4K83bdFwpK2jepuYqefXSilJKuTgtcqWUcnFa5LVm2A5wgTRv03GlrKB5m5pL5NVr5Eop5eL0jFwppVycFrlSSrk4jy5yEflQRLbVvTJFZFvd+zEicqre1/5lOSoi8ryIHK6X6ep6X/utiKSLyG4RucpmztNE5GURSRORZBH5VERC6t53umN7moiMqzuG6SLylO089YlIJxFZJSKpIpIiIo/VvX/Onwvb6v5M7ajLlVT3XqiIfCEie+v+28Z2TgAR6VHvGG4TkSIRedyZj299eo28joi8ChQaY6aLSAywxBjT23Ks74jI88BJY8wrZ7yfAMwHBgPtgZVAd2NMdbOH/H6unwD/M8ZUicifAYwxv3HGYwsgIt7AHuBKIAvYDEwyxuyyGqyOiLQD2hljvhWRYGALcD1wC2f5uXAGIpIJJBpjCuq99xJwzBjzp7q/LNsYY35jK+PZ1P0sHAaGAPfipMe3Po8+Iz9NahemvIXaQnQ1E4EFxphyY8x+IJ3aUrfKGLPCGFNV9+kGoKPNPA0wGEg3xuwzxlQAC6g9tk7BGHPEGPNt3cfFQCrQwW6qizIReK/u4/eo/cvI2YwBMowxzvT0+Y/SIq81Asg1xuyt914XEdkqIl+JyAhbwc7wSN2linfr/ZO0A3Co3vdk4Xx/wO8DPqv3uTMeW1c4jkDt5SlgALCx7q2z/Vw4AwOsEJEtIjK17r0oY8wRqP3LCYi0lu7cbuP7J3XOeny/4/ZFLiIrRWTnWV71z7Ym8f3fuCNAtDFmAPBLYJ6ItLKc9Z9AN6B/Xb5XT/9vZ9lUs1wva8ixFZGngSpgbt1bVo5tA1g7jhdCRIKAfwOPG2OKOPfPhTMYbowZCIwHponISNuBzkdEfIHrgI/r3nLm4/sdH9sBmpoxZuyPfV1EfIAbgUH1/p9yoLzu4y0ikgF0B5KaMOp5s54mIjOBJXWfZgGd6n25I5Dt4Ghn1YBjew9wLTDG1N2MsXVsG8DacWwoEWlBbYnPNcb8B8AYk1vv6/V/LqwzxmTX/TdPRD6l9vJVroi0M8Ycqbvun2c15A+NB749fVyd+fjW5/Zn5A0wFkgzxmSdfkNEIupueCAiXYE4YJ+lfKcztav36Q3AzrqPFwG3iYifiHShNuum5s53JhEZB/wGuM4YU1rvfac7tnU2A3Ei0qXurOw2ao+tU6i7j/MOkGqMea3e++f6ubBKRALrbsoiIoHAT6jNtgi4p+7b7gEW2kl4Tt/717mzHt8zuf0ZeQOceT0MYCQwXUSqgGrgIWPMsWZP9n0viUh/av+5nwk8CGCMSRGRj4Bd1F7CmGZ7xEqdNwA/4IvaDmKDMeYhnPPYUje65hFgOeANvGuMSbEcq77hwF3ADqkbJgv8P2DS2X4unEAU8Gnd770PMM8Y87mIbAY+EpH7gYPAzRYzfo+IBFA7aqn+MTzrnztno8MPlVLKxemlFaWUcnFa5Eop5eK0yJVSysVpkSullIvTIldKKRenRa6UUi5Oi1wppVzc/wcCVTTIBYj0cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#  Plot the northward heat transport in this model\n",
    "Rtoa = np.squeeze(diffmodel2.timeave['ASR'] - diffmodel2.timeave['OLR'])\n",
    "plt.plot(diffmodel2.lat, inferred_heat_transport(Rtoa, diffmodel2.lat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}