{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling spectral bands with `climlab`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a brief introduction to the `climlab.BandRCModel` process.\n",
    "\n",
    "This is a model that divides the spectrum into 7 distinct bands: three shortwave and four longwave.\n",
    "\n",
    "As we will see, the process works much like the familiar `climlab.RadiativeConvectiveModel`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About the spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The shortwave is divided into three channels:\n",
    "\n",
    "- Channel 0 is the Hartley and Huggins band (extreme UV, 200 - 340 nm, 1% of total flux, strong ozone absorption)\n",
    "- Channel 1 is Chappuis band (450 - 800 nm, 27% of total flux, moderate ozone absorption)\n",
    "- Channel 2 is remaining radiation (72% of total flux, largely in the visible range, no ozone absorption)\n",
    "\n",
    "The longwave is divided into four bands:\n",
    "\n",
    "- Band 0 is the window region (between 8.5 and 11 $\\mu$m), 17% of total flux.\n",
    "- Band 1 is the CO2 absorption channel (the band of strong absorption by CO2 around 15 $\\mu$m), 15% of total flux\n",
    "- Band 2 is a weak water vapor absorption channel, 35% of total flux\n",
    "- Band 3 is a strong water vapor absorption channel, 33% of total flux\n",
    "\n",
    "The longwave decomposition is not as easily related to specific wavelengths, as in reality there is a lot of overlap between H$_2$O and CO$_2$ absorption features (as well as absorption by other greenhouse gases such as CH$_4$ and N$_2$O that we are not representing)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example usage of the spectral model"
   ]
  },
  {
   "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": [
    "First try a model with all default parameters. Usage is very similar to the familiar `RadiativeConvectiveModel`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (1,) \n",
      "  Tatm: (30,) \n",
      "The subprocess tree: \n",
      "Untitled: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "col1 = climlab.BandRCModel()\n",
    "print(col1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out the list of subprocesses.\n",
    "\n",
    "We now have a process called `H2O`, in addition to things we've seen before.\n",
    "\n",
    "This model keeps track of water vapor. We see the specific humidity in the list of state variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "AttrDict({'Ts': Field([288.]), 'Tatm': Field([200.        , 202.68965517, 205.37931034, 208.06896552,\n",
       "       210.75862069, 213.44827586, 216.13793103, 218.82758621,\n",
       "       221.51724138, 224.20689655, 226.89655172, 229.5862069 ,\n",
       "       232.27586207, 234.96551724, 237.65517241, 240.34482759,\n",
       "       243.03448276, 245.72413793, 248.4137931 , 251.10344828,\n",
       "       253.79310345, 256.48275862, 259.17241379, 261.86206897,\n",
       "       264.55172414, 267.24137931, 269.93103448, 272.62068966,\n",
       "       275.31034483, 278.        ])})"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col1.state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The water vapor field is initialized to zero. The `H2O` process will set the specific humidity field at every timestep to a specified profile. More on that below. For now, let's compute a radiative equilibrium state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 730 steps, 730.4844 days, or 2 years.\n",
      "Total elapsed time is 1.9986737567564754 years.\n"
     ]
    }
   ],
   "source": [
    "col1.integrate_years(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Field([-0.00148377])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check for energy balance\n",
    "col1.ASR - col1.OLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEfCAYAAABvWZDBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7cElEQVR4nO3deXxU1f3/8dc7O2Ffw5IAAdl3CAi1tVC1avVbrbb90UXRWu2irVqtW2vVr9W61Lb2q9aqdW/rXrV133BpZYewEwJhXwUUwg75/P64N3WcTCATksxM8nk+HvOYmXO3zxzCfOaee+45MjOcc865eKUlOgDnnHOpyROIc865WvEE4pxzrlY8gTjnnKsVTyDOOedqxROIc865WvEE4pxrUJIKJT0vabMkk/RwWP7f1xHrrpA0OQFhuhrISHQALvlIiufmoEIzW1FfsaQaScOB04GHvV6q9TAwFLgJ2AAsS2g0rtY8gbhYzop6/wXgAuA+4P2oZZsbJKLUMRy4DpgMrEhkIMlIUjbB39NdZvbbqMXNgIMNH5WrLU8grgozezzyvaQMggTyYfSyxkxSSzPbkeg4IiVTTJIygXQz2xPHZnmAgK3RC+Lcj0sCfg3E1ZoCP5I0U9IuSTskvSNpQtR6PcP27eslfVPSHEm7JZVKOjdcp7ukZyRtDffzuKSWUft5ONxPR0mPStoiaaektySNqCbG/yfpg3CfuyRNlfT1GOtZuP/jwvXLgX+Gy7pKuiOMe5ukPZIWSrpSUnrEPq4HHgrfvhPuM7KN//rwfc8Yx6/S1n+omMLlRZL+IekjSXslLZH0izDhH1bE/o+XNCWsnw2S7pTUPGrdytgHSfqdpDXAHmBsuLyDpLslrZa0L3y+W1L7iH08DKwM314XUT/jI+OpYew1+uxhvE9LWhuutyH8Gz2lJsdxh+ZnIO5IPAZ8C3iG4IszG/gO8IakM8zsxaj1TwV+CNxD8Av0POBBSfuAm4G3gWuA0cD3CL6gvh/juK+G218PdAYuAt6TNM7M5leuJOnXwC/C9a8FKoCvAU9LusjM7o7abxFwJnA/8EhE+VDgDOAfBO31mcDJwC1AL+AH4XrPAV0IztZuBhaF5UfSxh8zJklfCeMpBe4gqI9xwP8SNKN9o4b7Hwl8Pdz/o8AE4KfAYEknmFlF1Pp/BXaHxzRgvaTWwH+Ao4AHgVnACOBHwJckjQnPmv4MzAF+H8b+XLjPRcShpp89TF5vh5vdS5C8OhDU6dHAS/Ec18VgZv7wxyEfwDkEXxbnRJR9LSy7IGrdDGAGUAYoLOsZrrsT6BGxbkeCJFEB/CxqP88B+4AWEWUPh/t5rnLfYfmocB+vRpSNDNe9OcbneR7YDrSMKLPwcXyM9ZtFHi+i/DGCNvsuMepqfIz1rw+X9YyxbAUwOaosZkxADsHF5/eAjKhll1Z3/BjHrNz/6VHld4blE2PEPjnGMW8Kl/04qvzCsPzGiLLKv4Xrq4nn4UPVSzyfHfhq+P6bif4/1Fgf3oTlauu7wA7g+bD5ooOkDkAbgmaWnkCfqG2eN7PKJgzMbDOwhODLP/ps4H2CX/o9Yxz7Ngu/IcL9zATeAI6X1CIs/g7Bl8cjkfGFMb4ItCT41Rqp2MzejD6Yme2uPJ6kLEntwv28RtAMXBQjxroSK6YTCK4lPAS0ifpsL4frfLmG+19iZs9Hld0SPn8txvp/MLMDUWVfI+hMcV9U+Z+Bj6rZT23F89k/CZ9PltSqDmNwIW/CcrU1gOBLeOMh1skDSiLeL4+xzjZgvZntjVEO0J6qYjV5LCT44ugBLAjjE7D4MPFFKom1UtiufhVwNkEzjaJWaXuIYxypWDENCJ8fPMR20Z+tOlXq0szWS/qYoHmuJvEUAjOiE4uZHZC0hOBssK7U+LOb2buSHiU4K/yOpOnAm8CTZrawDmNqsjyBuNoSwa/Obx9inflR76vronmorpvRX9Y1XU8EZyAnH2L/C6Le76pmvd8BPwGeJGiu2QTsJ/hivJWad0Y51P011f1fjBVT5Wf9OcE1hVjWHWFM1dV7dXXUUOL67GY2SdLtwFeAzwOXAb+QdImZ3VWfgTYFnkBcbS0F+gJTzKy8gY89AJgSo+wgn/byWQqcBKwys7gu0sZwFvCemU2MLJR0VIx1D5UkKruutiPiHhFJOQQX30trGM/S8HlnrCa3OA2MLpDUBWhN7DPGWJYD/SRlRJ6FhGdufePYT03E/dkt6FgxH7hNUhtgKnCLpLsjm0Jd/PwaiKutRwn+fn4Ta6Gkmjah1MYVkv77C1nSSOB44K2IZPZY+HxzZFfbiG06xXG8g0T9Ig+7uV4aY93K47eLsayy+ef4qPJLie//4msEZ0FXSapyHEnNFNUF+hD6STo9quzK8Pn5Gu7jeYIOEdE95s4Py/9Rw/3URI0/e3it6jP1amYfE3TwyCW4IO+OgJ+BuFoxs2ckPQRcFH6B/4vggmk+wcXpo4jdhl4XegCvSXqR4Jf7RQRdS38eEd90SdcBNwBzJD1N0LTRhaDX1leArBoe7xngB5KeJGhDzyPoZrwlxrrTCToF/EJSW4KeZ2VmNjXcdjHwv2EX0zKCZpWxBHVXI2a2U9LZBF/cSyQ9SHD20gboT9Dl+GsEPaYOZx7wuKT7CX7dTyDo1vsuQZNdTdxG0HX27vBvYTZBN97zCDpJ3FbD/RxWnJ/9bOBSSZVdfvcDXwROBJ4ys911FVdT5QnE1ZqZfU/SOwT3PVxN8IW8geA+gKvr8dAnEVyXuIGgi+0U4OdmNjcqvv+VNJPgvoZLgOYEv17nAxfHcbyfEfQ4+yZwGrCaoMdR5UXZyGOukvQ9gl/xfyLoSfYIMNXMDko6DfgjwTWVfcDrBF9q/44jHszsNUmjCS7uf5fgl/42gntOfgfMPcTmkWaFn+8mgnt0tgN3AddY1XtAqovlE0nHEPx7fBU4l6Bzxb3AdVbHd87H8dknEySyUwl+OBwkSNqXE3xGd4TkTYAuVYR3KU8ys5peWHeHoGDQzEfM7JxEx+JSk18Dcc45VyueQJxzztWKJxDnnHO14tdAnHPO1UqT6oXVoUMH69mzZ6LDaBA7d+6kefPmh1+xCfE6qcrrpCqvk6pmzpz5kZl1jC5vUgmkZ8+ezJgxI9FhNIjJkyczfvz4RIeRVLxOqvI6qcrrpCpJK2OV+zUQ55xzteIJxDnnXK2kdAKRdFI4lWWppKsSHY9zzjUlKZtAwgHy7iYYrnsg8C1JVUYWdc45Vz9SNoEAY4BSM1tuZvuAJwjGKXLOOdcAUrkXVjeCQe0qrQGOjl5J0gUEg/2Rl5fH5MmTGyS4RCsvL28yn7WmvE6q8jqpyuuk5lI5gcQaUK/KXZFmdh/hXM1FRUXWVLrneVfEqrxOqvI6qcrrpOZSuQlrDVAQ8T6fmk/jGZe71qzhmU2b2F9Ro9GtnXOuSUjlM5DpQB9JhcBaYCKHnp+7VirM+PP69czfuZMuWVmc36UL53fpQn6OT2bmnGvaUvYMJJx7+SKCKS4XEcwwtqCuj5MmMaeoiH8OHsyIFi24ceVKek6Zwhnz5/PG1q1U+FhizrkmKpXPQDCzl4GX6/s46RKndujAqR06ULZ7N39et46/bNjAPz76iD7NmvHDrl05p3Nn2mVm1ncozjmXNFL2DCRRCps145bevVkzbhyPDxhAx8xMLlu2jG4ffsi5ixczffv2RIfonHMNIqXPQBIpOy2N7+Tl8Z28PIrLy/nT2rU8vnEjD2/YwKgWLfhxt25M7NSJ3PT0RIfqnHP1ws9A6sCwFi24t18/1n3uc9zVpw97Kio4b8kSun34IZeWlrJk165Eh+icc3XOE0gdapWRwYXdujFv9GjeHT6ck9q14+61a+k/bRonFBfz8pYtftHdOddoeBNWPZDEsW3acGybNmzct4+/rF/PPWvXcsq8efTPzeXS/HzOysujmTdvOedSmJ+B1LO8rCyu6dGDsrFj+euAAeSmpfGDkhIKPvyQa8vK2LB3b6JDdM65WvEE0kAy09L4dl4eM0aN4t3hw/l869bctHIlPaZM4dzFi5lbXp7oEJ1zLi7ehNXAIpu3lu7axZ1r1vDQhg08vGEDx7dty8/y8zmxXTvSFGuoL+ecSx5+BpJAfXJzuatvX1aPG8dvCgtZuHMnX5k3j8HTp3P/unXsPngw0SE651y1PIEkgXaZmVwVXid5fMAActLSuKCkhO5TpvArv07inEtSnkCSSFZ4c+LMUaOYPHw4n2vVil+HY2/9qKSE5bt3JzpE55z7L08gSUgSX2zThheGDGHxmDFM6tyZB9evp+/UqZy1aBELdu5MdIjOOecJJNn1zc3lz/36sXzsWC7Oz+cfmzczePp0Tp83j6k+7pZzLoE8gaSIbtnZ3HHUUawcN47revTgvU8+YeysWRw3Zw5vbduG+R3uzrkG5gkkxbTPzOT6wkJWjh3L7b16sWjXLo4vLmbsrFm88NFHPlSKc67BeAJJUS0zMri8e3eWH3009/bty+b9+zl9/nyGTp/O4xs24B2AnXP1zRNIistJT+cHXbtSMmYMfx0wAEmctXgxZwH3rl3LXp/H3TlXTzyBNBIZ4VApxUVFvDB4MG2AHy1dSp+pU/nzunXs80TinKtjnkAamTSJr3bowN3A60OHkp+dzQ9LSug7dSoPrFvHfk8kzrk64gmkkRJwQrt2/HvECF4ZMoS8rCzOLymh37RpPLh+vScS59wR8wTSyEnipPbtmTJyJC8NGUL7zEzOW7KEAdOm8ciGDRzwROKcqyVPIE2EJL7Svj3TRo7kxcGDaZWRwTmLFzOwsteWd/91zsXJE0gTI4n/6dCBmaNG8fzgweSmpXHW4sUMmjaNv23c6InEOVdjnkCaKEmc1qEDs4qKeHbQoGAgx0WLGDJ9Os9u3ux3tjvnDssTSBOXJnFGx47MKSri6YEDEfD1BQs4etYs3t62LdHhOeeSmCcQBwSJ5OudOjF39Gge7t+fjfv2cVxxMV8uLmbmjh2JDs85l4Q8gbjPSJeY1LkzS8aM4fe9ezNrxw6KZs7k/y1YQMmuXYkOzzmXRDyBuJhy0tO5pKCA5WPH8qsePXhpyxYGTpvGD5csYZ3PkOicwxOIO4xWGRncUFjIsrFj+XG3bjy4YQNHTZ3KVcuWsW3//kSH55xLIE8grkbysrL4Y58+LBkzhjM7duS21avpNXUqt65axa6DPvavc02RJxAXl8JmzXhswADmFBVxTKtWXLV8OX2nTuXpTZu8669zTUxSJBBJBZLekbRI0gJJF4fl7SS9IWlp+Nw2YpurJZVKWiLpxMRF3zQNbdGCfw0dyrvDh9MxK4tvLlzISXPnstQvtDvXZCRFAgEOAJeZ2QBgLHChpIHAVcBbZtYHeCt8T7hsIjAIOAm4R1J6QiJv4o5t04bpI0dy51FHMWX7dgZPn861ZWXs9mYt5xq9pEggZrbezGaFr3cAi4BuwGnAI+FqjwCnh69PA54ws71mVgaUAmMaNGj3Xxlpafw0P5/FY8bwjY4d+fXKlQyaPp1/ffRRokNzztUjJVu7taSewHvAYGCVmbWJWLbNzNpKuguYYmaPh+V/AV4xs2di7O8C4AKAvLy8UU888UT9f4gkUF5eTosWLRJy7DnAH4CVwDHARUDnhETyWYmsk2TldVKV10lVEyZMmGlmRdHlGYkIpjqSWgDPApeY2XZJ1a4aoyxmJjSz+4D7AIqKimz8+PF1EGnymzx5Mon6rOOBH1dU8Ic1a7hhxQq+B/yyRw8uKyggOy1xJ72JrJNk5XVSlddJzSVFExaApEyC5PFXM3suLN4oqUu4vAuwKSxfAxREbJ4PrGuoWN3hZaWlcUX37iweM4avtGvHL8rKGDZ9Om/5+FrONRpJkUAUnGr8BVhkZr+LWPQiMCl8PQl4IaJ8oqRsSYVAH2BaQ8Xraq4gJ4dnBg/mlSFDOGDG8cXFXF5a6nO0O9cIJEsT1jHAWcA8SXPCsmuAW4CnJJ0HrAK+AWBmCyQ9BSwk6MF1oZl5t58kdlL79sxr04bLly3jjjVreP+TT/j7wIH0atYs0aE552opKRKImX1A7OsaAMdVs81NwE31FpSrc83S07m7b1++1LYt5y1ezIgZM3igXz++0alTokNzztVCUjRhuablzHD+kQG5uXxz4UJ+VFLi9404l4I8gbiE6NmsGe+PGMEVBQXcu24dR8+axaKdOxMdlnMuDjVOIJLGSrpe0quS5obDi3wo6WFJ50YOM+JcTWSmpXFr7968PGQI6/fto2jmTB5av97H1HIuRRw2gUiaJGke8B/gEiAXWApMBbYBRwMPAGvDZFJYf+G6xujk9u0pLiri6Fat+N6SJXx30SIfKt65FHDIi+iSioFOwKPA2cAci/HzUFJr4FTgO8ACSeea2ZP1EK9rpLpmZ/PGsGHcvHIlN6xYwdsff8w9ffrwtY4dEx2ac64ahzsDeQgoNLMrzWx2rOQBYGafmNlfzewrwDjg4zqO0zUB6RLX9uzJ9FGj6JyVxRkLFvDNBQvYuG9fokNzzsVwyARiZn8wsz3x7NDMis3stSMLyzVlI1q2ZNrIkdxUWMgLH33EwGnTeHzDBr824lyS8V5YLillpqVxTY8ezCkqol9uLmctXsyp8+axek9cv2ecc/UorhsJJWUBJwP9gJyoxWZmN9ZVYM4BDGjenPdHjOCutWu5ZvlyBk2fzm29enFB166kVT/YpnOuAdQ4gUjqCnwA9CQY+bbyf29ku4InEFfn0iUuzs/nq+3bc/6SJfxo6VKe2LSJB/r146jc3ESH51yTFU8T1u3AZqA7QfI4GuhFMJxIafjauXpT2KwZbwwbxgP9+jGnvJyhM2ZwYUkJC/0GROcSIp4E8gXgDj4dNr3CzFaY2a+AZ4A/1nVwzkWTxHldurBgzBi+2bEjD6xfz6Dp0/nSnDk8u3kzB3yUX+caTDwJpD2wzswqgJ1A5J3nbxPMI+Rcg+iWnc3DAwawZtw4flNYyLLdu/n6ggX0nDKFX69Y4V1/nWsA8SSQNUCH8PUy4MsRy8YA3j3GNbiOWVlc1aMHy8eO5fnBgxnYvDnXrlhBwYcf8u2FC/nPJ59491/n6kk8vbDeAb4IPA/8Gbhb0nBgP3BiWOZcQqRLnNahA6d16MCSXbu4Z+1aHt6wgb9v2sTwFi24sGtX8hMdpHONTDxnIL8E/gRgZn8CLiYYF6sLcBtwWZ1H51wt9MvN5c4+fVg7bhz39u3LATPOLynhG8BlpaWU7tqV6BCdaxRqlEDCM43xQHdJ2QBm9n9m9nkzG2lm18R7x7pz9a1FRgY/6NqVuUVFvDt8OEXAH9eupc+0aXxl7lxe2rKFCm/ecq7WDplAJLWR9DYwE3gSeA0olTS4IYJzri5I4tg2bbgOWDl2LNf16MGc8nJOnTePgdOmcf+6dezxCa2ci9vhzkB+RXC/xw0Eo+3+FEgH7qrnuJyrF12zs7m+sJCVY8fytwEDaJ6ezgUlJfScMoWbVq5kqw8j71yNHe4i+inAjWZ2S2WBpBLgVUktzWxHvUbnXD3JTEvjW3l5TOzUiXc+/pjbV6/ml2Vl/GblSs7r0oVL8/Pp2axZosN0Lqkd7gykJ/DvqLIPCO5E714fATnXkCTxpbZteWXoUOYWFXFmx47cs24dR02dyrcWLmTmDv+N5Fx1DpdAMoG9UWWVd2hl1304ziXOkBYteGTAAMqOPppLCwp4acsWimbO5Etz5vDKli1+P4lzUWpyH8j/RF00TyMYQPGrYe+s/zKzB+swNucSIj8nh9t79+aXPXpw37p13LlmDV+ZN4/BzZtzeUEB3+rUiaw0nwnBuZokkF9UU/6rqPcGeAJxjUbrjAx+3r07F+fn88SmTdy+ejXnLF7MNcuXc0l+Pj/q2pUWGXHNiOBco3K4v/7CBonCuSSWlZbG2Z07c1ZeHq9t3crtq1dzxfLl3LpqFZcVFHBht2608kTimqBD/tWb2cqGCsS5ZCeJk9q356T27ZnyySfcuHIl15SVcfvq1fwsP5+f5OfT2hOJa0K8Ide5WhjbujUvDR3K9JEj+Xzr1ly7YgU9PvyQ68vK2Ob3krgmIq4EImmSpFclLZS0POqxrL6CdC5ZFbVqxYtDhjBr1Ci+1LYtN6xcSc8pU7i2rIwtnkhcIxfPlLbXEtyRPh+YQ9Xuvc41WSNatuS5wYMpLi/n1ytX8uuVK/nDmjVc1K0bP8vPp2NWVqJDdK7OxdNgex5wp5ldWl/BOJfqhrVowdODBjG/vJybVq3i1lWr+L81a/hxt25cVlBAnicS14jEOyPhP+srEOcak8EtWvD3gQNZMHo0p3fowB2rV1M4ZQo/X7bMm7ZcoxFPAnkXGFZfgQBISpc0W9K/wvftJL0haWn43DZi3asllUpaIunE+ozLudoa0Lw5jw8cyKIxY/h6x47csXo1vcJpd8sPHEh0eM4dkcMN555W+QAuAc6VdLakDpHLItY5UhcDiyLeXwW8ZWZ9gLfC90gaCEwEBgEnAfdISq+D4ztXL/rm5vLogAHMLSpiQps2XLtiBb2nTuX/1qxhb0VFosNzrlYO96V/gGDK2v3AEmAw8BCwMaK88rGvmn3UiKR8gtF/H4goPg14JHz9CHB6RPkTZrbXzMqAUoJ52Z1LaoNbtOD5IUP4cMQIBjZvzk9LS+k/bRqPbdjAQR9ry6UYHWqAOEnXEwxRUiNmdkOtA5GeAX4DtAQuN7NTJX1sZm0i1tlmZm0l3QVMMbPHw/K/AK+Y2TMx9nsBcAFAXl7eqCeeeKK2IaaU8vJyWrRokegwkkqy1YkBM4D7gaUEQ1+fBxxDMNx1Q0i2OkkGXidVTZgwYaaZFUWXH+5O9OvrLaIIkk4FNpnZTEnja7JJjLKYic7M7gPuAygqKrLx42uy+9Q3efJkmspnralkrJMJwGVmPLt5M78sK+Pa3bsZ26oVvyksZHzbtofd/kglY50kmtdJzSXLnejHEIzuuwJ4AviSpMeBjZK6AITPm8L11wAFEdvnA+saLlzn6k6axDc6dWLB6NHc37cvq/fsYUJxMScVFzPL5yNxSexwF9EvlZQTzw4ljZR0UjzbmNnVZpZvZj0JLo6/bWbfBV4EJoWrTQJeCF+/CEyUlC2pEOgDTIvnmM4lm4y0NL7ftStLjz6a3/buzfQdOxg1cyZnLVrE6j17Eh2ec1Uc7gzkbGCFpFskVduFV1JbSWdJep1gxsJWdRTfLcAJkpYCJ4TvMbMFwFPAQuBV4EIzO1hHx3QuoZqlp3NZQQHLx47l6u7deXrTJvpNm8avysq8669LKoe7E30kcBZwGXCFpO3APGAzwVAmbYFeQO/w/ZPAQDNbUduAzGwyMDl8vQU4rpr1bgJuqu1xnEt2rTMyuLlXL37QtStXLV/OjStX8sD69dxUWMikzp1JU0NdancutkOegVjgUTMbBowDfg/sIEgaIwh6TL0PfA/oambnHknycM5V1SMnh78PHMh/RoygR04O31uyhKKZM5m8bVuiQ3NNXI3HwjKzqcDUeozFOXcI41q35j8jRvDkpk1cuXw5E4qLOb1DB27r1Ys+ubmJDs81QcnSC8s5VwOSmJiXx+IxY7i5sJA3t21j0PTp/Ky01OchcQ3OE4hzKahZejpX9+jB0jFjOKdzZ+5cs4ajwqFRDvjQKK6BeAJxLoV1zs7mvn79mF1UxIgWLfhpaSkjZ87k/Y8/TnRorgnwBOJcIzC0RQveGDaMfwwaxCcHDnDsnDmcvWgRG/b6vG+u/ngCca6RkMTpHTuyaMwYftG9O0+G94/c6c1arp54AnGukclNT+fXvXoxf/RoPte6NZd4s5arJ3ElEEnNJf1U0jOS3pHUJyyfKKl//YTonKuNPrm5vDxkyGeatc7yZi1Xh2qcQCQVAHOB2wnGnjqW4EZCCAYVvbzOo3POHZHIZq1f9ujBU2Gz1h9Wr/ZmLXfE4jkDuYNguJI+wCg+O6T6uwQJxTmXhHLT07mxsPC/zVqXLlvGyJkzmZ/owFxKiyeBnABcZ2arqDr3xlqgW51F5ZyrF9HNWj8FLiwp4RMfpNHVQjwJJItgHKxYWhNMa+ucS3KVzVoLRo/mTODedesYOG0az2/enOjQXIqJJ4HMBc6sZtnJwMwjD8c511BaZGRwITBl5Eg6ZmbytQULOGP+fNb6RXZXQ/EkkNuB8yTdz6fXOwZKuoFgKufb6zo451z9G92qFdNHjeLWXr14ZetWBkybxj1r11JhMWeJdu6/apxAzOw54MfAN4A3w+JHgUuAi8zs1TqPzjnXIDLT0riie3fmjx7N0a1aceHSpXx+9mwW7NyZ6NBcEounG29r4CGCi+UnAt8laLrKN7P76ic851xD6t2sGa8PHcqj/ftTsmsXI2bM4NqyMvYc9Ak/XVU1SiCSMoAtwJfNbKeZvWlmfzOz18ysugvrzrkUJImzOndm0ZgxTOzUiV+vXMnwGTP48JNPEh2aSzI1SiBmdgDYCPjPEOeaiI5ZWTw6YACvDR3K7ooKPj97Nj9ftozdfjbiQvFcRH8c+H59BeKcS05fbteOeaNH8/0uXfjt6tWMnDmTKX424ohjSltgBfBtSdOBF4D1RN1QaGYP1l1ozrlk0Sojgz/368fXO3bkvCVLOGb2bC4vKOCGnj3JSU9PdHguQeJJIHeHz90IhjKJZoAnEOcasRPatWP+6NFcvmwZt61ezT+3bOHh/v0Z06pVokNzCRBPE1bhYR696jw651zSaZWRwX39+vHq0KHsOHiQcbNmcdWyZd5Tqwmq8RmIma2sz0Ccc6nlxPBs5LLSUm6NOBsZ7WcjTYZPKOWcq7XWGRk80L8/rwwZwvbwbOR/V6zwoeKbiBqfgUgqo+oovJ9hZt6M5VwTdFL79swrKuInpaVct2IFr2zdymP9+3NUbm6iQ3P1KJ4zkHdjPOYDrcL9TK7r4JxzqaNNZiaPDRjA3wcMYPGuXQyfMYMH1q3DfEytRiueayDnxCqX1AZ4lU/Hx3LONWET8/I4pnVrzlm8mPNLSvjXli3c368fHbOyEh2aq2NHfA3EzD4mGIn3V0ccjXOuUSjIyeGNYcO4o3dvXtm6lSHTp/PKli2JDsvVsbq6iL4HyK+jfTnnGoE0iZ8VFDB91Cg6ZmXxlXnzuLCkhF3e3bfROKIEIilD0nDgemBBXQTknGtchrZowfSRI/lZfj73rFvHyBkzmL3Dx2BtDOIZzr1C0sHIB7CXYCbCo4BLjyQQSW0kPSNpsaRFksZJaifpDUlLw+e2EetfLalU0hJJJx7JsZ1z9SsnPZ07jjqKN4cNo/zgQcbOmsUf16zxC+wpLp6hTP6Xqt149wArgVfM7EhHV7sTeNXMvi4pC8gFrgHeMrNbJF0FXAVcKWkgMBEYBHQF3pTU18z83Ni5JHZc27bMKSri3CVLuLi0lDe3beOh/v1pn5mZ6NBcLcTTC+v6+gpCUiuCaXLPCY+1D9gn6TRgfLjaIwRdha8ETgOeMLO9QJmkUmAM8GF9xeicqxsdsrJ4cfBg/rh2LT9ftozhM2bwtwED+EKbNokOzcXpSK+BDJR0pqSuRxhHL2Az8JCk2ZIekNQcyDOz9QDhc6dw/W7A6ojt14RlzrkUIImL8/P5cORIctLSGD9nDjeuWMFBb9JKKfHciX4XkGFmPwzfnwE8CaQD2yWdYGbTjyCOkcBPzGyqpDsJmquqDSdGWcy/PEkXABcA5OXlMXny5FqGmFrKy8ubzGetKa+TqpKhTu4Efg/8asUK/rFiBdcAHRIYTzLUScowsxo9gGXA2RHv5wH/AIYArwP/qum+Yuy7M7Ai4v0XgJeAJUCXsKwLsCR8fTVwdcT6rwHjDnecUaNGWVPxzjvvJDqEpON1UlWy1ElFRYU9tG6d5b77rnX44AN76aOPEhZLstRJMgFmWIzv1HiasDoTTCqFpHyCC9i/MbN5wB+B0bVLYWBmG4DVkvqFRccBC4EXgUlh2SSCiawIyydKypZUCPQBptX2+M65xJLEOV26MHPUKLpmZXHKvHlcsWyZD8qY5OLphbUbaBG+/iKwHZgRvi8HWh5hLD8B/hr2wFoOnEtwjeYpSecBq4BvAJjZAklPESSZA8CF5j2wnEt5/Zs3Z+rIkVy6bBm3r17NlO3beXLgQLpkZyc6NBdDPAlkFnChpFXAhcAbZlb586CQYIrbWjOzOUBRjEXHVbP+TcBNR3JM51zyyUlP5099+3JMq1b8oKSEETNm8MTAgYxv2/bwG7sGFU8T1i+AsUAx0A+4MWLZ6XgTknOuDn23c2emjRpFm4wMjisu5tZVq6jwXlpJpcYJxIIeVt0J7rcoNLO5EYvvA66r49icc03coObNmT5qFF/v2JGrli/na/Pns23//kSH5UJx3QdiZjvNbKaZba8sk9TezF4ys5K6D88519S1zMjgiYED+eNRR/HK1q2MmjmTWT6WVlKIZyys8yX9POL9EElrgE2SZkjqXC8ROueaPEn8JD+f94YPZ78Zn5s1i/t9sqqEi+cM5CcEPbEq/Q74GLgEaE0wVpZzztWbsa1bM2vUKI5t04YLSkq4oKSEvd7VN2Hi6YXVHVgMIKk1QVfe083sZUlbgN/UQ3zOOfcZHbOyeGXoUK4rK+OmVauYV17Os4MH0827+ja4eM5A0oHKVP95gqFDJofvV/PpOFXOOVev0iV+3asXzw4axPydOymaOZN/f3KkA4K7eMWTQJYCp4SvJwL/MbNd4fuuwNa6DMw55w7njI4dmTpqFC3S05kwZw73rl3r10UaUDwJ5LfAJZI+Ar4N/F/EsgnA3JhbOedcPRrUvDnTRo7k+LZt+dHSpX5dpAHFMx/I38K70I8GppvZexGLNxKMT+Wccw2ubWYm/xwyhF+VlXHzqlXM37mTZwcNoqtfF6lX8VxEx8w+AD6IUe43ETrnEipd4qZevRjZsiWTFi1i1MyZPDtoEJ9r3TrRoTVacd1IKKm5pJ+Gc5e/I6lPWD5RUv/6CdE552ruzKjrIo9u2JDokBqteG4kLCC4znE7wfDpx/LpCLwTgMvrPDrnnKuFQeGovp9v3ZpJixdz1bJlPtthPYjnDOQOYC9B8hjFZ2cFfJcgoTjnXFJol5nJq0OH8sOuXbl19WrOmD+fHQcOJDqsRiWeBHICcJ2ZraLq9LFr8TnJnXNJJjMtjT/17ctdffrw0pYtHDN7Niv37El0WI1GPAkkC6huBLPWgA+R6ZxLShd268YrQ4eyas8eRvtNh3UmngQyFzizmmUnAzOPPBznnKsfJ7Rrx9RwfpEvzZnDI35x/YjF0433duAZSQB/C8sGSjoNOA/4ah3H5pxzdapfbi5TRo7kmwsWcM7ixSzauZObe/UiTTr8xq6KeCaUeg74McG85G+GxY8SjMZ7kZm9WufROedcHWuXmckrERfXJy5cyO6DBxMdVkqq8RlIOALvQ8BjwDiCwRO3EIyJ5bO7OOdSRmZaGvf06UPvnBx+vnw5a/bu5YXBg+mYlZXo0FJKjRKIpAyCZPE1M/snn56BOOdcSpLE5d270zMnh7MWL2bcrFm8PHRoosNKKTVqwjKzAwTjXfl5nnOuUfl6p068M2wY2w8eZNysWT4qbBzi6YX1OPD9+grEOecSZWzr1kwZOZKOmZlcDvxt48ZEh5QS4umFtQL4tqTpwAvAeqJuKDSzB+suNOecazi9mjXjPyNHMuHf/+Y7ixZRtmcP13TvjryHVrXiSSB3h8/dCIYyiWaAJxDnXMpql5nJ7cCjnTrxy7IyVu7Zwz19+pCRFte4s01GPAmksN6icM65JJEFPDZgAD1ycrh51So27tvH3wcOJDc9PdGhJZ14EshOoNzMfCAZ51yjpnBukW7Z2Vy0dCnHFxfzzyFDaJ+ZmejQksohz8skpUu6XtLHBL2wtkt6VlKbhgjOOecS6cfduvH0oEHM2rGDz/tAjFUcrmHvh8CvgFkEc6K/AJwG/L6e43LOuaRwZseOvDFsGBv27Qu6+ZaXJzqkpHG4BHI+cL+ZfcnMrjSzbwAXAt+V5LdsOueahC+0acP7w4eTBnxh9mwmb9uW6JCSwuESSC/g6aiyJ4F0oEe9ROScc0locIsWfDhyJPnZ2Zw4dy5PbdqU6JAS7nAJpAWwPaqsctyrljjnXBNSkJPD+yNGMKZVKyYuXMif161LdEgJVZPOzd0k9ap8EJyVVCkPl9WapEslLZA0X9LfJeVIaifpDUlLw+e2EetfLalU0hJJJx7JsZ1zrqbaZWby+tChnNK+PT8sKeE3K1diTXS+9Zp0432mmvLnY5TVqqO0pG7AT4GBZrZb0lPARGAg8JaZ3SLpKuAq4EpJA8Plg4CuwJuS+pqZj9XlnKt3zdLTeW7QIM5dvJhrysrYsn8/t/fu3eTuWj9cAjm3QaIIZADNJO0HcoF1wNXA+HD5I8Bk4EqCnmBPmNleoExSKTAG+LAB43XONWGZaWk8OmAAbTMzuWPNGrYdOMCf+/ZtUnetHzKBmNkjDRGEma2V9FtgFbAbeN3MXpeUZ2brw3XWS+oUbtINmBKxizVhWRWSLgAuAMjLy2Py5Mn19CmSS3l5eZP5rDXldVKV10lV8dbJGQQXih/csIGlGzbwS4K72ZuCeO5ErzfhtY3TCIZL+Rh4WtJ3D7VJjLKYjZBmdh9wH0BRUZGNHz/+iGJNFZMnT6apfNaa8jqpyuukqtrUyQRg5Jo1XFJayq1t2vD84MG0zEiKr9d6lSznWscDZWa22cz2A88BnwM2SuoCED5X9ptbAxREbJ9P0OTlnHMJcXF+Po/078+7H3/MccXFbNm/P9Eh1btkSSCrgLGSchVchToOWAS8CEwK15lEcCc8YflESdmSCoE+wLQGjtk55z7j7M6deXbwYIrLy5kwZw6b9u1LdEj1KikSiJlNJejtNQuYRxDXfcAtwAmSlgInhO8xswXAU8BC4FXgQu+B5ZxLBqd16MC/hgyhdPdujp09m7V79yY6pHqTFAkEwMyuM7P+ZjbYzM4ys71mtsXMjjOzPuHz1oj1bzKz3mbWz8xeSWTszjkX6YR27Xht6FDW7dvHsbNns2L37kSHVC+SJoE451xj8oU2bXhz2DC2HjjAF+bMoWTXrkSHVOc8gTjnXD0Z06oV7wwbxp6KCo6dPZv5jWwkX08gzjlXj4a3bMl7w4eTJjF+zhxm7dhx+I1ShCcQ55yrZwOaN+e94cNpnp7Ol+bMYer26DFqU5MnEOecawBH5eby/ogRtM/M5MvFxUz55JNEh3TEPIE451wD6Z6Tw+Thw+mYmcmX587lPymeRDyBOOdcAyrIyeHdESPonJXFiXPn8sHHHyc6pFrzBOKccw2sW3Y2k4cPp2tWFifNncv7KZpEPIE451wCdA2TSEFODifPncu7KZhEPIE451yCdMnO5p1hw+iRk8NX5s7lnW3bEh1SXDyBOOdcAnXOzuad4cMpzMnhlHnzUupMxBOIc84lWKesLN4ePpyeOTmcMncu/06R3lmeQJxzLgl0ysrirWHD6Jadzclz56bEfSKeQJxzLkl0yc7m7eHD6ZSZyYlz5zIjye9Y9wTinHNJpFt4TaR9ZiYnzJ3L7CQeO8sTiHPOJZmCnBzeHjaMVunpnFBczLwkHcXXE4hzziWhns2a8fbw4eSkpXFccTELd+5MdEhVeAJxzrkk1btZM94ZPpwMieOLi1mWZDMbegJxzrkk1ic3lzeHDWNfRQXHzZnD6j17Eh3Sf3kCcc65JDeweXNeHzaMbQcOcHxxMRv37Ut0SIAnEOecSwkjW7bk5aFDWbN3LycUF7N1//5Eh+QJxDnnUsUxrVvzwuDBLNm1i5PmzmX7gQMJjccTiHPOpZDj27XjmUGDmF1ezv/Mm8eugwcTFosnEOecSzH/06EDj/Xvz/uffMIZ8+ezr6IiIXF4AnHOuRQ0MS+P+/v147Vt2zh70SIOmjV4DBkNfkTnnHN14rwuXdi6fz9XLF9Ou6VLubtPHyQ12PE9gTjnXAr7effufLR/P7etXk37zExuLCxkb0UFz27ezO2rVrF49272VlSQnZZG/2bNuKJ7d87o2JHstCNvgPIE4pxzKe6WXr3YeuAAv165kpKdO3l12zYM2BFxgX1PRQVzdu7kgpISflhSwm979+b8rl2P6Lh+DcQ551KcJO7t25d+zZrx1Ecfsf3gwc8kj0jlBw+y/eBBLikt5dqysiM6ricQ55xrBB5cv57Ve/fWeP1dFRX8bvVq7l+3rtbH9ATinHMpbm9FBZcvW8auOLvz7gq3q2034AZNIJIelLRJ0vyIsnaS3pC0NHxuG7HsakmlkpZIOjGifJSkeeGyP6ohux0451ySeW7zZmrbibci3L42GvoM5GHgpKiyq4C3zKwP8Fb4HkkDgYnAoHCbeySlh9v8CbgA6BM+ovfpnHNNxm2rVlV7zeNwyg8e5NbVq2u1bYMmEDN7D9gaVXwa8Ej4+hHg9IjyJ8xsr5mVAaXAGEldgFZm9qGZGfBoxDbOOdfkLD7CeUKW7NpVq+2SoRtvnpmtBzCz9ZI6heXdgCkR660Jy/aHr6PLY5J0AcHZCnl5eUyePLnuIk9i5eXlTeaz1pTXSVVeJ1WlYp3U/NJ5bHsqKmr1mZMhgVQn1nUNO0R5TGZ2H3AfQFFRkY0fP75Ogkt2kydPpql81pryOqnK66SqVKyT7PfeY88RjIeVk5bG+GOPjXu7ZOiFtTFsliJ83hSWrwEKItbLB9aF5fkxyp1zrknq36zZEW3fLze3VtslQwJ5EZgUvp4EvBBRPlFStqRCgovl08Lmrh2Sxoa9r86O2MY555qcK7p3p2V6+uFXjKFlejpXFhQcfsUYGrob79+BD4F+ktZIOg+4BThB0lLghPA9ZrYAeApYCLwKXGhmld0MfgQ8QHBhfRnwSkN+DuecSyZndOwYs22/JhRuXxsNeg3EzL5VzaLjqln/JuCmGOUzgMF1GJpzzqWs7LQ0ftu7N5eUlsZ1M2FuuF1WLQdWTIYmLOecc0fo/K5d+VlBAbk1TAa5aWn8rKDgiAZUTOZeWM455+JwY2Eh3bOzuXzZMioIbhKM1iI9nTSok9F4PYE451wjcn7XrpzduTPPbd7MbatXs2TXLvZUVJCTlka/3FyuLCjgjI4da91sFckTiHPONTLZaWl8Ky+Pb+Xl1etxZAmYRzdRJG0GViY6jgbSAfgo0UEkGa+TqrxOqvI6qaqHmVXpqtWkEkhTImmGmRUlOo5k4nVSlddJVV4nNee9sJxzztWKJxDnnHO14gmk8bov0QEkIa+TqrxOqvI6qSG/BuKcc65W/AzEOedcrXgCcc45VyueQFKQpAJJ70haJGmBpIvD8naS3pC0NHxuG7HN1ZJKJS2RdGLioq8fh6iT2yUtljRX0j8ktYnYpknWScTyyyWZpA4RZU22TiT9JPzcCyTdFlHeqOvkiJiZP1LsAXQBRoavWwIlwEDgNuCqsPwq4Nbw9UCgGMgGCgmGwE9P9OdooDr5MpARlt/qdcLA8H0B8BrBjbUdmnqdABOAN4HscFmnplInR/LwM5AUZGbrzWxW+HoHsIhgXvjTgEfC1R4BTg9fnwY8YWZ7zayMYB6VMQ0adD2rrk7M7HUzOxCuNoVPZ7NssnUSLv49cAWfnQ66KdfJj4BbzGxvuKxyZtRGXydHwhNIipPUExgBTAXyLJixkfC5U7haN2B1xGZr+PSLpNGJqpNI3+PTyceabJ1I+iqw1syKo1ZrsnUC9AW+IGmqpHcljQ5Xa1J1Ei8fTDGFSWoBPAtcYmbbgxl+Y68ao6xR9t+OrpOI8l8AB4C/VhbF2LzR1wlBHfyCoGmvyqoxyhp9nYT/dzKAtsBYYDTwlKReNKE6qQ0/A0lRkjIJ/gP81cyeC4s3SuoSLu8CVJ6GryFo866UD6xrqFgbSjV1gqRJwKnAdyxs2Kbp1klvgrb8YkkrCD73LEmdabp1AsFnf84C04AKgkEVm0Sd1JYnkBSk4FTjL8AiM/tdxKIXgUnh60nACxHlEyVlSyoE+gDTGirehlBdnUg6CbgS+KqZ7YrYpEnWiZnNM7NOZtbTzHoSfEGONLMNNNE6CT0PfClcpy+QRTAib6OvkyPhTVip6RjgLGCepDlh2TXALQSn3ucBq4BvAJjZAklPAQsJmjAuNLOqU5Wlturq5I8EPWjeCJv4ppjZD5tynZjZy7FWbsp1AjwIPChpPrAPmBSerTaFOqk1H8rEOedcrXgTlnPOuVrxBOKcc65WPIE455yrFU8gzjnnasUTiHPOuVrxBOISLhwR9nCPFYmOs6FJ6inp+vCO6KQj6UxJGyXlRpStkPR41Hppkh6RVCHp/LDs0nCEZP8OSmF+H4hLBuOi3v+DYATU6yPK9jZYNMmjJ3Ad8AGwPLGhfFY49MfNwO1RN2jGWu8xgnuSzjGzR8NF9xLc4DkJeKiew3X1xBOISzgzmxL5XtJe4KPo8sZAUnbliK8pHsNpBAnuwUMcJxN4Avgq8G0ze6pymZntlvQocDmeQFKWnz66lCCpUNJfJW2WtFfSHElfi1rn+rC5q7+k1yTtlLRK0rnh8rMUTC5VHk4q1Dtq+xWSHpd0fjiB0B5JsyRNiBHPFyW9JWlHeJzXJA2OWmeypA8k/Y+k2WFi/HG47CJJH0raKuljSVMknRKx7XjgnfDtGxFNeePD5Sbp+qjj9QzLz4koe1jSGknjJP1H0m6CeWOQ1EHSnyStDet0saQLavhP8n3gVTPbGmuhpGzgOYIxyL4RmTwiPAEMlPS5Gh7TJRlPIC7pSSogGHJ7GHApwS/aWcCzCoYmj/Y08BLBfCgzCYaouJlgzoergHOBfsDfYmz7ReBnBCPWTiRoOntFUr+IeE4B3gLKge8C3yaYnOj9MNZIfQmGU/k/4MRwOwh+vT9A0LTz/4AZwL8knRwunwVcGL7+KUEz37iwPF6tCb6s/w6cDPxNUivg38ApBE2FpwD/BP4k6SeH2lmYHMYD71ezSjOCMaSOB043s+erWW8OsB04qcafxCWXRM9o5Q9/RD+AFcDjEe//AmwG2ket9wYwJ+L99QRDbZ8dUdaWYAyjLUCriPKfhuv2iDruPqB7RFlLYCvwWERZKfBWVCytCAbf+0NE2WSCUV2HH+bzphE0J78OvBBRPj6M8fgY2xhwfVRZz7D8nIiyh8Oy06LWvRbYA/SJKr8//BwZh4j36HCfJ1Tzb2fh43s1+Ld+H3g90X9z/qjdw89AXCo4CXgZ+ERSRuWDYErWYeGv6UiVk0ZhZtsIhrWfYhHzgwCLw+foM4YpZrYqYvsdBGcz4wAk9SEYEv2vUbHsAj4Ejo3a3wozmxP9gSSNkvQvSRsJEtx+4ASCM6O6dgD4V1TZSQRndWUx6rQ9wVSu1ekaPm+uZvk0giR0naQeh4ltc8T+XIrxBOJSQSfgbIIv2cjH7eHy9lHrb4t6v6+aMoCcqPKNMY6/kU9noauc5fEvMeI5NUYs66N3FjZzvQW0A34CfI5gEqNXY8RTFzZZ1RFkOxEku+jP8HS4PPpzRKqMsboL8UsJmq9aAm9JOtQMfrsJmrxcCvJeWC4VbCFo6ri1muV1OcFPXjVlayNiAbgaeDPGuvui3sca7vokgusS3zSzNZWFkfdT1MBegjkrIlX3pR8rhi0EZ2YXV7PNkkMcu7IO2la3gpkVS/oyQaJ8S9IXzSxWcm5HcLbiUpAnEJcKXiVoQlpgZrvr+VhjJRWY2WoASS0JLjC/FC5fQtDOP8jMbqnlMSoTxf7KAgWTGB1DMMFTpcpf+LF+oa8EBkeVnRJjveq8SnD2s8rMNh1u5SiVzX+9gP9Ut5KZzQg7BbwGvClpgplFJ4tCfIKmlOUJxKWCXxF8ybwn6S6CL/C2BF+gvczse3V4rI3A62EX2b0EN7s1B24EMDOTdCHwgqQs4CmCX9B5BE1Rq+yzM93F8ibBdYlHJd0BdAFuIJgELLJZuSRc73uStobxLAmvyzwB/FLBXO9TgC8A34rjc/6eoPfX+5J+T5AYmwP9gS+Y2WnVbWhmqyStBMYAj1e3XrjufySdSnAN63VJx4XXpZDUhqCX2m/jiNslEb8G4pJeeFG7iODu9JsJel/9iaDL7dt1fLh3gTvC4zxJ0N5/spmVRMTzMsH1g+YEXXFfI7i3ojPBhfRDMrMFwHeAHgTdXa8g6F78XtR6W4CLCLovvwtMB0aFi38D3BUufx4YQDDTXo2Y2ScECe9lgiT5GsFNgafx6f0nh/IkwTWfmhzrXYIu1QOAVyM6PZxC0OT3j5rG7ZKLz0joXEjBeFsfmNl3Ex1LsgtvwlwCjDezD2q5j1cIRhyoceJzycWbsJxzcTOzZZIeIjhzqtGZSCRJw4EJVL2O41KIN2E552rrWmB6nL3HKnUGzjWz0jqOyTUgb8JyzjlXK34G4pxzrlY8gTjnnKsVTyDOOedqxROIc865WvEE4pxzrlb+P6FKBxOb/cTBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot( col1.Tatm, col1.lev, 'c-', label='default' )\n",
    "ax.plot( col1.Ts, climlab.constants.ps, 'co', markersize=16 )\n",
    "ax.invert_yaxis()\n",
    "ax.set_xlabel('Temperature (K)', fontsize=16)\n",
    "ax.set_ylabel('Pressure (hPa)', fontsize=16 )\n",
    "ax.set_title('Temperature profiles', fontsize = 18)\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default this model has convective adjustment.  We can set the adjusted lapse rate by passing a parameter when we create the model.\n",
    "\n",
    "The model currently has no ozone (so there is no stratosphere). Not very realistic!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More reasonable-looking troposphere, but still no stratosphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### About the radiatively active gases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Band model is aware of three different absorbing gases: O3 (ozone), CO2, and H2O (water vapor). The abundances of these gases are stored in a dictionary of arrays as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'CO2': Field([0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038]),\n",
       " 'O3': Field([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       " 'H2O': Field([5.00000000e-06, 5.00000000e-06, 5.00000000e-06, 5.00000000e-06,\n",
       "        5.00000000e-06, 5.00000000e-06, 6.38590233e-06, 9.08848690e-06,\n",
       "        1.33273826e-05, 2.34389689e-05, 3.84220914e-05, 5.95564299e-05,\n",
       "        8.82144990e-05, 1.25843839e-04, 1.73951159e-04, 2.34088411e-04,\n",
       "        3.07840683e-04, 3.96815735e-04, 5.02635028e-04, 6.26926041e-04,\n",
       "        7.71315753e-04, 9.37425100e-04, 1.12686431e-03, 1.34122899e-03,\n",
       "        1.58209684e-03, 1.85102493e-03, 2.14954752e-03, 2.47917415e-03,\n",
       "        2.84138824e-03, 3.23764591e-03])}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col1.absorber_vmr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ozone and CO2 are both specified in the model. The default, as you see above, is zero ozone, and constant (well-mixed) CO2 at a volume mixing ratio of 3.8E-4 or 380 ppm.\n",
    "\n",
    "Water vapor is handled differently: it is determined by the model at each timestep. We make the following assumptions, following a classic paper on radiative-convective equilibrium by Manabe and Wetherald (J. Atmos. Sci. 1967):\n",
    "\n",
    "- the relative humidity just above the surface is fixed at 77% (can be changed of course... see the parameter `col1.relative_humidity`\n",
    "- water vapor drops off linearly with pressure\n",
    "- there is a small specified amount of water vapor in the stratosphere."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Putting in some ozone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.Dataset&gt;\n",
       "Dimensions:    (lat: 64, lev: 59, lon: 128, time: 12)\n",
       "Coordinates:\n",
       "  * lev        (lev) float64 0.2842 0.3253 0.3719 ... 849.5 959.0 1.004e+03\n",
       "  * lon        (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2\n",
       "  * lat        (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 82.31 85.1 87.86\n",
       "  * time       (time) float64 4.382e+04 4.384e+04 ... 4.412e+04 4.415e+04\n",
       "Data variables:\n",
       "    P0         float64 1.004e+05\n",
       "    date       (time) int32 19900116 19900214 19900316 ... 19901115 19901216\n",
       "    datesec    (time) int32 0 0 0 0 0 0 0 0 0 0 0 0\n",
       "    OZONE_old  (time, lat, lev, lon) float64 ...\n",
       "    OZONE      (time, lev, lat, lon) float64 ...\n",
       "Attributes:\n",
       "    Conventions:                     NCAR-CSM\n",
       "    Source:                          AMIP II (symmetric for APE project)\n",
       "    Written_By:                      olson\n",
       "    Date_Written:                    August 22 2003\n",
       "    Host:                            zen\n",
       "    Command:                         ncgen\n",
       "    history:                         Wed Jul 30 08:35:58 2008: ncrename -v OZ...\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-7fb73d55-5358-4307-9ab8-4ecc926e6f0a' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-7fb73d55-5358-4307-9ab8-4ecc926e6f0a' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>lat</span>: 64</li><li><span class='xr-has-index'>lev</span>: 59</li><li><span class='xr-has-index'>lon</span>: 128</li><li><span class='xr-has-index'>time</span>: 12</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-8a0d17be-318b-4e08-8e48-be25dde0bf95' class='xr-section-summary-in' type='checkbox'  checked><label for='section-8a0d17be-318b-4e08-8e48-be25dde0bf95' class='xr-section-summary' >Coordinates: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lev</span></div><div class='xr-var-dims'>(lev)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.2842 0.3253 ... 959.0 1.004e+03</div><input id='attrs-9a8151d1-065b-43ea-9874-2f6369a4e524' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-9a8151d1-065b-43ea-9874-2f6369a4e524' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-11a67b32-b43c-4b21-94d5-b9159662ad98' class='xr-var-data-in' type='checkbox'><label for='data-11a67b32-b43c-4b21-94d5-b9159662ad98' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>hybrid level at layer midpoints (1000*(A+B))</dd><dt><span>units :</span></dt><dd>hybrid_sigma_pressure</dd><dt><span>positive :</span></dt><dd>down</dd></dl></div><div class='xr-var-data'><pre>array([2.84191e-01, 3.25350e-01, 3.71936e-01, 4.24615e-01, 4.84153e-01,\n",
       "       5.51432e-01, 6.27490e-01, 7.13542e-01, 8.11039e-01, 9.21699e-01,\n",
       "       1.04757e+00, 1.19108e+00, 1.35511e+00, 1.54305e+00, 1.75891e+00,\n",
       "       2.00739e+00, 2.29403e+00, 2.62529e+00, 3.00884e+00, 3.45369e+00,\n",
       "       3.97050e+00, 4.57190e+00, 5.27278e+00, 6.09073e+00, 7.04636e+00,\n",
       "       8.16371e+00, 9.47075e+00, 1.10000e+01, 1.27894e+01, 1.48832e+01,\n",
       "       1.73342e+01, 2.02039e+01, 2.35659e+01, 2.75066e+01, 3.21283e+01,\n",
       "       3.75509e+01, 4.39150e+01, 5.13840e+01, 6.01483e+01, 7.04289e+01,\n",
       "       8.24830e+01, 9.66096e+01, 1.13155e+02, 1.32514e+02, 1.55122e+02,\n",
       "       1.81445e+02, 2.11947e+02, 2.47079e+02, 2.87273e+02, 3.32956e+02,\n",
       "       3.84573e+02, 4.42610e+02, 5.07601e+02, 5.80132e+02, 6.60837e+02,\n",
       "       7.50393e+02, 8.49521e+02, 9.58981e+02, 1.00369e+03])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lon</span></div><div class='xr-var-dims'>(lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>0.0 2.812 5.625 ... 354.4 357.2</div><input id='attrs-c2ab461c-b5a0-465a-b3eb-b416bee70b56' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-c2ab461c-b5a0-465a-b3eb-b416bee70b56' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-f040f216-7ba7-406a-8b8d-300eb02614ff' class='xr-var-data-in' type='checkbox'><label for='data-f040f216-7ba7-406a-8b8d-300eb02614ff' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>longitude</dd><dt><span>units :</span></dt><dd>degrees_east</dd></dl></div><div class='xr-var-data'><pre>array([  0.    ,   2.8125,   5.625 ,   8.4375,  11.25  ,  14.0625,  16.875 ,\n",
       "        19.6875,  22.5   ,  25.3125,  28.125 ,  30.9375,  33.75  ,  36.5625,\n",
       "        39.375 ,  42.1875,  45.    ,  47.8125,  50.625 ,  53.4375,  56.25  ,\n",
       "        59.0625,  61.875 ,  64.6875,  67.5   ,  70.3125,  73.125 ,  75.9375,\n",
       "        78.75  ,  81.5625,  84.375 ,  87.1875,  90.    ,  92.8125,  95.625 ,\n",
       "        98.4375, 101.25  , 104.0625, 106.875 , 109.6875, 112.5   , 115.3125,\n",
       "       118.125 , 120.9375, 123.75  , 126.5625, 129.375 , 132.1875, 135.    ,\n",
       "       137.8125, 140.625 , 143.4375, 146.25  , 149.0625, 151.875 , 154.6875,\n",
       "       157.5   , 160.3125, 163.125 , 165.9375, 168.75  , 171.5625, 174.375 ,\n",
       "       177.1875, 180.    , 182.8125, 185.625 , 188.4375, 191.25  , 194.0625,\n",
       "       196.875 , 199.6875, 202.5   , 205.3125, 208.125 , 210.9375, 213.75  ,\n",
       "       216.5625, 219.375 , 222.1875, 225.    , 227.8125, 230.625 , 233.4375,\n",
       "       236.25  , 239.0625, 241.875 , 244.6875, 247.5   , 250.3125, 253.125 ,\n",
       "       255.9375, 258.75  , 261.5625, 264.375 , 267.1875, 270.    , 272.8125,\n",
       "       275.625 , 278.4375, 281.25  , 284.0625, 286.875 , 289.6875, 292.5   ,\n",
       "       295.3125, 298.125 , 300.9375, 303.75  , 306.5625, 309.375 , 312.1875,\n",
       "       315.    , 317.8125, 320.625 , 323.4375, 326.25  , 329.0625, 331.875 ,\n",
       "       334.6875, 337.5   , 340.3125, 343.125 , 345.9375, 348.75  , 351.5625,\n",
       "       354.375 , 357.1875])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>lat</span></div><div class='xr-var-dims'>(lat)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>-87.86 -85.1 -82.31 ... 85.1 87.86</div><input id='attrs-d43b4b64-f418-46c9-bd60-5b8cdb25b5ea' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-d43b4b64-f418-46c9-bd60-5b8cdb25b5ea' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-83919486-ec24-4b90-afda-57a524b9eca7' class='xr-var-data-in' type='checkbox'><label for='data-83919486-ec24-4b90-afda-57a524b9eca7' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>latitude</dd><dt><span>units :</span></dt><dd>degrees_north</dd></dl></div><div class='xr-var-data'><pre>array([-87.8638, -85.0965, -82.3129, -79.5256, -76.7369, -73.9475, -71.1577,\n",
       "       -68.3678, -65.5776, -62.7873, -59.997 , -57.2066, -54.4162, -51.6257,\n",
       "       -48.8352, -46.0447, -43.2542, -40.4636, -37.6731, -34.8825, -32.0919,\n",
       "       -29.3014, -26.5108, -23.7202, -20.9296, -18.139 , -15.3484, -12.5578,\n",
       "        -9.7671,  -6.9765,  -4.1859,  -1.3953,   1.3953,   4.1859,   6.9765,\n",
       "         9.7671,  12.5578,  15.3484,  18.139 ,  20.9296,  23.7202,  26.5108,\n",
       "        29.3014,  32.0919,  34.8825,  37.6731,  40.4636,  43.2542,  46.0447,\n",
       "        48.8352,  51.6257,  54.4162,  57.2066,  59.997 ,  62.7873,  65.5776,\n",
       "        68.3678,  71.1577,  73.9475,  76.7369,  79.5256,  82.3129,  85.0965,\n",
       "        87.8638])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>4.382e+04 4.384e+04 ... 4.415e+04</div><input id='attrs-36d940a5-93b3-46f1-9e18-2803cb943d10' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-36d940a5-93b3-46f1-9e18-2803cb943d10' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-57ecf187-721f-4258-a6c8-a1b724ef12c2' class='xr-var-data-in' type='checkbox'><label for='data-57ecf187-721f-4258-a6c8-a1b724ef12c2' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>time</dd><dt><span>units :</span></dt><dd>day number (Jan 1 = 1)</dd><dt><span>calendar :</span></dt><dd>365_days</dd></dl></div><div class='xr-var-data'><pre>array([43816., 43844., 43875., 43905., 43936., 43966., 43997., 44028., 44058.,\n",
       "       44089., 44119., 44150.])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-854d155b-99cb-4a5a-a73c-a3be10457251' class='xr-section-summary-in' type='checkbox'  checked><label for='section-854d155b-99cb-4a5a-a73c-a3be10457251' class='xr-section-summary' >Data variables: <span>(5)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>P0</span></div><div class='xr-var-dims'>()</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-99f5e254-ddf5-431b-8e96-60c300f2a3f8' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-99f5e254-ddf5-431b-8e96-60c300f2a3f8' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-bdfcf6af-dc6b-4d45-a488-fc8c0d0aef36' class='xr-var-data-in' type='checkbox'><label for='data-bdfcf6af-dc6b-4d45-a488-fc8c0d0aef36' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>reference pressure</dd><dt><span>units :</span></dt><dd>Pa</dd></dl></div><div class='xr-var-data'><pre>array(100369.)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>date</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>int32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-477f8172-13f3-49c8-9c13-7cb5010083dc' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-477f8172-13f3-49c8-9c13-7cb5010083dc' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-5518b65d-83d8-44a7-958d-f4f19bf8b4a1' class='xr-var-data-in' type='checkbox'><label for='data-5518b65d-83d8-44a7-958d-f4f19bf8b4a1' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>current date as 8 digit integer (YYYYMMDD)</dd></dl></div><div class='xr-var-data'><pre>array([19900116, 19900214, 19900316, 19900415, 19900516, 19900615, 19900716,\n",
       "       19900816, 19900915, 19901016, 19901115, 19901216], dtype=int32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>datesec</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>int32</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-b61c460d-c9dc-4e1c-ac18-49b851c08399' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-b61c460d-c9dc-4e1c-ac18-49b851c08399' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-442015d2-be75-4d88-aa56-c28f69ddf1ab' class='xr-var-data-in' type='checkbox'><label for='data-442015d2-be75-4d88-aa56-c28f69ddf1ab' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>seconds of current date</dd><dt><span>units :</span></dt><dd>s</dd></dl></div><div class='xr-var-data'><pre>array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int32)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>OZONE_old</span></div><div class='xr-var-dims'>(time, lat, lev, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-810db05b-af14-4513-9acb-8190bb4b31cd' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-810db05b-af14-4513-9acb-8190bb4b31cd' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-80b963ed-c97c-49c8-85d8-80f3bfc29ed7' class='xr-var-data-in' type='checkbox'><label for='data-80b963ed-c97c-49c8-85d8-80f3bfc29ed7' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>OZONE</dd><dt><span>units :</span></dt><dd>FRACTION</dd></dl></div><div class='xr-var-data'><pre>[5799936 values with dtype=float64]</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>OZONE</span></div><div class='xr-var-dims'>(time, lev, lat, lon)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>...</div><input id='attrs-9224701c-a332-4fe1-b7ba-b6d4ef027879' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-9224701c-a332-4fe1-b7ba-b6d4ef027879' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-c2bb42dc-83e9-4e9e-a9c6-195303b8b9e5' class='xr-var-data-in' type='checkbox'><label for='data-c2bb42dc-83e9-4e9e-a9c6-195303b8b9e5' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>long_name :</span></dt><dd>OZONE</dd><dt><span>units :</span></dt><dd>FRACTION</dd></dl></div><div class='xr-var-data'><pre>[5799936 values with dtype=float64]</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-a170dcff-b4d1-4757-9f62-7f5ab536f5bf' class='xr-section-summary-in' type='checkbox'  checked><label for='section-a170dcff-b4d1-4757-9f62-7f5ab536f5bf' class='xr-section-summary' >Attributes: <span>(8)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'><dt><span>Conventions :</span></dt><dd>NCAR-CSM</dd><dt><span>Source :</span></dt><dd>AMIP II (symmetric for APE project)</dd><dt><span>Written_By :</span></dt><dd>olson</dd><dt><span>Date_Written :</span></dt><dd>August 22 2003</dd><dt><span>Host :</span></dt><dd>zen</dd><dt><span>Command :</span></dt><dd>ncgen</dd><dt><span>history :</span></dt><dd>Wed Jul 30 08:35:58 2008: ncrename -v OZONE,OZONE_old apeozone_cam3_5_54.nc\n",
       " </dd><dt><span>DODS_EXTRA.Unlimited_Dimension :</span></dt><dd>time</dd></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (lat: 64, lev: 59, lon: 128, time: 12)\n",
       "Coordinates:\n",
       "  * lev        (lev) float64 0.2842 0.3253 0.3719 ... 849.5 959.0 1.004e+03\n",
       "  * lon        (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2\n",
       "  * lat        (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 82.31 85.1 87.86\n",
       "  * time       (time) float64 4.382e+04 4.384e+04 ... 4.412e+04 4.415e+04\n",
       "Data variables:\n",
       "    P0         float64 ...\n",
       "    date       (time) int32 ...\n",
       "    datesec    (time) int32 ...\n",
       "    OZONE_old  (time, lat, lev, lon) float64 ...\n",
       "    OZONE      (time, lev, lat, lon) float64 ...\n",
       "Attributes:\n",
       "    Conventions:                     NCAR-CSM\n",
       "    Source:                          AMIP II (symmetric for APE project)\n",
       "    Written_By:                      olson\n",
       "    Date_Written:                    August 22 2003\n",
       "    Host:                            zen\n",
       "    Command:                         ncgen\n",
       "    history:                         Wed Jul 30 08:35:58 2008: ncrename -v OZ...\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#  Put in some ozone\n",
    "import xarray as xr\n",
    "\n",
    "ozonepath = \"http://thredds.atmos.albany.edu:8080/thredds/dodsC/CLIMLAB/ozone/apeozone_cam3_5_54.nc\"\n",
    "ozone = xr.open_dataset(ozonepath)\n",
    "ozone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Dimensions of the ozone file\n",
    "lat = ozone.lat\n",
    "lon = ozone.lon\n",
    "lev = ozone.lev\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taking annual, zonal, and global averages of the ozone data\n",
    "O3_zon = ozone.OZONE.mean(dim=(\"time\",\"lon\"))\n",
    "\n",
    "weight_ozone = np.cos(np.deg2rad(ozone.lat)) / np.cos(np.deg2rad(ozone.lat)).mean(dim='lat')\n",
    "O3_global = (O3_zon * weight_ozone).mean(dim='lat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaF0lEQVR4nO3da4xc533f8e9/brvLXS6vS5pc0iFVrRyRimU7W1m2CqEoU0iJbFMIoIQpnBKBHAKBEl8QIJCStn5ToU4R5NIiMkLYTulGsMzKSkWnjmOBttIGCSQvKdniRYxoUSaXXIlLSryTuztz/n0xZ2bOXJaXGe7O8JnfB1icc57znHOeGZC/ffaZ55wxd0dERLpDqt0NEBGR+aPQFxHpIgp9EZEuotAXEekiCn0RkS6SaXcDrmX58uW+bt26djdDROSWsmfPnlPuPlRb3vGhv27dOsbGxtrdDBGRW4qZ/bRRuYZ3RES6iEJfRKSLKPRFRLqIQl9EpIso9EVEusi8h76ZPWhmh8zssJk9Pt/XFxHpZvMa+maWBv4c+EVgA/BrZrZhPtsgItLN5nue/j3AYXd/E8DMngE2AwfmuR23pDfeOc+3f3SiUmBWWU0UWbyV2F21v7hM7Gxw3DXPF2/UnTexv6osUb9SfvXrJReNrtfoPI3ad63rlffeyPVmqZusU/36G5eXj61rU6Nya/g+Vb/umtc563Ubv+ZZX0Pi302pTt37aZVjGr0vpfet9jWlrP5ape1U8no1/2alOfMd+sPAscT2OPDR2kpmtg3YBvD+979/flp2C/jJ5AX++w8OA6CvQZBuVffLoPSLpqq88ssnlSqWpUplVr1dWqZTVlWWKu8z0qnEely38pOKy1KkU5BJpUinjEzayKZSxWU6RTZtZNIpcukUPdl4mUnRk0nTm0uzsCdDf0+Ghb0ZVi/qY9GC7Jy8f/Md+o1+VdfFl7tvB7YDjI6OKt5iD961iiP/5aFZ97t7+ZeB15Qnyyp1vGq7tJ4s98Q5ao9nlvN4zTWrrlfTvqtdr1HbSuepfZ0N29ewXoP2Nii72vUavbbaa1XO7Q2vjc/epkavido6zbYnLk8el2xD9XvSuA3eoO3lsrh+7eu/1jWjuutXzhl54/NHVeeO67gXz9WojMQ+96pzRJETxfWixL7InUJUOodTiI/NF5yCOzMzEYWoQCHy8k8+ipgpOJem81yYynNlJqIZIysG+N4X7r/pf+HMd+iPA2sT22uAE7PUlRuUHFKp2TPfTRGZM+7FcJ0pONOFiJlCRL7gzBSi8vZ0vrQslpfKivsblCXqTcXHVs7jTOUjCpEznS+Uy0rHTucrx5TKClHrfdX77xiakyGt+Q79HwIjZrYeOA5sAf7dPLdBRBoohWkpBKdLIVoKtEQIVodoMvC8Ur98nJfXS+fJJ4M3ccxMlFiPAzoZ7KX1uRjeTKeMbDwUk0unikMyGSOXTpHLpMllUvSkUyzIZcimjVwmLk+niuvlshS5dJpsxujJpOMhnBQ92cR6Jk1PNrGeSdHfk2FxX5ZUam47afMa+u6eN7PfBv4OSANfc/f989kGkU7iXuxFTs1EXMkXmJqJmMoXimX5qKo3WVwv1kn2MEt1iuWFRN1KnXK9uvNVh/dchGkxBFPlQM3GIZncLoVof3l/9b7KejFcs+kUmcRxVecpB3Zxu6dqf6oc0tm0kc1UAj49x2HbKeb9KZvu/h3gO/N9XZHrMVOIuDxT4PJ0gUvTxeXlmeK47JWZYhhfmSlwJQ7n5LK0fypRr1w/DvTycqbAlTh4b4ZybzNT/UFhsmywL1v+8LBUVtqfTfRWk8Gci3u75d5vIiTLvd1EyObS1cdkUqZZNx2m4x+tLFJrKl/gwpU8F6cKXJjKJ0I6Xx3YVev5Wcrj4+L1mUJzXd3yTIxsmt74z/be+M/53myawb5sXVnxz/tS/cqyFMpVy3S6qiwZ5rl0as6HBCQcCn2ZF7VBfWEqz8WpPOfj5YUr+esqvzCVv6FgThksyGXozaZZkCv+lNYXL8gmyjP05dL0xdvJ9d5scb0U1L2ZdFV45zLdMzQgtz6FvlwXd+fCVJ6zl2fKP+fi5ZlLM1XltXVuJKj7c2n6ezIM9GYY6MnQn8uwdukCBnri7Xgec38uzUBvlv44oBfkMvRlS+vp8npPJqXhBZEEhX4XmilEvHtxmlMXpjh9obh89+J0VXifqQn2s5dnrjoNLZ0yFvdlWdSXZbAvy5IFOdYt62ewL8PC3mxVaFfW08UAj7cX5DLqMYvMMYV+AEq98NMXpjl9cYrJ88VlKdBLy1MXpjgdh3sjKYPBvmxVeK9d0sfiBcXt6p9ccRnv68+l1aMWuQUo9DtcFDmnLkwxfuYyx9+7zPF4eeLMZSYTgT41yyyQRX1Zlg3kWD7Qwwfet5Bl/T0sH+iJy3Lxeg9L+3Ms7MnoA0GRwCn022ymEPH22SvlMD9+5jLj712qhPvZK3XT+gZ7MwwvWcCKhT3cvmKAoTjEl/X3sHxhD8v6i2G+tD9HLqOvTBCRCoX+PLk0nefgxHn2nzjLvuNnOXLqIsffu8zb565QO1S+fKCHNUv62Di8iAc2vo/hJX0ML+4rLxf2zs2DmEQkfAr9OXD28gwHTpwrB/y+E+d4c/JCOdyX9ue4fcUA9962rC7QVy/uozebbu8LEJFgKfRbFEXOy2+9y96j77H/+Dn2nTjLT09fKu9/32Avdw0P8tDPrWLj6kHuGl7EqkW9+tBTRNpCod+kwycv8Nzecf73K8c5cfYKAO9fuoC7hgf5ldG1bFw9yMbVixha2NPmloqIVCj0b8B7F6f59o9P8K29x/nRsTOkU8b9I8t54pfu5P6RoTn70gMRkZtFoX8d9h59j7/4+5/w/ddPMlNw7lw1yH946E4+9aHVrFjY2+7miYhcN4X+Vbx7cZo//NvX+ebYMZb159j6sXX88kfWsGH1YLubJiLSFIV+A1HkfHPsGH/43de5cCXPtvtv47ObRhjo0dslIrc2pViN/SfO8gd/vY9Xj53hnvVL+c8P38UdKxe2u1kiIjeFQj/h/70xyWd2jLGwN8uf/OrdPPyhYU2tFJGgKPRjf//Pk/zm18e4bXk/T3/moywb0FRLEQmPQh948dBJtv3PPdw+NMDTn/koS/pz7W6SiMic6PrQf/HQSbZ9fQ8jK4uBv3iBAl9EwtXVj2B89+I0n//mq9y+QoEvIt2hq0P/v8ZTMv/kVz+kwBeRrtC1ob/np+/yzA+P8ei/Ws8H3qcpmSLSHboy9POFiD/4632sXtTLZzeNtLs5IiLzpitD/3/tGef1t8/znz65kX7dZSsiXaQrQ/+5vePcsXKABzaubHdTRETmVdeF/sTZy/zwrff45AdX625bEek6XRf6/+fHEwB84u7VbW6JiMj867rQ//aPJ7hreJD1y/vb3RQRkXnXVaF/7N1L/OjYGT75QfXyRaQ7NR36ZrbWzH5gZgfNbL+ZfS4uX2pmL5jZG/FySeKYJ8zssJkdMrMHbsYLuBEv/vMkAP92gz7AFZHu1EpPPw/8rrvfCdwLPGZmG4DHgd3uPgLsjreJ920BNgIPAk+ZWbqVxt+ofzx8itWLejW0IyJdq+nQd/cJd98br58HDgLDwGZgR1xtB/BwvL4ZeMbdp9z9CHAYuKfZ69+oKHL+6c3TfPz25Zq1IyJd66aM6ZvZOuDDwEvASnefgOIvBmBFXG0YOJY4bDwumxcHJs5x5tIM992+bL4uKSLScVoOfTMbAL4FfN7dz12taoMyn+Wc28xszMzGJicnW20iAP9w+BQA9/2L5TflfCIit6KWQt/MshQD/2l3fy4ufsfMVsX7VwEn4/JxYG3i8DXAiUbndfft7j7q7qNDQ0OtNLHsH39ympEVA6wY7L0p5xMRuRW1MnvHgK8CB939jxO7dgFb4/WtwPOJ8i1m1mNm64ER4OVmr38josh55eh7/Mv1S+fjciIiHauVp43dB/w68JqZvRqX/T7wJWCnmT0KHAUeAXD3/Wa2EzhAcebPY+5eaOH61+2t0xc5fyXP3WsWzcflREQ6VtOh7+7/QONxeoBNsxzzJPBks9ds1mvHzwLwc8OL5/vSIiIdpSvuyP3RsbP0ZFLcsXKg3U0REWmrrgj9146fYePqQTLprni5IiKzCj4F84WIfcfP8cE1i9vdFBGRtgs+9N88dZHLMwU+qA9xRUTCD/23Tl0E4PYVGs8XEQk+9E+cuQzAqkV9bW6JiEj7BR/6E2evkMukWNafa3dTRETaLvjQP37mMqsX9ZJK6cmaIiLBh/7E2Ssa2hERiQUf+ifOXGb1YoW+iAgEHvr5QsQ7564wvFhP1hQRgcBD/+T5KSKHVerpi4gAgYf+2cszACzuy7a5JSIinSHo0L80XXxyc19uXr9/XUSkYwUe+nkA+nta+doAEZFwBB76cU8/q56+iAgEH/rq6YuIJAUe+sWe/gKN6YuIAKGH/pRCX0QkKezQL/f0NbwjIgLBh36eXCZFWg9bExEBAg/9qXxETybolygickOCTsR8FJHTl6GLiJQFnYgzeSeT1tCOiEhJ2KFfiMiqpy8iUhZ0Is5EruEdEZGEoBNxJh9peEdEJCHo0M9HGt4REUkKOhGnC05GoS8iUhZ0Is7kI3Ia3hERKQs69PNRRCYV9EsUEbkhLSeimaXN7BUz+5t4e6mZvWBmb8TLJYm6T5jZYTM7ZGYPtHrta5kpOFndkSsiUnYzEvFzwMHE9uPAbncfAXbH25jZBmALsBF4EHjKzOb08ZeFyMnouTsiImUthb6ZrQEeAr6SKN4M7IjXdwAPJ8qfcfcpdz8CHAbuaeX615KPXA9bExFJaLWn/6fA7wFRomylu08AxMsVcfkwcCxRbzwuq2Nm28xszMzGJicnm25cIYrU0xcRSWg69M3sE8BJd99zvYc0KPNGFd19u7uPuvvo0NBQs01UT19EpEYr3y5yH/ApM/sloBcYNLO/At4xs1XuPmFmq4CTcf1xYG3i+DXAiRauf00a0xcRqdZ0T9/dn3D3Ne6+juIHtN93908Du4CtcbWtwPPx+i5gi5n1mNl6YAR4uemWX4d8wUlryqaISNlcfI/gl4CdZvYocBR4BMDd95vZTuAAkAcec/fCHFy/TD19EZFqNyX03f1F4MV4/TSwaZZ6TwJP3oxrXo985KR1R66ISFnQYx+avSMiUi3o0NfsHRGRakGHfiFy0qbQFxEpCT/01dMXESkLOvQjd1IKfRGRsqBDX8M7IiLVgg19dydy1NMXEUkINvSj+Kk+6umLiFQEG/qFOPX1FbkiIhXBRmLkxdDX8I6ISEX4oa/hHRGRsmBDvzy8o9AXESkLNvSj+Lu8NLwjIlIRbOgXvNTTb3NDREQ6SLihX569o9QXESkJNvQ97umbxvRFRMqCDf3yzVnq6YuIlAUc+nFPv83tEBHpJMGGftzR1zx9EZGEYEM/ikpj+m1uiIhIBwk29OPRHX2QKyKSEG7oU3oMQ5sbIiLSQYIN/dLsHY3pi4hUBBz6GtMXEakVbOhrTF9EpF7Aoa8xfRGRWsGGfmlM33R7lohIWbChr9k7IiL1gg390vP0NaYvIlIRbuhr9o6ISJ2WQt/MFpvZs2b2upkdNLOPmdlSM3vBzN6Il0sS9Z8ws8NmdsjMHmi9+demefoiIhWt9vT/DPiuu/8scDdwEHgc2O3uI8DueBsz2wBsATYCDwJPmVm6xevPSk/ZFBGp13Tom9kgcD/wVQB3n3b3M8BmYEdcbQfwcLy+GXjG3afc/QhwGLin2etfS2mefirYASwRkRvXSiTeBkwCf2lmr5jZV8ysH1jp7hMA8XJFXH8YOJY4fjwuq2Nm28xszMzGJicnm2pcpaevvr6ISEkroZ8BPgJ82d0/DFwkHsqZRaP09QZluPt2dx9199GhoaGmGlc+sTJfRKSsldAfB8bd/aV4+1mKvwTeMbNVAPHyZKL+2sTxa4ATLVz/qip35Cr1RURKmg59d38bOGZmH4iLNgEHgF3A1rhsK/B8vL4L2GJmPWa2HhgBXm72+tduX3GpyBcRqci0ePzvAE+bWQ54E/gNir9IdprZo8BR4BEAd99vZjsp/mLIA4+5e6HF689KX5coIlKvpdB391eB0Qa7Ns1S/0ngyVaueb30dYkiIvWCndBY6ukr9EVEKoINfU3ZFBGpF2zoU/66xPY2Q0SkkwQb+pG+OUtEpE6woa/n6YuI1As29Cs9/fa2Q0SkkwQb+qU7cnV7lohIRcChX1ymNb4jIlIWbOhHrjF9EZFaAYd+canHMIiIVAQc+noMg4hIrXBDP9KjlUVEaoUb+vogV0SkTsChrw9yRURqBR/6egyDiEhFsKHvmr0jIlIn2NDX8I6ISL1gQ7+g2TsiInWCDf3y8I66+iIiZcGGvoZ3RETqBRz6xaWGd0REKgIOfT2GQUSkVrCh764PckVEagUb+hreERGpF2zoV6ZstrkhIiIdJNjQL8/eUeqLiJSFH/oa3hERKQs49IvLtEJfRKQs2NAvjekr80VEKoINfU3ZFBGpF2zo65uzRETqtRT6ZvYFM9tvZvvM7Btm1mtmS83sBTN7I14uSdR/wswOm9khM3ug9ebPTlM2RUTqNR36ZjYMfBYYdfe7gDSwBXgc2O3uI8DueBsz2xDv3wg8CDxlZunWmj87d8dM35wlIpLU6vBOBugzswywADgBbAZ2xPt3AA/H65uBZ9x9yt2PAIeBe1q8/qwK7hrPFxGp0XTou/tx4I+Ao8AEcNbdvwesdPeJuM4EsCI+ZBg4ljjFeFxWx8y2mdmYmY1NTk421b7INV1TRKRWK8M7Syj23tcDq4F+M/v01Q5pUOaNKrr7dncfdffRoaGhptoXxcM7IiJS0crwzi8AR9x90t1ngOeAjwPvmNkqgHh5Mq4/DqxNHL+G4nDQnIgiDe+IiNRqJfSPAvea2QIrflq6CTgI7AK2xnW2As/H67uALWbWY2brgRHg5Rauf1WRa7qmiEitTLMHuvtLZvYssBfIA68A24EBYKeZPUrxF8Mjcf39ZrYTOBDXf8zdCy22f1aFSMM7IiK1mg59AHf/IvDFmuIpir3+RvWfBJ5s5ZrXy93V0xcRqRH0Hbka0xcRqRZs6Bfn6be7FSIinSXY0HfdnCUiUifY0I8iPVZZRKRWuKHvrjtyRURqBBz6etiaiEitYEPf3UkF++pERJoTbCxG+iBXRKROwKGvefoiIrUCDn09hkFEpFawoe/q6YuI1Ak29CPdkSsiUifo0LeG39siItK9gg19d92RKyJSK9jQ1+wdEZF6wYa+bs4SEakXbCzq5iwRkXoBh76evSMiUivg0NeUTRGRWsGGvm7OEhGpF2zoq6cvIlIv6NDXmL6ISLWAQx/djysiUiPY0NcXo4uI1As49NHNWSIiNYKNRT1wTUSkXsChDylN3xERqRJs6LumbIqI1Ak29DV7R0SkXrCh72j2johIrWuGvpl9zcxOmtm+RNlSM3vBzN6Il0sS+54ws8NmdsjMHkiU/7yZvRbv+282x3dORZEeuCYiUut6evr/A3iwpuxxYLe7jwC7423MbAOwBdgYH/OUmaXjY74MbANG4p/ac95UegyDiEi9a4a+u/9f4N2a4s3Ajnh9B/BwovwZd59y9yPAYeAeM1sFDLr7P7m7A19PHDMn9MA1EZF6zY7pr3T3CYB4uSIuHwaOJeqNx2XD8XpteUNmts3MxsxsbHJysqkGFp+909ShIiLButkf5DaKWb9KeUPuvt3dR919dGhoqKmG3H/HED//M0uuXVFEpItkmjzuHTNb5e4T8dDNybh8HFibqLcGOBGXr2lQPmf+4yc2zOXpRURuSc329HcBW+P1rcDzifItZtZjZuspfmD7cjwEdN7M7o1n7fz7xDEiIjJPrtnTN7NvAP8aWG5m48AXgS8BO83sUeAo8AiAu+83s53AASAPPObuhfhUv0VxJlAf8Lfxj4iIzCMrTqbpXKOjoz42NtbuZoiI3FLMbI+7j9aWB3tHroiI1FPoi4h0EYW+iEgXUeiLiHQRhb6ISBfp+Nk7ZjYJ/LSJQ5cDp25yc0Ki92d2em9mp/fm6jrp/fkZd697pEHHh36zzGys0XQlKdL7Mzu9N7PTe3N1t8L7o+EdEZEuotAXEekiIYf+9nY3oMPp/Zmd3pvZ6b25uo5/f4Id0xcRkXoh9/RFRKSGQl9EpIsEGfpm9qCZHTKzw2b2eLvb0ynMbK2Z/cDMDprZfjP7XLvb1GnMLG1mr5jZ37S7LZ3GzBab2bNm9nr8b+hj7W5TpzCzL8T/p/aZ2TfMrLfdbZpNcKFvZmngz4FfBDYAv2Zm+hqtojzwu+5+J3Av8JjemzqfAw62uxEd6s+A77r7zwJ3o/cJADMbBj4LjLr7XUAa2NLeVs0uuNAH7gEOu/ub7j4NPANsbnObOoK7T7j73nj9PMX/tLN+QX23MbM1wEPAV9rdlk5jZoPA/cBXAdx92t3PtLVRnSUD9JlZBljAHH8dbCtCDP1h4FhiexwFWx0zWwd8GHipzU3pJH8K/B4Qtbkdneg2YBL4y3j46ytm1t/uRnUCdz8O/BHFbxGcAM66+/fa26rZhRj61qBM81ITzGwA+BbweXc/1+72dAIz+wRw0t33tLstHSoDfAT4srt/GLgI6PMywMyWUBxNWA+sBvrN7NPtbdXsQgz9cWBtYnsNHfyn1nwzsyzFwH/a3Z9rd3s6yH3Ap8zsLYpDgv/GzP6qvU3qKOPAuLuX/jJ8luIvAYFfAI64+6S7zwDPAR9vc5tmFWLo/xAYMbP1Zpaj+IHKrja3qSOYmVEckz3o7n/c7vZ0End/wt3XuPs6iv9mvu/uHdtbm2/u/jZwzMw+EBdtAg60sUmd5Chwr5ktiP+PbaKDP+TOtLsBN5u7583st4G/o/gp+tfcfX+bm9Up7gN+HXjNzF6Ny37f3b/TvibJLeR3gKfjztSbwG+0uT0dwd1fMrNngb0UZ8i9Qgc/jkGPYRAR6SIhDu+IiMgsFPoiIl1EoS8i0kUU+iIiXUShLyLSRRT6IiJdRKEvItJF/j9EBXceyQhRcgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot( O3_global*1E6, lev)\n",
    "ax.invert_yaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to create another instance of the model, this time using the same vertical coordinates as the ozone data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (1,) \n",
      "  Tatm: (59,) \n",
      "The subprocess tree: \n",
      "Untitled: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#  Create the column with appropriate vertical coordinate, surface albedo and convective adjustment\n",
    "col2 = climlab.BandRCModel(lev=lev)\n",
    "print(col2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Set the ozone mixing ratio\n",
    "col2.absorber_vmr['O3'] = O3_global.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 730 steps, 730.4844 days, or 2.0 years.\n",
      "Total elapsed time is 1.9986737567564754 years.\n"
     ]
    }
   ],
   "source": [
    "#  Run the model out to equilibrium!\n",
    "col2.integrate_years(2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x169c510d0>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEfCAYAAABbIFHdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABFpElEQVR4nO3dd3hUZdrH8e8dSAIkEHrooUOokQAqi0AUEMSCZe3dFbG+9rUhoKsi9ra6LquC2LFgV8BQLChFUJCmEKrSW5Ce+/3jOYFhnJRJZjKT5P5c17kyc86Zc34zMOee055HVBVjjDHGX0ykAxhjjIlOViCMMcYEZAXCGGNMQFYgjDHGBGQFwhhjTEBWIIwxxgRkBcIYU2wi0kxEPhCRjSKiIvKKN/7QY595s0RkagRimiBVjHQAUzJEJJgbXpqpala4spQ2IpIGDAZesc8lT68AnYAHgD+A3yKaxoSEFYjy4yK/58cBQ4AXgRl+0zaWSKLSIw0YDkwFsiIZJBqJSDzu/9Ozqvqo3+TKwMGST2VCwQpEOaGq432fi0hFXIH4zn9aWSYiVVV1Z6Rz+IqmTCISC1RQ1T1BvCwZEGCL/4Qgl2OijJ2DMEcQ52oRmSMif4rIThHJFJEMv/maeseXR4jI2SIyT0R2i8ivInKZN08TEZkgIlu85YwXkap+y3nFW04dERknIptFZJeITBGRo/LIeI6IfO0t808R+V5Ezgown3rLP8GbPxv4yJvWQEQe83JvFZE9IvKLiPxTRCr4LGME8LL3NNNbpu8x9hHe86YB1v+XY+35ZfKmdxWR90Vkk4jsFZElInK3V9AL5LP8viIy0/t8/hCRp0QkwW/e3OztReRxEVkD7AGO8abXFpHnRGS1iOzz/j4nIrV8lvEKsNJ7Otzn8+njm6eQ2Qv13r2874jIWm++P7z/o4MKsx5TeLYHYfy9CpwHTMBtGOOBC4BJInKGqn7oN//JwFDg37hfkFcAL4nIPuBB4CvgLqAbcDluA/SPAOv93Hv9CKAecB0wXUSOVdUFuTOJyL+Au735hwE5wOnAOyJynao+57fcrsCZwH+BsT7jOwFnAO/jjpfHAgOBUUBz4CpvvveA+ri9rQeBRd744hxjD5hJRE7y8vwKPIb7PI4F7sMd5vp7IZffBTjLW/44IAO4AeggIv1UNcdv/teA3d46FfhdRJKAb4GWwEvAXOAo4GrgeBHp7u31/AeYBzzhZX/PW+YiglDY9+4Vp6+8l72AK061cZ/p0cAnwazXFEBVbSiHA3ApbmNwqc+4071xQ/zmrQjMBlYA4o1r6s27C0jxmbcOrgjkADf7Lec9YB+Q6DPuFW857+Uu2xuf7i3jc59xXbx5Hwzwfj4AdgBVfcapN/QNMH9l3/X5jH8Vd8y8foDPqk+A+Ud405oGmJYFTPUbFzATUAl3cnc6UNFv2k15rT/AOnOXP9hv/FPe+HMDZJ8aYJ0PeNOu8Rt/rTf+fp9xuf8XRuSR55X8Ppdg3jtwqvf87Eh/h8rDYIeYjK8LgZ3AB97hhdoiUhuojjsM0hRo5feaD1Q19xADqroRWILbuPv/mp+B+6XeNMC6R6u3BfCWMweYBPQVkURv9AW4jcNY33xexg+Bqrhfnb7mq+pk/5Wp6u7c9YlInIjU9JbzBe7Qa9cAGUMlUKZ+uGP5LwPV/d7bp948/Qu5/CWq+oHfuFHe39MDzP+kqh7wG3c67mKFF/3G/wfYlMdyiiqY977d+ztQRKqFMIMJwA4xGV+puI3s+nzmSQaW+jxfHmCercDvqro3wHiAWvxVoEMSv+A2DCnAQi+fAIsLyOdraaCZvOPadwAX4w6jiN8sNfJZR3EFypTq/X0pn9f5v7e8/OWzVNXfRWQb7vBZYfI0A2b7Fw5VPSAiS3B7c6FS6PeuqtNEZBxur+4CEZkFTAbeUtVfQpjJYAXCHElwvxrPz2eeBX7P87qEMb9LG/03xoWdT3B7EAPzWf5Cv+d/5jHf48D1wFu4wykbgP24Dd/DFP4CjvzuL8nr+xUoU+57vQ13TD+QdcXMlNfnntdnVFKCeu+qeomIPAKcBPQEbgHuFpEbVfXZcAYtb6xAGF/LgNbATFXNLuF1pwIzA4w7yOGrZJYBA4BVqhrUSdAALgKmq+q5viNFpGWAefMrArmXdtbE5x4JEamEO7n9ayHzLPP+7gp0SCxI7fxHiEh9IInAe3yBLAfaiEhF370Ib8+rdRDLKYyg37u6CxcWAKNFpDrwPTBKRJ7zPVRpisfOQRhf43D/Jx4KNFFECnuIoyhuF5FDv3BFpAvQF5jiU6xe9f4+6Hspqs9r6gaxvoP4/aL2LgO9KcC8ueuvGWBa7uGZvn7jbyK479cXuL2YO0TkL+sRkcrid4lwPtqIyGC/cf/0/n5QyGV8gLvgwP+Ksyu98e8XcjmFUej37p0rOuJzVdVtuAsoquBOeJsQsT0Ic4iqThCRl4HrvA30x7gTko1wJ39bEvgYdiikAF+IyIe4X97X4S69vM0n3ywRGQ6MBOaJyDu4Qw/1cVc9nQTEFXJ9E4CrROQt3DHsZNxluJsDzDsLd9L9bhGpgbtya4Wqfu+9djFwn3cJ5grcYY9jcJ9doajqLhG5GLdhXiIiL+H2PqoDbXGX5J6Ou+KoID8D40Xkv7hf5xm4y16n4Q6pFcZo3KWlz3n/F37EXeZ6Be4ihNGFXE6BgnzvFwM3iUjuJbH7gd7AicDbqro7VLkMdplreR0IcJmrz7SLcFcc7cBdspqFuwz1HJ95mpL3pY1Tgax81tnHZ9wr3rg6uD2Ezbhj4l8B6XlkH4T71bkF2AusBj4Drvab7y+XWPpMqwI8gjt8tQe3Ib0DOCHQ5wJcgjtpvs9/ubhDLp97ubcBbwMNyfsy14CZvOkdgPHAWm9d63H3IwwDahbi31W9z7Qv7rDLbm8Zz+BzCbA37wjyuETXm14Hd3/LGtyGeA3uyrTafvPl93+hwMtcg3nvuHsixuKKwy7c/9H5uPMQ8ZH+XpW1IfeadmMiwrvL9hJVLeyJa5MPcY0yjlXVSyOdxZR+dg7CGGNMQFYgjDHGBGQFwhhjTEB2DsIYY0xApf4y19q1a2udOnVISEgoeOYI2bVrV1TnA8sYKpYxNCxjaOSXcc6cOZtUtU6+C4j0ZVTFHdLT0zUzM1OjWbTnU7WMoWIZQ8MyhkZ+GXFtbeW7fbVzEMYYYwKyAmGMMSYgKxDGGGMCKvUnqQPZv38/a9asYc+e6OgvPSkpiUWLitv4aOhUqlSJRo0aERsbG+koxpgoViYLxJo1a6hatSpNmzbFp4HQiNm5cydVqxa2Ic7wUlU2b97MmjVraNasWaTjGGOiWJk8xLRnzx5q1aoVFcUh2ogItWrVipq9K2NM9CqTBQKw4pAP+2yMMYVRZguEMaYEbdwIt94Kv/0W6SQmhKxAhEmFChVIS0ujQ4cOnH322Wzbtu3QtEcffZS2bdvSoUMHOnfuzLhx4wDo06cPbdq0IS0tjbS0NM4666wIpTcmSF9/DU89Ba1awaBB8PnnkJMT6VSmmKKmQIhImojMFJF5IjJbRLpHOlNxVK5cmXnz5rFgwQJq1KjBc889B8ALL7zApEmT+OGHH1iwYAHTp09HfdrDeu2115g3bx7z5s1jwoQJkYpvTHBOPx1WroR774W5c2HgQGjTBp58Enx+HJnSJWoKBK4Lw5GqmgbcSwi7NIy07t27s3btWgAefPBB/v3vf1OtWjXAXQJ7ySWXRDKeMaHRoAGMGOEKxeuvQ926cNNN0KgRXH01/PRTpBOaIEXTZa4KVPMeJ+H6Gi62G5ctY152dsEzBiEtMZEnW7Uq1LwHDx5k2rRpXHXVVezcuZOdO3fSokWLPOe/4IILqFy5MgD9+vXjkUceCUlmY0pMXBycd54b5s6F556DV16BF16Azp3hoovg/POhfv1IJzUFiJrmvkUkFdfPsOD2bHqo6so85h0CDAFITk5OHzNmDImJiYemJyUl0bJlSwD+uXo1P//5Z0izdqxShYcbN853nurVq9O+fXtWrVpF586dmThxIrt27aJDhw6sWrUq4GtOOukk/vWvf9GlS5eQ5g3k119/Zfv27YeeZ2dnH/EZRiPLGBqRyFhx+3aSp0whedIkqi1ejMbEsLVLF9b368em447joPejKJIZg1XaM2ZkZMxR1a75LqCg1vxCOQCTgQUBhtOAp4EzvfnOBiYXZpmBWnP95ZdfgmjvMDwSEhJUVXXbtm167LHH6lNPPaWqqo0aNdLffvst4Gt69+6ts2bNKpF8/p9RaW+ZMlpYxkJYvFj1nntUmzZVBdUqVVQvuED1s89U9++PjoyFUNozEm2tuapqX1XtEGCYCFwCvOfN+g5Qqk9S50pKSmL06NE8+uij7N+/nzvvvJNrr72WHTt2ALBjxw5efPHFCKc0pgS1aQP33w/Ll8OMGe6Q06efuhPbjRvD3XdT6fffI53SEF3nINYBvYGpwPHAsoimCaHOnTvTuXNn3nzzTa6++mqys7Pp1q0bsbGxxMbGcssttxya1/ccRO3atZk8eXKkYhsTXiLQs6cbnnoKPvkEXn4ZRo3iaFV33mLIEDj1VLB2wyIimgrElcBTIlIR2IN3jqG0yvY7Mf7RRx8denz77bdz++23/+U1U6dODXcsY6JTfDyccYYbVq8m6957aTZ5Mpx1FiQnw2WXwT/+Aflc4GFCL2ouc1XVr1U1XVU7q+rRqjon0pmMMRHQuDErL7kEsrLg44/h6KNh9Gho2RL69oUJE2D//kinLBeipkAYY8wRKlRwd2VPnAirVsF998GyZfD3vx86V0FWVqRTlmlWIIwx0a9hQxg2zJ3Y/vhj6N4dRo2C5s3hpJNcETlwINIpyxwrEMaY0iN3r+LDD2HFClc05s+HwYOhaVMYPhxWr450yjLDCoQxpnRq0gRGjnRNe7z/PnTo4C6fbdoUTjkFPvrI9iqKyQqEMaZ0q1jR7UF8/rlrbvyOO2D2bHd5bLNmrn0o26soEisQYRLO5r4/+OADOnXqRNu2benYsSMffPDBoWnDhg2jU6dOpKWl0b9/f9atC0mTVsaUDs2awQMPuJPa774L7du7k9u2V1EkViDCJFzNfc+fP59bb72ViRMnsnjxYj788ENuvfVWfvJayrztttv46aefmDdvHieffDL33XdfybxhY6JJbKy7pyLQXkXTpq5Z8pUBm3ozPqxAlIBQNvf96KOPctddd9GsWTMAmjVrxp133nmo1dfc5QLs2rXLuhc1xnev4r33oFMn+Ne/3PiBA935C7uvIqBoupM6PG68EebNC+0y09JcRyiFEOrmvhcuXMitt956xLiuXbse2kMBuPvuuxk3bhxJSUlkZmYW7j0ZU9bFxrqOjXI7N3rpJfjf/9yeRr16cOmldre2H9uDCJPdu3eTlpZGrVq12Lp1K/369UNVC/xF73uIKVBfEIGW4T/ugQceYPXq1VxwwQU8++yzoXlDxpQlKSnuCqisLHfJbLduh+/W7tMHXn0VQtxNQGlU9vcgCvlLP9Ryz0Fs376dgQMH8txzz3HDDTeQkJDA8uXLad68eZGW2759e2bPnk2nTp0OjZs7dy7t2rX7y7znn38+gwYNYuTIkUV+H8aUaRUrupPXp5wCa9fC2LFuz+Lii+G66+Dcc+GKK1wBKYeHa20PIsxC3dz3rbfeykMPPUSW18RAVlYWDz744KEWYZctO9wI7ocffkjbtm1D92aMKcsaNoS77nLNeUyd6i6dffVV1xZUx47w+OOwaVOkU5aosr8HEQVC2dx3WloaDz/8MKeccgr79+8nNjaW0aNHk5aWBsAdd9zBkiVLiImJISUlhRdeeKHE3qcxZYII9O7thqefhrfecnsVt9ziCsj558P110c6ZYmwAhEm4Wzu+4wzzuCMM84IOO3dd98tfEhjTP6SklyfFEOGwIIFrn/tcePg5ZdJ69gR7rnHnfQuo/1V2CEmY4wpjA4d4Pnn3bmKxx8nfvNmOOecw5fRbtwY6YQhFzUFQkTeEpF53pAlIvMinckYY/6ienW46Sa+HzfO3Zndvr3bk2jUCM47z52/8Ln5tTSLmgKhqueoapqqpgHvcrh/6qIuLyS5yiL7bIwJgQoV4OST4YsvYNEiuPpqd+d2RgakprqT2ps3RzplsURNgcgl7oL+s4E3irqMSpUqsXnzZtsQBqCqbN68mUqVKkU6ijFlR9u27pL6devcpbK1armT2g0awAUXwPTppXKvQqJtIyoivYDHVbVrPvMMweuzOjk5OX3MmDEkJib6TichIYEKFSqEPW9hFOYGuZJ08OBBdu3adUQBzc7OPuIzjEaWMTQsY2gUlDFh+XLqf/wx9b78koq7dpHdvDmrzj2Xjccfj5bQtim/jBkZGXPy284CbuNVUgMwGVgQYDjNZ57ngVsKu8z09HTNzMzUaBbt+VQtY6hYxtAoUxmzs1XHjFFNTVUF1ZQU1aefduPDLL+MwGwtYPtaooeYVLWvqnYIMEwEEJGKwBnAWyWZyxhjwiYhwd2NvWCB6xq1YUO44YbDzX1E8c130XYOoi+wWFXXRDqIMcaEVEyMa278m2/g66+hRw/XmVFKCtx5J2zdGumEfxFtBeJcinFy2hhjSoW//c01ErhwoWvS4+GHoXlzGDUqqhoJjKoCoaqXqqq1DWGMKR/atYPXXnNdEvTs6fYkWrZ0N+RFQR8VUVUgjDGmXOrUyd10N2OG64/immvcvRTTpkU0lhUIY4yJFj17unsmPvnEPc/IgNtug717IxLHCoQxxkQTETjpJHfYacgQePRR1x/F/PklHsUKhDHGRKPERHjhBfj4Y9iwwRWJ//ynRCNYgTDGmGg2aJC7h6JvXxg6tESLhBUIY4yJdrVrwwcfuGIxdCi8/HKJrNYKhDHGlAZxcTBhAvTr5+7MnjAh7Ku0AmGMMaVFpUpuT6JdO3jssbCvzgqEMcaUJlWqwPHHw88/Q05OWFdlBcIYY0qbzp1h1y5Yvjysq7ECYYwxpUlOjuvWFGDPnrCuygqEMcaUFqpw440wfjzcfz906BDW1VmBMMaY0kDV9R/xzDNw001w991hX2XFsK/BGGNM0e3a5Vp8ffZZd2L60ktd8xsl0I2x7UEYY0w0+u03uPlm1wPdVVdBhQrwv//BmDGu86ESEFV7ECJyPXAdcAD4RFVvj3AkY4wpOTk58OWX7jDSZ5+5onDWWXDdda4HuhLYa/AVNQVCRDKA04BOqrpXROpGOpMxxoSVqttTmDrVDZmZsG4d1KsH997rWnNt0CBi8YpUIESkHtAAqAxsAlao6r5iZrkaGKWqewFUdUMxl2eMMdHFvyBMnQpr17ppycnQpw+cdhqceaZrWiPCCl0gRKQr8A9gANDYb/I+EZmF60/6NVXdUYQsrYHjROQBYA9wq6rOKsJyjDEmOhSmIOQObdqU+CGkgoiq5j+DKwyPAr2An4HpwI/ARmA3UBNoBhwNnOC9bDTwmKru8VvWZKBegNXcDTwAfAX8H9ANeAtorgECisgQYAhAcnJy+pgxY0hMTCzE242M7OzsqM4HljFULGNolMqMOTlU+uMPEpYvJ3HFChKWLydpwQLiN20CYF+NGmxLSzs0/Nm4cdgLQn6fY0ZGxhxV7ZrvAlQ13wHYBTwJpBZi3krAOcAcYFhB8/u99nOgj8/z34A6Bb0uPT1dMzMzNZpFez5VyxgqljE0oj7jpk364xNPqD71lOo//qF69NGqiYmqbp/BDc2aqZ59turzz6suWqSak1PiMfP7HIHZWsD2tTCHmFqo6h+FmA91ewxvAW+JSHJhXuPjA+B4YKqItAbicOc3jDEmMnbvhkWL3P0HvsPvv5OWO0+tWtCxI1x2mfvbsSO0bw9Vq0YweGgUWCAKWxwCvG59kC95CXhJRBYA+4BLvCpnjDHhtX27O1fw66+wePHhQrBs2eEWUytVcs1s9+8PHTsyPyeHzhde6K44irJzB6FS1KuYOuHOSdQC/qOqf4hIS2C9qu4syjLVXQV1YVFea4wx+VKFTZtcAfj118PFIPfvJp+DFSLQooXbEzjnnMN7BS1buvsSPFunToX69Uv+vZSgoAqEiMQD44EzAAEU+Aj4A3dieilwR4gzGmNMwXJy3D0EeRWBnT6/XUWgSRO30T/jDPe3RYvDfxMSIvc+okiwexAPAH2Bi4BJgO9hpM+Aa7ACYYwJl/37YdWqIzf8uY+XLz+y+evYWGjWzG3we/Y8sgg0bQrx8RF7G6VFsAXiPOAeVX1dRCr4TVsBNA1JKmNM+fTnn7ByJTV/+AGWLIGVK48c1q07she1ypXdRr91azjppCP3Aho3hopR01hEqRTsp1cLWJTHtBjASrIxJjBV2Lbt8MY+K+uvBcA7F9Ap9zUVK0KjRpCS4rrZTElxv/5btXJFoH79MnuCOBoEWyBWAMfibmjz1x1YUuxExpjSKScH1q//60bfd9jpdw1L5cpuo5+SAunphx7/uGULRw0e7NohquB/sMKUlGALxDjgLhHJAt7zxqnX0N5NwIjQRTPGRJU9e2D1ancOwHfI3fivXg179x75murV3S/+Fi0O7wH4DrVrB9wD2D51qjtEZCIq2AIxGugMvAqM8cZ9jbuD+k1VfSaE2YwxJSUnBzZsOHLD718MNvi1nyniDvE0bgxdusDpp/+1AFSrFpn3Y0IiqAKhqgeBc0XkOeBEoC6wGfhcVaeFIZ8xJhR27Tq0wa8/aRJ89dVfi8E+vwaZExLcRr5JE1cAmjQ5cmjYMCpaHDXhE0xrrnHATOAOVf0SmBG2VMaYwlN1v+6zstwQ6Nf/5s2HZm8Drkeyhg3dhr5bN9e8tH8BqF7dTgCXc4UuEKq6T0Sa4Xp7M8aUlNy7gHMLwIoVRz5eudK1GeQrKenwhv6YY47Y8H+3di3HnnmmXQJqChTs/5BJQH8CX8VkjCkKVdiy5fBG378IZGW5Q0S+atZ0J3/bt4dBg9zjpk0PHxJKSspzdXunTrXiYAol2P8lzwDjRaQirvXV33HNbRyiqstDE82YMmTbtsC//nMf+1/+mXv1T6tWrnG43AKQO9jJX1MCgi0QuSeib8Zd1hqIXbRsyqddu0hctgx+/93dBZw7LF/uWgv1VbWqawaiWTPIyHAb/WbNDheA6tVLPr8xfoItEJeFJYUxpUVODqxZ4zb8ixcfLgKLF8OaNRzqnkvEbejbtIG//e3whj+3CNSoYSeATdQL9jLXseEKYkxUyc6GpUv/WgiWLDnyhHC1aq4IZGRAmzYsPHCA9mee6doDqlQpcvmNCQE7U2XKt+xsmDMHfvrpyEKwZs3heUTcL/82bVzn8m3busdt27qO5332BDZOnQodOpT42zAmHILtD+KlAmZRVb2iKEFEZARwJbDRG3WXqn5alGUZE9DBg7BwIfzwA3z/vRsWLjzcOmi1am6j7+0NHCoEtjdgyqlg9yCOx++qJaAmUBXY5g3F8YSqPlrMZRjjrFlzZDGYPfvw5aI1akD37q55iKOPhqOOKtNdRxpTFMGeg2gaaLyI9AJeAC4IQSZjgrdzpysAvgVh3To3LTYW0tJcp/JHH+2Gli2tGBhTAFH13yEo4oJErgQuUdWeRXz9COBSYAcwG7hFVbfmMe8QYAhAcnJy+pgxY0hMTCzKaktEdnZ2VOeD0pmx8po11J0yhTrTp5OQlYV4h4r+bNiQnW3bsiM1lR2pqWS3bImWUJtBpfFzjEaWMTTyy5iRkTFHVbsGnJhLVUMyACcAuwqYZzKwIMBwGpCMu4ciBte16UuFWW96erpmZmZqNIv2fKqlKOOaNaqPPaaanq4KqiKqvXqpDh+u+umnqps2RT5jlLOMoVHaMwKztYDta0iuYvLurL4UWJPffKrat5DL+y/wcfGTmTJh82Z4913S/v1vd7WRqutc5tFH4ZxzXI9jxpiQC/YqpkBtMMUBrXHdkQ4tahARqa+qv3tPT8ftWZjyat8+mDABXn8dvvgCDhwgrnFjGD4czjvP9UFsjAmrYPcgYvjrVUw7cb3LvamqU4uRZbSIpHnLzwKuKsayTGn2449w6aVub6FRI7jxRjj/fH7Yto0+GRmRTmdMuRHsVUx9wpQDVb0oXMs2pcTevXD//TBqFNSpA+++C4MHu74LAKZOjWQ6Y8qdUJ2DqKWqmwue05g8zJrlLkNduBAuvhieeMI1aW2MiZiYYGYWkStF5Daf5x1FZA2wQURmi0i9kCc0Zd9DD7lObbZtg48/hrFjrTgYEwWCKhDA9YBv11WP4+6evhFIAu4LSSpTfrz9Ntx1F5x1FixY4Dq/McZEhWAPMTUBFgOISBLQGxisqp+KyGbgoRDnM2XZokVw+eXQowe8+iqU0M1sxpjCCXYPogLgtWxGT9wVR1O956uBuqGJZcq8nTvhjDOgShW3F2HFwZioE2yBWAbkHgM4F/hWVf/0njcAtoQqmCnjrr3W9bfw5pvQsGGk0xhjAgj2ENOjwKsicglQA/i7z7QM4KdQBTNl2MSJ7pDSsGFw/PGRTmOMyUOw90G8LiKrgKOBWao63WfyeuDDUIYzZdDmzXDVVdC5M9xzT6TTGGPyEfR9EKr6NfB1gPHDQ5LIlG033OCKxBdf2HkHY6JcsPdB9BCRk32e1xKRN0TkZxF5VEQqhD6iKTPeftu1rTRsmNuDMMZEtWBPUo8C0n2ePwKcBCwFrgbuClEuU9asXQtDh7rOeu68M9JpjDGFEGyBSMV15oOIxAJnATep6pnA3cD5oY1nyoScHNf43t697uR0bGykExljCiHYcxCJuB7fALoDCRzut2Eu7kY6Y470+OMweTK8+CK0ahXpNMaYQgp2D2ItkHvweCCwQFU3eM9rAH8GfJUpvz76CG6/Hc48E/7xj0inMcYEIdgC8QbwoIhMAG4GxvtM64K7ka7U2Lp/f6QjlG1z58K557re38aOBZFIJzLGBCHYAjECeBiIx52wfsJnWmfgneIGEpFbRURFpHZxl5Wf65YuJX3OnNy+sk2orV4NJ58MtWu7vYiEhEgnMsYEKdgb5Q4CD+QxbXBxw4hIY6AfsKq4yyrI0dWq8dy6dXyzfTs9q1cP9+rKl1WrYMAA2LULvvkG6lkr8MaURsHuQQAgIp1E5DoRGZ7bB4SItBSRqsXM8wRwO3/t1jTkTq9dmyoxMby6fn24V1W+/PQTHHssrFvn9hw6dIh0ImNMEUkwh1hEJB533uEMQHAb8m6qOldE3gOWquodRQoicipwgqr+n4hkAV1VdVMe8w4BhgAkJyenjxkzhsTExKDX+QAwE3gXCOc9vdnZ2UXKV5JCkbH6jz/SYdgwDlauzE8PP8yu5s1DlM4pL59juFnG0CjtGTMyMuaoatd8F6CqhR5wjfVtxd3vUAfX9HcXb9qVwI8FvH4ysCDAcBrwPZDkzZcF1C5MpvT0dM3MzNSi+HzzZiUzU9/dsKFIry+souYrScXO+NprqrGxqu3bq65aFZJM/srF51gCLGNolPaMwGwtYPsa7H0Q5wH3qGu0z79ZjRVA0wKKUd9A40WkI9AMmC/uSpdGwFwR6a6qfwSZsdBOqF6denFxvLp+PWfUqROu1ZRtGzfCjTe6JjR69YIPPoAaNSKdyhgTAsGeg6gFLMpnWfFFCaGqP6tqXVVtqqpNgTW4PZOwFQeAijExnF+3Lp9s3sxmu+Q1OKquL4d27eCdd2DECJg0yYqDMWVIsAViBXBsHtO6A0uKF6fkXZSczH5V3t6woeCZjbN2LZx2Gpx3HjRr5u53GD7cWmc1powJtkCMA+4QkQs4fF5XRSQDuAl4KRShvD2JgCeoQ61zYiLtq1RhvF3NVLCDB+H5591ew+TJ8Oij8N13dqWSMWVUsAViNPAJ8CqHuxf9Gnfy+XNVfSaE2UqEiHBRvXp8u2MHv+3eHek40WvaNOjSBa65xt0Z/dNPcMstUMFaeDemrAqqQKjqQVU9F+gNPAaMAZ4GjlfVC8KQr0RcULcuArYXEcjKlXD22dCnD2zb5vp0mDIFWraMdDJjTJgV+iomEYnD3TZwh6p+CcwIW6oS1qhSJTKqV+fVP/7g3pQUxNoMgj//hIcfhtGjXRtKI0fCrbdClSqRTmaMKSGF3oNQ1X24S1EPhC9O5FyUnMxve/Ywc8eOgmcu6yZPhtRUuO8+GDwYliyBe++14mBMORPsOYhJQP9wBIm0M+rUoXJMTPk+zLR7t7unoV8/VwymT4c33oDGjSOdzBgTAcHeKPcMMF5EKgIfAL/j126Sqi4PTbSSVa1iRU6rXZs3N2zgiZYtiYspUjNVpVbismVw7bXwyy9w/fUwapTtMRhTzgW7FZwGNMb1BTEN1xf1Mr+h1LooOZktBw7w2ZYtBc9cVhw8CKNG0eWaa2DrVvj8c3j6aSsOxpig9yAupwRaWo2U/jVqUDc2llf/+IPTaoe1O4rosH49nH8+fPUVm3r3pu6770KtWpFOZYyJEsH2B/FKmHJEhYoxMZxXty7Pr1vH1v37qREbG+lI4fP113DOObBlC7z0Er80bUpdKw7GGB+FOsQkIhki8oGILBCRr0Xk+nAHi5SL6tVjnyoTNm6MdJTwUIUnnnD3NVSpAjNnwmWXWXegxpi/KLBAiMiJuDulewO7gObAkyJyV5izRUSXxETaVqlSNjsS2rHD3fR2881w6qkwezZ07hzpVMaYKFWYPYi7cM1pNFHVo3EnqccA/5QyeEeZiHBRcjIztm8nqyw1vfHrr3DMMfD++64NpXffhaSkSKcyxkSxwhSIVOAxVd0Jh/qlvh+oiisWZc4FyclAGWp6Y9Ik6N4dNmxwj2+5xQ4pGWMKVJgCURvw75fhd+9vmTyrmVKpEr2Tkhi/fn1uT3ilkyo8+SQMGAANG8KsWZCREelUxphSorD3QZTirWTRXFSvHkt272b2zp2RjlI0e/fCFVfATTe5vhu++8713WCMMYVU2ALxkYisyh1wHQcBfOo7XkRWFjWIiNwvIj+JyDwR+VJEGhR1WaFwVp06xIuUzpPVmzdD//7w8suuDaUJEyDKO1c3xkSfwtwHMTbsKZxHVHUYgIjcANwLDC2hdf9FUsWKnOo1vfFYixbElpamN5Ytg0GDXDPdr7/uen0zxpgiKLBAqOplJRFEVX2bUU0gCg5rXZSczDsbN/LFli2cXBrurJ4xw7W+KuL6bOjZM9KJjDGlmETTSVgReQC4GNgOZKhqwLvVRGQIMAQgOTk5fcyYMSSG4RDKAeAcoDrwPIf7WA1WdnZ2WPL5qjt5Mm1Hj2ZPvXr8/NBD7G7YMKjXl0TG4rKMoWEZQ6O0Z8zIyJijql3zXYCq5jsAXQqaJ8BrKgFtA4yfDCwIMJzmN9+dwMjCrCs9PV0zMzM1XD7etEnJzNQbli4t8jLCmU9zclQffFAVVHv1Ut28uUiLCWvGELGMoWEZQ6O0ZwRmawHb18IcWJ8uIh+KyAARyXd+EWni3WG9Ajg5QDHqq6odAgwT/WZ9HTizENnCblCtWtzQsCFPr13LJ5s3RzrOkQ4cgKFD4a673LmGL7+EmjUjncoYU0YU5iR1G9yNcROBHSLyHTAf2AjsBWrgmt/oDnTAFYdbVPX1YIKISCtVzW0u/FRgcTCvD6eHmzdn2rZtXLp4MT917Ur9+PhIR4LsbNfY3qefwp13wr/+BaXlRLoxplQozEnqtcDlInIHcBlwIq4/iMo+s60ApgN3AF94uy/BGiUibYAcYCURvILJX6UKFXijXTvS58zh4sWL+aJTJ2IieSfyjh2u17fZs+GFF+CqqyKXxRhTZhW6uW9V3QA87A2ISHXcuYbNqrq/uEFUNSoOKeUlNSGBp1q2ZMjSpTy2ejW3NWkSmSB79rgrlebMgffeczfBGWNMGBT5mISqblPVP0JRHEqLf9Svz5m1a3PXihXM3rGj4BeE2oEDcO65kJkJY8dacTDGhJUdtA6CiPDfNm2oHxfHeYsWsfPAgZJbeU4OXHklTJwIzzwDF1xQcus2xpRLQRUIEckRkYN5DAdEZLOITBKR/uEKHGk1YmN5LTWV5bt3c/2yEuyC+4UX4JVXYORIuO66kluvMabcCnYP4n5gNe4Kpldw5yPGes/XAK8CdYDPROQvl7mWFcdVr849KSmMXb+eN0qiraa9e+Ghh9yd0cOGhX99xhhDkH1SA3twVywNVNU9uSNFpDLwGa5QdAE+wXU09HGIckadYSkpTNm6laFLl3JMtWo0q1y54BcV1bhxsGYN/O9/1o+DMabEBLsHMRR4wrc4AKjqbuAJYKiq5uB6nOsUmojRqWJMDK+1a4cA5y9axP6cnPCsaP9+t/fQrZu7tNUYY0pIsAWiLhCbx7Q4DncgtAko8z91UypV4sU2bZi5YwfXLFtGTjjatfrhB1ixwnqBM8aUuGAPMc0GRojIt6qa26scXt8Nw73pACnAutBEjG5n163Lz7t28a+VKzmgypg2bagQyg35wYPub506oVumMcYUQrAF4v+AKcAKr8mNDbi9imOBP4ELvfla4tpTKhfub9aMWBGGZ2WxLyeHsW3bUjFUzV5U9P6JSvKSWmOMIcgCoapzRaQlcAtwNNAR1z/1Y8DjqrrZm+/eUAeNdvc2bUqcCHeuWMF+VV5LTQ1NJ0O5BWJ/ubkf0RgTJYLdg8ArAneFIUupd0dKCvExMdz822/sy8nhrfbtiS9ukajlndZZvrz4AY0xJghF2nqJSE0RGSQiF4nIQBGxNqY9NzVuzLOtWjFx82bOWLCAPbnnEIqqRQtITXXtLhljTAkKukCIyL+AtcBHuJvkPgHWisj9Ic5Wal3bsCEvtm7NZ1u2cOqCBewp+CX5O+ssmD4dSuKmPGOM8QTb1MaNuMNL44EMINX7Ox64S0RuCHXA0urKBg14qU0bJm/dyp1AdnFOMp91lmuLyfYijDElqCg3yj2lqleq6jRVXeL9vRJ4Grgm9BFLr0vr12d8aio/AX3nz2ft3r1FW1DHjnDUUe6GuV27QprRGGPyEmyBaIo7pBTIJ970IhGRR0RksYj8JCLve/1NlHrnJyczEliwaxdHzZ7NV1u3Br8QEXj6aVi92hUJY4wpAcEWiM24bkUDae9NL6pJQAdV7QQsBe4sxrKiSk9gVno6tWNj6Td/PqNWrgz+ruuePeHCC+GRR6AkW5E1xpRbwRaI94H7vauXYgFEpKKInAfcB7xb1CCq+qWq5h6onwk0KuqyolFqQgI/dOnC3+vU4c4VKxi8YAFbg723YfRoiI93zX2Ho1kPY4zxIcF0Hy0iVYFPgb8BB4EtQE2gAvA1cJKqZhc7lMhHwFuqOj6P6UOAIQDJycnpY8aMITExsbirDZvs7OxD+RRXZf8NJAMjcbedF1aD99+n9dNPs+y661h7Zuh6afXNGK0sY2hYxtAo7RkzMjLmqGrXfBegqkENuEb4TgZGA//F9QlxEl6xKeC1k4EFAYbTfOa5G7cNLXB5qkp6erpmZmZqNAuU79tt27ThN99opWnT9KV16wq/sJwc1ZNPVo2LU503L6wZo41lDA3LGBqlPSMwWwvYvhblTmrF9fMQdF8Pqto3v+kicgmu+JzgrafMOjYpibldu3LeL79w+ZIlTN22jSdbtqRGbF6N5XpE4KWXoHNn1z/1nDlQpUrJhDbGlCsFnoMooJvRv3Q7WtQgIjIA+Cdwqqr+WdTllCZ14+L4snNn7k1J4bX162k3axbvbdxY8Avr1HGdCC1ZAlddZecjjDFhUZg9iPtwh87D7VkgHpgkrrnsmao6tATWG1EVRBjZrBmDa9fm8iVLOHPhQs6qU4dnW7UiOS4u7xf27Qv33ee6IG3e3PVVbYwxIVRggVDVESWQA1UN5lxtmXNU1ar80KULj6xezcisLL7aupUnW7bkwuRkJK/+Je6+23UmdN99kJICl19esqGNMWVaiDotMKEQGxPDXSkpzO/albZVqnDx4sUM+vlnVu3JozUnEXjhBejf3x1q+vLLkg1sjCnTrEBEobYJCUw/6iieatmSadu20X7WLJ5as4YDgfq9jo2Fd96B9u3hjDPgu+9KPrAxpkyyAhGlKohwQ6NGLOjWjb9Vq8aNv/5K97lz+WHHjr/OXK0afPYZ1K8PAwfC3LklH9gYU+ZYgYhyzSpX5rNOnXi7XTv+2LePY+bO5dqlS9nmfxd2/fowZQokJblDTgsXRiawMabMsAJRCogIf69bl8Xdu3N9w4a8sG4dbX/4gTfWr+eI20WaNIGvvoK4ODjhBFi6NHKhjTGlnhWIUqRaxYo81aoVP6Sn07hSJc5ftIgTf/qJFbt3H56pRQu3J5GTA717256EMabIrECUQulVqzKzSxeeadmS73bsoMOsWTyxejUHc/cmUlNh2jR3lVOfPvDjjxHNa4wpnaxAlFIVRLiuUSN+6daNjOrVufm33+gxdy4/Z3ttJaamum5Kq1SB44+H77+PbGBjTKljBaKUa1ypEh917Mjrqaks37OHLnPmcM/y5ew4cABatnRFolYtd+d1Zmak4xpjShErEGWAiHBecjKLunXjvLp1eWDVKprOnMm/srLY3rChKxIpKXDiifDaa5GOa4wpJaxAlCG14+IYl5rKrC5d6JmUxLCsLJrOnMl9+/axfepU+NvfXK90Dz5oDfwZYwpkBaIM6lqtGh927Mjs9HR6JSUxPCuLlEWLGPn882w75xzXhtPQoXCgyI3vGmPKASsQZVh61apM7NiRuenpZNSowYg//qDmkCE8e/HF8OKL7G7cmL033OCa57A9CmOMn6A7DDKlz1FVq/J+hw6s3rOHz7ds4bPbbmN6airnffIJA194AZ55ho7JyWw45xzqXHgh0rWru0TWGFOuWYEoRxpXqsSVDRpwZYMG7G/Xjm+HDuWhrCwOTJzIsZ9/Tv9//xt5+mnWN2nClsGDaXDRRSSlp1uxMKacipoCISJ/B0YAqUB3VZ0d2URlW2xMDL2rV6d3WhqkpTGhd2/eTk5mx4QJtPnoI3o/+ywVn36arKZNWXnyydQ85RRS+/ShYn6dGBljypSoKRDAAuAM4D+RDlIe1Qb6pKbCsGEcuPtuflyxgvVvvEH9Dz6g57//TYVnn2V7YiKLjjmGA/360fzUU2nQtm2kYxtjwihqCoSqLgLy7j3NlJiKMTF0a9EC7rkH7rmHrevXs+jjjznw+ee0mDGDhpMnwz//SVbjxqzq1YvKAwaQOmgQiTVqRDq6MSaEoqZAmOhVIzmZHldcAVdcgebksHTePNZ+9BFVpkwh/d13SXjtNfZXqMD8Tp3YdsIJ1B00iDY9exJT0f57GVOaiZbg5Y0iMhmoF2DS3ao60ZtnKnBrfucgRGQIMAQgOTk5fcyYMSQmJoYhcWhkZ2dHdT4oesb9+/ax5ZdfqDRrFi1mzaLDsmUAbKlWjZ+OOopNXboQ37UrVRs0iFjGkmQZQ8MyhkZ+GTMyMuaoatd8F6CqUTUAU4GuhZ0/PT1dMzMzNZpFez7V0GVcv3q1znjuOZ0+eLCuq1NH1d1hoasaNNDpZ5+tM8eM0W2//x7RjOFkGUPDMoZGfhmB2VrA9tWOAZiQqtuoEXWvuQauuQbNyeHXefNY/emnVP7qKzp98glJb79NjgiL2rZlQ69eJJ14Im3796dSQkKkoxtj/ERNgRCR04FngDrAJyIyT1VPjHAsUwwSE0PLLl1o2aUL3HMP+/btY97UqWz57DNqTptGjzFjiP3Pf9gdF8fco45iZ58+JA8YQGs7f2FMVIiab6Gqvg+8H+kcJnzi4uJI69/f9ZkN7Ni2jR8//ZQ/J02i4fTpdHn4YXj4YbYnJrKkWzf2ZmTQ6KSTaHrUUUiMtQpjTEmLmgJhyp9q1avT/fzz4fzzAVi/ahXLPv2UnMmTafrddzTJzIR772V9rVr8dswxbGrXjnX169OgTZsIJzemfLACYaJGcpMmJA8dCkOHoqpkLVrEys8+o8JXX9H6u+/o8ckn8MgjrGrYkJU9elDhhBNoddJJ1GncONLRjSmTrECYqCQiNG3Xjqbt2sEtt5CTk8O7Y8dSa9UqKk+dSqdPPyXpnXcA+K1pU9b26EHcCSfQeuBAatavH+H0xpQNViBMqRATE0OtZs3oc9llMHw4B/bt45dvvmH9l1+SOGMGXd5/n8TXXwdgacuW/NGjB5VPOIHWAwaQVLduhNMbUzpZgTClUsW4ONplZNAuIwOAfXv3Mn/aNDZPmkS1GTPo+tZbVBk3jhwRlrRqxR9/+xtVcgtGrVoRTm9M6WAFwpQJcfHxdPa5QmrPn3/y47RpbJ00ieozZnD0+PFUevllDsbE8EubNmw89liqZGTQ+sQTSapTJ8LpjYlOViBMmVSpShWOGjgQBg4EYPeuXfz41VdsmzyZ6t9+y9Gvvkqll15yexgtW7K+Rw8q5xaMeoFagzGm/LECYcqFygkJHHXKKXDKKYArGHOnTmXblCkkff01Xd98kypjxwKwrEULfj/2WOIzMmjZvz+1GjWKZHRjIsYKhCmXKick0GXQIBg0CIA9u3e7Q1JTplDt669JnzCBhPHjAXeV1Lqjj6Zi794079eP5JYtIxndmBJjBcIYoFLlyhw1YAAMGAB4J71nzGDLpEkkfPstnT7+mKS33gJgdf36rOreHXr1okm/fjRq397u9DZlkhUIYwKIi4+nc9++0LcvAAf272fhzJmsnzKFyl9/TZupU6k9cSIAG2rWZHm3buzv2ZNt9euTc9xxxFSoEMn4xoSEFQhjCqFibCztjzuO9scdB0BOTg5Lf/yRtZMnU3HGDFrMmkWDL74AYNuNN7KsSxf+7NGDWn360LpXL+IqVw5Lrr27dzP7lVeo9cQTpKxcSfz+/eyNjWVlSgpbbr6Z9EsuIT5M6zZlnxUIY4ogJiaG1unptE5Ph3/+E83JYeXSpcwcN476S5fSaPZsuk2fDqNGsScujvnt27PtmGOo2rs3LU84gWq1axc7w/RHHiFt5Eg65eRQdffuQ+Mr79tH22XL2Hnzzey99Va+Hz6cXrfdVuz1mfLHCoQxISAxMaS0bcuK/v3p9eCDAKxfs4blU6awZ9o06vzwAz1efJHY558nR8Td7d29O7HHHUezvn2p16JFUOvLvP56uv/3vyTs3ZvnPLlFI33YMDJXrSLjmWeK/gZNuWQFwpgwSW7UiORLLoFLLgFg5/btzM/MZMe0aVSbOZMu771H4muvAbCmXj1WpadzsEcP6mdk0Lxbtzz7xJj+yCMFFgdfCXv30v2//2V6kya2J2GCYgXCmBJSNSmJroMHw+DBgOvPe+F337ExM5P4b7+l+cyZ1PvkEwB2JCSwrFMnso85hqRevWidkUGVpCT27t5N2siRhS4OuRL27iVt5Ej2XX99iN+VKcuirkCIyADgKaACMEZVR0U4kjFhERsXR/vevaF3bwB3HmPJElZNmcLBb76h/uzZHPXkk8Q88QQHYmJY1Lo1u2rUoP3+/UVan+TkMHvsWLD+NEwhRVWBEJEKwHNAP2ANMEtEPlTVXyKbzJjwk5gYUlJTSUlNheuuA2Drxo0sy8wke/p0qn//PUfNnEkF1SItv+ru3dR8/HH++M9/QhnblGFRVSCA7sCvqrocQETeBE4DrECYcqlGnTp0P/tsOPtsAHbHx1N5374iL6/JypX8EapwpswTLeKvkXAQkbOAAar6D+/5RcDRqnqd33xDgCEAycnJ6WPGjCExMbHE8xZWdnZ2VOcDyxgq4c7Y6/jjiSnGdzZHhE8//LDcf46hUNozZmRkzFHVrvm9Ptr2ICTAuL98G1T1ReBFgK5du2piYiJ9+vQJc7Simzp1alTnA8sYKuHOuDs2tlh7EHvi4oj27wvYv3WoFDdjtDUgswbw7WC4EbAuQlmMiTorU1KK9fpVxXy9KV+irUDMAlqJSDMRiQPOBT6McCZjosaWm29mZxGbzthZpQpbbr45xIlMWRZVBUJVDwDXAV8Ai4C3VXVhZFMZEz3SL7kELWLLsSpCV++mPWMKI6oKBICqfqqqrVW1hao+EOk8xkST+MqVmTd8OLvi44N63a74eOYNH05cpUphSmbKoqgrEMaY/PW67TZ+uPLKQheJXfHx/HDlldbMhgmaFQhjSqGMZ55hzv33syMhIc9zEjsrV2ZHQgJz7r/fGuozRWIFwphSqtdttxG/cSM/P/44i1u35s/4eHJE+DM+nsWtW/PzE09QadMm23MwRRZt90EYY4IQX7kyPYYOhaFDD42rArSNXCRThkTVndRFISIbgV3ApkhnyUdtojsfWMZQsYyhYRlDI7+MKapaJ78Xl/oCASAiswu6ZTySoj0fWMZQsYyhYRlDo7gZ7RyEMcaYgKxAGGOMCaisFIgXIx2gANGeDyxjqFjG0LCMoVGsjGXiHIQxxpjQKyt7EMYYY0LMCoQxxpiAor5AiMhLIrJBRBb4jEsTkZkiMk9EZotId59pd4rIryKyREROLKGMjUUkU0QWichCEfk/b3xNEZkkIsu8vzUikTOffI+IyGIR+UlE3heR6pHIl19Gn+m3ioiKSO1ozCgi13s5ForI6GjLGE3fGRGpJCI/iMh8L+NIb3xUfF8KyBhN35mAGX2mF/87o6pRPQC9gC7AAp9xXwIDvccnAVO9x+2A+UA80Az4DahQAhnrA128x1WBpV6W0cAd3vg7gIcjkTOffP2Bit74hyOVL7+M3vPGuCbgVwK1oy0jkAFMBuK9aXWjMGPUfGdwPUcmeo9jge+BY6Ll+1JAxmj6zgTM6D0PyXcm6vcgVHU6sMV/NFDNe5zE4V7nTgPeVNW9qroC+BXoTpip6u+qOtd7vBPXl0VDL89Yb7axwOBI5Mwrn6p+qa4PDoCZuB78Sjxffhm9yU8At3Nk97PRlPFqYJSq7vWmbYjCjFHznVEn23sa6w1KlHxf8ssYZd+ZvD5HCNF3JuoLRB5uBB4RkdXAo8Cd3viGwGqf+dZweCNTIkSkKXAUrponq+rv4L64QF1vtojl9Mvn63LgM+9xRD9H34wiciqwVlXn+80WNRmB1sBxIvK9iEwTkW5RmPFGoug7IyIVRGQesAGYpKpR933JI6OviH9nAmUM5XemtBaIq4GbVLUxcBPwP2+8BJi3xK7jFZFE4F3gRlXdkd+sAcaFPWde+UTkbuAA8Fok83lZDmX0Mt0N3Bto1gDjSjyj9zlWBGrgDkHcBrwtIhJlGaPqO6OqB1U1DfcLvLuIdMhn9qjLGC3fmQAZOxHC70xpLRCXAO95j9/h8G7SGtyxt1yNOLwrHVYiEov7Qr6mqrnZ1otIfW96fVyVj0jOPPIhIpcAJwMXqHegMhL58sjYAnesdL6IZHk55opIvSjKiJflPW+X/wcgB9dIWjRljLrvDICqbgOmAgOIou9LPhmj6jsTIONphPI7E86TKKEagKYceZJ6EdDHe3wCMMd73J4jT8Isp2ROUgswDnjSb/wjHHnSbXQkcuaTbwDwC1DHb3yJf455ZfSbJ4vDJ9yiJiMwFLjPe9watxsvUZYxar4zQB2guve4MjADt8GNiu9LARmj6TsTMKPfPMX6zoQtfAg/hDeA34H9uAp4BdATmOO92e+BdJ/578adnV+Cd9VGCWTsidtV+wmY5w0nAbWAKcAy72/NSOTMJ9+v3sYsd9wLkfoc88roN8+h/+zRlBGIA8YDC4C5wPFRmDFqvjNAJ+BHL+MC4F5vfFR8XwrIGE3fmYAZ/eYp1nfGmtowxhgTUGk9B2GMMSbMrEAYY4wJyAqEMcaYgKxAGGOMCcgKhDHGmICsQJiw8FqRLGjIinTOkiYiTUVkhIg0j3SWQETkTBFZLyJVfMZlich4v/liRGSsiOSIyJXeuJu8Vk5tu1JGVIx0AFNmHev3/H3cNfgjfMbtLbE00aMpMBz4GnejUtQQkYrAg8AjqvpnAfO9CvwduFRVx3mTXgD+ibtr++UwxzUlwAqECQtVnen7XET2Apv8x5cFIhKvXkuupTzDabgC9lI+64kF3gROBc5X1bdzp6nqbhEZB9yKFYgywXYFTcSISDMReU1ENorIXnGd2ZzuN88I73BUWxH5QkR2icgqEbnMm36R14FLtriOclr4vT5LRMaLyJVeRyl7RGSuiGQEyNNbRKaIyE5vPV/4NyInIlNF5GsROUVEfvQK3zXetOtE5DsR2SIi28R10DPI57V9gEzv6SSfQ219vOkqIiP81tfUG3+pz7hXRGSNiBwrIt+KyG5cXwqISG0ReV5E1nqf6WIRGVLIf5J/AJ+rqn/z+rnrjce153Qy8Hff4uDjTaCdiPQo5DpNFLMCYSJCRBrjmnzojGtd9FRcMxXvimuu2N87wCe4PgLmAC+JyIO4VkrvAC4D2gCvB3htb+BmXDMD5+IObX0mIm188gzCNe+QDVwInI/rcGeGl9VXa+Bp4BngRO914H59j8EdejkHmA18LCIDvelzgWu9xzfgDsMd640PVhJuY/wGMBB4XUSqAd8Ag3CH8gYBHwHPi8j1+S3M2/j3wbXnE0hl4EOgLzBYVT/IY755wA68hu1MKRfu9kJssEH1UJsw432e/w/YCNTym28SMM/n+Qhc20IX+4yrgWtqeTNQzWf8Dd68KX7r3Qc08RlXFdcJ1as+434FpvhlqQZswqfhO1yLmTlAWgHvNwZ3CPdLYKLP+D5exr4BXqPACL9xTb3xl/qMe8Ubd5rfvMOAPUArv/H/9d5HxXzyHu0ts18e/3bqDZcX4t96BvBlpP/P2VD8wfYgTKQMAD4FtotIxdwB101iZ+/XsK/cjllQ1a24pqBn6pH9biz2/vr/4p+pqqt8Xr8TtzdyLICItMI1Lf6aX5Y/ge9w3d76ylLVef5vSETSReRjEVmPK2D7gX64PZtQOwB87DduAG6vbEWAz7QWrsvJvDTw/m7MY/oPuCIzXERSCsi20Wd5phSzAmEipS5wMW4j6js84k2v5Tf/Vr/n+/IYB1DJb/z6AOtfz+HetHJ7LvtfgDwnB8jyu//CvMNQU4CawPVAD6Ab8HmAPKGwQVUP+o2riytm/u/hHW+6//vwlZsxrxPdy3CHl6oCU0Qkv57IduMOSZlSzq5iMpGyGXco4uE8poeys5XkPMat9ckCrhvOyQHm3ef3PFATyANw5wXOVtU1uSN97ycohL24psN95bVRD5RhM27P6v/yeM2SfNad+xnUyGsGVZ0vIv1xhXCKiPRW1UDFtyZub8OUclYgTKR8jjvEs1BVd4d5XceISGNVXQ0gIlVxJ3A/8aYvwR1nb6+qo4q4jtxCsD93hIi0Bv6G68ckV+4v9EC/sFcC/l1vDgowX14+x+29rFLVDQXN7Cf38Fxz4Nu8ZlLV2d5J9y+AySKSoar+xaAZ7pCUKeWsQJhIuRe3EZkuIs/iNtA1cBvI5qp6eQjXtR740ruEdC/uZq4E4H4AVVURuRaYKCJxwNu4X8DJuENFq1T18QLWMRl3XmCciDwG1AdGAqs48lDuUm++y0Vki5dniXde5E3gHnH9Hc8EjgPOC+J9PoG7emqGiDyBK3wJQFvgOFU9La8XquoqEVmJ64p0fF7zefN+KyIn484hfSkiJ3jnhRCR6rirvB4NIreJUnYOwkSEd9K4K+7u6gdxVy89j7sk9asQr24a8Ji3nrdwx9sHqupSnzyf4o7fJ+AuVf0Cd29BPdyJ6nyp6kLgAiAFdzno7bjLb6f7zbcZuA53ee80YBaQ7k1+CHjWm/4BkApcVNg3qarbcQXtU1wR/AJ309tpHL7/Ij9v4c65FGZd03CXHKcCn/tcVDAId0ju/cLmNtHLepQzZZq49p6+VtULI50l2nk3GS7B9V39dRGX8RnujvlCFzYTvewQkzEGAFX9TURexu35FGpPwpeIpAEZ/PU8iiml7BCTMcbXMGBWkFdf5aoHXKaqv4Y4k4kQO8RkjDEmINuDMMYYE5AVCGOMMQFZgTDGGBOQFQhjjDEBWYEwxhgT0P8DXp9PQj4J/MoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot( col1.Tatm, np.log(col1.lev/1000), 'c-', label='RCE' )\n",
    "ax.plot( col1.Ts, 0, 'co', markersize=16 )\n",
    "ax.plot(col2.Tatm, np.log(col2.lev/1000), 'r-', label='RCE O3' )\n",
    "ax.plot(col2.Ts, 0, 'ro', markersize=16 )\n",
    "ax.invert_yaxis()\n",
    "ax.set_xlabel('Temperature (K)', fontsize=16)\n",
    "ax.set_ylabel('log(Pressure)', fontsize=16 )\n",
    "ax.set_title('Temperature profiles', fontsize = 18)\n",
    "ax.grid()\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we include ozone we get a well-defined stratosphere. We can also a slight cooling effect in the troposphere.\n",
    "\n",
    "Things to consider / try:\n",
    "\n",
    "- Here we used the global annual mean Q = 341.3 W m$^{-2}$. We might want to consider latitudinal or seasonal variations in Q.\n",
    "- We also used the global annual mean ozone profile! Ozone varies tremendously in latitude and by season. That information is all contained in the ozone data file we opened above. We might explore the effects of those variations.\n",
    "- We can calculate climate sensitivity in this model by doubling the CO2 concentration and re-running out to the new equilibrium. Does the amount of ozone affect the climate sensitivity?  (example below)\n",
    "- An important shortcoming of the model: there are no clouds! (that would be the next step in the hierarchy of column models)\n",
    "- Clouds would act both in the shortwave (increasing the albedo, cooling the climate) and in the longwave (greenhouse effect, warming the climate). Which effect is stronger depends on the vertical structure of the clouds (high or low clouds) and their optical properties (e.g. thin cirrus clouds are nearly transparent to solar radiation but are good longwave absorbers)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (1,) \n",
      "  Tatm: (59,) \n",
      "The subprocess tree: \n",
      "Untitled: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "col3 = climlab.process_like(col2)\n",
    "print(col3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's double CO2.\n",
    "col3.absorber_vmr['CO2'] *= 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The radiative forcing for doubling CO2 is 1.310158 W/m2.\n"
     ]
    }
   ],
   "source": [
    "col3.compute_diagnostics()\n",
    "print('The radiative forcing for doubling CO2 is %f W/m2.' % (col2.OLR - col3.OLR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1095 steps, 1095.7266 days, or 3 years.\n",
      "Total elapsed time is 4.996684391891189 years.\n"
     ]
    }
   ],
   "source": [
    "col3.integrate_years(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Field([3.98517386e-07])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col3.ASR - col3.OLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Equilibrium Climate Sensitivity is 2.758433 K.\n"
     ]
    }
   ],
   "source": [
    "print('The Equilibrium Climate Sensitivity is %f K.' % (col3.Ts - col2.Ts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Ts: (1,) \n",
      "  Tatm: (30,) \n",
      "The subprocess tree: \n",
      "Untitled: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "col4 = climlab.process_like(col1)\n",
    "print(col4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'CO2': Field([0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038, 0.00038,\n",
       "        0.00038, 0.00038]),\n",
       " 'O3': Field([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "        0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
       " 'H2O': Field([5.00000000e-06, 5.00000000e-06, 5.00000000e-06, 5.00000000e-06,\n",
       "        5.00000000e-06, 5.00000000e-06, 6.38590233e-06, 9.08848690e-06,\n",
       "        1.33273826e-05, 2.34389689e-05, 3.84220914e-05, 5.95564299e-05,\n",
       "        8.82144990e-05, 1.25843839e-04, 1.73951159e-04, 2.34088411e-04,\n",
       "        3.07840683e-04, 3.96815735e-04, 5.02635028e-04, 6.26926041e-04,\n",
       "        7.71315753e-04, 9.37425100e-04, 1.12686431e-03, 1.34122899e-03,\n",
       "        1.58209684e-03, 1.85102493e-03, 2.14954752e-03, 2.47917415e-03,\n",
       "        2.84138824e-03, 3.23764591e-03])}"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col4.absorber_vmr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The radiative forcing for doubling CO2 is 4.421081 W/m2.\n"
     ]
    }
   ],
   "source": [
    "col4.absorber_vmr['CO2'] *= 2.\n",
    "col4.compute_diagnostics()\n",
    "print('The radiative forcing for doubling CO2 is %f W/m2.' % (col1.OLR - col4.OLR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1095 steps, 1095.7266 days, or 3.0 years.\n",
      "Total elapsed time is 4.996684391891189 years.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Field([-5.25654883e-07])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col4.integrate_years(3.)\n",
    "col4.ASR - col4.OLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Equilibrium Climate Sensitivity is 3.180993 K.\n"
     ]
    }
   ],
   "source": [
    "print('The Equilibrium Climate Sensitivity is %f K.' % (col4.Ts - col1.Ts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interesting that the model is MORE sensitive when ozone is set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}