{
 "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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Tatm: (30,) \n",
      "  Ts: (1,) \n",
      "The subprocess tree: \n",
      "top: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'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.        ]), 'Ts': Field([ 288.])}"
      ]
     },
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 730 steps, 730.4844 days, or 2 years.\n",
      "Total elapsed time is 1.99867375676 years.\n"
     ]
    }
   ],
   "source": [
    "col1.integrate_years(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.00150116])"
      ]
     },
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEfCAYAAABmsjC7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8XFX5x/HPN3vbdEvTpkvSfU33Ni1FBMsmFfGHIGoV\nAZVFEBQEREBAFFAQRFFEBERBkMqigmyySAG1+76nS9qme2lL23Rv8/z+uDc6ZGkyzSQzkzzv12te\nM3PumTvPnGnnyT333HNkZjjnnHP1lRLvAJxzzjUNnlCcc87FhCcU55xzMeEJxTnnXEx4QnHOORcT\nnlCcc87FhCcU51xcSJogaaakMkkm6TxJZ1U8jqg3JCy7Pp7xutp5QnH/Ff6nreutZ7zjTVSSxkq6\nXVLXeMeSqCR1Bl4ABHwHuACYHtegXL2lxTsAl1AuqPT8ROAy4BHg/UrbtjZKRMlpLPAD4GVgQ5xj\nSVQfA1oCN5nZGxWFktYBLYCD8QrMHTtPKO6/zOypyOeS0ggSypTK25oLSa3NbHe844iUaDFJyjaz\nsihf1jm83x5ZaGblwP6YBOYanXd5uXqRlCrpGklzJe2TtEvSm5JOqFTvv/3gki6QtFDSfknLJH0p\nrNNH0ouSPpS0U9LjklpW2s/zYZ97F0mTJO2QtEfS65KG1BDjhZKmhq/bI+nfkj5TqU52GN+Dks6U\nNEXSHuCZcHsPSQ9Imh/Gt0/SgvCzp0Ts5z7gV+HTGRFdhA9WbA+f51YT5weSXq5rTGGdDpLul7RK\n0kFJmyU9Kamgjt9fnduz0nd4Yfid7wfujqhTJOnlcD/7wza6WpIiPyfw60ptVBZuq3IOpZb4a/1u\nw3rnhtu2SdorabWk57zrNrb8CMUds/BH4nngM8Akgq6xlsBXgcmSPmVmb1V62ReBLsDDwE7gcuBp\nSUeAXxB0E90InAB8DdgFXFNpH6nAW8Bq4FagALgSeF/SGDNbERHjL4CrgZf43w/xF4CXJH3NzP5Q\nad8nAReG8T0OHArLxwBnAi8CK4Gs8HP/PHz/68J6fwI6hvu4DSgJy5dVacC6qzamMClNDd/vd8DS\nMJYrgNMkjTazjXXYf53bM/QVoBvwG4LEsC2M50TgTaAMeJCgW/Qcgu+1EPhG+PrLCdouso0OEaW6\nfreSziT4dzoLuBPYHcb/SaB7+LldLJiZ3/xW7Y0gMRjw1Rq2XxBu/3Kl8kxgEbAwomxIWPdDoHNE\neQFwGCgHLqu0nzeAvUBaRNnz4X7+WKnuiWH58xFlJ4VlN1eqK4Ifvg+AzLAsO6xrwPHVfNaWgKop\n/ytwAGgfUXZVuJ+iaurfF27LrWbbB8DLEc9ri+l3BD+OAyqV9wf2AQ/W4TuOpj0rvsN9QK9q9jWP\noLuqf0RZCsEfCQYcV1sbAWeF5edV877XH+N3+whwBGgd7/9TTf3mXV6uPr4CbAHekJRbcQNaA68A\ng1V1pNOfzWxTxRMzKwXWEPxIPVap7vsEJ2jzq3nveyKfmNn7wH+AT0tKD4vPJ0hWT1eKrwPBX7Ud\ngNGV9vtvM5tS+c3MbK9V/GJJmZJywn29DmQAI6qJMVaqxKTg/NYXCX48t1X6fNuB2QR/gddVXdqz\nwgtmVhJZIKk3MAx41syKI/ZTDvwkfHpOFPHUJprvdidBYjtHUmoMY3CVeJeXq49BQCeOPuIrj4+O\ndFpVTZ0dwMHwx6dyOQQ/DqsjysupvgtpMcHooa4ESWoQwb/x1dXUjYwvUnF1lSRlAt8n+CHrXU2V\n9kd5j/qqLqbuQCuCH+mafqh31XH/dW3Po8XTK7xfVM22irLq2u1YRfPd3g98CngC+KWk94HXgElm\ntr2mF7voeUJx9SFgLXDxUepU7n8/UkO9msor3udYiKA76qyj1Jlf6fneGur9huCczh+B2wmS6GGC\ncz23U/cBLkdbgKim/4/VxVTRJn8HflnD6w7XMaZo1dRGjanO362ZbZQ0AhgPnEbQXfYg8ENJnzSz\nOQ0ca7PhCcXVx3JgHPCemTXmdQMpwABgYaXyQoI+/IojouXAx4FlYdfaMQlHcX0ZeNXMLqy0rbqu\nrqMljYq/iHMI+vkr9pMDtI0irHUE3YStrOrAh2jVtT2PpuLIc3A12wor1YmFqL5bMztMMPDgLQgu\nPgWmEQwA+WIM42rW/ByKq48nCc5x/LC6jZIqdyfF0vcqvdeJBN0zr5pZxYihJ8P7uyOH9h5jfEeo\ndKQkqR3wrWrqVlyTkVPNtoruotMqlV9XueLRmNkB4FngFEkTqqsjqVMUu6xLex4tnhKCI4IvSOob\nsZ8Ugh9tCAYwxEqdv9vqhmgTJM+DVP8duWPkRyiuPp4g6Ju+UdLxBP3S2wlGbp1IcO5jWAO8735g\ntKRXgFfD97uK4ORrxY8XZjZZ0r3Ad4FBkv4CbCI4JzCW4C/cdrW9mZmVS/orcL6kJwkGC3QBLgE2\nE5zPiDQtvL9dUj5BF1Gxmc0m6KJaC9wnqRvBkcbJBH/ZR3ux4nXAccArkp4hmLrkCNCToCvobYJ2\nqU2d2rMOriIYJDBF0kMER2DnEHy+R8xs2tFeHI0ov9tnJLUiaI+1BKPnzicYTPFk5X27Y+cJxR0z\nMzNJEwl+RL4O3ELwb2ojMBN4oIHe+gjBX/g/J7iuIAP4F8Gw0uWVYrxB0lSC6yquJzii2kzw1/R3\nonjPKwiS5TkE1zqsDt9/BcGoosj3XCrpinD/DwPpBNdrzDazg5I+TdA21xL8mL9C8KO7OIp4MLNt\nko4j+FH9HHAewV/d64DJBNes1EWd27OWeN4Pj2xuJ7g+pAVB19Q11Hye55hF8d0+TjAi8etALsHQ\n9YXAZ8zsZVzMKBwJ6VxSkPQ8MMHMsuMdS1Pg7eliyc+hOOeciwlPKM4552LCE4pzzrmY8HMozjnn\nYqJZjfLKzc21nj17xjuMRrFnzx5atWoV7zASirdJVd4mVXmbVDVr1qwPzKxjbfWaVULp2bMnM2fO\njHcYjWLy5MmMHz8+3mEkFG+TqrxNqvI2qUrSmtpr+TkU55xzMeIJxTnnXEwkdUKRNEHBErIrJEUz\nRYRzzrkYS9qEEi6U82uCuaQKgS9JKjz6q5xzzjWUpE0oBBPArTCzVeHU6ZOAs+Mck3PONVvJPMqr\nGxC5DsI6gplXP0LSZcBlAHl5eUyePLlRgou3srKyZvNZ68rbpCpvk6q8TY5dMieUOjGzR4BHAIqK\niqy5DAf0oY9VeZtU5W1SlbfJsUvmLq/1BOs2VMgPy2LuwXXreH7LFg6VV17y3DnnXIVkPkKZAfST\n1IsgkUwkWKY1psrN+O3GjSzcs4cuGRlc2qULl3bpQn5WVqzfyjnnklrSHqGEa0RfBfwDWAI8a2aL\nYv0+KRJzi4r4+5AhjMzO5o41a+g5dSrnLlzIm9u3U+5zoTnnHJDcRyiY2asES5Y2qFSJs3JzOSs3\nl5J9+/jthg38btMm/vrBB/Rr0YLLu3blq507k5Oe3tChOOdcwkraI5R46dWiBXf36cO644/nqUGD\n6JieznUrV9JtyhS+tnQpM3btineIzjkXF0l9hBJPmSkpnJ+Xx/l5ecwrK+M369fz1ObN/GHTJkZn\nZ/PNbt2Y2KkTLVNT4x2qc841Cj9CiYHh2dk8PGAAGz72MR7s14/95eVcvGwZ3aZM4TsrVrBs7954\nh+iccw3OE0oMtUlL48pu3VgwZgzvjhjBhJwcfr1+PQOnT+f0efN4dds2P4nvnGuyvMurAUjipHbt\nOKldOzYfPMjvNm7kofXr+fSCBQxs2ZLv5OdzQV4eLbw7zDnXhPgRSgPLy8jg5h49KBk3jqcHDaJl\nSgrfKC6mYMoUbi0pYdOBA/EO0TnnYsITSiNJT0nhy3l5zBw9mndHjODjbdty15o19Jg6la8tXcr8\nsrJ4h+icc/XiXV6NLLI7bPnevTywbh2/37SJP2zaxGnt23Ntfj5n5OSQIsU7VOeci4ofocRRv5Yt\nebB/f0qPP56f9OrF4j17OHPBAobMmMGjGzaw78iReIfonHN15gklAeSkp3NjeJ7lqUGDyEpJ4bLi\nYrpPncptfp7FOZckPKEkkIzwYslZo0czecQIPtamDXeGc4ddUVzMqn374h2ic87VyBNKApLEJ9q1\n48WhQ1k6diwXde7M4xs30n/aNC5YsoRFe/bEO0TnnKvCE0qC69+yJb8dMIBV48ZxdX4+f926lSEz\nZvDZBQuY5vOGOecSiCeUJNEtM5Of9e3LmuOP5wc9evDezp2Mmz2bU+fO5e0dOzC/At85F2eeUJJM\nh/R0bu/VizXjxnFv794s2buX0+bNY9zs2bz4wQc+tYtzLm48oSSp1mlpXN+9O6uOO46H+/dn66FD\nfHbhQobNmMFTmzbhA46dc43NE0qSy0pN5Rtdu1I8dixPDxqEJC5YupQLgIfXr+dAeXm8Q3TONROe\nUJqItHBql3lFRbw4ZAjtgCuWL6fftGn8dsMGDnpicc41ME8oTUyKxP/l5vJr4I1hw8jPzOTy4mL6\nT5vGYxs2cMgTi3OugXhCaaIEnJ6Tw79HjuS1oUPJy8jg0uJiBkyfzuMbN3picc7FnCeUJk4SEzp0\nYOqoUbwydCgd0tO5eNkyBk2fzhObNnHYE4tzLkY8oTQTkjizQwemjxrFS0OG0CYtja8uXUphxagw\nH27snKsnTyjNjCQ+k5vLrNGj+duQIbRMSeGCpUsZPH06f9q82ROLc+6YeUJppiRxdm4us4uKeGHw\n4GBiyiVLGDpjBi9s3epX3jvnouYJpZlLkTi3Y0fmFhXxXGEhAs5btIjjZs/mnzt2xDs851wS8YTi\ngCCxnNepE/PHjOEPAwey+eBBTp03j0/Om8es3bvjHZ5zLgl4QnEfkSpxUefOLBs7lp/36cPs3bsp\nmjWLLy5aRPHevfEOzzmXwDyhuGplpaZyTUEBq8aN47YePXhl2zYKp0/n8mXL2OArSDrnquEJxR1V\nm7Q0ftirFyvHjeOb3brx+KZN9J02jRtXrmTHoUPxDs85l0A8obg6ycvI4Jf9+rFs7Fg+17EjPy0t\npfe0adyzdi17j/jcxs45TyguSr1atOCPgwYxt6iIE9q04cZVq+g/bRrPbdniQ42da+YSIqFIKpD0\njqTFkhZJujosz5H0pqTl4X37iNfcJGmFpGWSzohf9M3TsOxsXh42jHdHjKBjRgZfWLyYCfPns9xP\n3DvXbCVEQgEOA9eZWSEwDrhSUiFwI/C2mfUD3g6fE26bCAwGJgAPSUqNS+TN3Ent2jFj1Cge6NuX\nqbt2MWTGDG4tKWGfd4M51+wkREIxs41mNjt8vBtYAnQDzgaeCKs9AXw2fHw2MMnMDphZCbACGNu4\nUbsKaSkpfDs/n6Vjx/L5jh25c80aBs+YwcsffBDv0JxzjUiJ1u8tqSfwHjAEWGtm7cJyATvMrJ2k\nB4GpZvZUuO13wGtm9nw1+7sMuAwgLy9v9KRJkxrlc8RbWVkZ2dnZcXnvucAvgDXACcBVQOe4RPJR\n8WyTROVtUpW3SVUnn3zyLDMrqq1eWmMEU1eSsoEXgGvMbFeQQwJmZpKizn5m9gjwCEBRUZGNHz8+\nRtEmtsmTJxOvzzoe+GZ5Ob9Yt44frl7N14FbevTguoICMlPid1AczzZJVN4mVXmbHLuE6PICkJRO\nkEyeNrO/hMWbJXUJt3cBtoTl64GCiJfnh2UuQWSkpHBD9+4sHTuWM3Ny+H5JCcNnzOBtnx/MuSYr\nIRJK2J31O2CJmd0fsekl4KLw8UXAixHlEyVlSuoF9AOmN1a8ru4KsrJ4fsgQXhs6lMNmnDZvHtev\nWOFr3DvXBCVKl9cJwAXAAklzw7KbgbuBZyVdTNAl/wUAM1sk6VlgMcEIsSvNzIcVJbAJHTqwoF07\nrl+5kp+tW8f7O3fyTGEhvVu0iHdozrkYSYiEYmb/IlgGvTqn1vCau4C7GiwoF3MtUlP5df/+nNK+\nPRcvXcrImTN5bMAAPt+pU7xDc87FQJ0TiqRM4HiC60S6Ai2AD4BlwHtmtqpBInRNzuc6dmR0djYT\nFy/mC4sXc/mHH3J/nz60SPVLiZxLZrUmFEl9gWuA84G2QDmwE9gH5ABZgEmaBTwEPGlm3kHujqpn\nixa8P3Ikt5SU8NPSUv69cyd/LixkUKtW8Q7NOXeMjnpSXtKvCc5TjAF+FN5nmVkHM8s3s5ZAF+Bc\ngssP7gcWSTquYcN2TUF6Sgr39OnDq0OHsvHgQYpmzeL3Gzf6nGDOJanaRnl1Bcaa2XFm9nMzm2Vm\nhyMrmNlmM3vRzC4jSC6/AYY3ULyuCfpUhw7MKyriuDZt+PqyZXxlyRKfGt+5JHTUhGJm55jZ3KPV\nqVT/gJn9MryY0Lk665qZyZvDh/Ojnj3585YtFM6YwV+3bo13WM65KCTEdSjOQbD88K09ezJj9Gg6\nZ2Rw7qJFfGHRIjYfPBjv0JxzdRD1sOFwCvl+BCfjP8LM3otFUK55G9m6NdNHjeLe0lJ+uHo1b+/Y\nwQN9+3J+Xh6R0/E45xJLNMOGs4DHCS4urOl/tY/7dDGRnpLCzT16cE5uLhcvW8YFS5fyzJYtPNy/\nPwVZVf6Wcc4lgGi6vG4lmPfvIoKEchVwCfAvYCVwVqyDc25Qq1a8P3Ikv+jbl8kffsjgGTN4eP16\nyn0kmHMJJ5qE8jmCocMV879PM7Pfm9kngHkEC105F3OpElfn57NwzBjGtm7NFcuXc8rcuazw1SGd\nSyjRJJTuwKJwzqxDQOQVaI8DX4xlYM5V1qtFC94cPpzHBgxgblkZw2bO5MriYhbv2RPv0JxzRJdQ\nthFcKQ9QykevNcklmIrFuQYliYu7dGHR2LF8oWNHHtu4kcEzZnDK3Lm8sHUrh30WY+fiJppRXlOB\nkcDLBOuW3CGpNeF68ATnUpxrFN0yM/nDoEHc26cPv9u4kd9s2MB5ixbRLSODy7t25dKuXcnLyIh3\nmM41K9EcodxDMA0LwJ3APwnOqdwDrAKuiG1oztWuY0YGN/bowapx4/jbkCEUtmrFratXUzBlCl9e\nvJj/7NzpU7k410jqfIRiZjOBmeHj3cDnwhmIM81sVwPF51ydpEqcnZvL2bm5LNu7l4fWr+cPmzbx\nzJYtjMjO5squXcmPd5DONXF1OkKRNELSeZJOC5MI8N+pVjyZuIQyoGVLHujXj/XHH8/D/ftz2IxL\ni4v5PHDdihU+Osy5BlLbbMPtJP0TmAX8GfgHsELSkMYIzrn6yE5L4xtduzK/qIh3R4ygCPjl+vX0\nmz6dM+fP55Vt2/x6FudiqLYjlNuA44DbCS5cvJrgavgHGzYs52JHEie1a8cPgDXjxvGDHj2YW1bG\nWQsWUDh9Oo9u2MD+I76CtHP1VVtC+TRwh5ndYWavmdmDBFfKnxiO8HIuqXTNzOT2Xr1YM24cfxo0\niFapqVxWXEzPqVO5a80atvu0+c4ds9oSSk/g35XKKtZ/794QATnXGNJTUvhSXh4zR4/m7eHDGdm6\nNbeUlNB9yhSuXr6c1fv2xTtE55JObQklHThQqaxiLvFMnEtykjilfXteGzaM+UVFfK5jRx7asIG+\n06bxpcWLmbV7d7xDdC5p1GXY8GcqnYRPAQz4P0kjIiua2eOxDM65xjQ0O5snBg3irl69eGD9en67\nYQOTtmzh5Hbt+G5BARNycnz6fOeOoi4J5fs1lN9W6bkRzOnlXFLLz8ri3j59uKVHDx7ZsIEH1q3j\nzAULGNKqFdcXFPClTp3ISPG16ZyrrLaE0qtRonAuAbVNS+O73btzdX4+k7Zs4d7SUr66dCk3r1rF\nNfn5XNG1K9lpUa9R51yTddT/DWa2prECcS5RZaSkcGHnzlyQl8c/tm/n3tJSbli1invWruW6ggKu\n7NaNNp5YnPM15Z2rK0lM6NCBt0eMYMrIkRzXpg03l5TQc+pU7ly9mp2HD8c7ROfiKqqEIukiSa9L\nWixpVaXbyoYK0rlEM65tW14ZNowZo0bx8bZtuXX1anpMmcLtJSXs8GtZXDMVzZrytwI/BBYCc6k6\nnNi5ZqeoTRteGjqUObt3c8eaNfxwzRp+vm4d387P55r8fDqkp8c7ROcaTTQdvxcDD5jZdxoqGOeS\n1cjWrfnLkCHMKyvjzjVruHPNGn6xbh1XdevGtfn5dPS1WVwzEE2XVwfg7w0ViHNNwfDsbJ4bPJgF\nRUWc1aED96xdS6+pU7lh5Uo2HzxY+w6cS2LRJJR3+eiyv865GgzJzuaZwkIWjRnDZ3Nz+VlpKb2m\nTuW7K1eyzc+xuCaqtunrUypuwDXA1yRdKCk3cltEnXqRlCppjqSXw+c5kt6UtDy8bx9R9yZJKyQt\nk3RGfd/buYYwqFUrniosZMnYsZzXsSM/Ky2ldzgqrMxHhbkmprYkcBg4FN6WAUOA3wObI8orbrE4\nnr8aWBLx/EbgbTPrB7wdPkdSITARGAxMAB6SlBqD93euQfRv2ZInBw1iflERJ7drx62rV9Nn2jR+\ntW4dB8rL4x2eczFR20n5HxFMqdLgJOUTTJd/F3BtWHw2MD58/AQwGfheWD7JzA4AJZJWAGOBKY0R\nq3PHakh2Nn8bOpSpO3dyU0kJ316xgvvXreNHPXvy5bw8Un2uMJfEZAmyYp2k54GfAK2B683sLEkf\nmlm7cLuAHWbWTtKDwFQzeyrc9jvgNTN7vpr9XgZcBpCXlzd60qRJjfSJ4qusrIzs7Ox4h5FQEq1N\nDJgJPAosJ1gr4mLgBIL1IRpDorVJIvA2qerkk0+eZWZFtdVLiPkiJJ0FbDGzWZLGV1fHzExS1NnP\nzB4BHgEoKiqy8eOr3X2TM3nyZJrLZ62rRGyTk4HrzHhh61ZuKSnh1n37GNemDT/p1Yvx7dvX+vr6\nSsQ2iTdvk2NX20n5ayVlRbNDSaMkTYgyjhMIpsNfDUwCTpH0FLBZUpdwv12ALWH99UBBxOvzwzLn\nkk6KxOc7dWLRmDE82r8/pfv3c/K8eUyYN4/Zvh6LSyK1nZS/gOAcxd2SahwyLKm9pAskvUGwomOb\naIIws5vMLN/MehKcbP+nmX0FeIlgyWHC+xfDxy8BEyVlSuoF9AOmR/OeziWatJQULunaleXHHcd9\nffowY/duRs+axQVLllC6f3+8w3OuVrV1eY0iSCrXATdI2gUsALYSTL3SHugN9Amf/xkoNLPVMYrv\nbuBZSRcDa4AvAJjZIknPAosJRqJdaWZHYvSezsVVi9RUriso4JIuXbhn7VruLy3lha1bub6ggBsK\nCnzKfJewjnqEYoEnzWw4cDzwc2A3QRIZSXAC/X3g60BXM/tafZOJmU02s7PCx9vM7FQz62dmp5nZ\n9oh6d5lZHzMbYGav1ec9nUtEbdPS+HHv3iw77jjOzs3ljjVr6D99Or/fuJHyBBlM41ykOv+pY2bT\ngGkNGItzrho9srJ4prCQb3frxrUrV/L1Zcv41fr13N+nT6OcuHeurnw9FOeSxPFt2/KfkSN5ZtAg\nth06xMnz5nHOwoUs37s33qE5B3hCcS6pSGJiXh5Lx47lx7168daOHQyeMYNrV6zwdVhc3HlCcS4J\ntUhN5aYePVg+dixf7dyZB9ato284lcthn8rFxYknFOeSWOfMTB4ZMIA5RUWMzM7m2ytWMGrWLN7/\n8MN4h+aaIU8ozjUBw7KzeXP4cP46eDA7Dx/mpLlzuXDJEjYd8IVVXePxhOJcEyGJz3bsyJKxY/l+\n9+78ecsWBkyfzgPeDeYaSVQJRVIrSd+W9LykdyT1C8snShrYMCE656LRMjWVO3v3ZuGYMXysbVuu\n8W4w10jqnFAkFQDzgXsJpjo5ieDCRgjmuLs+5tE5545Zv5YteXXo0I90g13g3WCuAUVzhPIzgulV\n+gOj+egM2+8CJ8YwLudcDER2g93SowfPht1gvygt9W4wF3PRJJTTgR+Y2RqqLrq1HugWs6icczHV\nMjWVO3r1+m832HdWrmTUrFksjHdgrkmJJqFkEMzjVZ22BJM0OucSWOVusG8DVxYXs9PXt3cxEE1C\nmQ98roZtnwJm1T8c51xDq+gGWzRmDJ8DHt6wgcLp0/nb1q3xDs0luWgSyr3AxZIeJTghD1Ao6YcE\nK5feG+vgnHMNJzstjSuBqaNG0TE9nXMWLeLchQtZ7yft3TGqc0Ixs78A3wQ+D7wVFj8JXANcZWav\nxz4851xDG9OmDTNGj+ae3r15bft2Bk2fzkPr1/sU+S5qUV2HYmYPE5x8PwP4CkFXV364brtzLkml\np6RwQ/fuLBwzhuPatOHK5cv5+Jw5LNqzJ96huSRSp4QiKUPSXyWdZGZ7zOwtM/uTmf3DzHzRa+ea\niD4tWvDGsGE8OXAgxXv3MnLmTG4tKWH/EV8Q1dWuTgnFzA4Cp9W1vnMueUnigs6dWTJ2LBM7deLO\nNWsYMXMmU3bujHdoLsFFkyD+DYxrqECcc4mlY0YGTw4axD+GDWNfeTkfnzOH765cyT4/WnE1iCah\nXEcwyusqSfmSUiWlRN4aKkjnXPx8MieHBWPGcEmXLtxXWsqoWbOY6kcrrhrRJIEFQB/gAWANcBA4\nFHE7GPPonHMJoU1aGr8dMIA3hg1jz5EjnDBnDt9budLPrbiPSIui7o+oOuWKc64ZOT0nh4VjxnD9\nypX8tLSUv2/bxh8GDmRsmzbxDs0lgDonFDO7vQHjcM4liTZpaTwyYACf69iRS5Yt4/jZs/luQQG3\n9+xJVmpqvMNzceTnPZxzx+SM8Gjla507c09pKaNnzWLGrl3xDsvFUZ2PUCTdVksVM7M76hmPcy6J\ntE1L47GBAzmvY0cuLS7m+Nmzua1nT27u3p20FP97tbmJ5hzK7UfZVnFuxROKc83QhA4dWFBUxLdW\nrOAHq1fz2vbt/HHgQPq2bBnv0FwjimYur5TKNyAX+CqwEOjbQDE655JAu/R0/jhoEM8MGsTSvXsZ\nMXMmj23YgPmcYM1GvY5JzWy7mT0J/AH4dUwics4ltYl5ecwvKuK4Nm24tLiYcxYuZOtBv6qgOYhV\nJ+c8/jelvXOumSvIyuLN4cP5WZ8+vLZ9O0NnzOC1bdviHZZrYLFKKGcBvjqPc+6/UiSuLShgxujR\ndMzI4Mwyr31NAAAbQElEQVQFC7iyuJi9fjFkkxXNKK/HqynOAIYAQ4EfxCoo51zTMSw7mxmjRvH9\nkhLuX7eOt3fs4JnCQka2bh3v0FyMRXOEcgpwcqXbaGATwYqNd9UnEEntJD0vaamkJZKOl5Qj6U1J\ny8P79hH1b5K0QtIySWfU572dcw0rKzWVn/Xty1vDh1N25AjjZs/ml+vW+Qn7JiaaUV49zaxXpdsg\nM5tgZn+w+v/LeAB43cwGAsOBJcCNwNtm1g94O3yOpEJgIjAYmAA8JMkv0XUuwZ3avj1zi4r4ZE4O\nV69YwdkLF7Lt0KF4h+ViJCGuPJLUluCk/u8gWH/FzD4EzgaeCKs9AXw2fHw2MMnMDphZCbACGNu4\nUTvnjkVuRgYvDRnCL/r25fXt2xkxcybvf/hhvMNyMVDnhCLpbElfi3jeQ9IUSbvDrqrsesTRi+Ck\n/u8lzZH0mKRWQJ6ZbQzrbALywsfdgNKI168Ly5xzSUASV+fnM2XUKLJSUhg/dy53rF7NEe8CS2rR\nXCl/C/BcxPP7gXzgEeACgivpr69HHKOAb5nZNEkPEHZvVTAzkxT1vzZJlwGXAeTl5TF58uRjDDG5\nlJWVNZvPWlfeJlUlQps8APwcuG31av66ejU3E1wxHS+J0CZJy8zqdAO2AxPCxy2AfcDnw+eXACvr\nuq9q9t0ZWB3x/ETgFWAZ0CUs6wIsCx/fBNwUUf8fwPG1vc/o0aOtuXjnnXfiHULC8TapKlHapLy8\n3H6/YYO1fPddy/3Xv+yVDz6IWyyJ0iaJBJhpdfgtj+YcSlaYRAA+RnBU8Ub4fBnQNcpc9l9mtgko\nlTQgLDoVWAy8BFwUll0EvBg+fgmYKClTUi+gHzD9WN/fORdfkvhqly7MGj2arhkZfHrBAm5YuZLD\n5eXxDs1FIZour9XAx4F3CU6KzzKzinVAOwH1XRP0W8DTkjKAVcDXCM7xPCvpYoJVIr8AYGaLJD1L\nkHQOA1eamV8t5VySG9iqFdNGjeI7K1dyb2kpU3ft4s+FhXTJzIx3aK4OokkovwXuk3QOMAK4ImLb\n8QQ/7sfMzOYCRdVsOrWG+ndRz2tfnHOJJys1ld/0788JbdrwjeJiRs6cyaTCQsa3b1/7i11cRXMd\nygMEMwtPAb5uZo9GbG4N/D62oTnnmrOvdO7M9NGjaZeWxqnz5nHP2rWU+yiwhBbNEQpm9jTwdDXl\n34hZRM45FxrcqhUzRo/mkmXLuHHVKv6zcyd/GDiQ9unp8Q7NVSOa61D6Sxob8byFpJ9I+rukqxom\nPOdcc9c6LY1JhYX8sm9fXtu+ndGzZjF79+54h+WqEc0orweB8yKe3wVcRzC66+eSroxlYM45V0ES\n38rP570RIzhkxsdmz+ZRX7wr4USTUIYD/waQlAJcCHzPzEYDdxJePOiccw1lXNu2zB49mpPateOy\n4mIuKy7mgA8tThjRJJS2QMUKOSOB9sDz4fPJQO/YheWcc9XrmJHBa8OG8f3u3Xls40Y+MWcO6w8c\niHdYjugSymb+t278JwmujK+YTyub4HoQ55xrcKkSd/buzQuDB7Nwzx6KZs3i3zvreymcq69oEspL\nwE8k3Udw7iRyXq+hBBcjOudcozm3Y0emjR5NdmoqJ8+dy8Pr1/t5lTiKJqHcCLwMnEGQXH4cse3/\n+N80LM4512gGt2rF9FGjOK19e65YvtzPq8RRna9DMbM9wKU1bPtYzCJyzrkotU9P5+9Dh3JbSQk/\nXruWhXv28MLgwXT1KVsaVdQLbEnKlXSWpIsk5YRlWeHIL+eci4tUibt69+b5wYNZUFbG6Fmz+I+f\nV2lU0VzYKEn3Eixm9RLwONAz3Pwi8P2YR+ecc1H6XKXzKk9u2hTvkJqNaI4qbgKuAn4EHAcoYtvf\ngbNiGJdzzh2zweGsxR9v25aLli7lxpUrfTXIRhBNQrkE+JGZ/RiYXWnbCqBPzKJyzrl6yklP5/Vh\nw7i8a1fuKS3l3IUL2X3Yr25oSNEklG7A1Bq2HQRa1T8c55yLnfSUFH7Tvz8P9uvHK9u2ccKcOazZ\nvz/eYTVZ0SSU9cCQGrYNB0rqH45zzsXeld268dqwYazdv58xfhFkg4kmoTwH3CbphIgyk9Sf4ELH\nSTGNzDnnYuj0nBymheurnDJ3Lk/4yfqYiyah3A4sBd4DlodlzwELwud3xzQy55yLsQEtWzJ11ChO\nbNuWr4Yn633RrtiJZsXGfcB4glUb/wO8BcwgmGX4dDM72ADxOedcTOWkp/NaxMn6iYsXs+/IkXiH\n1STU6Up5SenAmcB8M/sj8McGjco55xpQekoKD/XrR5+sLL67ahXrDhzgxSFD6JiREe/QklqdjlDM\n7BDwLP+7kNE555KaJK7v3p3nCguZU1bG8bNnU7x3b7zDSmrRnENZBXRqqECccy4ezuvUiXeGD2fX\nkSMcP3s28+MdUBKLJqH8FPi+pI4NFYxzzsXDuLZtmTpqFB3T07ke+NPmzfEOKSnVebZh4BQgByiR\nNBXYCEQOjzAzuyiWwTnnXGPp3aIF/xk1ipP//W/OX7KEkv37ubl7dyTV/mIHRJdQPg4cArYSTLNS\neaoVH3vnnEtqOenp3As82akTt5SUsGb/fh7q14+0FJ9MvS6iWQ+lV0MG4pxziSAD+OOgQfTIyuLH\na9ey+eBBnikspGVqarxDS3jRTF+fKymrIYNxzrlEoHBtlV/368fft23jtHnz2HboULzDSnhHTSiS\nUiXdLmkHsBnYJekFSe0aJzznnIufb3brxnODBzN7924+7hNL1qq2I5TLgdsIpqu/j2BhrbOBnzdw\nXM45lxA+17Ejbw4fzqaDB4NhxWVl8Q4pYdWWUC4FHjWzU83se2Z2HnAl8BVJfkmpc65ZOLFdO94f\nMYIU4MQ5c5i8Y0e8Q0pItSWU3gQTQEb6M5AK9GiQiJxzLgENyc5myqhR5Gdmcsb8+Ty7ZUu8Q0o4\ntSWUbGBXpbLd4X3r2IfjnHOJqyAri/dHjmRsmzZMXLyY327YEO+QEkpdRnl1k9S74kZw1FKlPNx2\nzCR9R9IiSQslPSMpS1KOpDclLQ/v20fUv0nSCknLJJ1Rn/d2zrm6yklP541hw/h0hw5cXlzMT9as\nwXwKfKBu16E8X0P536opO6aB2pK6Ad8GCs1sn6RngYlAIfC2md0t6UbgRuB7kgrD7YOBrsBbkvqb\nmc9B7ZxrcC1SU/nL4MF8belSbi4pYduhQ9zbp0+zv6q+toTytUaJIpAGtJB0CGgJbABuIliDBeAJ\nYDLwPYKRZpPM7ADBVDArgLHAlEaM1znXjKWnpPDkoEG0T0/nZ+vWsePwYX7bv3+zvqr+qAnFzJ5o\njCDMbL2k+4C1wD7gDTN7Q1KemW0Mq20C8sLH3YCpEbtYF5ZVIekygkXAyMvLY/LkyQ3wCRJPWVlZ\ns/msdeVtUpW3SVXRtsm5BCeaH9+0ieWbNnELwdX2zVE0c3k1mPDcyNlAL+BD4DlJX4msY2YmKeqO\nSjN7BHgEoKioyMaPH1//gJPA5MmTaS6fta68TaryNqnqWNrkZGDUunVcs2IF97Rrx9+GDKF1WkL8\nvDaqRDk2Ow0oMbOt4WJefwE+BmyW1AUgvK8Yp7ceKIh4fX5Y5pxzcXF1fj5PDBzIux9+yKnNdKqW\nREkoa4FxkloqOKt1KrCE4Mr8iinxLwJeDB+/BEyUlCmpF9APmN7IMTvn3Edc2LkzLwwZwryyMk6e\nO5ctBw/GO6RGlRAJxcymEYwmmw0sIIjrEeBu4HRJywmOYu4O6y8iWJJ4MfA6cKWP8HLOJYKzc3N5\neehQVuzbx0lz5rD+wIF4h9RoEiKhAJjZD8xsoJkNMbMLzOyAmW0Lp33pZ2anmdn2iPp3mVkfMxtg\nZq/FM3bnnIt0ek4O/xg2jA0HD3LSnDms3rcv3iE1ioRJKM4515Sc2K4dbw0fzvbDhzlx7lyK9+6N\nd0gNzhOKc841kLFt2vDO8OHsLy/npDlzWNjEZyr2hOKccw1oROvWvDdiBCkS4+fOZfbu3bW/KEl5\nQnHOuQY2qFUr3hsxglapqZwydy7TdlWec7dp8ITinHONoG/Llrw/ciQd0tP55Lx5TN25M94hxZwn\nFOecayTds7KYPGIEHdPT+eT8+fyniSUVTyjOOdeICrKyeHfkSDpnZHDG/Pn868MP4x1SzHhCcc65\nRtYtM5PJI0bQNSODCfPn834TSSqeUJxzLg66hkmlICuLT82fz7tNIKl4QnHOuTjpkpnJO8OH0yMr\nizPnz+edHTviHVK9eEJxzrk46pyZyTsjRtArK4tPL1iQ1EcqnlCccy7OOmVk8M8RI+iZlcWn58/n\n30k6+ssTinPOJYBOGRm8PXw43TIz+dT8+Ul5nYonFOecSxBdMjP554gRdEpP54z585mZZFfUe0Jx\nzrkE0i08p9IhPZ3T589nThLN/eUJxTnnEkxBVhb/HD6cNqmpnD5vHguSZJZiTyjOOZeAerZowT9H\njCArJYVT581j8Z498Q6pVp5QnHMuQfVp0YJ3RowgTeK0efNYmeArP3pCcc65BNavZUveGj6cg+Xl\nnDp3LqX798c7pBp5QnHOuQRX2KoVbwwfzo7Dhzlt3jw2HzwY75Cq5QnFOeeSwKjWrXl12DDWHTjA\n6fPmsf3QoXiHVIUnFOecSxIntG3Li0OGsGzvXibMn8+uw4fjHdJHeEJxzrkkclpODs8PHsycsjI+\ns2ABe48ciXdI/+UJxTnnksxncnP548CBvL9zJ+cuXMjB8vJ4hwR4QnHOuaQ0MS+PRwcM4B87dnDh\nkiUcMYt3SKTFOwDnnHPH5uIuXdh+6BA3rFpFzvLl/LpfPyTFLR5PKM45l8S+2707Hxw6xE9LS+mQ\nns4dvXpxoLycF7Zu5d61a1m6bx8HysvJTElhYIsW3NC9O+d27EhmSuw7qDyhOOdckru7d2+2Hz7M\nnWvWULxnD6/v2IEBuyNO2O8vL2funj1cVlzM5cXF3NenD5d27RrTOPwcinPOJTlJPNy/PwNatODZ\nDz5g15EjH0kmkcqOHGHXkSNcs2IFt5aUxDQOTyjOOdcEPL5xI6UHDtS5/t7ycu4vLeXRDRtiFoMn\nFOecS3IHysu5fuVK9kY5fHhv+LpYDTtu1IQi6XFJWyQtjCjLkfSmpOXhffuIbTdJWiFpmaQzIspH\nS1oQbvul4jmswTnn4uwvW7dyrIOGy8PXx0JjH6H8AZhQqexG4G0z6we8HT5HUiEwERgcvuYhSanh\na34DXAr0C2+V9+mcc83GT9eurfGcSW3KjhzhntLSmMTRqAnFzN4DtlcqPht4Inz8BPDZiPJJZnbA\nzEqAFcBYSV2ANmY21cwMeDLiNc451+wsrec6Kcv27o1JHIkwbDjPzDaGjzcBeeHjbsDUiHrrwrJD\n4ePK5dWSdBlwGUBeXh6TJ0+OTdQJrqysrNl81rryNqnK26SqZGyTup+Kr97+8vKYfOZESCj/ZWYm\nKabzB5jZI8AjAEVFRTZ+/PhY7j5hTZ48mebyWevK26Qqb5OqkrFNMt97j/31OLGelZLC+JNOqncc\niTDKa3PYjUV4vyUsXw8URNTLD8vWh48rlzvnXLM0sEWLer1+QMuWMYkjERLKS8BF4eOLgBcjyidK\nypTUi+Dk+/Swe2yXpHHh6K4LI17jnHPNzg3du9M6NbX2itVonZrK9woKaq9YB409bPgZYAowQNI6\nSRcDdwOnS1oOnBY+x8wWAc8Ci4HXgSvNrGIYwzeBxwhO1K8EXmvMz+Gcc4nk3I4dOdZrJxS+PhYa\n9RyKmX2phk2n1lD/LuCuaspnAkNiGJpzziWtzJQU7uvTh2tWrIjq4saW4esyYjRRZCJ0eTnnnKun\nS7t25dqCAlrWMTm0TEnh2oKCmE4QmVCjvJxzzh27O3r1ontmJtevXEk5wUWLlWWnppICDTLbsCcU\n55xrQi7t2pULO3fmL1u38tPSUpbt3cv+8nKyUlIY0LIl3yso4NyOHWPWzRXJE4pzzjUxmSkpfCkv\njy/l5dVeOYZkCbAOcWORtBVYE+84Gkku8EG8g0gw3iZVeZtU5W1SVQ8zq3UoWLNKKM2JpJlmVhTv\nOBKJt0lV3iZVeZscOx/l5ZxzLiY8oTjnnIsJTyhN1yPxDiABeZtU5W1SlbfJMfJzKM4552LCj1Cc\nc87FhCcU55xzMeEJJQlJKpD0jqTFkhZJujosz5H0pqTl4X37iNfcJGmFpGWSzohf9A3jKG1yr6Sl\nkuZL+qukdhGvaZZtErH9OkkmKTeirNm2iaRvhf9WFkn6aUR5k26TmDIzvyXZDegCjAoftwaKgULg\np8CNYfmNwD3h40JgHpAJ9CKY8j813p+jkdrkk0BaWH6PtwmF4fMC4B8EF/rmNvc2AU4G3gIyw22d\nmkubxPLmRyhJyMw2mtns8PFuYAnQDTgbeCKs9gTw2fDx2cAkMztgZiUE68iMbdyoG1ZNbWJmb5jZ\n4bDaVP632mezbZNw88+BG4DIUTnNuU2uAO42swPhtoqVY5t8m8SSJ5QkJ6knMBKYBuRZsKIlwCag\nYiKfbkBpxMvW8b8flianUptE+jr/W4yt2baJpLOB9WY2r1K1ZtsmQH/gREnTJL0raUxYrVm1SX35\n5JBJTFI28AJwjZntClZEDpiZSWp2Y8Irt0lE+feBw8DT8YotXiLbhKANbiboCmy2qvm/kwbkAOOA\nMcCzknrHM8Zk5EcoSUpSOsF/iKfN7C9h8WZJXcLtXYCKw/b1BH3mFfLDsialhjZB0leBs4DzLewY\np/m2SR+CcwHzJK0m+NyzJXWm+bYJBEcef7HAdKCcYJLIZtEmseIJJQkpOBT5HbDEzO6P2PQScFH4\n+CLgxYjyiZIyJfUC+gHTGyvexlBTm0iaQHCu4P/MbG/ES5plm5jZAjPrZGY9zawnwQ/pKDPbRDNt\nk9DfCE7MI6k/kEEw43CTb5NY8i6v5HQCcAGwQNLcsOxm4G6CQ/WLCUbvfAHAzBZJehZYTNDlcaWZ\nVV3KLbnV1Ca/JBih82bYJTjVzC5vzm1iZq9WV7k5twnwOPC4pIXAQeCi8Gi2ObRJzPjUK84552LC\nu7ycc87FhCcU55xzMeEJxTnnXEx4QnHOORcTnlCcc87FhCcUFzfhTLe13VbHO854knSapNviHUdN\nJH1MUpmkThFlUyW9VU3de8Pv9Jbw+URJ6yS1aMyYXcPxhOLi6fhKt00EM+BGlp0Tt+gSw2lAwiYU\n4D7g4YjJFKtQ4JfA9cD1ZnZnuOlZYDfBlDCuCfALG13cmNnUyOeSDgAfVC5vSiSlElz/dbjWyg0X\ng4B0MztYz/2cQJD0L6zlvX4LXAJcZWa/rthmZuWSHgVukHSfmR2qTzwu/vwIxSWNsPtnctjFUibp\nFUmDKtWZKuktSZ9RsKjWPkkzJY2WlB52u2yWtE3So5HdLZIGhl0yl0j6laQPJO2R9KKkgkrvI0lX\nSlogab+kLZJ+K6ltRJ2scH+3SbpV0hqCq7D7SWol6ZcKFnraI2mDpL9J6hfx+ruB7wGpEV2A+8Nt\nE8Ln4yrFdXlY3jmibJOkx8JtxcAh4NRwW2tJP5O0RtJBSSsl3RAmgtpcAkw3sxU1fF8pwO+Bi4HL\nIpNJhEkEs2J/pg7v5xKcH6G4pCDpXOA54K/Al4FU4CbgPUnDIqbth2BRpDuAO4H9BN0yLxIsoHSQ\n4C/qYcBPgI1U7VL6ATAjrNc1rPeapOER0278HPhmeP82wQSCdwGFkj5hZuUR+/sGsIyga2c/waSd\nLYEWYZwbCCYivAqYImmAmW0Dfh2+/5eBj4f7itxvND5FMIvurcA2YIWkjLBNeoVxLCGYmuROoC3w\n/Vr2eQbwpxq2pRHM7Px5gmlMnqqukpltkLQSmAD8pbo6LonEe4Uvv/mt4gasBp6qpjyFYE2KVyuV\n5wAfEiyMVFE2leBHuyCi7AsEC0m9XOn1rxJMEljxfGBYbw7htERh+alh+fnh8/4EP+w3VNpfRb0J\n4fOs8PlqIKOWz55KsILgfuCKiPK7gcPV1J8Q7ntcpfLLw/LOEWWbCM5V5Faqe2n4OY6rVH4HsA9o\nd5R4e4Tvc0E126aG24xg7rDavvfngPnx/vfnt/rfvMvLJYPBBNOGPyUpreIG7CI4kjipUv1FZha5\nKNLS8P4fleot5aNTk1d4zsJfOgAze5tg5tnjw6IzAAFPV4rnPeBANfG8atWcr5B0vqQZknYSTDy4\ni2AiywHVxFRf75vZB5XKJhAsgTur0ud4gyAZHm1lwq7h/dYati8imKD0O5IG1xLb1oj9uSTmCcUl\ng4ohqU8T9P9H3k4DOlSqv6PS84NHKc+q5v0211BWsVJfRTzrKsVykCAhVI5nY6XnSPo88BQwF5gI\nHEfQJbWzhpjqq0oMBJ9jAFXb9L1we+XPEakixgM1bN8EnBJuf1vS0ZLkPoLuP5fk/ByKSwbbwvvr\n+N+PXaT9MX6/vBrKJleKZzywp5q6lf9qr25K74nAQjO7tKJAUkuCcxd1UfGZMyqV15QEqothG8G5\nna/U8JpVR3n/ijZoX1MFM1sl6RTgXYKk8gkzW1lN1RyCI0CX5DyhuGSwgODE9SD76KJIDeXzkn5S\n0e0l6VSCk+ZTwu1vEPxA55vZsS4p3JKgmyvSV6upd4BglFe6fXRY7ZrwfggfTbJnRhHD6wTdXjtq\n+KE/mpUE8R91mVwzK5Z0GvAO8E9JJ5nZmkrVehEkNpfkPKG4hGdmRyRdBTwX/hX/AsFfyJ0JRiUV\nm9mDMXzLXOAFSY8BXQhGeS0iGOKKmS2W9AvgEUlDgPcJfvi7E6zV/isz+08t7/E68AtJ9xAkqOMI\nTqiXVaq3OLz/roKrzw+b2WwzK5E0DbgtPAeznSAhdaPufk+wsuc7kn4GLCTosusL/B9whtWwmJSZ\n7ZE0i6OfZ6mou0jS6cA/CY5UTjKzDfDf63JGA/dEEbdLUH4OxSUFM/srwRKtOQRLuP6DYARULrFf\nkvWHBOuGPwn8imDU0qcif1zN7FrgWwTncJ4nWEL2eoLurpI6vMeDwE8Jhia/FO7n01TtQnsBeBS4\nNowjMlF9EZgNPESw4uBS4N66fkgzO0AwMu1J4EqCUW9/BM4nOOqpbYjyn4HTJWXW4b3mEQxmyCU4\nUqnoVhwPtAr35ZKcr9joXEjSQIJrMS6wGq6bcP8jKYdgOPdFZvb8Me7j9wRdh6fHNDgXF36E4pw7\nJma2Hbif4Gr+qIWzD3wJuCWWcbn48XMozrn6uAc4LKmTHWWCyBr0AL5lZtMaIC4XB97l5ZxzLia8\ny8s551xMeEJxzjkXE55QnHPOxYQnFOecczHhCcU551xM/D/ZrVbg84FW9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1168d7290>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'CO2': Field([ 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,\n",
       "         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,\n",
       "         0.00038,  0.00038,  0.00038,  0.00038,  0.00038,  0.00038]),\n",
       " 'H2O': Field([  5.00000000e-06,   5.00000000e-06,   5.00000000e-06,\n",
       "          5.00000000e-06,   5.00000000e-06,   5.00000000e-06,\n",
       "          6.38590233e-06,   9.08848690e-06,   1.33273826e-05,\n",
       "          2.34389689e-05,   3.84220914e-05,   5.95564299e-05,\n",
       "          8.82144990e-05,   1.25843839e-04,   1.73951159e-04,\n",
       "          2.34088411e-04,   3.07840683e-04,   3.96815735e-04,\n",
       "          5.02635028e-04,   6.26926041e-04,   7.71315753e-04,\n",
       "          9.37425100e-04,   1.12686431e-03,   1.34122899e-03,\n",
       "          1.58209684e-03,   1.85102493e-03,   2.14954752e-03,\n",
       "          2.47917415e-03,   2.84138824e-03,   3.23764591e-03]),\n",
       " 'O3': Field([ 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",
       "         0.,  0.,  0.,  0.])}"
      ]
     },
     "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": "markdown",
   "metadata": {},
   "source": [
    "We need to provide some ozone data to the model in order to simulate a stratosphere. As we did with the original column model, we will use the ozone climatology data provided with the CESM model.\n",
    "\n",
    "See here for more information, including some plots of the ozone data:\n",
    "<http://www.atmos.albany.edu/facstaff/brose/classes/ENV480_Spring2014/styled-5/code-3/index.html>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import netCDF4 as nc\n",
    "\n",
    "datapath = \"http://ramadda.atmos.albany.edu:8080/repository/opendap/latest/Top/Users/BrianRose/CESM_runs/\"\n",
    "endstr = \"/entry.das\"\n",
    "ozone = nc.Dataset( datapath + 'som_input/ozone_1.9x2.5_L26_2000clim_c091112.nc' + endstr )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Dimensions of the ozone file\n",
    "lat = ozone.variables['lat'][:]\n",
    "lon = ozone.variables['lon'][:]\n",
    "lev = ozone.variables['lev'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Taking annual, zonal, and global averages of the ozone data\n",
    "O3_zon = np.mean( ozone.variables['O3'],axis=(0,3) )\n",
    "O3_global = np.sum( O3_zon * np.cos(np.deg2rad(lat)), axis=1 ) / sum( np.cos(np.deg2rad(lat) ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGBlJREFUeJzt3X2sXPV95/H39859xBg/YOMYG7DbWCEOKhu4YUnTdLtL\nd0PTbB1ltchpk7WqqKxU2qZVpRYi7Ub9AykrVd0225IuG9KlSjbIokiwVbYJ6ySqoioQG7IKNqG4\nGD9hY4cH29i+TzPf/WOOr8fXnnvgjn3n+sz7JY3mzO+cM+d3B/z7nPP7nYfITCRJvamv2xWQJHWP\nISBJPcwQkKQeZghIUg8zBCSphxkCktTDDAFJ6mGGgCT1MENAknpYf7crUGbFihW5bt26bldDki4r\nO3bs+ElmrixbbsGHwLp169i+fXu3qyFJl5WI2Pt2lrM7SJJ6mCEgST3MEJCkHmYISFIPm/cQiIg7\nI+KFiNgdEffO9/YlSWfNawhERA34C+CXgI3AJyNi43zWQZJ01nwfCdwG7M7MlzJzAngE2DTPdZAk\nFeb7OoE1wP6WzweAfz7PdZCkt+3QsdP8r6f2MVjrY7C/j4HifbC/j6H+vguWD9aKef0z5tWar76+\n6PafNW1BXiwWEXcDdwNcf/31Xa6NpF526NgYf/6d3VzMx7EP1KI0PDauvoo/2nTTxdtoG/MdAgeB\n61o+ry3KzpGZDwIPAoyOjl7En15SlWQmk/Vkot5gcqrBZL3RnK5nc3qq0TIvW+Y35zU/J+OTdcan\nGoxN1otXMV2U/dy7V3ByfIrjY1McPz3J8bFJxiYbc653sy51Tk7UGezvY7i/j5HBGvWB2nTYTNTn\np+mb7xD4AbAhItbTbPw3A786z3WQVCIzqTdaGtiZDedUni2bLs8Zy81oiKdfOaNxbpyzndbG+ez8\n5jLj53xull1METDcX2N4oI/hgRrDAzWG+s9M97F22QjDK688d/5AH8P9NUYGawz3ny0fHuhjaKB2\n3vcNF8uf+e5udw3Nawhk5lRE/BbwTaAGfCUzd85nHaSFqtFIxqbO7oWebtkrHZ/+fHbedINYvI/X\nG0zOaJybDW79nAa5tSFu3TOeLBrdiaLsYnZ/nNHfFwzU+prdIUVXyEBrt0gtpstGBs/0t58tO7PM\n+evGOZ8Hp5cPBmb025/pihmo9Z2dV+tjaKDZFROxcPrr58O8jwlk5jeAb8z3dqV36kxXw9hUnbGJ\nogGeqnN6oj7dVXB6os741NnG+vTM7oTidXpmF8NEvWjwi++bajbQnRg8p7GLosE809id/bxoqP9s\nA1mb2UA2lx2cpXFuNpznNsTT27rAuq3b6fZer863IAeGpXZa95ZPtzSy5za6M+e1NsB1Tk80G/ML\n7l3P2AtvzHFv+Mye5UhLF8DIQI2hgRpLRgZ411VDzfKiG2GopYtg5Jyug7NdCed9V39tusGv9UXP\n7cHq4jAENC/qjeTE2CTHTp//On56qmW6+X5ifGp6b7l1z3uue8t9wdlGtb+P4cGzfbUjgzWWjgy0\naXT7in7fls9n+n8H+hianq5ND+4N9deoucery4QhoLet3sjpRvrCjXn7eW+NT83axzxQC5aMDHDV\nyABLile7veXzB+Bm2VsupgdrvdfXK70dhkCPmao3OD42NXsjfurC806MT8363YO1vqIR72fJyADX\nLB5iwzVXTjfqrQ38kpEBllxxdnpkoGYjLXWBIVABJ8YmOXRsjFfePM0rb45x6NhpXj0+1tKIT003\n8G+VNORD/X3nNNSrlwxz47sWn9+Az2jElxTdKZIuL4bAAjc2WefQsTEOvXmaV6bfzzb2h94cO28P\nvS9gxZVDLC0a6TVLh3nv6sXnN+IzXlfZkEs9xxDoosl6g1ePj03vxbc29mc+v35y4rz1rl40yOql\nw9xw9SI++FNXs3rpCNcuHeHaJcOsXjrCNYuHGKj5qAhJ5QyBS2iy3uCFwyfY//qp6b34Q8fGeKXY\ngz9yYuy8UxAXD/dz7ZIRVi8d5mfWLp1u2K9dOsy1S0Z415Jh99YlXTSGwEU0Nlnn/+1/k6f3vM7T\nL7/Ojr1vcGqiPj1/eKBvuoH/uQ0rphv41UuGWbN0hNVLR7hyyP8kkuaPLU4HTk1MsWPvGzy953We\n2vM6P9z/5vR57De+azH//ta1jK5bzvoVi1izdISlVwx4BoykBcUQeAfeGp/iqZdem270nzt4jKlG\nUusLbrr2KrZ88AZuW381H1i3jKVXDHa7upJUyhB4G14/OcFD33uJv/6HvZwYn2Kw1sfN1y3hP/6L\nn+K29Vdz6w3L7MaRdFmy5ZrFkRNj/I+/f4mvfn8fY1N1PnrTan7t9uu55fplDs5KqgRD4AIOHTvN\nX373n/j6D/YzVW+w6Z+t4Z5/+dO8+5rF3a6aJF1UhsAM33nhCPd87Rkmphp84pY1/OYvvJt1KxZ1\nu1qSdEkYAi22bt/PfY/9iPesWsx///StXLf8im5XSZIuKUOA5sNDvrhtN//1//4jH96wggd+7RYW\nDw90u1qSdMn1fAhkJv/p8ef46vf38Ylb1vBf/t3PeMsFST2j50PgmzsP89Xv7+M3Pryez330vV7M\nJamn9PQu78nxKf7of+/ixnct5g/vvNEAkNRzevpI4IvbXuTQsTH+/FffT79dQJJ6UM+2fP/46gke\n+t4e7hpdy603LO92dSSpK3o2BP7bt3czMlDjD++8sdtVkaSu6ckQeOPkBN987jCfuGUNV1851O3q\nSFLX9GQIPPbsQSbqDTbfdn23qyJJXdVzIZCZPPL0Pm6+binvXX1Vt6sjSV3VcyHwzL43ePHIW3zy\nA9d1uyqS1HU9FwKP//AVhgf6+NjN13a7KpLUdT0VApnJk7te5cMbVvoQGEmix0LguYPHOXRsjH+z\ncVW3qyJJC0JPhcC3dh2mL+CO9xoCkgQdhEBEXBcR34mIXRGxMyI+W5Qvj4gnI+LF4n1Zyzr3RcTu\niHghIj5yMf6Ad+LJXa/ygXXLWb7Ih8BLEnR2JDAF/H5mbgRuB+6JiI3AvcC2zNwAbCs+U8zbDLwP\nuBN4ICLm7UG9+147xY8Pn+Bf2xUkSdPmHAKZeSgznymmTwDPA2uATcDDxWIPAx8vpjcBj2TmeGbu\nAXYDt811++/UP/zTTwD4hfdcM1+blKQF76KMCUTEOuD9wFPAqsw8VMw6DJzZ9V4D7G9Z7UBRNi9+\n8PIbLF80yE+v9HnBknRGxyEQEVcCfwP8bmYeb52XmQnkHL7z7ojYHhHbjx492mkVAdix93VuvWGZ\nzwyQpBYdhUBEDNAMgK9l5mNF8asRsbqYvxo4UpQfBFov011blJ0nMx/MzNHMHF25cmUnVQTg6Ilx\nXn7tFKM3LCtfWJJ6SCdnBwXwEPB8Zv5Jy6wngC3F9Bbg8ZbyzRExFBHrgQ3A03Pd/juxY+/rAIyu\n87kBktSqk8tmPwR8GvhRRPywKPsc8AVga0R8BtgL3AWQmTsjYiuwi+aZRfdkZr2D7b9t219+g8H+\nPm5a4w3jJKnVnEMgM78HtOtgv6PNOvcD9891m3O1Y98b3Lx2CUP983ZGqiRdFip/xXBmsvvIW942\nWpIuoPIh8MapSU6MTbHuak8NlaSZKh8CL792EoB1K67ock0kaeGpfAjsLULgBo8EJOk8lQ+Bl39y\nir6AtctGul0VSVpwqh8Cr53k2qUjnhkkSRfQAyFwykFhSWqj8iGw77WTXH+1g8KSdCGVDoF6I3nj\n1CTXLB7qdlUkaUGqdAgcPz0JwJKRgS7XRJIWpkqHwJuGgCTNqtIhcMwQkKRZGQKS1MMMAUnqYYaA\nJPWwSofAmbODrjIEJOmCKh0Cb56aYHigj+EBbxkhSRdS6RB4a7zOlUOdPEFTkqqt0iEwMdXwxnGS\nNItqh0C9wWB/pf9ESepIpVvIiak6g7VK/4mS1JFKt5ATUx4JSNJsKt1C2h0kSbOrdAs5MdWwO0iS\nZlHpFtLuIEmaXaVbyHFDQJJmVekW0jEBSZpdpVvIiakGQ44JSFJblW4hJ+sN+mvR7WpI0oJV6RCo\nN6DWV+k/UZI6UukWspGJvUGS1F7HTWRE1CLi2Yj42+Lz8oh4MiJeLN6XtSx7X0TsjogXIuIjnW67\nTL2R1MLuIElq52LsJ38WeL7l873AtszcAGwrPhMRG4HNwPuAO4EHIuKS3uKz0Uj6+gwBSWqnoxCI\niLXALwNfbineBDxcTD8MfLyl/JHMHM/MPcBu4LZOtl+mnh4JSNJsOj0S+FPgD4BGS9mqzDxUTB8G\nVhXTa4D9LcsdKMrOExF3R8T2iNh+9OjROVeu3khqHglIUltzDoGI+BhwJDN3tFsmMxPId/rdmflg\nZo5m5ujKlSvnWkUaaXeQJM2mk2cvfgj4lYj4KDAMXBURXwVejYjVmXkoIlYDR4rlDwLXtay/tii7\nZBwYlqTZzflIIDPvy8y1mbmO5oDvtzPzU8ATwJZisS3A48X0E8DmiBiKiPXABuDpOde8vH40Eo8E\nJGkWl+Ip7F8AtkbEZ4C9wF0AmbkzIrYCu4Ap4J7MrF+C7QPQKDqhPBKQpPYuSghk5neB7xbTrwF3\ntFnufuD+i7HNMvUiBbxYTJLaq2wT2chmCNgdJEntVT8E7A6SpLYqGwLT3UGGgCS1VdkQODMwbHeQ\nJLVX2RDI6e6gLldEkhawyobAme4gxwQkqb3KhoDdQZJUrrIhYHeQJJWrbAjUPUVUkkpVNgS8bYQk\nlatuCBQpYAZIUnvVDQG7gySpVIVDoPnuk8Ukqb0Kh4DdQZJUproh4MViklSquiFgd5AklapwCBTd\nQV2uhyQtZJUPAW8bIUntVTYE8sy9gxwTkKS2KhsCDe8dJEmlKhwCzXePBCSpvQqHQJECZoAktVXZ\nEHBMQJLKVTgEHBOQpDKVDQHHBCSpXIVDwHsHSVKZyoeARwKS1F5lQ8CTgySpXOVDwNtGSFJ7lQ0B\nrxiWpHKVD4FwTECS2uooBCJiaUQ8GhE/jojnI+KDEbE8Ip6MiBeL92Uty98XEbsj4oWI+Ejn1W/P\ni8UkqVynRwJ/BvxdZt4I3Aw8D9wLbMvMDcC24jMRsRHYDLwPuBN4ICJqHW6/LbuDJKncnEMgIpYA\nPw88BJCZE5n5JrAJeLhY7GHg48X0JuCRzBzPzD3AbuC2uW6/TGP67CBTQJLa6eRIYD1wFPiriHg2\nIr4cEYuAVZl5qFjmMLCqmF4D7G9Z/0BRdp6IuDsitkfE9qNHj86pcunFYpJUqpMQ6AduAb6Ume8H\nTlJ0/ZyRzZY43+kXZ+aDmTmamaMrV66cU+W8bYQkleskBA4ABzLzqeLzozRD4dWIWA1QvB8p5h8E\nrmtZf21RdklM30Cusuc/SVLn5txEZuZhYH9EvKcougPYBTwBbCnKtgCPF9NPAJsjYigi1gMbgKfn\nuv0yHglIUrn+Dtf/beBrETEIvAT8Os1g2RoRnwH2AncBZObOiNhKMyimgHsys97h9tuqe3aQJJXq\nKAQy84fA6AVm3dFm+fuB+zvZ5tuVXiwmSaUq22N+5jqBmiEgSW1VNwQazXfHBCSpveqGgNcJSFKp\nyoaAt5KWpHKVDQHvHSRJ5SocAs13xwQkqb0Kh4BjApJUprIhkD5oXpJKVTYE7A6SpHIVDgEHhiWp\nTIVDoPnubSMkqb3KhkB6JCBJpSobAg0HhiWpVIVDoPluBkhSexUOAY8EJKlMZUMgPUVUkkpVNgQa\nDQeGJalMdUPAIwFJKlXhEPDeQZJUprIhkJlEeLGYJM2msiFQz7QrSJJKVDYEGulD5iWpTIVDIB0P\nkKQSlQ2BTM8MkqQylQ2BRiO9RkCSSlQ3BDwSkKRSFQ4BxwQkqUxlQyAz6bM/SJJmVdkQsDtIkspV\nOAQcGJakMhUOAW8ZIUllOgqBiPi9iNgZEc9FxNcjYjgilkfEkxHxYvG+rGX5+yJid0S8EBEf6bz6\n7XmKqCSVm3MIRMQa4HeA0cy8CagBm4F7gW2ZuQHYVnwmIjYW898H3Ak8EBG1zqrfXsN7B0lSqU67\ng/qBkYjoB64AXgE2AQ8X8x8GPl5MbwIeyczxzNwD7AZu63D7bSUODEtSmTmHQGYeBP4Y2AccAo5l\n5reAVZl5qFjsMLCqmF4D7G/5igNF2Xki4u6I2B4R248ePTqn+p15noAkqb1OuoOW0dy7Xw9cCyyK\niE+1LpOZSXOn/B3JzAczczQzR1euXDm3Cib0VXbYW5Iujk6ayV8E9mTm0cycBB4DfhZ4NSJWAxTv\nR4rlDwLXtay/tii7JBqZBHYHSdJsOgmBfcDtEXFFNM/FvAN4HngC2FIsswV4vJh+AtgcEUMRsR7Y\nADzdwfZn1RwTuFTfLknV0D/XFTPzqYh4FHgGmAKeBR4ErgS2RsRngL3AXcXyOyNiK7CrWP6ezKx3\nWP+2vE5AksrNOQQAMvPzwOdnFI/TPCq40PL3A/d3ss23K72BnCSVquzQaSaOCEhSieqGAF4sJkll\nKhsCjQZ2B0lSicqGgEcCklSusiHQ8IJhSSpV2RBIHyojSaUqHAKeIipJZaobAjgwLEllKhsCPk9A\nkspVNgS8WEySylU3BMD+IEkqUd0QyPRIQJJKVDYEwAMBSSpT2RBwTECSylU3BEifJyBJJaobAh4J\nSFKpaoeAKSBJs6puCOCD5iWpTHVDILE/SJJKVDcEMAMkqUxlQwDHBCSpVGVDwDEBSSpX3RDwSECS\nSlU3BPDJYpJUprIh0PDJYpJUqrIhkD5oXpJKVTcEwHsHSVKJyoYAPk9AkkpVNgR80LwklatsCDQy\nqZkCkjSryoZAveGYgCSVKQ2BiPhKRByJiOdaypZHxJMR8WLxvqxl3n0RsTsiXoiIj7SU3xoRPyrm\nfTEucQudmdQqG3GSdHG8nWbyfwJ3zii7F9iWmRuAbcVnImIjsBl4X7HOAxFRK9b5EvAbwIbiNfM7\nL6p6I71YTJJKlIZAZv498PqM4k3Aw8X0w8DHW8ofyczxzNwD7AZui4jVwFWZ+f3MTOCvW9a5JD68\nYSW33rCsfEFJ6mH9c1xvVWYeKqYPA6uK6TXA91uWO1CUTRbTM8svKCLuBu4GuP766+dUwf/8bzfO\naT1J6iUd95oXe/YX9frczHwwM0czc3TlypUX86slSS3mGgKvFl08FO9HivKDwHUty60tyg4W0zPL\nJUldNNcQeALYUkxvAR5vKd8cEUMRsZ7mAPDTRdfR8Yi4vTgr6D+0rCNJ6pLSMYGI+DrwC8CKiDgA\nfB74ArA1Ij4D7AXuAsjMnRGxFdgFTAH3ZGa9+KrfpHmm0Qjwf4qXJKmLIhf47TZHR0dz+/bt3a6G\nJF1WImJHZo6WLeflVJLUwwwBSephhoAk9bAFPyYQEUdpDj6/UyuAn1zk6lSFv82F+bu052/T3kL9\nbW7IzNILrRZ8CMxVRGx/O4Mivcjf5sL8Xdrzt2nvcv9t7A6SpB5mCEhSD6tyCDzY7QosYP42F+bv\n0p6/TXuX9W9T2TEBSVK5Kh8JSJJKVC4EIuLO4tGWuyPi3m7XZ6GIiOsi4jsRsSsidkbEZ7tdp4Um\nImoR8WxE/G2367KQRMTSiHg0In4cEc9HxAe7XaeFICJ+r/i39FxEfD0ihrtdp7moVAgUj7L8C+CX\ngI3AJ4tHXqp5Q7/fz8yNwO3APf425/ks8Hy3K7EA/Rnwd5l5I3Az/kZExBrgd4DRzLwJqNF8tO5l\np1IhANwG7M7MlzJzAniE5iMve15mHsrMZ4rpEzT/Ibd9uluviYi1wC8DX+52XRaSiFgC/DzwEEBm\nTmTmm92t1YLRD4xERD9wBfBKl+szJ1ULgTXA/pbPsz7GsldFxDrg/cBT3a3JgvKnwB8AjW5XZIFZ\nDxwF/qroKvtyRCzqdqW6LTMPAn8M7AMOAccy81vdrdXcVC0EVCIirgT+BvjdzDze7fosBBHxMeBI\nZu7odl0WoH7gFuBLmfl+4CTQ82NtEbGMZi/DeuBaYFFEfKq7tZqbqoVAu8dbCoiIAZoB8LXMfKzb\n9VlAPgT8SkS8TLML8V9FxFe7W6UF4wBwIDPPHDU+SjMUet0vAnsy82hmTgKPAT/b5TrNSdVC4AfA\nhohYHxGDNAdqnuhynRaE4rGeDwHPZ+afdLs+C0lm3peZazNzHc3/Z76dmZflXt3FlpmHgf0R8Z6i\n6A6aTw7sdfuA2yPiiuLf1h1cpgPmpY+XvJxk5lRE/BbwTZqj9V/JzJ1drtZC8SHg08CPIuKHRdnn\nMvMbXayTLg+/DXyt2LF6Cfj1Lten6zLzqYh4FHiG5pl3z3KZXjnsFcOS1MOq1h0kSXoHDAFJ6mGG\ngCT1MENAknqYISBJPcwQkKQeZghIUg8zBCSph/1/OTkTe+mCtnsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1168c2c10>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Tatm: (26,) \n",
      "  Ts: (1,) \n",
      "The subprocess tree: \n",
      "top: <class 'climlab.model.column.BandRCModel'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#  Set the ozone mixing ratio\n",
    "#   IMPORTANT: we need to flip the ozone array around because the vertical coordinate runs the wrong way\n",
    "#   (first element is top of atmosphere, whereas our model expects the first element to be just above the surface)\n",
    "#col2.absorber_vmr['O3'] = np.flipud(O3_global)\n",
    "\n",
    "#  Not anymore!\n",
    "col2.absorber_vmr['O3'] = O3_global"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 730 steps, 730.4844 days, or 2.0 years.\n",
      "Total elapsed time is 1.99867375676 years.\n"
     ]
    }
   ],
   "source": [
    "#  Run the model out to equilibrium!\n",
    "col2.integrate_years(2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11f279750>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEfCAYAAAC04jrjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VFX6wPHvmwBJCEgLhC5oAKkCQUAFBEUBEbDg6uqi\nu7a1odh7++murriKa9lVF10RFbsioGKhKSq9KSJdQlESeofk/f1xbnQYU2aSSe7M5P08zzyTuffO\nve+ZgfvOPefcc0RVMcYYYwqS4HcAxhhjopclCWOMMYWyJGGMMaZQliSMMcYUypKEMcaYQlmSMMYY\nUyhLEsaYiBGR/iIyR0R2iYiKyFAROSP/74Dt2nnLbvYzXlM8SxJxzvuPGOqjmd/xRisR6Soi94tI\nQ79jiVYiUh94BxDgBmAYMMvXoEypVfI7AFPmhgW97glcATwPzAhat7lcIopNXYH7gAnABp9jiVYn\nAFWBO1R1cv5CEckCUoADfgVmSs6SRJxT1bGBr0WkEi5JfB28rqIQkeqqutPvOAJFW0wiUk1Vd4X5\ntvre85bAhaqaB+yLSGCm3Fl1k/kdEUkUkREiskBE9orIDhH5VERODNru13plERkmIktEZJ+ILBOR\nP3rbHC0iH4jINhHZLiIvikjVoP287dVhNxCRcSKyVUR2i8jHItKukBgvEpFvvPftFpGvRGRQ0DbV\nvPieFpHTReRrEdkNvO6tP1JEnhSRRV58e0VksVf2hID9PAY85b2cHVA993T+eu91WgFxZovIhFBj\n8rapIyKPi8gqETkgIj+LyBgRaRLi9xfy5xn0HV7kfef7gEcCtukiIhO8/ezzPqPrRUQCywk8E/QZ\n7fLW/a5Nopj4i/1uve3O9tbliMgeEVkjIm9ZtWlk2ZWEOYz3H/9tYBAwDlctVRX4MzBVRAao6mdB\nbzsPaAD8B9gOXAm8KiK5wChcFc3twInAX4AdwIigfSQCnwFrgHuAJsA1wAwROU5VVwTEOAq4HhjP\nbyfXPwDjReQvqvq/oH33Ai7y4nsROOgtPw44HfgAWAkke+V+wjv+Td52rwF1vX3cC6z2li/73QcY\nugJj8hLNN97xRgM/eLFcBfQVkUxV3RjC/kP+PD1/AhoB/8ad7HO8eHoCnwK7gKdxVZJn4b7XNsBf\nvfdfifvsAj+jg4Qp1O9WRE7H/TudCzwE7PTiPw1o6pXbRIKq2qMCPXAnewX+XMj6Yd76C4KWJwHf\nAUsClrXztt0G1A9Y3gQ4BOQBVwTtZzKwB6gUsOxtbz+vBG3b01v+dsCyXt6yO4O2FdzJLBtI8pZV\n87ZV4PgCyloVkAKWvwfsB2oFLLvW20+XArZ/zFuXVsC6bGBCwOviYhqNO+G1ClreEtgLPB3CdxzO\n55n/He4Fmhewr4W4qqKWAcsScIlfgW7FfUbAGd7yoQUc9+YSfrfPA7lAdb//T8X7w6qbTLA/Ab8A\nk0UkLf8BVAcmAm3l9z183lDVTfkvVHUdsBZ34vlv0LYzcI2YjQs49j8CX6jqDGAmMFBEKnuLL8Ql\noFeD4quD+/VZB8gM2u9Xqvp18MFUdY/mn4VEkkSktrevj4EqQMcCYoyU38Ukrr3oPNwJMSeofFuA\nebhfyqEK5fPM946qrg5cICJHAR2AN1X1x4D95AEPey/PCiOe4oTz3W7HJauzRCQxgjGYIFbdZIK1\nBupRdE+ndA7v4bOqgG22Age8E0rwcnD/4dcELM+j4Oqb73G9ZhriEk9r3L/bNQVsGxhfoB8L2khE\nkoC7cCenowrYpFYRxyitgmJqCqTiTryFnXx3hLj/UD/PouJp7j1/V8C6/GUFfW4lFc53+zgwAHgZ\n+JeIzAA+Asap6pbC3mzCZ0nCBBPgJ+DSIrYJrs/OLWS7wpbnH6ckBFcVdEYR2ywKer2nkO3+jWsj\neQW4H5cYD+HaTu4n9I4dRU3KUtj/sYJiyv9MPgT+Vcj7DoUYU7gK+4zKU8jfrapuFJGOQG+gL66q\n6mngARE5TVXnl3GsFYYlCRNsOdAdmK6q5dmvPQFoBSwJWt4GVyeef+WyHOgBLPOqtUrE6710ATBJ\nVS8KWldQNVNRiSD/l2ttXL15/n5qAzXCCCsLV0WXqr/vHBCuUD/PouRfIbYtYF2boG0iIazvVlUP\n4RrnPwN3wyPwLa6TxHkRjKtCszYJE2wMrs3ggYJWikhwVU4k3RZ0rJ64qpFJqprfU2aM9/xIYDfV\nEsaXS9AVjYjUBIYXsG3+PQO1C1iXX1XTN2j5TcEbFkVV9wNvAieLSP+CthGRemHsMpTPs6h4VuN+\nuf9BRDIC9pOAOxGDa+SPlJC/24K6G+MS4gEK/o5MCdmVhAn2Mq6u93YROR5Xz7sF12OpJ64toUMZ\nHHcfkCkiE4FJ3vGuxTVQ5p+QUNWpIjISuAVoLSLvAptwdexdcb9EaxZ3MFXNE5H3gAtFZAyuQb0B\ncBnwM659INC33vP9ItIYVz3zo6rOw1UP/QQ8JiKNcFcEfXC/wMO9Qe4moBswUURexw1rkQs0w1XD\nfI77XIoT0ucZgmtxDelfi8izuCuls3Dle15Vvy3qzeEI87t9XURScZ/HT7heYxfiOhyMCd63KTlL\nEuYwqqoicj7uxHAJcDfu38lGYA7wZBkdOhf3S/wJXL/3KsCXuC6Sy4NivFVEvsH1+78Zd+XzM+5X\n7w1hHPMqXAI8C9cXf413/BW43jSBx/xBRK7y9v8foDLufoJ5qnpARAbiPpsbcSfoibgT6fdhxIOq\n5ohIN9yJ8hxgKO7XcRYwFXdPRShC/jyLiWeGdwVyP+7+hRRctdAICm83KbEwvtsXcT3xLgHScN2w\nlwCDVHUCJmLE6wFojG9E5G2gv6pW8zuWeGCfp4kka5MwxhhTKEsSxhhjCmVJwhhjTKGsTcIYY0yh\nYr53U1pamjZr1szvMEpk9+7dpKam+h1GqVk5okc8lAGsHOVh7ty52apat7jtYj5JNGvWjDlz5vgd\nRolMnTqV3r17+x1GqVk5okc8lAGsHOVBRNYWv5W1SRhjjCmCJQljjDGFsiRhjDGmUDHfJlGQgwcP\nkpWVxb590T33eo0aNVi6dGm5Hzc5OZnGjRtTuXLwvDPGGHO4uEwSWVlZVK9enWbNmhEwV3vU2blz\nJ9WrVy/XY6oqOTk5ZGVl0bx58+LfYIyp0OKyumnfvn3UqVMnqhOEX0SEOnXqRP1VljEmOsRlkgAs\nQRTBPhtjTKjisrrJGGPijips3gwrVsDy5e75kkugjKuNLUmUkcTERNq3b8+hQ4do3rw5r7zyCjVr\nuvlSfvzxR0aMGMGyZcuoUaMGGRkZPPXUUyxdupQhQ4Yc1lbw2GOP0bdv8IRnxpi4VFAiCHzeseO3\nbRMT4fjjLUnEqpSUFBYsWADAxRdfzDPPPMNdd93Fvn37GDhwII8//ji9e/emevXqTJ06lc2bNwPQ\ns2dPJkywOVOMiVvhJoJmzSAjwyWEjAxo0cI9N2sGVaqUebiWJMrB8ccfz6JFiwB47bXXOP744xk0\naBA7d7qZLfNv2586dapPERpjIkoVsrM5YskSWLs26hNBUeI+SYxYvpwFu3YVv2EYOlarxqgWLULa\nNjc3l88//5xLL70UgCVLlpCZmVno9jNmzKBjx46/vn7nnXc4+uijSxewMSbyvETA8uWFXhF0zt82\nIcGd8Fu0iMpEUJS4TxJ+2bt3Lx07dmT9+vW0bt2aU089NaT3WXWTMVEkhETwqwISwaK9e+lw9tlR\nnwiKEvdJItRf/JGW3yaxZ88e+vXrxzPPPMN1111H27ZtmTZtmi8xGWMKEJgI8k/+oSSC7t3dcxFX\nBFumToWWLcuzNBEX90nCb1WrVuVf//oXZ555JldffTUXXHABDz/8MBMnTqRXr14ATJ8+ndq1a/sc\nqTFxrAwTQbyzJFEOOnXqRIcOHXj99dcZNmwYEyZMYMSIEVx33XUkJSXRoUMHnnzySbKzs3/XJnH3\n3XczdOhQH6M3JkYUlAgCn7dv/23bghJBfjtBBUwERbEkUUZ2BTWWf/jhh7/+fcwxx/Dxxx//buym\n9PR0tgf+QzbGHM4SQbmzJGGMiS6qVN62DWbOtEQQBSxJGGPKXzFXBCdaIogaliSMMWWjJFVDGRnQ\nvTvLgRYDBlgiiAKWJIwxJZefCArqMVRMIvj1aqCARLB+6lRaeCMRGH9ZkjDGFK2MEoGJDZYkjDG/\nTwTBCSGURJCR4UYktUQQVyxJlJGyHCp8+/btDB8+nJkzZ6KqnHjiiTz11FPUqFGDtWvXctZZZ5GX\nl8fBgwcZPnw4V155ZbmW3UQpSwSmBCxJlJGyHCr80ksvpV27dowZMwaA++67j8suu4y33nqLBg0a\n8PXXX5OUlMSuXbto164dgwcPpmHDhmVbYBMdihuGOjgRHHnkb72GLBGYAliSKAeRHCp8xYoVzJ07\nlzfeeOPXZffeey8ZGRmsXLnysBFj9+/fT15eXuQKYqJDYVcEK1bQY+lS2L37t20DE8Gf/mSJwIQt\n/pPEiBHg/aKPmI4dYdSokDaN9FDh33//PR07diQxMfHXZYmJiXTs2JHvvvuOo48+mnXr1jFw4EBW\nrFjByJEj7SoiFhWRCIq6Ivi5b18a9e5ticBETPwnCZ/4OVR4kyZNWLRoERs2bODMM89k6NChpKen\nl2qfpgyoQk5O4cNQF1Y1dOGFhw86F5AIlk+d6pKEMRESNUlCRO4HLgc2e4vuVNVJpd5xiL/4I62s\nhgpv06YNCxYsIC8vj4SEBADy8vJYsGABbdq0OWzbhg0b0q5dO2bMmGGDBPqltIkgv3rIrgiMT6Im\nSXieUNXH/A4ikiI9VHhGRgadOnXioYce4t577wXgoYceonPnzmRkZJCVlUWdOnVISUlh69atfPnl\nl9xwww1lVj7D4YmgoHsJLBGYGBZtSSIuRXqo8NGjRzN8+PBf2yqOP/54Ro8eDcDSpUu56aabEBFU\nlZtvvpn27duXX2HjVaQSQbNmkJTkWzGMCZeoqt8xAL9WN/0F2A7MAW5S1a2FbHsFcAVAenp65rhx\n4w5bn3/vQbTLzc09rAG6PK1YsSJiw5Lv2rWLatWqRWRfftq1cye18vJIycoiZf36357Xr6dqVhaV\nAnoNaUIC+9LT2duokXs0bszeRo3Y06gR++rXR326Ioib78LKUeb69OkzV1W7FLdduSYJEfkMqF/A\nqruAb4BsQIEHgQaqeklx++zSpYvOmTPnsGVLly6ldevWpQ+4jAXPJ1GeIvkZTZ069dduvFGvoCsC\n7/nQ0qWHJYJfrwgCh5aI8iuCmPouimDlKHsiElKSKNfqJlXtW/xWICIvAKXr4mNMoKVL4fLLYcmS\ngquGMjJc99GTTjp8rKEoTATGlKeoaZMQkQaqutF7eRawpDT7U1VEpPSBxaFoqWIsN4cOwUUXwapV\nro0geNA5LxFY91Fjfi9qkgTwqIh0xFU3rQH+WtIdJScnk5OTQ506dSxRBFFVcnJySE5O9juU8jNq\nFMyZA6+/Duef73c0xsSUqEkSqjosUvtq3LgxWVlZv46HFK327dvny8k6OTmZxo0bl/txfbFiBdxz\nDwweDOed53c0xsScqEkSkVS5cuXDRlKNVlOnTqVTp05+hxG/8vLgssvcvQfPPgt2VWlM2OIySRgD\nwAsvwLRp7rlRI7+jMSYmJfgdgDFlYt06uOUWOPlk8AZXNMaEz5KEiT+qcNVVrlfTCy9YNZMxpWDV\nTSb+vPYaTJwIjz8ORx3ldzTGxDS7kjDx5Zdf4Prr3Uxr113ndzTGxDxLEia+XHcd7NwJo0eDT+Ni\nGRNPLEmY+PHBB/DGG3D33RA0t4YxpmQsSZj4sG2ba6zu0AFuu83vaIyJG9ZwbeLDLbfAzz/D+PE2\ncY8xEWRXEib2ff45/Pe/cPPN0KXYkY+NMWGwJGFi2+7dbgjwFi3g/vv9jsaYuGPVTSa23X03rF7t\nht9ISfE7GmPijl1JmNj1zTfw5JOuwbpXL7+jMSYuWZIwsWn/frjkEmjcGB55xO9ojIlbVt1kYtPf\n/uamJJ00CY44wu9ojIlbdiVhYs/ChfDwwzBsGAwY4Hc0xsQ1SxImthw65Ib+rl0bnnjC72iMiXtW\n3WRiyz//CXPnwptvQp06fkdjTNyzKwkTO157De68E845B4YO9TsaYyoESxImNrz1Flx0EfTsCWPG\n2ERCxpSTElU3iUh9oCGQAmQDq1X1QCQDM+ZX778PF1zg5oiYMAGqVvU7ImMqjJCThIh0AS4D+gFN\ng1YfEJHZwOvAWFXdGbkQTYU2cSL84Q+Qmem6u1ar5ndExlQoxSYJLzk8BvQCFgMTgPnAZmAvUBto\nDnQDHgEeEZFHgX+q6r4yittUBJ98Amef7Yb//vhjux/CGB+EciUxDXgBuEpVlxa1oYgkA0OAW3Ht\nHQ+WOsIykrVvH9/v2cNptWv7HYopyOefw5lnusmDJk+GmjX9jsiYCimUJHG0qm4KZWfelcMbwBsi\nkl6qyMrYfWvW8M7mzfxy4olUSbD2+6gybRoMGgQZGfDpp+6eCGOML4o9O4aaIAp4388leZ+I3CQi\nKiJpJXl/qAanpbE9N5cZ27eX5WFMuL76CgYOhGbN3NVEWpn+MzDGFKNEP6FFpIOIXCsi93k9nRCR\nDBGpXppgRKQJcBrwU2n2E4q+tWqRnJDA+Ozssj6UCdWsWW6YjYYNXYKoV8/viIyp8MJKEiKSJCJv\n4Rqu/wXci+sKC/AocFcp43kC156hpdxPsVITE+lbqxbjc3JQLfPDmeLMmwennQZ168IXX0CDBn5H\nZIwBJJwTpIg8BlwKXAN8CvwMdFHVeSJyOXC1qnYqUSAiQ4CTVfV6EVnj7bfAn/kicgVwBUB6enrm\nuHHjSnJIJgD/BEYDR5VoD6Wza9cuqsVBl87SliN1xQo63ngjuVWrMn/UKPbXrx/B6EIXD99HPJQB\nrBzloU+fPnNVtfj5flU15AewHrjG+zsRyAM6e6/7AluLef9nwJICHkOAb4Ea3nZrgLRQYsrMzNSS\n2rBvnzJlij60Zk2J91EaU6ZM8eW4kVaqcixerJqWptq4seqqVRGLqSTi4fuIhzKoWjnKAzBHQzjH\nhnvHdR2gsG6wCUBSMQmpb0HLRaQ97l6LheKGW2gMzBORrlrChvNQNEhKomv16ozPzuauI48sq8OY\nwvzwA5xyClSuDFOmQPPmfkdkjAkSbsP1auD4QtZ1BZaVJAhVXayq9VS1mao2A7JwVyhlliDyDU5L\nY9bOnWzcv7+sD2UCLV8OJ5/sxmD64gvX3dUYE3XCTRJjgNtF5EKgsrdMRaQPcAPwYiSDKw+DveGm\nJ+Tk+BxJBbJqlUsQhw65XkzHHON3RMaYQoSbJB4FJgKvAFu9ZV/i2ho+VtWnIhGUd0VRLn1T26Wm\n0iw5mfGWJMrH2rXQpw/s2QOffQZt2/odkTGmCGG1SahqLnC+iDyDG+ivHpCDSxDTyiC+MiciDK5T\nh+c3bmR3bi6piYl+hxS/srJcgtixw11BdOjgd0TGmGKEMwpsFeAq4HNVnQHMKLOoytngtDT+tX49\nn23dyhC7w7dsbNjgqphyctwVROfOfkdkjAlByNVN6uaLeAQ36mtc6VWjBjUSE+3u67KSleV6MW3c\n6EZzPe44vyMyxoQo3DaJpfhz31mZqpyQwIA6dfgwJ4dcu/s6smbMcHNBrF/v5oM4vrDOccaYaBRu\nkrgXuMe7ryGuDK5Th80HDzJrxw6/Q4kPqvDss66KqWZN+PZbN/WoMSamhHsz3W1ANWC+N3TGRg4f\nZ0lV9aQIxVau+teuTSURxufkcHyNGn6HE9v274drroHRo92Irq++CvaZGhOTwr2SyAW+xzVarwMO\necvyH3kRja4c1apcmV41ali7RGmtXw8nneQSxN13w/jxliCMiWHhdoHtXUZxRIXBaWmMWLGClXv3\ncnRKit/hxJ6vvoKhQ2HXLnjnHTf1qDEmptmUbAEGeXdff2hXE+F77jl3D0S1avDNN5YgjIkTYV1J\niEiv4rZR1eklD8dfR6Wk0C41lfE5OYxo0sTvcGLD/v20fOwxmDjRTRj06qtQq5bfURljIiTchuup\nFD8hUEzfsjy4Th3+8dNPbD14kFqVKxf/hopswwYYOpSGX38Nd9wBDz4Idse6MXEl3OqmPsDJQY9z\ngZdxc0CcEcng/DA4LY1c4KMtW/wOJbp9/TV06QKLFvHd/ffD3/9uCcKYOBRuw3Vh4zO9KyJPAIOA\nj0odlY+Oq16d9MqVGZ+dzQXp6X6HE53++1+4+mpo0gQmT2azteEYE7ci2XA9EfhDBPfniwQRBqWl\n8dGWLRzIi9kevWXjwAG46iq4/HLXSD17NrRr53dUxpgyFMkk0YoYvk8i0OA6ddiRm8v0bdv8DiV6\nbNrk7p7+z3/gttvcEBu1424YL2NMkHB7N11UwOIqQDvgUuDdSATlt1Nq1SIlIYHxOTn0tRMhzJrl\nurRu3QrjxsF55/kdkTGmnITbu+l/hSzfD7wBXF+qaKJE1cRETq1Vi/HZ2TyZkYE373bF9NJLcOWV\n0LAhzJwJxx7rd0TGmHIUbpIoaKb6far6cySCiSaD09IYn5PD4t276VCtmt/hlL+DB+GGG+CZZ6Bv\nX3cF4d1saIypOMLt3bS2rAKJNgO9aqbx2dkVL0ns3QvnnAMffQQ33wwPPwyVwv09YYyJB2E1XItI\nSxHpGvA6RUQeFpEPReTayIfnn/pJSXSrXr3izX29Zw8MHuwmB3r+eRg50hKEMRVYuL2bngaGBrz+\nG3AT0BB4QkSuiVRg0WBwWhqzd+5kw/79fodSPnbtckN7f/EF/O9/rqurMaZCCzdJHAt8BSAiCcBF\nwG2qmgk8BFwR2fD8Ndirg59QEa4mduyA/v3dTHKvvAIXFdSRzRhT0YSbJGoA+WfMTkAt4G3v9VTi\nbGrTtqmpNE9Ojv85JrZtg3793Oxxr78OF1zgd0TGmCgRbpL4Gcjw/j4NWKmq67zX1XCTEMUNEWFw\nnTp8tnUru3Nz/Q6nbGzZAqeeCnPnwltvwbnn+h2RMSaKhJskxgMPi8hjuLaItwLWtQdWRSqwaDE4\nLY39qnwajwP+ZWfDKafAokXw3ntw5pl+R2SMiTLhJonbgQlAP1zC+FvAusHA5AjFFTV61qhBjcTE\n+Ovl9MsvbvylH35wU4wOHOh3RMaYKBTufRK7gQK7vKjqCaUJREQeBIbgxn/6Bfizqm4ozT4joXJC\nAqfXqcOEnBxyVUmMh7uvN250VxBr1sCECe5vY4wpQLj3SSSISKWgZf1E5CYR6VjKWEaqagdV7Yi7\nWrm3lPuLmMF16rD54EFm7djhdyill5UFJ50EP/3kbpazBGGMKUK41U2vAy/mvxCRK3HzR4wEvhWR\nviUNRFUDz8CpFD8DXrnpX7s2lUR4e/Nmv0MpnbVrXYLYtAk++cT9bYwxRRDV0M/FIrIWd1/EOO/1\nSuBzXCP280B9Ve1T4mBE/oa792I70EdVCzwri8gVePdkpKenZ44bN66khwzZg8CXuBEOG0Ron7t2\n7aJaOQ35kbxxIx1vuIFKu3ax8NFH2dmmTcT2XZ7lKEvxUI54KANYOcpDnz595qpql2I3VNWQH8Be\noKf3dwau/aCD9/o0ILuY938GLCngMSRouzuAB0KJKTMzU8tD1r59Wm36dB20aFHE9jllypSI7atI\nq1erNm6sWquW6pw5Ed99uZWjjMVDOeKhDKpWjvIAzNEQzrHhDsqzA8gfCrS3lxQWea9zgeRiElKo\n1VGvApOA+8KMr8w0Skri/mbNuHnlSsZnZzM4Lc3vkEKzbZvrubRrF0yZAh1L23RkjKlIwm2TmAnc\nLiJnACNwJ/J8GUBWSQMRkRYBL4cAP5R0X2XlukaNaJeaynXLl7MnFm6uO3jQ3Rz344/wzjuWIIwx\nYQs3SdyKu5IYj7tquD9g3XnA16WI5RERWSIii3BVV1E3gVHlhASebdGCtfv387e1UT5quipccw18\n9pkbzfXkk/2OyBgTg8K9T2I50EJE6qhq8N1l1wObShqIqp5T0veWp541a3JRejoj163jovr1aVW1\nqt8hFeyxx+CFF+DOO+Evf/E7GmNMjAr3SgIAVc0RkWoicqSIVPaWLdZCeiPFm0ePPprUxESu+fHH\n/Ib26PLuu3DbbW4u6gcf9DsaY0wMCztJiMgZIjIP1011JW7MJkTkvyJSIYYPTa9Shb81b87n27bx\nZrTdOzF7NvzpT9Ctm5ufOqFEvwOMMQYI/47rM4EPgGzgtqD3rwYujlxo0e2vDRuSWa0aN6xYwY5D\nUTL47dq1MGgQ1K8PH3wAKSl+R2SMiXHh/sy8D3hJVU8DRgWtWwK0i0hUMSBRhH+3bMmmAwe4f80a\nv8OB7dvhjDNg3z6YOBHq1fM7ImNMHAg3SbQG3vD+Dq6M38pv91BUCMcdcQR/bdiQf2VlsWjXLv8C\nOXTItT/88IPr6tq6tX+xGGPiSrhJYgdQ2F1kzYAoq6Ave39r3pxalStz1Y8/kudXI/Yjj7ixmP7z\nHxuwzxgTUeEmiU+BO0SkZsAyFZEk4FrcYH8VSu3KlXn0qKOYuWMHL28qcQ/g0nn/fejZEy691J/j\nG2PiVrhJ4i6gPrAM+C+uyul2YAHQmMNvrqswLq5fnxOPOIJbV61iy8GD5Xvwbdtg/ny7gjDGlImw\nkoSqrgE64+Z7OBU3XlMv4Bugm0bBJEF+SBDh2ZYt2XrwIHeuKucZXKdPh7w86N27fI9rjKkQwrrj\nWkRqAJtV1eo1gnSoVo3rGjdmVFYWlzRoQNcjjiifA0+ZAsnJ0L17+RzPGFOhhHwl4c1Il4MbV8kU\n4P5mzahfpQpX/fgjueXViD1lCpxwAiQllc/xjDEVSshJQlUPAT/jqphMAY6oVIknMjKYt2sX/1y3\nruwPmJMDCxdCnxLP82SMMUUKt+F6LHBZWQQSL/5Qty7npKVx+6pVvF/WQ3Z884177tWrbI9jjKmw\nwp10aA1wgYjMxg3PsZGgm+pU9cUC3ldhiAhjWrfmpwULuGDpUqYnJdGlrNonano9kXfuLJv9G2Mq\nvHCTxDMA1XPEAAAgAElEQVTecyMgs4D1ClToJAFQNTGRD9u3p9vcuZyxeDHfZmZyZHKRk/aVTDtv\nFJTFi93sc8YYE2HhVjc1L+ZxVESji2HpVaowqUMH9uXlMXDRIraXxSCANWpA06YuSRhjTBkId9Kh\nKJ+OLbq0SU3lnXbt6L9oEUO/+45J7dtTOdJDd7drZ0nCGFNmQjpjiUgfEXnfm170SxEZXtaBxYtT\natXi+ZYt+WzrVq4qi0mK2rd3A/uV953expgKodgkISL9gM+Ak4DduCqlUSJyZxnHFjf+0qABdzVt\nyuhNm/jHTz9FduedO7sEMXNmZPdrjDGEdiVxJ/Al0FRVuwFNcOM23SYiUpbBxZMHmzfnj/Xqccfq\n1bzxyy+R2/EZZ0CtWvDss5HbpzHGeEJJEq2Bf6rqTgBVzQUeBKrjEoYJgYjwYqtW9KhRg4uXLuWr\n7dsjs+OqVeEvf3HzWm+okENnGWPKUChJIg0IHgN7o/dcoSYZKq3kxETea9uWJsnJDFm8mPWR2vFV\nV0FuLjz/fKT2aIwxQOhdYH2aTSf+pFWpwqT27QG4A8iJRINzRgb07++ShDVgG2MiKNQk8aGI/JT/\nAFZ7yycFLhcR6yIbghZVq/J+u3ZsAvovWkT2gQOl3+k118DGjW76UmOMiZBQ7pN4ucyjqIB61KzJ\nA8D/7d5NzwULmNyhA01Kc1d2//7Qpg3ccQcMGgSpqRGL1RhTcRWbJFT1L+URSEV0PPBJhw4MWryY\nE+fPZ3KHDhxT0pN7YiI895ybxvTee+Gf/4xorMaYiinCt/+WnIiMFJEfRGSRiLwXNI923OpVsybT\nOnZkf14ePebPZ/aOHSXfWY8e8Ne/wqhRMGdO5II0xlRYodxM1zncnYpIsogcE+bbPgXaqWoH4Edc\nu26F0LF6db7q1InqlSpx8sKFfLZlS8l39sgjkJ4Ol19ujdjGmFIL5UpiuoiMF5H+IlLk9iLS1LsT\nezVwRjiBqOpkb2IjcHNmNw7n/bEuo2pVvurUiWbJyQxcvJi3S3rDXc2a8NRTsGABPPFEZIM0xlQ4\nUtxYQiLSCHfz3IXADuBrYCGwGdgP1MIN1dEVaIdLEPep6mslDkrkQ+ANVR1byPorgCsA0tPTM8eN\nG1fSQ/lq165dVKtW7bBlO3GXUN8DNwCDSrJjVdrdcw+1Zs9mzujR7G1ctvm2oHLEongoRzyUAawc\n5aFPnz5zVbVLsRuqakgPoB5wG/AFbgynvIDHSuAloD9e4ilkH58BSwp4DAnY5i7gvaL2E/jIzMzU\nWDVlypQCl+8+dEgHLFyoTJmif1+zRvPy8sLfeVaWaq1aqt26qR48WLpAi1FYOWJNPJQjHsqgauUo\nD8AcDeEcG/JQ4ar6C/AP74HXsJwM5KhqSJXfqtq3qPUi8mdcNdUpXiEqpKqJiXzQrh1/+eEH7ly9\nms0HD/LY0UeTEM5QWY0awb//DeefDw8/DPfcU3YBG2PiVol7N6nqNlXdFGqCKI6I9AduBQar6p5I\n7DOWVU5IYEzr1gxv1IgnsrIYvHgxW8NtiD7vPLjgAnjgAZg9u2wCNcbEtajpAgs8jRs08FMRWSAi\n//E7IL8liPBkRgbPtGjB5K1b6Tx3LvPCnc/66aehQQMYNgz2VPjca4wJU1hJQkTyRCS3kMchEckR\nkU9F5LRwA1HVDFVtoqodvceV4e4jHokIVzdqxPSOHTmkygnz5jF648bi35ivVi343/9g2TK47bYy\ni9MYE5/CvZJ4EFiH69n0P1z7xMve6yzgFaAu8JGIhNUF1hSte40azMvMpGfNmly2bBmX/fAD+3Jz\nQ3vzKafAiBHuquKDD8o2UGNMXAk3SezDdXFtpqqXquqdqnoJ0BxYg0sWnYHJuMmKTATVrVKFjzt0\n+HWWuxPnz2f13r2hvfnhhyEzEy6+GFauLNtAjTFxI9wkcSXwhKruC1yoqnuBJ4ArVTUPN3Ndh8iE\naAIlivDQUUfxYbt2rNq3j85z5zIxJ6f4NyYnw9tvQ0ICnHMOhJpcjDEVWrhJoi5QuZB1VfhtEqJs\nwKY2LUNnpKUxNzOTZsnJnLF4MfeuXk1ucb2GmzWDsWNh4UK49tpyidMYE9vCTRJzgftFpEHgQhFp\nCNwH5I8qdyRgc2mWsaNSUpjZqRN/qV+fB9eupe/ChWzYv7/oN51+Otx9N7z4IoweXT6BGmNiVrhJ\n4nrcmEqrRGSKiLwhIlOAVUBD4DpvuwygxMNymNClJCYyulUrXmrVilk7dnDsnDlMKq766f77oW9f\nN1HRvHnlEqcxJjaFlSRUdR4uATyBG46jvff8T6CFqi7wtrtXVe+LcKymECLCnxs0YG5mJg2rVGHg\n4sXctGIFB/LyCn5DYiK89hrUqweDB8MGu+gzxhQs7JvpVDXH69V0iqq28Z7vUtUQWk9NWTomNZVv\nOnfm6oYNeTwrixPnz2dlYQ3UdevChx/Ctm0wZIjdaGeMKVCJ7rgWkdoiMlBEhonI6SJSO9KBmZJJ\nSUzkmZYteadtW1bs3UunOXN4/eefC9742GPh9ddh7lx3R3ZhVx7GmAor7CQhIg8B64EPcTfSTQDW\ni8iDEY7NlMLZdeuyoEsX2qWmcsHSpVz6ww/sKejmu0GD3FSn777rGrSNMSZAuMNyjMDdJDcW6AO0\n9p7HAneKyHVFvN2UsyOTk5nWsSN3Nm3KS5s20X3ePH4sqFppxAg37enDD7shPIwxxlOSm+meVNXL\nVXWaqi7zni8H/gVcHfkQTWlUTkjgb0cdxaT27dmwfz9d5s7lreBZ70TcbHZ9+8IVV8AXX/gTrDEm\n6oSbJJoBEwtZN9Fbb6JQ/zp1mNelC22qVuUP33/PiOXLD+/9VLkyvPUWtGrlGrJtaHFjDOEniRzc\nFKUFaeutN1GqaXIy0zt14rpGjXhy/XpOWrCAdfsCRlipWRM++cT1fBowAL7/3r9gjTFRIdwk8R7w\noNerqRKAiFQSkT8C/we8E+kATWRVSUjgyRYteKNNG5bs3k2n4JvvGjaETz91VxannQZr1vgWqzHG\nf+EmiTuABbheTXtF5GdgL/AqsBAb+TVm/KFePeZkZtIgKYmBixdzzpIlrM2/qjj6aJg8GXbvhlNP\nhcK60Bpj4l64d1zvBHoBg3F3XY8HHsfNS32Squ6KeISmzLSqWpXZnTvzUPPmfLRlC8fMmsX/rVnD\n3txcaN8eJk1yd2P36+duujPGVDglueNaVXWCqt7q9XK6TVUnqRY3BKmJRsmJidx15JH80LUrg+rU\n4b41a2gzezYfZGej3bvDe++5tonMTFcNZYypUIpNEsVMWfq7KUzLI2gTeU2Tk3mzbVs+P/ZYUhMS\nOHPJEgYsWsSyHj3gs8+gUiXXRnHhhRDchdYYE7cqhbDN/wF2lVBBnFyrFvO7dOGZ9eu5b80a2s6a\nxam1a/PHjz/m3NGjSXn0UfjoIxg5Ei65xN1jYYyJW8UmCVW9vxziMFGkckICI5o04Y/p6Tyxbh1v\nbN7MxVu2cNmpp3Lpccdx/8MPk37ZZTBmDDz3nN/hGmPKUIkG+DMVQ3qVKjxy9NGs6taNbzt35rpG\njZhYrx4N/v53rrzlFnYuWEBuhw40eukl2Lev+B0aY2KOJQlTLBGh6xFH8FhGBmu6d+fLzEySL7+c\nHq++yrhevWgxZgzrW7dm5n//y+6dO/0O1xgTQZYkTFgSRDihRg1GtWjB/IEDafr229z/6KMcPHCA\nEy6/HOrVY2b//i5h7Njhd7jGmFKyJGFKLEGEnjVr0vu442iyZg3z33+fuWeeSYtZszjh8suRevX4\n+rTTmPncc+yy+yyMiUmWJExEJFauTKchQ+j1+uvU/uUXFowfz+yhQzlq3jxOuPJKEtPT+faUU5j5\n73+zc8sWv8M1xoQoapKEiJwrIt9592V08TseU3KJlSrRcdAgTho7lrRNm1gwcSLfnnceTRct4oSr\nr6ZK/frM6tOHmU8/zfbsbL/DNcYUIWqSBLAEOBuY7ncgJnISK1Wi4+mn03vMGNI3bWLhxx/z9YUX\n0vj77zlh+HBSGjRgVq9ezBg1ii12k54xUSdqkoSqLlXVZX7HYcpOQmIix/brR++XXqL+hg0s/vRT\nZl58MQ2XL6fnDTdQvWFDZvXowbRHH+WXrCy/wzXGABJtQy6JyFTgZlWdU8Q2VwBXAKSnp2eOGzeu\nnKKLrF27dlGtWjW/wyi10pZDVdm6bBlVpk6l/YwZHLlhA7kJCcxv1441PXuS3LMn1dLTIxdwIeLh\n+4iHMoCVozz06dNnrqoWX7WvquX2AD7DVSsFP4YEbDMV6BLqPjMzMzVWTZkyxe8QIiKS5cjLzdUf\nv/5ap1x3nS7LyFAFVdDv2rTRKbfcoqvmzYvYsYLFw/cRD2VQtXKUB2COhnCODWXspohR1b7leTwT\neyQhgRbdu9Oie3d48klWLVrE2tdfJ23CBHqPHAkjR7KyWTOyBg6k/nnn0fLEE5GEqKk1NSbu2P8u\nE9WO6tCBPg8/TPvFi1m/fDlTH3yQ7Wlp9Pj3v2nVqxfrGzVi2iWXsOjjj8k9ZIMQGxNpUZMkROQs\nEckCjgcmisgnfsdkokujjAx63303nWfPZuu6dUx//HHWt2xJ97Fj6TBgADnp6Uw7/3xmv/02+20s\nKWMiolyrm4qiqu/h5tA2plhpDRvS64Yb4IYb2LFlC3Peegt57z0yx4+n2htvsK1aNeb06oUMGUK7\noUM5onZtv0M2JiZFzZWEMSV1RO3anPjXv3LCxx9TefNm5owdy5J+/Wg1cyYn/PWvJKenM+fEE5n+\nyCP8vHq13+EaE1MsSZi4kpSaSpcLL6TH229Ta/NmFk2axFcXX0zaTz/R6447SD/qKL5r25Zpt97K\nytmz0bw8v0M2JqpFTXWTMZGWWKkSHQYMgAED0Lw8ls+bR9abb1J30iRO8npKrWncmLX9+lFr6FBy\nK1f2O2Rjoo4lCVMhSEICLbp0oUWXLvDoo2xYtYrlb75J1QkTOP7ll6kyejSba9ZkxsknU+Wss2h/\n5plUjdKboIwpT1bdZCqkhkcdxUm3385xX37J3p9/5qvnnmNhZiYdJk+m27BhkJbG16ecwvRRo/hl\n/Xq/wzXGN5YkTIVXo3ZtTrziCirdfTcpmzcz7513mDN0KEcuWUKvG24grUkTFh97LFNvvZUVs2ZZ\nO4apUCxJGBOgSnIync8+m15jx9Jg40Z+nDGD6ddfT+LBg/QeOZKMbt3IatyY6RddxPx33+Wg3Y9h\n4py1SRhTCElIoGWPHrTs0QOAjatXs/ytt0ieOJHj3niDlFdeYUfVqszt0YPcgQNpfe651G7QwOeo\njYksu5IwJkQNmjen16230nXaNPKys/lm7FjmDxxIs7lzOfH666nRuDGLOnVi6h13sHLOHKuWMnHB\nkoQxJZBavTrdL7yQk958k3o//8ySKVOYfs01VN67l96PPMLRxx3HT02bMu3Pf2bue+9xwKqlTIyy\n6iZjSikhMZF2vXtD794ArF+xghVvv03KpEl0e+01kl9+me2pqczu0YO800+n1dlnU69xY19jNiZU\ndiVhTIQ1ysjgpNtvp+v06eRmZzNrzBgWnX46GXPm0PP660lr2pQl7dox9cYbWTp9Onm5uX6HbEyh\nLEkYU4ZSjziCrsOG0dOrllo2fTrTR4wAEXqNGkXrk05ic3o6X55zDt+++CI7t2zxO2RjDmPVTcaU\nE0lMpFXPnrTq2ROAzevXs+zdd0mYNIl2kydT8913OXjFFczv3Jkd/frR9OyzaXbssTapkvGV/esz\nxid1GzWix/DhnPDRR6Tm5LBg4kS+vPRSUrds4aSHHqJ5586u8XvYMOa8+SZ7d+/2O2RTAdmVhDFR\noHKVKnQ8/XQ4/XQA1i1bxqr33iPlo4847q23qDp2LHuSkvi2Wzf2DhjAUWedRdNWrXyO2lQEdiVh\nTBRq0qqVa/yeNg3Jzmb2uHHMPvdcGqxYQe877qDpMcewslkzpl5yCXPfeYf9e/f6HbKJU3YlYUyU\nS6lWjePOOw/OOw9UWbNgAT+9/z6pn35K97FjSX7pJXYnJ7Owa1c2d+rEuvR0mrRu7XfYJk5YkjAm\nlojQrFMnmnXqBA88wJ4dO5g9cSJ7Jk6k2dSpdJ0+HZ58ktVNm/JTnz6knnEGbQcMICU11e/ITYyy\nJGFMDKt6xBEc98c/wh//iObl8f6rr1Jz9WpSJ0+m6+uvk/Lyy+xJSmJW167sOe00jjzzTJq1bYuI\n+B26iRHWJmFMnJCEBGo2aULve+/luC+/hOxs5owbx5w//IF6q1fT+557aN6+PWubNGHahRfyzauv\nsnPbNr/DNlHOkoQxcSqlenW6nHcevcaModm6daxbvJgZDz3E5qOPJvO99+j+pz+RVLcu8487jim3\n3cbSr74iL4KDEu7fu5ev/v1vfmjZkr1JSeQlJLA3KYkfWrZk5n/+Y43tMcKShDEVRJN27eh5110c\nN20aVXJyWPDBB8y85BKqbt9On0cfpXWPHmTXrcuXQ4bw5bPPlmpGvukjR7K/bl063HQTxyxfTsqB\nAySoknLgAMcsX077G29kf926TB85MoIlNGXB2iSMqYCqpKTQcfBgGDwYgF/WrOHH8eNJ+OQTWk+b\nRp3x48m79lq+P+YYNvXuTa3TT6dN374kJScXu+8pw4fT9YUXSN2/v9BtqntXEZn33MOUn36iz1NP\nRaZgJuLsSsIYQ71mzehx3XWcMHEitbKzWTplCjNuvJFDSUn0eu45Og0axMHatZl10klMfeghVi5Y\nUOB8GdNHjiw2QQRK3b+fri+8YFcUUcyuJIwxh0moVInWvXvT2hv6fOeWLSydOJF9H33EkTNmcOQ9\n98A997Cufn1W9exJ5X79OOaMM0g94gg6PvBAyAkiX+r+/XR84AEODB9OlRCuVEz5irorCRHpLyLL\nRGSFiNzudzzGVHTVa9em67Bh9HrtNY5ct451S5Yw4+9/Z0ObNnSaOJETLruMGg0bsrFFC5LCTBD5\nJC+POS+/HOHITSREVZIQkUTgGWAA0Ab4o4i08TcqY0ygJm3b0vOOO+j2+edU3bKFJR9/zJfDh1M3\nO5ukQ4dKtM/qe/dS+/HHIxypiYSoShJAV2CFqq5S1QPAOGCIzzEZYwpRKSmJdv36cdKoUSSolmpf\nTdeujVBUJpJES/nFRpKIDAX6q+pl3uthQDdVvTZouyuAKwDS09Mzx40bV+6xRsKuXbuoVq2a32GU\nmpUjevhZhl4nn1yqRJEnwvQvvgDi47uA6C5Hnz595qpql+K2i8mGa1V9HngeoEuXLtrba2CLNVOn\nTiVWYw9k5YgefpZhb+XKpBw4UOL376tS5dfY4+G7gPgoR7RVN60HmgS8buwtM8ZEubVHHlmq9/9U\nyvebshFtSWI20EJEmotIFeB8YLzPMRljQrDlxhvZmZJSovfurFqVLTfeGOGITCREVZJQ1UPAtcAn\nwFLgTVX9zt+ojDGhyLz4YrSE83GrCF0uvjjCEZlIiKokAaCqk1S1paoerap/8zseY0xoklJSWHDf\nfexOSgrrfbuTklhw3312I12UirokYYyJXb1uuYVZl18ecqLYnZTErMsvp9ctt5RxZKakLEkYYyKq\nz1NPMffBB9mRmlpoG8XOlBR2pKYy98EHbXC/KGdJwhgTcb1uuYWkzZtZ/Pjj/NCyJXuSksgTYY83\nn8TiJ54gOTvbriBiQEzeJ2GMiX5JKSmccOWVcOWVvy6rChzjX0imBKLqjuuSEJHNQKzez58GZPsd\nRARYOaJHPJQBrBzl4UhVrVvcRjGfJGKZiMwJ5bb4aGfliB7xUAawckQTa5MwxhhTKEsSxhhjCmVJ\nwl/P+x1AhFg5okc8lAGsHFHD2iSMMcYUyq4kjDHGFMqShDHGmEJZkigjItJERKaIyPci8p2IXO8t\nry0in4rIcu+5VsB77hCRFSKyTET6+Rf9b4oox0gR+UFEFonIeyJSM+A9MVOOgPU3iYiKSFrAsqgq\nR1FlEJHh3vfxnYg8GrA8qsoARf6b6igi34jIAhGZIyJdA94TjeVIFpFZIrLQK8cD3vKY+j9eLFW1\nRxk8gAZAZ+/v6sCPQBvgUeB2b/ntwD+8v9sAC4EkoDmwEkiM4nKcBlTylv8jVsvhvW6CG55+LZAW\nreUo4rvoA3wGJHnr6kVrGYopx2RggLf8dGBqlJdDgGre35WBb4HusfZ/vLiHXUmUEVXdqKrzvL93\n4ubHaAQMAV72NnsZONP7ewgwTlX3q+pqYAXQFZ8VVg5Vnaxu/g+Ab3CzCEKMlcNb/QRwKxDYiyPq\nylFEGa4CHlHV/d66X7y3RF0ZoMhyKHCEt1kNYIP3d7SWQ1V1l/eysvdQYuz/eHEsSZQDEWkGdML9\n0khX1Y3eqk1Auvd3I2BdwNuy+O0kFhWCyhHoEuAj7++YKoeIDAHWq+rCoM2iuhxB30VLoKeIfCsi\n00TkOG+zqC4D/K4cI4CRIrIOeAy4w9ssasshIokisgD4BfhUVWP6/3hBLEmUMRGpBrwDjFDVHYHr\n1F2DxkQf5MLKISJ3AYeAV/2KLRyB5cDFfSdwr69BhamA76ISUBtX1XEL8KaIiI8hhqSAclwF3KCq\nTYAbgNF+xhcKVc1V1Y64K+muItIuaH3M/B8vjCWJMiQilXH/CV5V1Xe9xT+LSANvfQPcLxCA9bi6\n8XyNvWW+K6QciMifgTOAC73/DBBb5TgaVze8UETW4GKdJyL1idJyFPJdZAHvetUfs4A83MByUVkG\nKLQcFwP5f7/Fb1UxUVuOfKq6DZgC9CcG/48Xye9GkXh94Bq1xgCjgpaP5PBGrUe9v9tyeKPWKqKg\nUauIcvQHvgfqBi2PqXIEbbOG3xquo64cRXwXVwL/5/3dElelIdFYhmLKsRTo7f19CjA3Wr8LL666\nQE3v7xRgBu5HU0z9Hy+2nH4HEK8PoAfuMnMRsMB7nA7UAT4HluN6pNQOeM9duB4Py/B6efj9KKIc\nK7yTUf6y/8RiOYK2+TVJRGM5ivguqgBjgSXAPODkaC1DMeXoAcz1TqTfAplRXo4OwHyvHEuAe73l\nMfV/vLiHDcthjDGmUNYmYYwxplCWJIwxxhTKkoQxxphCWZIwxhhTKEsSxhhjCmVJwkScN5pqcY81\nfsfpJxHpKyJRe6e3iJwgIrtEpF7Asm9E5LMCth3pfad3e6/PF5EsEUkpz5hN2bAkYcrC8UGPTbhR\nVgOXneVbdNGhL9E9HMhjuHtffilsA3H+BdwM3KyqD3mr3gR24oY+MTGukt8BmPijqt8EvhaR/UB2\n8PJ4IiKJuOmADxW7cdnFIEBlVT1Qyv2ciEvkFxVzrOeAy4BrVfWZ/HWqmiciLwC3ishjqnqwNPEY\nf9mVhPGdV/Uy1ave2CUiE0WkddA234jIZyIySNxER3u9iWkyRaSyV+Xxs4jkiMgLgVUdInKMVx1y\nmYg8JSLZIrJbRD4QkSZBxxERuUZEFovIPhH5RUSeE5EaAdske/u7V0TuEZG1wAGghYikisi/xE2o\ns1tENojI+yLSIuD9jwC3AYkB1W/7vHX9vdfdg+K60lteP2DZJhH5r7fuR+AgbjgLRKS6iPxTRNaK\nyAERWSkit4Y48N9lwCxVXVHI95UAvARcClwRmCACjMONfjoohOOZKGZXEsZXInI2bjC394ALgETc\nENHTRaSD/jbkMrhJWx4EHgL24apEPsANfXAA98u3A/AwsJHfV+fcB8z2tmvobfeRiByrqrneNk8A\nV3vPn+MGZPsb0EZETlLVvID9/RU3vMIIL55fgKq4cXwexM2HkAZcC3wtIq1UNQd4xjv+BbihKMAN\nylcSA4DjgHuAHGCFiFTxPpPmXhxLgRNxn1sN3NAQRekHvFbIukq4EX/PBS5W1bEFbaSqG0RkJW6M\nr3cL2sbECL/HBbFH/D9wYyKNLWB5Am78p0lBy2sD23AT6eQv+wZ3Im4SsOwPuDGAJgS9fxKwNOD1\nMd5283FVQvnLT/GWX+i9bok7Wd8atL/87fp7r5O912uAKsWUPRE3+9o+4KqA5Y8AhwrYvr+37+5B\ny6/0ltcPWLYJV/efFrTt5V45ugUtfxDYizcoXSHxHukdZ1gB677x1ilwZwjf+1vAIr///dmjdA+r\nbjJ+aosbLnmsiFTKfwA7cL/4ewVt/52qBk7a8oP3/EnQdj9w+JDM+d5S7+wFoKqfA9m4+ndwv6AF\neDUonunA/gLimaQF1P+LyIUiMltEtuPmrNiBG/mzVQExldYMVc0OWtYfNyXo3KByTMYluKJmQ2vo\nPW8uZP13uGlebxCRtsXEtjlgfyZGWZIwfsrvXvkqrj498NEXN5pmoK1Brw8UsTy5gOP9XMiy/NnB\n8uPJCorlAO4kHxzPxqDXiMi5uBFZFwDnA91w1UHbC4mptH4XA64crfj9ZzrdWx9cjkD5Me4vZP0m\n4GRv/eciUlTi24urejMxzNokjJ9yvOeb+O0EFmhfhI+XXsiyqUHx9AZ2F7Bt8K/rgoZQPh9YoqqX\n5y8Qkaq4toBQ5Je5StDywk7sBcWQg2sr+VMh71lVxPHzP4NahW2gqqtE5GRgGi5RnKSqKwvYtDbu\nSs3EMEsSxk+LcY27rVX18XI43rki8nB+lZOInIJrWP7aWz8Zd9JtrKolnY61Kq6KKdCfC9huP653\nU2U9vIvoWu+5HYcnztPDiOFjXJXT1kJO3kVZiYv/qKI2UtUfRaQvbja2L0Skl6quDdqsOS5ZmRhm\nScL4RlVzReRa4C3v1/Y7uF+y9XG9cX5U1acjeMg04B0R+S/QANe76Ttcd01U9XsRGQU8L26u4hm4\nk3lT4DTgKVWdWcwxPgZGicg/cEmnG67ReVfQdt97z7eIu4v5kKrOU9XVIvItcK/XprEFl2QaEbqX\ncFOBThGRf+ImxEkCMoDBQD/9rTfXYVR1t4jMpeh2i/xtvxORU4EvcFcUvVR1A/x630gm8I8w4jZR\nyCSOu80AAAEcSURBVNokjK9U9T2gD65qYjSuEfoR3Al9VoQP9wBuTuExwFO43joDAk+YqnojMBzX\nJvI28D7ujuLNwOoQjvE08Cium+14bz8D+X311TvAC8CNXhyByec83AxzzwIv4hriR4ZaSFXdj+uR\nNQa4Btfb6xXgQtzVSXHdbd8AThWRpBCOtRDX4J+Gu6LIr9LrDaR6+zIxzGamM3FPRI7B3SswTAvp\n129+IyK1cV2TL1bVt0u4j5dw1XanRjQ4U+7sSsIYcxhV3QI8jrsrPGzeXex/BO6OZFzGH9YmYYwp\nyD+AQyJST4sY5K8QRwLDVfXbMojLlDOrbjLGGFMoq24yxhhTKEsSxhhjCmVJwhhjTKEsSRhjjCmU\nJQljjDGF+n9PTZ7mesNaqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11f184190>"
      ]
     },
     "metadata": {},
     "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Tatm: (26,) \n",
      "  Ts: (1,) \n",
      "The subprocess tree: \n",
      "top: <class 'climlab.model.column.BandRCModel'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "col3 = climlab.process_like(col2)\n",
    "print col3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Let's double CO2.\n",
    "col3.absorber_vmr['CO2'] *= 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The radiative forcing for doubling CO2 is 1.390281 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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1095 steps, 1095.7266 days, or 3 years.\n",
      "Total elapsed time is 4.99668439189 years.\n"
     ]
    }
   ],
   "source": [
    "col3.integrate_years(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  4.31207667e-07])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col3.ASR - col3.OLR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Equilibrium Climate Sensitivity is 2.790202 K.\n"
     ]
    }
   ],
   "source": [
    "print 'The Equilibrium Climate Sensitivity is %f K.' % (col3.Ts - col2.Ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "climlab Process of type <class 'climlab.model.column.BandRCModel'>. \n",
      "State variables and domain shapes: \n",
      "  Tatm: (30,) \n",
      "  Ts: (1,) \n",
      "The subprocess tree: \n",
      "top: <class 'climlab.model.column.BandRCModel'>\n",
      "   convective adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>\n",
      "   LW: <class 'climlab.radiation.nband.FourBandLW'>\n",
      "   H2O: <class 'climlab.radiation.water_vapor.ManabeWaterVapor'>\n",
      "   SW: <class 'climlab.radiation.nband.ThreeBandSW'>\n",
      "   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>\n",
      "\n"
     ]
    }
   ],
   "source": [
    "col4 = climlab.process_like(col1)\n",
    "print col4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'CO2': Field([ 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,\n",
       "         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,\n",
       "         0.00038,  0.00038,  0.00038,  0.00038,  0.00038,  0.00038]),\n",
       " 'H2O': Field([  5.00000000e-06,   5.00000000e-06,   5.00000000e-06,\n",
       "          5.00000000e-06,   5.00000000e-06,   5.00000000e-06,\n",
       "          6.38590233e-06,   9.08848690e-06,   1.33273826e-05,\n",
       "          2.34389689e-05,   3.84220914e-05,   5.95564299e-05,\n",
       "          8.82144990e-05,   1.25843839e-04,   1.73951159e-04,\n",
       "          2.34088411e-04,   3.07840683e-04,   3.96815735e-04,\n",
       "          5.02635028e-04,   6.26926041e-04,   7.71315753e-04,\n",
       "          9.37425100e-04,   1.12686431e-03,   1.34122899e-03,\n",
       "          1.58209684e-03,   1.85102493e-03,   2.14954752e-03,\n",
       "          2.47917415e-03,   2.84138824e-03,   3.23764591e-03]),\n",
       " 'O3': Field([ 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",
       "         0.,  0.,  0.,  0.])}"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col4.absorber_vmr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The radiative forcing for doubling CO2 is 4.421099 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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integrating for 1095 steps, 1095.7266 days, or 3.0 years.\n",
      "Total elapsed time is 4.99668439189 years.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ -5.32160527e-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": {
    "collapsed": false
   },
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}